© The Authors 2024

1. Introduction

Cyclotron resonant scattering features (CRSFs), also referred to as “cyclotron lines”, are characteristic absorption-line like features observed in the hard X-ray spectra of a number of accretion-powered X-ray binaries harboring a strongly magnetized neutron star (see Staubert et al. 2019, for a review). The X-ray spectra are expected to be formed in structures called accretion columns near the magnetic poles of the accreting neutron star (NS). Close to the NS, the accreted matter is funneled toward the magnetic poles due to the strong magnetic field, B ∼ 10¹² G. In the plasma of the accretion column, charged particles such as electrons experience Landau quantization, where their momenta perpendicular to the magnetic field lines are quantized. Inelastic scattering of photons off these electrons can result in the formation of an absorption-line like feature at an energy corresponding to the cyclotron energy ∼12 keV B/10¹² G.

The diagnostic power of such cyclotron lines lies in the possibility to infer the magnetic field strength at the region of their formation directly from the line centroid energy. In addition, cyclotron line parameters are sensitive to the other characteristics of the accreting plasma, such as temperature, density, and the velocity of the flow. Pulse-phase-resolved studies of CRSFs have further contributed to the investigation of the accretion geomerty of neutron stars and their emission characteristics (see, e.g., Maitra 2017). However, CRSFs are by no means observed in all accreting pulsars and the exact physical conditions necessary to form an observable CRSF are not yet well understood. Extensive observational campaigns over the last decades (e.g., Staubert et al. 2019; Pradhan et al. 2021) have found variability in CRSF parameters with a rotational phase of the pulsar as well as different types of correlations of the CRSF energy with X-ray luminosity. Different explanations for this behavior were proposed, involving, for example, the displacement of the line-forming region inside the accretion channel (e.g., Becker et al. 2012; Nishimura 2014), reflection from the NS surface (e.g., Poutanen et al. 2013), or a Doppler shift of photons during scattering in the accretion flow (Nishimura 2014; Mushtukov et al. 2015).

In sources that exhibit multiple cyclotron lines, the relation between those is of particular interest, raising the question about their origin as to whether they are individual features from different regions or harmonic lines. The spacing of the line energies allows for the region of their formation to be probed and it provides insights into the complex process of photon redistribution in resonant scattering events in an inhomogeneous medium. Kendziorra et al. (1992) were the first to claim CRSFs at 54 keV and possibly 27 keV for Vela X-1 (see Kretschmar et al. 2021, for an extensive review). Subsequent analyses confirmed the lower CRSF (see, e.g., Kretschmar et al. 1997; Kreykenbohm et al. 1999), whilst other studies (see, e.g., Orlandini et al. 1998) only found the higher CRSF. The phase-resolved analysis of an RXTE observation of Vela X-1 by Kreykenbohm et al. (2002) detected both lines and identified them as fundamental and second harmonic, which has been confirmed by many subsequent studies (see, e.g., Wang 2014; Fürst et al. 2014; Diez et al. 2022). Another prominent X-ray source with multiple cyclotron lines is 4U 0115+63. Its X-ray spectrum shows up to five cyclotron lines at energies of nearly integer multiples of 11 keV (Heindl et al. 1999; Ferrigno et al. 2009; Tsygankov et al. 2007), that are interpreted as a fundamental line and its higher harmonics (also, an additional fundamental line was reported at ∼16 keV, Iyer et al. 2015; Liu et al. 2020).

The focus of this work is the study of the presence and relation of the multiple cyclotron lines in GX 301−2 (alias 4U 1223−62). GX 301−2, discovered by McClintock et al. (1971), is a high-mass X-ray binary (HMXB) system that contains a NS rotating with a period of ∼700 s and accreting from the stellar wind of the B hypergiant Wray 977 (Bradt et al. 1977; Kaper et al. 2006) that is clumpy (Fürst et al. 2011; Roy et al. 2024). Under the assumption that the rotating vector model (see, e.g., Radhakrishnan & Cooke 1969) is applicable for GX 301−2 and that emission is dominated by the ordinary polarization mode, Suleimanov et al. (2023) determined the inclination and magnetic obliquity of the source using polarimetric measurements and report values of i_p = 135° ±17° and θ = 43° ±12°, respectively. Using the Gaia survey (Data Release 3, Gaia Collaboration 2016, 2023), Bailer-Jones et al. (2021) estimate a distance of GX 301−2 of $3.{55}_{-0.17}^{+0.19}$ kpc. GX 301−2 moves around its companion on an orbit with an eccentricity of 0.47 (Koh et al. 1997) and an orbital period of ∼41.5 d (see, e.g., Leahy 2002), showing rapid orbital decay (Manikantan et al. 2024). It shows regular pre-periastron flares, where the luminosity increases by roughly an order of magnitude compared to other orbital phases (Pravdo et al. 1995; Pravdo & Ghosh 2001; Islam & Paul 2014). Whilst the pre-periastron flare itself is regular, the timing and luminosity of the maximum cannot be predicted precisely, which complicates the scheduling of its observation. The X-ray pulsar is a well-established CRSF source with a line centered at ∼35 keV (Mihara 1995). The CRSF energy is strongly variable with the pulse phase between ∼30–40 keV (Kreykenbohm et al. 2004). Additionally, higher CRSF energies around 45–53 keV have been reported, leading to speculation about the luminosity dependence of the CRSF (La Barbera et al. 2005; Suchy et al. 2012) and whether changes in luminosity are responsible for the observation of higher CRSF energies.

The NUclear Spectroscopic Telescope ARray (NuSTAR, Harrison et al. 2013) observed the source twice between 2014–2015 at orbital phases of 0.67 and 0.87, where the luminosity was on the order of 3–4×10³⁶ erg s⁻¹, which is roughly one order of magnitude lower than during the periodic pre-periastron flares. Fürst et al. (2018) studied the two observations along with an archival Suzaku observation and report two nonharmonically spaced CRSFs, one around 35 keV and one at 50 keV. The existence of two separate lines was confirmed by Ding et al. (2021), using observations of GX 301−2 by the Hard X-ray Modulation Telescope Insight-HXMT (Zhang et al. 2020) at different orbital phases. Their systematic analysis includes data from about 40 pointed observations at different luminosities and finds a lower energy of the first CRSF of ∼26–35 keV and thus higher ratio of the centroid energies of the two CRSFs (≈1.63). We give an overview over past NuSTAR observations of GX 301−2 along with the determined cyclotron line energies in Table 1. Furthermore, we show the orbit of the NS around its companion along with the orbital phases of recent NuSTAR observations of GX 301−2 in Fig. 1.

Fig. 1. Orbital geometry of the system GX 301−2 based on Fig. 1 of Fürst et al. (2011). We show the NS orbit around its B hypergiant companion Wray 977 along with the orbital Swift (Krimm et al. 2013) count rate as determined from the long-term BAT light curve in blue. We further show the orbital phases of NuSTAR observations of GX 301−2 and the corresponding year in red. The 2022 observation is omitted as it was performed at nearly the same orbital phase as the 2014 observation. For the orbital parameters, we used the values presented by Kaper et al. (2006), that is, a sin i = 368.3°, e = 0.462, ω = 310.4°, and RWray 977 = 62 R⊙. We use an inclination of i = 55° (within the range given by Leahy & Kostka 2008). Epoch folding of the long-term Swift light curve yields an orbital period of 41.474 ± 0.006 d, which is in good agreement with the period stated by Doroshenko et al. (2010) as derived from Koh et al. (1997), including an orbital period derivative. We subsequently use the periastron passage TPA = 53531.65 MJD from Doroshenko et al. (2010). For the calculation of orbital phases, we take the orbital period derivative Ṗorb = −3.2 × 10−6 s s−1 (Doroshenko et al. 2010) into account.

Fürst et al. (2018) argue that unlike the first CRSF, the potential second harmonic does not show a clear phase variability. They propose a scenario where both CRSFs originate from different regions inside the same accretion column. In particular, one line is supposed to form at the bottom of the column and one at higher altitudes in a radiation-dominated shock, with the latter being subjected to relativistic boosting and therefore showing different pulse-phase variability. While the idea that in many accreting pulsars only one accretion column is visible for most of the rotational period was already discussed in the context of pulse profile modeling (Falkner 2018), in general, CRSFs with energies different from purely harmonic spacing are assumed to require the superposition of emission from two accretion columns and an asymmetric B-field or mass accretion (see, e.g., Liu et al. 2020).

For some HMXBs that show cyclotron lines, the centroid energy of the CRSF varies with the X-ray luminosity (see, e.g., Staubert et al. 2019, Fig. 9). Both positive and negative correlations have been observed for different sources. An analysis of a more recent NuSTAR observation of GX 301−2 at a luminosity of ∼0.3 × 10³⁷ erg s⁻¹ (in the observer’s frame, derived using d = 3.55 kpc in the energy range of 3–79 keV from Nabizadeh et al. 2019) resulted in higher centroid energies for both cyclotron lines compared to the observations studied by Fürst et al. (2018), where the source luminosity was on the order of 0.3–0.4 × 10³⁷ erg s⁻¹ (in the observer’s frame, derived using d = 3.55 kpc in the energy range of 5–50 keV). Nabizadeh et al. (2019) thus conclude a positive correlation between cyclotron line energy and luminosity of GX 301−2. The orbital phase-resolved analysis by Ding et al. (2021) showed a positive luminosity-dependence of both cyclotron line energies for luminosities below 10³⁷ erg s⁻¹ and a negative dependence above this threshold.

