The High Energy X-ray Probe (HEX-P): A New Window into Neutron Star Accretion

Accreting neutron stars (NSs) represent a unique laboratory for probing the physics of accretion in the presence of strong magnetic fields ($B\gtrsim 10^8$ G). Additionally, the matter inside the NS itself exists in an ultra-dense, cold state that cannot be reproduced in Earth-based laboratories. Hence, observational studies of these objects are a way to probe the most extreme physical regimes. Here we present an overview of the field and discuss the most important outstanding problems related to NS accretion. We show how these open questions regarding accreting NSs in both low-mass and high-mass X-ray binary systems can be addressed with the High-Energy X-ray Probe (HEX-P) via simulated data. In particular, with the broad X-ray passband and improved sensitivity afforded by a low X-ray background, HEX-P will be able to 1) distinguish between competing continuum emission models; 2) provide tighter upper limits on NS radii via reflection modeling techniques that are independent and complementary to other existing methods; 3) constrain magnetic field geometry, plasma parameters, and accretion column emission patterns by characterizing fundamental and harmonic cyclotron lines and exploring their behavior with pulse phase; 4) directly measure the surface magnetic field strength of highly magnetized NSs at the lowest accretion luminosities; as well as 5) detect cyclotron line features in extragalactic sources and probe their dependence on luminosity in the super-Eddington regime in order to distinguish between geometrical evolution and accretion-induced decay of the magnetic field. In these ways HEX-P will provide an essential new tool for exploring the physics of NSs, their magnetic fields, and the physics of extreme accretion.


INTRODUCTION 1.Accreting Neutron Stars
Accretion is a ubiquitous process in the Universe, from the formation of stars and planets to supermassive black holes at the center of galaxies.X-ray binary systems are composed of a compact object (CO), either a black hole (BH) or neutron star (NS), that accretes from a stellar companion.These systems can be further categorized based on the mass of the stellar companion: low-mass and high-mass (see Longair 2011 for a more detailed review).Figure 1 shows a simple schematic for accretion onto NSs in several types of binary systems.Low-mass X-ray binaries (LMXBs) have a ≲ 1 M ⊙ companion that transfers matter via Roche-lobe overflow to form an accretion disk around the compact object.High-mass X-ray binaries (HMXBs) have a more massive (≳ 5 M ⊙ ), early-type stellar companion and accretion onto the CO can occur in various ways, including periodic Roche lobe overflow (such as when the CO passes through the decretion disk of a Be type star), or capture of material ejected through dense stellar winds.We note that there is a rare class of LMXBs that accrete from the slow winds of a late-type stellar companion and exhibit properties akin to HMXBs (see Bozzo et al. 2022 and references therein) 1 .
NSs are the densest known objects with a surface in the Universe.The matter inside a NS exists in an ultra-dense, cold state that cannot be replicated in terrestrial laboratories.The only way to discern how matter behaves under these conditions is by determining the equation of state (EoS).The EoS sets the mass and radius of the NS through the Tolman-Oppenheimer-Volkoff equations (Tolman, 1934(Tolman, , 1939;;Oppenheimer and Volkoff, 1939) and astronomical measurements of NS masses and radii are therefore crucial for determining which theoretical EoS models are viable (see, e.g., Lattimer 2011).In particular, accretion onto NSs probes the properties and behavior of matter in the presence of a strong magnetic field (B ∼ 10 8−9 G for LMXBs and B ∼ 10 12 G for HMXBs; see Caballero and Wilms 2012).
In order to fully encapsulate the accretion emission from these systems, a broad X-ray passband is necessary (see §2 for LMXBs and §3 for HMXBs).Currently, NuSTAR (Harrison et al., 2013) is the only focusing hard X-ray telescope with a passband from 3-80 keV.Energies below 3 keV have to be supplemented with other X-ray telescopes, such as NICER (Gendreau et al., 2012), Swift (Gehrels et al., 2004), or XMM-Newton (Jansen et al., 2001), to construct a broad energy passband.While observations can be scheduled to occur during the same observing period, often the data are not strictly simultaneous due to the different orbits of the missions resulting in various degrees of overlap.Since NS binary systems are highly variable, sometimes on time-scales comparable to or shorter than the typical exposure itself, simultaneous observations over a broad X-ray passband are invaluable for studying these systems.
Figure 1.Simple schematic diagrams of several types of X-ray binary systems: (a) a LMXB system where the NS accretes from a stellar companion via Roche-lobe overflow, (b) a NS in a HMXB that accretes from material captured from stellar winds that are launched from the companion, and (c) a NS accreting from the decretion disk of a Be type stellar companion.

The High Energy X-ray Probe: HEX-P
The High Energy X-ray Probe (HEX-P; Madsen et al. 2023, in prep.) is a probe-class mission concept that offers sensitive broad-band X-ray coverage (0.2-80 keV) with exceptional spectral, timing and angular capabilities.It features two high-energy telescopes (HETs) that focus hard X-rays and one low-energy telescope (LET) that focuses lower-energy X-rays.
The LET consists of a segmented mirror assembly coated with Ir on monocrystalline silicon that achieves an angular resolution of 3.5 ′′ , and a low-energy DEPFET detector, of the same type as the Wide Field Imager (WFI; Meidinger et al. 2020) onboard Athena (Nandra et al., 2013).It has 512 × 512 pixels that cover a field of view of 11.3 ′ × 11.3 ′ .The LET has an effective passband of 0.2-25 keV, and a full frame readout time of 2 ms, which can be operated in a 128 and 64 channel window mode for higher count-rates to mitigate pile-up and faster readout.Pile-up effects remain below an acceptable limit of ∼1% for fluxes up to ∼100 mCrab in the smallest window configuration.Excising the core of the PSF, a common practice in X-ray astronomy, will allow for observations of brighter sources, with a typical loss of up to ∼60% of the total photon counts.
The HET consists of two co-aligned telescopes and detector modules.The optics are made of Nielectroformed full shell mirror substrates, leveraging the heritage of XMM-Newton, and coated with Pt/C and W/Si multilayers for an effective passband of 2-80 keV.The high-energy detectors are of the same type as flown on NuSTAR, and they consist of 16 CZT sensors per focal plane, tiled 4 × 4, for a total of 128 × 128 pixel spanning a field of view of 13.4 ′ × 13.4 ′ .This paper highlights interesting existing open questions about NSs and accretion in strong magnetic fields, and demonstrates HEX-P's unique ability to address these with the current best estimate (CBE) mission capabilities.All simulations presented here were produced with a set of response files that represents the observatory performance based on CBEs as of Spring 2023 (see Madsen et al. 2023, in prep.).The effective area is derived from raytracing calculations for the mirror design including obscuration by all known structures.The detector responses are based on simulations performed by the respective hardware groups, with an optical blocking filter for the LET and a Be window and thermal insulation for the HET.The LET background was derived from a GEANT4 simulation (Eraerds et al., 2021) of the WFI instrument, and the HET background was derived from a GEANT4 simulation of the NuSTAR instrument.Both assume HEX-P is in an L1 orbit.The broad X-ray passband and superior sensitivity will provide a unique opportunity to study accretion onto NSs across a wide range of energies, luminosity, and dynamical regimes.

Background
Persistently accreting NS LMXBs are divided into two types based upon characteristic shapes that are traced out in hardness-intensity diagrams and color-color diagrams (Hasinger and van der Klis, 1989): 'Z' sources and 'atoll' sources.Z sources trace out a Z-shaped pattern with three distinct branches: the horizontal, normal, and flaring branch (HB, NB, and FB, respectively), with the FB being the softest spectral state.They can be further divided into two subgroups, Sco-like and Cyg-like, based upon how much time they spend in the different branches.Sco-like sources spend little to no time in the HB and extended periods of time in the FB (> 12 hours), whereas Cyg-like sources have a strong HB and spend little time in the FB (Kuulkers et al., 1997).Atoll sources are either observed in the soft 'banana' state or hard 'island' state.A difference in mass accretion rate is thought to be the primary driver between the behavioral differences between atoll and Z sources.Atoll sources probe a lower range in Eddington luminosity (∼0.01-0.5 L Edd ) in comparison to the near-Eddington Z sources (∼0.5-1.0L Edd : van der Klis 2005).Further evidence that the average mass accretion rate is responsible for the observed difference between the two classes comes from observing transient NS LMXBs that have transitioned between the two classes during outburst (e.g., XTE J1701−462: Homan et al. 2010).
NS LMXBs occupy a number of spectral states which vary considerably in terms of the models and spectral parameters needed to describe the continuum emission (Barret, 2001).In the very hard state, the spectrum is dominated by Comptonization that can be modeled with an absorbed power-law component with a photon index Γ ∼ 1 (Ludlam et al., 2016;Parikh et al., 2017) or two thermal Comptonization components assuming two distinct populations of seed photons from different plasma temperatures (Fiocchi et al., 2019).In the hard state, the spectrum is dominated by a hard Comptonization component with Γ = 1.5-2.0 and a soft thermal component arising from a single temperature blackbody component or multi-color disk blackbody with a temperature ≲ 1 keV (Barret et al., 2000;Church and Balucińska-Church, 2001).In the soft state, the spectrum becomes thermally dominated with weakly Comptonized emission.Model choices for the thermal and Comptonization components in the soft state vary in the literature leading to two classical descriptions.The "Eastern" model, after Mitsuda et al. (1989), uses a multi-color disk blackbody in combination with a Comptonized blackbody component, while the "Western" model, after White et al. (1988), uses a single-temperature blackbody and a Comptonized disk component.Lin et al. (2007) devised a hybrid model for hard and soft state spectra based upon RXTE observations of two transient atoll systems that resulted in a coherent picture of the spectral evolution (e.g., the thermal components follow the expected L x ∝ T 4 relation).In this model, the soft state assumes two thermal components (i.e., a single-temperature blackbody and a disk blackbody) and weak Comptonization accounts for the power-law component.This hybrid double thermal model has been used in many NS LMXBs studies (e.g., Cackett et al. 2008Cackett et al. , 2009;;Lin et al. 2010), though not exclusively.For instance, a recent study using thermal Comptonization from a blackbody and a multi-color accretion disk blackbody (akin to the "Eastern" model) to model the X-ray spectra of the atoll 4U 1820-30 found good agreement with the observed jet variability in this system (Marino et al., 2023).In the absence of multi-wavelength data to support a choice of continuum model, it is difficult to ascertain which prescription of the spectra is appropriate.Due to the soft spectral shape in these states, the source spectrum typically becomes background dominated above 30 keV even when observed with NuSTAR.A broad X-ray passband with large effective collecting area and low X-ray background that can observe a large number of sources at various luminosity levels and spectral states is needed to further our understanding of these sources.
The accretion disks in these systems can be externally illuminated by hot electrons in the corona (Sunyaev et al., 1991) or from the thermal emission of the NS or boundary layer region (Popham and Sunyaev, 2001).We note that the exact geometry of the corona is unclear (see Degenaar et al. 2018 for some possible geometries), but X-ray polarization measurements with IXPE are beginning to shed light on coronal orientation and presence of boundary layer regions in some systems (e.g., Jayasurya et al. 2023;Ursini et al. 2023;Cocchi et al. 2023;Farinelli et al. 2023); building up our understanding of the accretion geometry when coupled with X-ray energy spectral studies.Regardless of the hard X-ray source, the disk reprocesses the hard X-ray emission and re-emits discrete emission lines superimposed upon a reprocessed continuum known as the 'reflection' spectrum (Ross and Fabian, 2005;García and Kallman, 2010).The spectrum is then relativistically blurred due to the extreme environment close to the NS, where accreting material reaches relativistic velocities as it falls into the NS's deep gravitational well (Fabian et al., 1989;Dauser et al., 2013).Studying reflection in NS LMXBs allows for properties of NSs and the disk itself to be measured, such as the NS magnetic field strength (Ibragimov and Poutanen, 2009;Cackett et al., 2009;Ludlam et al., 2017b), extent of the boundary layer (Popham and Sunyaev, 2001;King et al., 2016;Ludlam et al., 2021), and NS radius (Cackett et al., 2008;Ludlam et al., 2017aLudlam et al., , 2022)).Of particular interest is both the inner disk radius, R in , and dimensionless spin parameter, a = cJ/GM 2 (which is the mass normalized total angular momentum J of the CO); the latter is typically fixed in current studies due to limitations in data quality.The spin parameter has important consequences for both accreting NSs and BHs (see Connors et al. 2023, in prep., andPiotrowska et al. 2023, in prep., for HEX-P science with BH X-ray binaries and supermassive BHs, respectively).The spin sets the location of the innermost stable circular orbit (ISCO) where a higher spin corresponds to a smaller ISCO (Bardeen et al. 1972).Consequently, the position of the inner disk radius for higher spin values cannot be replicated by lower spin if the data is of sufficient quality to recover R in accurately.The majority of NSs in LMXBs have spin a ≲ 0.3 (Galloway et al., 2008;Miller et al., 2011) and the difference in the location of the ISCO decreases from 6 gravitational radii (R g = GM/c2 ) at a = 0 to 4.98 R g at a = 0.3 (see Fig. 2).Therefore, targeting an accreting NS with a higher spin, as indicated through a high measured spin frequency, can decrease the upper limit on NS radii obtained from reflection modeling by roughly 2-3 kilometers.
In order to capture the entire reflection spectrum (i.e., the low-energy O VIII and iron (Fe) L emission lines near 1 keV (Madej et al., 2014;Ludlam et al., 2018Ludlam et al., , 2021)), the Fe K emission lines at 6.4-6.97 keV, and the Compton backscattering hump above 15 keV), determine the appropriate continuum model, and measure the absorption column along the line of sight, a broad X-ray passband is necessary.Hence, the broad X-ray sensitivity provided by the LET and HETs on HEX-P is crucial for studying accreting NS LMXBs.

