CONFIRMATION OF A HIGH MAGNETIC FIELD IN GRO J1008−57

We fit two models frequently used in modeling accreting pulsar spectra to the phase-averaged continuum spectra: a power law with a high-energy cutoff, and an npex model consisting of two power laws with negative and positive spectral indices and an exponential cutoff (Makishima et al. 1999). We also included a Gaussian iron line and a low-energy blackbody component. For fits including data from Suzaku XIS, a second Gaussian component was needed to adequately fit the iron line complex. We used an updated version of the Wilms et al. (2000) absorption model (tbnew) as a neutral absorber with wilm abundances (Wilms et al. 2000) and vern cross-sections (Verner et al. 1996). For the power law with high-energy cutoff, we removed residuals due to the discontinuity at the cutoff energy with a Gaussian absorber tied to the cutoff energy (e.g., Coburn et al. 2002, and references therein). We allowed the normalizations to vary between all instruments being fit.

In contrast to the fits to Type I bursts reported by Kühnel et al. (2013), we found the npex model provides a better fit for all combinations of instruments despite having fewer free parameters, so we restrict our attention to this model for further analysis. Yamamoto et al. (2014) similarly found that the npex model provided the best fit to the Suzaku data from this giant outburst. Table 1 provides the best-fit values for the phase-averaged continuum parameters, and Figures 2–5 show the best fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Best-fit Phase-averaged npex Parameters for Fits with NuSTAR (N), NuSTAR and Swift (NX), Suzaku (S), and All Instruments (NXS)

Parameter	N	NX	S	NXS
NH (1022 cm−2)	$1.1^{+1.2}_{-1.1}$	1.06 ± 0.06	1.60 ± 0.03	1.45 ± 0.02
Γ1	0.15 ± 0.05	$0.27^{+0.06}_{-0.04}$	$0.94^{+0.02}_{-0.05}$	0.30 ± 0.01
A1 (10−1)	$1.76^{+0.26}_{-0.20}$	2.38 ± 0.15	$5.17^{+0.05}_{-0.15}$	$2.70^{+0.08}_{-0.10}$
Γ2	2	2	2	2
A2 (10−4)	$1.39^{+0.30}_{-0.18}$	$2.0^{+1.6}_{-0.4}$	$4.7^{+2.0}_{-0.4}$	1.81 ± 0.07
Efold (keV)	$8.44^{+0.14}_{-0.13}$	$8.34^{+0.38}_{-0.17}$	$7.47^{+0.16}_{-0.08}$	8.17 ± 0.05
E0, 1 (keV)	6.60 ± 0.02	6.60 ± 0.02	6.59 ± 0.02	6.58 ± 0.01
σ1 (keV)	0.29 ± 0.02	0.30 ± 0.02	0.28 ± 0.03	0.27 ± 0.02
Agauss, 1 (10−3)	6.8 ± 0.7	7.1 ± 0.4	$9.0^{+1.3}_{-1.4}$	6.5 ± 0.4
E0, 2 (keV)			$5.4^{+0.2}_{-0.1}$	$6.0^{+0.2}_{-0.3}$
σ2 (keV)			$0.07^{+0.03}_{-0.02}$	0.19 ± 0.02
Agauss, 2 (10−3)			$7.4^{+0.5}_{-0.2}$	$27^{+7}_{-5}$
kT (keV)	1.61 ± 0.03	$1.62^{+0.09}_{-0.07}$	$3.4^{+0.3}_{-0.2}$	0.42 ± 0.02
ABB (10−3)	14 ± 4	$8.7^{+1.7}_{-1.4}$	$79^{+5}_{-7}$	$3.9^{+0.3}_{-0.2}$
χ2/dof	632.6/536	881.8/692	1134.8/343	2424.9/1045
$\chi ^{2}_\mathrm{red}$	1.18	1.27	3.31	2.32

Notes. Errors are 90% C.L. In all fits the positive power-law index Γ₂ was free to vary but converged to the limiting value of 2. Normalization units are photons keV⁻¹ cm⁻² s⁻¹ for the power-law components (A_{1, 2}), photons cm⁻² s⁻¹ for the Gaussians (A_{gauss1, 2}), and $L_{39}/D_{10}^2$ for the blackbody (A_BB), where L₃₉ is the source luminosity in units of 10³⁹ ergs⁻¹ and D₁₀ is the distance to the source in units of 10 kpc.

Download table as:  ASCIITypeset image

