Cyclotron line evolution revealed with pulse-to-pulse analysis in the 2020 outburst of 1A 0535+262

We present a detailed analysis of the X-ray luminosity (Lx) dependence of the cyclotron absorption line energy (Ecyc) for the X-ray binary pulsar 1A 0535+262 during its 2020 giant outburst based on pulse-to-pulse analysis. By applying this technique to high cadence observations of Insight-HXMT, we reveal the most comprehensive Ecyc-Lx correlation across a broad luminosity range of ~(0.03-1.3)*10^38 erg/s. Apart from the positive and negative correlations between cyclotron line energy and luminosity at Lx~(1-3)*10^37 erg/s and ~(7-13)*10^37 erg/s, which are expected from the typical subcritical and supercritical accretion regimes, respectively, a plateau in the correlation is also detected at ~(3-7)*10^37 erg/s^-1. Moreover, at the lowest luminosity level (Lx<10^37 erg/s), the positive Ecyc-Lx correlation seems to be broken, and the pulse profile also occurs a significant transition. These discoveries provide the first complete view on the correlation between luminosity and the centriod energy of the cyclotron line, and therefore are relevant for understanding how accretion onto magnetized neutron stars depends on luminosity.


INTRODUCTION
X-ray pulsars (XRPs) are strongly magnetized rotating neutron stars (NSs) in binary systems, accreting matter from their massive (spectral type late O or early B) companion star (see Walter et al. 2015;Mushtukov & Tsygankov 2023;Neumann et al. 2023, for reviews).Because of the strong magnetic field ( ≳ 10 12 G), magnetic pressure dominates the ram pressure of the accreting gas within the relatively large Alfvén radius of ∼ 10 7−10 cm.In this case, the fully-ionized accreting plasma is funnelled by the magnetic field lines onto the magnetic poles of the NS, in which the kinetic energy of the infalling matter is converted to heat and radiation.A strong magnetic field in the vicinity of the NS surface could quantize the energy of electrons according to the Landau levels (see e.g.Isenberg et al. 1998; ★ E-mail: shuiqc@ihep.ac.cn (IHEP) † E-mail: szhang@ihep.ac.cn (IHEP) ‡ E-mail: wangpj@ihep.ac.cn (IHEP) Araya-Góchez & Harding 2000).Resonant scattering of photons in these electrons then leads to the generation of resonance features (generally appearing in absorption) in X-ray spectra of XRPs.These absorption-like features are called cyclotron resonant scattering features (CRSFs), and also called cyclotron lines (see Staubert et al. 2019, for a recent review).The centroid energy of the fundamental cyclotron line,  cyc , is directly related to the magnetic field strength of the line-forming regions as  cyc ≈ 11.6 ×  12 keV, where  12 is the magnetic field strength in units of 10 12 G, hence permits a direct measurement of the magnetic field strength near the NS surface.
Luminosity ( X ) dependence of the X-ray spectrum has been found in a number of XRPs (see e.g.Mihara et al. 1998;Reig & Nespoli 2013).In particular, two different types of luminosity dependence of the cyclotron line energy have been detected (Staubert et al. 2019).So far, a positive correlation between  cyc and  X has been observed in several sources, e.g.Her X-1 (Staubert et al. 2007), Vela X-1 (Fürst et al. 2014;La Parola et al. 2016), GX 304-1 (Yamamoto et al. 2011;Malacaria et al. 2015;Rothschild et al. 2017) and Cep X-4 (Vybornov et al. 2017), while a clear  cyc - X anti-correlation has been only found in V 0332+53 (Makishima et al. 1990;Tsygankov et al. 2006Tsygankov et al. , 2010;;Cusumano et al. 2016;Doroshenko et al. 2017;Vybornov et al. 2018) and 1A 0535+262 (Kong et al. 2021).In these two sources, a transition for correlation between positive and negative has been observed as well.In general, the negative  cyc - X correlation is detected in relatively high luminosity sources (e.g.V 0332+53), while the positive one is seen in lower luminosity sources (e.g.Her X-1).Accordingly, the observed bimodal behavior should be related to the mass accretion rate, which determines the structure of the accretion flow near the NS surface.At high mass accretion rates, the radiative pressure is strong enough to stop the falling matter near the polar cap, and hence supports an extended emission region above the NS surface, known as the accretion column (Basko & Sunyaev 1976).In this case, the accretion column grows higher with the increasing accretion rate.At low mass accretion rates, the gas deceleration is presumably dominated by Coulomb interactions (Burnard et al. 1991;Staubert et al. 2007;Becker et al. 2012), hence the characteristic emission height decreases with the increasing accretion rate.The transition between the aforementioned two different accretion regimes occurs at the so-called critical luminosity,  crit , which is associated with the onset of the accretion column (Basko & Sunyaev 1976;Becker et al. 2012;Mushtukov et al. 2015a).Since the cyclotron line energy is expected to reflect the characteristic emission height (Becker et al. 2012) and/or the bulk velocity of in-falling plasma (Mushtukov et al. 2015c), the  cyc - X relation provides an almost unambiguous technique to trace the different accretion regimes at work at the magnetic pole of an accreting neutron star (see e.g.Doroshenko et al. 2017;Vybornov et al. 2018;Kong et al. 2021).
In addition to the spectral properties, the angular pattern of the emitted radiation also depends on the mass accretion rate, resulting in the luminosity dependence of the pulse profile (Basko & Sunyaev 1976).At high luminosity ( X ≳  crit ), the emitted radiation primarily escapes through the walls of the accretion column, forming a so-called "fan beam" (Davidson 1973).At very low luminosity ( X ≪  crit ), the radiation is produced by hot spots or mounds at the polar cap and escapes roughly along magnetic field lines (Burnard et al. 1991;Nelson et al. 1993), resulting a "pencil beam".In the intermediate luminosity range ( X ≲  crit ), however, the emission pattern may be a hybrid combination of these two types (Blum & Kraus 2000).
In this study, we investigate the 2020 giant outburst of the transient accreting pulsar 1A 0535+262 discovered by the Ariel V satellite in 1975 with a pulsation period of ∼104 s (Rosenberg et al. 1975).This system contains an O9.7IIIe donor star and a strongly magnetized NS (Steele et al. 1998;Neumann et al. 2023).The eccentricity and period of the binary orbit are ∼0.47 and ∼110.3 days, respectively (Finger et al. 1996).The distance to this source is ∼2 kpc, which has been recently confirmed by Gaia (Bailer-Jones et al. 2018).Absorption line-like features at ∼50 and ∼100 keV were firstly detected with MirTTM/HEXE (Kendziorra et al. 1994), and have been confirmed by different missions and studies (see e.g.Kretschmar et al. 2005;Wilson & Finger 2005;Inoue et al. 2005;Ballhausen et al. 2017).Several studies have investigated the  cyc - X relation with intriguing results.Caballero et al. (2007) reported that there were no significant changes of  cyc with luminosity.On the contrary a possible positive  cyc - X correlation in a higher luminosity range was reported by Sartore et al. (2015).Based on the pulse-to-pulse technique, Klochkov et al. (2011) and Müller et al. (2013) revealed a positive  cyc - X correlation on the short time scale of the pulse period.Recently, Kong et al. (2021) revealed both positive and neg-ative  cyc - X correlations during the 2020 giant outburst based on the pulse-averaged spectroscopy derived from high-cadence observations of Insight-HXMT.Here, we focus on the luminosity dependence of cyclotron line energy during the 2020 giant outburst using pulseto-pulse analysis in order to study luminosity changes on timescales of the pulse period and therefore extending/refining the range of luminosity.We stress that this is significantly different approach than that used in Kong et al. (2021), and in fact new results are found.We introduce the Insight-HXMT observations and data reduction in Section 2, present the data analysis and results in Section 3, discuss these in Section 4 and finally summarize in Section 5.

