NuSTAR OBSERVATIONS OF THE MAGNETAR 1E 2259+586

For the timing analysis, we extracted photons in a circular region of 60'' radius around the nominal source position for the NuSTAR data, while for Swift we used a 22'' radius. The energy bands used for spectral analysis were 5–79 keV and 0.5–10 keV for NuSTAR and Swift, respectively. Since 1E 2259+586 is located in CTB 109, the SNR background must be considered. Sasaki et al. (2004) pointed out that no emission above 4 keV was observed in the Chandra data apart from the pulsar and the NuSTAR data agrees with this. Since in one of the two NuSTAR modules (FPMA) stray light is observable at the edge of the field of view, we excluded this region when defining the background extraction area. For NuSTAR, background was extracted from a circular region of radius 100'' nearby the source region. We also compared our results using these background regions to those obtained using the mission's background model nuskybgd (Wik et al. 2014), and confirmed that all results were compatible within the corresponding uncertainties.

Due to the good spatial resolution of the Swift XRT, background could be extracted from an annular region around the NS with an inner radius of 80'' and an outer radius of 200''. Before continuing we checked if any pile-up events were present in the data. We followed the standard Swift procedure¹⁵ and determined the count rate in a circular region of 472 radius (20 pixels) centered on the source to be above 0.5 counts s⁻¹, indicating that the central 1–2 pixel radius region of the bright core might suffer from pile-up. We removed this area for the subsequent analysis by substituting the circular source extraction region by an annulus with inner and outer radii of 5'' and 22'', respectively. We also re-generated the ARFs to correct the flux of the spectrum for the loss of counts due to the use of an annular region as well as the correction factor for the light curves.¹⁶

In the next step, we applied a barycentric correction to the selected source events using the multi-mission tool barycorr with the corresponding orbital and clock correction files at the Chandra position of α = 23^h01^m08295 and δ = +58°52'4445 (J2000.0) reported by Patel et al. (2001). Following the H-test method (De Jager et al. 1989), we searched for pulsations and determined the best period for the significant pulsations to be P(NuSTAR) = 6.97914(2) s and P(Swift) = 6.97915(2) s. The measured periods were in agreement with those obtained from the ephemeris of the Swift monitoring program of 1E 2259+586 (see Table 2).

The resulting pulse profiles are shown in Figure 3 for six energy ranges: 3–4 keV, 4–8.3 keV, 8.3–11.9 keV, 11.9–16.3 keV, 16.3–24.0 keV, and 24.0–79.0 keV. The bands up to 24.0 keV were chosen to coincide with those reported by Kuiper et al. (2006) for direct comparison. All pulse profiles, with the exception of the highest energy band (24–79 keV), are background-subtracted. The statistical significance of the observed pulsations in all energy bands is larger than 99% (p value¹⁷ smaller than 0.01). For the lower five energy bands, the pulse profiles agree generally with those reported by Kuiper et al. (2006); however, in the energy ranges between 11.9 and 24.0 keV slight differences are apparent. In the highest energy band (24–79 keV), significant pulsations are observable for the first time. We note, however, that due to limited statistics we did not subtract the background in order to avoid negative bin counts. The nominal background count rate in this case is 6.56 × 10⁻⁴ counts s⁻¹.

Zoom In Zoom Out Reset image size