The brightness of the source highlights systematic effects in joint fits between multiple instruments, producing poor goodness of fit. Only the NuSTAR-only fit has a reasonable goodness of fit, at $\chi ^2_\nu = 1.18$ for 536 degrees of freedom. There is substantial disagreement between the instruments below 10 keV (e.g., Figure 5). The NuSTAR and Suzaku observations are not simultaneous, and so some spectral evolution may have occurred between the two epochs. However, there is also disagreement even among the three XIS modules (Figure 4). This disagreement at low energies, driven primarily by XIS1, leads to differences in the best-fit blackbody temperature and power-law indices (Table 1). These discrepancies were also noted in this data set by Kühnel et al. (2013) and Yamamoto et al. (2014); these authors elected to excise several energy ranges from the XIS backside-illuminated detectors or exclude the data entirely. The fit for the blackbody temperature shows multiple minima in the Suzaku-only fit, for example, with the ∼3 keV temperature preferred to the 0.4 keV temperature suggested by the joint fit to all instruments. Similarly, our coarser binning and inclusion of all three XIS modules results in better fits to the iron line complex with two broadened Gaussians rather than the narrow 6.4, 6.67, and 7 keV lines fit by Kühnel et al. (2013) and Yamamoto et al. (2014). We are primarily interested in the spectral behavior at high energies, which is well above the folding energy E_fold and hence relatively insensitive to these parameters.

Next, we fit the data using the above continuum model and a multiplicative cyclotron scattering feature using a pseudo-Lorentzian optical depth profile (the XSPEC cyclabs model; Mihara et al. 1990). We initially confined our search to line centers above 50 keV, with fits initialized near the 75 keV value reported by Yamamoto et al. (2013). Table 2 reports the best-fit parameters, while Figures 2–5 compare the residuals to fits without a cyclotron line.

Table 2. Best-fit Phase-averaged npex Parameters with a Cyclotron Feature for Fits with Suzaku (S) and All Instruments (NXS)

Parameter	S	NXS
NH (1022 cm−2)	$1.62^{+0.04}_{-0.06}$	$1.43^{+0.03}_{-0.02}$
E0, cyc (keV)	$74^{+3}_{-2}$	$78^{+3}_{-2}$
cyclabs width (keV)	<5.89	$11^{+6}_{-4}$
cyclabs depth	$4.2^{+3.8}_{-3.3}$	$0.81^{+0.15}_{-0.13}$
Γ1	$1.03^{+0.12}_{-0.13}$	$0.29^{+0.03}_{-0.02}$
A1 (10−1)	$5.24^{+0.23}_{-0.29}$	$2.54^{+0.13}_{-0.17}$
Γ2	$1.10^{+0.04}_{-0.20}$	2
A2 (10−4)	$52.479^{+0.005}_{-7.573}$	$1.40^{+0.29}_{-0.19}$
Efold (keV)	8.9 ± 1.0	$8.55^{+0.31}_{-0.14}$
E0, 1 (keV)	$6.589^{+0.002}_{-0.010}$	$6.58^{+0.01}_{-1.44}$
σ1 (keV)	0.29 ± 0.03	0.27 ± 0.02
Agauss, 1 (10−3)	$9.3^{+1.3}_{-1.2}$	6.5 ± 0.4
E0, 2 (keV)	$5.36^{+0.17}_{-0.10}$	$5.6^{+0.4}_{-0.5}$
σ2 (keV)	$0.63^{+0.25}_{-0.14}$	$2.07^{+0.32}_{-0.27}$
Agauss, 2 (10−3)	$5.7^{+0.8}_{-1.9}$	$34^{+16}_{-9}$
kT (keV)	$2.90^{+0.99}_{-0.26}$	0.45 ± 0.02
ABB (10−3)	$43^{+31}_{-20}$	$4.2^{+0.5}_{-0.4}$
χ2/dof	1070.7/340	2266.0/1042
$\chi ^{2}_\mathrm{red}$	3.15	2.17

Note. Errors are 90% C.L. Normalizations units are as in Table 1.

Download table as:  ASCIITypeset image

Fits to the NuSTAR data alone do not provide strong constraints on the CRSF parameters, as there are degeneracies with the continuum modeling because the cyclotron line lies at the upper edge of the NuSTAR bandpass. However, there are clear residuals in the NuSTAR data above 70 keV, and NuSTAR-only fits are significantly improved by the CRSF (χ²/dof improves from 632.6/536 to 569.8/533). The best-fit Suzaku CRSF parameters are a reasonable match to those reported in Yamamoto et al. (2014) given the minor differences in analysis methods.