OBSERVATIONS AND DATA REDUCTION
In the present study, we analyze the data from the Hard X-ray Modulation Telescope (Insight-HXMT, Zhang et al. 2014Zhang et al. , 2020)), the first Chinese X-ray satellite launched on 2017 June 15 with a science payload that allows observations in a broad energy band (1-250 keV).Insight-HXMT include three collimated telescopes: the High Energy X-ray telescope (HE, NaI/CsI, 20-250 keV), the Medium Energy Xray telescope (ME, Si pin detector, 5-30 keV), and the Low Energy X-ray telescope (LE, SCD detector, 0.7-13 keV), working in scanning and pointing observational modes and Gamma Ray Burst (GRB) mode.For details about the Insight-HXMT mission, see Zhang et al. (2019); Liu et al. (2020); Cao et al. (2020); Chen et al. (2020).
Fig. 1 presents the light curve of 1A 0535+262 during its 2020 giant outburst from the Swift/BAT transient monitor (Krimm et al. 2013).The shaded areas represent the time intervals of the Insight-HXMT observations which well cover almost the entire outburst period.The Insight-HXMT Data Analysis Software v2.05, together with the current calibration model v2.061 and the standard Insight-HXMT Data Reduction Guide v2.052 , are adopted in the extraction of the event data, under a series of criteria as recommended by the Insight-HXMT team: (1) the elevation angle (ELV) is larger than 10 • ; (2) the geometric cutoff rigidity (COR) is larger than 8 GeV; (3) the offset for the pointing position is less than 0.04 • ; (4) and the good time internals (GTIs) are at least 300 s away from the South Atlantic Anomaly (SAA).In addition, we exclude the data of detector boxes affected by the Crab (see also Kong et al. 2021Kong et al. , 2022;;Wang et al. 2022), since the Crab is located close to 1A 0535+262.The backgrounds are produced from blind detectors, using the tools: LEBKGMAP, MEBKGMAP and HEBKGMAP, version 2.0.9 based on the standard Insight-HXMT background models (Liao et al. 2020a,b;Guo et al. 2020).The XSPEC v12.12.0 software package (Arnaud 1996) is used to perform the spectral analysis.

