Discovery of 40.5 ks Hard X-ray Pulse-Phase Modulations from SGR 1900+14 

The demodulation analysis was first applied to the Suzaku XIS data, using the 8–12 keV range where the HXC dominates. Although the XIS energy response is poorly calibrated at $\gtrsim 10$ keV, we included the 10–12 keV interval because it still contains usable HXC signals, and the calibration uncertainties does not affect timing studies. We used the data from XIS0 and XIS3 (front-illuminated CCDs), but not those of the XIS1 camera (back-illuminated CCD chips), because its background at $\approx 10$ keV is more than an order of magnitude higher than those of the other two cameras (Koyama et al., 2007). The maximum harmonic number of $Z_{m}^{2}$ was tentatively set at $m=2$, because any Fourier component with $m\gtrsim 3$ of the pulse profile would be strongly smeared out by the XIS time resolution. As above, $T$ was varied over the 10–100 ks range, with a step of 0.2 s to 1.0 s (depending on $T$). At each $T$, we scanned $A$ from 0 s to 1.5 s with a step of 0.1 s, $\psi$ from $0^{\circ}$ to $360^{\circ}$ with a $5^{\circ}$ step, and $P$ over the error range of Equation (2) with a step of $20~{}\mu$s.

Figure 5 (a1) presents the maximum value of $Z_{2}^{2}$ achieved at each $T$, when $A$, $\psi$, and $P$ are all optimized. After Paper I, this kind of plot is hereafter called a demodulation diagram (DeMD). The result reveals a clear peak at

where the error is estimated as the range where $Z_{m}^{2}$ decreases by 4.72 from the peak value (Paper I). As shown in the same figure with a dashed gray curve, this peak has $A\approx 1.3$ s, or $\approx\pm 25\%$ of $P$. For reference, the XIS0 and XIS3 data, when analyzed separately, consistently reveal the 40 ks peak. When the XIS1 data are included, the peak becomes somewhat higher, but the DeMD becomes noisier in the $T=10-20$ ks interval, presumably due to the background variations.

As indicated by a green arrow in Fig. 5 (a1), the XIS data give $Z_{2}^{2}=12.29$ before the demodulation. Therefore, the 40 ks peak in the DeMD, with $Z_{2}^{2}=26.36$, yields $\delta Z_{2}^{2}=14.07$, where $\delta Z_{m}^{2}$ denotes a relative increment in $Z_{m}^{2}$, and provides a measure of the pulse-significance increase owing to the demodulation (Paper I). As a fiducial value, an increment by $\delta Z_{m}^{2}=10$ means a decrease in the probability $Q$ by a factor of $0.67\times 10^{-2}(1+\delta Z_{m}^{2}/Z_{m}^{2})^{m-1}$, or approximately by two orders of magnitude.

In Fig. 6 (a) we compare two periodograms, both computed using the 8–12 keV XIS 0+3 data. The black one is before the demodulation, whereas the red one is after the demodulation employing Equation (4) and $T=40.0$ ks, together with $A$ and $\psi$ as given in the figure. The result visualizes that the demodulation selectively enhances $Z_{2}^{2}$ at $P\approx P_{\rm XIS}$, although it is not necessarily obvious whether this peak is statistically significant.

To the background-inclusive HXD-PIN data, we applied the demodulation analysis with the same procedure as for the XIS, except that the search step in $\psi$ is reduced to $3^{\circ}$ and $m=4$ is employed. The latter is because the HXC pulse profiles of 4U 0142+61 and 1E 1547.0$-$5408, though variable, often exhibit three to four peaks per cycle when demodulated, and hence $m=4$ has generally been found most appropriate (Makishima et al. 2019; Paper I). In addition, we tentatively chose an energy range of 16–50 keV. The derived DeMD is presented in Fig. 5 (a2), where a prominent peak with $\delta Z_{4}^{2}=24.07$ is found at

This $T$ is close to Equation (5), and the best-estimated amplitude of $A\approx 1.2$ s is consistent between the XIS and the HXD. The parameters characterizing this DeMD peak are summarized in Table 1 in comparison with those from the XIS.

The 16–50 keV HXD periodograms, before and after the demodulation, are shown in Fig. 6 (b) in black and red, respectively. The demodulation parameters employed in calculating the red periodogram are given in the figure. Thus, by correcting the data for the phase modulation with $T=43.5$ ks (Equation 6), the HXD pulses, which were undetectable in the raw data, have been clearly restored with $Z_{4}^{2}=40.64$, at a period of