Combining the NuSTAR data with Suzaku provides an independent confirmation of the line. In the joint fit, the line centroid moves to 78 keV and the best-fit width is 11 keV, matching the values obtained by Yamamoto et al. (2014) in their npex fits including XIS. The GSO cross-normalization changes from 1.19 relative to NuSTAR in the npex fit to 1.38 in the npex with cyclotron feature fit. All other cross-normalization constants remain constant within errors (Table 3). The cyclotron-line fit thus produces correctly the expected agreement between the normalizations of Suzaku PIN and GSO.

Both the NuSTAR and Suzaku data thus independently show evidence for a CRSF in the 70–80 keV range. Because the NuSTAR data do not cover the entire CRSF, the joint fit provides the best constraint on the line parameters, but the parameters are sensitive to the NuSTAR and Suzaku-HXD cross-calibration.

We assessed the significance of the detections using the method of posterior predictive p-values (ppp-values; Protassov et al. 2002). Briefly, we simulate many instances of a model without a cyclabs feature by folding the spectral model through instrumental responses; the exposure and binning are matched to the real data, and the data and background are perturbed to account for counting statistics. For each simulated data set, we then fit the null model and a test model with a cyclabs feature. For each simulated realization, we determine Δχ² between the two models and compare the distribution of Δχ² values for the simulations to the observed value.

If few of the simulated Δχ² values are as large as observed in the real data, this provides evidence for the CRSF. Rather than restricting the simulated model parameters to those of the best-fit null model, we use the Cholesky method to draw the simulated parameters from a multivariate Gaussian derived from the covariance matrix obtained in the fit to the null model (Hurkett et al. 2008).

We performed 10,000 simulations of the NuSTAR data alone as well as the joint NuSTAR, Swift, and Suzaku data. The line energy was allowed to vary in the 50–100 keV range and the line width between 1 and 30 keV. In all cases, the simulated Δχ² was less than the value observed in the real data, providing >3.9σ evidence for the existence of the line in each of the two fits. In most of the simulated cases, the best-fit depth of the line is zero, and so the two models are indistinguishable. The largest deviation in χ² was 17.0 (21.3) for the combined data sets (NuSTAR only), far smaller than the Δχ² values of 278.5 (62.8) seen in the real data. Based on the difference between the observed and simulated Δχ² distributions, it is clear that the CRSF detection is much more significant in the joint fit than when using NuSTAR alone.

Given the distribution of Δχ² in these simulations, it would be computationally unfeasible to simulate enough realizations to expect a Δχ² value near the true value and obtain a true ppp-value significance. We can obtain a simple estimate of the significance (and hence the number of simulations required to obtain that chance deviation) by summing the data and model counts in the ±1σ energy window around the best-fit 78 keV cyclotron line. Dividing the difference between the npex model without the cyclotron line and the data in this region by the statistical error allows us to estimate the level of chance fluctuation needed. The deviations in the NuSTAR data (which do not cover the full cyclotron line) are 1.8σ and 2.3σ when the modules are considered independently; the deviation in the GSO data taken alone is 8.0σ. We thus expect the GSO measurement to dominate the fit. (We do not correct for trials over energy because the high-energy line was previously reported in other observations.) If taken at face value, the statistical errors would require more than 8 × 10¹⁴ simulations to achieve deviations in χ² comparable to the observed values.

We considered whether systematic calibration uncertainties could be responsible for the observed feature. While the method of ppp-values provides a robust assessment of line significance (Hurkett et al. 2008), it is sensitive to false positives if systematic errors are present. If a line feature is due to inaccuracy in the instrumental responses or the modeled background, ppp-value tests will confirm its statistical significance but not its physical reality. The calibration of the NuSTAR responses in the 70–78 keV range is less certain than at lower energies due to the increasing faintness of astrophysical calibrators. However, measured deviations from a fiducial spectrum of the Crab Nebula are <15% from 70–78 keV (Madsen et al. 2014). Similarly, few-percent deviations of the Crab spectrum have been measured in Suzaku GSO spectra near 70 keV (Yamada et al. 2011). These effects are not large enough to produce the ∼30% deviation seen here, so we conclude that the feature is both significant and real.