In this work, we investigate whether the newly discovered CRSF is also present during the pre-periastron flare, where the source luminosity is on the order of a few 10³⁷ erg s⁻¹ as well as analyze its pulse-phase dependence. We present a spectral and timing analysis of a NuSTAR observation of the pre-periastron flare and compare the source behavior to previous observations at different luminosity states. The remainder of the paper is structured as follows: We describe the data reduction in Sect. 2. In Sect. 3 we describe the spectral analysis of the whole observation as well as three time-segmented spectra. In Sect. 4 we present our spectral analysis of eight phase bins, where we focus on the phase dependence of the continuum and cyclotron line parameters. We discuss the results of the preceding sections from a physical standpoint and elaborate on their physical implications in Sect. 5. In Sect. 6 we briefly summarize the key findings of our analysis.

2. Data acquisition and reduction

In the following we describe the data acquisition, data reduction methods, and setup for the spectral analysis. GX 301−2 was observed during one of its regular pre-periastron flares (Pravdo et al. 1995) on 2020 December 27 at an orbital phase of 0.97 (as calculated from orbital parameters given in the caption of Fig. 1) with NuSTAR’s two nearly identical focal plane modules (FPMs). We show a long-time Swift lightcurve of GX 301−2 as well as the light curve of the NuSTAR observation along with its hardness ratio in Fig. 2. The NuSTAR observation captured the pre-periastron flare, but a significant fraction of the data was lost due to a telemetry issue. The resulting effective exposure time is ∼48.5 ks per FPM.

Fig. 2. NuSTAR observation of GX 301−2 in the context of its long-term behavior. a: hardness ratio, (H − S)/(H + S), of NuSTAR observation with hard energy range 10–30 keV and soft energy range 4–10 keV. The dashed lines indicate the division into time segments studied individually in Sect. 3.2. b: light curve of NuSTAR observation. c: Swift light curve between MJD 59100 and MJD 59300 (Krimm et al. 2013). We show ∼5 orbital periods and indicate the pre-periastron flares with orange bars.

We use nupipeline from HEAsoft v. 6.30.1 with CALDB version 20220802 for standard reprocessing and screening, and treat both detectors individually. We select circular source and background regions with 120″ radius each. The time-resolution of the extracted light curves is 1 s. Our analysis requires us to resolve spectral features above 50 keV as accurately as possible. We therefore apply the rebinning criteria proposed by Kaastra & Bleeker (2016), which are optimized for the energy resolution of the detector rather than ensuring a certain signal-to-noise. We use the χ²-statistic to evaluate spectral models.

Misalignments of few tens of eV between the energy grids of the two NuSTAR focal plane modules with respect to the energy of the Fe K edge became apparent in a preliminary analysis. By introducing an additive gain shift to the energy grids of both detectors, respectively, and fitting both offset parameters g_FPMA/FPMB, the misalignments can be corrected for. All spectral analysis was performed with ISIS, version 1.6.2–51 (Houck & Denicola 2000). Unless otherwise noted, uncertainties are given at a 90% level. Luminosities are given in the observer’s frame unless stated otherwise.

3. Phase-averaged spectroscopy

GX 301−2 is highly variable on all timescales from the pulse period to beyond the orbital period. In this section, we analyze the phase-averaged spectrum and its variation on the timescale of hours, that is, the duration of the observation. Figure 2 shows that the actual pre-periastron flare (shown with orange bars in the bottom panel) is already preceded by significant flaring activity on timescales of ks that are also associated with spectral variability (see, e.g., Leahy & Kostka 2008; Evangelista et al. 2010). The figure shows significant fluctuations in the hardness ratio throughout the orbital phase, aligning with the fluctuations in the absorption column density throughout the source’s orbit, as reported in literature (Islam & Paul 2014). The spectral analysis of GX 301−2 therefore requires a solid understanding of its timing properties. In Sect. 3.1 we analyze the spectrum of the complete observation. In Sect. 3.2 we divide the observation into three time segments based on changes in flux level and hardness ratio and analyze spectral changes on the timescale of hours.

3.1. Time-averaged spectral analysis

We conduct the spectral analysis by making use of purely empirical models. An extensive list of commonly used empirical continuum models is presented by Müller et al. (2013) and Staubert et al. (2019). We use a powerlaw component in combination with a Fermi-Dirac cutoff (FDcut, Tanaka et al. 1986)

$$\begin{array}{cc}& \mathtt{continuum}(E)=\mathtt{powerlaw}(E)\times \mathtt{FDcut}(E)\hfill \\ \hfill & =K·E^{-\mathrm{\Gamma }}\times {[1+exp\left(\frac{E-E_{\mathrm{cutoff}}}{E_{\mathrm{fold}}}\right)]}^{-1},\hfill \end{array}$$(1)

where K is a normalization constant, Γ the photon index, and E_cutoff/fold the cutoff and folding energy, respectively. The model was also used by Kreykenbohm et al. (2004) to successfully describe the X-ray continuum of GX 301−2 during the pre-periastron flare as observed with RXTE and more recently by Nabizadeh et al. (2019) in their analysis of a NuSTAR observation of GX 301−2 from 2019. Fürst et al. (2018) on the other hand used the NPEX model (Mihara 1995) to describe the continuum in NuSTAR observations from 2014 and 2015.

For the observation in focus of this publication, multiple empirical continuum models are tested in a preliminary analysis. Only a powerlaw including a FDcut cutoff provides a good description of the data and is therefore used throughout this spectral analysis. In general, we acknowledge the variety of continuum models that have been successfully used to describe the X-ray spectrum of GX 301−2 and emphasize that possible correlations between continuum parameters and cyclotron lines (see, e.g., Müller et al. 2013) call for caution when comparing cyclotron line parameters resulting from analyses with different continuum models.

In order to account for photoelectric absorption in the interstellar medium and by neutral material near the X-ray source, we use a partial covering model (PCF, see, e.g., Ballhausen 2021; Diez et al. 2022),

$$\begin{array}{c}\hfill \mathtt{PCF}(E)=f_{\mathrm{PCF}}\times \mathtt{tbab}{s}_1(E)+(1-f_{\mathrm{PCF}})\times \mathtt{tbab}{s}_2(E),\end{array}$$(2)

where f_PCF is the partial covering fraction and tbabs is the absorption model of Wilms et al. (2000). Our absorption model therefore has three parameters, f_PCF and the hydrogen column densities N_H1/2 of the two underlying absorption components, for which we assume the cross sections by Verner & Yakovlev (1995) and the elemental abundances by Wilms et al. (2000). We model both CRSFs as absorption lines with Gaussian optical depth using the model gabs (see, e.g., Staubert et al. 2019). Its parameters are the cyclotron line energy E_CRSF, its width σ_CRSF, and its strength d_CRSF, from which the optical depth at E_CRSF can be derived as $\tau =d_{\mathrm{CRSF}}/(\sqrt{2\pi }{\sigma }_{\mathrm{CRSF}})$. For an in-depth discussion of CRSFs and their modeling, see, for example, Staubert et al. (2019). We further include components for the two prominent neutral iron emission lines Fe Kα and Fe Kβ as additive Gaussian emission lines using the model egauss. It is defined as

$$\begin{array}{c}\hfill \mathtt{egauss}(E)=\frac{A}{\sqrt{2\pi }\sigma }exp[-\frac{{(E-E_0)}^2}{2{\sigma }^2}]\end{array}$$(3)

where A is the normalization representing the area under the Gaussian curve, E₀ is the center of the Gaussian curve, and σ its width. The soft X-ray spectrum of GX 301−2 was studied in detail by Watanabe et al. (2003) using Chandra High-Energy Transmission Grating (HETG) observations, that show that the widths of the iron emission lines are smaller than the energy resolution of NuSTAR (∼400 eV at 6.4 keV, see Harrison et al. 2013), for which reason we do not expect to be able to constrain them. We therefore fix the line widths to 1 eV. We furthermore fix both line energies to their theoretical values as given in the X-Ray Data Booklet¹ (based on Bearden 1967; Krause & Oliver 1979), that is E_Fe Kα = 6.397 keV and E_Fe Kβ = 7.058 keV. The centroid energy of the Fe Kα line was determined as the average of the two Fe energy levels contributing to the observed emission line at 6.391 keV and 6.404 keV. The study of the fluorescence lines is outside the scope of this paper, we therefore refer to Fürst et al. (2011) for a detailed study of the fluorescence lines in the X-ray spectrum of GX 301−2.

Our preliminary analysis shows the presence of an excess at the lower end of the NuSTAR energy range at 3.5–5.0 keV, which cannot be described correctly with the continuum model. We therefore include a blackbody component bbody to account for the soft excess. Using the above model provides no satisfactory description of the soft X-ray spectrum at energies of 5–7 keV. We therefore include an additional Gaussian line (hereafter denoted by index SE for soft emission), with centroid energy of ∼6 keV and width below 1 keV. An extensive justification for this component is given in Appendix B. We note that this added emission component is far narrower than the blackbody component and cannot be compensated for by varying blackbody temperature and normalization.

To account for detector flux calibration uncertainties a multiplicative factor C_FPMB for the FPMB spectrum is introduced. This parameter is expected to be well within ∼3% of unity (see, e.g., Madsen et al. 2022). To account for a shift in the energy grids of both detectors we include additive gain shift components (g_FPMA and g_FPMB), that shift the energy grids according to their value. The total model can therefore be expressed as

$$\begin{array}{cc}& f(E)=\mathtt{detconst}\times [(\mathtt{continuum}(E)+\mathtt{bbody}(E))\hfill \\ \hfill & \times \mathtt{PCF}(E)\times \mathtt{gab}{s}_1(E)\times \mathtt{gab}{s}_2(E)\hfill \\ \hfill & +\mathtt{egaus}{s}_{\alpha }(E)+\mathtt{egaus}{s}_{\beta }(E)+\mathtt{egaus}{s}_{\mathrm{SE}}(E)].\hfill \end{array}$$(4)

