DEEP NuSTAR AND SWIFT MONITORING OBSERVATIONS OF THE MAGNETAR 1E 1841−045

3.1.1. Temporal and Spectral Properties of the Bursts

We performed a comprehensive search for bursts in the NuSTAR light curves. We extracted event time series, applied the barycenter correction using the position R.A. = 18^h41^m19343, decl. = −04°56'118 (J2000; Durant 2005), binned the light curves with a variety of bin sizes ranging from 0.01 to 10 s, and searched for time bins that contained more counts than expected above background including source persistent emission using Poisson statistics. The background was extracted in an interval 10 pulse periods long using the same extraction region, just prior to the time bin being considered. In total, we found seven time bins that are significantly above the mean level after considering the number of trials. Note that two of the seven significant bins turned out to be produced by one burst (burst 5; see below); hence, we found six bursts during our observations. The significance of the bursts is high (p < 10⁻¹⁰) but only on short timescales, e.g., ≲0.1 s. We list the burst times in Table 2. Note that bursts 5 and 6 were previously reported based on Fermi GBM detections (Collazzi et al. 2014; Pal'shin et al. 2014).

Table 2.  Deadtime-corrected Properties of NuSTAR-detected Bursts from 1E 1841–045

Burst	T0	ϕa	Tr	Tf	A	C1	C2	T90b	Nevt	Δpost−pre
	(day)		(s)	(s)	(cps)	(cps)	(cps)	(s)	(cts)	(cts)
1	0.35836545	${0.2266}_{-0.0002}^{+0.0003}$	${0.007}_{-0.001}^{+0.001}$	${0.047}_{-0.019}^{+0.017}$	${1740}_{-600}^{+760}$	${12}_{-2}^{+2}$	<16	${0.12}_{-0.05}^{+0.04}$	5	⋯
2	0.35837490	${0.2958}_{-0.0002}^{+0.0002}$	${0.006}_{-0.002}^{+0.002}$	${0.09}_{-0.01}^{+0.01}$	${700}_{-90}^{+100}$	${5.5}_{-1.4}^{+1.6}$	${5.9}_{-1.6}^{+1.8}$	${0.22}_{-0.03}^{+0.03}$	61	75 ± 41
3	0.60981692	${0.5583}_{-0.0001}^{+0.0002}$	${0.014}_{-0.001}^{+0.002}$	${0.011}_{-0.002}^{+0.003}$	${2010}_{-410}^{+460}$	${12}_{-2}^{+2}$	${5}_{-3}^{+3}$	${0.057}_{-0.007}^{+0.011}$	27	16 ± 40
4	7.27821493	${0.4463}_{-0.0002}^{+0.0002}$	${0.0053}_{-0.0014}^{+0.0016}$	${0.052}_{-0.02}^{+0.02}$	${840}_{-240}^{+300}$	${16}_{-2}^{+2}$	<3.4	${0.13}_{-0.05}^{+0.05}$	11	26 ± 40
5	8.62801288	${0.1075}_{-0.0001}^{+0.0001}$	${0.0090}_{-0.0003}^{+0.0003}$	${0.0184}_{-0.0004}^{+0.0004}$	${67000}_{-60000}^{+70000}$	${15}_{-1}^{+1}$	<3	${0.0631}_{-0.0007}^{+0.0007}$	22	211 ± 43
6	8.759684723	${0.83672}_{-0.00001}^{+0.00001}$	<0.0006	${0.0249}_{-0.0018}^{+0.0016}$	${8000}_{-1100}^{+1200}$	${14}_{-1}^{+2}$	<4	<0.059	31	7 ± 41

Notes. Parameters for the short-timescale light curve. T₀ is the burst arrival time and is days since MJD 56,540 (barycentric dynamical time). T_r,f are the rising and falling times for the burst light curves, A is the peak count rate, C_1,2 are the constant levels of the light curves before and after the bursts (see Equation (1)), T₉₀ is the time interval that includes 90% of the burst counts estimated with the exponential functions, N_evt is the number of events within T₉₀, and Δ_post−pre is the difference in numbers of photons contained in the pre- and the post-burst 200 s intervals.

^aSpin phase corresponding to T₀, where phase zero is defined at the pulse minimum (T_ref = 56,540.32899020 MJD), the same as that for the timing analysis in Section 3.2.2. ^bSince T₉₀ for the whole burst is not always well defined because the constants C_1,2 are different before and after the burst peak, T₉₀ values for the rising and the falling function were calculated separately and then summed to obtain that for the burst. When only an upper limit is available for T_r or T_f, we used the upper limit to calculate T₉₀ and show it without uncertainties.

Download table as:  ASCIITypeset image

We note that some of the bursts may not be fully sampled owing to high count rates. Specifically, the maximum count rate of NuSTAR detectors is ∼400 cps limited by a deadtime of 2.5 ms per event. In addition, an event that was detected in a very short elapsed time from the previous event is regarded as background and filtered out during the standard pipeline process. We investigated these effects by looking into the elapsed time for each event and found that those for events in two high-count time bins were significantly shorter than in other time intervals. We further reprocessed the observation data with a relaxed elapsed-time filter (see Madsen et al. 2015, for more details) and were able to recover an additional 58 events in an R = 120'' circular aperture in the 3–79 keV band for the two time bins combined. From this study, we find that the two high-count bins were actually connected in time (i.e., the gap between the two bins was filled by the recovered events) and the livetime of the detectors was ∼1/300 of the exposure. We also investigated other high-count time bins using a relaxed filter and were able to recover marginally additional events only for burst 6 (10 ± 8 events), the other GBM-detected burst. Below, results for burst properties are obtained with the data processed with the relaxed filter.