Simulated Science Cases for LMXBs
All simulations in this section were conducted in xspec (Arnaud, 1996) via the 'fakeit' command, which draws photons from a randomized seed distribution, with V07 of the HEX-P response files selecting an 80% PSF correction, assuming the data would be extracted from a 15 arcsec region for the HET and 8 arcsec region for the LET, as well as the anticipated background for the telescope at L1.The simulated spectra were grouped via grppha to have a minimum of 25 counts per bin to allow for the use of χ 2 statistics.The flux for each of the following simulations can be found in Table S1 in the Supplemental Materials.

Distinguishing between continuum models
To demonstrate HEX-P's ability to distinguish between different continuum models, we base our simulations on the models used to fit simultaneous NICER and NuSTAR observations of the accreting atoll 4U 1735−44 published in Ludlam et al. (2020).Model 1 is based on the hybrid double thermal continuum model of Lin et al. (2007) while Model 2 is based on the "Eastern" model of Mitsuda et al. (1989).Both models provide an adequate description of the data in the 0.3-30 keV energy band (the NuSTAR observations were background dominated above 30 keV) when neglecting the reflection spectrum, but data at higher energies would distinguish between these two continuum models.We take the continuum model and parameter values for Model 1 and Model 2 from Table 1 of Ludlam et al. (2020) as input for the HEX-P simulation; the exact values can be found in Table S2 of the Supplementary Materials.The overall models 2 are as follows: (a) Model 1: tbabs*(diskbb+bbody+pow) and (b) Model 2: tbabs*(diskbb+nthcomp).Note that we fix the absorption column index to a value between the ones used in Model 1 and Model 2, N H = 4 × 10 21 cm −2 , to focus on the difference in the spectral shape due to the model component choices rather than a difference in absorption column between the models.We chose an exposure time of 20 ks, roughly equivalent to the exposure time of the NuSTAR observation in the aforementioned study.
The simulated HEX-P spectra remain above the background in the HET bands up to 80 keV, though they dive into the LET background below 1 keV.This could be remedied by increasing the exposure time, though, notably, the spectrum above 30 keV is where the models diverge.The reduced χ 2 is ≤ 1.03 when fitting each simulation with its own input model.Attempting to fit each simulated spectrum with the opposing model description provides an inadequate fit with ∆χ 2 > 1000 for the same number of degrees of freedom (d.o.f.). Figure 3 shows the HEX-P spectrum for the Model 2 input of a thermal disk with a Comptonized blackbody from a boundary layer region.For direct comparison, we show a simulated spectrum using Model 1 (the double thermal continuum model with weak Comptonization) with parameter values from fitting the simulated spectrum of Model 2 below 30 keV, since these models perform equally blue).The models diverge above 30 keV where currently available missions quickly become background dominated for the same exposure time.To emphasize this point, the same S/N near 30 keV could be achieved with NuSTAR in an exposure time of 96.6 ks.The broad X-ray passband and improved sensitivity of HEX-P will differentiate the models.
well within this energy regime.This highlights that though these models are nearly identical below 30 keV, they diverge at higher energies inaccessible by current operating missions.The weakly Comptonized emission above 30 keV cannot be replicated by Model 2, and demonstrates the potential for a broadband X-ray mission like HEX-P to discern between competing continuum models that would otherwise equally describe the data below 30 keV.

