© ESO 2018

1. Introduction

Cyclotron resonant scattering features (CRSFs, or cyclotron lines) in the hard X-ray spectra of accreting neutron stars allow for the measurement of the magnetic field strength close to the neutron star surface. CRSFs have been detected in approximately 25 sources to date (Caballero & Wilms 2012) and the measurements show that these neutron stars possess strong magnetic fields of the order of 10¹² G (e.g., Trümper et al. 1978). These strong fields dominate the accretion geometry, leading to the formation of accretion columns above the magnetic poles in which most of the hard X-ray radiation is produced. To understand the physics within these accretion columns, as well as the emission geometry of the X-rays, a thorough understanding of the magnetic field configuration is therefore important.

CRSFs show up as broad absorption features in the hard X-ray spectrum of accreting pulsars. They are formed by scattering of photons off electrons quantized onto Landau-levels in the strong magnetic field close to the surface of the neutron star. As the Landau levels are directly correlated with the magnetic field strength, so is the CRSF energy, and we can estimate the magnetic field strength in the line-forming region following the so-called 12-B-12 rule,

$$B_{12}=E_{\text{CRSF}}\times (1+z)/\mathrm{11.6.}$$(1)

Here, B ₁₂ is the magnetic field in 10¹² G, E _CRSF is the cyclotron line energy in keV, and z is the gravitational redshift due to the mass of the neutron star.

In recent years, more and more evidence is being collected of the variability of CRSF energy as function of luminosity (see, e.g., Fürst et al. 2014 and references therein). The direction of the correlation seems to depend on luminosity: at very high luminosities (L ≳ 3 × 10³⁷ erg s⁻¹), an anti-correlation between CRSF energy and luminosity is evident e.g., V 0332+53; Tsygankov et al. 2010), while at lower luminosities a positive correlation has been found (e.g., Her X-1; Staubert et al. 2007). Different theories for this behavior have been proposed; for example, the formation of a shock in the accretion column (Becker et al. 2012), or reflection on the neutron star surface (Poutanen et al. 2013).

At the very lowest luminosities, below ∼5 × 10³⁶ erg s⁻¹, the situation becomes even more complicated: sources located there do not show a uniform behavior. For example, the CRSF energy of A 0535+26 seems to be constant as a function of luminosity (Caballero et al. 2007; Ballhausen et al. 2017), while in Vela X-1 a clear positive trend is observed (Fürst et al. 2014). Theoretical work explaining the different behavior at these luminosities is currently still missing.

While the CRSFs provide useful information about the magnetic field strength, measurements of the magnetic field geometry are not straight forward. To first order, the magnetic field can be assumed to be a simple symmetric dipole. This approach explains, for example, the apparent variation of the CRSF energy with pulse phase in GX 301−2, assuming that at different phases different magnetic field strengths are sampled (Suchy et al. 2012). However, it fails to explain the dependence of the CRSF energy in Cen X-3 (Suchy et al. 2008), for which a dipole magnetic field, off-set from the neutron star center, might provide a better explanation (Kraus et al. 1996).

Decomposition of pulse-profiles also often results in the requirement of an asymmetric dipole field (e.g., Sasaki et al. 2012 and references therein). In this approach, the energy-dependent pulse profile is modeled with phase-dependent contributions from the two accretion poles, under the influence of strong light-bending of the emitted photons by the gravitational field of the neutron star. For a good description of the observed pulse profile, this decomposition typically requires that the two poles not be exactly opposite of each other, but slightly offset (e.g., ∼25° in A 0535+26; Caballero et al. 2011).

In a detailed study of 4U 0115+63, the source known to have the most CRSFs (up to 5; Santangelo et al. 1999; Heindl et al. 1999), Iyer et al. (2015) find evidence that the lines can be separated into two pairs, with different fundamental energies. They argue that this separation could be explained if accretion happens on both poles, but each pole has a different magnetic field strength, due to an off-set magnetic field from the neutron star center.

Even more complicated magnetic field geometries are sometimes proposed for the recently discovered ultra-luminous pulsars (e.g., Bachetti et al. 2014). While the spin-up rate of NGC 7793 P13 seems to agree with a dipolar field of 1.5 × 10¹² G (Fürst et al. 2016), Israel et al. (2017) need to invoke a strong multi-polar component to explain the extreme luminosity of NGC 5907 ULX1.

To study the magnetic field geometry as well as the luminosity dependence of the CRSF energy at low luminosities, the high-mass X-ray binary (HMXB) GX 301−2 is an ideal target. It has a high pulsed fraction, a well-studied, strong CRSF and, due to it being only about 3 kpc away, also has a high flux at low luminosities. It was discovered in 1969 by balloon-borne experiments (Lewin et al. 1971; McClintock et al. 1971). Its pulse period is around ∼683 s (White et al. 1976) and the system has an orbital period of ∼41.5 d (Koh et al. 1997; Doroshenko et al. 2010) with a relatively large eccentricity of e = 0.47. This eccentricity gives rise to a regular pre-periastron flare, at which the neutron star overtakes the accretion stream (Leahy 2002; Leahy & Kostka 2008). Therefore, the average flux is strongly dependent on orbital-phase, varying in the range 2 − 20 × 10⁻⁹ erg cm⁻²s⁻¹ (Fig. 1). However, strong variability is present on top of that average luminosity, in particular outside the pre-periastron flare where accretion is dominated by direct accretion from the stellar wind. As an example, Fig. 1 shows the fluxes of two archival Suzaku observations, which differ significantly from the expected average flux at their respective orbital phases.

Fig. 1. Average light curve from the Burst Alert Telescope (BAT) of Swift, folded on the orbital period where orbital phase 1 is periastron. A strong flare shortly before periastron is clearly visible. We note that due to the eccentric orbit, the phases are not equidistant in time. The two NuSTAR observations are shown by the red crosses and the Suzaku observations by the blue diamonds. These data use the right-hand y-scale with the flux measured between 5 and 50 keV. The left and right y-scales are aligned based on the daily BAT count rate measured during the NuSTAR observations. We note the strong variability above and below the average flux in the pointed observations of NuSTAR and Suzaku.

Mihara et al. (1995) discovered a CRSF in the hard X-ray spectrum of GX 301−2 using data from the Ginga satellite (Makino & ASTRO-C Team 1987). They describe a very broad (σ ≈ 16 keV) and strong line around 35 keV. The line was later confirmed by Kreykenbohm (2004), using Rossi X-ray Timing Explorer (RXTE; Bradt et al. 1993) data. These authors also found that the line energy varies between ∼30 and 40 keV as a function of pulse phase. The variation was studied in more detail by Suchy et al. (2012), who explain the smooth sinusoidal variation of the energy as a function of pulse phase as being due to different viewing angles on a dipolar magnetic field. La Barbera et al. (2005) analyzed four BeppoSAX observations, and also found a CRSF, albeit at much higher energies, at around 45–53 keV, depending on orbital phase. These energies are consistent with the values reported by Doroshenko et al. (2010), who analyzed INTEGRAL data taken during the pre-periastron flare. We note that in all studies, the CRSF is very strong and removes a significant fraction of the continuum flux. Accurate modeling therefore depends on obtaining the correct continuum shape, and also requires good coverage of the energies above and below those dominated by the CRSF.