The two iron lines and the additional Gaussian line are not subject to absorption in our spectral model as defined in Eq. (4). The use of either unabsorbed or absorbed emission lines constitutes a simplification of the expected conditions, where the reprocessing of radiation within the ambient medium gives rise to the observed emission as well as absorption features simultaneously. A preliminary analysis shows that both fit statistic as well as line shape do not differ from a model in which the emission lines are also absorbed, we therefore use a model in which the lines are not absorbed.

In the remainder of this paper, the index 1 and term “lower line” denote the CRSF first claimed by Mihara (1995) around 35 keV; the index 2 and term “higher line” denote the CRSF discovered by Fürst et al. (2018) around 50 keV. Furthermore, unless otherwise stated, parameters by Fürst et al. (2018) used for comparison to our findings are results of their NPEX fit to the NuSTAR observation from 2015 (see their Table 2).

Our model with two cyclotron lines provides the most acceptable fit to the data among the tested models, with a goodness of fit statistics of χ²/d.o.f. = 720/477 = 1.51. We show the observed spectrum along with the best-fit model in the upper panel of Fig. 3. All free fit parameters are given in Table 2. We refer to this model as the whole observation fit (WOF) in the remainder of this paper. We further show χ residuals in the second panel as well as residuals of the fit with both cyclotron lines removed in succession in order to illustrate their impact on the spectral fit. The residuals of the best-fit model do not deviate significantly from the data and do not show a systematic structure, but rather stochastic variability. The only exception is an excess below 4 keV, that we attribute to the use of a simplified absorption model; this deviation does however not affect the analysis of the CRSFs.

Fig. 3. Continuum modeling of GX 301−2. a: spectra of the two NuSTAR FPMs in blue and purple with best fit model shown in black. Models with both CRSFs removed in turn and model without CRSFs are shown in green, red, and orange. b: χ residuals of best fit. c and d: residuals of best fit with both CRSFs removed in succession. e: best fit of alternative model with only one CRSF. f: best fit of alternative model without CRSFs. The spectra have been rebinned for plotting purposes.

The energy of the lower cyclotron line is $E_{\mathrm{CRSF}1}=37.4_{-1.0}^{+1.1}$ keV, which is roughly 3 keV higher than the value reported by Fürst et al. (2018) for observations at different orbital phase and luminosity. The width of the line is found to be ${\sigma }_{\mathrm{CRSF}1}=6.6_{-0.9}^{+1.3}$ keV, with a strength $d_{\mathrm{CRSF}1}=5.5_{-1.5}^{+2.4}$ keV, corresponding to an optical depth τ₁ = 0.33 at the centroid energy. The second cyclotron line has a line energy of $E_{\mathrm{CRSF}2}=51.5_{-1.0}^{+1.1}$ keV. This value is ∼2 keV higher than the values reported by Fürst et al. (2018). We further report a width of ${\sigma }_{\mathrm{CRSF}2}=5.1_{-0.8}^{+1.0}$ keV and a strength of $d_{\mathrm{CRSF}2}=7.7_{-2.2}^{+3.2}$ keV, corresponding to an optical depth τ₂ = 0.60 at E_CRSF2. It is evident that the lower CRSF is wider than the high-energy line. Since thermal broadening is proportional to the line energy, other broadening mechanisms are required to explain the observed behavior.

From an alternative fit model in which the Gaussian emission lines are subject to the partial covering absorption component, we determine the equivalent width of the absorbed Fe Kα emission line to be EW_Kα = 401 eV. This value is roughly a factor of 3 higher than that obtained by Fürst et al. (2018), who report a value of EW_Kα = 134 eV for their second NuSTAR observation, whose luminosity is lower by a factor of ∼8 compared to the observation studied here. During our spectral analysis we further tested models containing one and zero CRSFs. The residuals of these models are shown in Fig. 3e and f, and their CRSF parameters are given in Table 3. Both the residuals and fit statistic of the fit without CRSFs indicates that the model lacks a required component in the range 30–50 keV. The model with one CRSF contains a comparably narrow and shallow CRSF at 35 keV. The residuals indicate the presence of a second possibly also narrow CRSF at energies ∼50 keV. The inclusion of such a component however changes the parameters of the lower CRSF in an unexpected way: its centroid energy is shifted by 3 keV toward higher energies whilst nearly doubling in width and increasing sixfold in strength, which implies an increase in optical depth by factor 4. With a width of 5.1 keV, the second CRSF is also wider than expected from the residuals in Fig. 3, panel e. Investigation of the parameter space shows that no statistically favorable set of parameters that preserves the shallow and narrow nature of the lower cyclotron line exists. While the two features show very minor overlap at most, they nevertheless influence each other indirectly via changes in the overall continuum shape, that can lead to the observed unexpected change in CRSF parameters.

We further note that the fitted gain offsets are −32.2 eV and −61.2 eV for the two focal plane modules, that correspond to ∼0.75 and ∼1.5 NuSTAR energy bins, respectively. Grefenstette & Brightman (2022) quote expected gain offsets of ≈40 eV, which indicates that the found offsets are within the expected range. This is also consistent with Diez et al. (2023), who find an offset between NuSTAR and XMM-Newton of −87 eV and Ballhausen et al. (2020), who obtained gain offsets of −25 eV and −67 eV for the two NuSTAR FPMs, respectively. We therefore conclude that the obtained gain offsets are within the expected range.

We stress that while the residuals alone are a sufficiently strong argument for the inclusion of a second CRSF, we back up our finding by use of statistical means. The improvement in χ² statistic by including a second CRSF with three additional free parameters is found to be Δχ² = 100. We aim to show that the statistic improvement found above stems from the correct modeling of an actual observed absorption feature and is not merely a consequence of adding degrees of freedom to the fit model until it describes the spectrum sufficiently well. We therefore simulate spectra in order to show the statistical significance of the CRSF detection.

We simulate 100 000 observations of GX 301−2 based on the spectral model described in Eq. (4) with only CRSF1. The parameters are drawn from Gaussian distributions centered around the best-fit value as given in the left-most column of Table 2 and variance derived from the parameter uncertainty. We therefore make use of our best-fit model but remove CRSF2 manually. Alternatively, one could use the parameters of a model that contains only one CRSF ab initio. The statistical results of the two approaches may differ slightly as the parameters of CRSF1 change upon the inclusion of CRSF2. We argue however that the best description of the phase-averaged spectrum is obtained with CRSF1 as stated in the last column of Table 3, that is, in a fit with both CRSFs; we therefore base our statistical test on the best fit including both lines, where CRSF2 is removed manually. The spectra are simulated based on Poisson statistics.

We fit both models with one and two CRSFs respectively to the data. In the 100 000 simulated spectra and the corresponding fit models, we find no simulation showing an improvement of more than Δχ² = 100 in fit statistic upon the inclusion of CRSF2; we therefore detect the second CRSF with a statistical significance of > 99.999%.

Fürst et al. (2018) carried out a similar test and found one false positive out of 10 000 simulated and tested spectra, that corresponds to a false positive rate of ≤0.01%. It is noteworthy that the statistical test of the 2020 NuSTAR observation shows a much higher detection confidence, which is most likely due to the fact that the observation at hand has a higher luminosity of factor ∼8.

We lastly comment on the present statistical test, which is used both in this work as well as by Fürst et al. (2018). We assume that parameters are Gaussian-distributed, which is not necessarily the case. Moreover, drawing parameters independently from each other implies the assumption that they are uncorrelated, which is generally not justified for such a phenomenological model. As a result of this, we most likely sample a much larger parameter space than intended, which also contains spectra which do not agree with our best-fit model within uncertainty. Therefore, our statistical test should be understood as a lower limit of the significance of the presence of CRSF2 in the phase-averaged spectrum.

Apart from the solution presented here with two CRSFs, another solution with comparable fit statistic exists, in which σ_CRSF1 exceeds 10 keV. Cyclotron lines are expected to be broadened by a variety of causes: thermal motion of electrons and inhomogenities in the conditions of the line-forming region such as magnetic field strength, plasma temperature, and density. Lastly, the viewing angle on the accretion column is also expected to impact the observed line shape in the phase-averaged spectrum (see, e.g., Schwarm 2017), which introduces a pulse-phase variability as the viewing angle changes over the pulse period. All of the above can lead to wider CRSF profiles and potentially complex line shapes. Here, however, we argue that the “narrow-CRSF” solution for the phase-averaged spectrum is more physically motivated for the following reasons:

1. The wide CRSF removes ∼5% of the total counts, whereas the narrow CRSF in our best-fit model removes ∼1%. Although physically plausible cases of exceptionally influential CRSFs have been found before (Nakajima et al. 2010; Malacaria et al. 2023), Staubert et al. (2019, Table A.5) show that CRSFs with σ > 10 keV are more the exception than the rule. Furthermore, the table shows that such wide CRSFs are predominantly found for high-energy CRSFs with line energies E_CRSF ≳ 50 keV.

2. The residuals shown in Fig. 3f do not support the claim of a cyclotron line width > 10 keV as the absorption feature at ∼35 keV is at most ∼8 keV wide.

3. The existence of the solution with a broad CRSF parameters can be explained by correlations between continuum and CRSF parameters (see, e.g., Müller et al. 2013), where the CRSF models a part of the continuum and no longer describes the absorption feature caused by cyclotron resonant scattering.

We therefore continue with the narrow-CRSF solution described and shown above as the basis of physical interpretation of the phase-averaged spectrum.

3.2. Time-segmented spectral analysis