Neutron star radius constraints from relativistic reflection modeling
To demonstrate that HEX-P would improve the radius constraints obtained by reflection modeling of NSs, we base our simulations on the results of a recent investigation of Cygnus X-2 that utilized simultaneous data from NICER and NuSTAR (Ludlam et al., 2022).This source is of particular interest since it has a dynamical NS mass measurement of M NS = 1.71 ± 0.21 M ⊙ (Casares et al., 2010).The source was analyzed in the NB, HB, and the vertex between those branches.Fixing the spin at a = 0, Ludlam et al. (2022) found that the inner disk radius remained close to R ISCO .In the current literature it is customary to find studies in which the spin of the NS has been fixed while fitting models to the data given the highly degenerate nature of R in and a.We set up our simulations using the Cygnus X-2 data in the NB modeled with the latest public release of the self-consistent reflection model tailored for thermal emission illuminating the accretion disk, relxillNS (v2.2:García et al. 2022).The input parameters for the simulation3 are provided in Table 1.We leave R in at 1 R ISCO for all simulations while varying the input spin since the disk remained consistent with this value while the source was in the various branches.We chose three spin values as test cases: non-rotating (a = 0), the highest value expected for a NS LMXB (a = 0.3), and a = 0.17 which is an Table 1.Input model parameters for simulating reflection model data with HEX-P.Values are taken from Table 3: RNS1 of Ludlam et al. (2022) for Cygnus X-2 in the NB and updated using the public version of relxillNS (v2.2).Additionally, we show one case of the recovered model parameters from a simulated HEX-P spectrum after performing a Markov Chain Monte Carlo (MCMC) analysis with a burn-in of 10 6 and chain length of 5 × 10 5 to emulate a standard analysis of the data.Errors are reported at the 1σ confidence level, though the errors for R in and a are drawn from the bivariate distribution between the two due to their correlated nature.
Model A Fe 5.9 4.9 +0.7 approximation based on the spin frequency of the source (Mondal et al., 2018).We use an exposure time of 100 ks and investigate how well R in and spin can be recovered independently.
Given the degenerate nature of the parameters of interest and that each spectrum simulated is created via a randomized photon generation, the results obtained can vary with each iteration.We therefore simulate 10 4 spectra per spin value to determine the likelihood of constraining R in and a.The LET data were modeled in the 0.3-15 keV energy band while the HET data were modeled in the 2-80 keV band.We impose an upper limit on the spin (a ≤ 0.7) when fitting the simulated spectra so as to not surpass the rotational break-up limit of a NS (Lattimer and Prakash, 2004).We fit each spectrum and then perform error scans to ensure that each iteration reached a minimum goodness of fit of χ 2 /d.o.f.≤ 1.1 prior to obtaining values for inner disk radius and spin from each simulated spectrum, thus building up the posterior distribution shown in Figure 4.
Errors are drawn from the bivariate posterior distribution due to their correlated nature and reported at the 1σ confidence level.We find that for an input of a = 0, we recover a = −0.03+0.15  −0.17 and R in = 1.02 +0.07 −0.02 R ISCO .For an input of a = 0.17, we recover a = 0.17 +0.09 −0.16 and R in = 1.02 +0.05 −0.02 R ISCO .Lastly, for an input of a = 0.30, we recover a = 0.27 ± 0.13 and R in = 1.02 +0.05 −0.02 R ISCO .The ability to recover the input parameters improves at higher spin as the relativistic effects become stronger.Recall that the position of the inner disk radius for higher spin values cannot be replicated by lower spin, thus if the data is of sufficient quality to recover R in accurately this will become apparent.For example, an input of R in = 1 R ISCO for a spin of a = 0 corresponds to 6 R g .This value can be mimicked by a higher spin and R in beyond 1 R ISCO (e.g., 6 R g is equal to 1.2 R ISCO for a = 0.3).However, an input of R in = 1 R ISCO for a = 0.3 corresponds to 4.98 R g and cannot be replicated by a lower spin so long as the upper limit recovered when fitting the data returns a value of R in < 1.2 R ISCO .We note that Cygnus X-2 is a bright LMXB (∼0.7 Crab) and would thus need to be observed in a configuration that would reduce pile-up effects in the LET.However, we obtain consistent results at the 1σ level when performing the same exercise with just the HETs, which do not suffer from pile-up, though the 3σ level can no longer rule out a non-rotating NS as is the case when the LET is included.Hence the combined observing power of the HETs and LET provide an improvement over the current capabilities of existing missions.
Figure 5 shows the improvement provided by HEX-P for reflection modeling of spinning NS LMXBs in comparison to simultaneous NICER and NuSTAR results with fixed spin.Additionally, the current best constraints from gravitational wave events of binary NS mergers and pulsar light curve modeling demonstrate how the various methods can work in concert to narrow down the allowable region for the EoS.Each method has its own underlying systematic uncertainties, so the firmest EoS constraints will require multiple measurement approaches.Utilizing the dynamical NS mass estimate of Cygnus X-2 (M NS = 1.71 ± 0.21 M ⊙ : Casares et al. 2010), Ludlam et al. (2022) reported an upper limit on the radius of R NS = 19.5 km for M NS = 1.92 M ⊙ and R NS = 15.3 km for M NS = 1.5 M ⊙ .The lower limit on spin and upper limit on R in return a conservative upper limit on the NS radius.From the HEX-P simulations we calculate an upper limit of R NS = 16.7 km for M NS = 1.92 M ⊙ and R NS = 13.2 km for M NS = 1.5 M ⊙ for an input of a = 0.3.This demonstrates the power of HEX-P to study rotating NSs with reflection studies.For comparison, we conducted the same simulation with NuSTAR response files.For the highest spin value of a = 0.3, NuSTAR recovers a = 0.67 +0.03 −0.54 and R in = 1.02 +0.16 −0.02 R ISCO .The conservative upper limit on the NS radius from the location of the inner disk radius is larger than 6 R g and therefore does not provide an improvement over current reflection studies that assume the spin is fixed at a = 0. Performing this test with the combination of NICER and NuSTAR marginally improves the recovered spin value to a = 0.68 +0.02 −0.50 , but with a larger upper limit on R in = 1.02 +0.18 −0.02 R ISCO , which is still > 1 km larger than the constraints obtained from HEX-P (i.e., R NS = 18.5 km for M NS = 1.92 M ⊙ and R NS = 14.4 km for M NS = 1.5 M ⊙ ).
We demonstrate that these results do not depend on drawing from a posterior distribution of simulated spectra by conducting Markov Chain Monte Carlo (MCMC) analysis on one of the simulated HEX-P spectra with an input spin of a = 0.3.The results of this test are shown in Table 1 with the inclination fixed at the median value from optical and X-ray studies (Orosz and Kuulkers, 1999;Ludlam et al., 2022).The lower limit on spin of a = 0.1 and upper limit on R in = 1.12 R ISCO still provides an improvement over current NICER and NuSTAR analyses, which fix a = 0.The MCMC analysis finds an upper limit of R NS = 18 km for M NS = 1.92 M ⊙ and R NS = 14.1 km for M NS = 1.5 M ⊙ .
For completeness we note that the Kerr metric is used to describe the space-time close to the NS in the relativistic reflection model.As the angular momentum (i.e., spin) increases, the NS can become oblate, thus causing deviations in space-time from the Kerr metric due to an induced quadrupole moment.The exact induced deviation depends upon the EoS of the NS, but this has an effect < 10% (Sibgatullin and Sunyaev, 1998) in the anticipated range of spin parameters for NS LMXBs (a ≲ 0.3).Given that our 1σ errors for R in are larger than this deviation from the Kerr metric, the HEX-P spectra would not be sensitive enough to measure this effect.However, our conservative upper limits on NS radius are still at larger radii than the deviation in the R ISCO and hence are not in conflict even if we are insensitive to an induced quadrupole moment.Therefore, reflection studies of accreting NS LMXBs with HEX-P would provide improved upper limits on NS radii in comparison to current studies.