The error is estimated in the same way as in Equation (5).

The HXD is a low-background but non-imaging instrument, and we are analyzing its data without subtracting the background which amounts to $\approx 91\%$ of the 16–50 keV events. Therefore, we must examine whether the pulse-phase modulation is an artifact caused by background variations. In the first 1/3 of the present observation, the spacecraft was in such orbits as to pass through the South Atlantic Anomaly, where the background is higher and more variable than in the rest (Kokubun et al., 2007). We hence divided the HXD data into three disjoint time portions with comparable durations, and applied the demodulation analysis to them individually, with the modulation period fixed at $T=43.5$ ks because it cannot be constrained. Then, the three portions gave $(A,\psi)$=(1.0 s, 213^∘), (1.4 s, 207^∘), and (1.1 s, 216^∘) in this order, together with $P$ which agrees with Equation (7). The parameters are thus similar among the three time portions, without correlation to the background behavior. Considering further the XIS vs. HXD similarities in $T$ and $A$, we conclude that the HXD results are not much affected by the background variations. We also infer that the 43.5 ks phase modulation of the HXD pulse is a coherent phenomenon, because the three portions, each covering about one modulation cycle, indicate consistent values of $\psi$.

As given in Appendix A, the chance probability $Q_{\rm HXD}$, for a value of $Z_{4}^{2}\geq 40.64$ to arise via statistical fluctuations, is estimated as $Q_{\rm HXD}\approx 4\%$. Since it is considered reasonably low, the phase modulation is likely to be real, rather than due to statistical fluctuations. The selection of $m=4$ is examined and justified in Appendix B.

For reference, we tried expanding the search range of $T$ up to 200 ks, considering the long data length (242 ks) of NuSTAR. However, no additional DeMD peaks were found.

As seen so far, the XIS and HXD data suggest that the HXC pulses of SGR 1900+14 on this occasion were phase-modulated with $T=40-44$ ks and $A\approx 1.2$ s. However, several inconsistencies and puzzles still remain within the HXD data, and between the XIS and HXD data. They are;

1. S1:

   When the lower energy bound of the HXD data is lowered from 16 keV, the $T\approx 43.5$ ks DeMD peak gradually diminishes, down to $Z_{4}^{2}=33.09$ in 12–50 keV.

2. S2:

   The $T=43.5$ ks peak of the 16–50 keV HXD data is accompanied by $\psi=216^{\circ}$ (Table 1), but it changes to $\psi\approx 60^{\circ}$ if using, e.g., the 12–18 keV band instead. The latter is not consistent, either, with that from the XIS ($\psi\approx 160^{\circ}$).

3. S3:

   The values of $T\approx 40.0$ ks indicated by the XIS (Equation 5) and $T\approx 43.5$ ks by the HXD (Equation 6) appear somewhat discrepant.

4. S4:

   As in Fig. 6 (b), the optimum pulse period from the demodulated HXD data falls on the lowest end of the conservatively estimated error range of $P_{\rm XIS}$.

Among the above issues, [S1] must be taken most seriously, because it is specific to the HXD data, and hence is free from the limited XIS time resolution. It is on one hand consistent with the pulse non-detection (using $m=2$) in Fig. 3 (b). On the other hand, it is puzzling, because expanding the energy range from 16–50 keV to 12–50 keV should normally increase $Z_{4}^{2}$ (see an argument in the next subsection). Therefore, we infer that the pulse coherence degrades when a broader energy range is used, as seen in Paper I. We return to this issue after analyzing the NuSTAR data.

We applied the same demodulation analysis to the NuSTAR data, and obtained the DeMDs in panels (b1) and (b2) of Fig. 5, together with the periodograms in Fig. 6 (c). The lower energy bound of 6.0 keV is chosen to approximately coincide with the SXC vs. HXC cross-over energy (Tamba et al., 2019), and is made somewhat lower than that for the XIS (8 keV), because the NuSTAR effective area at these energies decreases towards lower energies, whereas that of the XIS behaves in the opposite sense. The 6.0–20 keV DeMD reveals a strong peak with $Z_{4}^{2}=66.51$ ($\delta Z_{4}^{2}=29.35$) at