The upper panel of Fig. 2 shows a significant temporal change in the hardness ratio of the 4–10 keV and 10–30 keV bands. This variation could be due to a change of the continuum shape, that in turn could influence the modeling of the CRSFs. In order to better understand the temporal variation of the spectrum that precedes the pre-periastron flare and the change in absorption during the flare, we divide the observation in three time segments marked as X, Y, and Z (TimeSegment X, TS Y, and TS Z). The segmentation is carried out in consideration of changes in both the hardness ratio and flux level of the source as indicated by the vertical line in the upper and middle panel of Fig. 2. Alternatively to such a temporal division, one can argue that a division only with regard to hardness ratio is more appropriate, but such an approach lacks the ability to gain an understanding of the interesting temporal evolution of the source at the onset of the pre-periastron flare. Time segment Z coincides with the flaring event, where the source reaches its peak luminosity (see, e.g., Kreykenbohm et al. 2004, their Fig. 1, or the lower panel of Fig. 2).

In order to identify potential changes in the spectral components between the three time segments, we start by inspecting the ratio between the time segment spectra and that of the whole observation, which are shown in Fig. 4. Deviations from a horizontal line indicate possible spectral changes during the time segment with respect to the full observation. Most notably, the three ratios do not match in shape or trend over the whole spectral range. Above 10 keV, especially TS X and TS Y show significant curvatures of concave and convex nature with respect to the expected horizontal line, respectively. The most profound differences however are present in the energy band up to 10 keV, most notably around the prominent iron line at 6.4 keV. This behavior is a potential indicator of a change in absorption during the observation.

Fig. 4. Ratios between spectra of the three time segments versus the spectrum of the whole observation. The dashed line indicates the behavior expected in the case of no spectral change.

To assess the spectral changes among the time-segmented spectra that are indicated by the ratios shown in Fig. 4, we apply the continuum model (Eq. (4)) to the individual spectra extracted for the three time segments. The χ residuals of the fits are shown in Fig. 5, the fit parameters are given in Table 2.

Fig. 5. χ residuals of best fits to time-segmented spectra.

Apart from few outliers at ∼6 keV the residuals of the fits to the time-segmented spectra show no systematic deviation from the data. We find that the absorption model parameters change during TS Y as N_H, 1 increases by a factor of 1.7 whilst N_H, 2 decreases by 25%. The partial covering fraction decreases over the whole observation from 0.64 to 0.46. The continuum parameters such as photon index, cutoff and folding energy also differ between the time segments, that is expected due to the observed curvatures in Fig. 4. The blackbody temperature increases gradually over the observation. For TS X and TS Y, the upper limit of its contribution to the spectrum cannot be constrained; for TS Y we find a relative norm that exceeds the value obtained for the whole observation by a factor of ∼60, which most likely massively overestimates the actual contribution of the blackbody component, as we expect its relative contribution with respect to the continuum to be well below a factor of 10. The contribution of both iron lines increases over the duration of the observation, the normalization of the Gaussian soft component increases from TS X to TS Y and remains constant within uncertainty from TS Y to TS Z.

We now turn to the parameters in the focus of this study, those of the CRSFs. The line energy of the lower line remains constant within uncertainty intervals and is consistent with the value determined for the whole observation in all time segments. The line widths and strengths during TS X and TS Z are in agreement with those width obtained for the whole observation, for TS Y the width increases to 10 keV and its strength increases to 14 keV. For CRSF2, we also find that the line energies are consistent with the ones obtained in Sect. 3.1. Both the width and strengths found in the time-segmented spectra are also consistent with the WOF. We note that due to the reduced number of counts in the segmented spectra the CRSF parameters are generally less well-constrained as in the fit of the whole observation.

In conclusion, we find that both the continuum shape and absorption are variable on the timescales of the observation. The cyclotron line parameters however vary only negligibly and are mostly consistent with the WOF values as shown by the above time-segmented analysis. The above claim is also supported by a preliminary analysis of the whole observation in the hard X-ray band from 10–79 keV, which poses the advantage of being mostly unaffected by changes in absorption; it further confirms that the excess below 4 keV in Fig. 3b does not affect the analysis of the cyclotron lines. We therefore conclude that it is justified to conduct the subsequent phase-resolved analysis for the whole observation together instead of a time-segmented manner as in Sect. 3.2.

4. Phase-resolved spectroscopy

After analyzing the phase-averaged spectrum of GX 301−2, we now turn to the analysis of the phase-resolved spectrum. Using epoch folding (see, e.g., Leahy et al. 1983), we determine the pulse period of the NS to be P = 671.62 ± 0.05 s based on NuSTAR light curves with 1 s binning. We show both the pulse profile for different energy bands as well as the hardness ratio in Fig. 6. The pulse profile shows two peaks. We define the phase zero as the minimum between the main and weaker pulse. The main peak is at a phase of φ ∼ 0.6–0.9 and a weaker peak at φ ∼ 0.1–0.3. A similar pulse profile structure was also obtained by Alonso-Hernández et al. (2022, their Fig. 3), who analyzed the 0.5–78 keV pulse profile of GX 301−2.

Fig. 6. Behavior of GX 301−2 over the pulse phase. Top panel: pulse profiles for different energy bands. Bottom panel: hardness ratio (H − S)/(H + S) with soft energy range 4–10 keV and hard energy range 10–30 keV.

To further illustrate GX 301−2’s behavior over pulse phase, we show a phase energy map, where count rates are shown as a function of both pulse phase and energy, in Fig. 7. The peak at φ ∼ 0.2 is most prominent at energies exceeding 10 keV and less pronounced at energies below ∼6 keV. The main peak at φ ∼ 0.7, on the other hand, is also present in the soft X-ray band. Both peaks show a splitting into two components toward lower energies. At energies above 35 keV a phase shift of the peaks is present. To further investigate the pulsations and especially whether the Fe Kα shows weaker pulsations than other soft components, we show the pulsed fraction as a function of energy in Fig. 8.

Fig. 7. Phase energy map for energies between 3.5–60 keV. For each energy bin, we show (N − ⟨N⟩)/(Nmax − Nmin).

Fig. 8. Pulsed fraction as a function of energy. We define the pulsed fraction as (Nmax − Nmin)/(Nmax + Nmin), where Nmax/min are the maximum and minimum rates of the pulse profile in a given energy bin.

It is evident that the pulsed fraction is highly variable with energy. In the 3.5–20 keV range it shows a slight decrease from 0.4 to 0.3. Toward higher energies the pulsed fraction increases significantly up to 0.9 at 50 keV. Most notably, the Fe Kα line shows far weaker pulsations than other soft components as the pulsed fraction decreases from 0.4 to ∼0.1 at the energy of the iron emission line. Such a behavior was previously reported by Fürst et al. (2011, Fig. 7 therein) based on a XMM-Newton observation of GX 301−2 during the pre-periastron flare.

The main focus of our phase-resolved analysis is to assess the presence of the two cyclotron lines and if present determine their pulse-phase variability. For this purpose, we fit our model as in Eq. (4) to eight phase bins that are extracted from the data of the whole observation.

We choose the phase bins by considering both the pulse profile and hardness ratios, where we aim to divide bins wherever significant changes in either quantity are present, as marked in Fig. 6 with vertical dashed lines. The phase-resolved background spectra contain very few counts due to their respective short exposure times. As the background is not expected to vary with pulse phase, we use the phase-averaged background spectra also for the phase-resolved analysis to improve statistics. We freeze parameters that are not expected to vary with pulse phase and show little to no phase dependence in a preliminary analysis. They include the detector constant, gain shift, blackbody temperature, Gaussian emission lines, and the Gaussian soft emission component. We show the phase dependencies of all free parameters in Fig. 9 and the residuals of the fits in Fig. 10.

Fig. 9. Phase dependencies of all free parameters in the phase-averaged analysis. Red represents the phase-averaged parameter values. Gray represents the normalized pulse profile. Arrows indicate where uncertainty intervals hit their predefined range limits.

Fig. 10. χ residuals of fits to the phase-resolved spectra. We use the same color-coding of NuSTAR focal plane modules as in above figures.

We find that CRSF1 is present in all eight phase bins, as residuals of fits excluding it show deviations from the model around 35 keV. For the higher CRSF, the picture is not as clear. At its high energy, background events and the low count rate complicate the investigation of its presence. We consider as CRSF2 as detected only if the following three criteria are fulfilled:

1. The residuals indicate a missing model component between 45–60 keV (see Table 4, second column).

   Table 4.

   Summary of the search for CRSF2 in the phase-resolved spectra.

2. The fitted energy, width, and strength of CRSF2 be constrained within reason² or a fit with E_CRSF2 and σ_CRSF2 as in the phase-averaged analysis in Sect. 3.1 and free d_CRSF2 provides a good description of the data (see Table 4, third column).

3. A statistical test supports the detection claim of CRSF2 with at least 3σ confidence (see Table 4, fourth through seventh column).

We note that in other publications the absence of a CRSF is often statistically substantiated by showing that the strength of the fitted CRSF is consistent with zero (see, e.g., König et al. 2020). This approach is however not feasible in this case due to the low signal-to-noise ratio at energies exceeding ∼50 keV: even if no CRSF is present, the included gabs component models either part of the continuum and cutoff or narrow continuum troughs of stochastic nature. In either case, the fitted strength is not consistent with zero. However, in those cases the improvement in χ² is sufficiently small to utilize the same statistical test as used for the phase-averaged analysis in Sect. 3.1.

We employ the above described strategy for each phase bin and summarize our findings in Table 4. Based on our established criteria, we detect CRSF2 with at least 3σ significance in phase bins 1, 2, and 4. We stress that the non-detection in the other phase bins can either be due to the abscence of the spectral feature during those time or insufficient count rates to detect the spectral feature. We therefore only claim the detection of CRSF2 in phase bins 1,2, and 4, but not that CRSF2 is absent in the other phase bins.