Background
NSs in HMXBs are typically highly magnetized, with a dipolar magnetic field strength of ∼10 12 G.The pressure of the TeraGauss magnetic field disrupts the inflow of the accreting matter, channelling it onto the magnetic poles (Davies and Pringle, 1981).There, the kinetic energy of the matter is released as radiation in a highly anisotropic manner.Since the NS rotates, the sources are visible as X-ray pulsars.Studying this emission is crucial to understanding the physics of plasma accretion and the interaction of radiation with high magnetic fields that are orders of magnitude stronger than those achievable in laboratories on Earth.Due to their extreme surface magnetic field strength, quantum effects take place at the site of emission.In particular, the electron motion in the direction perpendicular to the magnetic field lines becomes quantized.This affects the electron scattering cross section, leading to a resonance at an energy that is proportional to the magnetic field strength (Meszaros, 1992).As a result, the spectra of X-ray pulsars (XRPs) can contain absorption line-like features called cyclotron resonant scattering features (CRSFs), or cyclotron lines, at an energy of where z g is the gravitational redshift (typically 0.3 for standard NS mass and radius) and B 12 is the NS (polar) magnetic field strength in units of 10 12 G, while n is the number of Landau levels involved4 .Cyclotron lines represent the most direct way to probe the NS magnetic field on or close to the NS surface.
Up to now, cyclotron lines have been confirmed in the spectra of about 40 sources (Staubert et al., 2019).
Figure 6 shows a schematic description of the formation of cyclotron harmonics.
When observed at typical outburst luminosity of ≳ 10 36 erg s −1 , the spectra of XRPs have been described with phenomenological models (a power-law with high-energy cutoff) modified by interstellar absorption, an Fe Kα emission line around 6.4 keV, low-temperature blackbody components and, at times, a component called the 10 keV feature (Mushtukov and Tsygankov, 2022, and references therein).However, at lower luminosity accretion regimes, the spectrum often exhibits a transition to a double-hump shape (Tsygankov et al., 2019a) which has recently been explained in terms of a low-energy thermal component peaking at ∼5 keV and a high-energy Comptonized component peaking at ∼35 keV representing the broadened red wing of the cyclotron line (Sokolova-Lapa et al., 2021;Mushtukov et al., 2021).A dip between the two components was sometimes interpreted as a cyclotron line (see, e.g., Doroshenko et al., 2012, for the case of X Per), but is likely a continuum feature, as was shown for sources with a known fundamental cyclotron line at higher energies after the transition to a low-luminosity state (see, e.g., Tsygankov et al., 2019a).These low luminosity states of accretion characterized by two-hump spectra are generally associated with braking of the accretion flow in the NS atmosphere (Mushtukov et al. 2021, Sokolova-Lapa et al. 2021, although a model assuming an extended collisionless shock instead was recently proposed by Becker and Wolff 2022).Rapid deceleration of the plasma in the atmosphere proceeds mainly by Coulomb collisions and leads to the formation of an overheated outermost layer with electron temperatures of ∼30 keV.Such high temperatures significantly enhance Comptonization and together with resonant redistribution lead to the formation of the high-energy excess in the spectra around the cyclotron line.When the mass-accretion rate onto the poles is increased, the pressure of the emitted radiation starts affecting the dynamics of the accretion flow.The radiation becomes capable of decelerating the flow above the surface, reducing direct heating of the atmosphere and giving rise to a radiation-dominated shock (Basko and Sunyaev, 1976).The onset of this regime is typically associated with a critical luminosity (L crit ∼10 37 erg s −1 : Basko and Sunyaev, 1976;Becker et al., 2012;Mushtukov et al., 2015) and growth of the accretion column as an emitting structure.Modeling emission from accretion columns is a very challenging problem due to its dynamic nature, multi-dimensional geometry, and the necessity of including general relativistic effects.On the other hand, studying the low-luminosity regime with emission from the heated atmosphere of the polar cap allows easier access to fundamental parameters of accreting NSs (see § 3.2.3 for more details).
Cyclotron lines exhibit variations with accretion luminosity, pulse phase, and over secular time scales.Variation of the cyclotron line centroid energy E cyc with X-ray luminosity L X can be exploited to gain insight into the accretion regime.For example, low-luminosity XRPs often show a positive correlation between E cyc and L X , while high-luminosity XRPs show a negative correlation (see, e.g., Klochkov et al. 2011).The two trends are assumed to be separated by the critical luminosity and thus are related to the transition between the different accretion regimes described above.On the other hand, studying the pulse-phase dependence of the E cyc (as well as the line depth and width) can shed light upon the geometrical configuration of the NS (Nishimura, 2011).Moreover, some sources show evidence of a pulse phase-transient cyclotron line; i.e., a CRSF that only appears at certain pulse phases (Klochkov et al., 2008;Kong et al., 2022).Finally, secular variation of the cyclotron line has been attributed to a geometric reconfiguration of the polar field or to a magnetic field burial due to continued accretion (Staubert et al., 2020;Bala et al., 2020).
One of the biggest challenges for the detection and characterization of cyclotron lines is determining a robust broadband continuum model.The centroid energy of the line as well as its other parameters can show significantly different best-fit values when modelled with different continua.Due to this model-dependence, the very presence of a cyclotron line has been questioned at times (Di Salvo et al., 1998;Doroshenko et al., 2012Doroshenko et al., , 2020) ) as has its luminosity dependence (Müller et al., 2013).A similar argument also holds for the continuum best-fit values.The "true" continuum model can only be inferred if data cover a sufficiently wide energy passband to constrain absorption affecting the softer X-ray energies and the various spectral components at harder X-ray energies (Sokolova-Lapa et al., 2021;Malacaria et al., 2023a).Moreover, properly constraining the broadband continuum is necessary when multiple cyclotron line harmonics are present, as in the cases of 4U 0115+63 (Heindl et al., 2004) and V 0332+53 (Pottschmidt et al., 2005).In addition, observing the broadband continuum from accreting XRPs can help constraining physical parameters of the NSs thanks to the recent development of physically motivated spectral models (e.g., Farinelli et al. 2016;Becker and Wolff 2022, and references therein).This is especially important at low-luminosity, where the accretion flow free-falls on the NS surface and our understanding of the physical mechanisms is more uncertain concerning, e.g., how the flow transitions from radiation-dominated to gas-dominated, or the detailed production of seed photons from cyclotron, bremsstrahlung, and blackbody mechanisms at the site of impact.HEX-P will bring crucial contributions also in this field, as it will be able to probe broadband spectral emission from intrinsically dim sources whose required observing exposure with current facilities would be prohibitive.
Last but not least measuring phase-dependent spectral components with the highest sensitivity is required in order to develop consistent physical models for the emission process that include X-ray polarization properties.Recent IXPE observations of accreting pulsars like Cen X-3, Vela X-1, EXO 2030+375 and others provided new constraints as they showed low, energy-dependent polarization degrees of <10% (Tsygankov et al., 2022;Forsblom et al., 2023;Malacaria et al., 2023b).Radiative transfer models of highly magnetized plasmas predict considerably higher polarization degrees of 50-80% (Meszaros et al., 1988;Caiazzo and Heyl, 2021;Sokolova-Lapa et al., 2023).Initial candidates for explaining the difference highlight the potential importance of the temperature structure of the plasma in the accretion column, of angle-dependent beaming, and of the propagation of X-rays in the magnetosphere for lowering the polarization degree.Furthermore it has been demonstrated that the X-ray continuum and cyclotron line profile can be expected to directly show subtle, complex imprints of polarization effects (Sokolova-Lapa et al., 2023).This is an intriguing possible reason for the occasionally observed and notoriously difficult to constrain "10 keV feature" mentioned above, as well as for the fact that some accreting pulsars do not show cyclotron lines.
Advancing our understanding of HMXBs in general and cyclotron line sources in particular therefore requires broadband spectra with high sensitivity throughout the relevant energy passband, with a limited background as provided by focusing X-ray facilities, and medium spectral energy resolution (Wolff et al., 2019).Thanks to its leap in observational capabilities, HEX-P will be capable of overcoming all abovementioned challenges and push our understanding of accretion onto XRPs forward.In the next section, we simulate a few exemplary cases highlighting the gains enabled by HEX-P with respect to currently available X-ray facilities for which the investigated cyclotron lines centroid energy, accretion regime, or the necessary exposure time prevent a comprehensive study of the physical mechanisms at work.

Simulated Science Cases for HMXBs
The following simulations show HEX-P's potential to tackle different cyclotron line science cases using several example sources.Each source has been chosen to represent a specific science case, either because it allows us to observe a certain accretion regime and luminosity range or due to its specific cyclotron line properties, such as the existence of n > 1 harmonics.We begin with the prototypical persistent cyclotron line sources Cen X-3 and Vela X-1, which show moderate to high fluxes and luminosities, and allow for exquisitely detailed and sensitive parameter constraints with HEX-P.Then we focus on transient accreting pulsars for which HEX-P will provide unprecedented access to low fluxes, allowing study of extreme luminosity regimes, e.g., for GX 304-1 in quiescence or for super-Eddington outbursts of extragalactic sources like SMC X-2.Similar to § 2.2, all simulations were performed via the 'fakeit' command in xspec (Arnaud, 1996) or ISIS (Houck and Denicola, 2000), employing version v07 of the HEX-P response files, selecting an 80% PSF correction, assuming a 15 ′′ extraction region for the HETs and 8 ′′ for the LET, as well as taking the expected background at L1 into account.The flux for each of the following simulations can be found in Table S1 in the Supplemental Materials.

