Evidence for a nearly orthogonal rotator in GX 301–2 with phase-resolved cyclotron resonant scattering features

Cyclotron resonant scattering features (CRSFs) are the absorption features in the X-ray spectra of strongly magnetized accretion neutron stars (NSs), which are probably the most reliable probe to the surface magnetic fields of NSs. The high mass X-ray binary GX 301–2 exhibits a very wide, variable and complicated CRSF in the average spectra, which should be two absorption lines based on NuStar and Insight-HXMT observations. With the Insight-HXMT frequent observations, we performed the phase-resolved spectroscopy and confirmed two cyclotron absorption lines in the phase-resolved spectra, with their centroid energy ratio ∼ 1 . 6 − 1 . 7 in the super-critical luminosity case. A major hindrance in understanding those CRSFs is the very poorly constrained magnetic inclination angle, which is also a fundamental property of a NS and key to understanding the emission characteristics of a pulsar. Comparing the phase-resolved CRSF with simulated X-ray spectra, the magnetic inclination angle is found to be ≳ 70 ◦ , i.e., nearly orthogonal between the NS’s spin and magnetic axies. The implications of an orthogonal rotator and magnetic structure evolution in the accreting X-ray binary are also discussed


INTRODUCTION
Neutron stars (NS) are ideal astrophysical laboratories for testing theories for dense materials, gravity and strong magnetic field physics (Özel & Freire 2016).Rotating NS, namely pulsar, emits radio, X-ray and multiwavelength pulsed radiation, and provides a tool to probe the NS characteristics (Manchester 2004).The magnetic inclination as an angle between the magnetic and spin axis is a key parameter for understanding the radiation mechanism, but is not well understood yet (Beskin et al. 1993;Novoselov et al. 2020), especially for the accretion X-ray pulsars, which we still have no clear observational information (Biryukov & Abolmasov 2021).Cyclotron resonant scattering features are the absorption features in strongly magnetized accretion neutron stars (Truemper et al. 1978;Staubert et al. 2019), which can show the accretion processes near the magnetic poles of neutron stars and the magnetic field geometry (Lamb et al. 1973;Harding & Daugherty 1991;Schönherr et al. 2007;Nishimura 2011;Mushtukov & Tsygankov 2022).
The cyclotron resonant scattering feature of GX 301-2 around 40 keV was first discovered in Ginga observations (Makishima K. 1992), and then confirmed by other missions (Kreykenbohm et al. 2004; La Barbera et al. 2005;Doroshenko et al. 2010).However, the absorption feature is very broad, and its centroid energy varies wildly from ∼ 30 − 55 keV in different observations (Suchy et al. 2012), which makes the physical origin as a longstanding puzzle (Ding et al. 2021).With high resolution spectroscopy by NuSTAR, it has been argued that the reported broad feature is in fact from two separate line spawning regions located at the NS surface and 1.4 km above, generating two cyclotron absorption lines at  1 ∼ 35−40 keV and  2 ∼ 50 − 55 keV respectively (Fürst et al. 2018).The recent Insight-HXMT data confirmed two absorption lines in the average spectra, with a line ratio of ∼ 1.7 based on the multiple observations (Ding et al. 2021).
GX 301-2 has been observed by the Insight-HXMT for multiple times from August 2017 to December of 2020, showing the variable X-ray light curves over time and orbital phases (Ding et al. 2021).Phase-average spectroscopy based on Insight-HXMT shows that the centroid energies of both two CRSFs vary with X-ray luminosity, discovering the presence of critical luminosity and accretion state transition in GX 301-2 (Ding et al. 2021).
Modeling the structure of accretion column in accreting X-ray NS is of fundamental importance as it would answer how the materials are transferred into radiation around the NS, and could potentially unveil the inner structure and evolution of NS.However, such a task is till now a formidable problem due to the coupling of radiation transfer and magnetic fluid dynamical processes in the extremely strong magnetic field in NS atmosphere.There have been several attempts to model the radiation from NS accretion column (Schönherr et al. 2007;Nishimura 2011;Nishimura 2015;Schwarm et al. 2017;Nishimura 2008Nishimura , 2014)).To reduce the intractable complexity in solving the time-dependent radiation transfer and magnetic hydrodynamical equations, a common strategy in these works is to adopt an analytic super-critical accretion column solution, which assumes a magnetically confined ideal fluid impact on the NS magnetic polar, producing radiation dominated shock (Nishimura 2008;Becker 1998;Arons 1992;Basko & Sunyaev 1976).The velocity distribution along the vertical axis results in the Doppler boosting on the generated CRSFs, thus producing quite different absorption structure when changing the viewing angle (Nishimura 2015(Nishimura , 2022)).
In this work, the Insight-HXMT observation and the data reduction of the phase-resolved spectrum are described in Section 2.Then, the simulation process is presented in Section 3, including the spectral modelling and the generation of the simulated phase-resolved spectra.In Section 4, the results of both observational and simulated phase-resolved spectroscopic data are presented.The dependency of CRSF properties on pulse phases and magnetic configuration, which provides a unique method to constrain the magnetic inclination of NS, is discussed in Section 5. Finally, we summarize our results in Section 6.

