Confirmation of the presence of a CRSF in the NICER spectrum of X 1822-371

X 1822-371 is an eclipsing binary system with a period close to 5.57 hr and an orbital period derivative $\dot{P}_{\rm orb}$ of 1.42(3)$\times 10^{-10}$ s s$^{-1}$. The extremely high value of its $\dot{P}_{\rm orb}$ is compatible with a super-Eddington mass transfer rate from the companion star and, consequently, an intrinsic luminosity at the Eddington limit. The source is also an X-ray pulsar, it shows a spin frequency of 1.69 Hz and is in a spin-up phase with a spin frequency derivative of $7.4 \times 10^{-12}$ Hz s$^{-1}$. Assuming a luminosity at the Eddington limit, a neutron star magnetic field strength of $B = 8 \times 10^{10}$ G is estimated. However, a direct measure of $B$ could be obtained observing a CRSF in the energy spectrum. Analysis of \textit{XMM-Newton} data suggested the presence of a cyclotron line at 0.73 keV, with an estimated magnetic field strength of $B=(8.8 \pm 0.3) \times 10^{10}$ G. Here we analyze the 0.3-50 keV broadband spectrum of X 1822-371 combining a 0.3-10 keV NICER spectrum and a 4.5-50 keV \textit{NuSTAR} spectrum to investigate the presence of a cyclotron absorption line and the complex continuum emission spectrum. The NICER spectrum confirms the presence of a cyclotron line at 0.66 keV. The continuum emission is modeled with a Comptonized component, a thermal component associated with the presence of an accretion disk truncated at the magnetospheric radius of 105 km and a reflection component from the disk blurred by relativistic effects. We confirm the presence of a cyclotron line at 0.66 keV inferring a NS magnetic field of $B = (7.9\pm 0.5) \times 10^{10}$ G and suggesting that the Comptonized component originates in the accretion columns.


Introduction
The low-mass X-ray binary system (LMXB) X 1822-371 (4U 1822-37) is a persistent eclipsing source with an orbital period of 5.57 hr, hosting an accreting X-ray pulsar with a spin frequency close to 1.69 Hz (Jonker & van der Klis 2001) and in a spin-up phase with a derivative of ν = (7.39 ± 0.03) × 10 −12 Hz s −1 (Mazzola et al. 2019).X 1822-371 belongs to the class of accretion disc corona (ADC) sources (White & Holt 1982), with an inclination angle between 81 • and 84 • (Heinz & Nowak 2001).The distance to this source was estimated to be between 2-2.5 kpc by Mason & Cordova (1982) using infrared and optical observations.The 0.1-100 keV unabsorbed luminosity is 1.2 × 10 36 erg s −1 , adopting a distance of 2.5 kpc (Iaria et al. 2001b).Recently, using Gaia data, Arnason et al. (2021) estimated a distance to X 1822-371 of 6.1 +1.6  −2.7 kpc, which implies that the observed luminosity should be a factor by six larger than that estimated in the literature.
The most recent orbital ephemeris of the source X 1822-371 was reported by Anitra et al. (2021), who suggested that the orbital period derivative is Ṗorb = (1.42 ± 0.03) × 10 −10 s s −1 adopting a quadratic ephemeris.This high orbital period derivative cannot be explained by a conservative mass transfer at the accretion rate inferred from the observed source luminosity.A highly non-conservative mass transfer is needed at a rate up to seven times the Eddington limit for a neutron star (NS) mass of 1.4 M ⊙ (see Burderi et al. 2010;Bayless et al. 2010;Iaria et al. 2011;Mazzola et al. 2019).In this scenario, the intrinsic luminosity produced by the source is likely at the Eddington limit (a few 10 38 erg s −1 ), two orders of magnitude higher than the observed luminosity of 10 36 erg s −1 .
Assuming the neutron star (NS) is in spin equilibrium, Jonker & van der Klis (2001) showed that the NS magnetic field is close to B ≃ 8 × 10 10 G for an intrinsic luminosity of 10 38 erg s −1 while B ≃ 2 × 10 16 G for an intrinsic luminosity of 10 36 erg s −1 .Finally, they suggested a B value close to 10 12 G, discussing the spectral models reported in the literature up to that time, which yields an intrinsic luminosity of the source of 10 37 erg s −1 given the measured Ṗpulse /P pulse value.
However, a direct measure of the magnetic field strength B is possible only by observing a cyclotron resonance scattering feature (CRSF) in the spectrum.Sasano et al. (2014), studying the Suzaku broadband spectrum of X 1822-371, suggested the presence of a CRSF at 33 keV and a corresponding value of B ≃ 3 × 10 12 G.Iaria et al. (2015) showed that the B-value proposed by Sasano et al. (2014) is not consistent with the spin-up regime of the source.By studying the XMM spectrum of X 1822-371, Iaria et al. (2015) identified a possible CRSF at 0.73 keV corresponding to a value of B = (8.8± 0.3) × 10 10 G.The authors suggested that the observed Comptonized component could be produced in the accretion column onto the NS magnetic caps, and estimated a magnetospheric radius of ∼ 250 km.Iaria et al. (2013) suggested a geometry of the source consisting of an optically thick corona close to the neutron star iden-tified by a Comptonization component in the source spectrum and an optically thin extended corona with an optical depth of τ = 0.01.The optically thin corona was introduced for three reasons: 1) since the intrinsic luminosity is at the Eddington limit while the observed one is a hundred times weaker, an optically thin extended corona with τ = 0.01 would scatter along the line of sight one percent of the intrinsic luminosity (the source is observed with an edge-on geometry); 2) the optically thin corona must be extended because the observed eclipses are partial; 3) the NS spin pulsation is observed.If the corona were extended and optically thick, the pulsation would be diluted by photon scatterings in the corona while assuming the existence of an optically thin corona with τ = 0.01 the travel time of the photons through the optically thin corona (≃ 2 s) is comparable to the spin period (≃ 0.59 s) and part of the pulsed photons can reach the observer after scattering.
Bak Nielsen et al. (2017) proposed that the seed-photons spectrum (composed of a blackbody component plus a power law with a cut-off energy at 10 keV and photon index Γ = 1) is scattered by an extended optically thin corona with an optical depth ranging from 0.01 to 3. According to the authors, their model eliminates the optically thick Comptonizing region and explains the system only with an optically thin (or moderately thick; τ ≃ 0.01 − 3) corona surrounding the X-ray pulsar.
Recently, Anitra et al. (2021) showed that the broadband X-ray spectrum of the source can be fitted with a model composed of: a thermal component associated with an accretion disk emission, a Comptonized component and a reflection component from the accretion disk.The emission from the central region of the system is obscured from direct view by the outer edge of the disk.However, it is scattered along the line of sight by an optically thin corona with optical depth of 0.01.
In this work, we analyze the broadband spectrum of X 1822-371 combining a NICER spectrum with an exposure time of 11 ks and NuSTAR spectrum with an exposure time of 96 ks.We confirm the presence of a cyclotron absorption line close to 0.7 keV by taking advantage of the large effective area of NICER between 0.5 and 2 keV.We confirm that the model discussed by Anitra et al. (2021) fits the broadband spectrum well and suggest that the Comptonizing corona is placed at the accretion column over the magnetic caps of X 1822-371.