Both absorption column densities show little variation with pulse phase and are evenly distributed around the value obtained in the phase-averaged analysis. The partial covering fraction traces the pulse profile. The blackbody component is also variable over the pulse phase and decreases with pulse phase apart from two outliers in phase bins 1 and 3. Apart from one outlier in phase bin 4, the folding energy is consistently below the phase-averaged value and is variable in the range 5.5–6.5 keV. To investigate why the phase-resolved values are not evenly distributed around the phase-averaged value we remind that the continuum and CRSF1 parameters in Sect. 3.1 changed significantly upon the inclusion of CRSF2. In a fit of the phase-averaged spectrum with only CRSF1 present, we find a folding energy of $5.{812}_{-0.022}^{+0.049}$ keV; we find that the phase-resolved values for E_fold are distributed evenly around it. The cutoff energy is consistently above the phase-averaged value and varies with pulse phase, for which a similar argument as for the folding energy can be made. Its lower confidence limit cannot be constrained in phase bin 4, that is most likely the result of an interplay with CRSF2, whose strength is significantly overestimated in this phase bin and also cannot be fully constrained. The powerlaw flux varies with pulse phase as expected; the photon index Γ also shows significant pulse-phase variability. For CRSF1 we find that the line energy is generally lower than that found in the phase-averaged analysis and increases from ∼32 keV to ∼41 keV from phases 0.4–1.0. The width of CRSF1 is lower than the value obtained in the phase-averaged analysis and is consistent with a value of ∼4 keV. The strength is also notably weaker; the line appears to be shallow, but continuously detectable. For CRSF2, the line energy increases over the three phase bins during which it was detected. The widths scatter around the value determined in the phase-averaged analysis. The strength in phase bins 1 and 2 is consistent with the phase-averaged value, in phase bin 4 however its upper limit cannot be constrained and the parameter value itself most likely significantly overestimates its actual contribution due to correlations with the continuum parameters, in particular the cutoff and folding energies. In summary, whilst CRSF1 is detectable over the whole pulse period, CRSF2 can only be detected with high statistical significance during a fraction of the pulse period.

5. Discussion

In the discussion of our findings we focus on the properties of the detected cyclotron lines. In Sect. 5.1, we discuss the phase-variability of the cyclotron lines. In Sect. 5.2 we discuss the relation between the two cyclotron lines and their potential independence from each other. In Sect. 5.3 we interpret our findings from a physical point of view and propose scenarios for the formation of both lines. In Sect. 5.4 we discuss the luminosity-dependence of both CRSFs. In Sect. 5.5 we discuss an alternative interpretation of the hard X-ray spectrum of GX 301−2 by considering emission wings of cyclotron lines.

5.1. Phase-variability of the CRSFs

In this section we aim to compare the phase dependencies of the cyclotron line energies we observed with those obtained by Fürst et al. (2018). We show the cyclotron line energies as a function of pulse phase for both this work as well as Fürst et al. (2018) in Fig. 11.

Fig. 11. Phase dependencies of both cyclotron lines as determined in this work (blue circles) and by Fürst et al. (2018, black squares). We also show our phase-averaged values in red.

For the lower line we find that our obtained values are in very good agreement with those obtained for the previous observation. For the higher line, we find that our obtained values for phase bins 1 and 2 are in good agreement whereas the value for phase bin 4 is a few keV higher in comparison.

In order to further assess the pulse-phase dependency of the newly discovered line it is desirable to determine its dependency at different orbital phases. Ding et al. (2021) conduct an orbital phase-resolved analysis using Insight-HXMT and detected both cyclotron lines at different orbital phases. Furthermore, the low count rates at energies of 50 keV and higher significantly impede both the analysis of the higher cyclotron line and its fitting using commonly used algorithms. Observations with longer exposure time and more sensitive X-ray observatories can solve this problem to some extent.

5.2. Harmonicity of the CRSFs

We discuss the possibility of a harmonic relation between the two CRSFs, in order to discern their relative dependence. Other HMXBs such as 4U 0115+63 (see, e.g., Heindl et al. 1999; Ferrigno et al. 2009) and V0332+53 (Tsygankov et al. 2006) are known to exhibit multiple harmonics of a CRSF in their X-ray spectra, with the harmonic lines exhibiting nearly integer multiples of the fundamental line energy. Some anharmonicity of the lines is expected from relativistic theory, as given by the resonance condition (see, e.g., Eq. (1) of Araya & Harding 1999)

$$\begin{array}{c}\hfill ħ{\omega }_n=\frac{\sqrt{1+2n(B/B_{\mathrm{c}}){sin}^2\theta }-1}{{sin}^2\theta }{m}_{\mathrm{e}}{c}^2,\end{array}$$(5)

where n is the number of a cyclotron harmonic (n = 1 corresponds to the fundamental, which is considered the first harmonic), B_c ≈ 4.413 × 10¹³ G is the critical field given by the Schwinger limit, θ is the viewing angle with respect to the magnetic field axis, and m_ec² is the electron rest energy. For the smaller magnetic fields discussed here, with B/B_c ∼ 0.1, the ratio between the second harmonic and the fundamental resonance energy is expected to be between 1.9–2.0, regardless of observer inclination.

Here, we find a ratio between the line energies of ∼1.4, making a harmonic relation with CRSF1 as the fundamental line unlikely. We nevertheless discuss a scenario where both lines are harmonically related, as follows.

The two observed lines could be higher harmonics of a fundamental line at ∼17 keV, whose second and third harmonic lines would be at 34 keV and 51 keV, respectively. In this case, a ratio of ∼1.4 between ℏω₃ and ℏω₂ is expected (see Eq. (5)). However, in the many analyses of the X-ray spectrum of GX 301−2, such a CRSF was never found. This is corroborated by our analysis, as the residuals displayed in Fig. 3b do not indicate that the description of the data between 15–20 keV lacks a required model component. The non-detection of a CRSF at 17 keV could however be attributed to photon spawning, that would lead to a less pronounced fundamental line. Excitation of electrons to higher Landau levels, n ≥ 2, most likely leads to radiative decay to the closest lower level, n − 1. The energy of the emitted “spawned” photon is equal to the energy difference between the Landau levels. This leads to an excess of photons near the fundamental line energy (see, e.g., Schwarm et al. 2017a,b), that could potentially fill up the otherwise observable apparent absorption feature or even lead to an excess at this energy. Given both NuSTAR’s energy resolution of 400 eV below ∼40 keV (Harrison et al. 2013, their Fig. 8) as well as the high luminosity of the observation, we expect that even if photon spawning plays a significant role in the filling of the fundamental line, there should be a detectable shallow absorption or weak emission feature at the fundamental line energy, as it is highly unlikely that the process of photon spawning exactly compensates the absorption feature³. No such features were detected in this observation, or have been reported previously. In addition, this scenario is inconsistent with sources such as 4U 0115+63 and V0332+52, where photon spawning provides such an insignificant contribution that the clear classical ratio of harmonics is preserved to a high degree, together with a prominent fundamental line in the spectra (see, e.g., Ferrigno et al. 2009). This is made even more puzzling by the lower magnetic field of these sources which is expected to lead to higher efficiency of this process.

It is noteworthy that multiple anharmonic CRSFs have already been claimed at least once before the year 2018 for GX 301−2: Orlandini et al. (2000) found that the BeppoSAX spectrum of GX 301−2 requires two presumably independent CRSFs. They found time-variable line energies around 20–25 keV and 45–55 keV (Orlandini et al. 2000, their Fig. 1), where the temporal variability occurs on the order of tens of days, that is most likely related to the orbital cycle of the NS and the consequent change in mass accretion rate. Their higher line can be identified as the 50 keV line discovered by Fürst et al. (2018). However, Orlandini & Fiume (2001) later list the parameters of the higher CRSF as E_CRSF = 49.5 ± 1.0 keV and σ_CRSF = 17.9 ± 2.5 keV, where at least the latter parameter raises concerns on whether an actual cyclotron resonance scattering feature was detected. Such wide CRSFs have not been observed in GX 301−2 with NuSTAR observations and it is rather unlikely that the Doppler-broadened absorption feature would be so wide. To further illustrate this, we refer to Meszaros & Nagel (1985, their Eq. (1)), who give an angle-dependent formula for the width of CRSFs, which can be expressed as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{CRSF}}=2.55\mathrm{keV}·\frac{E_{\mathrm{CRSF}}}{30\mathrm{keV}}·\sqrt{\frac{k{T}_{\mathrm{e}}}{1\mathrm{keV}}}·cos\theta \end{array}$$(6)

where k is the Boltzmann constant, T_e the plasma temperature, and θ the observer’s inclination to the magnetic field axis. One can see that the width stated by Orlandini & Fiume (2001) would require a plasma temperature of at least ∼18 keV, that significantly exceeds the theoretically expected plasma temperature of ∼5 keV (see, e.g., Becker & Wolff 2007) as well as experimental findings such as ∼4–10 keV for A 0535+262 (see Ballhausen et al. 2017). Moreover, for most sources exhibiting electron cyclotron lines the width of the features are well below ∼12 keV (see, e.g., Staubert et al. 2019, Table A.5). The preliminary analysis of a more recent NuSTAR observation of GX 301−2 (McCulloch et al., in prep.) shows that both CRSFs are present in the spectrum at energies consistent with values stated by Fürst et al. (2018).

We further note that consequent anharmonically spaced cyclotron lines have been observed in the spectra of other accreting X-ray pulsars and were nonetheless interpreted as harmonics. For example, a second cyclotron line was observed in the spectra of Cep X-4, where a known fundamental line at ∼30 keV (Mihara et al. 1991) was supplemented by a harmonic at ∼45 keV (from a Suzaku observation of the 2014 outburst, see Jaisawal & Naik 2015) and ∼55 keV (from a NuSTAR observation of the same outburst, see Vybornov et al. 2017). Similarly, in Cen X-3, a second harmonic at ∼47 keV was recently claimed by Yang et al. (2023) from an Insight-HXMT observation, in addition to a previously discovered line at ∼28 keV (Audley 1997; Santangelo et al. 1998). A fundamental and a second harmonic were found at ∼44 keV and ∼73 keV in the spectrum of MAXI J1409−619 (Orlandini et al. 2012).