Insight-HXMT
Insight-HXMT has three main payloads: the Low Energy X-ray Telescope (LE, Chen et al. 2020), the Medium Energy X-ray Telescope (ME, Cao et al. 2020), and the High Energy X-ray Telescope (HE, Liu et al. 2020).HE contains 18 cylindrical NaI(Tl) / CsI(Na) detectors with a total detection area of 5000 cm 2 in the energy range of 20-250 keV.The time resolution is 25 s, and its energy resolution at 60 keV is 19% (Zhang et al. 2020).ME is composed of 1728 Si-PIN detectors with a total detection area of 952 cm 2 covering the range of 5-30 keV.Its time resolution is 280 s, and the energy resolution at 20 keV is 14%.LE uses Swept Charge Device (SCD) with a total detection area of 384 cm 2 covering 1-15 keV.It has a time resolution of 1 ms and the energy resolution at 6 keV is 2.5%.The typical Field of Views (FoVs) for the three detectors are 1.1 (Zhang et al. 2020).Insight-HXMT has carried out a series of performance verification tests by observing blank sky, standard sources, showing good calibration state and estimation of the instrumental background (Li et al. 2020).
Insight-HXMT satellite conduct intensive recurring pointing observations on GX 301-2 between August 3   , 2017 and June 3   , 2020 (from MJD 57968 to 59003), with typical net exposure time in each pointing ∼ 3 ks.These observations cover almost all orbital phases of the binary orbit.The observations (information of OBSIDs) used in this paper are the same as those in Ding et al. (2021), and they are listed in Table 1, along with the observation date, exposure, count rate, pulse period, orbital phase, and luminosity.In the following science analysis, we filtered data with the recommended settings in HXMT Data Reduction Guide: (1) pointing offset angle < 0.1 • ; (2) pointing direction above Earth > 10 • ; (3) geomagnetic cut-off rigidity value > 8 GeV; (4) time since SAA passage > 300 s and time to next SAA passage > 300 s; (5) for LE observations, pointing direction above bright Earth > 30 • .The Insight-HXMT Data Analysis (HXMTDAS) V2.04 was used to produce high level products.We made the barycentric correction of the data using the tool hxbary in HXMTDAS.Well-calibrated energy bands of LE, ME and HE were used in the following spectra analysis (3 − 8.5 keV, 10 − 30 keV and 30 − 70 keV, respectively, Li et al. 2020).ME data in 21 − 24 keV was ignored due to the significant increase of calibration uncertainties in this energy range.

Phase-resolved spectral analysis
The photon arrival time for binary motion is first corrected using the latest ephemeris provided by Doroshenko et al. (2010).For each Insight-HXMT observation, we derive the spin period of the neutron star of GX 301-2 using efsearch which folds the light curve and finds the best period by searching for the maximum chi-square.The error is then calculated with   =

√
2  /( √  ) provided by Larsson (1996) because of the uneven sampling rate, where  is a simulation constant of 0.469,  is the total data points,  is the sinusoidal amplitude, and  is the time duration.The pulse period and error are shown in Table 1.To improve accuracy, pulse periods shown in Table 1 are determined using combined data where observations taken within two days are merged together.In general, a period of approximately 685 seconds would result in only four pulses being observed during a continuous exposure time of around 3 ks.However, the typical total time span of each Insight-HXMT observation is approximately 9 ks, with the 3 ks exposure time resulting from the selection of Good Time Intervals (GTI).This means that each of our observation covers more than 10 pulses, which is sufficient to obtain a reasonably accurate pulse period.Pulse profiles can then be obtained.Here we show an example of pulse profiles of P010130900103 in different energy bands covering the energy range that used in this paper, see Figure 1.More detailed pulse information can be found in Ding et al. (2021).The observed spectra of the accretion X-ray pulsar would vary in different pulse phases, which can be used to infer the accretion process of the emission region and geometric characteristics.In each pointing observation, we calculate the corresponding pulse phase of each corrected time based on the derived pulse period, setting initial phase to be the periastron passage (Doroshenko et al. 2010).By dividing the phase into eight equal width parts (see Section 4.1), the observation time can also be divided into the corresponding eight parts, thereby the corresponding phase-resolved spectra can be extracted, with net exposure time of ∼ 150 − 300 s.For each phase-resolved spectrum, we use the HXMT software to regenerate the new background spectrum and response matrix.
Different phenomenological models with a power law modified by cut-off at higher energy can describe the broad band spectrum of GX 301-2 very well (Endo et al. 2002;Fürst et al. 2018;Ding et al. 2021).At lower energy the spectrum gets heavily absorbed due to the presence of a variable neutral absorber (Fürst et al. 2011), with the Hydrogen column density in the range of 10 23 −2×10 24 cm −2 along the observer's line of sight throughout the orbit, and thus indicating accretion mode from clumpy stellar wind.To explain the low-energy part of the spectrum, previous studies applied a partial covering component in addition to the usual galactic absorption (La Barbera et al. 2005;Mukherjee & Paul 2004).In 6 − 7 keV, GX 301-2 shows many fluorescent lines (Fürst et al. 2011;Ji et al. 2021) including the bright and prominent Fe K shell emission lines (Endo et al. 2002).
In this work, the K and K lines are fitted with Gaussian and are frozen at 6.4 and 7.0 keV respectively.Their widths are frozen at 10 eV.
The spectral fitting in this paper is conducted with Xspec v12.12.0g (Arnaud 1996).We conducted a phenomenological fitting similar to Ding et al. (2021) with the spectral model in Xspec's terminology as Constant*TBabs*(TBpcf*CRSFs*Continuum+Gaussian+Gaussian), where CRSF is represented using the gabs model.