The error is about half that in Equation (6), reflecting the NuSTAR data length which is about twice longer than that with Suzaku. The corresponding DeMD peak is also recognized in the 20–60 keV result at a consistent $T$, although it is considerably less conspicuous, and the amplitude, $A\approx 1.2$ s, is somewhat larger than that in 6–20 keV, $A\approx 0.7$ s.

The reality of the above results was examined in several ways. First, like in the HXD case, we applied the same analysis to the 1st and 2nd halves of the 6–20 keV NuSTAR data, this time allowing $T$ also to vary. Then, the two halves yielded fully consistent DeMDs, in terms of $T$, $A$, and $\psi$. Next, like in § 3.3, we changed $R_{\rm acc}$, to find that the DeMD peak again becomes highest at $R_{\rm acc}\approx 50^{\prime\prime}$, while the modulation parameters ($T$, $A$, and $\psi$) depend little on $R_{\rm acc}$, at least from $35^{\prime\prime}$ to $100^{\prime\prime}$. Therefore, the result is not likely an artifact caused by the stray light from GRS 1919+105. Finally, as given in Appendix A, the peak with $Z_{4}^{2}=66.51$ has a post-trial chance probability of $Q_{\rm NuS}\approx 1\%$, which is even lower than $Q_{\rm HXD}$. We hence conclude that the phase modulation in the 6–20 keV NuSTAR data is real.

From these DeMDs, we presume that the Suzaku and NuSTAR data recorded the same phenomenon. The value of $T$ indicated by NuSTAR is in fact consistent with that of the XIS ($\approx 40.0$ ks; Equation 5), However, $T$ could be inconsistent between the HXD ($\approx 43.5$ ks; Equation 6) and NuSTAR ($\approx 40.5$ ks). If so, the problem [S3] in § 4.2.3 is unlikely to be an artifact due to the limited time resolution of the XIS, and must be regarded as inherent to the HXD data.

Even ignoring for the moment this HXD vs. XIS (plus NuSTAR) discrepancy in $T$, the NuSTAR data themselves are subject to the following two puzzles.

1. N1:

   When the two NuSTAR energy ranges used in Fig. 5, 6–20 and 20–60 keV, are added together, the 40 ks DeMD peak decreases to $Z_{4}^{2}=44.97$, which is higher than that in 20–60 keV (33.47) but much lower than in 6–20 keV (66.51).

2. N2:

   As in Table 1, we find $\psi=3^{\circ}$ and $\psi=246^{\circ}$, in 6–20 and 20–60 keV, respectively. The latter, calculated for $T=38.7$ ks, becomes $\psi\approx 160^{\circ}$ if using $T=40.5$ ks. Therefore, the initial phase is $\sim 180^{\circ}$ off between the two energy ranges.

Evidently, [N1] and [N2] are of the same nature as [S1] and [S2], respectively. In particular, [N1] (as well as [S1]) is puzzling, because periodic signals with a constant pulse profile and insignificant background should satisfy a relation as (Paper I)

where $N_{\rm tot}$ is the total number of signal photons, and PF is the pulsed fraction.

Substituting $\psi(E)$ into Equation (10) and setting $S(E)=0$, we applied the EDMP1 correction to the 12–50 keV HXD data. Starting from initial guesses of $E_{1}=14.0$ keV, $E_{2}=16.0$ keV, and $\Delta\psi=180^{\circ}$, as suggested by [S2] and Fig. 10 (a1,a2), we trimmed these parameters (as well as $P$, $A$, and $\psi_{0}$), so as to maximize the DeMD peak. Then, $Z_{4}^{2}$ has increased (Table 1), and yielded the parameters as in Table 2. In contrast, assuming $\Delta\psi\sim-180^{\circ}$ did not increase $Z_{4}^{2}$. Therefore, like in the NuSTAR case (Fig. 10 b), $\psi(E)$ is thought to inrease by $\sim 180^{\circ}$ from 13.5 keV to 15.8 keV.

The 12–50 keV HXD data were further analyzed via the EDPV2 scheme, by activating $S(E)$. We optimized $R$ and $E_{3}$, as well as the EDPV1 parameters which changed to some extent. The EDPV1/2 parameters determined in this way are also given inTable 2. The DeMD peak further increased, and $R$ turned out to be negative (see a later discussion in § 5.3.3).

Figure 11 (a) superposes DeMDs from the 12–50 keV HXD data, derived under three conditions; via the simple demodulation (dashed black), the EDPV1 modeling (orange), and the EDPV2 scheme (red). Their basic features are summarized in Table 1. The progressively more complex timing corrections have produced the following four noticeable effects.