Pulse-to-pulse Analysis
We produce pulse-amplitude-resolved spectra by applying the pulseto-pulse technique elaborated in Klochkov et al. (2011) and Müller et al. (2013).This approach has recently employed to investigate the luminosity dependence of the CRSF energy in Cepheus X-4 and V 0332+53 using NuSTAR data (see Vybornov et al. 2017Vybornov et al. , 2018)).The technique enables the analysis of the spectral variability on the time scale of a pulse period.Furthermore, by segmenting a single time-averaged spectrum into multiple spectra at different luminosity levels, this method expands the luminosity range and increases the number of data points, providing detailed information regarding the luminosity dependence of the CRSF energy.
In our pulse-to-pulse analysis, we extract the net light curve from LE (1-10 keV), ME (10-30 keV) and HE (30-120 keV), respectively.The arrival times of photons are corrected to the solar system barycentre with the HXMTDAS tool hxbary, and the effects of binary orbital modulation are corrected using the ephemerides provided by Camero-Arranz et al. (2012).Based on the binary-and barycentrecorrected data, we measure the pulse periods and produce pulse profiles (see Wang et al. 2022;Hou et al. 2023, for details about the timing analysis).Next, we measure the number of counts from ME during each pulse period, as its energy range occupies an intermediate position and typically encompasses the peak flux in the spectrum of an X-ray Pulsar (XRP).The number of counts summed over a pulse represents the pulse amplitude and is thus used as a measure of the pulse brightness.According to the measurement of the pulse amplitude in each pulse, we can build a frequency distribution of pulses as a function of the pulse amplitude (Fig. 2).Based on the distribution of the pulse amplitude, the pulses can be divided into two bins at different brightness levels.As shown in the bottom panel of Fig. 2, the pulse profiles in different amplitude bins indeed have significantly different count rates.The selection of the number of the bins mainly depends on the available counting statistics of the energy spectra, which should ensure that the cyclotron line energy can be well constrained in the spectral fitting.We stack the GTIs of pulses belonging to the same pulse-amplitude bin and subsequently extract the energy spectra from LE, ME, and HE, respectively, using the stacked GTIs.We defined a single observation as the combination of exposures within an one-day period.However, low-luminosity spectra exhibit poor counting statistics.To bolster the robustness of spectral fittings during low-luminosity intervals, we moderately relaxed the criteria for combining exposures during these stages.Specifically, for the outburst rising phase, exposures from MJD 59,163 to 59,165 are amalgamated into a single observation, while for the fading phase, exposures from MJD 59,204 to 59,207 are likewise merged into one observation.Following this, the pulse-to-pulse analysis is executed for each individual observation.We note that moderate changes in the combination strategy have small impact on the spectral outcomes.We have tried to build histograms of the pulse amplitude with a longer-time period, and the obtained spectral results are similar to the those from the original 1-day combination.However, constructing the amplitude histograms solely based on luminosity may introduce bias, since a time drift of the cyclotron line energy has been reported (see Kong et al. 2021).This suggests that there could be intrinsic differences in the spectra between the rising and fading phases.
In the spectral analysis, the energy bands adopted for the spectral fittings are 2-7 keV for LE, 8-30 keV for ME and 30-120 keV for HE.Considering the current accuracy of the instrument calibration, we add 0.5%, 0.5%, and 1% systematic errors to the energy spectra for LE, ME, and HE, respectively.Following the recommendation of the Insight-HXMT calibration group, the spectra are rebinned as follows: (1) LE: channels 0-579 and 580-1535 are respectively grouped with 5 and 10 bins in each group; (2) ME: channels 0-1023 are grouped with 2 bins in each group; (3) HE: channels 0-255 are grouped with 2 bins in each group.With this grouping strategy, the spectra have at least 100 counts per energy bin, resulting in spectral data that closely resemble the Gaussian distribution.Additionally, the background spectra of the three in-struments of Insight-HXMT are dependent on background models, which are mainly based on a large amount of the blank sky observations, so the spectral data of the background are not Poissonian but Gaussian.As a result,  2 statistics is used in the spectral analysis.In the present study, we use the same spectral model: con-stant×tbabs× mgabs×(bbodyrad+bbodyrad+Gaussian+cutoffpl), as that used in the time-averaged spectral analysis (Kong et al. 2021), to fit each energy spectrum.In the fitting model, inter-calibration constants are required among the instruments.The constant value of the LE is fixed at 1, while the constant values of the ME and HE are allowed to vary.mgabs is a multiplicative absorption model with a Gaussion profile, where  ′ () is the spectrum modified by mgabs,  is the photon energy,  cyc is the cyclotron line centroid energy,  1 and  1 are the optical depth and the width of the line, respectively,  is the centroid energy ratio of the first harmonic line and fundamental line, and  2 and  2 are the optical depth and the width of the first harmonic line, respectively.To study the properties of the fundamental CRSF line energy, we fix the hydrogen column density ( H ), the ratio of the first harmonic line and fundamental line of CRSFs (), the width of the fundamental line of CRSFs ( 1 ), the width of the first harmonic line of CRSFs ( 2 ) and the width of the iron emission line ( Fe ) to 0.59 × 10 22 cm −2 , 2.3, 10 keV, 5 keV and 0.3 keV, respectively (see also Kong et al. 2021).We find each spectral fitting with the model obtains a satisfactorily reduced  2 .To compute the uncertainties at the 68% confidence level,we employ the Markov Chain Monte Carlo (MCMC) method using the Goodman-Weare algorithm with 8 walkers and a total length of 40,000 (Goodman & Weare 2010), and the initial 2000 elements are discarded as a burn-in period.We find that the autocorrelation length is typically 1000 elements, so the net number of independent samples in the parameter space we have is of the order of 10 4 .In order to further test the convergence of the MCMC chain, we compare the one-and two-dimensional projections of the posterior distributions for fitting parameters from the first and second halves of the chain, and we find no significant differences (see Fig. B1).In Fig. 3, we display the representative pulse-amplitude-resolved spectra and present the corresponding best-fit spectral parameters of each model in Table 1.
The best-fitting parameters of all spectra in this study are provided in Table A1.