Spectral Modelling
A prominent accretion mound, inside which is a sinking region, is expected to be formed above  ∼ 1 × 10 37 erg s −1 (Arons 1992) in which the line-forming region is expected to be located around the accretion mound.In this work, we consider a line-forming region around an accretion mound with a two-dimensional structure taking into account gravitational light bending under the assumption of two symmetrically located polar caps like Nishimura (2022).
Assuming a radiation dominated accretion column, we conduct a Monte Carlo (MC) simulation to deal with the radiation transfer process inside the accretion column.We adopt a model derived by Arons (1992) for an accretion mound shape as shown in Figure 2.For 0 ≤  ⊥ <   , the equality between the ram and radiation pressure yields the nominal location of the shock to be at the height where  ⊥ is the cylindrical radius with respect to the magnetic axis,   =   /   is the accretion luminosity of one pole where  is the gravitational constant,   is Eddington luminosity. ∥ ∼  ⊥ ∼ 1 is a good approximation over most of the mound's volume (Arons 1992).We consider a typical supercritical luminosity case with  = 2 × 10 37 erg s −1 , implying a shock height of ∼ 1 × 10 3 m (Arons 1992).In the critical luminosity, the mound is expected to be as high as   = 10 3 m (Basko & Sunyaev 1976).
The line-forming region is assumed to be a narrow region around an accretion mound with a height of ∼ 1 × 10 3 m, thus we actually compute the cyclotron line in altitude from ∼ 10 m to ∼ 1.1×10 3 m in an accretion column with the surface cyclotron energy being assumed as 75 keV, such that the re-constructed phase-averaged spectra would produce two cyclotron lines at ∼ 33 keV and ∼ 55 keV respectively (coincide with observations, see Ding et al. 2021).The strength of the surface magnetic field can be about 6.5 × 10 12 G when gravitational redshift is taken into account.
The bulk velocity is expected to be decelerated due to resonant scattering in the line-forming region as shown in Figure 2. We assume the bulk velocity profile in the line-forming region is given by where  is the distance from the surface of a NS along the column axis,    () = − √︁ 2    /(   + ) is free-fall velocity at a height , () is decelerated from −   (  ) to − 0 in the line-forming region where  0 = 0.1 is the velocity at the surface of an accretion mound, making use of an approximate velocity profile  ∝  ∝ √  (Becker & Wolff 2007).Here,  is the speed of light.
We consider a beam pattern that have a beam with a peak intensity of  1 = 60 • , a pencil-like beam, or a beam pattern composed of two beam patterns that have a beam with a peak intensity of  1 = 60 • and with a peak intensity of  2 = 90 • , a fan beam.Here we consider the angle-dependence of intensity of the continuum spectra, i.e., beam patterns, as follows (Kraus 2001): Here,  is the angle between the direction of photon's propagation and B-field.We use 60 • for parameter  1 and 90 • for  2 respectively and fix 30 • for parameter  as in previous theoretical papers (Nishimura 2015), since a pencil-like plus fan beam is considered in this work.The radiation spectrum energy distribution (SED) is modeled by NPEX, which performs well both in the phase-averaged spectra (Ding et al. 2021) and in the phase-resolved spectra (see Section 4.1), being produced with two cut-off power-law components.
As previous study has confirmed the lack of variability for the continuum parameter throughout the NS orbital motion, we set these pa- rameters to be the mean values observed in previous phase-averaged study.Specifically, the low energy power-law component (Wien tail) is assumed to be a thousand time stronger than the hard component, while the photon index and high energy cut-off are frozen at 1.2 and 7.0 keV respectively.Eight derived spectra, which obtained by uniformly distributing cos  from 0 to 1, are shown on the right panel of Figure 3.

Generation of the simulated phase-resolved spectra
A photon emitted at an angle  with respect to the radial direction would escape to infinity at a larger angle Ψ due to general relativistic (GR) light bending near a compact object (Beloborodov 2002): where   is the Schwarzschild radius of the NS, while  is the NS radius.We assume that the tidal stream from Wray 977, and the generated transient disk still enables significant alignment between NS spin and binary orbital angular momentum, thus the observer inclination with respect to the NS should be close to the binary inclination ( ≃ 66 • , Kaper et al. 2006), implying a geometric relation connecting the observed spin phases (Φ) and the magnetic inclination of NS (see Figure 4 for the definition of these angles): Θ AC here is the angle between accretion column and NS spin axis, being assumed to be identical to the magnetic inclination; Φ 0 is an arbitrary phase offset.To deal with the problem of accretion columns being occluded by the NS, we assumed that as the accretion column crossing the limb of the star, disappearing from our line-of-sight, it quickly darkens and the polar cap at the opposite side will become dominant (i.e.we did not take the height of the column into account).Consequently, each combination of phase and Θ AC can be used to calculate the corresponding cos .Then, if a specific NS inclination is assumed, we could re-construct both the phase-averaged and phaseresolved spectra from the corresponding photon emission angle in any specific spin phase.Specifically, we assume the NS inclination is 66 • .Same as the observational phase-resolved spectra, we set eight uniformly distributed phases.Additionally, we set Θ AC to be a arithmetic sequence ranging from 0 to 90 with an interval of 5.This allows us to generate a total of 8 * 19 spectra corresponding to different combinations of phases and Θ AC .For each spectrum, we calculate the intensity by dividing its phase into eight equal width parts again (which gives a total of 64 uniformly distributed phases by now), and combining them with the Θ AC of the spectrum to calculate cos , determining which simulated spectrum (i.e. on the right panel of Figure 3) each combination corresponds to.Subsequently, we average the eight simulated spectra to obtain the final simulated phase-resolved spectrum for each combination of phase and Θ AC , with a total simulated phase-resolved spectra number of 8 * 19.
We then utilized the generated spectra by convolving them with authentic Insight-HXMT ancillary and response matrix, so as to show how these model spectra would appear when being observed with Insight-HXMT.To be more specific, we fit the simulated spectra with NPEX, while keep adding extra Gaussian absorption until no apparent structure is observable.After that, we loaded the ancillary and response file of Insight-HXMT and used the command fakeit built-in Xspec to produce fake spectra.The fake spectra are renormalized so as to produce energy flux ∼ 1 × 10 −8 erg cm −2 s −1 in 10 − 70 keV, being comparable to radiation dominated state's value if consider the source distance.The generated HE spectra of 17 detectors are then combined using the same method provided in HXMTDAS v2.04.The resulting spectra are again fitted with NPEX models to derive the CRSF energy ratio.Such ratio, together with the centroid energy, help us compare the theoretical line profile with the observational one in the supercritical luminosity scenario.On the left panel of Figure 3, we show the centroid energy ratio of the two CRSFs fitted from each simulated spectrum in a given photon emission angle.The result shows that the photon emitted at  ≲ 40 • (cos  > 0.65,  is the emission angle between photon and magnetic field) should be dominant in our observed spectrum.