Although it was reported that the Fermi-detected bursts are likely from the magnetar 1E 1841–045 (Pal'shin et al. 2014), the localization was not unambiguous. In order to localize the bursts better and see whether they are really produced by the magnetar 1E 1841–045, we used the NuSTAR data to measure the position of bursts 2, 3, 5, and 6, which had sufficiently many counts for such an analysis. We projected their 3–79 keV event distributions collected for 2 s onto R.A. and decl. axes and fit the projected profiles with a Gaussian plus constant function. The results for the burst location offsets from the 1E 1841–045 position were ΔR = 9'' ± 2'', 8'' ± 2'', 9'' ± 3'', and 5'' ± 2'' for bursts 2, 3, 5, and 6, respectively, where the quoted uncertainties are purely statistical (1σ). Note that aspect reconstruction accuracy of NuSTAR is 8'' (90%), and so the measured positions are consistent with that of the magnetar 1E 1841–045.

In order to characterize the properties of the bursts, we fit the short-term (∼10 s) light curves around the bursts with exponentially rising and falling functions

$$\begin{eqnarray}F(t)=\left\{\begin{array}{ll}{{Ae}}^{(t-{T}_{0})/{T}_{{\rm{r}}}}+{C}_{1} & t\lt {T}_{0},\\ (A-{C}_{2}){e}^{-(t-{T}_{0})/{T}_{{\rm{f}}}}+{C}_{1}+{C}_{2} & t\geqslant {T}_{0},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where A is the amplitude, T₀ is the burst peak time, T_r is the rising time, T_f is the falling time, and C_1,2 are constants (see Gavriil et al. 2011, for a different model). Since there are only 10–120 counts detected within a 1 s window around each burst, we extracted events from the whole detector and used a maximum likelihood optimization without binning. The peak count rates are very high, and the T₉₀ durations are ≲1 s. We present the results in Table 2 and show the bursts morphologies in Figure 1. Some bursts occurred within 1 s of each other. Note that we do not find a clear increase in the tail emission except for in burst 5.

Zoom In Zoom Out Reset image size

The spectrum of a burst can provide information on the burst mechanism. Therefore, we extracted ∼0.2 s spectra in circles with radius R = 120'' centered at the source position for bursts 2 and 5 and fit the spectra with single-component models; we tried both a blackbody and a power-law model. In order to remove the persistent emission, we extracted background counts in the same region as we used for the source but in a 1 ks pre-burst interval. We did not attempt to fit the burst data with a multicomponent model because there were too few photons collected during the 0.2 s intervals. We grouped the spectrum to have 1 count per spectral bin and used lstat (Loredo 1992) in XSPEC 12.8.1 g. The NuSTAR bandpass is not sensitive to the relatively low N_H of the source, and we set it to 2.05 × 10²² cm⁻², which we obtained by joint fitting of the soft-band spectra with a broken power-law model (see Section 3.2.3). We use this value and the tbabs absorption model in XSPEC throughout this paper unless noted otherwise. The burst spectra can be described with a power-law model having Γ = 1–2 or a blackbody model with kT = 3–5 keV. The results of the fits are shown in Table 3. Note that we did not detect the high-temperature blackbody component (kT = 13 keV; Collazzi et al. 2014), probably because of the low statistics at high energy.

Table 3.  Best-fit Parameters of the Burst Spectra for a ∼0.2 s Interval around the Burst Peak

Burst	Γ	3–79 keV Fluxa	lstat/dof
		(10−8 erg s−1 cm−2)
2	1.6(3)	3.5(1.2)	37/46
5	0.96(24)	1300(600)	53/56
Burst	kT	LBBb	lstat/dof
	(keV)	1038 erg s−1
2	3.3(3)	1.5(3)	48/46
5	4.8(5)	400(100)	59/56

Notes.

^aDeadtime-corrected flux. ^bDeadtime-corrected bolometric luminosity for an assumed distance of 8.5 kpc.

Download table as:  ASCIITypeset image

3.1.2. Burst Tails

Short X-ray bursts from magnetars can exhibit long emission tails, lasting for hours (An et al. 2014). In order to search for long tails, we binned the light curves into 20 s bins and compared counts in 10 bins before and after a burst, excluding the burst bin. We then compared the pre- and post-burst counts and found that the difference is significant only for burst 5 (Δ_post−pre = 211 ± 43 counts for a 200 s time interval). We show the light curves for the bursts in Figure 2 and report the results in Table 2. We performed the same study on different timescales (e.g., 2 s and 50 s) and found the same results; the difference is significant only for burst 5.

Zoom In Zoom Out Reset image size

In order to search for spectral evolution after burst 5, we extracted spectra within a radius R = 120'' in the tail of burst 5 in several time intervals excluding the burst (T > T₀ + 0.5 s; see Table 2). Each interval had ≳100 events above the persistent emission plus background. We then fit the spectra with a blackbody and a power-law model. We show the results of the power-law fit in Figure 3. The spectral shape did not change significantly over ∼2 ks of tail emission, and the 3–79 keV flux decay is well described with a power-law function, having a decay index of 0.45 ± 0.10. Note that the decay index was measured after taking into account the covariance between the flux and the spectral index in the fit. A similar decay index is measured for the luminosity when using a blackbody model. Note that An et al. (2014) measured the flux evolution index in the tail to be 0.8–0.9 for 1E 1048.1–5937, much steeper than our measurement for 1E 1841–045. The count enhancement at later times seems to be spiky and might have been caused by undetected activity there. Since the other bursts do not have significant tail emission, we were not able to measure their spectral evolution.

Zoom In Zoom Out Reset image size

The tail spectrum of burst 5 is soft compared to the burst spectrum and may be similar to the persistent emission. In order to see whether the spectra of the burst tails are significantly different from the persistent source emission, we studied the persistent emission in a pre-burst interval in which approximately the same number of events was collected as in the tail spectra. We extracted a persistent spectrum in a 100 s pre-burst time interval and modeled it with a single-component model, either a power law or a blackbody. The persistent spectrum over the short interval was well described with single-component models, and we find that the best-fit parameters are Γ = 2.31 ± 0.24 or kT = 1.58 ± 0.13 keV, similar to the tail spectrum. We show the persistent levels with blue dashed lines in Figure 4.

Zoom In Zoom Out Reset image size

Since the spectral shape did not change significantly over the tail interval, we fit the combined (∼1.8 ks) tail spectrum to a power-law or a blackbody model after removing the persistent emission. The spectrum is well fit by a power-law model (χ²/degrees of freedom (dof) = 122/146) with a photon index Γ = 2.2 ± 0.2 (green lines in left panel of Figure 3), similar to the 100 s persistent spectrum. A blackbody model also fits the data with kT = 2.1 ± 0.2 keV (χ²/dof = 137/146). We tried to fit the combined tail emission with the model for the persistent spectra obtained below (Section 3.2.3) using the same fit parameters except for the normalization constants (Table 4). The model fits the spectrum well (χ²/dof = 171/147 with the null hypothesis probability p = 0.09). However, we see a trend in the fit residuals, which suggests that the tail spectra may be slightly softer than the persistent one (see Figure 3, right panel).

Table 4.  Phenomenological Spectral Fit Results for 1E 1841–045

Phase	Dataa	Energy	Modelb	NH	kT	Γsc	Ebreak/Fsd	Γh/βe	Fhf	LBBg	χ2/dof
		(keV)		(1022 cm−2)	(keV)		(keV/ )
0.0–1.0	N, S, X, C	0.5–79	BB+BP	2.05(3)	0.491(5)	1.95(1)	13.5(2)/⋯	1.24(1)	5.88(6)	1.64(5)	8060/7930
0.0–1.0	N, S, X, C	0.5–79	BB+2PL	2.49(5)	0.443(9)	2.82(8)	⋯/1.53(6)	0.97(4)	4.70(6)	1.15(9)	7931/7930
Pulsed	N, X, C	0.5–79	PL	2.05h	⋯	⋯	⋯	1.83(3)	1.4(1)	⋯	1236/1749
Pulsed	N	3–79	PL	2.05	⋯	⋯	⋯	1.81(3)	1.4(1)	⋯	1128/1591
Pulsed	N	5–79	PL	2.05	⋯	⋯	⋯	1.73(4)	1.4(1)	⋯	806/1140
Pulsed	N	10–79	PL	2.05	⋯	⋯	⋯	1.47(7)	1.5(2)	⋯	349/510
Pulsed	N	15–79	PL	2.05	⋯	⋯	⋯	1.28(12)	1.6(3)	⋯	174/279
Pulsed	N	20–79	PL	2.05	⋯	⋯	⋯	1.0(2)	1.5(4)	⋯	113/159

Notes.

^aN: NuSTAR; S: Swift; X: XMM-Newton; C: Chandra. ^bBB: Blackbody; PL: Power law; BP: Broken power law; and 2PL: Two power laws in XSPEC. ^cPhoton index for the soft power-law component. ^dBreak energy for the BB+BP fit or soft power-law flux in units of 10⁻¹¹ erg s⁻¹ cm⁻² in the 3–79 keV band for the BB+2PL fit. ^ePhoton index for the hard power-law component. ^fFlux in units of 10⁻¹¹ erg cm⁻² s⁻¹. The values are only the power-law (hard power-law) flux in the 3–79 keV band for the BP (PL, 2PL) model. ^gBlackbody luminosity in units of ${10}^{35}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ for an assumed distance of 8.5 kpc (Tian & Leahy 2008). ^hN_H for the pulsed spectral analysis was frozen.

Download table as:  ASCIITypeset image

In order to study the persistent emission, we removed the burst intervals using time filters. We used 20 s windows centered at the burst peak times for all the bursts except for burst 5, for which we used a 2 ks window because of the tail. For the soft-band spectrum below 10 keV, we used the Chandra (only Obs. ID 730 due to pileup in Obs. ID 6732), XMM-Newton, and Swift data (see Table 1), the same observations as used by A13. Although these observations were taken long before the NuSTAR observation, we show below that the source emission properties have been stable over 10 yr (see also A13). For the NuSTAR data, we used a circular aperture with R = 60'' for the source and an annular aperture with inner and outer radii of R = 60'' and R = 100'' for the background, respectively. Note that ∼15% of the source events fall in the background region because of the NuSTAR point-spread function (An et al. 2014). However, not all the ∼15% of the source flux is subtracted as background in the spectral fitting since we scale the area of the background region to that of the source for spectral fits, and ∼10% of source events will be lost during background subtraction. We take into account this effect using a normalization constant.

3.2.1. Timing Analysis

3.2.2. Pulse Profiles and Pulsed Fraction

The pulse profile of 1E 1841–045 is known to change with energy (Kuiper et al. 2006; An et al. 2013). In particular, A13 found that the pulse shape in the 24–35 keV band is different from those in the adjacent energy ranges, which suggested the existence of a spectral feature. However, no firm conclusion could be made owing to limited statistics. We investigate this here with much better statistics.

We produced pulse profiles for individual observations and aligned them with the template. Backgrounds were extracted from an annular region with inner and outer radii of 60'' and 100'' and were subtracted from the pulse profile of the source region. We verified that the pulse shape has not changed significantly over the ∼300 days of NuSTAR observations (Table 1) by comparing the pulse profile of individual observations with the combined profile in several energy bands (e.g., 3–6 keV, 6–10 keV, 10–15 keV, and 15–79 keV). We show the combined pulse profiles in several energy bands in Figure 5. Note that we find double-peaked structure similar to that seen by A13, e.g., in the 17–33 keV band (see Figure 5). However, with the much better statistics we have now, we find that the pulse shape does not change suddenly but instead changes gradually with energy. This does not support the existence of a narrow spectral feature, as suggested by A13.

Zoom In Zoom Out Reset image size

We find that the pulse shape at higher energies becomes more complicated, sometimes showing three peaks (e.g., 33–38 keV in Figure 5; a new peak seems to appear at phase ∼0.5). In order to see whether the triple-peaked structure at high energies (≳25 keV) is significant, we fit each pulse profile with a harmonic function in which the number of harmonics contained was varied between zero (constant) and three. In the fit, we calculate the χ² value and the F-test probability by adding higher-order harmonic functions one by one. From this study, we find that the pulse profiles below 38 keV are generally well fit with the sum of two harmonics, and the others with a single harmonic function; adding one more harmonic to these is not statistically required (99% confidence). We further fit the pulse profiles with a sum of the first two harmonics plus a fifth harmonic because the triple-peaked structure is best described with a fifth harmonic. The F-test probability for adding the fifth harmonic shows an improvement to the fit with 99.7% confidence in the 45–55 keV band and with 98.6% confidence in the 33–38 keV band. However, these may not imply a significant detection when considering the number of double-peaked profiles we have. Therefore, we conclude that the double-peaked structures are statistically significant but the triple-peaked ones are only marginally so.

The pulsed fraction of the source has previously been measured to be increasing with energy (Kuiper et al. 2006; A13). Note that A13 did not attempt to measure the area pulsed fraction defined by

$$\begin{eqnarray}&&{{PF}}_{\mathrm{area}}=\displaystyle \frac{{\displaystyle \sum }_{i}({p}_{i}-{p}_{\mathrm{min}})}{{\displaystyle \sum }_{i}{p}_{i}},\end{eqnarray} \tag{ 2 }$$

where p_i is counts in ith phase bin and p_min is the counts from the lowest bin in the pulse profile, due to insufficient statistics, and that PF_area is known to be biased upward in the low counts regime. Since we now have much better statistics, we measured the area pulsed fractions in different energy bands and show them in Figure 6. Even with good statistics, the area pulsed fraction is known to be biased upward (see the appendix), and so we show two alternative measures of the pulsed fraction as well. The first alternative is to fit the pulse profile with a harmonic function and use the best-fit function to calculate the pulsed fraction (PF_fit), which will remove the bias caused by incorrect identification of the baseline level by selecting the minimum phase bin. We find that the pulse profiles are well described with two harmonics having χ²/dof ∼ 1 (Figure 5) except for the lowest-energy band, for which the two-harmonic function yields an unacceptable fit, with p = 10⁻⁴. The large χ² in this case is mostly due to a sharp step in the pulse profile at phase ∼0.2, but the overall pulse shape agrees with the best-fit function well. Furthermore, rebinning the phase can make the fit acceptable. Therefore, we used two harmonics for the fit function, calculated PF_area using the best-fit parameters, and show the pulsed fractions in Figure 6 (blue squares). We also calculated the rms pulse amplitude defined by

$$\begin{eqnarray}&&{{PF}}_{\mathrm{rms}}=\displaystyle \frac{\sqrt{2{\displaystyle \sum }_{k=1}^{6}(({a}_{k}^{2}+{b}_{k}^{2})-({\sigma }_{{a}_{k}}^{2}+{\sigma }_{{b}_{k}}^{2}))}}{{a}_{0}},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray*}&&{a}_{k}=\displaystyle \frac{1}{N}\underset{i=1}{\overset{N}{\displaystyle \sum }}{p}_{i}\mathrm{cos}(\displaystyle \frac{2\pi {ki}}{N}),{b}_{k}=\displaystyle \frac{1}{N}\underset{i=1}{\overset{N}{\displaystyle \sum }}{p}_{i}\mathrm{sin}(\displaystyle \frac{2\pi {ki}}{N}),\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\sigma }_{{a}_{k}}^{2}+{\sigma }_{{b}_{k}}^{2}=\displaystyle \frac{1}{{N}^{2}}\underset{i=1}{\overset{N}{\displaystyle \sum }}{\sigma }_{{p}_{i}}^{2}{\mathrm{cos}}^{2}(\displaystyle \frac{2\pi {ki}}{N})+\displaystyle \frac{1}{{N}^{2}}\underset{i=1}{\overset{N}{\displaystyle \sum }}{\sigma }_{{p}_{i}}^{2}{\mathrm{sin}}^{2}(\displaystyle \frac{2\pi {ki}}{N})\end{eqnarray*}$$