Constraining magnetic field geometry from cyclotron line harmonics
As discussed earlier, cyclotron lines result from transitions between quantized Landau energy levels of electrons in the presence of a magnetic field (Meszaros, 1992).Measuring the energy, width, and strength of cyclotron lines in HMXBs is crucial for understanding the underlying physics of these systems.The energy of cyclotron lines provides insights into the line-forming region in the accretion column.The width of cyclotron lines provides valuable information about the geometry, temperature distribution, and plasma conditions within the accretion column (see e.g., Becker et al. 2012;Staubert et al. 2014 and references therein).Broadening of cyclotron lines can be influenced by factors such as electron thermal motion, turbulence, and relativistic effects.In the last few decades, comprehensive analyses of HMXBs have directly unveiled many important properties of these systems (see, e.g., Staubert et al. 2019;Pradhan et al. 2021 and references therein).In addition to the fundamental, higher harmonics of the cyclotron line are also sometimes present in the X-ray spectrum.They arise due to higher-order interactions between X-rays and the Landau levels.By detecting and analyzing fundamental cyclotron lines and their higher harmonics, one can determine the magnetic field strength in different regions of the accretion column.Furthermore, the strength or intensity of fundamental cyclotron lines and their higher harmonics relative to the fundamental line provide insight into the scattering efficiency and the fraction of scattered photons, improving our understanding of the emission processes and properties of the compact object (e.g., see Alexander and Meszaros 1991;Schwarm et al. 2017 and references therein).
Given the limited passband and / or higher background of instruments like RXTE and Suzaku, there have been few studies to date of cyclotron lines higher harmonics in HMXBs (see Table 3 of Pradhan et al., 2021).More recent instruments like NuSTAR and Insight/HXMT provide a significant improvement, but as we look into the future, these studies will be revolutionized by a mission with a broad passband coupled with very low background like HEX-P.In order to test this, we simulated a 20 ks HEX-P spectrum of Cen X-3 in xspec.Due to its bright nature and that Cen X-3 is a strong hard X-ray emitter (see Table S1 in the Supplemental Material), we concentrate on the HETs only which are devoid of pile-up although we have verified that the results below are consistent when including the LET.Cen X-3, an eclipsing HMXB, consists of an O star and a pulsar with a rotational period of ∼4.8 s (Chodil et al., 1967;Giacconi et al., 1971;Schreier et al., 1972).The X-ray spectrum of this source exhibits multiple emission lines (Fe, Si, Mg, Ne: Iaria et al. 2005;Tugay and Vasylenko 2009;Naik and Paul 2012;Aftab et al. 2019;Sanjurjo-Ferrín et al. 2021), and the X-ray emission has been detected beyond 70 keV.The fundamental cyclotron line occurs near 29 keV (Tomar et al., 2021) with a possible n 2 harmonic at ∼47 keV (Yang et al., 2023) as seen from recent HXMT results.With the passband of HEX-P extending up to 80 keV and significantly less background, we investigated the possibility of detecting an n 3 harmonic in the source.
The baseline model for our simulation is NPEX (Equation 2; Mihara 1995;Makishima et al. 1999) and the values of model parameters are taken from Table 1 (ObsID P010131101602) of Yang et al. (2023).The show the residuals obtained when we set the strength of fundamental, n 2 harmonic, and n 3 harmonic to zero, respectively.We re-binned the data visually and provide arrows for the cyclotron line components to aid with clarity.For comparison, we provide a simulated NuSTAR spectrum for using the same model and exposure time which is shown in panel (e).The NuSTAR data are strongly background dominated at these flux levels, which makes constraining the higher energy harmonics difficult even at much longer exposure times (see text for more details).functional form of NPEX is given by where Γ 1 and Γ 2 are the negative and positive power-law indices, respectively.Γ 2 , which approximates a Wien hump, was fixed at 2.0, Γ 1 = 0.8, the folding energy E fold =6.9 keV, and the absorption column N H = 2.06× 10 22 cm −2 .
The CRSFs are modeled with a multiplicative model of a line with a Gaussian optical depth profile with energy E i , strength d i , and width σ i .The fundamental line, i = 1, is at ∼29 keV with d 1 ∼ 1.4 and σ 1 = 7.6 keV.The n 2 harmonic line, i = 2, is at ∼47 keV, with d 2 ∼ 2.3 and σ 2 = 9.7 keV.The 2-75 keV flux is ∼1× 10 −8 erg cm −2 s −1 .In order to showcase the unparalleled accuracy of HEX-P in measuring magnetic fields correspondent to cyclotron line energies above 50 keV, we incorporated a representative example into our model: we include an n 3 harmonic line at 70 keV (i = 3), with the width (σ 3 ) and strength (d 3 ) set to match those of the n 2 harmonic line.By doing so, we investigate the ability of HEX-P to effectively detect these parameters in the n 3 harmonic line, from which we can constrain the magnetic field strength and details about the magnetic field structure.
A 20 ks observation HEX-P is able to constrain the energy of the n 3 harmonic well, to within 7% (see Fig. 7).We also applied an F-test to calculate the probability of chance improvement, with and without the n 3 harmonic line, and found that the probability of chance improvement when adding the n 3 harmonic line is below 0.06%, confirming the robustness of the detection.We reiterate here that the clear improvement of HEX-P over its predecessors in terms of broadband coverage and sensitivity makes it possible to investigate higher harmonic features -which have clearly eluded NuSTAR (see e.g., Tomar et al., 2021).
Since cyclotron lines result from the resonant scattering of photons by electrons whose energies are quantized into Landau levels by the strong magnetic field, and the quantized energy levels of the electrons are harmonically spaced, one would naively expect the energy of the n > 1 harmonics to be the integer multiples of the fundamental line.It has however been observed in various X-ray pulsars that this is generally not the case (see Yang et al. 2023;Orlandini et al. 2012 and discussions within).The discrepancy can be understood by assuming a difference in line formation mechanisms or in line-forming regions for the fundamental and higher harmonics.In the former case, the fundamental would primarily arise from resonant scattering, while the higher harmonics involve additional effects like multiple scattering and photon spawning (see e.g., Nishimura 2003;Schönherr et al. 2007).This, however, can not be the only reason because some systems show an an-harmonic spacing larger than predicted by this effect.One way to explain this difference is to take into account that the optical depths of the fundamental and the higher harmonics can be different if they are formed at different heights above the NS.The higher harmonics could be closer to the NS surface and the fundamental line situated at a height with weaker magnetic field strength (see, e.g., Fürst et al. 2018).Another possibility is to consider a displacement of the magnetic dipole, which would also explain the energy difference of the two lines if the lines originate from the different poles of the NS (Rodes-Roca et al., 2009).Therefore, a significant phase dependence of the strength of the fundamental and higher harmonics is expected, which is not possible to probe with the current data sets.Here, too, HEX-P can provide revolutionary new capabilities though (see section 3.2.2).
Finally, note that while the n 2 harmonic of Cen X-3 at 47 keV technically falls within the energy range covered by NuSTAR, meaningful constraints could not be derived with NuSTAR (see panel (e) of Fig. 7) due to the dominant background in that energy range even when the exposure was set to 500 ks.However, with HEX-P, the parameters associated with the cyclotron line can be much better constrained, while also alleviating the degeneracies with the continuum.To emphasize the latter point, we present a contour plot in Figure 8 depicting the relationship between the fundamental cyclotron line energy and the cutoff energy for both HEX-P and NuSTAR data for the same model and exposure of 20 ks for Cen X-3.The plot clearly demonstrates that HEX-P provides significantly more stringent constraints compared to NuSTAR.Such accurate measurements in short exposure times will enable detailed pulse phase-resolved spectra to explore the accretion physics in the accretion column in unprecedented detail, as discussed in the following section.
Note that Cen X-3 is highly variable, the X-ray flux varies by up-to two orders of magnitude even outside eclipses.Therefore, we also fit the actual NuSTAR data in bright flux state (ObsID 30101055002; exposure 21 ks), which has a flux an order of magnitude more than our simulations, and find the harmonics were not detected in the NuSTAR spectrum (see, (Tomar et al., 2021)).In order to make a comparison of HEX-P vs NuSTAR for this bright state of Cen X-3, we simulated an HET spectrum using the NuSTAR model in bright state, while keeping the parameters of harmonics as above.We find that, even for this bright state of Cen Figure 8. Contour plots (68%, 90% and 99% confidence levels) for cutoff energy versus the fundamental cyclotron line energy for Cen X-3 using the same model and exposure for NuSTAR and HEX-P.As evident from the plot, HEX-P will provide better constraints on the continuum and cyclotron line parameters; thereby mitigating degeneracies between the two.Note that the cutoff energies are slightly shifted to aid visualization.
X-3, NuSTAR will be able to obtain the same signal-to-noise as HEX-P only if the exposure is increased by 20 times (i.e., 420 ks).