Observational spectrum
As it is currently not practical to generate continuum from first principle in the context of accreting NS radiation simulation, we tested three different empirical continuum models first.highEcut is a high energy cut-off model, forcing the part above cut-off energy to decay exponentially: However, due to the discontinuous derivative, it will produce an artificial cusp at  cut , which was modeled in this work with a Gaussian absorption.newHcut (Burderi et al. 2000) uses polynomials to smooth the part around  cut , but introduced extra parameters like smooth range Δ, which is fixed to 5 keV here.NPEX (Mihara 1995) is a combination of two cut-off power law models with common  cut , having a soft component (photon index fixed to 2) and a hard component, which generally outperforms other two models when being applied to averaged spectrum (Ding et al. 2021).Besides these empirical continuum model, we also include a theoretical model for comparison, namely Becker-Wolff self-consistent cyclotron line model (i.e.bwcycl in XSPEC Becker & Wolff 2007;Ferrigno et al. 2009).Upon comparing the phase-resolved fitting results of these four models, we note that NPEX model outperforms the others, exhibiting more stable and normal parameters and yielding better chi-square values, especially when comparing with the derived parameters from highEcut and newHcut.On the other hand, bwcycl also performs good with its CRSF parameters very similar to those from the NPEX model, except for some phase-resolved spectra where the derived parameters of CRSFs and bwcycl differ significantly from those of other phaseresolved spectra in the same observation.Therefore, considering the quite stable results compared to others, NPEX has been chosen as our final continuum model.Besides, the final results of the centroid energy ratio of all the observations derived from the other three models are very similar the NPEX results, i.e., scattered below the critical luminosity and concentrated above it, with the best-fit ratios within the range of ∼ 1.6−1.8 as shown in Figure 5.
The fitted spectral parameters of the phase-resolve spectral analysis for P010130900103 is presented in Figure 6, which is a super-critical case of  = 4.2 × 10 37 erg/s.To confirm the existence of two CRSFs, residuals of models containing different number of gabs are shown on the right panel, with the corresponding reduced chi-square shown on the right side outside the panel, where the grey lines show the residuals from no gabs model, and the purple lines show the residuals from one gabs model.See the caption in Figure 6 for a detailed description of the figure elements.Combining both the reduced chisquare and the residual structure, it is evident that the addition of each CRSF component noticeably improves the fitting results of the spectrum.A lower luminosity observation with  = 6.6 × 10 36 erg/s is shown in Figure 7 for comparison.Although the fittings are still acceptable in lower luminosity cases, the results for several phases are not as good as in the super-critical observations.However, it is still evident that the inclusion of two CRSF models is necessary for the fitting in most phases, even though sometimes a second CRSF component only marginally improves the results.
We present the evolution of the fitted parameters with their pulse phases in Figure 8.For the higher luminosity case of P01030900103, the energy of the second CRSF remains almost constant, with a slight decrease around the pulse minimum and an increase at the beginning of the second pulse, which resembles the result of the 2015 outburst observed by NuSTAR (Fürst et al. 2018), even though these variations are not significant considering the error bars.The CRSF at lower energy exhibits a subtle decreasing followed by increasing with phase, and the continuum parameter Γ shows an inverse correlation with phase amplitude.All these results indicate a similarity in spectrum compared with the 2015 outburst, although the correlations presented by Fürst et al. ( 2018) is more pronounced.Meanwhile, the parameters of the lower luminosity observation showing on the right side of Figure 8 exhibit different behaviors.The CRSF at higher energy decreases at the second pulse, and the CRSF at lower energy decreases when the source is around its pulse minimum.The photon index shows a positive correlation with phase amplitude.These different behaviors may suggest that spectrum below and above critical luminosity could have different properties, although they could also be attributed to the observational precision due to lower brightness.It is worth noting that the 2015 NuSTAR observation analyzed by Fürst et al. (2018) has a luminosity of around 2.8 ×10 36 erg/s, which resembles our lower luminosity case.These fitting results are also shown in Table 2 with their chi-square and degree of freedom presented.In addition, observation P010130900107 is also included in the table for comparison, which has a luminosity of 1.8 × 10 37 erg/s.The errors reported in this paper are generated with Markov Chain Monte Carlo (MCMC) method (50 walkers, 200000 steps), and are at 90% confidence level except for special statements.
To avoid the influence of inaccurate fitting on the results, we do not consider phases with fitted absorption line energies greater than 70 keV, which is beyond the energy range of our spectra.After that, the best-fit ratio ( 2 / 1 ) of two cyclotron line centroid energies is   Table 2. Best-fit parameters for three observations in different spin phases, where P010130900103 has high luminosity of 4.2 × 10 37  /, P010130900107 has median luminosity of 1.8 × 10 37  /, and P010130900401 has low luminosity of 0.7 × 10 37  /. (1) (2)  derived by Orthogonal Distance Regression (ODR) so as to consider errors in both variables.The samples are shown in Figure 9, where the upper two panels are for the super-critical cases and the lower two panels are for the low luminosity cases.We can see that the upper ones are more concentrated with smaller errorbars, while the lower ones have large error bars and are scattered.