is the Fourier power produced by the noise in the data, p_i is the number of counts in ith bin, N is the total number of bins, ${\sigma }_{{p}_{i}}$ is uncertainty in p_i, and n is the number of Fourier harmonics included, in this case, n = 6 (see Archibald et al. 2014, for more details). We show the results in Figure 6 (red diamonds).

Zoom In Zoom Out Reset image size

We note that the small-scale features in Figure 6 change when we use a different energy resolution. For example, a sudden jump sometimes appears at ∼30 keV, similar to that seen by A13 (their Figure 3), if we use different energy bins. However, the overall trend is similar; we do not see any rapid increase in the pulsed fraction with energy above 10 keV. This is different from previous reports (Kuiper et al. 2006) of pulsed fraction increasing with energy and approaching 100% at 100 keV.

3.2.3. Phase-averaged Spectral Analysis

Since it is possible that the source has different emission properties during the bursting periods, we compared the spectral properties of individual observations. In this study, we did not include the soft-band spectra because they were taken at much earlier epochs. We jointly fit all the NuSTAR spectra in the 6–79 keV band with an absorbed broken power-law model in order to minimize the effect of the blackbody component, which is negligible above 6 keV, and found that the spectral shapes for the six NuSTAR observations (see Table 1) are consistent with one another.