The large width of the CRSF also resulted in relatively large uncertainties of the centroid energy of the line, making a detailed study of the behavior of the line energy as function of luminosity difficult. Suchy et al. (2012) see an indication for an anti-correlation between energy and luminosity, but the trend is not significant at the 90% level. Combining the results by Suchy et al. (2012) with the measurements by La Barbera et al. (2005), a positive correlation seems to be present but it is very likely that these values are dominated by systematic differences in the analyses.

In this paper we analyze two observations of GX 301−2 taken with the Nuclear Spectroscopic Telescope Array (NuSTAR; Harrison et al. 2013). The NuSTAR data provide unprecedented energy resolution and signal-to-noise ratios (S/N) above 10 keV and are therefore ideally suited to studying the CRSF region. In Sect. 2 we describe the data reduction and calibration. In Sect. 3 we present phase-averaged spectroscopy and compare the spectrum to archival Suzaku data. In Sect. 4 we discuss the pulse profiles and perform phase-resolved spectroscopy. We discuss a physical interpretation of the results and draw conclusions in Sect. 5.

2. Data analysis and calibration

The two NuSTAR observations are separated by roughly one year and were performed at different orbital phases (see Table 1). The first observation was performed just after apastron, during an expected intermediate flux state of the source. However, by chance the observations took place at a much lower luminosity than expected. The second observation was obtained during a later orbital phase and closer to the pre-periastron flare, which resulted in the desired higher fluxes (Fig. 1). As far as we know, these data have not been published elsewhere yet.

NuSTAR has two identical instruments, Focal Plane Module A and B (FPMA and FPMB, respectively). We extract data from both modules separately, using the standard pipeline nustardas v1.7.1, as distributed with HEASOFT v6.20 and calibration files 20170222. Source data from both observations were extracted from a circular region with a radius of 40′′, centered on the brightest spot in the respective sky images of each module. The extraction region size was chosen to optimize the S/N at energies above 45 keV. A larger region introduces more background photons than source photons at the highest energies, due to the shape of the point spread function (PSF). The influence of the smaller region on the S/N at low energies is negligible. Background data were extracted from a circular region with a radius of 140′′, located at the opposite end of the field-of-view to avoid contamination from source photons. The spectra were rebinned following the energy resolution of the detectors, by combining 2, 3, 5, 8, 16, 18, 48, and 72 raw channels for energies above 3, 10, 15, 20, 35, 45, 55, and 65 keV, respectively. Additionally we required at least 20 counts per bin and modeled the data in the 3–78 keV energy range.

Suzaku observed GX 301−2 twice, as described by Suchy et al. (2012). Here we reanalyze the longer of the two observations, ObsID 4030444020, which was taken on 2009 January 05 (Table 1) and took place at post-periastron phase (0.38). The data were obtained with the X-ray Imaging Spectrometer (XIS; Koyama et al. 2007) with CCD cameras 0, 1, and 3 operating in normal mode in a 1/4 window sub-array option, along with data from the PIN instrument from the Hard X-ray Detector (HXD; Takahashi et al. 2007). The exposure times for each instrument were 61.8 ks for XIS and 51.0 ks for PIN (Table 1).

The Suzaku data were reprocessed and extracted using the standard aepipeline as distributed with HEASOFT v6.19 and calibration packages HXD-20110913, and XIS-20160607, and the screened XIS events were filtered to exclude times of telemetry saturation. The XIS images were first extracted with xselect, and then further corrected for thermal attitude uncertainties with aeattcor2. For both 3 × 3 and 5 × 5 editing modes, the source data were extracted from circular regions with 120′′ radii, with 20′′ pile-up exclusion radii centered on the PSFs. The background regions, circular with 120′′ radii, were extracted at the furthest locations from the point source within the windows.¹ We generated the XIS source and background event files, images, and spectra using xselect. The energy and ancillary responses were generated with xisrmfgen and xissimarfgen.

Using hxdpinxbpi, we extracted the dead-time corrected PIN source spectrum and the total background spectrum and applied the ae_hxd_pinhxnome5_20080716.rsp PIN response file for the specific calibration epoch. The total PIN background spectrum is composed of the sum of cosmic X-ray background (CXB; ∼5% of the PIN background flux) and non X-ray modeled background (NXB; ∼95% of the PIN background flux).

Data analysis was performed with the Interactive Spectral Interpretation System (ISIS; Houck et al. 2000) v1.6.2. All uncertainties are reported at the 90% confidence level unless otherwise noted.

3. Phase-averaged spectroscopy

3.1. NuSTAR observation 2

As the second observation (ObsID 30101042002) was the brightest observation overall with a higher S/N, we begin our spectral modeling with these data and later apply our preferred models to the other NuSTAR observation, as well as to archival Suzaku data. The spectrum of GX 301−2 has been well studied in X-rays so we can start with a relatively sophisticated spectral model, including a partially covering absorber (with column density N_H and covering fraction f). We model the absorption with an updated version of the tbabs model (Wilms et al. 2000), using the wilm abundance vector and cross-sections presented by Verner et al. (1996). We add Fe Kα and Fe Kβ fluorescence lines (with line center energy E_Kα/β, width σ_Kα/β, and equivalent width EW_Kα/β). However, due to the proximity to the Fe K-edge, we cannot constrain the width of the Fe Kβ line, and we therefore fix it to be narrow (σ_Kβ = 10⁻⁶ keV). We apply different phenomenological models to describe the continuum emission, as no physical model was available that describes this emission at the relatively low luminosity of GX 301−2 during our observations (∼2 × 10³⁶ erg s⁻¹).

To model the previously observed CRSF in the hard X-ray spectrum, we use a multiplicative absorption model with a Gaussian optical depth profile (model gabs in XSPEC, described through its energy E_CRSF, width σ_CRSF, and strength d_CRSF²). From previous investigations at similar luminosities (e.g., Suchy et al. 2012) we expect the CRSF energy to be around 35 keV, and we initialize the model with this energy. We also add a multiplicative constant to the model for FPMB, to allow for slight differences in absolute flux calibration with respect to FPMA (CC_B).

This baseline model can thus be written as

$$\begin{array}{c}\hfill \text{CC}\times (f\times N_{\mathrm{H}}+[1-f])\times (\text{CONT}\times {\text{CRSF}}_1+\mathrm{Fe}\mathrm{K}\alpha +\mathrm{Fe}\mathrm{K}\beta )\end{array}$$(2)

where CONT is the respective continuum model.