Constraining accretion column emission geometry through pulse phase-resolved spectroscopy
Vela X-1 is an archetypical HMXB.It has been well studied with all current and past X-ray missions (see Kretschmar et al., 2021, for a recent review).While not being exceptionally luminous (∼10 36 erg s −1 ), its close distance (d = 1.9 kpc) and X-ray eclipses make it an ideal system to study NS magnetic fields and their interaction with the stellar wind of the companion (i.e., the mass donor).The X-ray spectrum of Vela X-1 shows two CRSFs, the fundamental at around 25 keV and the n 2 harmonic around 55 keV (see Fürst et al., 2014, and references therein).These line energies are ideally covered by the HEX-P energy band.However, in Vela X-1, the fundamental line around 25 keV is much weaker and shallower than the n 2 harmonic line at 55 keV.In fact, the 25 keV line is so weak that there has been a long-standing discussion in the literature about its existence, which could only be settled once NuSTAR data were available (Fürst et al., 2014;Kreykenbohm et al., 2002;Maitra and Paul, 2013).The cyclotron lines in Vela X-1 therefore represent a good test case to explore HEX-P's sensitivity to broad and shallow spectral features.
We performed simulations to study how well HEX-P can measure the energy, width, and depth of both CRSFs, in particular as a function of rotational phase of the NS.Variations of CRSF parameters as a Figure 9. Simulated Vela X-1 spectrum for the bin covering phases 0-0.1 (see Fig. 10 for the spectral values in that bin).The n 2 harmonic cyclotron line is clearly visible as a strong dip at the highest energies as indicated by the arrow.function of phase constrain the magnetic field geometry and emission geometry of the accretion column (see, e.g.Iwakiri et al., 2019;Liu et al., 2020).
We base our simulations on the spectral fits of NuSTAR data presented by Diez et al. (2022).In particular, the continuum is modeled by a powerlaw with an exponential cutoff at high energies (using the model "FDcut" in xspec, Equation 3; Tanaka 1986), modified by neutral absorption column at low energies: where Γ is the power-law index, E cut is the cutoff energy, and E fold is the folding energy.In addition to the continuum model used by Diez et al. (2022) we also include soft X-ray emission lines at 6.4 keV (Fe Kα), 2.4 keV (S Kα), 1.8 keV (Si Kα), 1.4 keV (Mg Lα), and 0.9 keV (Ne IX), representing various atomic features in the spectrum (Diez et al., 2023).A simulated spectrum is shown in Figure 9.
As mentioned in §3.2.1, CRSFs are modeled using a multiplicative line model with a Gaussian optical depth profile described by its energy E i , its strength d i , and its width σ i .Here the subscript i denotes either the fundamental line around 25 keV (i = 1) or the n 2 harmonic line around 55 keV (i = 2).The width σ 1 is set to 0.5 × σ 2 (Diez et al., 2022).Instead of relying on the (rather uncertain; see, e.g., Maitra et al. 2018) current knowledge of the phase-resolved behavior of the CRSF parameters, we simulate that the four most relevant parameters (E 1 , E 2 , d 1 , d 2 ) vary sinusoidally with random phase shifts to each other.This approach demonstrates the power of HEX-P to resolve small changes in any of the parameters, even for relatively weak lines.In particular, we assume that the fundamental line E 1 varies by about ±3.5 keV as a function of phase, while the n 2 harmonic line energy E 2 varies by ±5 keV.The strength d 1 and d 2 vary by ±0.3 keV and ±5.0 keV, respectively.
In addition, we also allow the absorption column of the partial absorber to vary, with ±5 × 10 22 cm −2 around the average value of 32.1 × 10 22 cm −2 .While we do not necessarily expect that the absorption column will vary significantly as function of pulse phase, this variability highlights the capabilities of the HEX-P/LET to measure small changes in absorbing columns on time-scales as short as 5 ks.These variations have been observed in time-resolved spectroscopy and allow us to study the physical properties, like density and clump sizes of the accreted medium (Diez et al., 2023).At the same time, the changes in N H do not influence the high energy spectrum where the CRSFs are present.
For our simulations, we assume that LET will be operated in a fast readout mode to avoid pile-up and with negligible deadtime.We simulate spectra with 50 ks exposure time each, and split the data into 10 phase-bins so that each one has an exposure time of 5 ks.We simulate 100 individual spectra for each phase-bin and calculate the standard deviation (SD) of that sample as uncertainties for each parameter.We additionally calculate the 90% and 10% percentile of the 90% uncertainties of each realization as a realistic estimate of the expected uncertainties in a real 50 ks observation.Using multiple realizations for each phase bin allows us to avoid issues with any one particular realization and shows that our uncertainties are properly sampled over the Poisson noise inherent to counting statistics.The results are shown in Figure 10.
All parameters can be very well reconstructed, with the fundamental energy having uncertainties < 0.5 keV for phases where its strength is > 0.4 keV (i.e., phases 0.5-1.0).Even for weaker lines (despite having larger uncertainties of about ±1 keV), the line can be still significantly detected.This is a significant improvement over NuSTAR, for which Fürst et al. (2018) found that the energy was basically unconstrained for line strengths < 0.5 keV.
The n 2 harmonic line energy can be very well constrained (with uncertainties <2 keV) for central energies of the line < 57 keV.Due to the exposure time of only 5 ks per phase-bin and the low cutoff energy of ∼25 keV, the line appears at the very edge of the useful range of the spectrum.Therefore, higher energies become more difficult to constrain, as most of the line is outside the useful passband; however, we still find typical constraints of ±3-4 keV, even for a simulated line energy of 61 keV.
These results represent a significant improvement over previous phase-resolved studies of Vela X-1.For example, using RXTE, Kreykenbohm et al. (1999) found variations of the CRSF energies for both the fundamental and the n 2 harmonic line, but could only use 10 phase bins and still had average uncertainties of 2-3 keV for the fundamental line, and 5-10 keV for the n 2 harmonic.Recently, Liu et al. (2022) published results obtained on Vela X-1 with Insight/HXMT.Using 16 phase-bins for a ∼100 ks exposure, they found average uncertainties comparable to our HEX-P simulations for the n 2 harmonic line, but with significantly larger uncertainties (a few keV) for the fundamental line.We also note that Liu et al. (2022) found the n 2 harmonic line at energies between 40-50 keV, i.e., at much lower energies than simulated here and therefore in a part of the spectrum with a much higher signal-to-noise ratio.
HEX-P observations would therefore be a big step forward in being able to obtain phase-resolved spectroscopy of bright HMXBs, where we can constrain the CRSF parameters, in particular the energy, much better than with existing instruments in shorter exposure times.This reduction in exposure time means that we can either observe more sources in less time or slice the phase-resolved spectra finer to obtain a more detailed look at the emission and magnetic field geometry of the accretion column (Nishimura, 2003;Schönherr et al., 2007;Schwarm et al., 2017).

Constraining the surface magnetic field strength from quiescent observations of Be X-ray binaries
Details regarding the formation of cyclotron lines in the spectra of accreting NSs in HMXBs are highly debated in the literature.The observed positive and negative correlations of the cyclotron energy with luminosity are often interpreted as a change of the dynamics of the accretion process in different accretion Figure 10.Absorption column and cyclotron line variability in Vela X-1, calculated from simulated phaseresolved HEX-P spectra using the results of Diez et al. (2022) with an average flux of 6.7×10 −9 erg s −1 cm −2 in the 3-80 keV band.The parameters from top to bottom are: absorption column N H , fundamental line energy (E 1 ), fundamental line strength (d 1 ), n 2 harmonic line energy (E 2 ), and n 2 harmonic line strength (d 2 ).Note that the fundamental line width is set to half of the n 2 harmonic width.The black uncertainties show the standard deviation for a sample of 100 simulations per phase bin.The orange dashed histogram indicates the expected 90% uncertainties on each individual realization.The green line indicates the input values for each phase bin.For details about the simulation set-up see § 3.2.2 and Table S3 in the Supplemental Material.
regimes, distinguished by a critical luminosity where a radiation-dominated shock forms in the column above the NS surface.For high-luminosity regimes, when it is assumed that L X > L crit , several mechanisms for the formation of cyclotron lines have been suggested that explain the observed negative line energy versus luminosity correlation.These scenarios involve different locations for the line production: above the NS surface in an accretion column that is growing with increasing luminosity (see, e.g., Becker et al., 2012) or in the illuminated atmosphere of the NS at lower magnetic latitudes, where the column radiation is reprocessed (as suggested by Poutanen et al. 2013;see, however, Kylafis et al. 2021).For lower luminosities, L X ≲ L crit , the formation of cyclotron lines is usually attributed to comparatively lower heights in the accretion column.In this case, a positive line energy versus luminosity correlation can be explained by the formation of a collisionless shock (see, e.g., Rothschild et al., 2017;Vybornov et al., 2017) that is moving closer to the NS surface for higher luminosities or by the redshift due to bulk motion of the accretion flow (Nishimura, 2014;Mushtukov et al., 2015).However, to apply these models in a consistent way and to distinguish between their predictions, it is required to know the surface field at the magnetic pole of the NS.This can then serve as a reference to estimate, for example, the characteristic height of the line-forming region in the column or the velocity of the accretion flow near the surface.
Recent evidence of accretion in Be X-ray binaries in quiescence, i.e., with L X ≪ L crit , and the simple emission region geometry expected in this case (i.e., a hot spot on the NS surface) together represent a unique opportunity to observe cyclotron lines at the energy corresponding to the surface value of the magnetic field.From current observational examples it seems that for this accretion state to occur, the magnetic field of the source might have to be sufficiently high, B ≳ 5 × 10 12 G, and the spin period might have to be comparatively long, ≳ 100 s (as, e.g., for GX 304−1 and GRO J1008−57: Tsygankov et al., 2019b;Lutovinov et al., 2021).The low flux level, however, makes the detection of cyclotron lines at high energies very challenging.GX 304−1 is one of the first sources which unambiguously exhibited stable quiescent accretion (Rouco Escorial et al., 2018) with a cyclotron line known from outburst observations (Mihara et al., 2010;Malacaria et al., 2015).The corresponding flux level in the 2-10 keV energy band, ∼0.4 mCrab, is below the sensitivity of Insight/HXMT.With NuSTAR we can access the high-energy emission, but are typically unable to constrain the turnover of the second hump of the characteristic double-hump spectrum (see §3.1) and the cyclotron line (Tsygankov et al. 2019b, Zainab et al., in prep.).
The best way to probe the above mentioned science cases is to observe the source during the quiescent accretion regime, i.e., at a low-luminosity state that is typically below the detection threshold for X-ray all sky monitors and accessible only with pointed observations subsequent to an X-ray outburst.Such a regime is also of importance since the simplified physics of plasma stopping at the nearly-static NS atmosphere allows for more detailed modeling of emission processes.The tenuous flow of matter stops at the NS atmosphere by Coulomb collisions, resulting in a high temperature gradient from ∼30 keV at the top down to ∼2 keV in the lower layers, separated by only ∼10 m (see, e.g., Sokolova-Lapa et al., 2021).The intrinsic emission is mainly produced by magnetic bremsstrahlung (for sufficiently low magnetic fields, cyclotron photons are produced as well from collisional excitations of electrons moving in bulk; Mushtukov et al., 2021), which is then modified by Compton scattering.This regime of accretion allows, for the first time, to combine modeling of the temperature and density structure of the emission region with a joint simulation of the continuum and cyclotron line formation (Mushtukov et al., 2021;Sokolova-Lapa et al., 2021).
In this way, existing physically-motivated models principally provide access to information about the field strength at the poles of accreting highly magnetized NSs.However, the limitations of current high-energy missions do not permit us to constrain the corresponding cyclotron lines in the spectra.HEX-P will open a In order to demonstrate HEX-P's capabilities to measure the surface magnetic field strength of a NS in quiescence, we simulate observations of the quiescent state of the Be X-ray binary GX 304−1, the first system for which a transition to the two-component spectrum was observed by NuSTAR (Tsygankov et al., 2019b).We use the physical model polcap presented by Sokolova-Lapa et al. (2021), which describes emission for the low-luminosity accretion regime.Here, we adopt an updated parametrization of the model (using the same underlying pre-calculated spectra; E. Sokolova-Lapa, priv. comm.).The parameters are the mass flux, Ṁ/πr 2 0 , where Ṁ is the mass-accretion rate and r 0 is the polar cap radius; the cyclotron energy, E cyc , corresponding to the polar magnetic field strength, B; and the normalization given in terms of r 2 0 /D 2 , where D is the distance to the source.We set the normalization and the mass flux based on the flux level and the shape of the spectrum as observed previously by NuSTAR.The magnetic field strength is set to B = 5.35 × 10 12 G, derived from the centroid energy of the cyclotron line observed during the previous outburst (∼50 keV, Jaisawal et al., 2016) and corrected by the gravitational redshift near the surface of the NS (z g ≈ 0.24, assuming standard NS parameters).The corresponding spectrum calculated with the polcap model and the exact values of the parameters used for the simulations are shown in Figure 11.Due to internal averaging over the emission angles to obtain the total flux in the NS rest-frame (see details in Sokolova-Lapa et al., 2021), the cyclotron line in the corresponding spectrum is located at ≈60 keV.Taking into account the gravitational redshift, the cyclotron line in the observed spectrum is therefore expected at ≈48 keV.We simulate a 60 ks observation for all three HEX-P instruments, combining the polcap model with tbabs (Wilms et al., 2000) to account for interstellar absorption (i.e., tbabs*polcap), fixing N H to the Galactic value of 1.1 × 10 22 cm −2 in the direction of the source.The resulting luminosity in the 1-80 keV range is ∼8 × 10 33 erg s −1 .
We compare the HEX-P spectra against a NuSTAR observation of the same exposure.We first fit both, the simulated HEX-P and existing NuSTAR observations, using a model that includes two independent Comptonized components (comptt;Titarchuk, 1994), tbabs*(comptt 1 + comptt 2 ).This model, with or without a multiplicative Gaussian-like cyclotron line (gabs), is commonly used to describe the twocomponent spectra of low-luminosity states (see, e.g., Tsygankov et al., 2019a;Lutovinov et al., 2021;Doroshenko et al., 2021).Similarly to earlier analyses (Tsygankov et al., 2019b), we obtain a good description of the NuSTAR data with this model, with χ 2 red = 136.93/107= 1.28 for the best fit.For the simulated HEX-P data, the same model consistent of two absorbed Comptonized components provides a formally satisfactory fit (χ 2 red = 322.64/275= 1.17), however, the residuals are flat only at low and intermediate energies, but indicate a dip at around 40-60 keV, where the cyclotron line is expected from the underlying physical model (see Figure 11).The best fit for a model which includes an additional Gaussian absorption line at high energies to describe the cyclotron line, provides the line's centroid energy of 46.5 +4.1 −2.7 keV , width 8.2 +3.1 −2.1 keV , and strength 54 +41 −21 keV, with the χ 2 red = 261.23/276= 0.95.The centroid energy corresponds (within its uncertainties) to the redshifted cyclotron line from the physical model, which is expected to be at ≈60/(1 + z g ) keV ≈ 48 keV, assuming the standard NS parameters to estimate the redshift, z g = 0.24. Figure 12 shows the resulting spectra and the corresponding best-fit models for the real NuSTAR and the simulated HEX-P observations.This example shows that even at this exposure, which is relatively low for the quiescent state of Be X-ray binaries, we can constrain the cyclotron line energy and thus the surface magnetic field to within ∼20%.This capability is crucial for systems which have never been observed at high-luminosity states, that is, for which cyclotron lines have never been observed in the spectra.This class of Be X-ray binaries accreting in quiescence is actively growing, e.g., through NuSTAR follow-up of Be X-ray binary candidates found during the first eROSITA All Sky Survey (see, e.g., Doroshenko et al., 2022).More similar discoveries are expected as deeper X-ray survey data from eROSITA become available.The spectral shape of the latter sources, in particular, the location of the high-energy excess associated with the energies below the cyclotron resonance, suggests a similar magnetic field strength as in GX 304−1.The search for high-energy cyclotron lines, which makes it possible to unambiguously constrain the surface magnetic field of an accreting neutron star, requires superior sensitivity at hard X-rays, as will be provided by HEX-P.
For a discussion of NS science that can be done with HEX-P for B-fields in excess of ∼10 14 G (i.e., magnetars), we refer the reader to Alford et al. 2023, (in prep.).