Jaisawal & Naik (2015) invoked a possible theoretical explanation for the nonharmonic ratio of cyclotron lines observed in Cep X-4. Nishimura (2005) and Schönherr et al. (2007) showed that if the magnetic field strength varies within the line-forming region, the ratio of subsequent harmonics deviates from 1.9, and instead, the ratio increases for magnetic fields that decrease with the height of the line-forming region and vice versa (see also Kumar et al. 2022). According to Schönherr et al. (2007), spawned photons play a major role in this process, providing more efficient filling of the fundamental line.

Cyclotron line energies have been found to vary both with time and luminosity. For Vela X-1, Fürst et al. (2014) find a positive correlation of the second harmonic line whilst the fundamental line shows a negative correlation for luminosities below ∼7 × 10³⁶ erg s⁻¹, followed by a flattening up to ∼10³⁷ erg s⁻¹. This implies that there is a direct proportionality between line energy ratios and luminosity, and therefore that the ratios of cyclotron line energies are also not necessarily constant over time. More recently, Diez et al. (2022) also report variability in the cyclotron line energies of both the first and second harmonic, but could not confirm an overall positive correlation of the cyclotron energies with luminosity. On the other hand, Liu et al. (2022) find a positive luminosity correlation for the fundamental line using Insight-HXMT observations, whilst the second harmonic appears to be uncorrelated with luminosity. Based on 11 observations taken over 1.5 years, they report a decrease in the line energy ratio from ∼2.0 to ∼1.7.

Despite the discrepancies, the aforementioned studies all agree on varying line energy ratios over time, furthermore providing a range of ω₂/ω₁ ∼ 1.5–1.7. Although this is slightly higher than our obtained line energy ratio of ∼1.4, they are significantly lower than the ∼1.9 expected from simple theoretical considerations. We display the ratio between the cyclotron line energies E_CRSF2/E_CRSF1 of GX 301−2 as stated in recent publications in Table 5. It is evident that our line energy ratio is the lowest to date, while the other reported values lie within the range observed in other sources. Our obtained value is in agreement with the two line energy ratios reported by Fürst et al. (2018) within uncertainties. As indicated by the examples from literature above, a non-integer line energy ratio alone is no convincing argument that two cyclotron lines are not harmonically related. Below, we therefore discuss possible physical scenarios for both the harmonic and anharmonic case.

5.3. Line-forming regions and physical interpretations

In consideration of the results of the phase-averaged analysis in Sect. 3, most notably the non-integer line energy ratio, and the phase-resolved analysis in Sect. 4, which showed that CRSF2 is only detectable during a fraction of the pulse phase, we here discuss possible formation scenarios for the two CRSFs in both the harmonic and anharmonic case. We further debate multiple physical scenarios which could result in this behavior.

Since CRSFs serve as probes for the magnetic field strength of the NS, the presence of two cyclotron lines in GX 301−2’s spectrum implies that if independently formed, they must form in regions with different magnetic field strengths. From Staubert et al. (2019), Eq. (1) therein, we can estimate the field strengths to be B₁ ≈ 4.2 × 10¹² G, and B₂ ≈ 5.8 × 10¹² G, in the corresponding line-forming regions. In the following, we briefly discuss possible explanations for the observed CRSFs and where they are formed. First, we revisit Fürst et al. (2018)’s interpretation that the high-energy line is formed at the NS surface whereas the lower line is formed at a height of ∼1.5 km. Fürst et al. (2018) argue that no phase-variability is expected for CRSF2 in this case, which is corroborated by their analysis. However, we would like to note that due to the strong dependence of the cyclotron resonance profile with respect to the angle between observer inclination and the magnetic field axis, features formed at the surface of the NS are still expected to show significant phase-variability (see, e.g., Schwarm et al. 2017a, Fig. 5). Our analysis instead shows CRSF1 present during all pulse phases and CRSF2 only present in phase bins 1, 2, and 4. However, this discrepancy can be resolved, given the following geometry: the orientation of the accretion column could be such that its bottom is obstructed for phases 0 and 5–8 (neglecting phase bin 3), resulting in the non-detection of CRSF2 there. CRSF1 however is expected to form at an altitude of 1.5 km above the surface and can either be observed directly without obstruction by the NS or aided by relativistic light bending around the NS (see, e.g., Falkner 2018). In conclusion, there indeed exists a scenario where two independently formed cyclotron lines could have formed within the same accretion column, and be in agreement with the observed pulse-phase variability.

Further, we consider the possibility that the two lines form within different accretion columns situated at the magnetic poles of the NS. In this picture, the observed phenomenon can be understood based on arguments related to the system geometry or magnetic field structure. For pulsars, we expect that the position of the accretion columns does not coincide with the rotational axis of the NS and are separated by an offset angle χ. For such inclined rotators, the magnetospheric radius decreases with increasing offset angle (Bozzo et al. 2018). This could result in different mass-accretion rates onto the poles of the NS. In this case, the CRSFs may form at different heights within the two columns. In such a scenario, these two independently formed CRSFs can be observed as anharmonic lines with differing pulse-phase variability.

Similarly, even with symmetric mass accretion rates, an explanation for the observed cyclotron line behavior can be found. Magnetic field structures can be more complex than the often assumed ideal dipole field. For example, distorted dipole fields that are invoked to explain asymmetric pulse profiles as shown in Fig. 6 (Kraus et al. 1995), could lead to the observed CRSFs. If present, such an asymmetric magnetic field would naturally lead to differing magnetic field strengths in line-forming regions within the accretion columns, which subsequently leads to the formation of CRSFs with differing centroid energies. Jaisawal & Naik (2015) suggest that a significant phase shift in CRSF parameters is expected in this case.

A distorted dipole field has already been invoked for other systems. For example, Liu et al. (2020) detected two fundamental CRSFs at 12 keV and 16 keV in the X-ray spectrum of 4U 0115+63, and based on the phase-variability of the line equivalent widths, argued that the lines each form in different accretion columns, in regions with magnetic field strengths of ∼1.1 × 10¹² G and ∼1.4 × 10¹² G. The small difference in line energy makes it difficult to disentangle the higher harmonics of each fundamental line Liu et al. (2020, Sect. 4.1).

The results of our analyses of GX 301−2 presented here can be explained by any of the aforementioned scenarios and further research is required to definitively determine whether the two cyclotron lines form within the same accretion column and if so, whether they are harmonically related.

5.4. Cyclotron line energies in historic context

Since the first claim of a cyclotron line in the spectrum of GX 301−2 by Mihara (1995) nearly 30 years ago, many subsequent analyses reported cyclotron lines at energies in the range 30–55 keV. In this section we collate many CRSF energy values for GX 301−2 from literature and investigate possible explanations for the differing values in the literature. We further aim to determine whether a possible luminosity-correlation of the CRSF line energy can account for them. Many cyclotron line sources are known to show such a correlation (see, e.g., Staubert et al. 2019). For GX 301−2 such a dependency was suggested by La Barbera et al. (2005).

As the luminosity and flux estimates stated in relevant publications naturally depend on the spectral model used and the energy range of the instrument, we resort to different means to estimate the source brightness. Both the Swift (Krimm et al. 2013) and the RXTE ASM (Levine et al. 1996) observatories have been monitoring GX 301−2 almost daily for more than 25 years (for an extensive analysis of the ASM light curve of GX 301−2 see Leahy 2002). A method to allow comparison between the observed count rates from the two observatories is given in Appendix A. We take the Swift count rate as a rough luminosity estimate and consider the variation of CRSF energy with respect to it. One caveat of this method is the strong assumption that the source spectrum remains constant over time, which is generally not justified. We note however that Staubert et al. (2014, 2016) employ the same method in their analysis of the long-term decay of the cyclotron line energy in Her X-1. We show the obtained data points in Fig. 12.

Fig. 12. Cyclotron line energy as a function of the Swift count rate. For the data points from this work we use the CRSF line energies as obtained from the phase-averaged analysis in Sect. 3.1. Where required, the ASM count rate was converted as Swift count rate as described in Appendix A. The uncertainties are given on a 90% level. The semi-transparent data points with dashed error bars do not satisfy the criteria described in Sect. 5.4, as they are either too wide or the line energy is too poorly constrained. The two gray areas represent regions where we find the cyclotron line energies to be believable. Data are from Ding et al. (2021), Nabizadeh et al. (2019), Fürst et al. (2018), Suchy et al. (2012), Doroshenko et al. (2010), La Barbera et al. (2005), Coburn et al. (2002), and Orlandini & Fiume (2001).

Figure 12 shows that no simple luminosity-correlation can be inferred. Moreover, we note that CRSFs have been found both within the known ranges of 30–40 keV and 50–55 keV as well as between the two intervals, where we note that CRSF energies are generally less well constrained. Furthermore, it appears that observations with NuSTAR and Insight-HXMT can constrain the CRSF line energy much better compared to observations with other instruments. Our analysis above showed correlations between CRSF and continuum parameters which raise doubts whether all data points shown in Fig. 12 correctly model an absorption component interpreted as a cyclotron line. Therefore, we only consider CRSF values that fulfill the following criteria:

1. The ratio between uncertainty in line energy and line energy ΔE_CRSF/E_CRSF must be at most 20%.

2. The width of the CRSF must be reasonably small. We estimate an upper limit of the CRSF width by using Eq. (6) for θ = 0, E_CRSF as stated in the literature, and kT_e = 6.50 keV as obtained in Sect. 3.1 for the folding energy, that serves as a rough estimate of the plasma temperature. We discard all CRSF values whose width exceeds this value by more than 20%.