Since there is no clear evidence that the source spectral shape has varied during the NuSTAR observations, we used all the NuSTAR observations for the phase-averaged spectral analysis. Furthermore, we used the soft-band data as well in this study since the shape of the soft-band spectrum is also known to be stable (Zhu & Kaspi 2010; Dib & Kaspi 2014). We tied all the model parameters between NuSTAR, Swift, XMM-Newton, and Chandra except for the cross-normalization factors. The normalization constant for NuSTAR FPMA (Obs. ID 30001025002) was set to be 0.9 as a reference in order to account for the source contamination in the background region (Section 3.2.2). To fit the data, we grouped the spectra to have at least 20 counts per bin. We used an absorbed blackbody plus broken power law and an absorbed blackbody plus a double power law to fit the 3–79 keV NuSTAR data and the 0.5–10 keV soft-band data. We present the results in Table 4 and the spectra in Figure 7.

Zoom In Zoom Out Reset image size

We note that the spectral parameters we report in Table 4 are slightly different from those reported previously (A13). In order to see whether the difference is due to the updated calibration, we analyzed the same data set that A13 used (Obs. ID 30001025002) and were able to reproduce their results except for the flux. The flux we measure is lower by ∼15% than what A13 reported, because of a ∼15% increase in the NuSTAR effective area from CALDB 20131007.¹⁷ The hard-band component is much better constrained with the new long exposures, and thus the new results we report are more accurate. We note that the new parameters in Table 4 are not inconsistent with the data A13 used, providing acceptable fits to the data with χ²/dof = 2965/2878 and 2914/2878 for the blackbody plus broken power law and the blackbody plus double power law, respectively.

We conducted a phase-resolved spectral analysis for 10 phase intervals to study distinct features in the pulse profiles (see Figure 5 for pulse profiles). We did not use the Swift XRT or XMM-Newton MOS data in this study because their temporal resolutions were insufficient. The Chandra and XMM-Newton pn data were phase-aligned with the NuSTAR data by correlating the pulse profiles.