Constraining super-critical accretion and the critical luminosity via cyclotron line evolution in extragalactic sources
A correlation between cyclotron line energy and X-ray luminosity has been reported for a handful of X-ray pulsars (see Staubert et al., 2019, and above).The dependence of the line energy on luminosity can be attributed to changes in the height of the accretion column and can be positive or negative for low and high X-ray luminosity, respectively (Becker et al., 2012).However, in a couple of sources a secondary effect has been reported: for equal levels of luminosity, the energy of the cyclotron line can be different (up to 10% change) between the rise and the decay of an outburst or between different outbursts (e.g., V 0332+53, Cusumano et al. 2016 andSMC X-2, Jaisawal et al. 2023).Doroshenko et al. (2017) proposed that this effect is likely caused by a change of the emission region geometry (e.g., different combinations of height and width of the accretion column), while an alternative explanation is that it is due to accretion-induced decay of the NS's magnetic field (Cusumano et al., 2016;Jaisawal et al., 2023).
Given this complex observational behavior, it is critical to detect CRSFs in more transient systems and to study their evolution over luminous outbursts.Although a significant population of HMXBs is found in the Milky Way, interpreting results from observing its members is often hampered by their uncertain distance and large foreground absorption.In that sense nearby galaxies like the star-forming Large and Small Magellanic Clouds (LMC and SMC) offer a unique laboratory to complement our studies of luminous Galactic HMXBs.Sources in these galaxies have well determined distances of ∼50 kpc (LMC) or ∼60 kpc (SMC) and low Galactic foreground absorption (∼10 20 cm −2 ) making them ideal targets for spectral and temporal studies during major outbursts.Based on past observations and recent statistics, an outburst that peaks above 2 × 10 38 erg s −1 occurs in the Magellanic Clouds every few years (e.g.Vasilopoulos et al., 2014;Koliopanos and Vasilopoulos, 2018;Maitra et al., 2018;Vasilopoulos et al., 2020).Such outbursts are brighter than the Eddington luminosity for a typical NS, more precisely, they are brighter than the critical luminosity5 , L crit ∼10 37 erg s −1 , for a typical accretion column (see §3.1).They are therefore called super-Eddington or super-critical outbursts.SMC X-2 is a good example to demonstrate HEX-P's capabilities to detect and follow the evolution of a CRSF during a bright outburst.This BeXRB pulsar has exhibited two outbursts which reached luminosities above the Eddington limit over the last decade, one in 2015 and another in 2022.Both outbursts were followed up with Swift/XRT, as well as with NuSTAR ToOs (three observations in 2015 and one in 2022), covering a broad range in luminosity of (2 − 6) × 10 38 erg s −1 (see Fig. 13, left panel).When comparing the CRSF energy for the two outbursts at the same luminosity level, the line energy was about 2 keV higher during the 2022 outburst than observed previously in 2015 (see Fig. 13, right panel).This corresponds to a difference in magnetic field strength of B ∼ 2 × 10 11 G, akin to the difference reported by Cusumano et al. (2016) for V 0332+53.
In order to test HEX-P's capabilities to detect CRSF features for extragalactic sources and determine if the data would be sufficient to distinguish trends with luminosity, we performed several simulations based on the spectral properties of SMC X-2 as measured by NuSTAR (Jaisawal et al., 2023).For all simulations we used both HET and LET detectors.For the continuum we used an absorbed (tbabs with N H = 10 20 cm −2 ; Wilms et al., 2000) cut-off power-law model (i.e.cutoffpl) and a cyclotron absorption feature with the gabs model.We limited our simulations to the super-critical regime and a bolometric  (Jaisawal et al., 2023).Note that the NuSTAR points on the trend 1 line are from the 2015 outburst, while the point on the trend 2 line was measured during the 2022 outburst; hence different cyclotron line energies can be measured even when the source is at nearly the same luminosity during outburst indicating an inherent difference between outbursts.We demonstrate HEX-P's capabilities to discern between different possible trends in cyclotron line energy as a function of luminosity in these types of systems in as little at 10 ks, though greatly improved for 50 ks exposures.luminosity of 1.5-6×10 38 ergs s −1 which yields an absorbed flux range of 2-8×10 −10 ergs cm −2 s −1 (1-50 keV).We simulated spectra where the CRSF energy varied with luminosity assuming two arbitrary inverse linear energy dependencies (trend 1 and trend 2 in Fig. 13).We then fit the spectra and derived the best-fit parameters for the line with corresponding uncertainties.The results are shown in Figure 13 (right panel).HEX-P can constrain the CRSF parameters (here the central line energy) at least two times better than NuSTAR and with an exposure time (10 ks) that is only between half and a fourth that of NuSTAR, as well as discern between different trend inputs.For a longer exposure time of 50 ks the line parameters can be measured with unparalleled accuracy at the distance of the SMC.Thus, HEX-P will be able to detect and follow cyclotron line evolution for super-critical outbursts in the Magellanic Clouds and thereby test competing physical models (i.e.similar to V 0332+53 Cusumano et al., 2016;Vybornov et al., 2018).Note that the flux regime for which these trends are simulated (see Table S1) and the science cases of the previous subsections suggest that trends in CRSF energy can also be probed at sub-critical through super-critical luminosities in galactic sources.For science that can be executed with HEX-P regarding pulsations and cyclotron lines in extragalactic ultraluminous X-ray (ULX) sources which reach even higher luminosities of > 100× their Eddington limit, we refer the reader to Bachetti et al. (accepted).
In addition, it is expected that accreting XRPs crossing the critical luminosity should show a reversal of the cyclotron line energy versus luminosity correlation.To date only marginal evidence has been observed for such a trend reversal and in only a couple of sources (Doroshenko et al., 2017;Malacaria et al., 2022).As mentioned above, the Magellanic Clouds host comparatively many transient accreting pulsars displaying high-luminosity outbursts.HEX-P's ability to constrain the cyclotron line energy with unprecedented accuracy across orders of magnitude in luminosity (and thus over distinct accretion regimes) will provide us with the opportunity to trace more and well defined trend reversals.This will allow us to constrain the critical luminosity for a given source which, in turn, places constraints on physical parameters of the observed system, such as the NS magnetic field strength or source distance (Becker et al., 2012).