Luminosity Dependence of the CRSF Centroid Energy
Based on spectral fittings of the pulse-amplitude-resolved spectra, we calculate the X-ray luminosity ( X ) in the 2-150 keV energy range, assuming the source is at a distance of 2 kpc.To achieve this, we extend the energy response of the model from 0.1 to 200 keV using the  command in XSPEC.Subsequently, the model fluxes in the 2-150 keV energy range and their corresponding uncertainties are computed for the three instruments using the   command in XSPEC when the MCMC chains are loaded.This allows for the derivation of the flux distribution from the MCMC chains.The luminosities reported here represent average values derived from the three instruments, and the associated errors are calculated using standard error propagation formula.It is worth noting that the reported luminosity is not corrected for absorption.However, this does not affect the luminosity dependence of  cyc , since the value of  H is fixed in the fitting.Additionally, we find that the systematic difference between the values of unabsorbed flux and that linked to absorption is less than 1%, primarily due to the relatively small distance of the source.The left panel of Fig. 4 displays the luminosity dependence of the CRSF centroid energy, with orange squares and blue circles representing the rising and fading phases of the outburst, respectively.As a result of the pulse-to-pulse analysis, the number of data points shown in the left panel of Fig. 4 is twice as large as the result from the pulse-averaged spectral analysis (see Kong et al. 2021), enabling us to study the  cyc - X relations in greater detail.Our analysis confirms the previously observed anti-correlation between  cyc and  X at high luminosity when ( X ≳ 7 × 10 37 erg s −1 ).However, at lower luminosity, we find a more complex relationship between  cyc and  X .Specifically, in the luminosity range of ∼ (3-7) × 10 37 erg s −1 ,  cyc is not correlated with  X .In the narrow luminosity range of ∼ (1-3) × 10 37 erg s −1 ,  cyc is positively correlated with  X during both the rising and fading phases.At the lowest luminosity level ( X ≲ 1 × 10 37 erg s −1 ), the positive  cyc - X correlation appears to break, and  cyc shows an increasing trend with decreasing  X during both the rising and fading phases.Such a complex behaviour is reported here for the first time.
It is also important to note that the measured  cyc during the rising phase is systematically lower than that during the fading phase, a behavior that has been reported in the pulse-averaged spectral analysis by Kong et al. (2021).Similarly, line energy differences between the two outburst stages were detected in V 0332+53 during its 2015 giant outburst (Cusumano et al. 2016;Doroshenko et al. 2017).However, in the case of 1A 0535+262, contrary to what reported in literature,  cyc during the rising phase is smaller than during the fading phase.Nevertheless, the trends in evolution between the two phases are nearly parallel.
To study changes in the  cyc - X relation by combining data points from the rising and fading phases, we follow the method of Doroshenko et al. (2017), Vybornov et al. (2018), and Kong et al. (2021), assuming that  cyc undergoes a linear time drift (  cyc ) during the outburst.Specifically, we assume that each observation probes an "instantanous" changes in the cyclotron line energy.However, the influence of the time variation is negligible within individual observations.Then, for each data point present in the left panel of Fig. 4,  cyc is adjusted by subtracting a factor of  cyc ×Δ, where Δ represents the difference between the reference time (MJD 59,159) and the time of observation.We then use a broken power-law model with three break luminosities to fit the data points with the linear energy drift of  cyc taken into account.This model has nine free parameters: four powerlaw indices (Γ 1 , Γ 2 , Γ 3 , and Γ 4 ), three break luminosities ( 1 ,  2 , and  3 ), normalization (), and the linear time drift of line energy,  cyc .The model fitting is performed by using the MCMC method with the Goodman-Weare algorithm with 8 walkers and a total length of 40,000.The uncertainty in  cyc at each data point is accounted for using a heteroscedastic Gaussian likelihood, and the uncertainty in  X is accounted for by incorporating a sampling process during the execution of the MCMC chain.The luminosity dependence of  cyc , corrected for  cyc , is presented in the right panel of Fig. 4, with the best-fitting result using the broken power-law model plotted as the red dashed line.The model yields luminosity break values of  1 = (1.14±0.12)×10 37 erg s −1 ,  2 = (2.89±0.35)×10 37 erg s −1 , and  3 = (7.47 ± 0.26) × 10 37 erg s −1 .These transitional luminosity levels suggest the division of the  cyc - X relation into four zones (Zones I, II, III, and IV) with different types of the  cyc - X correlation.For the linear time drift,  cyc , we obtain an increase rate of 0.046±0.004keV day −1 for cyclotron line energy, consistent with the result of Kong et al. (2021).We also find that the best-fitting statistics of the broken power-law model with only one break at  3 improve from  2 = 219.62 for 73 degrees of freedom (d.o.f.) to  2 = 117.7 for 69 d.o.f. for the presented three-breaks broken power-law model.This indicates that the data preferred the model with three breaks at ∼ 5.8 confidence based on the F-test.Performing a similar comparison between broken power-law models with two breaks at  1 and  3 and three breaks, the F-test gives ∼ 4.9 confidence for the three breaks.This implies that the data prefer that the  cyc - X relation changes at  2 at ∼ 4.9 confidence.It is also important to note that the estimated confidence level can be affected by the choice of fitting model and the assumption that the CRSF centroid energy exhibits a linear drift during the outburst.