1. 1.

   The DeMD peak became higher, from $Z_{4}^{2}=33.09$ to $Z_{4}^{2}=41.57$ (EDPV1), and further to $Z_{4}^{2}=47.98$ (EDPV2). Thus, [S1] was mostly solved.

2. 2.

   From the 12–50 keV HXD data, the EDPV1 correction deduced $\psi_{0}=153^{\circ}$, which agrees with $\psi=160^{\circ}$ (Table 1) of the 8-12 keV XIS data. The EDPV2 scheme further enhanced the agreement. This means that [S2] was mitigated.

3. 3.

   The peak centroid of the HXD DeMD evolved from $T\approx 43.5$ ks to $T\approx 42$ ks (EDPV1), and finally to $T\sim 41$ ks (EDPV2) which is fully consistent with those from the XIS and NuSTAR. Therefore, [S3] was solved rather unexpectedly, although its mechanism is unclear.

4. 4.

   The best pulse period with the HXD was at first described by Equation (7), but it has increased to 5.21003 s, which agrees well with $P_{\rm XIS}$. Therefore, [S4] was also solved automatically.

To elucidate the item 4 above, Fig. 12 (a) compares four pulse periodograms, all from the 12–50 keV HXD data but derived in four different ways as explained in the caption. It visualizes that the series of timing corrections have not only increased the pulse significance, but also brought the HXD pulse period in a full agreement with $P_{\rm XIS}$, thus solving [S4].

We applied the same EDPV1 modeling to the 6–60 keV NuSTAR data. Figure 10 (c) depicts how the pulse significance varied with $\Delta\psi$, when $T=40.5$ ks is fixed but all the other parameters ($P$, $A$, $\psi_{0}$, $E_{1}$, and $E_{2}$) are allowed to vary. The data give a clear constraint as $\Delta\psi=162^{\circ}\pm 15^{\circ}$, where the error has been determined in the same way as before. The result also reveals a pair of subsidiary peaks which are $\pm 180^{\circ}$ off the central peak, but they are lower by $\approx 5$ in heights, so that their occurrence probability is each $\sim 1/6$ of that of the central peak. This difference can be explained in the following way: if we select $\Delta\psi$ from either side peak, the coherence in $\Psi$ between the $E\leq E_{1}$ and $E\geq E_{2}$ regions becomes the same as the case with $\Delta\psi=162^{\circ}$, but the coherence in $\Psi$ must be lost between $E_{1}$ and $E_{2}$, causing a decrease in $Z_{4}^{2}$. We hence adopt the central peak, in agreement with Fig. 10 (b) and our conclusion from the HXD data. Including this $\Delta\psi$, the EDPV1 parameters from NuSTAR are given in Table 2.

We analyzed the 6–60 NuSTAR data further employing the EDPV2 corrections, and obtained the optimum parameters as given also in Table 2. Compared with the HXD results, $\Delta\psi$ and $E_{1}$ can be regarded as the same within errors, whereas $E_{2}$ and $E_{3}$ are higher. i Unlike the HXD case, $R$ became marginally positive, and its absolute value is considerably smaller. Therefore, some of the EDPV2 parameters are considered to change with time, just as $A$ varies on time scales of months to years (Makishima et al., 2019).

Figure 11 (b) presents the DeMDs from the 6–60 keV NuSTAR  data, calculated under the same three conditions as for the HXD data. Their basic properties are again given in Table 1. The DeMD peak has become significantly higher, from $Z_{4}^{2}=44.97$ with the simple demodulation (dashed black) to $Z_{4}^{2}=59.04$ (EDPV1; light green), and further to $Z_{4}^{2}=67.40$ (EDPV2; blue), though still smaller than that ($Z_{4}^{2}=66.51$) derived in 6–20 keV via the simple demodulation. The final value of $\psi_{0}=6^{\circ}$ agrees with $\psi=3^{\circ}$ obtained originally in 6–20 keV. In the course of these improvements, the best values of $T$, $A$, and $P$ have remained relatively unchanged. Thus, the EDPV2 modeling on the NuSTAR data has solved [N2], and at least partially [N1] as well.

So far, we used the start energy of 8 keV in the $S(E)$ formalism (equation 12), but this is somewhat arbitrarily. We hence tried changing it to 6 keV or 10 keV, to find that the case with 8 keV is slightly preferred by the NuSTAR data. (The 12–50 keV HXD data are almost insensitive.) This justifies our use of 8 keV as the start energy, both for Suzaku and NuSTAR.