Simulated spectrum
It has been argued that the complex CRSFs in GX 301-2 may not result from the presence of two separate line forming regions, but entail multiple line forming regions (Ding et al. 2021).Specifically, recent theoretical models could naturally produce such absorption profiles (Nishimura 2015;Schwarm et al. 2017).Using NPEX model as the input spectra and a fan-like input beaming pattern, we can successfully reproduce the observed structures, which are presented on the right panel of Figure 3 with different emission angles.To illustrate how they appear when observed by Insight-HXMT, they are also convolved with authentic Insight-HXMT ancillary and response matrix.
They are then fitted with NPEX model plus two CRSF components, where the NPEX parameters are fixed to the values used to generate the spectra as discussed in Section 3.1.Figure 10 and Tabel 3 show the fitted results.
The resulted ratio in phase-averaged spectra is in good agreement with observation (∼ 1.6, but has no apparent variation with magnetic inclination), while the re-constructed phase-resolved spectra show strong dependency on magnetic inclination.In Figure 11, we show the fitting results of the simulated spectra, where   is the difference between the ratio of two cyclotron lines from the simulated spectra and the observed ratio obtained in Figure 5.The result shows that the observed photon emission angle is seriously restricted due to the magnetic geometry of the NS and the GR effect.We found that to obtain the CRSF structure in the phase-resolved spectral observations, the NS magnetic inclination should take special values.Assuming different magnetic inclination angles, the phase-resolved spectral simulations and fittings are performed to derive the cyclotron line centroid energy ratio.The results imply that the magnetic inclination of the NS should have the solution of ≳ 70 • , i.e. a nearly Table 3. Best-fit parameters with uncertainties for the MC simulated spectra as shown in Figure 10.The NPEX parameters are not shown in the table because they are fixed to the values used to generate the simulated spectra, i.e. the photon index and high energy cut-off are 1.2 and 7.0 keV respectively.orthogonal rotator, with 68% confidence level.The centroid energy ratio between the two CRSFs in GX 301-2 is the result of a nearorthogonal magnetic pole relative to the spin axis of the NS.