Relationship Between the CRSF energy and Pulse Profiles
The Insight-HXMT data have been utilized for providing detailed timing analyses of 1A 0535+262 during its 2020 giant outburst, as reported by Reig et al. (2022); Wang et al. (2022); Ma et al. (2022).Theoretical models suggest that the beam pattern changes across different accretion regimes (e.g.Basko & Sunyaev 1976;Becker et al. 2012).Therefore, analyzing the correlation between  cyc and  X , along with the pulse profile, can offer valuable insights into state transitions.In general, high-energy photons mainly originate from the polar caps through direct emission, while low-energy photons include contributions from direct emission, thermal components, and reflection on the NS surface (see e.g.Hu et al. 2023;Li et al. 2023).
Consequently, the pulse profile in the high-energy range is more effective in tracing beam pattern transitions.Here, we average the pulse profiles in the high-energy bands of 30-120 keV within each pulse bin according to the pulse-to-pulse analysis performed in Section 2, resulting in the amplitude-resolved pulse profiles (referred to as pulse profiles hereafter).The left panel of Fig. 5 displays the luminosity dependence of the pulse profile, with the color bar representing the normalized intensity.Based on the evolution of the pulse-averaged pulse profiles, Wang et al. (2022) suggested a transition luminosity of nearly 1.1 × 10 37 erg s −1 .We note that our results from the pulseto-pulse analysis are consistent with those of Wang et al. (2022), as we find a transition luminosity close to this value (see Figure 5).Interestingly, such a transition luminosity aligns with  1 , which represents the lowest break luminosity in the  cyc - X correlation.By examining the luminosity dependence of the pulse profile and that of  cyc together, we find a co-evolution between these two aspects.At luminosities below ∼  1 = (1.14 ± 0.12) × 10 37 erg s −1 ,  cyc exhibits to be anti-correlated with  X , and the pulse profile in the high energy bands (30-120 keV) appears as a single broad peak.In the luminosity range of ∼ (1 − 3) × 10 37 erg s −1 ( 1 ≲  X ≲  2 ), where the  cyc - X correlation is positive, the pulse profile changes from a single-peak structure to a double-peak one.When  2 ≲  X ≲  3 , the pulse profile exhibits a significant double-peak structure, and  cyc remains relatively constant.In Zone IV, the pulse profile still shows a double-peak structure, but the two peaks in one cycle have significantly different intensities, and  cyc shows an anti-correlation  .Left: Luminosity ( X ) dependence of the cyclotron line centroid energy ( cyc ) derived from the pulse-to-pulse analysis, where the uncertainties are reported at the 68% confident level (1).The data points from the rising and fading phases are plotted as orange squares and blue circles, respectively.Right: Luminosity dependence of cyclotron line centroid energy with the linear energy drift of  cyc = 0.047 keV day −1 taken into account (symbols are the same as in the left panel).The red dashed line represents the best fitting result with the broken power-law model which has three break luminosities ( 1 ,  2 and  3 ).These three vertical dashed lines and shaded ares are the best fitting break luminosities and corresponding uncertainties, respectively.These transitional luminosities suggest the division of the  cyc - X relation into four zones (Zones I, II, III, and IV) with different types of the  cyc - X correlation.
with luminosity.To clearly illustrate the evolution of the pulse profile across different zones, we present example pulse-amplitude-averaged profiles for the four zones in the right panel of Fig. 5.

DISCUSSION
We have studied the luminosity dependence of the cyclotron line energy of 1A 0535+262 during its 2020 giant outburst.Using highcadence observations from Insight-HXMT, we performed the pulseto-pulse analysis to reveal the comprehensive scenario of the  cyc - X relation in a broad luminosity range of ∼ (0.03-1.3) × 10 38 erg s −1 .Specifically, through the use of a large number of data points obtained from the pulse-to-pulse analysis, we were able to fit the  cyc - X correlation with a broken power-law model.Based on the fitting, we identified three transitional luminosity levels ( 1 ,  2 , and  3 ) with a confidence level of ∼ 5.These transitional luminosity levels suggest the division of the  cyc - X relation into four zones (Zones I, II, III, and IV) with different types of the  cyc - X correlation.Our results confirm the recently reported anti-correlation between  cyc and  X in the high luminosity range (Zone IV, i.e.  X ≳  3 ) from pulse-averaged spectral analysis (see Kong et al. 2021).However, we found that the  cyc - X relation at the lower luminosities ( X ≲  3 ) is more complex, and the pulse profile also shows significant changes.
In Zone III ( 2 ≲  X ≲  3 ),  cyc remains roughly constant, and the pulse profile in the high energy bands (30-120 keV) displays a double-peak structure.In Zone II ( 1 ≲  X ≲  2 ),  cyc shows positive correlation with  X , and the pulse profile exhibits a transition between double-peak and single-peak structure.At the lowest luminosity level (Zone I, i.e.  X ≲  1 ),  cyc seems to be anti-correlated with  X , and the pulse profile becomes a pure single broad peak.
The variation of  cyc with pulse phase is generally believed to be due to the changing viewing angle under which the emission regions are seen (Staubert et al. 2007).A phase-resolved analysis of this source with Insight-HXMT data have been conducted by Kong et al. (2022).They showed that the transition of  cyc - X correlation at the high luminosity level (∼ 6 × 10 37 erg s −1 ) identified in our analysis was also valid when certain pulse phases were analyzed.However, the remarkable phase dependency of  cyc variations above ∼ 6 × 10 37 erg s −1 were observed (see Fig. 3 of Kong et al. 2022), which means the phase variations of  cyc might potentially affect the obtained  cyc - X correlation based on the phase-averaged spectral fitting in Zone IV.Nevertheless, it appears that the influence of this phase dependency is relatively weak at luminosities in Zones I, II and III.To provide a more specific discussion, in this study, we disregard this effect for a pulse-to-pulse analysis with higher counting statistics and focus on studying the pulse-amplitude-resolved spectra averaged over the pulse profile.