We binned the NuSTAR and soft-band instrument spectra to have at least 20 counts per spectral bin, and we froze the cross-normalization factors to those obtained with the phase-averaged spectral fit. We fit the spectra with the two models that we used for fitting the phase-averaged spectrum: an absorbed blackbody plus broken power-law and an absorbed blackbody plus double power-law model. We find that both models explain the data well, having χ²/dof ≲ 1.003 for dofs of 1413–1966. The spectra vary with spin phase, having harder power-law spectra when the flux is high. However, the detailed variation depends on the spectral model used. We show the results in Figure 8.

Zoom In Zoom Out Reset image size

In order to see whether there is a spectral feature that shifts with energy, as was seen in SGR 0418+5729 (Tiengo et al. 2013), we produced an energy-phase image using 25 phase bins and 40 energy bins in the 3–79 keV band. We first divided counts in each pixel in the energy-phase 2D map by the phase-integrated counts in the same energy bin, and we present it in the left panel of Figure 9, which shows similar structures to the energy-resolved pulse profiles (Figure 5). We then divided the map further by the energy-integrated counts in the same phase bin (Figure 5) in order to have a better contrast, and we find no clear phase-dependent feature in the image (Figure 5 right). We also tried different binning and found the same results.

Zoom In Zoom Out Reset image size

We measured the pulsed spectrum in the 0.5–79 keV band in order to see whether it is significantly harder than the phase-averaged spectrum as seen in other hard-X-ray bright magnetars (Kuiper et al. 2006). We grouped the spectra to have at least 200 counts per spectral bin and subtracted the spectrum in the phase interval 0.9–1.1 (the DC level) from the phase-averaged spectra obtained in Section 3.2.3. We then jointly fit the broadband spectrum with a power-law model, letting the cross-normalization constants vary. We used lstat and χ² statistics and found that they give consistent results.

A simple power-law model with a photon index of 1.8 fits the 0.5–79 keV data well (reduced χ²/dof < 1). In order to verify the measurements for the pulsed spectrum above 15 keV (Γ = 0.72 ± 0.15; Kuiper et al. 2006), we restricted the fit range to high energies. As the lower-energy cutoff is increased, the power-law index decreases, consistent with spectral hardening. The results are summarized in Table 4.

We also estimated pulsed fractions in the hard band using the spectra (defined as the ratio of pulsed and total flux densities) in order to compare with those in Figure 6. We fit the total and the pulsed spectra (≳15 keV) to single power-law models. The power-law index of the total spectrum above 15 keV is 1.19 ± 0.02 (1.13 ± 0.02 above 20 keV), slightly smaller than what we obtained using the absorbed blackbody plus broken power-law model (see Table 4). The power-law index of the pulsed spectrum above 15 keV is 1.28 ± 0.12 (1.0 ± 0.2 above 20 keV), similar to that of the total spectrum. This suggests that the pulsed flux does not rapidly increase in the hard band, as also seen in Figure 6.

A13 found that the properties of the persistent hard-X-ray emission of 1E 1841–045 were consistent with the coronal outflow model proposed by Beloborodov (2013). The model envisions an outflow of relativistic electron-positron pairs created by electric discharge near the neutron star. The outflow fills the active "j-bundle"—a bundle of closed magnetospheric field lines that carry electric current (Beloborodov 2009). The pair plasma flows out along the magnetic field lines and gradually decelerates as it scatters the thermal X-rays. It radiates most of its kinetic energy in hard-X-rays before the e^± pairs reach the top of the magnetic loop and annihilate. The magnetic dipole moment of 1E 1841–045 is estimated from its spin-down rate, μ ≈ 7 × 10³² G cm³, assuming the neutron star radius to be 10 km. Similar to A13 and Hascoët et al. (2014), we assume a simple geometry where the j-bundle is axisymmetric around the magnetic dipole axis. However, instead of assuming that the j-bundle emerges from a polar cap, its footprint is allowed to have a ring shape. This possibility was introduced in Vogel et al. (2014), because the NuSTAR data for 1E 2259+586 favored a ring footprint over a polar cap. The model has the following parameters: (1) the power ${L}_{j}$ of the e^± outflow along the j-bundle, (2) the angle ${\alpha }_{\mathrm{mag}}$ between the rotation axis and the magnetic axis, (3) the angle ${\beta }_{\mathrm{obs}}$ between the rotation axis and the observer's line of sight, (4) the angular position ${\theta }_{j}$ of the j-bundle footprint, and (5) the angular width ${\rm{\Delta }}{\theta }_{j}$ of the j-bundle footprint. In addition, the reference point of the rotational phase, ${\phi }_{0}$, is a free parameter, since we fit the phase-resolved spectra.

We follow the method presented in Hascoët et al. (2014) and explore the whole parameter space by fitting the phase-averaged spectrum of the total emission (pulsed+unpulsed) and phase-resolved spectra of the pulsed emission. We use five equally spaced phase intervals. NuSTAR data are fitted above 10 keV, where the coronal outflow has to account for most of the observed emission.

Figure 10 shows the map of the p-value of the fit in the ${\alpha }_{\mathrm{mag}}$-${\beta }_{\mathrm{obs}}$ plane. The acceptable model is clearly identified in this map.¹⁸ It has ${\alpha }_{\mathrm{mag}}=0.25\pm 0.15,{\beta }_{\mathrm{obs}}=1.0\pm 0.2$, ${\theta }_{j}=0.24\pm 0.02$, and ${\rm{\Delta }}{\theta }_{j}/{\theta }_{j}\gt 0.26$, consistent with a polar cap. The corresponding magnetic flux in the j-bundle is (2.5 ± 0.4) × 10²⁶ G cm². The power dissipated in the j-bundle is ${L}_{j}=(6\pm 1)\times {10}^{36}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$. Most of the released energy is radiated in the MeV band (peaking at ∼6 MeV) and is not seen to NuSTAR. Using the obtained best-fit model for the hard-X-ray component, we have investigated the remaining soft-X-ray component. The procedure is similar to that in A13 and Hascoët et al. (2014); we freeze the best-fit parameters of the outflow model and fit the spectrum in the 0.5–79 keV band using the NuSTAR, Swift, Chandra, and XMM-Newton data. As in A13, we find that the 0.5–79 keV spectrum is well fitted by the sum of two blackbodies (which dominate below 10 keV) and the coronal outflow emission (which dominates above 10 keV). The cold and hot blackbodies have luminosities L_c = (2.2 ± 0.1) × 10³⁵ erg s⁻¹, L_h = (9.8 ± 1.3) × 10³⁴ erg s⁻¹ and temperatures kT_c = 0.45 ± 0.01 keV, ${{kT}}_{{\rm{h}}}=0.75$ ± 0.02 keV. Note that these values are different from those we obtained with the phenomenological models in Table 4.