Phase-resolved spectra from Insight-HXMT
With frequent Insight-HXMT observations, we performed phaseresolved spectral analysis of GX 301-2, confirming the existence of two CRSFs in all the phases independent of the continuum spectrum models.The centroid energy ratio ( 2 / 1 ) between the two CRSFs with the X-ray luminosity is plotted (Figure 5), discovering the ratio  2 / 1 ∼ 1.6−1.7 in the super-critical state ( ≳ 1.8×10 37 erg s −1 ).The absorption structures in the phase-resolved spectra are modeled by Gaussian absorption components while the continuum is described with the empirical models, i.e., the Negative Positive Exponential (NPEX, Mihara 1995) model.The centroid energy ratio has a large scatter from 1.2 − 2 in the low luminosity, while reaches a narrow value range of ∼ 1.6 − 1.7 in high luminosity, and this value is well consistent with the observed ratio in the averaged spectra (Ding et al. 2021).This is remarkable because it means that the fixed energy ratio observed previously in phase averaged spectra is not due to averaging between different absorption profile, but is the consequence of the appearance of similar spectral structures in different spin phases.Such dependency appears most prominently in high luminosity state, with the critical luminosity possibly being the turning point (∼ 1.8 × 10 37 erg s −1 , Ding et al. 2021).It is currently not clear whether this is predominantly due to the more stable radiation dominated accretion column structure above the critical luminosity, because Insight-HXMT is not capable of fully resolving the CRSF structure in the phase-resolved spectra at low luminosity (≲ 10 37 erg s −1 ).For observations with energy ratio close to unity where it is not easy to separate one line from the other, the uncertainties of the fitted absorption energies at low luminosity are large compared to those above the critical luminosity, resulting in larger errorbars in the energy ratio as shown in Figure 5. Thus, the scattering and decrease of energy ratio could be the consequence of degeneration between two Gaussian absorption components in low signal-to-noise ratio.
We may note that the first CRSF is in the energy range of ∼ 30− ∼ 45 keV considering the uncertainties as shown in Table 2.For a first CRSF at 45 keV, the upper limit of 70 keV chosen in our spectra analysis can cover only an energy ratio of 1.5, which may seriously affect the final results.However, the HE detector of Insight-HXMT can reach 250 keV at its maximum (e.g.Ma et al. 2021), and has been used for CRSF research up to 100 keV (Liu et al. 2022).Consequently, the spectra fitting of 3-100 keV has been tested.The energy ratio remains stable compared to the 3-70 keV fitting results for high luminosity observations, and fluctuated for low luminosity observations.Most importantly, the centroid energy of the second CRSF do not exceed 70 keV even for the low luminosity cases (not include the uncertainties), and the energy ratios of two CRSFs for all observations exhibit similar characteristics: scatter in the low luminosity, but concentrate in the high luminosity, with the best-fit ratio of the super-critical state derived to be 1.64±0.02.However, the parameters derived from 3-100 keV is unstable with larger uncertainties, and the spectra above 70 keV is noisy, following the averaged spectra analysis performed in Ding et al. (2021), we keep the 3-70 keV results.
To check if the results are affected by different separated phases, we recalculated the centroid energy ratios in phase-resolved spectra with different pulse phases, i.e., 0.15-0.4,0.4-0.6,0.6-0.9,0.9-1.15,which separate the pulse phases roughly based on the first peak, pulse minima and second peak.The results are shown in Figure 12.Despite minor differences existed as compared with Figure 5, the plot is still consistent with it, i.e., scattered below the critical luminosity and concentrated above it, and the best-fit ratio is around 1.67, which is within the error bar of our 8-phases results.We need to point out that only ratios with more than four successfully fitted phase-resolved spectra in each observation are shown in Figure 5 for more reliable We could observe these photons with a given magnetic axis (y-axis) in different spin phases (x-axis).The simulated spectra in Figure 3(b) are integrated into each phase bin in this map.After assuming a specific magnetic inclination, we are able to use the simulated spectra to re-construct phase-averaged and phase-resolved spectra.(b) The generated input pulse profile when magnetic inclination is 85 • .(c) Fitting results (ODR) of simulated phase-resolved spectra with different magnetic inclinations, where   is the difference between the ratio of two cyclotron lines from the simulated spectra and the observed ratio obtained in Figure 5.(d) Similar with (c) except the ratio from the simulated spectra is calculated using the average instead of ODR.All errors indicated in this figure are in 68% confidence level (1).Spectral simulations imply that the NS is a nearly orthogonal rotator with a magnetic inclination Θ AC ≳ 70 • .results, and ratios fitted from less phases as shown in Figure 12 may raise inaccuracy.
Combining multiple spectra may increase the confidence level of the results, thus we select observations with similar luminosity, and combine their phase-resolved spectra for each of the eight phases.The final results are shown in Figure 13.On the left panel, we combine three observations of luminosity ∼ 0.35 × 10 37 erg/s, and the right panel contains the combined results of four observations with their luminosity ∼ 1.0 × 10 37 erg/s.Before combining the spectra, the energy uncertainties for different phases in these selected observations were quite large, and the energy distribution was scattered.Due to the low luminosity, some of these observations only had one or two successfully fitted phases, resulting in randomly distributed ratios ranging from ∼ 1.1 to ∼ 1.5.However, as shown in Figure 13, the fitting results are significantly improved with much smaller uncertainties after combination, and the data points become more con-centrated.The fitted ratios are both around 1.3, much lower than the best-fit ratio of the super-critical state.The enhanced spectra results again demonstrate the significant difference in the energy ratio above and below the critical luminosity.