Cyclotron line behaviour at 𝐿
Since the work of Basko & Sunyaev (1976), it is generally accepted that the accretion mode of an accreting pulsars changes between suband supercritical regimes with the mass accretion rate, , hence with luminosity.The transition between the two accretion regimes occurs at the critical luminosity,  crit , which is determined by the properties of the neutron star (e.g., magnetic field strength, Basko & Sunyaev 1976;Becker et al. 2012;Mushtukov et al. 2015a,b).In the supercritical regime, where  X ≳  crit , the radiation pressure dominates the gas deceleration, leading to the formation of an accretion column whose height could be positively correlated with the mass accretion rate.In this case, the energy of CRSF forming in accretion column (Tsygankov et al. 2007) or at the stellar surface due to the reflection and reprocessing of X-ray photons (Poutanen et al. 2013;Lutovinov et al. 2015) is expected to be anti-correlated to luminosity.
At X ≲  crit , the accreting gas is decelerated in the atmosphere of a neutron star (Zel'dovich & Shakura 1969) or at some height either via Coulomb interactions (Staubert et al. 2007;Becker et al. 2012), or through a collissionless shock above the stellar surface (Langer & Rappaport 1982;Bykov & Krasil'Shchikov 2004).These scenarios require the positive correlation between the cyclotron line energy and accretion luminosity.In the case of accretion flow deceleration in the atmosphere of a neutron star, the positive correlation is due to the Doppler effect in the accretion channel above the polar caps (see e.g., Mushtukov et al. 2015c).In the case of Coulomb interactions or collisionless shock above the stellar surface, the positive correlation is due to the expected variations of emission region height above the surface (see Discussion in Rothschild et al. 2017).
The transition of X-ray pulsars between sub-and supercritical states is expected to be associated with the changes of a beam pattern forming close to the neutron star surface (Gnedin & Sunyaev 1973).Beam pattern variations change the typical angle between the field direction and the momentum of photons leaving the line-forming region.It can affect the line centroid energy (Nishimura 2014) and possibly complicates the transition between the positive and negative correlation at  X ∼  crit .
As presented in Figs 4 and 5, the cyclotron line energy is positively correlated to  X in Zone II, where  X ∼ (1-3) × 10 37 erg s −1 .
Conversely, the sign of the correlation turns negative when  X ≳  3 (Zone IV).These two luminosity ranges could potentially signify sub-critical and supercritical regimes, respectively.Before the 2020 outburst of 1A 0535+262, the occurrence of both positive and negative  cyc - X correlations in an outburst was a rare phenomenon, previously observed only in V 0332+53 (Makishima et al. 1990;Tsygankov et al. 2006;Cusumano et al. 2016;Doroshenko et al. 2017;Vybornov et al. 2018).It is worth noting that previous studies of 1A 0535+262 at luminosities below  3 yielded inconclusive results; some claimed a positive correlation (Sartore et al. 2015;Kong et al. 2021), while others reported the absence of any correlation (Tsygankov et al. 2006;Caballero et al. 2007Caballero et al. , 2013;;Müller et al. 2013).The complexity of these findings may be attributed to the relatively narrow luminosity range of the positive correlation in this source, making it challenging to detect.Thanks to the high-cadence observational strategy and a pulse-to-pulse technique, the positive  cyc - X correlation can now be robustly observed by Insight-HXMT.It is important to note that the observed luminosity range of the positive  cyc - X correlation in this study is consistent with that reported by Sartore et al. (2015) using INTEGRAL/SPI observations.Moreover, we observe that  cyc appears to exhibit a plateau in an intermediate luminosity range (Zone III) between Zones II and IV as determined by fitting with the broken power-law model (see Fig. 4).This plateau of the  cyc - X correlation observed in the luminosity range of  2 - 3 constitutes a novel observational phenomenon, potentially offering new perspectives for theoretical studies to explore the complex physical processes within the accretion column.A possible explanation for the plateau could be that when  1 ≲  X ≲  2 , the source enter a mixed regime.In this mixed regime, it is possible that both Coulomb interactions and radiative pressure play crucial roles in decelerating the infalling gas within the accretion column.Consequently, the radiative height of the CRSF might exhibit relative stability within this luminosity range.While in Zone II, Coulomb interactions dominate the deceleration, resulting in the positive correlation between  cyc and  X .In this scenario, the so-called critical luminosity might be located in the Zone III.