Observations and data analysis
Neutron Star Interior Composition Explorer ( NICER, Gendreau et al. 2016) observed X 1822-371 on 2022 Jun 04 for a total exposure time of 11 ks (obsid.5202780101).The primary instrument of NICER is the X-ray Timing Instrument (XTI), which is an array of 56 photon detectors that operate in the 0.2-12 keV energy range.We reduced the NICER data using the pipeline tool nicerl2 in NICERDAS v10 available with HEASOFT v6.31 and adopting the standard filters; the used calibration database version was xti20221001.During the observation the Focal Plane Modules (FPM) 63 and 43 were noisy; the tool nicerl2 rejected all the events collected by FPM 63 and generated GTIs for FPM 43 after data screening.In order to create the spectrum, we utilized the nicerl3-spect tool to exclude FPM 14 and 34 due to elevated detector noise.The 3-10 keV light curve was extracted using the tool nicerl3-lc.The source spectrum and the scorpeon background file were extracted using the tool nicerl3-spect and setting the options bkgformat=file.We added the systematic error suggested by the NICER calibration team to the spectrum1 , and grouped the data using the ftool ftgrouppha applying an optimal rebinning grouptype=optmin to have at least 25 counts per energy bin (groupscale=25, Kaastra & Bleeker 2016).
NuSTAR observed X 1822-371 from 2022 February 15 at 08:21:09, to February 17 at 14:26:19 (ObsID 30701025002), with a duration of 200 ks and a total exposure time of 96 ks.To produce the light curve, the source and the background spectra of FPMA and FPMB, we utilized the Nupipeline and Nuproducts scripts in HEASOFT.A circular region with a radius of 100 arcsec was chosen to extract the events from the source.The background events were extracted from a circular region within the Field of View (FOV) that is devoid of source photons.This circular region was selected from the same detector where the source is located and has a radius of 100 arcsec.After checking that the FPMA and FPMB spectra are similar in the 4-50 keV energy band we summed the two spectra using the ftool addspec2 , the summed spectrum was rebinned as done for the NICER spectrum.

Update of the orbital ephemeris
We produced the NuSTAR light curve extracting the FPMA events in the 5-40 keV energy range and applying the barycentric corrections to the event arrival times using the ftool barycorr.We folded the light curve, adopting a reference epoch T fold =59625 MJD and a reference period of P fold =0.2321107 days (see Fig. 1).To estimate the eclipse arrival time T ecl we used the procedure proposed by Burderi et al. (2010) finding T ecl =59625.19445(48)MJD/TDB.Adopting a reference orbital period of P 0 = 0.232109571 days and a reference eclipse time T 0 =50353.08728MJD, we found that the number of orbital cycles N is 39947 and the delay associated with the eclipse ar- rival time is 2259(41) s.To update the orbital ephemeris, we included the eclipse arrival time shown above with those reported by Mazzola et al. (2019) and Anitra et al. (2021).The delays as a function of orbital cycles reveal that the source exhibits stable expansion over 44 years of data (see top panel of Figure 2).We fitted the delays as a function of cycles adopting the quadratic model y = a + bN + cN 2 , where a = ∆T 0 is the correction to the adopted T 0 , b = ∆P is the correction to the adopted P 0 and c = 1/2P 0 Ṗ allows us to estimate the orbital period derivative Ṗ.We obtained a χ 2 (d.o.f.) of 47(32) (similarly to what obtained by Mazzola et al. 2019).The best-fit quadratic curve (red) and the corresponding residuals in units of σ are shown, respectively, in the top and bottom panels of Figure 2.
The updated orbital ephemeris are: T ecl = 50353.08740(13)MJD/TDB + 0.232109559(6)N+ where the first and the second term represent the refined values of the reference epoch T 0,orb and orbital period P 0,orb , respectively.The third term allows us to estimate an orbital period derivative of Ṗorb = 1.426(26) × 10 −10 s s −1 .We added a cubic term in order to test the presence of a second derivative of the orbital period.By fitting the delays with a cubic function we found a χ 2 (d.o.f.) of 43(31), that is the addition of a cubic term is not statistically significant (σ ∼ 1.3).
Finally, we verified whether a gravitational quadrupole coupling produced by tidal dissipation (Applegate & Shaham 1994) could be detectable in our data as done by Mazzola et al. (2019).To this aim, we added a sinusoidal term to the quadratic model (hereafter LQS model).By fitting the delays with the LQS model we obtained a χ 2 (d.o.f.) of 34(29), indicating that the addition of the sinusoidal term is not statistically significant with a detection of 2σ.
The folded NuSTAR/FPMA light curve, imposing as null phase the passage at the eclipse and adopting the orbital period at the times of the observation, is shown in fig. 3

(black curve).
We applied the barycentric correction to the 3-10 keV NICER light curve using the ftool barycorr and folded it using the eclipse arrival time (T e = 59734.51859MJD/TDB) and the orbital period predicted by the ephemeris shown in eq. 1 (P orb = 20054.36448s).We show the folded light curve in fig.
3 (red curve), the orbital modulation of the two light curves is similar along the orbital period.We searched for the NS spin frequency in the 5-40 keV NuSTAR/FPMA light curve binned at 0.01 s.We corrected the event arrival times for the binary orbital motion using a sin i = 1.006(5) lt-s (Jonker & van der Klis 2001), the eclipse arrival time shown above and the orbital period at the time of the observation.

Spin-period search
We used the ftool efsearch of the XRONOS package, adopting the start time of the observation as reference time and a resolution of the period search of 10 −6 s.We explored around a period of 0.591116 s, estimated using eq.(2) in Mazzola et al. (2019), and subsequently, we fitted the peak of the corresponding χ 2 curve with a Gaussian function.We assumed that the centroid of the Gaussian function is the best estimate of the spin period.We found that the spin period is 0.591124781 s, the χ 2 peak associated with the best period is 135 (see the left panel in  Fig. 4), and the probability of obtaining a χ 2 value greater than or equal to the χ 2 peak by chance, having 15 degrees of freedom, is 2.2 × 10 −21 for a single trial.Considering the 1000 trials in our research, we expect almost 2.2 × 10 −18 periods with a χ 2 value greater than or equal to χ 2 peak.This implies a detection significance with σ ≃ 8.7.
At the light of this result we searched for the spin period in the NuSTAR obsid.30301009002 discussed by Anitra et al. (2021).The FPMA light curve was rebinned at 0.01 s and we selected the events in the 5-40 keV energy range.After applying the barycentric correction to the event arrival times we corrected them for the orbital motion adopting the eclipse arrival time of T ecl = 58234.1536(5)MJD/TDB (see Anitra et al. 2021) and the orbital period at the time of the observation.Using the ftool efsearch, we explored around a period P spin of 0.591472 s obtaining a spin period of 0.59145871 s.The χ 2 peak associated with the best period is 59 (see the right panel in Fig. 4), and the probability of obtaining a χ 2 value greater than or equal to the χ 2 peak by chance, with 15 degree of freedom, is 3.7 × 10 −7 for a single trial.Considering the 1000 trials in our research, we expect almost 3.7 × 10 −4 periods with a χ 2 value greater than or equal to the χ 2 peak.This implies a detection significance with σ ≃ 3.4.
Finally, we searched for periodicity in the NICER data without success.We searched in the 2-10 keV and 5-10 keV energy band because the pulse amplitude increases with increasing energy, it is less than 0.5% between 2 and 6 keV and less than 1% between 6 and 10 keV (see Jonker & van der Klis 2001).We believe that the low effective area of NICER above 5 keV and the short exposure time prevent us from finding a statistically significant signal.
We refined the long-term spin period derivative using the spin period best-values, which span 25 yrs, shown in Tab. 1.We fitted the P spin values with a linear function (see the top panel of Fig. 5) and associated the post-fit errors to the parameters.We found a dependence of P spin on time given by: where the time t is given in MJD and the spin period derivative is Ṗspin = −2.64(2)× 10 −12 s s −1 (ν = 7.55(4) × 10 −12 Hz s −1 ) confirming that X1822-371 is spinning up.Propagating the errors on the P spin from eq.2 and using the eclipse arrival times from the two NuSTAR observations, we obtained the errors on the P spin .We found that P spin = 0.59112(2) s and P spin = 0.59144(2) s for the NuSTAR observations taken in 2022 and 2017, respectively.
Fluctuations are present in the data (see bottom panel in Fig. 5) as already discussed by Bak Nielsen et al. (2017) although the decreasing trend is evident.We show the pulse profiles in Fig. 6.By fitting with a sinusoidal function we found that the pulse amplitude is 1.8% for both the observations.Finally, we note that the pulse profiles are not perfectly sinusoidal and their shape is similar to that shown by Jonker & van der Klis (2001) for the energy band 9.4-22.7 keV.