Consistent with the very high absorption column and prominent iron lines, GX 301−2 is known to also show a Compton shoulder at the red wing of the Fe Kα line (Watanabe et al. 2003; Fürst et al. 2011). However, this shoulder extends to only about 6.2 keV and is within the energy resolution of NuSTAR, completely subsumed by the Gaussian model for the Fe Kα line. We therefore do not include this shoulder in our spectral modeling. Small residuals at the red end of the line might be related to this shoulder, or to the unidentified absorption feature found around 5.6 keV in XMM-Newton data by Fürst et al. (2011).

We apply different phenomenological models typically used to describe the continuum of magnetized, accreting neutrons stars: cutoffpl, FDcut (Tanaka et al. 1986), NPEX (Mihara et al. 1998), and HIGHECUT. For the latter, we have smoothed the artificial kink at the cutoff energy with an additional multiplicative Gaussian line. The continuum parameters are, where applicable, the flux ℱ between 5–50 keV, the cutoff energy E_cut, the folding energy E_fold, the photon index Γ, and the normalization of the second component A₂ in the NPEX model.

None of the continuum models provide an acceptable fit, see Table A.3 in the appendix, where we also give the best-fit parameters of all tested models (including a gain-shift component as described below). We found that the best fit is obtained by using the NPEX model, resulting in ${\chi }_{\mathrm{red}}^2=1.27$ for 457 degrees of freedom (d.o.f). Still, as shown in Fig. 2(b) even this model leaves very wavy residuals between 20 and 50 keV and requires a very broad CRSF (σ > 10 keV). The wavy residuals may indicate that the CRSF shape is not captured correctly by the one gabs component (cf. Cep X-4, Fürst et al. 2015) or that a second line is present.

Fig. 2. Spectrum of NuSTAR observation 2. Panel a: data from NuSTAR/FPMA (FPMB not shown for clarity), together with the best-fit NPEX model (black). We also show the two components of the NPEX model separately, indicating the soft component with a dashed line and the harder component (for which the photon-index is fixed to −2) with a dotted line. Additionally, we show the model with each CRSF removed in turn (orange and green) and with both removed (magenta). Panel b: residuals in terms of χ for the best-fit model with only 1 CRSF. Panel c: residuals to the best-fit model with 2 CRSFs Panel d: residuals in terms of data-to-model ratio without the 50 keV CRSF. Panel e: residuals without the 35 keV CRSF. Panel f: residuals without both CRSFs. In this panel the non-symmetric shape of the residuals can be clearly seen.

We have therefore added a second multiplicative Gaussian line, CRSF₂, also described by a gabs model. This approach results in a significantly improved fit with ${\chi }_{\mathrm{red}}^2=1.05$ for 454 d.o.f (i.e., Δχ² = 103 for three additional parameters). We find a strong line around 50 keV, as well as a much shallower feature around 35 keV.

However, in this fit small discrepancies between the FPMA and FPMB spectra become visible, beyond the small normalization difference of ∼3%. This discrepancy appears mainly as a small energy off-set between the two modules, in particular around the Fe Kα line. We therefore add an ad-hoc gain shift (GS) to the FPMA energies, for which we find a best-fit intercept value of ∼20 eV, that is, half a bin of the native NuSTAR binning. We note that this is a constant, energy-independent shift applied to all model channels within ISIS. This small shift is only necessary due to the extremely high statistics in our data and is about a factor of 20 smaller than the full width at half-maximum (FWHM) energy resolution of the detectors (Harrison et al. 2013). It is also below the typical absolute uncertainties of the NuSTAR gain, which are on the order of 40 eV in the iron line region (Madsen et al. 2015). We note that the applied gain shift provides limits on how well the line energies can be reconstructed on an absolute scale, even though the statistical uncertainties on, for example, the Fe Kα energy are much smaller.

While the NPEX continuum provided the statistically best-fit overall, using two independent gabs components to model the complex region between 30 and 60 keV resulted in a statistically good fit for the FDcut and HIGHECUT models as well, but not for the cutoffpl model (see Table A.4). This is in agreement with previous hard X-ray studies of GX 301−2, which also found a good description of the data only with more complex continuum models (e.g., La Barbera et al. 2005; Suchy et al. 2012). For all models with statistically acceptable fits the parameters of the two CRSF lines agree: the stronger one is located at ∼50 keV, while the weaker one is located at ∼35 keV. We note, however, that in the FDcut model, we had to fix the strength of the high-energy CRSF at the upper limit of 30, as it was unconstrained in the fit. We list the best-fit parameters for all models in Table 2.

As can be seen in Fig. 2 the CRSFs show significant overlap, due to their relatively large widths. They remove over half of the continuum flux between 25 and 70 keV, raising the question of how realistic the estimated continuum is. We therefore tried to fit other continuum models without CRSFs, including a broken power-law model as well as a comptt model. However, we do not find another statistically acceptable description of the very sharp turnover around 25 keV, together with the relatively straight continuum above that energy, modulated by wavy structure and the hardening above 60 keV.

The relative width σ/E of the CRSF is in line with what is expected from theoretical models and thermal line broadening (Schwarm et al. 2017), with values around 15% for both lines. The similar relative width of both lines is consistent with them being formed in the same plasma and gives us an indication of their physical origin.

We also checked whether or not the CRSF is better described by a single line with a Lorentzian optical depth profile (using the CYCLABS model in XSPEC), however, we do not find an acceptable fit with only one line. This finding can be understood due to the fact that the lines are broad, and the difference between a Gaussian and Lorentzian optical depth profile is negligible in this regime.

While the spectrum of GX 301−2 is rather complex at soft energies, due to the strong absorption and large contribution from fluorescence lines, no such features, besides the CRSF, are known above 10 keV. That is, the shape of the hard X-ray continuum, caused by Comptonization of seed photons in the hot electron gas of the accretion column, is expected to be very smooth. Any deviation from such a smooth continuum should therefore only be caused by effects from the magnetic field, i.e., CRSFs.

To test the assumption that the continuum above 10 keV is smooth, we applied the compmag model (Farinelli et al. 2012). This model implements a physical model of the accretion column based on theoretical calculations by Becker & Wolff (2007). With this physical continuum model alone, we do not find a good fit. However, when adding a thermal component in the form of a multi-temperature black-body (diskbb), we find a very good description of the data, if, and only if, we also include two CRSFs (${\chi }_{\mathrm{red}}^2=1.02$ for 451 d.o.f). The parameters of the CRSFs are consistent with the ones we find using the NPEX model (see Table A.4). Using only one CRSF does not lead to an acceptable fit (Table A.3). The parameters we find for the compmag model produce a continuum that is very smoothly rolling over at high energies, which is very closely mimicked by the NPEX model.