Cyclotron line behaviour at 𝐿
At the lowest luminosity level (Zone I,  X ≲ 1 × 10 37 erg s −1 ), the early researches reported less changes of the cyclotron line energy (e.g.Terada et al. 2006;Caballero et al. 2007Caballero et al. , 2013)).Additionally, a detailed analysis of a deep NuSTAR observation at a low luminosity of ∼ 7 × 10 34 erg s −1 by Tsygankov et al. (2019b) excluded a positive correlation at low luminosities of this source.In this study, we observe a break in the positive  cyc - X correlation when  X <  1 .Furthermore, the pulse profile displays a single broad peak in this luminosity range (see Fig. 5).It is important to note that a similar transition of the correlation during the 2011 outburst of 1A 0535+262 have been reported by Sartore et al. (2015) using INTEGRAL/SPI observations, who claimed that the transition occurred between 6.8 × 10 36 and 1.2 × 10 37 erg s −1 , likely corresponding to the Coulomb-stopping limit luminosity (Becker et al. 2012).They also detected changes in the pulse profile during the transition (see Fig. 12 of Sartore et al. 2015).Our correlation fitting suggests a transition luminosity of ∼ 1.14 × 10 37 erg s −1 , consistent with their result.In comparison with this previous work, we used more observations and a novel technique to robustly confirm such a transition at a ∼ 5 confidence level.
In the sub-critical regime, if the final deceleration of the in-falling gas occurs via Coulomb interactions, ones expect there could be a Coulomb-stopping limit luminosity (see Becker et al. 2012, for de-tails),  coul ≈ 1.17 × 10 37 × ( * /10 12 G) −1/3 erg s −1 , where  * is the magnetic field strength at the stellar surface.For 1A 0535+262, the detected  cyc is around 46 keV, so we can estimate  coul at ∼ 0.74 × 10 37 erg s −1 .We note that the estimated  coul is close to the lowest transition luminosity,  1 = (1.14 ± 0.12) × 10 37 erg s −1 .However, it is not very conclusive what happens when  X ≲  coul .A possible scenario is that the gas deceleration occurs via passage through a gas-mediated shock near the NS surface (Langer & Rappaport 1982), forming a small accretion mound (Mukherjee et al. 2013).In this case, the emission region approaches the stellar surface quickly as  X is reduced and the radiative beam pattern consists of pencil component only (Burnard et al. 1991;Nelson et al. 1993).When  X decreases, the change in the sign of the  cyc - X correlation from positive to negative has also been observed in GRO J1008-57 during its 2017 outburst (Chen et al. 2021).However, the pulse profile in GRO J1008-57 does not exhibit significant changes during the transition of the  cyc - X correlation (Ge et al. 2020).Consequently, Chen et al. (2021) proposed variations in the accretion channel geometry to explain the observed anti-correlation between  cyc and  X at lower luminosity levels.This is attributed to the fact that a reduction in the length of the accretion channel could potentially lead to a smaller redshift of the line energy by increasing the local X-ray energy flux, even as the total luminosity decreases.
Another explanation of the negative correlation at low luminosity level is related to the specific structure of NS atmosphere.If the accretion flow is stopped by the Coulomb collisions in the upper atmospheric layers, one would expect an inverse temperature profile with the overheated optically layers and relatively cold underlying atmosphere (Zel'dovich & Shakura 1969).At very low mass accretion rates, this atmospheric structure results in the appearance of two components of broadband energy spectra: the low-energy one related to the thermal emission of the atmosphere and the high-energy one due to emission and further reprocessing of cyclotron photons due to magnetic Compton scattering and absorption (see, e.g., Mushtukov et al. 2021;Sokolova-Lapa et al. 2021).This spectral transition was reported in a few accreting pulsars, including 1A 0535+262 (Tsygankov et al. 2019a,b).However, in this study, there is no evidence of an additional high-energy component, as reported by Tsygankov et al. (2019b) at low luminosities below 10 36 erg s −1 , appearing in spectra.We suggest this may because all Insight-HXMT observations on this source were conducted at luminosities beyond this level.According to the numerical simulations based on the Monte Carlo approach, the cyclotron absorption feature tend to be shifted towards higher energies with respect to the cyclotron resonance energy in the atmosphere (Mushtukov et al. 2021).Thus, the negative correlation between the cyclotron line energy and luminosity observed at low mass accretion rates can be related to the changes in the atmospheric structure.
An alternative explanation could be that when the accretion rate is extremely low ( X ≲ 10 −4  crit ), one would expect no gas-mediated shock above the NS and the optical thickness of the accretion channel becomes low even for cyclotron photons and the scattering feature is formed in the NS atmosphere (see Section 3.2 in Mushtukov et al. 2015c).In this situation, there would be no redshift of the line energy, and the line energy is expected to exhibit an anti-correlation with the luminosity (Chen et al. 2021).However, the transitional luminosity,  1 (∼ 10 37 erg s −1 ), observed in 1A 0535+262 is at least three orders of magnitude higher than what is anticipated in this scenario.

Cyclotron line energy drift
As already mentioned in Sec. 3,  cyc in the rising phase is systematically lower than that in the fading phase, although the behaviors of the  cyc - X relation from the two stages are similar (see also Kong et al. 2021).The difference of  cyc between the two outburst stages could be corrected by introducing a phenomenological increase rate of 0.046 keV d −1 to the cyclotron line energy.Differences in  cyc from the two stages of an outburst have been also found in the 2015 outburst of V 0332+53 (Cusumano et al. 2016;Doroshenko et al. 2017;Vybornov et al. 2018).However, in the case of V 0332+53,  cyc in the rising phase of the outburst is larger than in the fading phase.Cusumano et al. (2016) suggest that this decrease in  cyc might be due to a rapid dissipation of the external magnetic field of the NS, affected by diamagnetic screening from the accreted matter.Since we observe an increase of  cyc , the scenario proposed by Cusumano et al. (2016) cannot be applied in this case.An alternative interpretation invoked by Doroshenko et al. (2017) is that the energy drift of the CRSF was related to the different emission geometries because of different magnetosphere sizes in the rising and fading phases presumably associated with changes in the internal structure of the accretion disk.This scenario suggests that in comparison to V 0332+53, the geometrical evolution in 1A 0535+262 is opposite: the inner radius of the disc is larger and the polar cap footprint of the magnetic field is smaller during the rising phase compared to those of the fading phase.Consequently, there is a higher accretion column in the supercritical regime of the rising phase.Furthermore, if the observed  cyc in the sub-critical regime is affected by the Doppler effect (see Mushtukov et al. 2015b), one could expect lower  cyc values in the sub-critical regime of the rising phase.This is because, at a certain luminosity level and a specific height above the polar cup, the local radiation pressure upon the accreting matter is larger for a smaller radius of the hot spot, thereby enhancing the braking effect on the falling material.However, the reasons behind the different inner disk radius evolutions in the two sources remain unclear and are beyond the scope of our current study.

CONCLUSION
Based on the pulse-to-pulse analysis of Insight-HXMT high cadence observations of 1A 0535+262 during its giant outburst in 2020, we have detected three transitional luminosity levels in the luminosity dependence of the cyclotron line energy ( cyc ) for the first time in a Type II outburst.We find both positive and negative correlations between cyclotron line energy and luminosity at  X ∼ (1-3) × 10 37 erg s −1 and  X ∼ (7-13) × 10 37 erg s −1 , respectively, during both rising and fading phases of the outburst.These two luminosity ranges are generally expected to correspond to the typical sub-and super-critical accretion regimes, respectively.However, an unexpected plateau in the correlation between the CRSF energy and luminosity is observed in the range of ∼ (3-7) × 10 37 erg s −1 .We suggest that the plateau of line energy might result from a combination of Coulomb interactions and radiative pressure decelerating the accreting gas within the accretion column.Additionally, below the transitional luminosity of ∼ 10 37 erg s −1 , the positive  cyc - X correlation appears to break, and the pulse profile simultaneously transitions to a pure single broad peak.We propose there are several models might account for the  cyc behavior at  X ≲ 10 37 erg s −1 .These findings offer the first comprehensive view of the luminosity dependence of CRSF energy across a broad luminosity range of ∼ (0.03-1.3) × 10 38 erg s −1 , providing valuable insights for further theoretical work.
Table A1.Best-fitting parameters and goodness of fit of all pulse-amplitude-resolved spectra. cyc and  1 are the centroid energy and optical depth of the fundamental cyclotron line, respectively,  2 is the optical depth of the harmonic line, Γ and  fold are the spectral index and cutoff energy of the cutoffpl model,  1 and  2 are temperatures of the two bbodyrad models, respectively,  Fe is the centroid energy of the iron line, and  X is the luminosity in the 2-150 keV energy range.In our analysis of the low-pulse-amplitude spectrum observed on MJD 59159 we find that the lower-temperature blackbody component is not well constrained.Similarly, the harmonic cyclotron line and the iron emission line in the low-pulse-amplitude spectrum from MJD 59204-59207 observations are also poorly constrained.Therefore, we fixed these parameters to those of the nearest in time spectra where they are well constrained.The uncertainties are reported at a 68% (1) confidence level.  .One-and two-dimensional projections of the posterior probability distributions, and the 0.16, 0.5, and 0.84 quantile contours derived from the MCMC analysis for each free physical parameter.To test the convergence, we compare the one-and two-dimensional projections of the posterior distributions from the first (blue) and second (orange) halves of the chain, and we find no large differences.This illustration corresponds to the spectral fitting of the high-pulse-amplitude bin on MJD 59167.The best-fitting parameters of this example are summarized in Table 1.

Figure 1 .
Figure 1.Swift/BAT light curve of 1A 0535+262 during the 2020 giant outburst with time time resolution of one point per day.The shaded areas indicate the time intervals of Insight-HXMT observations.

Figure 2 .
Figure 2. Top: distributions of counts from ME in individual pulse (i.e. of pulse amplitudes) for the observation on MJD 59183, the red dashed line indicates the boundary of the amplitude bins, and colored shaded areas represent the counting ranges of low and high amplitude bins, respectively.Bottom: averaged ME pulse profiles of higher (blue) and lower (gray) pulseamplitude bins whose counting ranges are presented in the top panel.

Figure 3 .
Figure3.Unfolded pulse-amplitude energy spectra (the top panels) and the corresponding fitting residuals (the two bottom narrow panels) of MJD 59167 and MJD 59199.The spectral analysis is performed in the energy bands of 2-7 keV, 8-30 keV and 30-120 keV for LE, ME and HE, respectively.The middle panels show the best-fitting residuals with model: cons×tbabs× mgabs×(bbodyrad+bbodyrad+Gaussian+cutoffpl), while the bottom panels display the residuals without the CRSF line model (mgabs).
Figure 4. Left: Luminosity ( X ) dependence of the cyclotron line centroid energy ( cyc ) derived from the pulse-to-pulse analysis, where the uncertainties are reported at the 68% confident level (1).The data points from the rising and fading phases are plotted as orange squares and blue circles, respectively.Right: Luminosity dependence of cyclotron line centroid energy with the linear energy drift of  cyc = 0.047 keV day −1 taken into account (symbols are the same as in the left panel).The red dashed line represents the best fitting result with the broken power-law model which has three break luminosities ( 1 ,  2 and  3 ).These three vertical dashed lines and shaded ares are the best fitting break luminosities and corresponding uncertainties, respectively.These transitional luminosities suggest the division of the  cyc - X relation into four zones (Zones I, II, III, and IV) with different types of the  cyc - X correlation.

Figure 5 .
Figure 5. Left: Luminosity dependence of pulse profiles in the energy band of 30-120 keV during the whole outburst.The color bar displays the intensity of the pulse profile normalized in the 0-1 range.These three vertical dashed lines are the break luminosities ( 1 ,  2 and  3 ) shown in the left panel of Fig. 4, and the shaded areas are the corresponding uncertainties at 1 confidence level.These transitional luminosities suggest the division of the  cyc - X relation into four zones (Zones I, II, III, and IV).Right: Example pulse-amplitude-averaged profiles for the four distinct zones.These example pulse profiles for the four zones are from observations of MJD 59202, 59196, 59188 and 59173, respectively.
Figure B1.One-and two-dimensional projections of the posterior probability distributions, and the 0.16, 0.5, and 0.84 quantile contours derived from the MCMC analysis for each free physical parameter.To test the convergence, we compare the one-and two-dimensional projections of the posterior distributions from the first (blue) and second (orange) halves of the chain, and we find no large differences.This illustration corresponds to the spectral fitting of the high-pulse-amplitude bin on MJD 59167.The best-fitting parameters of this example are summarized in Table1.

Table 1 .
Best-fit Parameters of the representative Pulse-amplitude-resolved Spectra.The spectral analysis is performed in the energy bands of 2-7 keV, 8-30 keV and 30-120 keV for LE, ME and HE, respectively.The 0.5%, 0.5%, and 1% systematic errors are added to the energy spectra for LE, ME, and HE, respectively.Uncertainties are reported at the 68% confident level (1) using the MCMC technique, with length of 40,000.