The NICER spectrum
Initially, we fitted only the NICER spectrum in the 0.3-10 keV energy band with XSPEC v12.13.0c.We adopted the TBabs component accounting for the ISM abundances (Wilms et al. 2000) and set the photoelectric cross section reported by Verner et al. (1996).We adopted a model composed of a Comptonized component (Comptb in XSPEC, Farinelli et al. 2008) and two Gaussian emission lines at 6.4 keV and 6.96 keV.We kept fixed the parameter logA at the value of 8 in the Comptb component.The parameter logA is called the illumination factor, where 1/(1 + A) is the fraction of the seed-photon radiation directly seen by the observer, whereas A/(1 + A) is the fraction upscattered by the Compton cloud.By fixing logA = 8 we are assuming that all the seed-photons are upscattered.Hence, the adopted model is Model A: TBabs*(Comptb+gauss+gauss); this simple model does not give a good fit of the spectrum; we found a χ 2 (d.o.f.) value of 453(144).To improve the fit we added a blackbody component (bbodyrad in XSPEC) to the model (Model B: TBabs*(bbodyrad+Comptb+gauss+gauss)) finding a χ 2 (d.o.f.) value of 258(142).The fit is still statistically unsatisfactory, since large residuals are present below 1 keV, where the data-to-model ratio is larger than 5% (see top panel in Figure 7).A possible explanation for those residuals below 1 keV is that there might be abundances of oxygen and/or iron deviating from solar values.Therefore, we replaced the photoelectric absorption model Tbabs with the Tbfeo model, allowing the oxygen and iron abundances to vary 3 .We fitted the spectrum adopting Model C: TBfeo*(bbodyrad+Comptb+gauss+gauss).
We obtained a χ 2 (d.o.f.) value of 113( 140), the abundances of oxygen and iron were 3.6 ± 0.4 and > 10 times the solar abundances, respectively.Regardless of the plausibility of the abundance values associated with oxygen and iron, we still notice that the ratio exhibits residuals larger than 5% around 0.7 keV (see the second panel from the top in Figure 7).We replaced the Tbfeo component in the model with Tbabs and added two absorption edges (edge in XSPEC) with threshold energies fixed at 0.56 keV and 0.71 keV.The two absorption edges are associated with the neutral K-shell O and neutral L-shell Fe.This choice implies that the neutral iron abundance estimated by the edge at 0.71 keV is different from that estimated by the neutral iron K-shell edge at 7.1 keV with Tbabs component, making the model non-selfconsistent.The sole purpose of this model is to address the residuals below 1 keV.We fitted the spectrum with Model D: TBabs*edge*edge(bbodyrad+Comptb+gauss+gauss).We found a χ 2 (d.o.f.) value of 99( 140).The ratio corresponding to this model is shown in the third panel from the top in Figure 7.We observe that close to 0.7 keV the model deviates from the data by 5%.Since the introduced systematic error is less than 2%, we conclude that this model, besides being non-self-consistent, fails in its purpose of accurately modeling the spectrum below 1 keV.
To take into account a potential misrepresentation of the Solar Wind Charge Exchange (SWCX) background during the observation we fitted the spectrum adopting Model E: TBpcf*TBabs*(ga OK +ga Ne9 +Comptb+gauss+gauss), where the components ga OK and ga Ne9 are two Gaussian emission lines with energy fixed at 0.53 keV and 0.92 keV and associ-3 https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/workshops/NICER-CalStatus-Markwardt-2021.pdfated with neutral oxygen and Ne ix, respectively.Initially, we let the widths of the two lines free to vary and obtained a χ 2 (d.o.f.) value of 116(138).However, the best-fit values of the widths of the K-shell O and the Ne ix lines are 60 ± 15 eV and 280 +100 −50 eV; the broadening of the two lines suggests that the Gaussian components try to fit the continuum emission.In fact, SWCX lines are typically modeled as delta functions (see e.g., Fujimoto et al. 2007), so this model must be rejected.Then, we fitted the spectrum keeping fixed the widths of the ga OK and ga Ne9 components at 5 eV (hereafter Model F).We found a χ 2 (d.o.f.) value of 136(140) but large residuals persist below 1 keV as shown in the bottom panel of Figure 7.
Moreover, we investigated the possibility that the residuals below 1 keV are due to a poor modeling of the SCORPEON background spectrum.Therefore, we extracted the source spectrum and the background model (using the nicerl3-spect tool and setting bkgformat=script) to simultaneously fit the source spectral model and the SCORPEON background model4 .To fit the source spectrum and background together we adopted the pgstat statistic.However, the fit in the 0.3-15 keV energy range leaves the broad residuals below 1 keV unchanged with a ratio larger than 4%, suggesting that the residuals observed at 0.7 keV are not related to a mismodeling of the background.
Before proceeding with the investigation of the X 1822-371 spectrum below 1 keV, we note that the models from C to F (summarized in Table 2) show a reduced χ 2 less than 1 but the corresponding data-to-model ratios exceed 5% below 1 keV (see Figure 7).The reduced χ 2 values should not lead to the conclusion that we found a good fit because they are influenced by the systematic error we added to the spectrum.On the other hand, ratios exceeding 2% below 1 keV indicate the need to add an additional component to obtain a good fit of the spectrum below 1 keV5 .
Analyzing a XMM-Newton observation of X 1822-371, Iaria et al. (2015) suggested the presence of a cyclotron line at 0.7 keV, and we decided to investigate this scenario.We added to Model B a gabs component at 0.7 keV to take into account the possible presence of a cyclotron line.(hereafter Model G: gabs*TBabs*(bbodyrad+Comptb+gauss+gauss)).We obtained a χ 2 (d.o.f.) value of 85(139) with a ∆χ 2 of 172.9 with respect to Model B and the residuals below 1 keV disappear.We show the unfolded model and the corresponding ratio in Figure 8, the best fit values of the parameters associated with this model are shown in the third column of Table 3.
To verify whether the gabs component at 0.7 keV is statistically significant we used Monte-Carlo simulations as shown by Bhalerao et al. (2015).We adopted Model B as our null hypothesis and simulated 1000 fake spectra.Then we fitted each fake spectrum with Model B and Model B plus gabs (i.e.Model G) annotating the ∆χ 2 value obtained for the addition of gabs.Since the cyclotron line adds three free parameters, we expect that the histogram of ∆χ 2 values follows a χ 2 distribution with three degrees of freedom.We show the results in Figure 9.The highest ∆χ 2 obtained in our simulation is 15.8, significantly lower than ∆χ 2 = 172.9obtained in real data.From our simulation the probability to observe a ∆χ 2 of 172.9 by chance is 7.8 × 10 −35 that corresponds to a significance of ∼ 12 σ.
We repeated the same procedure for Model C, D and F. We show the observed ∆χ 2 due to the addition of the gabs compo- Best-fit values adopting the systematic error suggested by the NICER calibration team (third column) and a systematic error of 1% (fourth column).The errors are at 90% c.l.
nent and its statistical significance in the fourth and fifth column of Table 2.The addition of the gabs component to the tested models is always larger than 3 sigmas except for Model D for which we find the significance is only 0.1 σ.So, except for Model D, the addition of a gabs component significantly improves the quality of the fit.