Zoom In Zoom Out Reset image size

To investigate the spectrum of 1E 1841–045 in the monitoring observations, we extract spectra for each observation using a 10-pixel (24'') long strip centered on the source. An annulus of inner radius 75 pixels and outer radius 125 pixels centered on the source was used to extract background spectra. The spectra were grouped to have a minimum of 20 counts per bin. The spectra were fitted with a photoelectrically absorbed blackbody with an added power-law component, using the tbabs(bbody+pow) model in XSPEC 12.8.1, with N_H held fixed at 2.05 × 10²² cm⁻², the value we obtained in Section 3.2.3. No significant change in source 1–10 keV flux was observed over the monitoring period of ∼3 yr (χ²/dof = 49/59), including the NuSTAR-observed bursting period. The result is shown in Figure 11(a).

Zoom In Zoom Out Reset image size

We also searched all the Swift observations for bursts by binning the source region light curves above 1 keV into 0.01 s, 0.1 s, and 1.0 s bins. The counts in each bin were compared to the mean count rate of its GTI, assuming the Poisson distribution. We found no significant bursts in the Swift XRT data. Note that the Swift observations did not cover the NuSTAR-detected burst times presented in Table 2.

In order to derive the timing solution and to search for glitching activity, we barycentered the source events using the location of 1E 1841–045, R.A. = 18^h41^m1934, decl. = −4°56'112. We then extracted times of arrival (TOAs) using a maximum likelihood (ML) method (see Livingstone et al. 2009; Scholz et al. 2012). The ML method compares a continuous model of the pulse profile to the profile obtained by folding a single observation. These TOAs were fitted to a pulse arrival model (e.g., the quadratic function in Section 3.2.1) using the TEMPO2 (Hobbs et al. 2006) pulsar timing software package. We find that a timing model consisting of ν and $\dot{\nu }$ does not fit the data well (Figure 11(b) and "Fit 1" in Table 5), as was already observed in the previous RXTE monitoring observations (e.g., see Dib & Kaspi 2014). We need to include 12 frequency derivatives in order to achieve an acceptable fit (i.e., χ²/dof ∼ 1) with Gaussian residuals (Figure 11 (c) and "Fit 2" in Table 5). Note that the timing solutions presented in Table 5 are valid only over the time interval of the monitoring campaign.

Table 5.  Timing Parameters for 1E 1841–045

R.A. (J2000)	18h41m1934
Decl. (J2000)	−4°56'112
MJD range	55,795–56,799
Epoch (MJD)	56,300
Fit 1
ν(s−1)	0.084 806 860 6(9)
$\displaystyle \frac{{d}^{1}\nu }{{{dt}}^{1}}$ (s−2)	−2.9121(8) × 10−13
rms (s)	4.82
χ2/dof	2889.05/53
Fit 2
ν(s−1)	0.084 806 897(4))
$\displaystyle \frac{{d}^{1}\nu }{{{dt}}^{1}}$ (s−2)	−2.985(12) × 10−13
$\displaystyle \frac{{d}^{2}\nu }{{{dt}}^{2}}$ (s−3)	6.8(23) × 10−22
$\displaystyle \frac{{d}^{3}\nu }{{{dt}}^{3}}$ (s−4)	4.6(8) × 10−28
$\displaystyle \frac{{d}^{4}\nu }{{{dt}}^{4}}$ (s−5)	−8.7(16) × 10−35
$\displaystyle \frac{{d}^{5}\nu }{{{dt}}^{5}}$ (s−6)	−2.5(5) × 10−41
$\displaystyle \frac{{d}^{6}\nu }{{{dt}}^{6}}$ (s−7)	6.0(12) × 10−48
$\displaystyle \frac{{d}^{7}\nu }{{{dt}}^{7}}$ (s−8)	1.2(3) × 10−54
$\displaystyle \frac{{d}^{8}\nu }{{{dt}}^{8}}$ (s−9)	−3.4(7) × 10−61
$\displaystyle \frac{{d}^{9}\nu }{{{dt}}^{9}}$ (s−10)	−4.5(10) × 10−68
$\displaystyle \frac{{d}^{10}\nu }{{{dt}}^{10}}$ (s−11)	1.4(3) × 10−74
$\displaystyle \frac{{d}^{11}\nu }{{{dt}}^{11}}$ (s−12)	9.5(22) × 10−82
$\displaystyle \frac{{d}^{12}\nu }{{{dt}}^{12}}$ (s−13)	−3.3(8) × 10−88
rms (s)	0.43
χ2/dof	49.34/42

Notes. All errors are TEMPO2-reported 1σ errors.

Download table as:  ASCIITypeset image

Motivated by the apparent "kink" in the residuals of the simple spin-down model around MJD 56,100, we attempted to fit a glitch at the epoch of the kink. However, the data are better fit using a model with four frequency derivatives (rms residual of 0.97 s) versus a glitch model (rms residual of 0.99 s). Therefore, we do not need to invoke a sudden event to explain the measured TOAs. We present our best timing solutions in Table 5.

Magnetars often show bursting behavior in the X-ray band, which may be caused by instabilities inside the neutron star or its magnetosphere (Thompson & Duncan 1995; Lyutikov 2003; Woods & Thompson 2006). The time profiles and spectra of bursts show significant diversity. Woods et al. (2005) suggested that there are two types of magnetar bursts, one having significant tail emission and the other having orders of magnitude smaller tail emission. The authors attributed the former to crustal activity and the latter to magnetospheric activity. Lenters et al. (2003) found a strong correlation between bursts and tail energies in the magnetar SGR 1900+14.