Figure 12 (b1) shows the folded pulse profiles from the HXD data, processed through the EDPV1 scheme. Compared to the results with the simple demodulation in Fig. 7 (b2), the 16–50 keV pulse amplitude became smaller for some reasons, but the profiles became less energy dependent. They exhibit a relatively symmetric pair of horn-like peaks at $\Phi/2\pi\approx\pm 0.3$, which are similar to those found in the XIS profile in 8–12 keV (and to a lesser extent in 6–10 keV). The HXD profiles also show another pair of weaker peaks at $\Phi/2\pi\approx-0.1$ and $\approx+0.5$, and the four peaks are spaced by $\approx 1/4$ cycles. Furthermore, the 12–18 keV profile is seen to lag behind that in 16–50 keV. This “soft lag” presumably demanded the negative value of $R$. (Due to the insufficient XIS time resolution, we cannot draw any conclusion about XIS vs. HXD time lags.)

Figure 12 (b2) shows the HXD profiles corrected for the EDPV2 effect, using the parameters in Table 2 which dictate that the HXD pulse phase is most advanced at $E=(8+E_{3})/2=25.5$ keV, by $\Delta\Phi=R\,(E_{3}-8)/4=-38^{\circ}$ ($-0.11$ pulse cycles). As a result, the profiles have become mostly free from the soft lag, and exhibit the four-peak structure more clearly than before. As an assuring fact, such a four-peak profile was observed from SGR 1900+14 during its Giant Flare in 1988 (Feroci et al., 2001), and from 4U 0142+61 when demodulated (Makishima et al., 2014). On the other hand, this raises a suspicion that we might be biased towards particular EDPV2 solutions that selectively enhance the $m=4$ power. So, in Appendix B, we repeated the analysis by changing $m$, and removed this concern.

The NuSTAR pulse profiles, originally in Fig. 8, became as in Fig. 12 (c), after processed through the EDPV2 modeling. (We skip showing the EDPV1 results, because they are rather similar.) The profiles (in particular that in 10–17 keV) still depend on the energy, but they commonly exhibit a narrow peak at $\Phi/2\pi\approx-0.1$ which are nearly in phase across the entire entire energy. On both sides of this main peak, we observe secondary peaks at a separation by 1/3 to 1/4 cycles, just like in the demodulated HXD profiles (see also Appendix B).

The increase in the pulse significance, from the simple demodulation to the EDPV1 modeling, is $\delta Z_{4}^{2}=8.48$ and $14.07$, with Suzaku  (12–50 keV) and NuSTAR (6–60 keV), respectively. The implied decrease in $Q$ is a factor of $2.7\times 10^{-2}$ (Suzaku) and $1.9\times 10^{-3}$ (NuSTAR). Since these values are not too small, we may not readily exclude the chance origin of these improvements, when the increase in the parameters ($\Delta\psi$, $E_{1}$, and $E_{2}$) are considered. Unlike the case of the simple demodulation, it is not easy, either, to conduct any simulation studies, using the actual or Monte-Carlo-simulated data.

In spite of these limitations, we regard the EDPV1 improvements (Suzaku and NuSTAR altogether) as real, for several reasons. First, the energy-dependent changes in $\psi$ are directly visible from the data (§ 5.1; Fig. 10). Second, the EDPV1 corrections have not only increased $Z_{4}^{2}$, but also solved (at least partially) the puzzles [S1]-[S4], [N1], and [N2]. Third, the Suzaku and NuSTAR solutions agree within errors on $\Delta\psi$ and $E_{1}$, even though they differ in $E_{2}$. Such a coincidence would not easily happen if the pulse enhancements were simply due to chance fluctuations. Finally, a very similar phenomenon has already been confirmed in 1E 1547.0$-$5408 (Paper I) at broadly similar energies, from 10 to 30 keV.