Combined fit of the NICER spectra
To further investigate the possible presence of overabundances of neutral oxygen and neutral iron with respect to the solarabundances values, that is Model C and Model D, we combined the aforementioned NICER spectrum with five additional NICER spectra obtained from five different pointings of the source.
The five NICER observations (obsid.from 5202780102 to 5202780106) have an exposure time of 1.4 ks, 4.3 ks, 5.9 ks, 2.4 ks and 1.7 ks, respectively.The observations were conducted from June 23 to June 26, 2022.During the observation the Fo-cal Plane Module (FPM) 63 was noisy, and during the screening of the data the tool nicerl2 rejected all the events collected by FPM 63.Moreover, in order to produce the spectrum we removed FPM 14 and 34 using nicerl3-spect tool because of increased detector noise.
The source spectrum and the scorpeon background file were extracted using the tool nicerl3-spect and setting the options bkgformat=file.We added the systematic error suggested by the NICER calibration team to the spectrum and grouped the data using the ftool ftgrouppha by applying an optimal rebinning grouptype=optmin to have at least 25 counts per energy bin (groupscale=25, Kaastra & Bleeker 2016).
To fit the six spectra we included to the analyzed models a constant component (const) to account for flux variations of the source in different observations.Additionally, we left the blackbody normalization and the temperature of the seed photons (kT s ) free to vary for all six spectra.Finally, we added an edge at 2.1 keV to address the calibration feature at 2 keV present in the NICER spectra (see Figure 8).
Initially, we fitted the six spectra using Model B. We found a χ 2 (d.o.f.) of 1219(844).Below 1 keV, the model deviates from the data by 10%, as shown in the top-left panel of Figure 10.By adding a gabs component at 0.7 keV we obtained a χ 2 (d.o.f.) of 594(841) and a ∆χ 2 of 625.We show the ratio corresponding to this model in the bottom-left panel of Figure 10, the residuals below 1 keV disappear.
To verify whether the gabs component at 0.7 keV is statistically significant we used Monte-Carlo simulations.We adopted Model B as our null hypothesis and simulated 1000 fake spectra.We show our simulations is the left panel of Figure 11, the largest value of ∆χ 2 obtained from the simulations is 23.2 that is significantly smaller than ∆χ 2 = 625 obtained in real data.The probability to observe a ∆χ 2 of 625 by chance is 1.6 × 10 −132 corresponding to 24.5 σ.The best-fit parameters of the gabs component are E = 0.726 ± 0.008 keV, σ = 0.106 ± 0.012 keV and Strength= 0.064 +0.011  −0.009 .We conclude that the addition of the gabs component to Model B is highly significant.
By fitting the data using Model D We found a χ 2 (d.o.f.) of 657(842).The depths of the edges at 0.56 keV and 0.71 keV are 0.31 ± 0.03 and 0.20 ± 0.02, respectively.Below 1 keV, the model deviates from the data by at least 5%, as shown in the topcenter panel of Figure 10.Adding a gabs component at 0.7 keV we obtained a χ 2 (d.o.f.) of 592(839) and a ∆χ 2 of 65.We show the ratio corresponding to this model in the bottom-center panel of Figure 10, the residuals below 1 keV disappear.The depth of the edge at 0.56 keV and 0.71 keV are < 0.17 and < 0.12, respectively.The best-fit parameters of the gabs component are E = 0.74 +0.02 −0.03 keV, σ = 0.09 ± 0.02 keV and Strength= 0.04 ± 0.02.
Similarly, we assessed the statistical significance of incorporating the gabs component in this instance.We adopted Model D as our null hypothesis and simulated 1000 fake spectra.We show our simulations in the center panel of Figure 11, the largest value of ∆χ 2 obtained from the simulations is 15.7 that is smaller than ∆χ 2 = 65 obtained in real data.The probability to observe a ∆χ 2 of 65 by chance is 2.1 × 10 −11 corresponding to 6.6 σ.We conclude that the addition of the gabs component to Model D is highly significant.
Finally, we fitted the six NICER spectra adopting Model C, where we adopted the Tbfeo component leaving free to vary the O and Fe abundance.We found a χ 2 (d.o.f.) of 713(842); the abundances of O and Fe are 3.4 ± 0.2 and 11.5 ± 1.3 times the corresponding solar-abundance values.However, the ratio in the top-right panel of Figure 10 shows that the model devi-  −0.014 keV and Strength= 0.07 ± 0.03.Likewise, we evaluated the statistical significance associated with the inclusion of the gabs component in this case.We adopted Model C as our null hypothesis and simulated 1000 fake spectra.We show our simulations is the right panel of Figure 11, the largest value of ∆χ 2 obtained from the simulations is 17.2 that is smaller than ∆χ 2 = 119 obtained in real data.The probability to observe a ∆χ 2 of 119 by chance is 5.7 × 10 −23 corre-sponding to 9.8 σ.We conclude that the addition of the gabs component to Model C, is highly significant.So, the inclusion of the gabs component is always necessary in models B, C, and D, as the addition of this component exhibits a significance exceeding 5 sigma for all the three examined models.