We searched for a cyclotron line at half the energy of the 78 keV line reported in Section 3.2. The NuSTAR data enable a more sensitive search than is possible with PIN: the combined NuSTAR data have a signal-to-noise ratio (S/N) of 135 keV⁻¹ at 40 keV in these data, compared to 60 for PIN. (Data from both instruments are strongly source-dominated, given the brightness of the outburst.) No obvious residuals are apparent in the phase-averaged NuSTAR data near ∼38 keV (Figure 2), consistent with previous non-detections in phase-averaged data. Some residual structure is present in the Suzaku-PIN data below 40 keV (Figure 4). Using PIN data, Yamamoto et al. (2014) reported a possible fundamental with $E_0 = 36.8^{+1.1}_{-0.7}$ keV, optical depth at line center of $0.06^{+0.08}_{-0.03}$, and width of $11.1^{+7.2}_{-10.2}$ keV in their Suzaku-only fit, but concluded it is not statistically significant. A double-cyclotron line NuSTAR–Suzaku–Swift joint fit with the fundamental restricted to half of the 90% error limits for the 78 keV does fit a line depth at 39.8$^{+0.5}_{-1.2}$ keV. It has depth 1.0$^{+0.7}_{-0.2}$ and width 6.0$^{+8.4}_{-4.5}$ keV. However, the improvement in χ² is modest, only 6.3 for three additional free parameters. A NuSTAR-only fit to a line at this position results in a fundamental with depth consistent with zero. The 90% CL upper limit on the optical depth at 39 keV is 0.04. The possible 39 keV fundamental fit by the Suzaku data is thus disfavored by the more sensitive NuSTAR data. A broader NuSTAR search from 34–40 keV (at half of the line centroid identified by the independent NuSTAR and Suzaku fits) similarly yielded line depths consistent with zero.

We also performed a generic search for lower-energy lines by stepping a cyclabs feature through a 2 keV grid of energies over the 10–60 keV range. We used the NuSTAR data only, as in the joint fit the residuals show dips (due to response differences highlighted by the brightness of the source, Figure 5(c)) not present in the NuSTAR-only fit (Figure 2(c)). For speed, the continuum parameters were frozen in the initial search. Only one trial (at 26 keV, with Δχ² of 5.4) fit a line depth greater than zero. In this case the best fit line width was 1 keV, at the narrowest limit, and there is a known response calibration feature at this energy, so we do not consider this a reliable detection. Over all energies, the largest 90% CL upper limit on the optical depth was 0.09 at 52 keV. Accordingly, we can rule out a lower-energy fundamental with greater depths in the phase-averaged data.

Because cyclotron line intensities and energies may vary with pulse phase (e.g., Fürst et al. 2013), we split the NuSTAR observation into phase bins of roughly constant S/N and conducted spectral fits on each.

We barycentered the NuSTAR event data with barycorr using the DE200 solar system reference frame (Standish 1982) and applied a correction for light-travel time across the binary orbit using the ephemeris of Kühnel et al. (2013). Figure 6 shows the NuSTAR data folded at the best-fit spin period of 93.57434 s.

Zoom In Zoom Out Reset image size

The phase-resolved data were well-fit by the npex spectral model. We fixed the positive power-law index to 2, consistent with the values obtained in the phase-averaged fits. The photon indices show a correlation with intensity (Figure 6), while the folding energy shows a mild secular increase throughout the pulse. No obvious deficits are present in the residuals at lower energies.

We attempted to observe phase evolution of the 80 keV CRSF. Using four phase bins, Yamamoto et al. (2014) found only a slight dependence of the CRSF energy on the pulse data using the Suzaku data. We froze the width and depth of the CRSF to the best-fit values from the joint NuSTAR–Suzaku–Swift fit, but left the energy free to vary. However, the relatively limited NuSTAR energy coverage of the line does not permit firm constraints on the line energy.

We also performed a harmonic CRSF fit to the phase-resolved data. We fit for a fundamental line in the 34–47 keV range and froze the ∼80 keV harmonic width and depth to the best-fit phase-averaged values. In all phases, the line depth was consistent with zero, the linewidth unconstrained, and the improvement in χ² < 2. The 90% C.L. line depth upper limits were in the range 0.09–0.3.

We repeated the generic grid search for low-energy CRSFs in the phase-resolved spectra in case phase-dependent intensity was pushing a fundamental below detectability in the phase-averaged fit. As in the phase-averaged case, the additional CRSF fit widths pegged at the minimum value of 1 keV, fit depths that were consistent with zero, and/or were associated with the small known response feature at 26 keV. Over all energies, the largest upper limits on the line depth were 0.2–0.4 and occurred in the 50–60 keV range.

Accordingly, our phase-resolved fits rule out a phase-dependent fundamental CRSF below 70 keV with depth greater than one third of the depth of the 80 keV CRSF.