The case for the EDPV2 corrections, where we further incorporate $S(E)$, is more subtle. Probably it is not very significant for NuSTAR, because $\Delta Z_{4}^{2}=5.66$ from EDPV1 to EDPV2 is rather small, in agreement with the fact that $R=0$ is marginally excluded. Actually, the apparent soft lag seen in Fig. 8, between the 6–10 and 10–17 keV profiles, was mostly rectified (though not shown) by the EDPV1 corrections, and a minor hard lag remained, which demanded $R>0$ in EDPV2. In contrast, on the HXD data, the EDPV2 scheme (with $\Delta Z_{4}^{2}=6.41$ over the previous step) is considered more effective, because it removed the 12–18 keV versus 16–15 keV soft lag (Fig. 12) which was left over by the EDPV1 step. This agree with the positive value of $R$ in the EDPV2 solution. Moreover, the 12–50 keV profile has become sharper by the EDPV2 step.

As developed through a series of our studies, the phase modulation of the HXC pulses, which is now confirmed in SGR 1900+14, can be described using the basic dynamics of an axisymmetric rigid body. Namely, we identify $T=40.5$ ks with the slip period of the NS in SGR 1900+14, which is axially elongated by $\epsilon\approx P/T=1.3\times 10^{-4}$ and performs free precession. The deformation can be ascribed to the stress by extreme $B_{\rm t}$ which reaches $\sim 10^{16}$ G.

To be more concrete, consider two Cartesian frames with a common origin (Fig.18 of Makisima et al. 2021); an inertial frame $\Sigma=(\hat{X},\hat{Y},\hat{Z})$ with the unit vector $\hat{Z}$ parallel to ${\bf L}$, and $\Sigma^{*}=(\hat{x}_{1},\hat{x}_{2},\hat{x}_{3})$ fixed to the NS. As before, $\hat{x}_{3}$ is identified with the NS’s symmetry axis. The triplet $(\Phi,\alpha,\Psi)$ provides the three Euler angles, which specify the instantaneous attitude of $\Sigma^{*}$ relative to $\Sigma$. While $\alpha$ is constant, $\Phi$ and $\Psi$ both vary. Every time $\Phi$ completes its one cycle (in the period $P_{\rm pr}$, or one pulse), $\Psi$ advances by $2\pi\epsilon\cos\alpha=2\pi(P_{\rm pr}/T)$ (§ 1). Assuming $\cos\alpha\approx 1$, $\Psi$ returns to its initial value in the slip (or beat) period $T$, which comprises $T/P_{\rm pr}$ precession cycles and $(T/P_{\rm pr})+1$ rotations around $\hat{x}_{3}$.

We further assume that the X-ray emission pattern at each energy is constant when described in the $\Sigma^{*}$ coordinates, and the changes of $\Sigma^{*}$ relative to $\Sigma$ are responsible for all observed variability, including the pulsation and its phase modulation. When the emission pattern is symmetric around $\hat{x}_{3}$, and hence independent of $\Psi$, the pulsation will be strictly periodic, like what we observe for the SXC. If instead the emission breaks symmetry around $\hat{x}_{3}$ (and hence in $\Psi$), the pulse-peak phase becomes dependent on $\Psi$, as modeled by Equation (4) in the simplest case, and actually exhibited by the HXC from SGR 1900+14.

In short, the pulses from a rotating NS become phase-modulated when the following three symmetry breakings all take place; (i) $\alpha\neq 0$, (ii) $\epsilon\neq 0$, and (iii) an asymmetric emission pattern around $\hat{x}_{3}$. The HXC of the relevant magnetars is thought to satisfy all these conditions, whereas the SXC only (i) and (ii). These clear distinctions between the SXC and HXC, in their timing behavior and their spectral shapes, suggest a fundamental difference in their origins.

The above general conditions can be satisfied by a specific geometry given in Paper I, which also affords an explanation of the EDPV effects. Suppose that the SXC is emitted by a region symmetric around $\hat{x}_{z}$, whereas the HXC, arising via, e.g., the two-photon process, has a conical beam pattern around $\hat{x}_{z}$. The cone brightness is assumed to vary with $\Psi$ (the cone azimuth) due, e.g., to the presence of local multipoles which break the symmetry around $\hat{x}_{z}$. Then, the pulse-phase modulation up to $A\sim P/4$ can be explained in a semi-quatitative way (Paper I). Furthermore, let the directional vector $\hat{\xi}$ represent the generatrix along which the HXC is brightest. If $\hat{\xi}$ moves in $\Psi$ with energy, due to some strong-field physics such as proton cyclotron resonances (Paper I), the EDPV1 effects can be explained.