To test the significance of the detection of a second CRSF we run Monte Carlo simulations. We simulate 10 000 spectra, based on our best-fit with only one CRSF using the NPEX continuum (for the sake of computing speed, we do not use the physical compmag model; however, as NPEX produced an almost identical continuum shape, the results will also hold up for the physical modeling). For each simulated spectrum, we draw the model parameters from a Gaussian distribution around the best-fit values. We choose the width of the distribution to match the 90% uncertainties of the spectral parameters and apply Poisson statistics before folding it through the NuSTAR responses. We then apply the same binning schemes as for the real data and fit the simulated spectra with an NPEX continuum model, once with one CRSF and then again with two CRSFs. We compare the χ ² between these two models, which allows us to sample how likely it is to find an improvement in terms of χ ² as large as in the real data, which translates to a false-positive rate of an improvement of that strength. We find that the models show an improvement as large as or larger than the real data only for 1 of 10 000 simulated spectra. This translates to a false-positive rate of ≤0.01% and means that the feature in the real data is statistically significant at ≥99.9%.

Another test of significance is to apply the Akaike Information Criterion (AIC; Akaike 1974), which takes the degrees of freedom in each model into account. In particular we use ΔAIC, the difference in AIC between the fit with one CRSF and the one with two CRSFs. From that we find a false-positive rate of only 10⁻¹⁹. While this is much lower than inferred from the simulations, the AIC probability at these extreme ends of the distribution should be treated with some caution. Furthermore, we expect the false-positive rate in our simulations to decrease further if running more simulated spectra; however, the required numbers are computationally prohibitive.

La Barbera et al. (2005) used BeppoSAX data to perform a detailed study of the GX 301−2, and did not find evidence for a secondary line in BeppoSAX data in apparent contradiction to our results. In fact, they do not detect any significant line during their observation at a similar orbital phase (0.62–0.65) to that of the NuSTAR data presented here. However, the authors state that this non-detection is likely due to a lack of signal in the data.

To check if BeppoSAX could have picked up on the second line during a theoretical observation in a similar state as the NuSTAR observation with high enough S/N, we simulated 10 000 BeppoSAX spectra of the MECS, HPGSPC, and PDS instruments, based on responses and background data as published in the CALDB from May 1999. We use exposure times of 100 ks in MECS and 50 ks in HPGSPC and PDS, the latter being lower due to the rocking of the instruments. We use our best-fit NPEX model with two CRSFs as input, again drawing randomly from a Gaussian distribution around the best-fit values. We then compare the improvement between a model with two lines and a model with one line, based on the AIC. We find that in about 50% of all simulated spectra the fit is significantly improved (ΔAIC > 5) when using two lines. This result shows that the second line is at the edge of the BeppoSAX sensitivity, and is consistent with the non-detection in the existing data, which show a lower S/N than our simulated spectra (which also implicitly assume perfect knowledge of the responses). In comparison to INTEGRAL and Suzaku, BeppoSAX had the best detectors to measure the line profile, given its continuous coverage and good energy resolution. We can therefore confidently say that the non-detection of a secondary line in both Suzaku and INTEGRAL can be ascribed to the lack of data quality in these instruments (and see Sect. 3.3 for a comparison to archival Suzaku data).

The NPEX model essentially consists of two cutoffpls with the same folding energy, but with one component having a positive and the other a negative photon-index. We froze the index of the negative component to Γ = −2 (we note that in the definition of cutoffpl ℱ ∝ E ^−Γ ). Therefore, one of these components dominates the hard energies while the other component dominates the soft energies. We calculate the flux for each of these components separately, which gives us another indicator for the spectral hardness. These two components are plotted in Fig. 2a. We refer to these fluxes as “hard” and “soft”, respectively.

3.2. NuSTAR observation 1

Having established a good description of the spectrum of observation 2 in the previous section, we use the same models for the fainter data of observation 1. We use the same binning scheme as for observation 2 and find again an acceptable fit only when using two independent CRSFs. The statistical quality of the fit for the NPEX model is worse than in observation 2 and not significantly better than the FDcut or highecut models (Table A.2). All three models resulted in statistically acceptable fits, and we find that the CRSF parameters are independent of the chosen continuum model (Table 2). The gain-shift intercept is larger in observation 1 than in observation 2, with a best-fit value around 36 eV, which is still below the uncertainty of the energy calibration of NuSTAR.

We also tried to model this observation with one CRSF only, but find that the improvement when adding a second line is statistically significant for all continua (Table A.1). The smallest improvement is found for the highecut model, which still gives us a ΔAIC value of 12.3, which corresponds to a significance of > 99%.

We find that the parameters of the CRSFs are consistent with the parameters found in observation 2, in particular with respect to the energy, and therefore we do not find evidence for a luminosity dependence of the CRSF energy. However, as the flux is only about 30% lower compared to observation 2, any correlation is likely hidden within the uncertainties.

In Fig. 3 we show the spectra of both observations. The increased flux of observation 2 mainly occurs below ∼40 keV, while above, the spectra are very similar. This is reflected in the similar spectral parameters we find in both observations, where the only significant difference is the absorption column and the covering fraction.

Fig. 3. As in Fig. 2, but for observation 1. Panel a: additionally to the NuSTAR/FPMA data and best-fit NPEX model, we show the spectrum of observation 2 in gray for comparison. Panel b: residuals in terms of χ for the best-fit model with two CRSFs. Panel c: residuals in terms of data-to-model-ratio without the 50 keV CRSF. Panel d: residuals without the 35 keV CRSF. Panel e: residuals without both CRSFs.

3.3. Comparison to Suzaku

Suzaku observed GX 301−2 twice, as described by Suchy et al. (2012). Here we re-analyze the longer of the two observations, ObsID 4030444020, taken on 2009 January 05 (Table 1). Suchy et al. (2012) found a very good description of the spectrum with just one CRSF at around 35 keV and our goal is to see if the best-fit model we find for NuSTAR, which requires two CRSFs, is also able to describe the Suzaku data. As Suchy et al. (2012) discuss, the spectral parameters are somewhat dependent on the combination of XIS detectors used, due to calibration difference between the front-illuminated chips (XIS0 and XIS3) and back-illuminated chip (XIS1). As our main focus is on the hard X-ray spectrum, we only use data from XIS1 and PIN here for simplicity.

We re-binned the XIS data within ISIS following the scheme described by Nowak et al. (2011), which follows the energy resolution of the detector, requiring at least a S/N of 8 per bin. We restricted the energy range of XIS to be between 3 and 10 keV to reduce the influence of the very complicated soft X-ray spectrum on the statistical quality of the fit. We re-binned the PIN data to an S/N of 5 per bin and used them between 15 and 80 keV. Fluxes are given relative to PIN and we added a cross-calibration constant for the XIS data to account for differences in the absolute flux calibration between these detectors.

To model the data, we start with the best-fit NPEX model for NuSTAR observation 2. As the S/N of PIN above 50 keV is very low, we chose to fix the high-energy CRSF to the parameters found in NuSTAR. Due to the much higher spectral resolution of XIS compared to NuSTAR, we additionally added a Ni Kα line at around 7.5 keV and a Ca Kα line at 3.7 keV.