Combined fit of the NICER spectrum with two XMM/EPIC-pn spectra
To further enhance the robustness of our detection of the cyclotron line at 0.7 keV and to better investigate whether there are two edges associated with the K-shell of neutral oxygen and the L-shell of neutral iron below 1 keV, we reanalyzed two XMM/EPIC-pn spectra combining them with the NICER spectrum.The obsid.0784820101 was already analyzed by Anitra et al. (2021).The authors analyzed the EPIC-pn spectrum in the 2-10 keV energy band suggesting calibration issues in the spectrum associated with the presence of a silicon edge at 1.8 keV.The source was observed by the XMM-Newton observatory on 2017 March 3 for a duration of 69 ks.During the observation the EPIC-pn camera was operating in Timing Mode.We reduced the XMM-Newton data using the Science Analysis Software (SAS) v21.0.0.Initially, we extracted the EPIC-pn events between 10 and 12 keV to verify the presence of particle flaring background6 .The light curve shows a strong flaring activity after the first 30 ks from the beginning of the observation.To mitigate the effects of particle flaring, we selected the time intervals with a rate lower than 1.5 c/s in the energy range 10-12 keV.
We show the source and background EPIC-pn light curve in the 0.6-10 keV energy range cleaned from the particle flaring background in Figure 12 (left panel).The source light curve clearly shows the orbital modulation and ranges between 40 c s −1 and 90 c s −1 , two partial eclipses are present at 14000 s and 34000 s from the beginning of the observation.
The EPIC-pn spectrum and background were extracted selecting RAWX between 26 and 47 and RAWX between 3 and 5, respectively.The EPIC-pn spectrum has an exposure time of 41 ks, it was grouped using the ftool ftgrouppha by applying an optimal rebinning grouptype=optmin to have at least 25 counts per energy bin (groupscale=25, Kaastra & Bleeker 2016).To fit the data, we adopted the range of 0.6-10 keV for the EPIC-pn spectrum.
The obsid.0111230101 was already analyzed by Iaria et al. (2015); the authors suggested the presence of an absorption line, discussed as cyclotron absorption line, close to 0.7 keV by analyzing the EPIC-pn and RGS spectra in the 0.6-10 keV and 0.35-2 keV energy band, respectively.The observation was carried out on 2001 March 7 for a duration of 54 ks.During the observation the EPIC-pn camera was operating in Timing Mode.We reduced the XMM-Newton data using the Science Analysis Software (SAS) v21.0.0.
To mitigate the effects of particle flaring background, we selected the time intervals with a rate lower than 1.5 c/s.We show the source and background EPIC-pn light curve in the 0.6-10 keV energy range cleaned from the particle flaring background in Fig. Figure 12 (right panel).The source light curve shows the orbital modulation and ranges between 30 c s −1 and 70 c s −1 , two partial eclipses are present close to 20000 s and 40000 s from the beginning of the observation.
The EPIC-pn spectrum and background were extracted selecting RAWX between 31 and 44 and RAWX between 3 and 5.The EPIC-pn has an exposure time of 41 ks, it was grouped as shown above.
The two XMM observations are distant in time, we produced the hardness ratios (HRs) using the energy band 0.6-2 keV and 2-10 keV to check a possible change in the shape spectrum.The HRs are shown in Fig. 13.We observe that the HR value is approximately 2 for the observation carried out in 2001 and approximately 1.7 for the one carried out in 2017.
Since the NICER spectrum alone may suggest that the addition of two edges at 0.56 and 0.71 keV does not require the presence of an additional cyclotron line at 0.7 keV although we note that the data-to-model ratio below 1 keV is approximately 5%, we decide to verify that this continues to be true when we perform a combined fit of the NICER and the two EPIC-pn spectra.At 0.7 keV, the effective area of EPIC-pn is approximately 900 cm 2 , which is only 30% lower than that of NICER, which is around 1250 cm 2 .Given that the NICER spectrum has a systematic error larger than 1% and an exposure time of 11 ks, while both EPIC-pn spectra have an exposure time of 41 ks, we can compare the statistics of the three spectra at 0.7 keV without the risk of biasing the fit toward the NICER spectrum.
Initially, we fitted the NICER spectrum (0.3-10 keV), the 0.6-10 keV EPIC-pn spectrum corresponding to obsid.0784820101 (hereafter PN1 spectrum), the 0.6-10 keV EPIC-pn spectrum corresponding to obsid.0111230101 (hereafter PN2 spectrum) using Model B, to which we added a Gaussian component in the Fe-K region of the spectrum because PN1 and PN2 spectra require it.We added a constant to accommodate the different normalization between the three spectra.The constant value was held fixed at 1 for the NICER spectrum, while it was allowed to vary for PN1 and PN2 spectra.To take into account the systematic feature present close 2 keV in the NICER spectrum (see Figure 8) we added an absorption edge with energy threshold fixed at 2.1 keV; we kept fixed to zero the depth of the edge for the other spectra.In addition, since the three spectra were extracted from non-simultaneous observations, we allowed the normalization of the blackbody component, the seed-photon temperature and the energy index α of the Comptb component to vary independently.We forced the value of the electronic temperature of the PN1 spectrum to that of the NICER spectrum.We left to vary independently the energies and the normalizations of the two narrow Gaussian line at 6.4 and 6.96 keV for PN1 spectrum and kept fixed at 70 eV their widths.Finally, we added two Gaussian emission lines at PN1 and PN2 spectrum at 1.18 keV and 1.343 keV with width fixed at 5 eV.To take into account of the systematic at 2 keV in the PN2 spectrum we added a Gaussian component with null width and negative normalization (hereafter we call this model Model H).
The N H value is in units of 10 20 cm −2 .
χ 2 (dof) of 494(356) with a ∆χ 2 = 843.The ratio is within 2% as shown in the bottom-left panel of Figure 14.The best-fit parameters are shown in the third column of Table 4.The best-fit parameters of the gabs component are E 0 = 0.728 ± 0.005 keV, σ 0 = 0.097 ± 0.008 keV and Strength= 0.059 ± 0.009, these results are consistent with those obtained by fitting the NICER spectrum alone.
To test the significance of adding the gabs component, we simulated 1000 spectra using the Monte-Carlo method as described in the previous section.The histogram of the ∆χ 2 values obtained from the simulation is shown in the left panel of Figure 15 and it is well-fitted by a χ 2 with three d.o.f.(red curve).We found that the largest value of ∆χ 2 obtained from the simulation is 15.9, the probability to observe a ∆χ 2 of 849 by chance is 5.2 × 10 −180 corresponding to a statistical significance of 29 σ.
We started from Model H adding to it two absorption edges at 0.56 keV and 0.71 keV.By fitting the three spectra we found a χ 2 (dof) of 615(357), the data deviate from the model up to 5% below 0.6 keV as shown in the top-center panel of Figure 14.We added to the model a gabs component at 0.7 keV (hereafter Model L).The best-fit improves significantly; we found a χ 2 (dof) of 484(354) with a ∆χ 2 = 131.We show the residuals in the bottom-center panel of Figure 14.The best-fit parameter are shown in the fourth column of Table 4.The best-fit parameters of the gabs component are E 0 = 0.749 +0.012  −0.027 keV, σ 0 = 0.075 ± 0.013 keV and Strength= 0.033 +0.013 −0.009 , these results are consistent with those obtained by adopting Model I.
Also for this model, we assessed the significance of adding the gabs component as done previously.The histogram of the simulated spectra is shown in the middle panel of Figure 15.We found that the largest value of ∆χ 2 obtained from the simulation is 16.7, the probability to observe a ∆χ 2 of 131 by chance is 8.6 × 10 −26 corresponding to a statistical significance of 10.5 σ.Our test indicates that the addition of a component representing a cyclotron line at 0.7 keV is necessary, even though two edges are included in the model at 0.56 and 0.71 keV.
Finally, instead of adding two edges at 0.56 and 0.71 keV we used a self-consistent absorption model.Specifically, we replaced TBabs with TBfeo, allowing the abundances of oxygen and iron to vary freely, as suggested by the NICER team.By fitting the spectra we found a χ 2 (dof) of 658(357); also in this case the data deviate from the model up to 5% below 0.6 keV (top-right panel of Figure 14).We added to the model a gabs component at 0.7 keV (hereafter Model M).We found a χ 2 (dof) of 489(354) with a ∆χ 2 = 169.The best-fit parameters are shown in the fifth column of Table 4, the residuals in the bottom-right panel of Figure 14).The best-fit parameters of the gabs component are E 0 = 0.744 +0.011  −0.018 keV, σ 0 = 0.081 +0.016 −0.013 keV and Strength= 0.040 +0.014 −0.010 , these results are consistent with those obtained by adopting Model I and Model L. In this case, we find that the largest value of ∆χ 2 obtained from the simulation is 17.2, the probability to observe a ∆χ 2 of 169 by chance is 5.2 × 10 −34 corresponding to a statistical significance of 12 σ (right panel of Figure 15).Since the component associated with a 0.7 keV cyclotron line is required for all three described models, we discuss a model that requires a cyclotron line at low energies, with abundances of O and Fe fixed to solar values (Model G) in the subsequent sections of the work.
However, since the reduced-χ 2 associated with the NICER spectrum alone is much smaller than one and attributable to an overestimation of systematic errors in the 0.3-10 keV range (see the third column in Table 3), we extracted the spectrum by assigning a 1% systematic error to our data.Fitting the NICER spectrum adopting Model G, we obtained a χ 2 of 117( 138) with the ratio between 0.3 and 10 keV below 2%.The best-fit values shown in the fourth column of Table 3 indicate that the spectral shape remains unchanged when transitioning from the systematic error suggested by the NICER calibration team to a 1% systematic error, while the reduced-χ 2 increases from 0.61 to 0.87.