The CRSF values that meet the criteria are shown in Fig. 12 as opaque data points with solid error bars, while values not meeting the criteria are shown as semi-transparent data points with dashed error bars. Our classification mostly removes CRSF values that are between the two well-known “bands”. A clear correlation with luminosity in terms of BAT count rate can however still not be inferred from our collected and filtered data points.

5.5. Models including cyclotron line emission wings

In the above sections, we discussed spectral models utilizing two cyclotron line components. There are however physical models that can potentially explain the observed X-ray spectrum of GX 301−2 using only one such feature. Schönherr et al. (2007) studied cyclotron line shapes in detail and state that CRSFs often exhibit emission wings. If emission wings are present in the spectrum, the neglect to model them can lead to wavy residuals, from which additional cyclotron lines could erroneously be identified. Recently, it was shown that for low mass-accretion rates, emission wings of a cyclotron line formed by resonant redistribution in the hot NS atmosphere can significantly modify a wide range of the X-ray continuum (Sokolova-Lapa et al. 2021; Mushtukov et al. 2021). This process has likely led to erroneous interpretation of a dip in the spectrum of Be X-ray binary X Persei as a cyclotron line (Coburn et al. 2001). Similar dips were shown to be present at intermediate energies in low-luminosity observations of GX 304−1 (Tsygankov et al. 2019b) and A 0535+262 (Tsygankov et al. 2019a), for which cyclotron lines are well-known at higher energies (see, e.g., Malacaria et al. 2015; Ballhausen et al. 2017). At higher mass-accretion rates, relevant for GX 301−2, the formation of emission wings likely has less dramatic manifestation, but redistribution during Compton scattering in a warm plasma should still contribute to their formation (Schönherr et al. 2007; Schwarm et al. 2017a). A similar nature can be attributed to the high-energy hump below the cyclotron line in the spectra of A 0535+262 (Ballhausen et al. 2017). In addition, it was recently shown by Loudas et al. (2024) that pronounced blue wings of cyclotron lines can be formed by bulk Comptonization in the radiative shock.

Here, we investigate the two possibilities that the complex structure of the continuum of GX 301−2 is created by either the 35 keV cyclotron line with a blue emission wing situated at ∼40 keV or the 50 keV CRSF with a red emission wing also at ∼40 keV. We therefore fit a model of the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f(E)=& \mathtt{detconst}\times \mathtt{PCF}(E)\times \mathtt{gabs}(E)\hfill \\ \hfill & \times (\mathtt{powerlaw}(E)\times \mathtt{FDcut}(E)+\mathtt{egauss}(E))\hfill \end{array}\end{array}$$(7)

to the 10–79 keV spectrum of GX 301−2. The parameters of the partial covering absorption component are set to those determined in Sect. 3.1. The spectra along with residuals are shown in Fig. 13, the fit parameters given in Table 6.

Fig. 13. Spectral fits of models including either the 35 keV cyclotron line with the blue emission wing or the 50 keV cyclotron line with the red emission wing. Top panel: 10–79 keV NuSTAR data and models. Lower panels: residuals of the shown models.

For the model of the 35 keV CRSF, the gabs component naturally describes the CRSF and the egauss component its blue wing. The centroid energy of the blue emission wing cannot be constrained, we therefore tie it to the CRSF parameters via E_Wing = E_CRSF + σ_CRSF. The second column of Table 6 shows that the area of the blue wing is consistent with zero and its width is unexpectedly small. Along with the fact that the residuals (Fig. 13, middle panel) still show a significant dip at 50 keV, we conclude that this model cannot describe the X-ray continuum of GX 301−2 satisfactorily.

The model including the 50 keV CRSF and its red wing on the other hand provides a better description of the data: it can constrain both emission wing and CRSF parameters and the residuals show no significant deviations from the data; we especially stress that no dip at 35 keV, the energy of the prominent cyclotron line of GX 301−2, is evident. Also in terms of the statistical evaluation of goodness of fit it provides an acceptable result. The strength of the CRSF of 36 keV is larger than found in the analyses in Sect. 3, which itself is however no ground for concern, as the models are built on different underlying physical ideas (multiple CRSFs versus one CRSF including emission wings). We note that the width of the CRSF of 12.7 keV is also high in comparison and at least raises some concern whether an actual absorption feature due to cyclotron resonant scattering is being modeled (see, e.g., Staubert et al. 2019, Table A.5 for typical CRSF widths). In order to follow up on this interesting and promising alternative interpretation of the observation better constraints on the high-energy emission of GX 301−2 are required.

6. Conclusion

In summary, our analysis of the NuSTAR observation of GX 301−2 has shown that

1. Two absorption features, most likely CRSFs, are present in both the phase-averaged and -resolved spectra. The line spacing is not harmonic and no potential fundamental line around 17 keV has been found. Both CRSFs are present throughout the whole observation.

2. The lower cyclotron line is present during the whole pulse period; the higher cyclotron line only during phases φ ∼ 0.2–0.6. The pulse-phase dependency of the lower line is in very good agreement with results by Fürst et al. (2018). In the phase-resolved analysis, both lines are generally more narrow and weaker than in the phase-averaged spectra.

3. Although the ratio between the line energies of 1.38 does itself not suggest a harmonic relation between the two lines, we cannot rule out this possibility on the grounds of multiple studies of other sources, where ratios below 2 have been found. We note however, that the ratio found here is the lowest reported for GX 301−2 to date.

4. Based on observations of GX 301−2 by various observatories over the last 20 years, no clear correlation of the cyclotron line energy with Swift count rate can be inferred. While there are also publications in which values for E_CRSF between 40–50 keV have been found, these have to be considered in more detail as it is possible that they describe two CRSFs with one model component. Only the analysis of archival data can allow to draw conclusions about these certain values.

5. An alternative interpretation, in which the hard X-ray spectrum is modeled with a powerlaw continuum and 50 keV CRSF with a red emission wing provides a good description of the data and calls for further investigation of its applicability to the X-ray spectrum of GX 301−2.

Considering the long history of the spectral study of GX 301−2, we note that further significant progress in disentangling cyclotron lines and continuum effects will require noticeably better constraints of the high-energy emission above ∼60 keV at intermediate flux levels. New missions such as the probe-class mission concept high energy X-ray probe (HEX-P) are expected to take on this particular challenge for accreting neutron stars in HMXBs (Ludlam et al. 2023).

Acknowledgments