The X-ray bursts from 1E 1841–045 in 2013 September have very short rise and fall times, with T₉₀ of 0.01–0.6 s, and hard spectra (Γ = 1–2 or kT ∼ 3–5 keV). The blackbody temperature we measured for burst 5 (4.8 ± 0.5 keV; see Table 3) is consistent with that of the colder blackbody measured with the GBM (kT_l = 5.3 ± 0.2 keV; Collazzi et al. 2014).

Kumar & Safi-Harb (2010) reported detection of emission lines at 27, 40, and 60 keV in the 2010 May burst spectrum with Swift BAT, although these are argued against later by Lin et al. (2011). Interestingly, all the lines are in the NuSTAR band and could be detected by NuSTAR if they appeared again. However, we do not see evidence of line emission. Therefore, we estimate 90% upper limits on any 27 keV Gaussian line flux to be 0.24 photons cm⁻² s⁻¹ and 0.11 photons cm⁻² s⁻¹ in the brightest burst spectrum for the blackbody and the power-law continuum models, respectively. Note that we are not able to compare our results with those of Kumar & Safi-Harb (2010) since they did not present the line flux.

An extended tail is reliably detected only in energetic burst 5, and we find a hint of a tail in burst 6. Note that these two bursts were also observed by the GBM and had significant flux above the NuSTAR band (Collazzi et al. 2014); they are the two most energetic bursts in our sample. Thus, the tail brightness and the burst energy we measure for the 1E 1841–045 bursts seem to agree qualitatively with the correlation reported for SGR 1900+14 and SGR 1806–20 (Lenters et al. 2003; Göğüş et al. 2011). A similar trend was also seen in the recent bursts from 1E 1048.1–5937 (An et al. 2014). We note, however, that the large energy seen by GBM in bursts 5 and 6 indicates a tail-to-burst ratio (E_tail/E_burst ∼ 2 × 10⁻²) that is much lower than in 1E 1048.1–5937 (E_tail/E_burst ∼ 5–60). This confirms the known diversity of tails of magnetar bursts (Kaspi et al. 2004; Woods et al. 2005). Whether the tail-to-burst ratio follows a bimodal or a random distribution is not yet clear, and further investigation is required.

In contrast to the burst tails observed by NuSTAR for 1E 1048.1–5937, the burst tail in 1E 1841–045 shows no clear correlation between spectral hardness and flux. This correlation was also absent in some of the long-term (months to years) flux relaxation of other magnetars (e.g., An et al. 2012).

The flux evolution in the tail of burst 5 followed a power-law decay with a decay index of 0.45 ± 0.10. The flux decay is similar to the tails of bursts from SGR 1900+14 (decay index of 0.43–0.7; Lenters et al. 2003). A significantly faster decay was observed for burst tails in 1E 1048.1–5937 (decay index of 0.8–1; An et al. 2014). It is possible that the tail contains many unresolved weaker bursts that affect the observed flux decay. It is still unclear what controls the resulting decay rate and why it is significantly different in 1E 1048.1–5937 and 1E 1841–045.

Furthermore, we find that the tail spectra in 1E 1841–045 are similar to (or slightly softer than) the persistent emission. In contrast, the tail spectra in 1E 1048.1–5937 were significantly harder than its persistent emission. This further contributes to the diversity of magnetar bursts. For instance, 1E 2259+586 exhibited bursts with and without kilosecond-long tails, and the tail emission was observed to soften with decreasing flux (Kaspi et al. 2004). The observed diversity of magnetar bursts is not well explained by current theoretical models.

It is known that the pulse profiles of magnetars can look significantly different in different energy bands (den Hartog et al. 2008). This is also true for 1E 1841–045 (Figure 5). An interesting feature of the pulse profile is the double-peaked structure in the narrow band of ∼24–35 keV. A13 found this in the previous NuSTAR observation and suggested that it may be caused by a spectral feature. The new long observation confirms the double-peaked structure. It shows, however, that the change in the pulse profile is not as sharp as was suggested previously, and its shape may be more complicated, as a hint of another peak appearing between the two peaks is seen in several energy bands (e.g., in the 33–35 keV profile in Figure 5).

The pulsed fraction (and the rms pulse amplitude) increases with photon energy below ∼8 keV (Figure 5). We found no significant increase in the pulsed fraction above 8 keV, in contrast with previous measurements (Kuiper et al. 2006). This conclusion is insensitive to the choice of energy bins in Figure 5, which can affect the measured values in individual bins, but not the general trend. Our spectral analysis also suggests that the pulsed and steady components of the flux have similar power-law photon indices at high energies >15 keV (see Section 3.3 and Table 4). In the coronal outflow model, this implies that one of the angles, ${\alpha }_{\mathrm{mag}}$ and ${\beta }_{\mathrm{obs}}$, is small, or the emission region (${\theta }_{j}$) is broad.

We have measured the phase-averaged, phase-resolved, and pulsed spectra of the source. Our best-fit parameters are slightly different from those reported by A13. The use of the new CALDB for the NuSTAR data may have some effect on the obtained spectral shape; however, we showed in Section 3.2.3 that this effect is small.

We do not see any significant change in the source flux among the observations taken over 1 yr despite the fact that bursts were detected in some observations but not in the others. We further compared the NuSTAR-measured spectrum with that reported previously (Γ = 1.32 ± 0.11 and L_{10−100 keV} = 3.0 × 10³⁵ erg s⁻¹ for an assumed distance of 6.7 kpc; Kuiper et al. 2006) and find that our measurement (Γ = 1.37 ± 0.01 and L_{10−100 keV} = 3.0 × 10³⁵ erg s⁻¹) is fully consistent with the previous values, suggesting that the source hard-band spectrum has been stable over 10 yr. The same trend in the soft band (1–10 keV) is seen in the Swift monitoring data (see Figure 11). Stability in the soft-band pulsed flux has been reported by Dib & Kaspi (2014) and Zhu & Kaspi (2010) for a longer period. We note that Kumar & Safi-Harb (2010) reported a possible increase in the soft-band flux for 1E 1841–045 (<10 keV) associated with the source burst activity. However, Lin et al. (2011) showed that there is no significant change in the soft-band flux over a 1400-day interval, including the burst period reported by Kumar & Safi-Harb (2010).