We find that this model describes the data very well, with similar parameters for the low energy CRSF to those we found in the NuSTAR data (Table 2). However, the strength of the CRSF is somewhat stronger, and the feature is removing significant flux from the implied continuum (Fig. 4). From the Suzaku data alone our model therefore might appear somewhat unphysical, as the CRSFs dominate the spectral shape above ∼20 keV. However, as we have seen that this model is clearly required by the NuSTAR data, the fact that we also find an acceptable description to the much brighter Suzaku data is a further indication that we use the correct description. We note that the line at 35 keV found by Suchy et al. (2012) has very similar parameters to the line we find here, indicating that in their model the 50 keV feature was adequately modeled by the different continuum used.

Fig. 4. Same as Fig. 2 but for the Suzaku spectrum. Panel a: XIS 1 data are shown in red, PIN data are shown in blue. Panel b: residuals in terms of χ for the best-fit model with 2 CRSFs. Panel c: residuals in terms of data-to-model-ratio without the 50 keV CRSF. Panel d: residuals without the 35 keV CRSF. Panel e: residuals without both CRSFs.

4. Phase-resolved spectroscopy

4.1. Pulse profiles

For the timing and pulse-phase-resolved analysis, we shifted all time information to the solar barycenter using the task barycorr. We also corrected for the binary motion, using the ephemeris by Koh et al. (1997). To find the pulse period, we folded the events and light curve on test periods in the vicinity of the known 685 s period and searched for the folded profile with the largest deviation from a constant (i.e., using the epoch-folding technique, Leahy et al. 1983). To determine the uncertainty on the period, we simulated 1000 light curves based on the pulse profile and performed the same search. The uncertainty is the 1σ interval of the distribution of the resulting periods. To include the strong pulse-to-pulse variations of GX 301−2, we added a noise term to these simulated light curves to obtain the same variance as in the original light curve (Davies 1990). We find P_{Obs. 1} = 686.61 ± 0.03 s and P_{Obs. 2} = 685.83 ± 0.05 s, indicating significant torque transfer onto the neutron star between these two epochs. The measured periods are in good agreement with those measured at the respective times by Fermi/GBM³ (Finger et al. 2009).

The pulse profiles of observation 2 in three energy bands are shown in Fig. 5, together with the pairwise hardness ratios of these profiles. A very strong energy dependence is evident, that is, the main peak (around phase ϕ = 0.75) becomes significantly narrower at higher energies, while the relative strength of the secondary peak (around ϕ = 0.2) increases. This evolution agrees well with the pulse profiles observed by RXTE, BeppoSAX, and Suzaku (Kreykenbohm 2004; La Barbera et al. 2005; Suchy et al. 2012). The energy dependence is strongest at low energies, while above ∼12 keV it is much less pronounced. We therefore chose to combine the high-energy data into one profile to increase the S/N.

Fig. 5. Panel a: pulse profiles of observation 2 in three different energy bands: 3–6 keV (blue), 6–12 keV (yellow), and 12–78 keV (red). The profiles were normalized by subtracting the respective mean rate and dividing by the respective standard deviation to allow for easier comparison of the shape between energy bands. The lower panels show the hardness ratio between each pair of these pulse profiles. The hardness ratio is calculated as (H − S)/(H + S), where H is the profile in the harder band, and S the profile in the softer band. Panel b: hardness ratio between the 3–6 keV and 6–12 keV profiles. Panel c: hardness ratio between the 6–12 keV and 12–78 keV profiles. Panel d: hardness ratio between the 3–6 keV and 12–78 keV profiles. The dashed lines indicate the phases selected for phase-resolved spectroscopy.

To highlight the energy dependence of the pulse profile we also show a color-coded map of the energy-resolved profiles of observation 2 (Fig. 6). The energy bins were chosen so that each band contains 7500 photons. From this binning it follows that the highest energy resolution is obtained around the Fe Kα line at 6.4 keV, where the highest count rates are measured. The different behavior as a function of energy between the two main peaks can be clearly seen, with the secondary one (around phase ϕ = 0.2) being much weaker at low energies, but becoming similarly strong to the primary one at higher energies.

Fig. 6. Color-coded energy resolved pulse profile map. The x-axis indicates phase (repeated once for clarity) and the y-axis energy in keV. The color indicates the relative strength in each energy band according to the scale on the right, normalized by (r − 〈r〉)/σ(r), where r is the rate in each bin, < r> is the average count-rate over each energy band, and σ(r) is the corresponding standard deviation.

The energy-resolved pulse profiles of observation 1 (not shown) appear to be qualitatively the same; the energy dependence in particular is very similar. The gap between the primary and secondary peaks is somewhat narrower in observation 1, but still clearly discernible.

To perform a phase-resolved analysis, we split the data into seven phase bins, indicated in Fig. 5. These bins were selected to roughly cover intervals of constant hardness, and in particular to provide a very good S/N across the two peaks of the pulse profile where the spectrum seems to change only marginally as a function of phase. We label the phase bins with letters A–G, where phase-bin B corresponds to the first peak and phase-bin F corresponds to the second peak.

4.2. Phase-resolved spectral analysis

To describe the spectra of the seven selected phase bins, we use the NPEX model, including two CRSF lines and the iron lines, as well as partially covering absorption. This is the same model as was found to describe the phase-averaged spectrum best (see Table 2). However, due to the reduced S/N in each spectrum, we had to fix the width σ of the CRSF lines, as well as the width and energy of the Fe Kα and Fe Kβ lines to the best-fit values of the phase-averaged spectrum. After preliminary fits, we also decided to fix the folding energy E_fold to the phase-average value, as we did not find a significant change of this parameter as a function of phase and found that it has some degeneracy with the photon-index Γ. Therefore, we artificially describe all changes in the continuum shape by variations of Γ and the ratio of the hard and soft continuum component.

4.2.1. Observation 2

The results of the phase-resolved fits are shown in Fig. 7 for observation 2. As can be seen in the figure there is strong variability in all parameters. The flux of the hard and soft power-law follow very well the corresponding soft and hard pulse profiles (see Fig. 5). In particular, the second peak in the hard flux (bin F) is much narrower than in the softer band.

On top of this variation, the photon index also changes significantly, varying from around 0.6 in bin B to almost 1.2 in bin D (i.e., the deep minimum in the pulse profile). It is interesting to note that the hardness ratios do not take on extreme values during phase-bin D, indicating a complex interplay between the photon-index and the fluxes of the two components of the NPEX model.

Fig. 7. Panel a: flux between 3 and 40 keV in keV s−1 cm−2. Shown in black is the “soft” component, in red the “hard” component (see Sect. 3 for a description of the components). Panel b: covering fraction. Panel c: absorption column in 1022 cm−2. Panel d: photon index. Panel e: energy of CRSF1. Panel f: strength of CRSF1. Panel g: energy of CRSF2. Panel h: strength of CRSF2. Panel i: iron line flux in ph s−1 cm−2.

Fig. 8. Same as Fig. 7, but for observation 1.