As mentioned in Paper I, the EDPV2 effect represented by $S(E)$ is more difficult to interpret. To see this, let $\Pi_{3}$ be the plane defined by $\hat{Z}\|{\bf L}$ and $\hat{x}_{3}$, which rotates around $\hat{Z}$ with a period $P_{\rm pr}$. Then, $\hat{\xi}$ as defined above, rotates relative to $\Pi_{3}$ with a period $T$, in which it crosses $\Pi_{3}$ twice. At every crossing, the pulse arrival-time delay $\delta t$ will change its sign. When averaged over a cycle, we expect $\langle\delta t\rangle\approx 0$, or $S(E)\approx 0$, as long as Equation (4) provides a good approximation. However, if the time profile of $\delta t$ is much deviated from a sinusoid, and asymmetric between $\delta t>0$ and $\delta t<0$, we may expect $S(E)\neq 0$. Thus, we tentatively regard $S(E)$ as a modeling artifact, which would vanish when we improve Equation (4). Such attempts were already made in Paper I, but only very preliminarily.

So far, we have adopted the interpretation of mangetars as isolated NS powered by magnetic energies (Mereghetti, 2008). Some alternative models however describe them as NSs with ordinary dipole magnetic fields as $B_{\rm d}\sim 10^{12}$ G, powered by mass accretion from fossil or fallback disks around them (e.g. Benli and Ertan, 2016). In fact, infrared observations provided evidence for such disks around some magnetars, including 4U 0142+61 in particular (e.g. Wang et al., 2006).

Based on the disk-accretion scenario, Grimani (2021) argued that the 55 ks modulation in 4U 0142+61 can be explained when the source is hidden periodically by the disk if it is in a Keplerian rotation and is free-precessing. However, the disk is so distant (Appendix C) that the emission region would look point-like, and the disk is not highly ionized. Then, the modulation must get stronger towards lower energies because of the increasing photo-absorption, contradicting to the general absence of phase modulation in the SXC pulses. Therefore, this scenario would work for neither 4U 0142+61, nor SGR 1900+14 which is in a similar condition. Some other mechanisms must be sought for if the observed phenomenon it to be explained by the disk-accretion scenario.

Regardless of details of the phase-shift production, the EDPV effects in the two objects must be explained. Taking SGR 1900+14 for example, the 12 keV pulses at a particular phase in $\Psi$ need to arrive by $\sim 1$ s earlier than expected, whereas the 20 keV pulses later by a similar amount (Fig. 10). Such a sharp energy dependence is incompatible with the broadband nature of X-rays from accretion columns, wherein the $\sim 12$ keV and $\sim 20$ keV photons must behave in positively correlated ways. In contrast, our scenario based on the strong-field physics (§ 6.2.2) can explain the essential timing properties of SGR 1900+14 including its EDPV behavior. Therefore, the X-rays produced via accretion, if any, should contribute little to the overall X-ray emission.

We also examined how a circum-stellar disk affects the NS dynamics by forced precession. This may take place through two channels; one is via the accretion torque, wherein the NS can be spherical but needs to be accreting. The other is via direct gravity; the NS needs to be deformed, but the accretion is not required. As given in Appendix C, the former mechanism predicts a long period of forced precession as $\sim 2\times 10^{3}$ yr, whereas the latter is even longer by many orders of magnitude. Thus, the rigid-body dynamics of these magnetars are not affected by the disks around them.

Although we have argued against the accretion scenario, we do not mean that NSs with $B_{\rm d}\gtrsim 10^{14}$ G cannot become accreting sources, or cannot reside in binaries. In fact, the binary X-ray pulsar X-Persei, accreting from the companion’s stellar winds, has been found to have $B_{\rm d}\sim 1\times 10^{14}$ G (Yatabe et al., 2018). As another interesting case, the gamma-ray binary LS 5039 is likely to harbor a non-accreting magnetar, whose magnetic energy is released via interactions with the primary’s stellar winds to drive the remarkable non-thermal activity (Yoneda et al., 2020, 2021). Yet another example is so-called Central Compact Objects (CCOs), rather inactive NSs found at the center of some supernova remnants (e.g. Esposito et al., 2019). They are thought to have weak $B_{\rm d}$ but intense $B_{\rm t}$, and the latter sustains their activity. Thus, a fair fraction of NSs may be born as magnetars in a broad sense (Nakano et al., 2015), and reside in various environments. As suggested by the present study, some, if not all, of them might have $B_{\rm t}\sim 10^{16}$ G. It would hence be an interesting future work to classify magnetized NSson the $(B_{\rm d},B_{\rm t})$ plane.