The new measurements of the phase-resolved spectrum agree with the previous NuSTAR observation. The spectrum is harder near the pulse peaks, with smaller Γ_h and greater E_break. We also note a possible hardening of the phase-averaged spectrum above >15 keV. However, we cannot reliably measure the spectrum curvature in the hard-X-ray band owing to the statistical noise.

We have searched for the spectral feature that was suggested as a possible explanation for the change in the pulse profile near 30 keV (A13). We find that the phase-resolved spectra are all well fit by a blackbody plus broken power law or a blackbody plus double power law; no emission or absorption line is required or hinted at by the fit residuals. We also searched for a spectral feature that shifts with rotational phase, as was seen in SGR 0418+5729 (Tiengo et al. 2013), but did not find such a feature (Figure 9).

The more accurate phase-resolved spectrum obtained from the 350 ks NuSTAR observation provides a new opportunity to apply the coronal outflow model. We found that the model still fits the observed hard-X-ray spectrum, and it does so in a small region of the parameter space (Figure 10), which allowed us to estimate the size of the active j-bundle, the dissipated power in the magnetosphere, and the angles between the magnetic axis, the rotation axis, and the line of sight. The increase in exposure by a factor of ∼6 excluded at more than the 3σ level the second solution for the outflow model that was found in A13.

Interestingly, adding the new free parameter ${\rm{\Delta }}{\theta }_{j}$ does not introduce new acceptable solutions (with a p-value above 10⁻³) and does not significantly affect the best solution—the best fit shows that the footprint of the j-bundle on the star is a broad ring, hardly distinguishable from a polar cap. This contrasts with 1E 2259+586, where a thin ring with ${\rm{\Delta }}{\theta }_{j}/{\theta }_{j}\lt 0.2$ is clearly statistically favored over a polar cap. This diversity may point to different distributions of the crustal magnetic stresses that are responsible for the magnetospheric twisting.

The soft-X-ray component (below ∼10 keV) is well fitted by the sum of two blackbodies. The best-fit model is similar to that found in A13. The cold blackbody covers a large fraction of the star area, ${{\mathcal{A}}}_{c}\approx 0.42{{\mathcal{A}}}_{\mathrm{NS}}$, where ${{\mathcal{A}}}_{\mathrm{NS}}$ is the area of the neutron star with an assumed radius R_NS = 10 km. The emission area of the hot blackbody, ${{\mathcal{A}}}_{{\rm{h}}}\approx 0.024{{\mathcal{A}}}_{\mathrm{NS}}$, is comparable with the area of the outflow footprint ${{\mathcal{A}}}_{j}=\pi {\mathrm{sin}}^{2}{\theta }_{j}\simeq 0.014\ {({\theta }_{j}/0.24)}^{2}$. This is consistent with the coronal outflow model, where the footprint of the j-bundle is expected to form a hot spot, as some particles accelerated in the j-bundle flow back to the neutron star and bombard its surface.

Figure 10 suggests that the coronal outflow correctly describes the hard-X-ray source as a decelerating e^± outflow ejected from a discharge zone near the star. The outflow parameters inferred from the fit of the phase-resolved data have rather small statistical uncertainties, and the results reveal a puzzling feature of 1E 1841–045. Using Equation (48) in Beloborodov (2009), one can infer the discharge voltage in the active j-bundle: Φ ≈ 10¹¹ ψ⁻¹ V, where Ψ is the twist implanted in the magnetosphere, which does not exceed ψ_max ∼ 3 (Parfrey et al. 2013). Owing to our refined constraint on ${\theta }_{j}$ and the strong dependence of Φ on ${\theta }_{j}$ (${\rm{\Phi }}\propto {\theta }_{j}^{4}$), the inferred voltage is one order of magnitude higher than that given in A13. The high voltage is surprising in two ways: (1) it exceeds the expected threshold for e^± discharge (Beloborodov & Thompson 2007) by at least a factor of 10, and (2) it implies a short timescale for ohmic dissipation of the magnetospheric twist t_diss ≈ 0.1 ψ² yr (Equation (50) in Beloborodov 2009). Thus, without continued energy supply, the j-bundle is expected to untwist on a year timescale or faster, which is not observed—the persistent X-ray emission from 1E 1841–045 has been stable for at least one decade.

This puzzle is related to a more general question: why are some magnetars transient and others persistent? The magnetospheric activity should be fed by magnetic energy pumped from the star by its surface motions. The surface motions are caused by the crust yielding to accumulating magnetic stresses inside the magnetar. Energy supply to the twisted magnetosphere may or may not be intermittent, depending on the mechanism of the crustal motions. Beloborodov & Levin (2014) have recently shown that the crust can yield through a thermoplastic instability that launches a slow wave resembling the deflagration front in combustion physics. The thermoplastic wave rotates the crust and twists the magnetosphere. This mechanism naturally triggers the outbursts observed in transient magnetars, yet it is still unclear how the quasi-steady state observed in 1E 1841–045 is formed. It could be formed by frequently repeated energy supply to the j-bundle combined with some form of feedback on its untwisting rate.

We did not detect any changes in the 1–10 keV source flux during the 3 yr of Swift XRT monitoring since 2011. In particular, the flux did not change significantly during the bursting period, which is consistent with the stability of persistent emission in the NuSTAR data. Our measurements do not agree with the observation that the source persistent emission properties changed due to bursts (Kumar & Safi-Harb 2010). Note that Lin et al. (2011) also observed no change in the persistent properties of the source after bursts. The flux stability also implies that the source seems to be in a perpetually bursting state, and perhaps this is a difference between the classical bright anomalous X-ray pulsars monitored with RXTE and "transient" sources, i.e., perhaps the bursting indicates that there is more heat dissipating and hence more magnetic activity in this magnetar, consistent with the higher L_X in "quiescence."

The pulse timing behavior of 1E 1841–045 is known to be noisy, and Dib & Kaspi (2014) had to use five derivative terms to fit the timing data. Our timing solution requires 12 frequency derivatives. Note that although Dib & Kaspi (2014) used a much longer data span (13 yr), they break the data into smaller pieces that span 1–3 yr, similar to ours. The Swift monitoring data require more frequency derivatives than the RXTE data do owing to the large kink at MJD 56,100, which may indicate a glitch. We investigated the possibility that the kink is caused by a glitching activity but did not find any evidence of a glitch during the monitoring period of 3 yr including the bursting period.