Combined NICER-NuSTAR spectra
We fitted the combined 0.3-10 keV NICER spectrum and the 4-50 keV NuSTAR spectrum using Model G where we left free to vary the parameter logA of the Comptb component.To take into account the systematic feature present close 2 keV in the NICER spectrum (see Figure 8) we added an absorption edge with energy threshold fixed at 2.1 keV; we kept fixed to zero the depth of the edge for the NuSTAR spectrum.In addition, since the two spectra are extracted from non-simultaneous observations, we allowed the energy index α of the Comptb component to vary independently between the two spectra.We added a constant to account for the different normalization between the two spectra.The constant value was kept fixed at 1 for the NICER spectrum, while it was free to vary for the NuSTAR spectrum.Finally, we kept fixed at 70 eV the widths of the two lines in the Fe-K region.
Fitting the spectrum, we obtained an unacceptable fit with a χ 2 (dof) of 1038(324).In the NuSTAR spectrum, a feature at 9.6 keV is evident, which we modeled with an absorption edge.We set the depth of the edge to zero for the NICER spectrum.The addition of the absorption edge improves the fit; we obtained a χ 2 (dof) of 905(322).However, large residuals were still present in the Fe-K region of the spectrum.For this reason, we added an additional Gaussian component to the model.The model (hereafter Model 1) is: We found a χ 2 (d.o.f.) value of 469(319), the added Gaussian emission line has a centroid of 6.13±0.08keV and a σ of 730±90 eV.To better constrain this broad line, we left the width of one of the two narrow lines free to vary and constrained the amplitude of the second narrow line to vary with the first one.The obtained fit was statistically equivalent; we found a χ 2 (dof) of 469(318).The best-fit parameters are reported in the third column of Table 5, the corresponding residuals are shown in the bottom-left panel of Figure 16.Analyzing the residuals, it is evident that the NuSTAR spectrum (in red) is not well modeled above 15 keV.We added a power-law component with a cut-off at low energies to fit the residuals at high energies.We fixed the energy of the cut-off component (expabs in XSPEC) to the seed-photon temperature where po indicates the power-law component.By adding the absorbed power-law component we obtained a χ 2 (d.o.f.) of 302(316) and the residuals above 15 keV disappeared.We show the best-fit parameters in the fourth column of Table 5, the unfolded spectrum and the residuals are show in the top-left and middle-left panel of Figure 16.
In the unfolded spectrum we show the blackbody, the Comptonized component and the power-law in red, blue and green color, respectively.The narrow emission lines are in orange color and, finally, the broad Gaussian line is shown in magenta.
The temperature of the thermal emission is 0.173±0.007keV, the Comptonized component has a seed-photon temperature of 1.34 ± 0.05 keV, an electron temperature of 4.22 ± 0.07 keV, a αindex < 0.13 and 0.22 ± 0.06 for NICER and NuSTAR spectrum respectively.The logA parameter is 0.20 ± 0.10 implying that (39 ± 6)% of the seed photons is directly seen by the observer while the remaining part is scattered in the Comptonizing cloud.The narrow emission lines are at 6.39 ± 0.02 keV and 6.91 ± 0.03 keV with a width of 70 +40 −30 eV.The first line is associated with an emission line of neutral iron, the second line is marginally compatible with the presence of Fe xxvi ions.
The broad emission line has an energy of 6.19 +0.10 −0.34 keV and a width of 530 +90 −250 eV.The photon-index of the power-law component is 1.96 ± 0.05.The gabs component at low energies, that is the cyclotron line, has a centroid at 0.724 ± 0.014 keV, a width of 0.10 ± 0.02 keV and an optical depth of 0.25 ± 0.10.Finally, an absorption edge at 9.62 ± 0.09 keV is observed.
The unabsorbed fluxes in the 0.1-100 keV energy range are 1.3 × 10 −9 erg s −1 cm −2 , 2.6 × 10 −11 erg s −1 cm −2 and 1.4 × 10 −10 erg s −1 cm −2 associated with the Comptonized component, the blackbody and the power-law component, respectively.Assuming a distance to the source of 2.5 kpc, the total unabsorbed luminosity in the 0.1-100 keV energy range is 1.1 × 10 36 erg s −1 (6.6 × 10 36 erg s −1 for a distance to the source of 6.1 kpc, recently estimated by Arnason et al. (2021).Hereafter, we assume that the distance to the source is 6.1 kpc.
Following a scenario describing the presence of cyclotron line at 0.7 keV we changed the gabs component with a cyclabs component.We obtained a statistically equivalent fit with a χ 2 (d.o.f.)= 304(316); the best-fit values associated with the cyclabs component are: E = 0.69 ± 0.02 keV, width= 0.12 ± 0.03 keV and depth= 0.31 ± 0.04.These values are compatible with those reported by Iaria et al. (2015).
Because we are observing a broad Gaussian line in the Fe-K region of the spectrum we explored a scenario in which a reflection component from the accretion disk is present as suggested by Anitra et al. (2021).To model the spectrum we changed the broad Gaussian line with a reflection component blurred by relativistic effects.We assumed that the incident radiation comes from the Compton cloud.The model (hereafter Model 3) is: where rfxconv is the component that takes into account the reflection component from the accretion disk while rdblur is the component that takes into account the relativistic blurring.We The errors are at 90% confidence level.The values in square brackets were kept fixed during the fit.The intensity of the lines is in units of photons/cm 2 /s.fixed the outer radius R out at 3000 gravitational radii after checking that the χ 2 value was not sensible to this parameter.The inclination angle was fixed at 82.5 • as reported in literature (Mason & Cordova 1982).
We show the best-fit values in the fifth column of Table 5, the unfolded spectrum and the residuals are shown in the top-right and bottom-right panels of Figure 16.In the unfolded spectrum of Model 3 the red, blue, magenta and green component are the blackbody, the Comptonized component, the reflection component and the power-law.The two narrow emission lines in the Fe-K region are in orange.

Discussion
Using the NuSTAR data (obsid.30701025002), we updated the orbital ephemeris of X 1822-371 by adding one eclipse arrival time.We find a Ṗorb of 1.426(26) × 10 −10 s s −1 , compatible with the values present in literature (see, e.g., Anitra et al. 2021).As discussed before, the addition of a cubic term or a sinusoidal term to the quadratic model does not improve significantly the fit.
We estimated the spin period of the source during the two NuSTAR observations discussed above and confirmed the presence of a linear long-term decrease of the spin period finding Ṗspin of −2.64(2) × 10 −12 s s −1 , compatible with the value reported by Bak Nielsen et al. (2017).
The analyzed energy spectrum was obtained by combining the 0.3-10 keV NICER spectrum (obsid.5202780101) with the 4.5-50 keV NuSTAR spectrum.The continuum emission is composed of a blackbody component plus a Comptonized component and a power-law component.The blackbody component has a temperature of 0.158 ± 0.011 keV; by assuming a distance to the source of 6.1 kpc we infer an emission radius of 43 +8 −5 km.The seed-photon temperature and the electron temperature of the Comptonized component are kT 0 = 1.32 +0.03 −0.05 keV and kT e = 4.34 ± 0.07 keV, respectively.The illumination factor logA is 0.35 ± 0.12 and the α value, the energy index of the Comptonization spectrum, is 0.19 +0.05  −0.07 and 0.35±0.07for NICER and NuSTAR spectrum, respectively.We find that the percentage of seed-photon radiation directly seen by the observer is 31 ± 6% suggesting that the Comptonizing corona is geometrically compact and does not completely surround the inner regions of the system.The seed-photon emission radius is 2.4 ± 0.2 km.
The gabs component is interpreted as cyclotron absorption line (see also Iaria et al. 2015), with a centroid energy of E 0 = 0.66 ± 0.03 keV.The relation between the energy of the fundamental harmonic and the magnetic field strength is where m is the NS mass in unit of solar masses and R 6 is the NS radius in units of 10 6 cm.The NS mass was estimated to be between 1.61 ± 0.15 and 1.70 ± 0.13 M ⊙ (Iaria et al. 2015).Assuming a NS radius of 10 km and the value of E 0 , we estimate that the magnetic field strength is B = (7.9± 0.5) × 10 10 G and B = (8.1 ± 0.5) × 10 10 for a NS mass of 1.61 ± 0.15 M ⊙ and 1.70 ± 0.13 M ⊙ , respectively.These values of the magnetic field strength are compatible with the estimate done by Jonker & van der Klis (2001) who calculated a value of B ∼ 8 × 10 10 G for a luminosity of the source of 10 38 erg s −1 .In the following we discuss our results for a NS mass of 1.61 M ⊙ (similar results are obtained for a NS mass of 1.70 M ⊙ ).
We observe the weakest magnetic field strength ever inferred from an electron cyclotron line.Since the formation of a cyclotron line requires that the magnetic field is B > (v/c) 2 B crit , where v/c is the velocity of electrons (in units of the speed of light) perpendicular to the field lines and B crit = 4.4 × 10 13 G, we expect that the electron velocity component perpendicular to the field lines is of the order of 10% the speed of light to obtain a threshold value of B of the order of 8 × 10 10 G.
To estimate the expected luminosity from the source we used the Ghosh & Lamb (1979) equation, which links the spin-period derivative, its magnetic field and its luminosity: where I 45 is the NS moment of inertia in units of 10 45 g cm 2 , R 6 is the NS radius in units of 10 6 cm, m is the NS mass in units of solar masses and P is the spin period of the source, that is, 0.591124781(13) s.The parameter n(ω s ) is the dimensionless accretion torque that is a function of the fastness parameter ω s .
When ω s < 0.95, we can use the following approximate expressions (Ghosh & Lamb 1979): (4) Iaria et al. (2015) showed that ω s is between 0.063 and 0.083 for a NS mass of X 1822-371 between 1 and 3 M ⊙ .Assuming a magnetic field strength of B = (7.9± 0.5) × 10 10 G, a NS mass of 1.61 M ⊙ , a NS radius of 10 km and a spin-period derivative Ṗspin of −2.64(2) × 10 −12 s s −1 we infer an intrinsic luminosity of 1.5 × 10 38 erg s −1 (corresponding to a mass accretion rate of 1.13 × 10 −8 M ⊙ yr −1 ) that is very close to the Eddington luminosity of ∼ 2 × 10 38 erg s −1 .The intrinsic luminosity is 20 times larger than the observed luminosity (6.6×10 36 erg s −1 ).This discrepancy is likely due to the high inclination angle of the source, between 81 • and 84 • (see Iaria et al. 2013;Anitra et al. 2021) which causes the emission from the central part of the system to be obscured by the outer edge of the disk.Although the discrepancy between the intrinsic luminosity of the source and the observed luminosity has decreased for a source distance of 6.1 kpc, the scenario described above still holds true.As the ratio between intrinsic and observed luminosity is approximately 20, the optical depth of the ADC should be around 0.05.Therefore, the ADC remains optically thin, and all the previous conclusions remain unchanged.The flux associated with the blackbody component corrected by a factor ∼ 20 gives a blackbody radius of 190 +40 −20 km and a seed-photon radius of 10.7±0.9 km compatible with an emission coming from the NS.We hence associate the blackbody component to the emission from the accretion disk, its radius indicating the inner radius of the disk.We do not correct this value for the inclination angle since we assume that the optically thin corona covers the whole accretion disk.We infer the magnetospheric radius R m using the formula in Sanna et al. (2017), where G is the gravitational constant, M is the NS mass, µ the NS magnetic dipole moment, and Ṁ is the mass accretion rate.The factor ϕ is given by the following equation  (Sanna et al. 2017), where µ 26 , R 6 and Ṁ−9 are the NS magnetic moment in units of 10 26 G cm 3 , the NS radius in units of 10 6 cm and the mass accretion rate in unit of 10 −9 M ⊙ yr −1 .The mean molecular weight κ is 0.615 for fully ionized matter, while the parameter α for a standard Shakura-Sunyaev disk model is set to 0.1 (Shakura & Sunyaev 1973).Adopting a NS mass of 1.61 M ⊙ , a B-field strength of (7.9 ± 0.5) × 10 10 G and a NS radius of 10 km, we obtain that R m = 105 ± 5 km, that is smaller than the inner radius of the accretion disk, as expected for an X-ray pulsar.
We observe a reflection component from the accretion disk with a ionization parameter of log(ξ) = 1.4 +0.4  −1.4 , compatible with the value obtained by Anitra et al. (2021).This value implies a relatively low ionization state of matter in the disk; the iron line is dominated by the Fe I-XX Kα transition (Fabian et al. 2000).The inner radius at which the reflection component is produced is 62 +21 −13 gravitational radii, that corresponds to 150 +50 −30 km for a NS mass of 1.61 M ⊙ .This value is compatible with the inner radius of the accretion disk of 190 +40 −20 km derived above.The Comptonized component is associated with a opticallythick corona.By using eq. 10 in Farinelli et al. (2008) we find that the optical depth τ is 24 and 17 for the NICER and NuSTAR spectrum assuming a spherical geometry and 12 and 8 assuming a slab geometry.
From these findings, we can conclude that although there is a magnetosphere in place and the matter being accreted is channeled along the magnetic field lines, a portion of the matter has sufficient kinetic energy to overcome the magnetic barrier imposed by the magnetosphere.As a result, this matter penetrates the magnetosphere (possibly due to Rayleigh-Taylor and Kelvin-Helmholtz instabilities) and falls directly onto the surface of the NS.This phenomenon is likely favored by the large accretion rate, reaching the Eddington limit, and the magnetic field being relatively weak with a strength of 7.9 × 10 10 G.
We therefore suggest that the Comptonized component is formed at the accretion column onto the magnetic caps and the seed-photons come from all the NS surface, in agreement with the observation that 31% of the seed-photons do not interact with the Comptonized region.
If we assume that the broadening of the cyclotron line (σ 0 = 0.17 ± 0.03 keV) has thermal origin, we estimate a temperature of 13 keV not consistent with the electron temperature of 4.3 keV obtained from the fit.A possible explanation is that we are observing a gradient of the B-field strength along the accreting column.In this hypothesis we can infer the height of the column considering that the magnetic dipole moment µ = BR 3 does not change, and then −∆B/B = 3∆R/R.By imposing ∆B as the sigma of the gabs component we find that ∆R ≃ 0.9 km, assuming a NS radius of 10 km.
Assuming a fully ionized hydrogen plasma (n e = n p ) in the accretion column, the proton density is given by n p ∼ 10 21 Ṁ18 /(βS 10 ) cm −3 (see eq. 1 in Mushtukov et al. 2019), where Ṁ18 is the mass accretion rate in units of 10 18 g s −1 , S 10 is the cross section of the accretion column in units of 10 10 cm 2 and β = v/c with v the local velocity of plasma and c the speed of the light.Since the accreting matter slows down from a free-fall velocity v f f to about 0.1v f f in the radiation dominated shock at the top of accretion column we assume v = 0.1v f f , that is β ≃ 0.07.The cross section of the accretion column, i. e. the area on the star over which the accretion occurs, is S ≃ R 3 R −1 m cm 2 (Pringle & Rees 1972), where R is the NS radius.For a NS radius of 10 km we obtain S ≃ 8.3 × 10 10 cm 2 .Using the mass accretion rate estimated from eq. 3 we find that n p ≃ 1.2 × 10 21 cm −3 .Since n e = n p we can estimate the thickness l of the Comptonized region in the accretion column using τ = σ T n e l where σ T is the Thomson cross section.In the reasonable assumption that the accretion column has a slab geometry we can adopt an optical depth τ = 8 estimated above.We find l ≃ 0.1 km, that is the Comptonizing region has a thickness that is 1/8 of the height of the accretion column.
Finally, Jonker & van der Klis (2001) observed that the pulse fraction increases going from 5 keV to 20 keV, this can be explained in our scenario considering that the spin modulation is associated with the accretion column from which upscattered photons escape.The observed pulsation is due to the Comptonized photons leaving the accretion column which are modulated by the NS rotation.In this scenario, we do not observe a spin modulation at energies below 5 keV because this range is dominated by the soft thermal components corresponding to the inner disk and the seed-photons that come from all over the NS surface.
Our model suggests that the Comptonization component primarily originates from the accretion column.From our best-fit to the data, we find that the reflection amplitude is f re f l = 0.66 +0.08 −0.10 (a comparable value of 0.61 +0.04 −0.05 was obtained by Anitra et al. 2021).Typical values of the reflection amplitude Ω/2π for NS-LMXB atoll sources range from 0.2-0.3(see, e.g., Di Salvo et al. 2015Salvo et al. , 2019;;D'Aì et al. 2010;Egron et al. 2013;Matranga et al. 2017;Marino et al. 2019) and suggest the presence of a spherical corona in the inner part of the accretion disk (see Fig. 5 in Dove et al. 1997).The large reflection amplitude could be due to the peculiar geometry, in our scenario the radiation incident the disk (truncated at the magnetospheric radius) comes from a compact region located at the surface of the neutron star.Iaria et al. (2013) proposed that the Fe i and the Fe xxvi emission lines are produced in the photoionized surface of the accretion disk at a distance from the NS of 2×10 10 cm and < 3.7×10 10 cm, respectively.However, the ionization parameter ξ in the inner region of the accretion disk, where the reflection component originates, is close to 25, which is too low to produce Fe xxvi ions.We suggest that these narrow lines are produced at large distance from the central source, possibly in a disk atmosphere with a large gradient in temperature and/or photoionization, although further investigations are necessary to address their origin.
The absorption edge at 9.62 ± 0.09 keV is not compatible with a Fe xxvi transition expected at 9.278 keV in a rest frame.If we interpret this feature as blue-shifted we derive a wind speed of 11,000 km/s.However, Anitra et al. (2021) suggested that this feature could have a systematic origin.
Finally, we observe the presence of a power-law component never observed before in this source.This component looks similar to the power-law hard tail often observed in the X-ray spectra of bright LMXBs, both in Z-sources (see, e.g., Di Salvo et al. 2000;Iaria et al. 2001a;Di Salvo et al. 2001) and in atoll sources in their soft states (see, e.g., Piraino et al. 2007, and references therein).When significantly detected, these components usually show a power-law shape with a photon index ∼ 2−3 contributing to a few percent of the total flux from the source, and sometimes a correlation with the radio emission has been observed (e.g., Homan et al. 2004;Migliari et al. 2007).The origin of these features is still debated, but it is probably related to the presence of electrons with a non-thermal velocity distribution (e.g., Di Salvo et al. 2006).Considering that X 1822-371 is an X-ray pulsar, we tentatively suggest that non-thermal electrons moving along the magnetic field lines towards the NS polar caps may be responsible for the observed power-law hard tail in this source.

Conclusions
We analyzed a 0.3-50 keV broadband spectrum of X 1822-371 combining a NICER spectrum (0.3-10 keV) with a NuSTAR spectrum (4.5-50 keV).Our best-fit to the spectrum gives clear evidence of an absorption cyclotron line with energy close to 0.66 keV, confirming the detection reported by Iaria et al. (2015).
From the temporal analysis we find that: -The orbital period is expanding at a constant derivative of Ṗorb = 1.426(26) × 10 −10 s s −1 , the NS is spinning up with a long-term spin-period derivative of Ṗspin = −2.64(2)× 10 −12 s s −1 (with P spin = 0.59112(2) s), the observed absorption cyclotron line close to 0.66 keV corresponds to a magnetic field strength of B = (7.9±0.5)×10 10 G for a NS mass of 1.61 M ⊙ , using the Ṗspin and B values we find that the luminosity of the source is 1.5×10 38 erg s −1 and the magnetospheric radius R m = 105 ± 5 km, because the observed 0.1-100 keV luminosity is a factor 20 smaller than 1.5×10 38 erg s −1 and the inclination angle of the system is between 81 • and 84 • , we suggest the presence of an extended optically thin corona (with optical depth of 0.05) that scatters along the line of sight about 5% of the radiation produced in the inner region of the system.
By correcting the observed flux by a factor 20 we associate the thermal component with a thermal emission from the accretion disk finding an inner radius of 190 +40 −20 km.The reflection component is produced at an inner radius of 150 +50 −30 km, that is compatible with the inner radius of the accretion disk.
The Comptonized region is associated with the accretion column at the magnetic caps.We find that the equivalent spherical radius of the seed-photons emission region is 10.7 ± 0.9 km, that is very similar to the NS radius; we find that 31% of these soft photons do not interact with the Comptonizing region.
The broadening of the cyclotron line can be explained with a gradient of the magnetic field B associated with a height of the accretion column of ∼ 0.9 km.We estimate an average electron density of 1.2 × 10 21 cm −3 in the accretion column.
We observe a reflection component with an amplitude of 0.66 +0.08 −0.10 , about a factor of ∼ 2 larger than the typical values of 0.2 − 0.3 observed in the NS-LMXBs.A possible explanation is that the incident radiation to the disk comes from the photons leaving the accretion column with a geometry roughly similar to the lamp-post geometry observed in AGNs.
Finally, we suggest that the two observed narrow lines associated with the presence of neutral iron and Fe xxvi ions form in a region of the system, possibly a disk atmosphere, likely at a large distance, ∼ 2 × 10 10 cm, from the central source.

Fig. 1 .
Fig. 1.NuSTAR/FPMA folded light curve assuming as epoch T fold = 59625 MJD and an orbital period of 0.2321107 days.The period is divided into 512 bins.For clarity, two orbital phase cycles are shown.

Fig. 2 .
Fig. 2. Delays vs. orbital cycles for the quadratic ephemeris (top panel).We used the eclipse arrival times shown by Mazzola et al. (2019) (see tab. 1 in their work), the one obtained by Anitra et al. (2021) and that shown in the text.Residuals in units of σ (bottom panel).

Fig. 5 .
Fig.5.Long-term spin period change (top panel) and corresponding residuals in units of µs.The linear decrease of the pulse period over time suggests that the NS is in a spin-up regime.
Fig. 6. 5-40 keV NuSTAR/FPMA folded light curve of the obsid.30701025002 (left panel) and obsid.30301009002 (right panel).The curves are folded to have 16 bins per period.Two periods are shown fro clarity.The red curves represent the best-fit model using a sinusoidal function.

Fig. 7 .
Fig. 7. 0.3-10 keV ratios for Models B, C, D and F. The systematic error is that suggested by the NICER calibration team.The red horizontal lines indicate a deviation of the model from the data of 2%.

Fig
Fig. 8. 0.3-10 keV unfolded spectrum using Models G (top panel); the blackbody, the Comptonized component and the emission lines at 6.4 keV and 6.96 keV are shown with the red, blue and orange colors, respectively.The corresponding data-to-model ratio is in bottom panel.
Fig.9.Monte-Carlo simulations results for testing the gabs significance for Model B. The blue histogram shows ∆χ 2 values obtained in the simulations, and the red curve is a χ 2 distribution with 3 d.o.f..The ∆χ 2 = 172.9attained in actual data is significantly higher than values attained in simulations (by ∼ 12σ).

Fig. 10 .
Fig.10.Data-to-model ratio of the six NICER (0.3-10 keV) spectra with respect to the different adopted models, with or without the gabs component.The adopted colors for each spectrum are shown in the legend in which we indicate the last three digits of the corresponding observation ObsId.In the top panels we show the ratio corresponding to Model B, Model C, and Model D, respectively.In the bottom panels we show the ratio of the same models including a gabs component at 0.7 keV.The black horizontal lines indicate a deviation of the model from the data corresponding to 2%.

Fig. 14 .
Fig.14.Ratios combining the 0.3-10 keV NICER (black), PN1 (red) and PN2 (green) spectrum.In the top panels the residuals corresponding to Model I, Model L, and Model M removing the gabs component.In the bottom panels the residuals associated with the best-fit model shown in Table4.The orange horizontal lines indicate a deviation of the model from the data larger than 2%.

Fig. 16 .
Fig. 16.Unfolded spectrum of Model 2 (top-left panel) using a 0.3-10 keV NICER/XTI spectrum (black color) combined with a 4-50 keV NuSTAR/(FPMA+FPMB) spectrum (red color), the colors of the components are described in the text.Residuals adopting Model 2 and adopting Model 1 are shown in the middle-left panel and bottom-left panel, respectively.Unfolded spectrum and residuals of Model 3 in top-right and bottom-right panel, respectively.

Table 2 .
Adopted models to fit the NICER spectrum Obtained adding to the model a gabs component at 0.7 keV.b Statistical significance for the addition of the gabs component at 0.7 keV.c,d,e Best-fit values of the energy, width and strength of the gabs component, errors are at 90% confidence level.f The widths of the neutral O and Ne ix emission lines were left free to vary.g The widths of the neutral O and Ne ix emission lines were kept fixed at 5 eV. a

Table 3 .
Best-fit values of the parameters associated with Model G

Table 5 .
Best-fit parameters combining NICER and NuSTAR spectra