As previously pointed out by Kuiper et al. (2006), there is a gradual change in pulse morphology observable with increasing energy. The double-peak profile at low energies (below ∼8 keV) evolves into a less pronounced double peak or possibly a single-peak profile for intermediate energies (about 8–16 keV) and back to a double-peak structure at higher energies. A possible way to describe this development is that while for energies below ∼8 keV, the pulse near phase 0.9 (peak A") tends to dominate, it becomes less significant for increasing energy with respect to the second peak ("peak B") near phase 1.45. In the intermediate energy range (8–16 keV), both peaks are approximately equally pronounced, while peak B starts to slightly dominate at energies above ∼16 keV. A difference between NuSTAR and RXTE profiles is apparent in the energy bands of 11.9–16.3 keV and 16.3–24.0 keV. While Kuiper et al. (2006) observe a clearly dominating peak B in this energy range, our analysis shows only a small lead for this pulse over peak A.

In order to quantify the energy dependence of the pulse morphology, we decompose the pulse profiles in terms of their Fourier components. This decomposition is shown in Figure 4. We find six harmonics to be sufficient, and we included these in the plots. We observe that the ratio of the second harmonic to the first harmonic changes with energy. For energies below about 8 keV, the second harmonic dominates over the first, while in the intermediate range of 8–16 keV the first harmonic takes the lead, before the second harmonic dominates again for 16–24 keV. If the main contribution originates from the first harmonic number, this indicates that there is either one single peak in the pulse profile or a double peak with maxima of similar magnitude, while a leading second harmonic number points toward a double peak with pulses of different height and/or shape. This Fourier decomposition therefore reflects the change in dominance of one peak over the other in the double-peak morphology or, likewise, a possible change from a two-peak structure to one peak and back. For the NuSTAR band (E > 24 keV), the first and third harmonics dominate (2, 4, and 6 contribute at the 10% level). When interpreting this power spectrum (lower right plot of Figure 4) one should keep in mind that background was not subtracted and the peak-to-valley ratio is not very large. Here the first and (small) second harmonic take care of the wide peak B (from phase 1–1.8, including the small feature in the main minimum), while the third harmonic describes the narrower, subdominant peak A.

Zoom In Zoom Out Reset image size

As previously reported by Kuiper et al. (2006) as well as Gavriil et al. (2002), a small pulse-like feature can be observed in the main minimum around phase 0.60–0.65. This feature is the main reason for the presence of higher orders of harmonics (n ⩾ 3) in the Fourier decomposition (see Figure 4) of the pulse profile.

We also note that there is a hint of a small feature in the intermediate minimum between peak A and B around phase 1.2, best seen in the lowest two energy bands. With the current statistics, however, we do not consider this feature to be significant.

Table 3. Mean Peak Position for Peaks A and B in Various Energy Bands

Energy Band	Position	Error	Position	Error
(keV)	Peak A	Peak A	Peak B	Peak B
3.0–4.0	0.93	<0.01	1.39	<0.01
4.0–8.3	0.94	<0.01	1.38	<0.01
8.3–11.9	0.94	$\hphantom{<}0.01$	1.35	0.01
11.9–16.3	0.94	$\hphantom{<}0.01$	1.35	0.01
16.3–24.0	0.93	$\hphantom{<}0.01$	1.36	0.01

Download table as:  ASCIITypeset image

Another observation is a possible shift of the peak B position in phase toward lower values with increasing energy. Dashed lines near phase 0.9 and 1.45 have been added to Figure 3 to guide the eye for alignment comparison of peak A and peak B. In order to test this pulse-shift hypothesis we determined the mean value of the peak position for each peak in various energy bands by considering the unbinned group phase between phase 0.75–1.1 for peak A and 1.2–1.55 for peak B. In addition, we crosschecked our results by fitting a Gaussian to the peaks. The results for the mean and its error are listed in Table 3 and indicate that peak A remains at the same phase position of 0.937(17) throughout all energy ranges (0.933(3) for E ≲ 8 keV, 0.940(16) for E ≳ 8 keV). For peak B, there is a slight tendency to change position toward lower phase values for higher energies. Below ∼8 keV, peak B can be found at 1.384(3), above these energies it is located at 1.349(17). We note that the shift is not very significant at the given statistics, and more data are needed to confirm or rule out the observed tendency.

In addition, we used the ephemeris-folded profiles to determine the pulsed fraction (PF) as a function of energy for our observations. A pulse profile that changes with energy or time makes it challenging to determine the PF of the source accurately and different methods are commonly used in the literature (see, e.g., Archibald et al. 2007, A. M. Archibald et al. 2014, in preparation for a discussion). All of them have their advantages and disadvantages and therefore there are some caveats to keep in mind. The area PF is probably the most natural and physically meaningful definition. It is defined as the difference between the pulsed flux and the constant flux integrated over a full phase cycle and can be calculated according to

$$\begin{equation} {\rm PF}_{\rm area}=\frac{\frac{1}{N}\displaystyle \sum \nolimits _{j=1}^{N} p_{j}-p_{\rm min}}{\frac{1}{N}\displaystyle \sum \nolimits _{j=1}^{N} p_{j}}, \end{equation} \tag{ 1 }$$

where p_j is the number of events in the jth phase bin, N is the total number of bins, and p_min is the minimum flux. Note that the determination of the true minimum and its error is the main challenge and the biggest caveat when using PF_area, because both noise and binning show a tendency to bias the values of the area PF upward. In Figure 3 we included the results for p_min as the minimum off-pulse value determined from the Fourier representation (horizontal magenta line). The results for the area PF using the p_min determined in this way are shown in Figure 5. Note that the NuSTAR energy bands here (black squares) correspond to those presented in Kuiper et al. (2006) for the RXTE Proportional Counter Array (PCA, five bands in the range from 3.0–24.0 keV) and HEXTE data (14.8–27.0 keV). The Swift data (blue diamonds) covers the two lowest RXTE PCA bands in Kuiper et al. (2006) plus an additional band from 0.3–3.0 keV. Background was subtracted in all energy bands.

Zoom In Zoom Out Reset image size

The small inset plot in the same figure displays the standard PF (peak-to-peak PF), conventionally defined as ratio of peak flux minus minimum flux to the total flux

$$\begin{equation} {\rm PF}_{\rm P2P}=\frac{F_{\rm max}-F_{\rm min}}{F_{\rm max}+F_{\rm min}}. \end{equation} \tag{ 2 }$$

This quantity is extensively used in the literature since it is straightforward to calculate, although the difficulty in determining the maximum and minimum flux values accurately is biased in a way similar to the area PF. Furthermore, the PF_P2P cannot provide information on the total energy since the peak width does not feed into the flux variation. The PF_area and PF_P2P derived from our data agree with each other within uncertainties and show the same trend: the PF increases dramatically with increasing energy. For energies around 10 keV, PF_area is 60% ± 12%, while for energies around 20 keV the fraction reaches 96% ± 14%.

Following the method presented in Gonzalez et al. (2010; see also the forthcoming publication by A. M. Archibald et al. 2014, in preparation for details), we additionally used a root-mean-square (rms) estimator as a measure for the deviation of the pulsed flux from its mean defined as

$$\begin{equation} {\rm PF}_{\rm {\rm rms}}=\frac{ \sqrt{2\displaystyle \sum \nolimits _{k=1}^{6}\left(\left(a_{k}^{2}+b_{k}^{2}\right)-\left(\sigma _{a_{k}}^{2}+\sigma _{b_{k}}^{2}\right)\right)}}{a_{0}}, \end{equation} \tag{ 3 }$$

where the Fourier coefficients a_k and b_k are

$$\begin{equation} a_{k}=\frac{1}{N}\displaystyle \sum \limits _{j=1}^{N}p_{j}\cos \left(\frac{2\pi kj}{N}\right), \end{equation} \tag{ 4 }$$

$$\begin{equation} b_{k}=\frac{1}{N}\displaystyle \sum \limits _{j=1}^{N}p_{j}\sin \left(\frac{2\pi kj}{N}\right) \end{equation} \tag{ 5 }$$

and $\sigma _{a_{k}}$, $\sigma _{b_{k}}$ are the uncertainties of a_k and b_k, respectively,

$$\begin{equation} \sigma _{a_{k}}^{2}=\frac{1}{N^{2}}\displaystyle \sum \limits _{j=1}^{N}\sigma _{p_{j}}^{2}\cos ^{2}\left(\frac{2\pi kj}{N}\right), \end{equation} \tag{ 6 }$$

$$\begin{equation} \sigma _{b_{k}}^{2}=\frac{1}{N^{2}}\displaystyle \sum \limits _{j=1}^{N}\sigma _{p_{j}}^{2}\sin ^{2}\left(\frac{2\pi kj}{N}\right). \end{equation} \tag{ 7 }$$

p_j is again the number of events in the jth phase bin, $\sigma _{p_{j}}$ is the uncertainty of p_j, and N is the total number of bins. The number of Fourier harmonics included is n = 6. Note that ${\rm PF}_{\rm {\rm rms}}$ is the rms deviation of the flux from its average and can be defined in the Fourier domain. It does not cover the full range of values from 0–1 without additional (pulse shape dependent) renormalization, i.e., it is not a PF in the conventional sense, even though it is often referred to as the rms PF in the literature. This definition yields a very robust measure for the pulsed power that is less sensitive to noise and therefore significantly less biased than the area PF and more meaningful than the peak-to-peak PF. Figure 6 shows the rms variation for NuSTAR (black squares) and Swift data (blue diamonds). Note that the energy bands again correspond to those presented in Kuiper et al. (2006) for the RXTE PCA (five bands in the range from 3.0–24.0 keV) and HEXTE data (14.8–27.0 keV) and represent background-subtracted data in all energy bands. Since the definition of ${\rm PF}_{\rm {\rm rms}}$ differs from those of PF_area and PF_P2P as mentioned above, ${\rm PF}_{\rm {\rm rms}}$ does not have the same values as the other definitions in each energy band, but the same trend is observable in all cases: the PF increases dramatically with increasing energy. Without renormalization we determined that the ${\rm PF}_{\rm {\rm rms}}$ is 32% ± 9% for energies around 10 keV, while for energies around 20 keV it reaches 71% ± 15%.

Zoom In Zoom Out Reset image size

3.2.1. Phase-averaged Spectral Analysis

In this section, we use the Swift and NuSTAR observations to derive a phase-averaged spectrum. The same extraction regions as defined in Section 3.1 were used and pile-up for the Swift data was corrected for as described above.

We started with fitting the Swift data only in order to compare the results to previous soft X-ray measurements (Zhu et al. 2008) as well as to obtain a best-fit N_H to be used for the combined of Swift and NuSTAR data sets. All extracted spectra were grouped to have at least 50 counts per bin using the grppha tool. We fit the resulting spectra with XSPEC v12.8.1 using an absorbed BB plus PL model (tbabs*(bbody+powerlaw)) as well as an absorbed double BB (tbabs*(bbody+bbody)). The Swift spectrum is a little softer (Γ = 4.1(1)) than previously reported (Γ = 3.75(3); Zhu et al. 2008), but consistent within uncertainties with these earlier results. Our best-fit result yields N_H = 1.10(6) × 10²² cm⁻² in the case of the BB plus PL and N_H = 0.60(4) × 10²² cm⁻² for the double-BB model. The first value agrees well with estimates of Durant & van Kerkwijk (2006) obtained from fitting individual absorption edges of O, Fe, Ne, Mg, and Si in the XMM-Newton Reflection Grating Spectrometer spectra. The latter value is in agreement with Zhu et al. (2008) and consistent with the best-fit N_H value for CTB 109, N_H = (0.5–0.7) × 10²² cm⁻², as measured by Sasaki et al. (2004).

In a next step we fit the spectra from 0.5–79 keV by using the NuSTAR and Swift data sets together. All model parameters were tied except for the cross-normalization factors, which were set to 1.0 for NuSTAR's FPMA and left to vary for FPMB and Swift. Spectra were rebinned to have at least 50 counts per bin. We started out with fitting an absorbed double BB (tbabs*(bbody+bbody)) but the fit was poor, yielding a χ² per degree of freedom (dof) of 1448/827. The fit improved slightly by using an absorbed BB plus PL (tbabs*(bbody+powerlaw)), with χ² over dof of 1152/827. Adding an additional component improved the fit significantly and therefore a more complex model (BB+double PL or BB+broken PL, χ²/dof = 689/825 or 747/825, respectively) than suggested by previous data was favored. The F-test supports this conclusion, yielding a probability significantly smaller than 0.05. All spectral fit results can be found in Table 4, while Figure 7 shows our best-fit model for the spectra of NuSTAR and Swift.

Zoom In Zoom Out Reset image size

Table 4. Best-fit Parameters for 1E 2259+586 for Phase-averaged and Phase-resolved Analysis

Phase	Dataa	Energy		NHc	kT	Γs/kTd	$E_{\rm break^e}$	Γhf	Fluxg	Ratioh	χ2/dof
		(keV)	Fit Modelb	(× 1022 cm−2)	(keV)	(-/keV)	(keV)
0–1	S	0.5–10.0	BB+BB	0.60(4)	0.35(2)	0.71(6)	⋅⋅⋅	⋅⋅⋅	1.9(2)	1.00(1)	331/537
0–1	S	0.5–10.0	BB+PL	1.10(6)	0.42(2)	4.1(1)	⋅⋅⋅	⋅⋅⋅	2.1(1)	1.71(9)	402/537
0–1	S/N	0.5–79.0	BB+BB	0.60(4)	0.423(3)	1.78(4)	⋅⋅⋅	⋅⋅⋅	2.25(5)	0.2(1)	1448/827
0–1	S/N	0.5–79.0	BB+PL	1.10(6)	0.38(1)	3.91(3)	⋅⋅⋅	⋅⋅⋅	2.2(1)	2.5(3)	1152/827
0–1	S/N	0.5–79.0	BB+2PL	1.10(6)	0.410(7)	4.34(3)	⋅⋅⋅	0.4(1)	2.36(7)	0.31(2)	689/825
0–1	S/N	0.5–79.0	BB+bPL	1.10(6)	0.409(9)	4.08(3)	11.5(3)	1.2(1)	3.0(1)	2.2(2)	747/825
Peak A	S/N	0.5–79.0	BB+2PL	1.10(6)	0.416(6)	4.57(7)	⋅⋅⋅	1.0(3)	0.68(3)	0.29(4)	716/803
Peak B	S/N	0.5–79.0	BB+2PL	1.10(6)	0.399(5)	4.41(5)	⋅⋅⋅	0.3(1)	0.73(3)	0.53(4)	645/812
Off-pulse	S/N	0.5–79.0	BB+2PL	1.10(6)	0.404(6)	4.87(9)	⋅⋅⋅	1.0(4)	0.48(3)	0.37(6)	879/1031
Peak A	S/N	0.5–79.0	BB+bPL	1.10(6)	0.397(5)	4.37(4)	10.8(5)	1.38(1)	1.51(4)	0.78(2)	776/803
Peak B	S/N	0.5–79.0	BB+bPL	1.10(6)	0.386(5)	4.26(5)	10.9(4)	0.9(2)	1.46(2)	0.85(1)	696/812
Off-pulse	S/N	0.5–79.0	BB+bPL	1.10(6)	0.393(5)	4.66(5)	9.6(4)	1.4(4)	1.43(3)	0.47(1)	911/1031

Notes. Numbers in parentheses indicate the 1σ uncertainty in least significant digits. Cross-normalization factors were used when data from NuSTAR and Swift were combined. These factors were set to 1.0 for module A of NuSTAR and to 1 for Swift if no NuSTAR data were included. Cross-normalization factors were frozen to values obtained for phase-averaged spectral analysis, when performing phase-resolved spectral analysis simultaneously for NuSTAR and Swift. Fluxes are unabsorbed fluxes measured using XSPEC's cflux model. ^aNuSTAR 5–79 keV (N), Swift 0.5–10 keV (S). ^bBlackbody (BB), power law (PL), broken power law (bPL). ^cN_H was frozen to the value obtained from the phase-averaged spectral analysis for the phase-resolved analysis. ^dSoft photon index Γ_s if BB+PL model is used. Blackbody temperature of second (hot) BB if BB+BB is fit. ^eBreak energy for the broken power law (where applicable). ^fPhoton index of the hard power-law component. ^gUnabsorbed flux in units of 10⁻¹¹ erg cm⁻² s⁻¹ from 2–10 keV for Swift and 2–79 keV if NuSTAR data are included for all model components (BB+BB, BB+PL, BB+bPL) or two PL components (BB+2PL). ^hRatio of flux in the 2–10 keV range for Swift and 2–79 keV range if NuSTAR is included. We show $F_{\rm BB(\rm hot)}/F_{\rm BB(\rm cold)}$ (BB+BB), F_PL/F_BB (BB+PL), F_bPL/F_BB (BB+bPL) and $F_{\rm PL(\rm hard)}/F_{\rm PL(\rm soft)} ({BB}+2{PL})$.

Download table as:  ASCIITypeset image

3.2.2. Phase-resolved and Pulsed Spectral Analysis

3.2.3. Spectral Fits with the Coronal Outflow Model

Similar to what was done by An et al. (2013) for NuSTAR observations of 1E 1841−045, we apply the coronal outflow model of Beloborodov (2013a) to the observations of 1E 2259+586. In this model the hard X-ray component is produced by a decelerating e^± outflow in a twisted magnetic loop (the "j-bundle"). The relativistic e^± flow is produced near the star with a Lorentz factor γ ≳ 10³ via continual discharge. As the outflow expands and fills the j-bundle, it experiences a radiative drag (deceleration) due to resonant scattering of thermal photons (with energies in the keV range) around the NS. Beloborodov (2013a, 2013b) showed that the Lorentz factor of the flow decreases proportionally to the local magnetic field

$$\begin{equation} \gamma \sim 100 \frac{B}{B_Q}, \end{equation} \tag{ 8 }$$

where $B_Q = m_e^2 c^3 / \hbar e \simeq 4.44 \times 10^{13}$ G. Close to the star, where B ≳ 10¹³ G, the energy of scattered photons,

$$\begin{equation} E_{\rm sc} \sim \gamma (B/B_Q) m_e c^2 \sim 50 (B/B_Q)^2 \ \mathrm{MeV} \sim 5 \gamma ^2 \ \mathrm{keV} \,, \end{equation} \tag{ 9 }$$

is large enough to immediately convert to e^± pairs. At larger radii, where B ≲ 10¹³ G, the scattered photons escape, and the flow radiates almost all its kinetic energy before it reaches the magnetic equator, where the e^± pairs annihilate. The predicted hard X-ray spectrum has an average photon index of −1.5 (this follows from Equations (8) and (9)) and cuts off in the MeV band. The hard X-rays are beamed along the magnetic field lines, and the observed spectrum varies greatly depending on the line of sight (see Figure 7 in Beloborodov 2013b).

To fit the NuSTAR data, we make the simple assumption of an axisymmetric j-bundle. The poloidal magnetic field is assumed to be close to the dipole configuration, and the magnetic dipole moment is fixed to the value inferred from spin-down, μ_sd ≃ 5.9 × 10³¹ G. The model is described by five parameters: (1) the angular position θ_j (magnetic colatitude) of the j-bundle footprint, (2) the angular width Δθ_j of the j-bundle footprint, (3) the power L of the e^± outflow along the j-bundle, (4) the angle α_mag between the rotation axis and the magnetic axis, and (5) the angle β_obs between the rotation axis and the observer's line of sight. In addition, the reference point of the rotational phase, ϕ₀, must be introduced as a free parameter when fitting the phase-resolved spectra. Note that Δθ_j = θ_j would describe a polar cap. In this paper, we allow the j-bundle to have a ring-shaped footprint, which corresponds to Δθ_j < θ_j.

To test observations against the model, we follow the two-step method proposed by Hascoët et al. (2014). We first explore the whole parameter space by fitting the phase-averaged spectrum of the total (pulsed+unpulsed) emission together with three phase-resolved spectra of the pulsed emission. We only consider data above 16 keV where the observed hard component starts to dominate over the soft component as apparent from Figure 7. The results are shown in Figure 8. We find that the model can fit the data. However, the fit does not give strong constraints on the parameters of the model; a large region in the α_mag − β_obs plane is allowed. The only significant constraint derived from this analysis is that the polar-cap shape of the j-bundle footprint does not give a good fit; we find that a ring-shaped footprint is statistically preferred with 0.4 ≲ θ_j ≲ 0.75 and Δθ_j/θ_j ≲ 0.2 (at the 1σ level).

Zoom In Zoom Out Reset image size

At the second step, we freeze the best-fit parameters for the outflow model and fit the phase-averaged spectrum in the 0.5–79 keV range using both NuSTAR and Swift data. We test different models for the soft component, taking into account the extension of the outflow contribution to low energies. Among these models BB_tail is a phenomenological modification of the BB, where the Wien tail is replaced by a PL smoothly connected to the Planck distribution (see Hascoët et al. 2014). The energy E_tail, at which the PL starts, is a free parameter. All results are summarized in Table 5.

Table 5. Parameters of the Best-fit Model for the Soft X-Ray Component (NuSTAR and Swift Phase-averaged Spectra)

Model	NH	kT1	kT2	Etail	L1a	L2a	Γ	νLνa	χ2/dof	p Value
	(1022 cm−2)	(keV)	(keV)	(keV)				@ 1 keV
BB+BB	0.44(2)	0.430(4)	1.22(2)	⋅⋅⋅	1.7(2)	0.073(3)	⋅⋅⋅	⋅⋅⋅	632/497	3.8 × 10−5
BB+PL	1.17(4)	0.43(1)	⋅⋅⋅	⋅⋅⋅	0.6(1)	⋅⋅⋅	4.28(5)	3.3(3)	526/497	0.18
BBtail	0.52(2)	0.401(6)	⋅⋅⋅	2.59(5)	1.13(9)	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	539/498	0.098
BB1+ BBtail, 2	0.62(4)	0.33(3)	0.65(9)	4.3(6)	1.1(1)	0.3(2)	⋅⋅⋅	⋅⋅⋅	522/496	0.20

Note. ^aIn units of 10³⁵(D/D₀)² erg s⁻¹, where D₀ = 4.0 kpc is the distance inferred for 1E 2259+586 (Tian et al. 2010).

Download table as:  ASCIITypeset image

The NuSTAR pulse profiles of 1E 2259+586 generally agree with those presented previously in Kuiper et al. (2006) based on RXTE data. A change in peak predominance occurs with increasing energy at the boundary between the high-energy end of the soft component (where thermal emission from the star is the dominant source) and the beginning of the hard component (where the e^± outflow becomes the dominant source).

The two pulses of the double-peak profile seen in 1E 2259+586 are separated by half a period and have similar peak fluxes. This might be an indication that the NS is seen from the same magnetic latitude¹⁸ (in absolute value) every half period (assuming that the emission pattern around the NS is relatively axisymmetric). This is possible if (and only if) α_mag ≈ π/2 (corresponding to a nearly orthogonal rotator) or β_obs ≈ π/2. Note that no specific line of sight is a priori favored and β_obs ≈ π/2 is the most probable observer angle. We also note that the center of pulse B shifts from 1.38 to 1.35(2) for energies above 8 keV, while the initially dominant peak A remains centered at the same phase independent of energy. The phase shift of peak B is, however, not significant with current statistics, and more data are needed to confirm or rule out the observed trend. The e^± outflow and its emission are probably not perfectly axisymmetric, which can easily result in an energy-dependent phase shift of the peak maxima. Modeling such a minor shift with a small number of parameters is difficult and would require a different, more complicated fitting procedure.

Both the area and the peak-to-peak PF as well as ${\rm PF}_{\rm {\rm rms}}$ show an increase with energy. The area PF is about 60% for energies around 10 keV and rises to close to 100% at 20 keV, ${\rm PF}_{\rm {\rm rms}}$ is around 30% and 70% at these energies, respectively. The Swift and NuSTAR values for area and peak-to-peak PF are consistent with those reported by Patel et al. (2001) as peak-to-peak PF for E = 0.5–7.0 keV of 35.8% ± 1.4%. The use of different analysis techniques precludes a direct comparison of the PFs we determined and those established in some previous publication (spectral PF of ≳ 43% for energies below 10 keV as reported by Kuiper et al. 2006). We can compare, however, the general behavior for the PF versus energy observed for 1E 2259+586 with that for other magnetars. Kuiper et al. (2006) reported an increase of PF with energy for 1E 1841−045, 4U 0142+61 and 1RXS J1708−4009. An et al. (2013) observed the same tendency for 1E 1841−045, although the trend was not quite as pronounced (24% ± 4% rms PF at 20 keV and 41% ± 18% at 80 keV). For 1E 2259+586 we clearly see a dramatic increase in PF following the general trend observed for other magnetars. Such an increase in the PF is naturally expected in the coronal outflow model: photons of higher energy are more strongly beamed along the magnetic field lines as they are produced where the flow moves with a higher Lorentz factor (see Equation (9)).

We found that the spectral parameters of Swift and NuSTAR for the phase-averaged spectrum (see Section 3.2.1) agree well with those reported by Patel et al. (2001) and Zhu et al. (2008). These results support the fact that especially the soft-band spectrum (below ∼10 keV) has been stable over a period of ∼13 yr from Chandra observations in 2000 (Patel et al. 2001) to 2013 despite an antiglitch (Archibald et al. 2013) and glitches (Kaspi et al. 2003; Icdem et al. 2012), that temporarily altered the spectral parameters.

With the new data we were able to test if an additional PL component is required and find that it is. Therefore 1E 2259+586 exhibits the general tendency for magnetars to get harder at X-ray energies above 10 keV. We parameterized the NuSTAR plus Swift spectrum with a BB plus double PL model as well as a BB plus a broken PL. This phenomenological parameterization is a convenient way to characterize the presence of a hard X-ray component; the physical model of a coronal outflow is discussed in Section 4.3 below.

We note that the results of our analysis support the anticorrelation of spectral turnover (Γ_s − Γ_h) with magnetic field B and spin-down $\dot{\nu }$ suggested by Kaspi & Boydstun (2010). The authors demonstrated that there is a trend suggesting that those magnetars with the highest (lowest) B or $\dot{\nu }$ have the smallest (largest) spectral turnover. In Figure 9, we summarize the results of Kaspi & Boydstun (2010) for pulsed (green triangles) and total fluxes (blue circles) and add our new findings for 1E 2259+586 (red star) which further strengthens their case. Errors shown are 1σ uncertainties. It is interesting to point out that the observed spectral turnover for energetic rotation-powered pulsars (RPPs) is practically zero despite their much lower magnetic fields. Therefore Kaspi & Boydstun (2010) concluded that the production mechanism for X-rays above 10 keV in RPPs differs significantly from the one in magnetars. We also find a value for the hardness ratio, i.e., F_h/F_s (flux ratio of hard to soft spectral component in the 2–79 keV band), that agrees well with the observed correlation of hardness ratio and characteristic age that has been inferred from the spin-down rate, as reported by Enoto et al. (2010). They found a decrease in hardness ratio with increasing characteristic age or, likewise, an increase of F_h/F_s with magnetic field B.

Zoom In Zoom Out Reset image size

Furthermore, we investigate the phase-resolved spectra for NuSTAR and Swift and find that the PL indices were slightly harder for peak B than for peak A, while the breaking point for the energy is located around 10 keV for both pulses. We also examine pulsed spectra extracted from the NuSTAR observation, but due to low statistics we are not able to rule out a single PL. More high-energy data are needed in order to compare the results to those of Kuiper et al. (2006).

We comment further that our analysis did not find any spectral lines or absorption features in the phase-averaged and phase-resolved spectra that could be interpreted as cyclotron lines (Tiengo et al. 2013). We also searched for short bursts in the data, but in contrast to 1E 1048.1−5937, for which recently six bright X-ray bursts have been reported in NuSTAR data (An et al. 2014), we did not observe a similar behavior for 1E 2259+586.

We find that the phase-resolved spectra of 1E 2259+586 are consistent with the model of Beloborodov (2013a), although the available data do not allow us to derive strong constraints on the parameters of the model. In contrast to 1E 1841−045 (An et al. 2013), 4U 0142+61, and 1RXS J1708−4009 (Hascoët et al. 2014), a broad range of α_mag and β_obs is so far allowed for 1E 2259+586. The reason for this degeneracy is twofold: (1) the statistics of NuSTAR data of 1E 2259+586 are not as good as for the other objects; (2) the magnetic dipole moment of 1E 2259+586 is rather low, one order of magnitude smaller than for 1E 1841−045. As a result, the boundary of the radiative zone (where B ∼ 10¹³ G) is closer to the star. This leads to a larger allowed range of the j-bundle latitudes and to more flexible predictions of the model. Better statistics (with longer exposure times), extension of observations to the MeV range (where the emission of the e^± outflow peaks), or polarization measurements would help to break the degeneracy of α_mag and β_obs.

The footprint location of the j-bundle is better constrained, with 0.4 ≲ θ_j ≲ 0.75 and Δθ_j/θ_j ≲ 0.2 (at the 1σ level). This means that a ring-like footprint located at rather high colatitudes¹⁹ is statistically preferred to a polar cap footprint covering low colatitudes. This contrasts with the results obtained for 1E 1841−045 (An et al. 2013), 4U 0142+61, and 1RXS J1708−4009 (Hascoët et al. 2014), where a polar cap footprint provided a good fit. If confirmed with future, higher-statistics data, it would point to diversity in the crustal motions responsible for the twist of the magnetosphere.

We also analyzed the soft component, taking into account the extension of the outflow emission to low energies. The results are not very sensitive to the choice of parameters for the coronal outflow, because all good fits have similar extensions to low energies, which cut off below ∼5 keV. We find that the models BB+PL and BB+BB_tail provide equally good fits to the soft component; the fit with a modified blackbody (BB_tail) is only slightly worse but still acceptable, while the two-blackbody model (BB+BB) is statistically unacceptable. Even though BB+PL provides a good fit, it is physically problematic because the PL component is brighter than the BB (at all frequencies) and its energy content diverges at the low-energy end of the spectrum (in case of no cut off). It likely overpredicts the (sub) keV flux and the hydrogen column density N_H. Our physically preferred model is BB+BB_tail, where the cold BB corresponds to emission from most of the NS surface and the hot modified BB to a hot spot on the star (Hascoët et al. 2014). Such a spot may be expected at the footprint of the j-bundle, as some particles produced in the e^± discharge bombard the footprint and heat it (Beloborodov 2009). The area of the cold and hot thermal components are $\mathcal {A}_1 \approx 0.7$ $\mathcal {A}_{\rm NS}$ and $\mathcal {A}_2 \approx 0.02$ $\mathcal {A}_{\rm NS}$, respectively, where $\mathcal {A}_{\rm NS}$ represents the surface area of a NS with a typical radius of R_NS = 10 km. These results are similar to those obtained for 1E 1841−045, 4U 0142+61, and 1RXS J1708−4009 (Hascoët et al. 2014). Note however that "BB_tail only" also provides an acceptable fit for the soft component in 1E 2259+586, and the parameters of the hot spot in this source are not strongly constrained from these data.