The CRSF around 35 keV (CRSF₂) varies significantly in energy as well as in strength as a function of pulse phase. On the other hand all energies of CRSF₁ (around 50 keV) are consistent with a constant, except the measurement in phase bin D. As this bin is taken during the pulse minimum, we do not have data above 50 keV, and therefore the lower CRSF energy might be a systematic effect. To explore this further, we froze the energy of the CRSF in this phase bin to the value from the phase-averaged fit, i.e., 49.6 keV. With this, we find a statistically slightly worse fit of Δ χ ² = 11, while all other parameters did not change significantly.

The variation of the parameters of the low-energy CRSF (energy and strength) are only single peaked, in contrast to the double-peaked pulse profile. In particular, the extrema do not line up with any of the peaks in the pulse profile, but seem to be shifted to later phases.

4.2.2. Observation 1

We performed the same phase-resolved analysis for observation 1. As the energy resolved pulse-profiles are very similar between the observations, we used the same definition of phase bins. We also applied the same model. However, due to the overall lower statistics in observation 1, we additionally had to fix the energy of the high-energy CRSF to the phase-averaged value of 50.6 keV. As discussed for observation 2, the energy of this CRSF does not vary significantly as a function of phase, with a possible exception being during the pulse profile minimum.

We find a similar behavior of the spectral parameters as a function of phase as in observation 2, albeit at a much reduced statistical quality. In particular, E_CRSF,2 shows the same single peaked variation, with the extrema phase shifted with respect to the main peaks in the pulse profiles. The variation of the photon index seems to be somewhat different, as phase-bin A now shows a relatively low Γ instead of a high one. This might be related to the fact that in this bin the hard power-law flux is lower than that of the soft component, while the opposite is true for this bin in observation 2. This behavior might be due to a slight model intrinsic degeneracy between the photon-index Γ and the flux ratio of the two NPEX components. However, we tested for such a degeneracy in the phase-averaged spectrum by calculating a confidence contour between Γ and the flux ratio and did not find any indications for its existence.

5. Discussion and conclusion

We have presented an in-depth analysis of the hard X-ray spectrum of GX 301−2, as observed with NuSTAR in 2014 and 2015. We find evidence for two CRSFs, one at 35 keV and the other at 50 keV. Both features are broad and overlap strongly. However, thanks to NuSTAR’s unprecedented spectral resolution above 10 keV, the lines can be disentangled. Together they remove a significant fraction of the continuum flux.

Our results therefore confirm and reconcile the previous measurements of the CRSF. In particular the measured high values of the CRSF energy by La Barbera et al. (2005) are in agreement with our high-energy line, while the line detected around 35–40 keV, in, for example, RXTE and Suzaku (Kreykenbohm 2004; Suchy et al. 2012), corresponds to our lower-energy line. Due to the lack of high-energy coverage (e.g., in the case of Suzaku/PIN) or lower-energy resolution (as in the BeppoSAX data), it has so far not been possible to disentangle both lines in a single data set.

By applying our model to the archival Suzaku data, we show that those data are fully consistent with the two CRSFs proposed here. Furthermore, we find no significant change in energy as a function of luminosity, despite covering about a factor of five in flux between NuSTAR observation 1 and Suzaku. Previous indications of a strong dependence of the line energy on luminosity are therefore very likely due to not being able to disentangle both lines.

5.1. The location and nature of the CRSF

Measuring two or more CRSFs in an accreting HMXB is not unusual (e.g., 4U 0115+63 with up to 5 lines; Heindl et al. 1999). In this case, the lines are interpreted as harmonics of the fundamental line, that is, they are produced through scattering with electrons in higher Landau levels. The line energies should therefore be harmonically related at multiples of the fundamental energy. Here, we find that the two lines are not related at the expected factor of two. Instead, we find a factor of ∼1.4. While relativistic effects can cause the ratio to deviate slightly from two (Pottschmidt et al. 2005; Müller et al. 2013), the difference is clearly much larger than expected (see, e.g., Mészáros 1992).

It is possible that the low-energy CRSF is not the fundamental line, but the first harmonic. In that case we would expect a fundamental line around 17 keV, which gives harmonics at ∼34 keV and ∼51 keV, as observed. However, the spectra contain no indication for a feature around that energy, for any continuum model. Adding another CRSF at this energy with a fixed width of 3 keV results in a very small equivalent width of (−70±110) eV which is consistent with zero.

The fundamental line can often be much shallower than the harmonics due to photon-spawning (Schwarm et al. 2017), in particular in a two-column geometry where it can almost vanish at all observation angles (Falkner et al. 2018). However, it seems unlikely that this effect would result in a perfectly filled line, without any emission wings or other deviations from a smooth spectrum.

Another possibility would be that we observe only one line, but with a significant deviation from the Gaussian shape. Such deviations are theoretically predicted, in particular in the fundamental line (Schönherr et al. 2007; Schwarm et al. 2017), and have been observed in other sources, for example, Cep X-4 (Fürst et al. 2015), V 0332+53 (Doroshenko et al. 2017, but see Pottschmidt et al. 2005 for further discussion), and 4U 1626−26 (Iwakiri et al., in prep.). However, in those cases the secondary absorption features were much weaker compared to the main line. In GX 301−2, the lower-energy feature is only about a factor four weaker than the high-energy one and more easily visible as a stand-alone line.

It is also possible that both CRSFs are formed in the same accretion column, but at different locations within the column. If we assume that the 50 keV line is produced at the surface of a standard neutron star with a radius of 10 km and the field behaves like a dipole, we can estimate the altitude at which the lower-energy line is formed. This assumption implies an altitude of about 1.4 km, taking into account the reduction of the gravitational redshift z as well. If both lines are formed within a single column, this approach would therefore require two line-forming regions: in a shock about 1.4 km above the surface, and at the surface of the neutron star itself (cf. Poutanen et al. 2013).

According to the calculations by Becker et al. (2012), GX 301−2 is in the gas shock regime, where they postulate that the CRSF is formed at the surface of the neutron star. However, this is only the case for a canonical neutron star of mass 1.4 M_⊙. For a slightly more massive neutron star of 1.8 M_⊙, similar to what is postulated for Vela X-1 (Quaintrell et al. 2003), GX 301−2 shifts into the Coulomb-radiation-dominated deceleration regime. In this regime, we can calculate the expected shock height using the formula

$$\begin{array}{c}\hfill h=1.48\times {10}^5{L}_{37}^{-5/7}{B}_{12}^{-4/7}\text{cm},\end{array}$$(3)

where L₃₇ is the luminosity in 10³⁷ erg s⁻¹ and where we assume 10 km for the radius of the neutron star (see Eq. (51) in Becker et al. 2012). Using the values we measure in observation 2 (L₃₇ = 0.2782 ± 0.0017, $B_{12}=5.{58}_{-0.14}^{+0.15}$, corrected for the gravitational redshift at the neutron star surface for a neutron star with a mass of 1.8 M_⊙) we find h = 1.48 ± 0.03 km (uncertainties are given only by propagating the statistical uncertainties of the measurements). This value is in very good agreement with the estimate from the CRSF energy.