CONCLUSION
We have presented open questions about NSs and accretion onto strongly magnetized sources, and demonstrated HEX-P's unique ability to address these questions.The broad X-ray passband, improved sensitivity, and low X-ray background make HEX-P ideally suited to understand accretion in X-ray soft sources, down into low accretion rate regimes.In particular for LMXBs, we have shown that HEX-P observations will discriminate between competing continuum emission models that diverge above 30 keV; an energy band inaccessible to current focusing X-ray telescopes for these soft spectrum sources.Additionally, leveraging the improved sensitivity and broad X-ray passband, HEX-P will achieve tighter constraints on NS radius measurements through reflection modeling.These measurements are complementary and independent to other methodologies, and will narrow the allowed region on the NS mass-radius plane for viable EoS models for ultra-dense, cold matter, providing fundamental physics information in a regime inaccessible to terrestrial laboratories.
For the case of NSs in HMXBs, HEX-P will vastly improve our ability to: 1) detect multiple cyclotron line features in a single observation, aiding in our understanding of the magnetic field strength in different regions of the accretion column and the poorly known physical mechanisms behind the formation of higher harmonic features; 2) obtain detailed phase-resolved spectra to track the dependence of CRSFs on pulse phase in order to explore the geometric configuration of the accretion column emission in unprecedented details; 3) constrain the surface magnetic field strength of accreting NSs at the lowest accretion regime and characterize the continuum spectral formation to provide updated physically-motivated spectral models for future observations; and 4) identify CRSFs in extragalactic sources and follow the evolution of their line energy as a function of luminosity in order to distinguish between competing theories regarding changing emission region geometry or accretion-induced decay of the NS B-field.
HEX-P will provide a new avenue for testing magnetic field configuration, emission pattern, and accretion column physics close to the surface of NSs, as well as enhance our understanding of extreme accretion physics in both LMXBs and HMXBs.
California Institute of Technology, under a contract with NASA.E.S.L. and J.W. acknowledge partial funding under Deutsche Forschungsgemeinschaft grant WI 1860/11-2 and Deutsches Zentrum für Luftund Raumfahrt grant 50 QR 2202.ACKNOWLEDGMENTS R.M.L. and A.W.S. would like to thank Dr. Kazuhiko Shinki for providing consultation on statistical methods of §2.2.2.The authors thank the reviewers for their detailed comments that enhanced the science cases presented within the manuscript.

SUPPLEMENTAL MATERIAL
Here we provide the input values used for creating the various HEX-P simulations that are presented in the accompanying manuscript.Additionally, we specify the flux levels of the simulations since the introduction mentions a pile-up limit for the LET in units of mCrab.

Figure 2 .
Figure 2. Position of the innermost stable circular orbit (ISCO) in units of gravitational radii versus the dimensionless spin parameter a of the CO.The hatched region indicates where the NS surpasses the rotational limit and would break apart.The horizontal dot-dashed lines indicate the corresponding values of 10 km and 12 km radius for a canonical NS mass of 1.4 M ⊙ .The horizontal dotted lines indicate the same radius values for a more massive NS of 2.0 M ⊙ .

Figure 3 .
Figure 3.Comparison of simulated 20 ks HEX-P (LET+2 HETs) observation for two different continuum models typically used to describe NS LMXB soft state spectra: a double thermal continuum model with weak Comptonization (Model 1: black) and an accretion disk with blackbody Comptonization (Model 2:blue).The models diverge above 30 keV where currently available missions quickly become background dominated for the same exposure time.To emphasize this point, the same S/N near 30 keV could be achieved with NuSTAR in an exposure time of 96.6 ks.The broad X-ray passband and improved sensitivity of HEX-P will differentiate the models.

Figure 4 .
Figure 4. Posterior distribution of a and R in for 10 4 iterations of simulating HEX-P spectra for three spin values: (a) a = 0.0, (b) a = 0.17, and (c) a = 0.3.The 1σ, 2σ, and 3σ contours are shown.The white points indicate the highest probability value.The dotted grey line indicates a constant line of 6 R g which corresponds to the ISCO for a = 0.The distribution tightens at high values of spin as relativistic effects become stronger.However, the distributions show that the data are able to trace out unique values of inner disk radii by trending towards lower values of R g at higher a.

Figure 5 .
Figure 5. Mass and radius constraints from reflection modeling compared to NS gravitational wave events and NICER pulsar light curve modeling.The darker solid orange region indicates the improved radius constraints for reflection modeling of HEX-P data based on Cygnus X-2 for the case of high a.Note that both R in and a are free parameters when fitting the HEX-P data, while the NuSTAR plus NICER analysis fixes a = 0 (lighter solid orange region in panel a; Ludlam et al. 2022).The solid gray region indicates where causality is violated (i.e., the sound speed within the NS exceeds the speed of light).Pulsar light curve modeling of NICER data for PSR J0740+6620 is indicated in teal (dashed lines: Miller et al. 2021; solid lines: Riley et al. 2021).Maroon indicates the results for light curve modeling of PSR J0030+0452 (dashed lines: Miller et al. 2019; solid lines: Riley et al. 2019).The black dotted line denotes the mass-radius constraints from the combined GW170817 (Abbott et al., 2019) and GW190425 (The LIGO Scientific Collaboration et al., 2020) signatures using a piece-wise polytropic model as reported in Raaijmakers et al. (2021).Confidence contours correspond to the 68% and 95% credible regions.Panel (b) shows select EoS models from Lattimer and Prakash (2001) to demonstrate the behavior of different internal compositions on the mass-radius plane.

Figure 6 .
Figure 6.Simplified depiction of the formation of cyclotron lines in the spectra of accreting highly magnetized neutron stars due to photon scattering off electrons, whose motion perpendicular to the field lines is restricted by Landau quantization.The latter leads to strong cyclotron resonances in the Compton scattering cross section, broadened by the electron thermal motion parallel to the magnetic field.The centroid energies of the cyclotron lines in observed spectra approximately correspond to the redshited energies of the resonances (with some correction for the scattering redistribution effects during the lines' formation).The gravitational redshift is dependent on the height of the line-forming region, h, above the surface of the neutron star.

Figure 7 .
Figure 7. Simulated 20 ks HET spectra for Cen X-3, as well as the best-fit model are shown in panel (a).Panels (b), (c), and (d)show the residuals obtained when we set the strength of fundamental, n 2 harmonic, and n 3 harmonic to zero, respectively.We re-binned the data visually and provide arrows for the cyclotron line components to aid with clarity.For comparison, we provide a simulated NuSTAR spectrum for using the same model and exposure time which is shown in panel (e).The NuSTAR data are strongly background dominated at these flux levels, which makes constraining the higher energy harmonics difficult even at much longer exposure times (see text for more details).

Figure 11 .
Figure 11.Model used to simulate the observation of the Be X-ray binary GX 304−1 in quiescence.The vertical orange line indicates the center of the cyclotron line.The cyclotron line and the continuum are calculated together using polarized radiative transfer simulations (Sokolova-Lapa et al., 2021).The low-energy "thermal" hump and the high-energy hump below the cyclotron line (the red wing of the cyclotron line) are typically observed as a "two-hump" spectrum from Be X-ray binaries in quiescence.

Figure 12 .
Figure12.NuSTAR observation (left, a) and simulated HEX-P (right, a) spectra of a 60 ks observation of the Be X-ray GX 304−1 in quiescence and corresponding best-fit models.Both panels b show residuals for the best-fit models.For the simulated HEX-P spectra, we also show the residuals for the case when the cyclotron line strength is set to zero (left, c).The high-energy cyclotron line can be constrained from the HEX-P data, with the centroid energy 46.5 +4.1 −2.7 keV.

Figure 13 .
Figure 13.Left: The 2015 and 2022 super-critical outbursts of SMC X-2 observed by Swift/XRT.The vertical lines mark the epochs of four NuSTAR observations of the system.Right: The cyclotron line dependence on luminosity.The star symbols mark the measured CRSF values for SMC X-2 for the NuSTAR observations(Jaisawal et al., 2023).Note that the NuSTAR points on the trend 1 line are from the 2015 outburst, while the point on the trend 2 line was measured during the 2022 outburst; hence different cyclotron line energies can be measured even when the source is at nearly the same luminosity during outburst indicating an inherent difference between outbursts.We demonstrate HEX-P's capabilities to discern between different possible trends in cyclotron line energy as a function of luminosity in these types of systems in as little at 10 ks, though greatly improved for 50 ks exposures.

Table S2 .
Input spectral parameters for the two continuum models used to demonstrate HEX-P's capabilities to distinguish continuum shapes in NS LMXBs in section 2.2.1.The int type parameter of nthcomp in Model 2 is set to 0 so that the seed photons come from a single-temperature blackbody.