We thank the anonymous referee for their very helpful comments that allowed us to improve the manuscript. We acknowledge support from Deutsche Forschungsgemeinschaft grants WI1860/11-2 and WI1860/19-1. The material is based upon work supported by NASA under award number 80GSFC21M0002. The authors thank Lorenzo Ducci, Victoria Grinberg, and Robbie McCulloch for their collaboration and insight on GX 301−2 and the XMAG collaboration, especially Kent Wood, for their helpful comments and productive discussions. This research has made use of ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT (https://www.sternwarte.uni-erlangen.de/isis/), of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center, and of NASA’s Astrophysics Data System. 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). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

Appendix A: RXTE/ASM-Swift cross calibration

In order to derive daily light curves for GX 301−2 for the past ∼3 decades, we use both Swift (Krimm et al. 2013) and RXTE/ASM (Levine et al. 1996) light curves binned to 1 d time resolution⁴. In order to combine the measurements from both observatories we aim to derive a simple method to relate their count rates. This is possible as there is a ∼6 year time window in which both satellites observed the source almost daily. In total, we find 1900 days with coverage by both observatories between MJD 53416–55925. As the individual data points of the daily light curves can be strongly influenced by background radiation and noise they are not suited to derive a relation between the two detectors. We therefore calculate the cumulative sum of all data points for both light curves as

$$\begin{array}{c}\hfill {\mathrm{CC}}_{\mathrm{ASM}}(T)={\displaystyle \sum _{i=1}^T}{C}_{\mathrm{ASM}}(i),\end{array}$$(A.1)

where C_ASM(i) is the count rate observed with the All-Sky-Monitor on day MJD = i + 53415 and T ∈ [1, 2509]. We show the obtained cumulative sum of the ASM as a function of the cumulative sum of the BAT monitor in Fig. A.1, together with a linear fit of the form

$$\begin{array}{c}\hfill {\mathrm{CC}}_{\mathrm{ASM}}({\mathrm{CC}}_{\mathrm{BAT}})=m·{\mathrm{CC}}_{\mathrm{BAT}}+t.\end{array}$$(A.2)

The best-fit parameters are given in Table A.1. The χ² fit statistic does not take the uncertainties of the input values into account, that are used to calculate the expected values. However, as the uncertainties of CC_BAT are very small in comparison to the corresponding quantity of the ASM, it is sufficient to not consider them for the fitting algorithm. We find that the fit describes the data very well, and that there are only very small deviations from the best fit.

Fig. A.1. Cross calibration between Swift and RXTE/ASM. Black represents the cumulative sum of the ASM count rates as a function of the cumulative sum of the BAT count rates. Red shows the fitted linear curve.

Using Eq. (A.2), we can determine the BAT count rate from measured ASM count rates via

$$\begin{array}{c}\hfill C_{\mathrm{BAT}}=m^{-1}·C_{\mathrm{ASM}}.\end{array}$$(A.3)

The parameter t is not required for this conversion as it corresponds to the - theoretically zero - offset between the two cumulative sums and bears no meaning for individual count rates. We can now use the above equation to estimate BAT count rates from ASM count rates for those MJDs where no BAT observation took place.

In order to judge how well this cross calibration works, we show both BAT and ASM-converted BAT count rates for the interval MJD 54000–54300 in Fig. A.2.

Fig. A.2. Light curve comparison between the BAT and converted ASM. Blue shows the light curve of GX 301−2 as observed with Swift. Red represents the ASM light curve converted to BAT-comparable count rates.

We note that the general orbital luminosity-dependency of GX 301−2 is apparent in the ASM-converted BAT rate, whereas there are a non-negligible number of outliers where the converted rate is significantly higher or lower than the BAT count rate. Overall, it is evident that the procedure can be used to roughly estimate the corresponding Swift count rate from an observed ASM count rate; however it is also apparent that there are some data points in Fig. A.2 for which the values determined by the two observatories differ significantly. This is not surprising since instrument effects are not taken into account in the procedure and since ASM and Swift cover different energy bands of 2–10 kev and 15–50 keV, respectively. In general, the light curve obtained with Swift shows far fewer outliers and in comparison shows a more periodic light curve that more closely resembles the expected periodic X-ray emission that stems from the orbital characteristics of the system (see, e.g., Islam & Paul 2014, their Fig. 1). This is consistent with the harder energy band being less affected by processes outside of the accretion column-like flaring or absorption. Therefore, the Swift light curve is in our case favorable over the ASM light curve. In conclusion, we obtained a light curve combined from the two observatories, which roughly resembles the X-ray luminosity of the source.

Appendix B: Gaussian soft emission

In Sect. 3.1 we introduced a Gaussian soft emission component in our spectral model. In this appendix, we briefly explain why this parameter is required. Our preliminary spectral analysis includes the application of different phenomenological empirical models such as FDcut, cutoffpl, and NPEX, the former two also in combination with blackbody components. Applying a comparably simple model of the form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f(E)=\mathtt{detconst}\times (& \mathtt{PCF}(E)\times \mathtt{continuum}(E)\hfill \\ \hfill & +\mathtt{egaus}{s}_{\alpha }(E)+\mathtt{egaus}{s}_{\beta }(E))\hfill \end{array}\end{array}$$(B.1)

Fig. B.1. χ residuals of the cutoffpl+bbody fit with and without an added Gaussian soft component. The component has a center position of ∼6.1 keV, a width of ∼0.9 keV and a normalization of ∼0.07 photons s−1 cm−2. The color scheme for NuSTAR’s focal plane modules is as in previous figures.

either with or without a gain shift shows very strong deviations from the data at energies of 6–7 keV. The structure of the residuals indicate an incorrect modeling of the iron line complex, which can partly be solved by allowing the Fe Kα line energy to shift to ∼6.1 keV. Such a shift is not expected from a physical point of view; further it does not fully solve the modeling issues. In order obtain a better empirical description of the spectrum in the range 6–7 keV, we introduce a Gaussian emission component centered at ∼6 keV. We show the residuals of a cutoffpl+bbody fit with and without the Gaussian emission component in the soft X-ray band in Fig. B.1.

We note that the additional soft emission component is very prominent and contributes around 20% of the total model counts at its centroid energy. Whilst the X-ray spectrum of GX 301−2 is known to exhibit a Compton shoulder (see, e.g., Watanabe et al. 2003, Fig. 1) at energies of ∼6.2 keV, the residuals shown in the upper panel of Fig. B.1 do not indicate that the Compton shoulder is responsible for the observed deviations from the model. We argue that both calibration issues and a complex continuum shape in the soft X-ray band lead to deviations from the model, which therefore needs to include the empirical Gaussian emission component. Furthermore, NuSTAR’s energy resolution of ∼400 eV at 6.4 keV (see Harrison et al. 2013) does not suffice to model the soft emission in detail. To model both the soft and the hard X-ray band simultaneously, coordinated contemporaneous observations of NuSTAR and XMM-Newton or Swift/XRT are required. In the present work, we see the inclusion of the additional soft Gaussian feature as a means to account for the complex soft X-ray emission structure of GX 301−2, which is not of interest in this work and does not influence our investigation of the cyclotron lines.

Appendix C: Detection significance of the 35 keV CRSF in the phase-resolved spectra

Based on the fact that the strength of the 35 keV CRSF in the phase-resolved spectral fits (see Fig. 9) is predominantly low and under ∼5 keV, the question arises with which significance the line is detected and whether even spectral fits without it can provide a good description of the data. We therefore employ the same statistical test as for CRSF2 as in Sect. 4 in order to determine the detection confidence for each phase bin, respectively. Analogous to Table 4, we give the results of this statistical test in Table C.1.

It is evident that the 35 keV CRSF can be detected with at least 3σ significance during the whole pulse phase. We consequently argue that CRSF1 is present throughout the whole pulse phase as found in many previous studies (see, e.g., Kreykenbohm et al. 2004; Fürst et al. 2018)

All Tables

All Figures

Fig. 1. Orbital geometry of the system GX 301−2 based on Fig. 1 of Fürst et al. (2011). We show the NS orbit around its B hypergiant companion Wray 977 along with the orbital Swift (Krimm et al. 2013) count rate as determined from the long-term BAT light curve in blue. We further show the orbital phases of NuSTAR observations of GX 301−2 and the corresponding year in red. The 2022 observation is omitted as it was performed at nearly the same orbital phase as the 2014 observation. For the orbital parameters, we used the values presented by Kaper et al. (2006), that is, a sin i = 368.3°, e = 0.462, ω = 310.4°, and RWray 977 = 62 R⊙. We use an inclination of i = 55° (within the range given by Leahy & Kostka 2008). Epoch folding of the long-term Swift light curve yields an orbital period of 41.474 ± 0.006 d, which is in good agreement with the period stated by Doroshenko et al. (2010) as derived from Koh et al. (1997), including an orbital period derivative. We subsequently use the periastron passage TPA = 53531.65 MJD from Doroshenko et al. (2010). For the calculation of orbital phases, we take the orbital period derivative Ṗorb = −3.2 × 10−6 s s−1 (Doroshenko et al. 2010) into account.

In the text

Fig. 2. NuSTAR observation of GX 301−2 in the context of its long-term behavior. a: hardness ratio, (H − S)/(H + S), of NuSTAR observation with hard energy range 10–30 keV and soft energy range 4–10 keV. The dashed lines indicate the division into time segments studied individually in Sect. 3.2. b: light curve of NuSTAR observation. c: Swift light curve between MJD 59100 and MJD 59300 (Krimm et al. 2013). We show ∼5 orbital periods and indicate the pre-periastron flares with orange bars.

In the text

Fig. 3. Continuum modeling of GX 301−2. a: spectra of the two NuSTAR FPMs in blue and purple with best fit model shown in black. Models with both CRSFs removed in turn and model without CRSFs are shown in green, red, and orange. b: χ residuals of best fit. c and d: residuals of best fit with both CRSFs removed in succession. e: best fit of alternative model with only one CRSF. f: best fit of alternative model without CRSFs. The spectra have been rebinned for plotting purposes.

In the text

	Fig. 4. Ratios between spectra of the three time segments versus the spectrum of the whole observation. The dashed line indicates the behavior expected in the case of no spectral change.
In the text

	Fig. 5. χ residuals of best fits to time-segmented spectra.
In the text

	Fig. 6. Behavior of GX 301−2 over the pulse phase. Top panel: pulse profiles for different energy bands. Bottom panel: hardness ratio (H − S)/(H + S) with soft energy range 4–10 keV and hard energy range 10–30 keV.
In the text

	Fig. 7. Phase energy map for energies between 3.5–60 keV. For each energy bin, we show (N − ⟨N⟩)/(Nmax − Nmin).
In the text

	Fig. 8. Pulsed fraction as a function of energy. We define the pulsed fraction as (Nmax − Nmin)/(Nmax + Nmin), where Nmax/min are the maximum and minimum rates of the pulse profile in a given energy bin.
In the text

	Fig. 9. Phase dependencies of all free parameters in the phase-averaged analysis. Red represents the phase-averaged parameter values. Gray represents the normalized pulse profile. Arrows indicate where uncertainty intervals hit their predefined range limits.
In the text

	Fig. 10. χ residuals of fits to the phase-resolved spectra. We use the same color-coding of NuSTAR focal plane modules as in above figures.
In the text

	Fig. 11. Phase dependencies of both cyclotron lines as determined in this work (blue circles) and by Fürst et al. (2018, black squares). We also show our phase-averaged values in red.
In the text

Fig. 12. Cyclotron line energy as a function of the Swift count rate. For the data points from this work we use the CRSF line energies as obtained from the phase-averaged analysis in Sect. 3.1. Where required, the ASM count rate was converted as Swift count rate as described in Appendix A. The uncertainties are given on a 90% level. The semi-transparent data points with dashed error bars do not satisfy the criteria described in Sect. 5.4, as they are either too wide or the line energy is too poorly constrained. The two gray areas represent regions where we find the cyclotron line energies to be believable. Data are from Ding et al. (2021), Nabizadeh et al. (2019), Fürst et al. (2018), Suchy et al. (2012), Doroshenko et al. (2010), La Barbera et al. (2005), Coburn et al. (2002), and Orlandini & Fiume (2001).

In the text

	Fig. 13. Spectral fits of models including either the 35 keV cyclotron line with the blue emission wing or the 50 keV cyclotron line with the red emission wing. Top panel: 10–79 keV NuSTAR data and models. Lower panels: residuals of the shown models.
In the text

	Fig. A.1. Cross calibration between Swift and RXTE/ASM. Black represents the cumulative sum of the ASM count rates as a function of the cumulative sum of the BAT count rates. Red shows the fitted linear curve.
In the text

	Fig. A.2. Light curve comparison between the BAT and converted ASM. Blue shows the light curve of GX 301−2 as observed with Swift. Red represents the ASM light curve converted to BAT-comparable count rates.
In the text

	Fig. B.1. χ residuals of the cutoffpl+bbody fit with and without an added Gaussian soft component. The component has a center position of ∼6.1 keV, a width of ∼0.9 keV and a normalization of ∼0.07 photons s−1 cm−2. The color scheme for NuSTAR’s focal plane modules is as in previous figures.
In the text