We note however that in this regime we expect a positive correlation of the CRSF energy with luminosity, as the shock should be increasing in altitude with increased mass accretion rate. Such a correlation has been observed in Vela X-1 and Cep X-4 at a similar luminosity (Fürst et al. 2014; Vybornov et al. 2017). However, a correlation of the strength as observed in Cep X-4 is within the uncertainties of our values, and we therefore cannot rule it out nor confirm its presence.

5.2. Modeling of the phase dependence of the CRSF energy

Further information about the accretion and emission geometries can be obtained from the phase dependence of the CRSF energy. Here we concentrate on observation 2 due to the better constraints on all spectral parameters.

The 35 keV-line shows a very significant, almost sinusoidal variation as a function of phase (Fig. 7g). This variation could indicate that we sample different heights of the accretion column with different intrinsic magnetic field strengths over the rotational phase of the neutron star. However, to have such a clear variation, this means that at each phase we can only observe a very small region of the column, which could only be realized if the emission profile is very narrow and altitude dependent. It is unclear how such an emission profile would physically be produced.

Within the accretion column, the in-falling material can obtain relativistic velocities of about 0.3–0.5 c (Basko & Sunyaev 1976). This would result in the emitted radiation being strongly boosted towards the neutron star surface, and respectively experiencing a significant red-shift for the observer. As at different phases the angle of our line-of-sight towards the velocity vector changes, we observe different boosting factors and therefore different energies.

In this picture the amplitude of the phase-dependence of the CRSF energy strongly depends on the velocity in the line-forming region. If the 50 keV line is formed close to the surface, its velocity is basically zero, in agreement with the observed variation of energy against phase, which is compatible with a negligible phase-dependence of that line. On the other hand, the 35 keV CRSF is formed much higher in the column, experiencing a strong relativistic redshift. Here we present a simple model based on this idea. More details for the model setup will be described in a forthcoming publication (Falkner et al. 2018).

The model features a single accretion column with a negligible radius extended only in height. The observable energy,

$$\begin{array}{c}\hfill E=E{}^{\prime }\sqrt{1-\frac{R_{\mathrm{S}}}{R}}\frac{\sqrt{1-{\beta }^2}}{1+\beta cos{\eta }^{\star }},\end{array}$$(4)

is related to the intrinsic energy, E′, taking into account the gravitational redshift at the radius R from the central mass with the Schwarzschild radius R_S = 2GM_NS/c² and the boosting due to the local bulk velocity β = v/c, which depends on the emission angle η^⋆ with respect to the magnetic field (in the rest frame of the column).

In the simple picture of a cylindrical accretion column with negligible radius, η^⋆ can be identified with the emission angle in the rest frame of the accretion column α (Beloborodov 2002)

$$\begin{array}{cc}\hfill cos{\eta }^{\star }& =cos\alpha =1-(1-cos\Psi )(1-\frac{R_{\mathrm{S}}}{R}),\hfill \end{array}$$(5)

accounting for general relativistic light bending based on the Schwarzschild metric, where Ψ is the apparent emission angle.

The apparent emission angle is related to the geometrical setup by

$$\begin{array}{c}\hfill cos\Psi =cosicos{\mathrm{\Theta }}_{\mathrm{AC}}+sinisin{\mathrm{\Theta }}_{\mathrm{AC}}cos(\mathrm{\Phi }-{\mathrm{\Phi }}_{\mathrm{AC}}),\end{array}$$(6)

where i is the observer inclination, Φ the rotational phase, and Φ_AC and Θ_AC are the phase offset and the polar angle of the accretion column, respectively (Falkner et al. 2018).

The intrinsic energy $E_{\text{CRSF}}^{\prime }(h)$ in the rest-frame of the emitting plasma, which is emitted at a certain height h above the surface can be related to the observed surface energy E_CRSF(0), where we assume β = 0, by

$$E_{\text{CRSF}}^{\prime }(h)=E_{\text{CRSF}}(0)\frac{R_{\text{NS}}^3}{{(R_{\text{NS}}+h)}^3\sqrt{1-\frac{R_{\text{S}}}{R_{\text{NS}}}}}·$$(7)

Equation (7) accounts for the gravitational redshift and the decrease of the dipolar magnetic field. Further, E_CRSF(0) is a direct tracer of the surface magnetic field. Here we assume that the observed surface energy E_CRSF(0) is the phase-averaged observed energy of the high-energy CRSF line, 49.6 keV, and therefore that this line does not vary significantly as a function of phase.

Due to the many free parameters and degeneracies in this simple model, we perform MCMC calculations to estimate the parameter space using the emcee routine in ISIS. This routine is based on the parallel “simple stretch” method presented by Foreman-Mackey et al. (2013).

Figure 9 shows the probability distribution of the parameters and Table 3 gives the resulting parameter constraints. We used 50 free walkers for each free parameter and 100 000 iterations. The underlying model assumes that the two CRSFs originate from the same accretion column. Furthermore we fix Φ_AC to 0.562, which was found to be well constrained in a preliminary analysis and does not show any correlations to other parameters.

Fig. 9. Parameter probabilities for the model described in Sect. 5.2 obtained from Monte Carlo simulations (EMCEE). Parameter correlations are shown as color maps, where black corresponds to the highest probability. Purple, magenta, and orange correspond to the 99%, 90% and 68% probability levels. In the i − ΘAC space there are two solutions visible, which are divided by ΘAC = −i + 180° (green dashed line), where i is the inclination and ΘAC is the polar angle of the accretion column. The bottom panels show the one dimensional parameter probabilities, where red and blue histograms correspond to the individual solutions.

This simple model provides an excellent description of the observed phase-dependence of the CRSF energy, as shown in Figure 10. It also provides some limits on the geometry of the system; for example it indicates an inclination i > 20° and a velocity β between 0.2 and 0.4 c. These velocities are in good agreement with theoretical calculations. We also find a column height of around 1 km, which is very similar to the height we estimated independently from the shock height in the accretion column model by Becker et al. (2012).

Fig. 10. Phase dependence of the CRSF energies during NuSTAR observation 2, together with the predicted energies by our simple accretion column model in red. The energy of the higher-energy line has been fixed to 49.6 keV in the model.

We note that in a further step we could assume a certain velocity profile within the accretion column; for example, the shock model as discussed by Becker et al. (2012). In this case, the velocity below the shock would be much slower, and the emission region must be constrained to be around 1 km. This would lead to better constraints of the geometric parameters. However, a detailed calculation of this model is beyond the scope of the this paper.

Of course, a possible solution could also include contributions from both accretion columns. However, we do not consider this case here, as we cannot constrain the relative contribution of each column to the observed flux with the available data. In such a model, we would need to make sophisticated assumptions about the altitude-dependent emission profile in each column, which is not necessary in the present setup.

5.3. Summary