Magnetic inclination angle
The CRSF is highly anisotropic due to the uneven velocity distribution and magnetic field strength inside the accretion column (Becker 1998;Nishimura 2015).Consequently, the cyclotron absorption lines can provide important information about the material distribution near the NS magnetic polar region.For GX 301-2, two reported cyclotron absorption lines  1 ∼ 30 − 42 keV and  2 ∼ 48 − 56 keV show the similar relationships of the centroid energy versus X-ray luminosity (Ding et al. 2021), suggesting a state transition at a critical luminosity of ∼ 1.8 × 10 37 erg s −1 .To further investigate the Centroid energy ratios in phase-resolved spectra between the two CRSFs of GX 301-2 based on all analyzed Insight-HXMT observations.The results are derived with four pulse phases of 0.15-0.4,0.4-0.6,0.6-0.9,0.9-1.15.Elements can be referred in Figure 5. properties of the two CRSFs and their relation to the magnetized NS, we model the accretion structure to study the simulated phaseresolved spectra.Based on the aforementioned observational results, we generate and fit the simulated spectra to study the physical origin of the CRSFs.
Our findings provide strong constraints on the magnetic inclination angle in an accreting neutron star.Previously, by assuming a circular beam along the axis of the dipole magnetic field, the magnetic inclination angles of radio pulsars could be measured with polarization data (Tauris & Manchester 1998).Without regard to the accretion torque, the evolution of magnetic inclination would be dominated by pulsar radiation spin down and inner viscous dissipation.Recent theoretical results indicate that isolated NS naturally evolves into two classes: near-aligned and near-orthogonal rotators -with typical pulsars falling into the latter category (Lander & Jones 2018;Novoselov et al. 2020).This scenario was supported by the increasing separation between the main pulse and interpulse in Crab pulsar, implying the growth of the magnetic inclination angle (Lyne et al. 2013).
However, the magnetic inclination angle in accreting X-ray pulsars is still poorly understood both in observations and theories.Preliminary studies deduced the magnetic inclination distribution from pulse profiles by assuming black-body radiation from a hot spot on the polar cap region (Annala & Poutanen 2010;Leahy & Matsuoka 1990).For accreting NSs in binary systems, the torque injection through either transient or prolonged accretion disc further complicates the physical picture.Magnetic angle evolution study in accreting NSs showed that the magnetosphere-disc interaction could orthogonalize the magnetic pole, and the magnetic inclination is greater than 55 • with all reasonable parameter combination (Biryukov & Abolmasov 2021).By predicting a termination of inclination evolution after the disc-spin alignment stage, the calculation implies that the NS magnetic inclination evolution would be again dominated by inner (possibly bulk) viscosity once the disc is aligned with spin, which is likely to drive the NS into an orthogonal rotator (Lander & Jones 2018;Biryukov & Abolmasov 2021).Our observations confirm that a slowly rotating and magnetized accreting NS in a binary system could indeed take an orthogonal configuration.
The radiation spectral properties (including the phase-resolved CRSF) of the NS accretion system GX 301-2 can be described by a radiation dominated accretion column of a ∼ 500 m shock height on the magnetic poles with a surface magnetic field strength of  ≃ 8 × 10 12 G.The NS magnetic field generally can be derived by  =  cyc (1+) 11.57keV 10 12 G, where  cyc is cyclotron line energy,  ∼ 1.2 is the gravitational redshift near the NS surface (Staubert et al. 2019), which determines  ∼ (4 − 5) × 10 12 G for GX 301-2.This simple calculation should be an under-estimate because apart from the general relativistic redshift effect, the bulk motion in accretion column and real line production regions could significantly lower the resonant energy from the value produced by the corresponding local magnetic field (Nishimura 2022).Moreover, the width of CRSF is mainly determined by the bulk velocity gradient in our simulation and the depth of the absorption profile is influenced by beaming pattern, which can explain the wide absorption feature in GX 301-2 and other accreting X-ray pulsars (Staubert et al. 2019).On the other hand, a less collimated (e.g.fan-like pattern) beaming pattern along with a smaller gradient would result in a broader and shallower line profile, providing a possible explanation for the non-detection of CRSFs in many accreting X-ray pulsars.
These rotating NSs in accreting systems could be continuous gravitational wave (GW) sources (Bildsten 1998).It has already been shown that for an axially symmetric rotating fluid neutron star to emit GW, the rotation axis and the magnetic field axis should not be aligned to each other (Bonazzola & Gourgoulhon 1996).An orthogonal rotator, due to its large toroidal  field, is likely to have significant quadrupolar distortion, thus would act as the optimal configuration for GW emission (Cutler 2002).The characteristic frequency of GW emissions from GX 301-2 would be around 3 mHz, which has the possibility of detection using future space-based GW telescopes like LISA (detailed calculations would be beyond the scope of this paper, would be resolved in next work).GW signals could provide insights into the magnetic inclination angle of GX 301-2 and other rotating neutron stars in accreting systems.

CONCLUSIONS
As a summary, phase-resolved spectral analysis is performed on the observation data of GX 301-2.Two CRSFs are detected in all phases, with their centroid energy ratio ∼ 1.6 − 1.7 in the super-critical case, while it scattered in the range of 1.2−2 when the source luminosity is lower.Based on the observed results, we simulate the phase-resolved X-ray spectra to determine the NS magnetic inclination.Our simulation shows that in order to obtain results similar to the observed data, the inclination angle needs to be greater than 70 degrees.This new method can be applied to other NS X-ray binaries and enable us to possibly derive the magnetic inclination angle distribution in different accreting NS systems.Thus it would provide the opportunity to probe the evolution of magnetic inclination which may influence the evolution of pulsar structure, including the relaxation process in the differentially rotating core and the complex interplay between magnetism and relativistic hydrodynamics, with the Equation of State (EoS) also playing a role (Lander & Jones 2017).Future X-ray polarization observations on these accreting X-ray pulsars are expected to provide new independent measurements of the magnetic inclination and the geometry of X-ray binary, e.g., recent IXPE polarization measurement suggested the magnetic inclination of Cen X-3 to be ∼ 17 • , corresponding to a near aligned rotator (Tsygankov et al. 2022).Therefore, the X-ray polarization combined with our methods will measure the magnetic inclination angle for more X-ray pulsars, probing the whole picture of magnetic structure evolution in accretion binaries, then should be interesting for NS physics.

Figure 1 .
Figure 1.Pulse profiles of P010130900103 in different energy bands.

Figure 2 .
Figure 2. Schematic of the line-forming region around an accretion mound.The region represented by oblique dot line denotes the line-forming region while the darker shadow region denotes the continuum-forming region.  is the radius of the accretion column.Solid arrows denote the velocity of bulk motion with the velocity  = /,  is the light speed.Wavy arrows denote diffusive flux emerging from the continuum-forming region.At  = 2 × 10 37 erg s −1 , the bulk velocity in the line-forming region is decelerated significantly by radiation pressure due to resonant scattering.

Figure 3 .
Figure3.The simulated spectra and the centroid energy ratios of different photon emission angles.(a) Centroid energy ratio between the two cyclotron absorption lines (y-axis), fitted from each simulated spectrum in a given photon emission angle (, x-axis), which is the angle between the photon's motion direction and the magnetic field when the photon interacts with the line forming region.E cyc is the centroid energy of the higher energy line in two observed cyclotron absorption lines.The dashed line is the best fitted ratio from the observations, which will be discussed later.(b) Simulated spectra produced in different emission angles.The absorption core shifts with the local emission angle in ∼ 40 − 60 keV due to Doppler boosting and anisotropic beaming pattern.