The very high data quality of the NuSTAR observations of GX 301−2 has allowed us to identify two CRSFs, one at 34 keV and one at 51 keV, where the secondary line feature is statistically significant at ≥99.9%. At the same time, we rule out with high confidence a fundamental line at 17 keV, and therefore the possibility that the two observed lines are harmonically related. We speculate on different origins for the two lines, but current available theoretical models do not yet allow to draw firm conclusions. We argue that a likely scenario is a symmetric magnetic field (and accretion) geometry, in which the CRSFs are formed at two different altitudes above the neutron star: the low-energy line at about 1–1.4 km above the surface in a shock, and the high-energy line at the surface of the star. Both lines are therefore interpreted as the fundamental line corresponding to the magnetic field strength at the respective line-forming altitude.

Based on this model, we calculated the expected energy variance of the lines solely due to the changing viewing angle with phase and the boosting in the relativistic velocities of the in-falling material. Here we do not make any assumptions about the configuration and emission profile of the accretion column. This simple model provides a very good description of the phase dependence of the CRSF energy and gives some constraints on the geometry of the neutron star and its magnetic field. Further observations of the phase dependence at different luminosities and theoretical calculations are needed to confirm this model.

Acknowledgments

We thank the referee for their helpful comments. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). We would like to thank John E. Davis for the slxfig module, which was used to produce all figures in this work. The Swift/BAT transient monitor results were provided by the Swift/BAT team. This research has made use of data obtained from the Suzaku satellite, a collaborative mission between the space agencies of Japan (JAXA) and the USA (NASA). L.N. acknowledges support by ASI/INAF grant I/037/12/0 and PRIN-INAF 2014 grant “Towards a unified picture of accretion in High Mass X-Ray Binaries”.

References

Appendix A: Fit parameters for all phase-averaged models

Tables A.3 and A.4 give the best-fit parameters for all tested continuum models, using one and two CRSFs, respectively. The values for the NPEX, FDcut, and HIGHECUT models given in Table A.4 are the same as in Table 2, but repeated for convenience. We note that in particular the cutoffpl model did not lead to an acceptable fit, and the given uncertainties on the fit parameters are therefore likely underestimated and should be regarded with care. In this model we also had to allow for a secondary, very-high-absortion column (N_H,2) as well as fix the energy of the Fe Kβ line to 7.05 keV.

All Tables

All Figures

Fig. 1. Average light curve from the Burst Alert Telescope (BAT) of Swift, folded on the orbital period where orbital phase 1 is periastron. A strong flare shortly before periastron is clearly visible. We note that due to the eccentric orbit, the phases are not equidistant in time. The two NuSTAR observations are shown by the red crosses and the Suzaku observations by the blue diamonds. These data use the right-hand y-scale with the flux measured between 5 and 50 keV. The left and right y-scales are aligned based on the daily BAT count rate measured during the NuSTAR observations. We note the strong variability above and below the average flux in the pointed observations of NuSTAR and Suzaku.

In the text

Fig. 2. Spectrum of NuSTAR observation 2. Panel a: data from NuSTAR/FPMA (FPMB not shown for clarity), together with the best-fit NPEX model (black). We also show the two components of the NPEX model separately, indicating the soft component with a dashed line and the harder component (for which the photon-index is fixed to −2) with a dotted line. Additionally, we show the model with each CRSF removed in turn (orange and green) and with both removed (magenta). Panel b: residuals in terms of χ for the best-fit model with only 1 CRSF. Panel c: residuals to the best-fit model with 2 CRSFs Panel d: residuals in terms of data-to-model ratio without the 50 keV CRSF. Panel e: residuals without the 35 keV CRSF. Panel f: residuals without both CRSFs. In this panel the non-symmetric shape of the residuals can be clearly seen.

In the text

Fig. 3. As in Fig. 2, but for observation 1. Panel a: additionally to the NuSTAR/FPMA data and best-fit NPEX model, we show the spectrum of observation 2 in gray for comparison. Panel b: residuals in terms of χ for the best-fit model with two CRSFs. Panel c: residuals in terms of data-to-model-ratio without the 50 keV CRSF. Panel d: residuals without the 35 keV CRSF. Panel e: residuals without both CRSFs.

In the text

	Fig. 4. Same as Fig. 2 but for the Suzaku spectrum. Panel a: XIS 1 data are shown in red, PIN data are shown in blue. Panel b: residuals in terms of χ for the best-fit model with 2 CRSFs. Panel c: residuals in terms of data-to-model-ratio without the 50 keV CRSF. Panel d: residuals without the 35 keV CRSF. Panel e: residuals without both CRSFs.
In the text

Fig. 5. Panel a: pulse profiles of observation 2 in three different energy bands: 3–6 keV (blue), 6–12 keV (yellow), and 12–78 keV (red). The profiles were normalized by subtracting the respective mean rate and dividing by the respective standard deviation to allow for easier comparison of the shape between energy bands. The lower panels show the hardness ratio between each pair of these pulse profiles. The hardness ratio is calculated as (H − S)/(H + S), where H is the profile in the harder band, and S the profile in the softer band. Panel b: hardness ratio between the 3–6 keV and 6–12 keV profiles. Panel c: hardness ratio between the 6–12 keV and 12–78 keV profiles. Panel d: hardness ratio between the 3–6 keV and 12–78 keV profiles. The dashed lines indicate the phases selected for phase-resolved spectroscopy.

In the text

Fig. 6. Color-coded energy resolved pulse profile map. The x-axis indicates phase (repeated once for clarity) and the y-axis energy in keV. The color indicates the relative strength in each energy band according to the scale on the right, normalized by (r − 〈r〉)/σ(r), where r is the rate in each bin, < r> is the average count-rate over each energy band, and σ(r) is the corresponding standard deviation.

In the text

Fig. 7. Panel a: flux between 3 and 40 keV in keV s−1 cm−2. Shown in black is the “soft” component, in red the “hard” component (see Sect. 3 for a description of the components). Panel b: covering fraction. Panel c: absorption column in 1022 cm−2. Panel d: photon index. Panel e: energy of CRSF1. Panel f: strength of CRSF1. Panel g: energy of CRSF2. Panel h: strength of CRSF2. Panel i: iron line flux in ph s−1 cm−2.

In the text

	Fig. 8. Same as Fig. 7, but for observation 1.
In the text

Fig. 9. Parameter probabilities for the model described in Sect. 5.2 obtained from Monte Carlo simulations (EMCEE). Parameter correlations are shown as color maps, where black corresponds to the highest probability. Purple, magenta, and orange correspond to the 99%, 90% and 68% probability levels. In the i − ΘAC space there are two solutions visible, which are divided by ΘAC = −i + 180° (green dashed line), where i is the inclination and ΘAC is the polar angle of the accretion column. The bottom panels show the one dimensional parameter probabilities, where red and blue histograms correspond to the individual solutions.

In the text

	Fig. 10. Phase dependence of the CRSF energies during NuSTAR observation 2, together with the predicted energies by our simple accretion column model in red. The energy of the higher-energy line has been fixed to 49.6 keV in the model.
In the text