Figure 4 .
Figure 4. Main axes and angles of an accreting NS.The plane determined by the spin axis and observer is indicated with the black dashed line.The angles shown are Ψ (apparent photon emission angle), Θ AC (magnetic inclination),  (photon interacting angle with the line forming regions of accretion column),  (NS inclination) and Φ − Φ 0 (angle between the planes characterising the initial and instantaneous spin phase).

Figure 5 .
Figure5.Centroid energy ratios in phase-resolved spectra between the two CRSFs of GX 301-2 based on all analyzed Insight-HXMT observations.The ratios are derived from phase-resolved spectra in each short pointing exposures, with typical net exposure time of ∼ 3ks.Only ratios derived from more than four successfully fitted phase-resolved spectra in each observation are shown.The ratios converge to the value of ∼ 1.65 implied in the results of phase averaged spectroscopy (also see Figure13inDing et al. 2021 for details).The best-fit ratio, 1.69±0.03, is derived from the super-critical state with higher luminosity (≳ 1.8 × 10 37 erg s −1 in 5 − 70 keV, i.e. the vertical line).

Figure 6 .
Figure 6.Phase-resolved spectra derived in a super-critical observation.The ObsID is P010130900103 with  = 4.2 × 10 37 erg/s.Only results for the first three phases are presented, while the complete figure can be found in the appendix.Left panels: Fitting spectra with pulse profiles inserted, highlighting the corresponding pulse phases in each panel.The model used here is Constant*TBabs*(TBpcf*gabs*gabs*Continuum+Gaussian+Gaussian).Right panels: Residuals corresponding to the left panels with the same colors.Residuals from model with only one gabs are shown in each panel with purple lines, and their values are all increased by 6, while residuals from model with no absorption are shown in each panel with grey lines, and their values are all increased by 12.The reduced chi-square value is shown on the right side for each model.The broad band spectral shape is stable throughout all pulse phases.Almost all complications are due to the cyclotron (hard band) and neutral stellar wind (soft band) absorption.

Figure 7 .
Figure 7. Phase-resolved spectra derived in a lower luminosity observation.The ObsID is P010130900401 with  = 6.6 × 10 36 erg/s.Elements are the same as Figure 6.Only results for the first three phases are presented, while the complete figure can be found in the appendix.

Figure 8 .
Figure 8.The evolution of the best-fit parameters with pulse phases.The parameter explanations and their units are the same as described in Table2.ObsID P010130900103 results are shown on the left panels, which is a super-critical case, and obsID P010130900401 results with lower luminosity are shown on the right panels.Pulse profiles of the source in the HE band (30-60 keV) are shown in both sub panels with light blue lines.

Figure 10 .
Figure10.MC simulated spectra convolved with Insight-HXMT response matrix and background viewing from different angles.We assume the exposure is 150 s, and the gravitational light bending effect is not considered.We generated spectra in 10 − 70 keV range and fitted them with the same observational model without neutral absorption from interstellar medium and stellar wind, where  in the sub panels is equal to cos .

Figure 11 .
Figure11.Fitting results to simulated spectra after combining with Insight-HXMT ancillary and response matrix.(a) A heat map shows the photon interacting angle (cos ) with the line forming regions.We could observe these photons with a given magnetic axis (y-axis) in different spin phases (x-axis).The simulated spectra in Figure3(b) are integrated into each phase bin in this map.After assuming a specific magnetic inclination, we are able to use the simulated spectra to re-construct phase-averaged and phase-resolved spectra.(b) The generated input pulse profile when magnetic inclination is 85 • .(c) Fitting results (ODR) of simulated phase-resolved spectra with different magnetic inclinations, where   is the difference between the ratio of two cyclotron lines from the simulated spectra and the observed ratio obtained in Figure5.(d) Similar with (c) except the ratio from the simulated spectra is calculated using the average instead of ODR.All errors indicated in this figure are in 68% confidence level (1).Spectral simulations imply that the NS is a nearly orthogonal rotator with a magnetic inclination Θ AC ≳ 70 • .

Figure 13 .
Figure 13.Cyclotron line centroid energies of the pulse-phase resolved spectral analysis derived from two combined spectra.The combined observations and fitted cyclotron line centroid energy ratio ( 2 / 1 ) are shown on each panel.

Figure A1 .
Figure A1.Phase-resolved spectra derived in a super-critical observation.This is a completed version of Figure 6 with all 8 phases.The ObsID is P010130900103 with  = 4.2 × 10 37 erg/s.Left panels: Fitting spectra with pulse profiles inserted, highlighting the corresponding pulse phases in each panel.The model used here is Constant*TBabs*(TBpcf*gabs*gabs*Continuum+Gaussian+Gaussian).Right panels: Residuals corresponding to the left panels with the same colors.Residuals from model with only one gabs are shown in each panel with purple lines, and their values are all increased by 6, while residuals from model with no absorption are shown in each panel with grey lines, and their values are all increased by 12.The reduced chi-square value is shown on the right side for each model.The broad band spectral shape is stable throughout all pulse phases.Almost all complications are due to the cyclotron (hard band) and neutral stellar wind (soft band) absorption.

Figure A2 .
Figure A2.Phase-resolved spectra derived in a lower luminosity observation.This is a completed version of Figure7with all 8 phases.The ObsID is P010130900401 with  = 6.6 × 10 36 erg/s.Elements are the same as FigureA1.

Table 1 .
List of observations analyzed in this paper.The observation date, start MJD, HE exposure, count rates of three detectors, pulse period, orbital phase, and luminosity are also listed.