© The Authors 2026

1. Introduction

Accretion-powered X-ray pulsars (XRPs) are highly magnetized neutron stars (NSs) in binary systems accreting from a companion, typically an OB supergiant or a Be star (for a recent review see Weng & Ji 2024). Due to the star’s strong magnetic field (usually a few 10¹² G), accretion proceeds along the magnetosphere, and the plasma is channeled onto the magnetic poles (Ghosh & Lamb 1979). This process heats the magnetic polar regions, creates hotspots, and, under certain conditions, produces accretion columns, where the gravitational energy of the relativistic inflowing plasma is released primarily via hard X-ray emission (Basko & Sunyaev 1975, 1976; Arons et al. 1987).

A substantial fraction of these sources shows absorption-line-like features in the 10–100 keV energy range of their spectra, the so-called cyclotron lines or cyclotron resonant scattering features (CRSFs; for a review, see Staubert et al. 2019). These features originate from a quantum-mechanical process called magnetic scattering. The perpendicular momentum of charged particles in a strong magnetic field is quantized, resulting in discrete energy states, known as Landau levels (Landau 1930). Resonant scattering of photons with inflowing electrons in the vicinity of the NS surface imprints the fundamental CRSF in the continuum (Trümper et al. 1977; Truemper et al. 1978), with several XRPs also featuring harmonics (Truemper et al. 1978; Santangelo et al. 1999; Fürst et al. 2018; Yang et al. 2023). The energy difference between two successive Landau levels (hereafter referred to as the cyclotron energy) in the nonrelativistic case is

$$\begin{array}{c}\hfill E_{\mathrm{c}}\equiv \frac{ħeB}{m_{\mathrm{e}}c}\approx 11.6\left(\frac{B}{{10}^{12}\mathrm{G}}\right)\mathrm{keV},\end{array}$$(1)

where B is the magnetic field strength in the line forming region, ℏ is the reduced Planck constant, m_e is the electron’s rest mass, c is the speed of light, and e is the elementary charge. As E_c is linearly proportional to the magnetic field, the detection of CRSFs provides a direct and robust means of measuring the magnetic field strength in XRPs and, more broadly, serves as a diagnostic of the physical conditions near the neutron star surface Gnedin & Sunyaev (1974).

Detailed observations of CRSFs in XRPs have revealed significant variability both with pulse phase and X-ray luminosity, L (for reviews, see Staubert et al. 2019; Mushtukov & Tsygankov 2022). Pulse-phase modulations indicate anisotropic emergent spectra (e.g., Maniadakis et al. 2025), whereas variations with L (accretion rate) hint at changes in the emission region (e.g., Shui et al. 2024). Evidence of long-term CRSF variability has also been observed in Her X-1 (Staubert et al. 2017; Xiao et al. 2019), but the underlying physical mechanism remains unclear.

Of particular interest is the dependence of the observed CRSF centroid energy, E_CRSF, on the X-ray luminosity, L. Such behavior has been reported in several XRPs, including V0332+53 (Tsygankov et al. 2006, 2010), GX 304−1 (Klochkov et al. 2012; Malacaria et al. 2015; Rothschild et al. 2017), and GRO J1008−57 (Chen et al. 2021). High-luminosity sources (L ≳ 10³⁷ erg s⁻¹) typically present a negative correlation between E_CRSF and L, whereas in low-luminosity sources (L ≲ 10³⁷ erg s⁻¹) the trend reverses, with a positive correlation emerging. A complete transition from positive to negative correlation has recently been identified in pulse-to-pulse analyses of 1A 0535+262 by Shui et al. (2024), further reinforcing this distinction among sources. The leading physical explanation is the emergence of an accretion column once the luminosity exceeds a critical value, L^* ∼ 10³⁷ erg s⁻¹ (Basko & Sunyaev 1976; Mushtukov et al. 2015a), signaling a transition in the accretion regime (Becker et al. 2012); from a surface-dominated hotspot regime (L ≲ L^*) to an accretion-column regime (L ≳ L^*).

For supercritical (accretion-column-dominated) sources, it has been shown analytically (Basko & Sunyaev 1976; Becker & Wolff 2007) and through first-principles radiation-magnetohydrodynamic (RMHD) simulations (Klein et al. 1996; Zhang et al. 2022, 2023, 2025a) that the braking of the infalling plasma is caused by an extended radiation-mediated shock (RS) forming in the magnetically confined, optically thick, accretion column. Both analytical studies (Basko & Sunyaev 1976) and Monte Carlo (MC) approaches (Loudas et al. 2024b) support the RS as the likely formation site of CRSFs in high luminosity XRPs. This scenario successfully predicts the observed anti-correlation between E_CRSF and L: As the accretion rate increases, the accretion column grows in height, while the magnetic field strength drops with altitude. Thus, the cyclotron energy in the line-forming region shifts to lower values at higher luminosities (Burnard et al. 1991; Loudas et al. 2024a). An alternative scenario, proposed by Poutanen et al. (2013), attributes CRSF formation to surface reflection of hard X-rays emitted from the accretion column, but detailed radiative-transfer simulations challenge the viability of this interpretation (Kylafis et al. 2021).

At low luminosities (hot-spot-dominated emission), the outgoing radiation is dynamically unimportant, the accreting plasma reaches the NS surface in free fall, and the formation of CRSFs is believed to occur within the NS atmosphere. Even though no settled theory exists explaining the deceleration of the falling material, two leading mechanisms are found in the literature: 1) braking of plasma in the atmosphere by Coulomb collisions (Zel’dovich & Shakura 1969; Miller et al. 1987; Sokolova-Lapa et al. 2021) and 2) emergence of a collisionless shock (CS) where the infalling plasma is thermalized upon passing through the shock (Bisnovatyi-Kogan & Fridman 1970; Langer & Rappaport 1982; Bykov & Krasil’Shchikov 2004). The former model considers spectral formation in the NS atmosphere and can explain the observed double-hump in the hard X-ray spectra of low-luminosity XRPs (Mushtukov et al. 2021; Sokolova-Lapa et al. 2021), but it predicts zero variability of the cyclotron line centroid with luminosity and thus cannot reproduce the observed trends of cyclotron lines.

Conversely, the CS scenario has been successful in explaining the positive correlation between E_CRSF and L in several sources (Staubert et al. 2007; Vybornov et al. 2017; Rothschild et al. 2017; Roy & Sharma 2025). In that case, spectral formation takes place between the NS surface and the shock discontinuity. Since the magnetic field strength (and effectively E_CRSF) increases as the altitude decreases, and the CS height decreases with accretion rate (luminosity L) (Bykov & Krasil’Shchikov 2004), a positive correlation between E_CRSF and L arises (see also Becker et al. 2012). However, first-principles kinetic simulations are required to address whether the formation of such a shock is feasible under conditions encountered above the NS atmosphere and determine its properties. Once the shock dynamics is understood, subsequent radiative transfer calculations need to be conducted, followed by comparison with observed spectra and pulse profiles.

Mushtukov et al. (2015b) pointed out that for subcritical luminosities (L ≲ L^*) even though the propagation of the hotspot’s radiation through the mildly relativistic infalling plasma cannot halt accretion, it can significantly modify the velocity and density profiles above the NS surface. This effect was verified by robust RMHD simulations by Markozov et al. (2023). In this scenario, the formation of the CRSF takes place above the hotspot. Upward-propagating X-rays undergo resonant scattering with the infalling electrons, producing a redshifted CRSF relative to the local cyclotron energy due to the Doppler effect between the laboratory frame (hotspot) and the moving electron’s bulk frame. This model predicts a positive correlation between E_CRSF and L without postulating the presence of a CS. A qualitative explanation follows. Radiation pressure acts to decelerate the accretion flow, influencing the degree of redshift in the CRSF. As luminosity (and hence radiation pressure) increases, the accretion-flow velocity right above the hotspot decreases, reducing the Doppler shift and shifting the CRSF centroid to higher energies (closer to the local E_c). Although the model shows promise and has been successful in qualitatively explaining observations of the X-ray pulsar GX 304−1 (Mushtukov et al. 2015b), no quantitative calculations of cyclotron line formation and variability have yet been carried out within this framework.

In this paper, we investigate the aforementioned model by conducting robust MC radiative transfer simulations to obtain the angle-resolved spectrum for a typical X-ray pulsar, focusing on the properties of the emergent CRSF, and studying its dependence on accretion rate. We aim to quantitatively test the following predictions made by Mushtukov et al. (2015b):

• (i)

  E_CRSF is positively correlated with the luminosity L.

• (ii)

  The line width, σ_CRSF, is positively correlated with L.

• (iii)

  Both E_CRSF and σ_CRSF vary significantly over the pulse phase.

The structure of this paper is organized as follows. In Sect. 2, we introduce the physical model, detailing the velocity and density profiles, along with the structure of the accretion channel. In Sect. 3, we provide a detailed description of our MC code. In Sect. 4, we present the simulated spectra for a typical X-ray pulsar and analyze their primary features. In Sect. 5, we compare our MC results with observational data from GX 304−1 and discuss the degeneracies inherent in the methodology. In Sect. 6, we comment on the results and draw our conclusions.

2. The model

In this section, we outline the underlying features of the model developed by Mushtukov et al. (2015b), which applies to sources in the subcritical luminosity regime (L ≲ L^*).

2.1. Structure of the accretion funnel

We considered magnetospheric accretion onto the magnetic pole of a highly magnetized NS. The strong magnetic field confines the accretion flow into a funnel, and the plasma moves along the field lines only. The impact of infalling plasma heats the NS atmosphere at the base of the funnel, forming a hotspot that produces an outgoing photon beam. For simplicity, we adopted a cylindrical geometry for the filled accretion column, although more realistic accretion flow models favor bow-shaped or annular hotspots (Basko & Sunyaev 1976; Romanova et al. 2004; Zhu 2025). We injected seed photons isotropically upward from the base of the cylinder with a blackbody spectrum. We further assumed a uniform magnetic field B aligned with the magnetic axis, i.e., $\mathit{B}=B\widehat{z}$. This approximation is justified because the vertical extent of the line-forming region and the size of the funnel are, by construction, much smaller than the NS radius.

The interaction between upward-moving photons and the infalling electrons through Compton scattering exerts a drag (radiation-pressure) force on the falling plasma, altering its velocity and density profiles in a nontrivial manner. In the subcritical luminosity regime considered here, the plasma reaches the surface with a nonzero terminal velocity before being completely thermalized via Coulomb collisions in the atmosphere (Zel’dovich & Shakura 1969). Conservation of momentum in the nonrelativistic regime for a stationary flow dictates

$$\begin{array}{c}\hfill (\mathit{u}·\mathrm{\nabla })\mathit{u}=\frac{1}{\rho }{\mathit{G}}_{\mathrm{rad}},\end{array}$$(2)

where ρ and u are the plasma density and bulk velocity, respectively, while G_rad is the radiation-pressure force, proportional to the radiation flux times the opacity. In this equation, we neglected the gas pressure and gravity terms, as the radiation-pressure force is the dominant one in the vicinity of the NS surface. Mushtukov et al. (2015b) solved Eq. (2) and derived the following analytical expression for the velocity profile β(h)≡|u(h)|/c in the center of the column, considering isotropic emission from the hotspot, accounting for the dependence of the outgoing radiation flux on the height h above the NS surface, and considering motion along the field lines only,

$$\begin{array}{c}\hfill \beta (h)={[{\beta }_{\mathrm{ff}}^2-({\beta }_{\mathrm{ff}}^2-{\beta }_0^2)(1-\frac{2}{\pi }{tan}^{-1}\left(\frac{h}{d}\right))]}^{1/2}.\end{array}$$(3)

Here, d is the radius of the hotspot, estimated by the following power-law scaling (Mushtukov et al. 2015b)

$$\begin{array}{c}\hfill d\approx 6\times {10}^3{\mathrm{\Lambda }}^{-3/8}{L}_{37}^{9/35}{B}_{12}^{-3/14}{m}^{-81/140}{R}_6^{39/35}\mathrm{cm},\end{array}$$(4)

where Λ = 0.5 is a geometrical parameter related to the coupling of the magnetosphere to the disk (Ghosh & Lamb 1978), L₃₇ ≡ L/10³⁷ erg s⁻¹, B₁₂ ≡ B/10¹² G, m = M/M_⊙, and R₆ ≡ R/10⁶ cm. In Eq. (3), β_ff denotes the free-fall speed, divided by c, of the accreting plasma at the neutron-star surface, while β₀ represents the terminal velocity at the stellar surface (h = 0). The latter is related to the accretion luminosity through energy conservation, assuming that the kinetic energy of the accreting plasma is instantaneously converted into radiation via thermalization at the hotspot and Compton scattering above it. This yields a scaling relation of the form

$$\begin{array}{c}\hfill {\beta }_0={\beta }_{\mathrm{ff}}\sqrt{1-\frac{L}{L^{\ast }}},\end{array}$$(5)

with the critical luminosity L^* defined as the limiting luminosity at which the infalling material is fully decelerated (β₀ = 0) by radiation pressure at the stellar surface; above this threshold, an RS is expected to rise above the hotspot (Basko & Sunyaev 1976).

To ensure that the model is self-consistent, the velocity profile described by Eq. (3) necessitates a nonuniform density profile that obeys the conservation of mass

$$\begin{array}{c}\hfill \mathrm{\nabla }·(\rho \mathit{u})=0.\end{array}$$(6)

Given the cylindrical geometry and a constant accretion rate (Ṁ = const), the electron number density profile reads (Mushtukov et al. 2015b)

$$\begin{array}{c}\hfill n_{\mathrm{e}}(h)\approx 2\times {10}^{19}{L}_{37}^{3/5}{B}_{12}^{1/2}{\beta }^{-1}(h){\mathrm{cm}}^{-3}.\end{array}$$(7)

Finally, following Mushtukov & Tsygankov (2022) (see Sect. 4.3 therein), we estimate the effective temperature of the hotspot using the Stefan-Boltzmann law

$$\begin{array}{c}\hfill k_{\mathrm{B}}T=k_{\mathrm{B}}{\left(\frac{L}{2{\sigma }_{\mathrm{SB}}{S}_D}\right)}^{1/4}\approx 6.3L_{37}^{3/20}{\mathrm{\Lambda }}^{7/32}{m}^{13/80}{R}_6^{-19/40}{B}_{12}^{1/8}\mathrm{keV},\end{array}$$(8)

where σ_SB is the Stefan-Boltzmann constant and S_D denotes the hotspot area, given by Mushtukov et al. (2015b)

$$\begin{array}{c}\hfill S_D\approx 3\times {10}^9{\mathrm{\Lambda }}^{-7/8}{m}^{-13/20}{R}_6^{19/10}{B}_{12}^{-1/2}{L}_{37}^{2/5}{\mathrm{cm}}^2.\end{array}$$(9)

The factor of two in the denominator of Eq. (8) accounts for the existence of two antipodal hotspots.

The approximate power-law scaling adopted for the hotspot effective temperature captures the expected positive correlation with luminosity. However, the absolute normalization of the temperature remains uncertain at the level of a factor of ∼2, depending on the accretion geometry (e.g., the shape of the landing region, stochastic variability, and system inclination) and on the microphysics of the thermalization process in the NS atmosphere (e.g., CS versus Coulomb collisions). Deriving a robust prescription for temperature would require detailed RMHD simulations (such as those of Markozov et al. 2023) and is therefore beyond the scope of this work. Nevertheless, we emphasize that our conclusions regarding the CRSF properties and their dependence on luminosity and/or viewing angle are insensitive to the assumed injected photon spectrum, and hence to the specific choice of the hotspot (blackbody) temperature.

2.2. CRSF formation

As the accreting plasma descends toward the surface, upward-propagating photons undergo Compton scattering with the infalling electrons. For typical values of the electron density and hotspot size, the perpendicular Thomson optical depth of the accretion funnel lies in the optically thin regime,

$$\begin{array}{c}\hfill \tau \sim n_{\mathrm{e}}{\sigma }_{\mathrm{T}}d\lesssim 0.2,\end{array}$$(10)

such that most photons escape without scattering. However, magnetic scattering includes resonant processes in which photons with energies near the cyclotron resonance experience a strongly enhanced cross section (up to ∼10⁶σ_T, Canuto et al. 1971; Ventura 1979; Harding & Daugherty 1991; Sina 1996; Mushtukov et al. 2016; Loudas et al. 2021) and are effectively trapped within the medium. These resonant photons dominate the radiative deceleration of the infalling plasma (Mushtukov et al. 2015b): they undergo multiple resonant scatterings until they exchange sufficient energy with the electrons to shift out of resonance and escape the funnel. This process naturally produces an absorption-line-like feature accompanied by two broad, emission-like wings, reflecting photon redistribution in energy rather than true photon absorption, since the total number of photons is conserved in scattering.

In this study, we assumed that the magnetic field strength is well below the quantum electrodynamics (QED) critical value,

$$\begin{array}{c}\hfill B_{\mathrm{cr}}=\frac{m_{\mathrm{e}}^2{c}^3}{eħ}\approx 4.413\times {10}^{13}\mathrm{G},\end{array}$$(11)

such that the dimensionless field strength, b ≡ B/B_cr ≪ 1, and we also considered the infalling electron’s temperature T_e to be substantially smaller than the local cyclotron energy, that is, k_BT_e ≪ E_c. In short, we restricted our analysis to systems that satisfy

$$\begin{array}{c}\hfill k_{\mathrm{B}}{T}_{\mathrm{e}}\ll E_{\mathrm{c}}\ll m_{\mathrm{e}}{c}^2,\end{array}$$(12)

which is suitable for most XRPs in the hotspot regime (Becker & Wolff 2005). Under these conditions, the vast majority of electrons occupy the fundamental Landau level (n = 0), allowing us to safely limit our treatment to transitions between the ground state and the first excited Landau level (0 → 1 → 0) (Nobili et al. 2008; Loudas et al. 2021).

In the infalling electron’s rest frame, the relativistic cyclotron resonance energy is given by

$$\begin{array}{c}\hfill E_{\mathrm{res}}=m_{\mathrm{e}}{c}^2\frac{2b}{1+\sqrt{1+2b{sin}^2{\theta }_{\mathrm{rf}}}}=\frac{2E_{\mathrm{c}}}{1+\sqrt{1+2b{sin}^2{\theta }_{\mathrm{rf}}}},\end{array}$$(13)

where θ_rf is the angle between the photon momentum and the magnetic field in the electron frame. For small values of b, one would expect the CRSF centroid energy to appear at E_res ≈ E_c. Due to the bulk motion of the infalling plasma, this resonance energy appears Doppler shifted when viewed from the laboratory frame; thus, it is not photons with energy ≈E_res (in the lab frame) that are in resonance, but those that the Doppler boost to the electron frame brings them in resonance, that is

$$\begin{array}{c}\hfill E\approx \frac{E_{\mathrm{res}}}{\gamma (1+\beta cos\theta )}<E_{\mathrm{res}},\end{array}$$(14)

where θ ∈ (0, π/2) is the angle between the photon momentum and the magnetic field in the lab, and γ = (1 − β²)^−1/2. Therefore, the CRSF will not form at E_res but will appear redshifted. As Eq. (14) indicates, the amount of redshift depends on the plasma’s velocity in the region where the interaction takes place, as well as on the photon’s direction of propagation.

Higher luminosities correspond to lower infall velocities (see Eq. (5)), leading to higher CRSF centroid energies. This results in a positive correlation between the luminosity and the CRSF centroid energy (Mushtukov et al. 2015b). In addition, for a given velocity, the Doppler redshift decreases with increasing viewing (propagation) angle θ, implying that the CRSF centroid energy increases with θ. However, this qualitative picture is based on the single-scattering approximation. Multiple scatterings could induce modulations in the shape and location of the emergent CRSF. Obtaining accurate predictions of the model requires full radiative-transfer calculations.

3. Monte Carlo code

The relativistic MC radiative transfer code employed in this work is based on the implementation developed by Loudas et al. (2024b), to which the reader is referred for a detailed description. Here, we summarize the main parts of the code and describe the improvements with respect to the version of Loudas et al. (2024b). Radiative transfer is treated using the forced first collision technique, which ensures efficient sampling of photon–electron interactions. Photon packets are injected in the accretion funnel and are followed by the MC code until they escape. A key advantage of the present implementation is its flexibility, as it supports arbitrary density and velocity profiles within the accretion channel, enabling direct tests of different physical scenarios. Moreover, it incorporates an adaptive step size algorithm that increases efficiency, while maintaining accuracy in regions of steep gradients. Figure 1 shows an illustration of the system’s geometry and the numerical propagation of a photon packet.

Fig. 1. Illustration of the MC code’s geometry, coordinates, and calculation of the line-of-sight optical depth. A photon packet located at r0 is moving in direction $\widehat{n}$. We computed the accumulated optical depth by dividing the line of sight into a sequence of points r0 using an adaptive step-size algorithm. We then sampled a target optical depth, τrnd, and determined the position of the next scattering event via linear interpolation between the discretized points. We assumed a uniform magnetic field because the height of the line-forming region is much smaller than the NS radius.

3.1. Accretion channel and plasma

We modeled the accretion channel as an axisymmetric cylindrical column of radius d and semi-infinite length, with its base located at the NS surface. We considered the magnetic field to be homogeneous along the z axis, that is, $\mathit{B}=B\widehat{z}$, because the height of the line-forming region is much smaller than the NS radius. The electron density, n_e(r), evaluated at r = (x, y, z = h) is given by Eq. (7), while the bulk velocity of the infalling plasma is taken to be $\mathit{u}(\mathit{r})=-\beta (\mathit{r})c\widehat{z}$, with β(r) computed by Eq. (3). Even though both density and velocity profiles are expected to differ across the column (Markozov et al. 2023), for the purposes of this work, we utilized the analytical 1D profiles derived by Mushtukov et al. (2015b). We treated infalling electrons as a cold plasma, an assumption justified by observations (Becker & Wolff 2005); thus we set T_e = 0 in our simulations.

3.2. Seed photons

We emitted photon packets isotropically upward, with their initial positions sampled homogeneously from the base of the accretion channel. In our coordinate system, this base corresponds to the region ℛ = {(x, y, z = h)∈ℝ³ : x² + y² ≤ d², h = 0}. We initialized each photon packet with a specified weight, w = N (typically 10⁶), and assigned it an initial energy E drawn from a blackbody distribution at temperature T (given by Eq. 8),

$$\begin{array}{c}\hfill n(E)dE\propto \frac{E^2}{e^{E/k_{\mathrm{B}}T}-1}dE.\end{array}$$(15)

3.3. Optical depth calculation and photon propagation

We carried out the photon packet propagation as follows. We considered a photon packet with weight w, located at position r₀ and propagating in direction $\widehat{n}=(cos\phi sin\theta ,sin\phi sin\theta ,cos\theta )$ with energy E. Here, θ is the angle relative to the magnetic axis, while φ is the azimuthal one. The cumulative optical depth along the line of sight is defined as

$$\begin{array}{c}\hfill {\tau }_n({\mathit{r}}_{\mathbf{0}},\widehat{n},E)={\int }_{\mathcal{C}}{n}_{\mathrm{e}}(\mathit{r}){\sigma }_{\mathrm{eff}}(E_{\mathrm{rf}},\mathit{r},{\widehat{n}}_{\mathrm{rf}})ds,\end{array}$$(16)

where the trajectory 𝒞 is a straight line extending from r₀ to the cylindrical funnel’s boundary and ds is the infinitesimal scalar path length along 𝒞. The effective cross section σ_eff accounts for the bulk motion of the plasma and is defined as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{eff}}(E_{\mathrm{rf}},\mathit{r},{\widehat{n}}_{\mathrm{rf}})=(1+\beta (\mathit{r})cos\theta )\sigma (E_{\mathrm{rf}},{\theta }_{\mathrm{rf}}),\end{array}$$(17)

where σ is the magnetic scattering cross section expressed in terms of quantities evaluated at the electron rest frame (see Sect. 3.4) and β ≥ 0, by construction. The subscript “rf” denotes quantities measured in the electron rest frame. The Lorentz transformations relate the rest-frame quantities to the laboratory-frame ones as follows,

$$\begin{array}{c}\hfill cos{\theta }_{\mathrm{rf}}=\frac{\beta +cos\theta }{1+\beta cos\theta },{\phi }_{\mathrm{rf}}=\phi ,\end{array}$$(18)

and

$$\begin{array}{c}\hfill E_{\mathrm{rf}}=\gamma E(1+\beta cos\theta ),\end{array}$$(19)

where γ = (1 − β²)^−1/2.

To determine the location of the next scattering event, we computed numerically τ_n using an adaptive step-size algorithm to handle regions where the electron density n_e(r) and/or β(r) vary rapidly. We discretized the integration path into a sequence of points r₀, r₁, …, r_n along 𝒞, where r_n is the point of the line-of-sight’s intersection with the column boundary. Given a point r_i, we determined the subsequent point r_i + 1 as follows. First, we calculated the local mean free path

$$\begin{array}{c}\hfill {\lambda }_i\equiv \overline{\lambda }({\mathit{r}}_i,\widehat{n},E)={\left[n_{\mathrm{e}}({\mathit{r}}_i){\sigma }_{\mathrm{eff}}(E_{\mathrm{rf}},{\mathit{r}}_i,{\widehat{n}}_{\mathrm{rf}})\right]}^{-1}.\end{array}$$(20)

We then proposed a tentative step size equal to 1% of the local mean free path

$$\begin{array}{c}\hfill |{\mathit{r}}_{i+1}-{\mathit{r}}_i|=0.01{\lambda }_i.\end{array}$$(21)

To ensure linearity over this step, we evaluated λ_i + 1 at the new location and enforced the relative error condition

$$\begin{array}{c}\hfill \frac{|{\lambda }_{i+1}-{\lambda }_i|}{{\lambda }_i}<0.01.\end{array}$$(22)

If condition (22) was violated, we rejected the step and redefined r_i + 1 closer to r_i by halving the step size (backtracking), that is by increasing the spatial resolution locally until the condition was satisfied,

$$\begin{array}{c}\hfill {\mathit{r}}_{i+1}\to {\mathit{r}}_i+\frac{1}{2}({\mathit{r}}_{i+1}-{\mathit{r}}_i).\end{array}$$(23)

Once the step was accepted, the accumulated optical depth was updated as

$$\begin{array}{c}\hfill {\tau }_{i+1}={\tau }_i+\frac{|{\mathit{r}}_{i+1}-{\mathit{r}}_i|}{{\lambda }_i},\end{array}$$(24)

with the initial condition τ₀ = 0. This process continues until the boundary is reached (or τ > 20 for efficiency), yielding the total optical depth τ_n. A fraction we^−τ_n of the photon packet escapes the domain, and its energy and direction are recorded.

The remaining fraction underwent scattering with an electron in the column. The probability to scatter after covering a distance equivalent to optical depth τ is given by a truncated exponential distribution in τ, and thus we sampled a target optical depth τ_rnd using

$$\begin{array}{c}\hfill {\tau }_{\mathrm{rnd}}=-ln[1-\xi (1-e^{-{\tau }_n})],\end{array}$$(25)

where ξ ∈ [0, 1) is a uniform random number; ξ = 0 implies scattering without moving forward at all, while ξ → 1 indicates that the photon packet will scatter right before it reaches the boundary. We then located the index i such that τ_i ≤ τ_rnd < τ_i + 1. We found the precise location of the next scattering event, r_new, via linear interpolation along the segment,

$$\begin{array}{c}\hfill {\mathit{r}}_{\mathrm{new}}=\frac{{\mathit{r}}_i({\tau }_{i+1}-{\tau }_{\mathrm{rnd}})+{\mathit{r}}_{i+1}({\tau }_{\mathrm{rnd}}-{\tau }_i)}{{\tau }_{i+1}-{\tau }_i},\end{array}$$(26)

where |r_i − r₀|≤|r_new − r₀|< |r_i + 1 − r₀|.

Once the exact location of the scattering event was specified, we updated the weight of the remaining photon packet as follows,

$$\begin{array}{c}\hfill w_{\mathrm{new}}=w(1-e^{-{\tau }_n}).\end{array}$$(27)

This formula ensures conservation of the total number of photons. If w_new was larger than 1, we carried out the scattering event as described below; otherwise, we proceeded to the injection and propagation of a new photon packet (Sect. 3.2).

We treated the scattering process in three steps. First, we Lorentz-transformed the photon-packet energy and propagation direction from the laboratory frame to the electron rest frame using the local bulk velocity evaluated at the scattering location r_new (Eqs. (18) and (19)). Next, we sampled a new photon propagation direction $\widehat{n}{\prime }_{\mathrm{rf}}$ in the electron rest frame employing the differential scattering cross section (see Sect. 3.4); specifically, we drew θ_rf′ from the cross section and sampled φ′_rf isotropically. We then computed the post-scattering photon energy $E{\prime }_{\mathrm{rf}}$ via energy conservation (Eq. (8) in Loudas et al. 2021). Finally, we applied the inverse Lorentz transformation to return to the laboratory frame and updated the photon-packet properties,

$$\begin{array}{c}\hfill E,\widehat{n},{\mathit{r}}_0,w⟵E^{\prime },{\widehat{n}}^{\prime },{\mathit{r}}_{\mathrm{new}},w_{\mathrm{new}},\end{array}$$(28)

after which the photon was propagated to the next interaction.

3.4. Prescription of resonant Compton scattering

In this work, we focus on the formation of the CRSF, so we need to adopt resonant magnetic scattering cross sections. Fully relativistic QED cross sections have been extensively studied (Harding & Daugherty 1991; Sina 1996; Nobili et al. 2008; Gonthier et al. 2014; Mushtukov et al. 2016), but they turn out to be cumbersome and the primary bottleneck in radiative transfer codes (see, e.g., Araya & Harding 1999; Schwarm et al. 2017; Kumar et al. 2022). To tackle this issue, Loudas et al. (2021) offers approximate, but accurate, resonant magnetic scattering cross sections that are applicable when b ≪ 1 and E ≪ m_ec². Given that the systems of our interest lie in this regime, we adopt the resonant Compton scattering prescription developed by Loudas et al. (2021).

In this study, we neglected polarization effects; therefore, we employed the polarization-averaged differential cross section. In the electron’s rest frame, and integrated over the azimuthal angle φ_rf′, it is expressed as

$$\begin{array}{c}\hfill \frac{d\sigma }{dcos\theta {\prime }_{\mathrm{rf}}}=2\pi \frac{3\pi r_{0}c}{16}L(E_{\mathrm{rf}},E_{\mathrm{res}})(1+{cos}^2{\theta }_{\mathrm{rf}})(1+{cos}^2\theta {\prime }_{\mathrm{rf}}),\end{array}$$(29)

where θ_rf and θ′_rf denote the incident and scattered photon angles relative to the magnetic field. The azimuthal averaging implies that, in our MC implementation, the azimuthal scattering angle is sampled isotropically; this is justified in our regime, where the resonant contribution dominates the spectral formation and is effectively independent of φ′_rf (see, e.g., Ventura 1979). Here, L(E_rf, E_res) is the resonance (Lorentz) line profile,

$$\begin{array}{c}\hfill L(E_{\mathrm{rf}},E_{\mathrm{res}})=\frac{\mathrm{\Gamma }/2\pi }{{(E_{\mathrm{rf}}-E_{\mathrm{res}})}^2/{ħ}^2+{(\mathrm{\Gamma }/2)}^2},\end{array}$$(30)

centered on the relativistic cyclotron energy, E_res (see Eq. (13)). For the relativistic cyclotron transition rate Γ, we used the first-order corrected expression (Loudas et al. 2021)

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{4}{3}\frac{m_{\mathrm{e}}c{e}^2}{{ħ}^2}{b}^2(1-2.7b).\end{array}$$(31)

The magnetic scattering cross section, σ, is the integral of the differential cross section (29) over θ′_rf, namely

$$\begin{array}{c}\hfill \sigma (E_{\mathrm{rf}},{\theta }_{\mathrm{rf}})=2\pi \frac{\pi r_{0}c}{2}L(E_{\mathrm{rf}},E_{\mathrm{res}})(1+{cos}^2{\theta }_{\mathrm{rf}}).\end{array}$$(32)

We employed this expression in the calculation of the effective cross section (17).

3.5. Treatment of the boundary

Although most photons diffuse away from the base, where they are emitted, and escape, a small fraction of them return to the boundary (base). In this work, we do not model the physics of the NS atmosphere located below the base of the column (z < 0). So, when a photon packet reaches the base or when the line of sight intersects the base of the cylinder, a special treatment is required. If the wave packet crosses the base at z = 0, two distinct boundary conditions are considered, “reflection” or “thermalization.”

In the reflection scheme, the packet undergoes geometric reflection, retaining both its energy and weight. This is analogous to a scattering-dominated, stationary atmosphere. The main body of the paper uses this boundary condition. In the thermalization scheme, the wave packet is instantaneously absorbed and reemitted isotropically upward with a new energy sampled from a blackbody distribution (Eq. (15)), while its weight w remains unchanged. We examine this boundary condition in detail in Appendix A.

4. Results

In this section we present the results of our MC simulations for a typical XRP in the subcritical regime, employing the reflection scheme (for a comparison with results obtained using the thermalization scheme, we refer the reader to Appendix A). We assumed the seed photon distribution to be isotropic in the upward hemisphere and uniformly emitted from the base, described by a blackbody spectrum with temperature T given by Eq. (8). We adopted a fiducial value of L^* = 10³⁷ erg s⁻¹ for the critical luminosity; the remaining fixed, dimensionless parameters of the model are listed in Table 1. The local cyclotron energy corresponding to the chosen magnetic field strength at h = 0 was E_c = 25.55 keV. The only remaining free parameter was the relative accretion luminosity, L/L^*, which we varied over the range (0, 1) in order to study its impact on the emergent spectral features.

In each simulation run, we injected 10⁶ seed photon packets of equal initial weight w = 10⁶ to ensure high-fidelity spectra. Throughout this work, we present the resulting spectra in the form of dN/dE as a function of the photon energy E and, where indicated, polar angle θ. We normalized the angle-averaged spectra to unity (i.e., ∫(dN/dE) dE = 1) and scaled the angle-resolved ones according to their relative contribution to the total photon flux.

4.1. Spectral formation in the accretion channel

To illustrate the primary spectral features and the underlying physical picture, we first performed a simulation for a representative luminosity of L = 0.6 L^*. The corresponding emergent spectrum is shown in Fig. 2. As it can be seen there, the resulting spectrum exhibits the expected characteristics outlined in Sect. 2.2. The total emergent spectrum closely follows the seed blackbody at energies well below and well above the cyclotron energy, reflecting the low optical depth experienced by photons outside the resonance. The computed spectrum also features a prominent CRSF, resembling an absorption line, followed by a broad blue wing (bump). Photons in resonance with the quantized electrons experience a strongly enhanced opacity and are effectively trapped, undergoing multiple Compton scatterings and gaining on average energy in the process through Doppler boosting until their energies shift out of resonance, resulting in the formation of the CRSF. This energy redistribution of the scattered photons naturally gives rise to the pronounced blue wing (for a comparison with the case where the thermalization scheme is implemented, see Fig. A.1; there, in contrast, the blue wing is suppressed, with the red wing dominating).

Fig. 2. Representative emergent spectrum for an XRP with a fiducial luminosity L = 0.6L*. The solid black line represents the seed blackbody spectrum with temperature kBT = 5.9 keV. The solid red line shows the total normalized emergent spectrum. The dashed blue and green lines correspond to the angle-selected spectra observed at polar angles in the ranges (0, π/2) and (π/2, π), respectively. The vertical dashed gray line indicates the cyclotron energy at Ec = 25.55 keV. The spectrum features a redshifted, but prominent CRSF around E ∼ 20 keV followed by a bump. The angle-selected spectra are scaled relative to their contribution to the total photon flux.

The CRSF appears redshifted with respect to the local cyclotron energy. The origin of this shift is related to the Doppler effect between the lab frame and the electron’s bulk rest frame, as described qualitatively in Sect. 2.2. We discuss this effect in more detail in the following paragraphs.

Furthermore, we note that the dominant contribution to the total flux arises from photons escaping at angles in the range (0, π/2), owing to the low optical depth values away from resonance and the reflection (or thermalization) of photons crossing the base of the funnel. Yet, a minor fraction of photons escape toward the NS surface. The spectrum of this latter population does not exhibit a prominent absorption feature at the energies corresponding to the CRSF in the total spectrum, since photons propagating with cos θ < 0 are typically out of resonance at those energies. Instead, an absorption feature appears at higher energies, where photons are in resonance with the electrons, but it is strongly smeared out by the contribution of red- and blue-wing photons propagating in different directions. Escaped photons directed toward the NS will likely be reflected off the surface; however, their contribution to the total spectrum is negligible. We therefore safely neglected this photon population and the associated modeling of surface reflection outside the funnel.

To further illuminate the formation of the CRSF, we investigated the angular dependence of the emergent spectrum, which is relevant for phase-resolved studies of CRSFs. Fig. 3 shows the total spectrum in the range (0, π/2) in black, along with decomposed into five narrow viewing angle bins spectra in various colors. We note that the choice of the number of angular bins was arbitrary and does not affect the conclusions. The angle-dependent spectra were scaled based on their relative contribution in the photon flux. For convenience, the position of the classical cyclotron energy is also marked by a dashed vertical gray line.

Fig. 3. Angle-resolved spectra for a typical XRP source with luminosity L = 0.6L*. The solid black line corresponds to the angle-averaged spectrum in the range (0, π/2). The dashed colored lines show spectra extracted in five angular bins within this range. Each spectrum has been normalized to its angular flux contribution. The vertical gray line marks the local cyclotron energy. With increasing cos θ, the CRSF centroid energy shifts to lower values and the blue wing becomes more pronounced.

Evidently, the CRSF profile and its wings depend strongly on the viewing angle, exhibiting two primary characteristics. First, the centroid energy of the CRSF decreases with increasing cos θ. Specifically, as cos θ approaches zero (corresponding to directions perpendicular to the z-axis), the redshift relative to E_c becomes smaller. This behavior arises from the Doppler effect discussed in Sect. 2.2. A photon propagating at an angle θ with respect to the z-axis is in resonance with the infalling electrons when its energy is approximately E_c/γ(1 + β cos θ) (see Eq. (14)); thus, larger values of cos θ correspond to lower resonant energies, naturally producing the trend shown in Fig. 3.

Second, the blue wing becomes increasingly prominent and shifts to higher energies as cos θ increases. This behavior is a direct consequence of the implementation of the reflection scheme: Seed photons undergoing resonant scattering acquire their maximum energies when they are backscattered and therefore Doppler beamed toward the stellar surface (i.e., at angles near cos θ ≈ −1). Upon reflection from the base of the funnel, these high-energy photons reemerge at viewing angles near cos θ ≈ 1, manifesting as prominent blue wings.

Finally, we find that the angle-averaged spectrum features a broader, asymmetric, but shallower CRSF compared to those in angle-resolved spectra. The superposition of CRSFs formed at different viewing angles–each centered at a slightly different energy–produces a wider and shallower feature due to the partial overlap of individual absorption profiles with the wings of neighboring lines. At low luminosity, this is very prominent (see Appendix B).

4.2. Variability due to the accretion luminosity

Having studied the angular dependence of the CRSF properties and the physical origin of the associated spectral features, we now turn to their dependence on the accretion luminosity. At higher luminosities, the braking of the infalling plasma becomes more efficient; consequently, the Doppler effect, which causes the cyclotron energy to appear redshifted, is reduced. As a result, we expect the centroid energy to lie closer to the local value E_c at higher luminosities (see also Sect. 2.2).

We performed a total of 11 simulations, varying the accretion luminosity while keeping all other model parameters fixed. Figure 4 displays the angle-averaged emergent spectra for different luminosities, with the inset highlighting the energy range around the CRSF. All spectra are normalized to unity. As anticipated, the CRSF centroid energy shifts systematically toward higher energies with increasing luminosity, establishing a positive correlation between L and E_CRSF (for a quantitative characterization, see Sect. 4.3), in agreement with the predictions of Mushtukov et al. (2015b).

Fig. 4. Angle-averaged emergent spectra for luminosities spanning the range L/L* ∈ (0, 1). The color gradient from dark blue to yellow corresponds to increasing luminosity. The inset highlights the spectral region around the CRSF. The vertical gray dashed line denotes the local cyclotron energy. The positive correlation between L and the CRSF energy ECRSF is evident.

In addition, the blue wing becomes increasingly pronounced and the CRSF deepens with increasing luminosity. The overall spectrum also hardens. Higher luminosities correspond to larger plasma number densities (and thus optical depth), as Eq. (7) suggests, due to (i) an increased accretion rate and (ii) enhanced compression associated with a reduced terminal velocity. The resulting increase in optical depth naturally accounts for the deepening of the CRSF, the strengthening of the blue wing, and the observed spectral hardening.

4.3. Phase-dependent variability of the line features

Having demonstrated the strong variability of the CRSF with luminosity (Sect. 4.2) and viewing angle (Sect. 4.1), we now quantify the dependence of the line centroid energy E_CRSF and width σ_CRSF on these two parameters. Figure 5a displays the CRSF centroid energy as a function of luminosity for five angular intervals spanning the range (0, π/2), together with the values measured from the angle-averaged spectra in black. We defined E_CRSF as the centroid of the CRSF after subtracting the underlying continuum (see Appendix C). A clear positive correlation between E_CRSF and luminosity is present for all viewing angles, as well as for angle-averaged spectra.

Fig. 5. Dependence of the CRSF centroid energy ECRSF (panel a) and width σCRSF (panel b) on luminosity and viewing angle θ (measured with respect to the magnetic axis). Results are shown for five angular bins spanning the range θ ∈ (0, π/2), with Ec = 25.55 keV. Solid black lines denote quantities extracted from angle-averaged spectra. Both ECRSF and σCRSF exhibit a positive correlation with luminosity for all viewing angles. At a fixed luminosity, ECRSF (σCRSF) decreases (increases) toward smaller viewing angles, indicating substantial variability over the pulse cycle.

At a fixed luminosity, the centroid energy varies systematically with viewing angle: smaller angles relative to the z-axis (i.e., face-on views of the hotspot) consistently yield lower values of E_CRSF, indicating an anticorrelation with cos θ. An explanation is given in Sect. 4.1. The amplitude of this angular variation is ∼30% at low luminosities and decreases as the luminosity increases. From Eq. (14), the resonant energy range of photons moving away from the star spans [E_CRSF(θ = 0),E_CRSF(θ = π/2)], that is

$$\begin{array}{c}\hfill \frac{{\mathrm{\Delta }E|}_{\mathrm{CRSF}}}{E_{\mathrm{res}}}\approx \frac{1}{\gamma }-\frac{1}{\gamma (1+\beta )}={[\gamma (1+1/\beta )]}^{-1}.\end{array}$$(33)

Substituting β with the terminal velocity β₀, we get ΔE|_CRSF/E_res ≈ 29% for β₀ = 0.5, which corresponds to L ≪ L^*. At higher luminosities, β₀ decreases (Eq. (5)), leading to a smaller ΔE|_CRSF, in agreement with the simulation results.

Figure 5b is analogous to Fig. 5a, but shows the CRSF width. While several definitions of line width exist, depending on the assumed profile, here we characterized σ_CRSF as the full width at half minimum (FWHM) of the CRSF after subtracting the underlying continuum; for details we refer the reader to Appendix C. We find that the line width increases with luminosity for all viewing angles, although its absolute value depends substantially on θ. The widths corresponding to the angle-averaged spectra (black curve) are systematically higher than those in angle-dependent spectra, as the width in the former case is dominated by the superposition of CRSFs forming at different angles (see Figs. 3 and B.1). In addition, we find that σ_CRSF is positively correlated with cos θ, in agreement with the predictions of Mushtukov et al. (2015b).

This behavior can be easily understood qualitatively as follows. Higher luminosities produce stronger radiative drag, resulting in a steeper velocity gradient in the accretion flow (see Eq. (3)). Photons traversing the funnel therefore sample a broader range of electron bulk velocities, leading to increased Doppler broadening of the CRSF and a positive correlation between σ_CRSF and L. For a given velocity profile (i.e., fixed luminosity), the magnitude of this effect depends on the propagation direction relative to the flow (in this case, the z-axis): smaller values of cos θ suppress Doppler shifts (see Eq. (14)), while the maximum broadening occurs for photons propagating nearly parallel to the magnetic (flow) axis, cos θ ≈ 1.

Taken together, these results imply that both the CRSF centroid energy and width increase with luminosity, such that higher values of E_CRSF are accompanied by broader lines. Consequently, in studies of phase-averaged spectral evolution, this model predicts a correlated increase of the CRSF energy and width with luminosity. Conversely, the trend reverses for short-timescale, phase-resolved variability. As shown above, E_CRSF anticorrelates with cos θ, whereas σ_CRSF correlates positively with cos θ. For a given luminosity, therefore, the model yields an anticorrelation between the CRSF energy and width. Over a pulse cycle, the line energy (width) reaches its minimum (maximum) value when the accretion funnel is viewed face-on, i.e., along the magnetic-field direction, and its maximum (minimum) when it is seen edge-on.

5. Comparison with data: An application to the X-ray pulsar GX 304−1

Following Mushtukov et al. (2015b), we tested our model by applying it to the observed CRSF variability in the X-ray pulsar GX 304−1, a source exhibiting well-established E_CRSF − L and σ_CRSF − L positive correlations (Klochkov et al. 2012). We adopted the CRSF centroid energy and width measurements reported by Mushtukov et al. (2015b), who analyzed INTEGRAL (Winkler et al. 2003) observations of GX 304−1 obtained during the 2012 type-II outburst (Yamamoto et al. 2012).

To facilitate a direct comparison with the data, we performed an additional set of MC simulations employing the reflection scheme, as described in Sect. 4, but, following Mushtukov et al. (2015b), we took the classical cyclotron energy to be $E_{\mathrm{c}}^{\mathrm{GX}}=66.5\mathrm{keV}$, corresponding to a dimensionless magnetic field strength of b ≈ 0.13, and adopted a critical luminosity of L_GX^* = 2.7 × 10³⁷ erg s⁻¹. The remaining dimensionless NS parameters used in this exercise are listed in Table 2. The photon injection protocol was identical to that utilized in runs shown in Sect. 4.

Figures 6a and b show the GX 304−1 data (red and black circles) together with our simulation results (colored lines) for E_CRSF (left panel) and σ_CRSF (right panel) as functions of luminosity, for three representative narrow viewing-angle intervals. We measured E_CRSF and σ_CRSF from the simulated spectra as described in Sect. 4.3. Regarding the data, black circles correspond to CRSF parameters inferred using a Gaussian absorption profile (GABS), while red circles denote results obtained using a Lorentzian profile (CYCLABS). Further details concerning the data analysis can be found in Sect. 3 of Mushtukov et al. (2015b).

Fig. 6. Comparison of observed cyclotron line variability in GX 304-1 with our model. Panel a (b) shows the CRSF centroid energy (width). Red and black circles represent the observed ECRSF and σCRSF values of GX 304-1, obtained using the CYCLABS model (Lorentzian profile) and the GABS model (Gaussian profile), respectively, to fit the line feature (Mushtukov et al. 2015b). Green, blue, and purple lines correspond to the results derived from our MC simulations for three different viewing angle intervals. The model predictions, for specific viewing angles, agree sufficiently well with the data over an order of magnitude in luminosity.

The model successfully reproduces the observed trends and, for fixed viewing angles, provides good agreement with the data over nearly an order of magnitude in luminosity. In particular, viewing angles corresponding to cos θ ∈ (0.0 − 0.1) (purple curves) simultaneously match both the CRSF centroid energy and width, whereas larger values of cos θ are inconsistent with one or both observables. Since the viewing angle is primarily governed by the system geometry (neglecting gravitational light bending, which is not included here), this result favors configurations in which the accretion funnel (and thus the hotspot) is viewed predominantly edge-on over the pulse cycle. Such configurations may correspond either to a near-orthogonal rotator (α ∼ 70° −90°) observed at low inclination (i ∼ 0° −20°), or to a near-aligned rotator (α ∼ 0° −20°) observed at high inclination (i ∼ 70° −90°). Here, i denotes the angle between the line of sight and the rotation axis, and α is the inclination of the magnetic axis relative to the rotation axis. The latter scenario is the one proposed by Mushtukov et al. (2015b), based on phase-averaged modeling of the CRSF variability, and is consistent with our findings.

Nevertheless, it would be premature to use this comparison to place quantitative constraints on the physical parameters of GX 304−1. Our analysis does not include phase averaging over the pulse cycle, nor does it account for general-relativistic light bending. Additional limitations arise from parameter degeneracies inherent to our modeling approach, including the assumed values of $E_{\mathrm{c}}^{\mathrm{GX}}$, L_GX^*, and other NS parameters. For instance, adopting a higher value of $E_{\mathrm{c}}^{\mathrm{GX}}$ would favor larger values of cos θ in order to maximize the Doppler redshift, potentially altering the inferred geometry. Likewise, increasing β_ff would steepen the predicted E_CRSF − L correlation. A comprehensive comparison with the data therefore requires treating all relevant model parameters as free and performing full phase-averaged calculations, which we defer to future work.

6. Discussion and conclusions

In this work, we investigated quantitatively the origin of the CRSF variability observed in subcritical (L ≲ L^*) XRPs by means of resonant scattering in the accretion funnel above the hotspot, following the model proposed by Mushtukov et al. (2015b). We focused on two key spectral properties: the positive correlations of the CRSF centroid energy, E_CRSF, and line width, σ_CRSF, with the X-ray luminosity, L.

We developed a relativistic MC code, based on the framework described by Loudas et al. (2024b), to perform detailed radiative-transfer calculations in the accretion funnel and derive emergent, angle-dependent spectra. For simplicity, we assumed a semi-infinite cylindrical accretion channel above the hotspot, centered on the magnetic axis. Seed blackbody radiation was injected isotropically upward from the hotspot. The accretion flow was taken to be cold, while the plasma density and velocity profiles were computed analytically, capturing the effects of radiation pressure.

The primary model parameters governing our calculations were the critical luminosity L^* and the surface magnetic field strength B. Other free parameters – such as the NS mass, NS radius, and free-fall speed – were fixed to fiducial values. We varied the accretion luminosity L over a wide range to study its impact on the CRSF properties. As a proof-of-concept, we applied our model to observational data of the X-ray pulsar GX 304−1.

Our main conclusions can be summarized as follows:

1. Resonant scattering of hotspot radiation by the accretion flow imprints a prominent CRSF, followed by a broad blue wing (Fig. 1). The CRSF is always redshifted relative to the local cyclotron energy, E_c, due to the Doppler effect; the magnitude of the redshift decreases both at higher luminosities (Fig. 4) and at larger viewing angles relative to the magnetic axis (Fig. 3).

2. The CRSF centroid energy and width are positively correlated with luminosity for all viewing angles (Fig. 5), reflecting changes in the accretion-flow velocity profile, although their absolute values vary across viewing angles. With increasing luminosity, the blue wing becomes increasingly pronounced, the CRSF deepens, and the overall spectrum hardens due to the increase in the optical depth (Fig. 4).

3. The CRSF properties exhibit strong variability over the pulse cycle. At a fixed luminosity, E_CRSF (σ_CRSF) decreases (increases) with cos θ; the centroid energy reaches its maximum when the funnel is viewed edge-on due to the suppression of the Doppler shift. In contrast, the line width is maximized when the funnel is viewed face-on, as photons sample a broader range of bulk velocities. Thus, in phase-resolved studies, E_CRSF anticorrelates with σ_CRSF during the pulse cycle (Fig. 5).

4. Angle-averaged spectra feature broader, asymmetric, but shallower CRSFs compared to angle-resolved spectra (i.e., at a given pulse phase), owing to the superposition of CRSFs formed at different viewing angles (Fig. 3).

5. When applied to the observational data of the X-ray pulsar GX 304−1, the model successfully reproduces the observed variability of the CRSF centroid energy and width over nearly an order of magnitude in luminosity (Fig. 6). In particular, configurations in which the accretion channel is predominantly viewed edge-on are favored, in line with Mushtukov et al. (2015b), providing tentative constraints on the viewing geometry of the system (see Sect. 5).

Our findings are consistent with the qualitative predictions of Mushtukov et al. (2015b), supporting the hypothesis that the accretion channel constitutes the line-forming region and that the flow velocity profile drives the observed variability. However, unlike those predictions, which rely solely on the behavior of the resonant optical depth, our calculations self-consistently incorporate key radiative-transfer effects, such as photon energy redistribution during scattering, which is responsible for the formation of the CRSF wings.

While the model captures the essential physics governing CRSF variability in subcritical XRPs and successfully reproduces well-established observational trends, it relies on a number of simplifying assumptions that introduce limitations and motivate future extensions. First, we assumed a simplified cylindrical, filled accretion channel; however, more realistic magnetospheric flows point to bow-shaped or annular footprints (see, e.g., Zhu 2025, and references therein), which, nevertheless, we do not expect to alter our conclusions. Second, the funnel was extended to infinity, neglecting the global magnetospheric accretion picture, and assuming a constant magnetic field strength. Neither assumption is physical; however, as we argued in Sect. 2.2, the accretion flow is optically thick to resonant scattering; therefore, the line-forming region lies above the hotspot and is relatively compact. Considering a dipole field (B(r)∝r⁻³) and a vertical size of the line-forming region comparable to the hotspot size (Eq. (4)), the expected variation due to the magnetic field gradient is ≲(3 − 6)%.

We employed analytical, one-dimensional velocity and density profiles; however, RMHD simulations (Markozov et al. 2023) reveal more complex profiles, as the effect of radiation pressure force is less effective near the edges of the funnel, where radiation can easily escape sideways. This subtlety might affect the exact CRSF shape and slope of predicted correlations, but it would not change our conclusions. In addition, we treated the accretion flow as cold (i.e., k_BT_e ≪ E_c). We expect the electron’s temperature to modify the CRSF shape through thermal broadening. So, this approximation is valid, so long as the Doppler effect and the velocity gradient govern the CRSF width. In Sect. 4, we showed that for a typical cyclotron line energy of E_c ≈ 25 keV, σ_CRSF is always greater than 2 keV. Therefore, for electron temperatures of k_BT_e ≲ 2 keV, thermal effects will likely not be important.

Our model does not include a detailed description of the NS atmosphere beneath the accretion funnel (i.e., z < 0) or the associated physical processes. To assess the sensitivity of our results to the treatment of the lower boundary, we considered two limiting cases: one in which photons reaching the base undergo geometric reflection (reflection scheme), conserving their energy, and another in which photons get absorbed and reemitted as blackbody radiation (thermalization scheme). We find that our primary conclusions regarding the CRSF variability are robust against the choice of boundary condition. Additionally, as discussed in Sect. 4.1, we neglected the reflection of downward-escaping photons off the NS surface outside the accretion funnel, as their contribution to the total emergent spectrum is negligible. Nevertheless, a more physically motivated treatment that self-consistently accounts for the structure, plasma dynamics, and radiative transfer in the NS atmosphere – as well as the surrounding stellar surface – would further strengthen the predictive power of the model. Such an extension is beyond the scope of the present work and is left for future explorations.

Another simplification in our calculations concerns the treatment of the hotspot emission. We assumed isotropic, upward blackbody emission, uniformly from the base of the accretion funnel. In contrast, more detailed hotspot emission models predict anisotropic, Comptonized spectra (see, e.g., Meszaros & Nagel 1985; Sokolova-Lapa et al. 2021). To assess the sensitivity of our results to this assumption, we tested alternative emission prescriptions, including pencil-like patterns and point-source emission centered at the origin, (0, 0, 0). In all cases, the resulting variations in the CRSF centroid energy and width were at ≲10% level. This indicates that the CRSF variability explored in this work is not strongly sensitive to the assumed hotspot emission pattern, although we expect such effects to play a critical role in the modeling of pulse profiles, which we will explore in future work.

In this work, we performed polarization-independent calculations. With the current (IXPE Weisskopf et al. 2022; XL-Calibur Abarr et al. 2021; XPoSat¹) and upcoming (eXTP Zhang et al. 2025b) X-ray polarimetry missions, polarization data are becoming increasingly available. X-ray polarization holds strong constraining power regarding the system’s geometry (Poutanen et al. 2024); thus, incorporating polarization into our radiative transfer scheme is a vital next step. We plan to address this topic in detail in our forthcoming work.

Finally, we neglected general relativistic effects in this analysis, which can be important for pulse-profile modeling and thorough comparisons with observations (see, e.g., Markozov & Mushtukov 2024). Such effects may influence our results through: (i) gravitational redshift of the CRSF energy and width, as well as the inferred luminosity, and (ii) gravitational lensing by the NS. While gravitational redshift may shift the absolute values of the CRSF parameters, we expect this shift to be nearly identical in all cases, since the extent of the line-forming region above the hotspot is small. Thus, it should not affect the variability of the CRSF dramatically. Similarly, gravitational light bending will modify the observed range of viewing angles and introduce contributions from the antipodal hotspot. We expect these effects to become important when constructing the predicted pulse profiles.

An alternative interpretation of CRSF variability in subcritical XRPs is provided by the CS scenario, which identifies the line formation site as the region between the hotspot and the shock front. In this framework, the observed variations in the CRSF energy and width are attributed to changes in the local magnetic field strength at the line-forming region, which are tied to the shock altitude. When applied to the GX 304−1 system, this model successfully reproduces the luminosity dependence of the CRSF properties in phase-averaged spectra (Rothschild et al. 2017). Because bulk motion effects in the post-shock region are negligible, the CS model does not predict pulse-resolved CRSFs to exhibit strong variability with the pulse phase (i.e., viewing angle), in contrast to the scenario explored in this work. This fundamental difference between the two models offers a potential observational discriminant that could help distinguish between competing interpretations and advance our understanding of cyclotron line formation in subcritical XRPs. However, such a test requires the development of a robust model for the shock structure, followed by self-consistent radiative-transfer calculations.

In summary, building on the framework proposed by Mushtukov et al. (2015b), we have demonstrated a plausible physical origin of key observational trends in subcritical X-ray pulsars, including the luminosity- and phase-dependent variability of cyclotron lines, by means of resonant scattering in the accretion funnel. Our results highlight the central role of bulk motion in shaping the CRSF properties and provide new insights into the dynamics and geometry of the accretion flow.

Data availability

The results presented in this paper were obtained using a proprietary code developed by the authors. The data supporting the figures, as well as the code, are available from the corresponding author upon reasonable request.

Acknowledgments

We thank the referee for their careful reading of the manuscript and for their constructive comments that helped improve the clarity of this work. PF thanks Chryssi Koukouraki and Orestis Zoumpoulakis (Institute of Astrophysics, Foundation for Research and Technology-Hellas) for their technical assistance. NL acknowledges support from the Venture Forward Fellowship of the Department of Astrophysical Sciences at Princeton University. NL thanks Lizhong Zhang for valuable discussions on accretion column dynamics, and Anatoly Spitkovsky for insightful discussions on the underlying physics of the collisionless shock scenario as an accretion-halting mechanism above the hotspot.

References

Appendix A: Boundary effects

In this appendix, we investigate the effect of the boundary condition at the base of the accretion column on the CRSF properties. To do so, we repeated the calculations presented in Sect. 4, using the thermalization scheme instead of the reflection scheme. Specifically, photons hitting the base of the accretion column are not reflected back; rather, they are absorbed and reemitted as seed photons, following the initial blackbody distribution. We expect this change to modify the shape of the line wings, but not the centroid position, as the scattering mechanism remains the same. In particular, we expect an enhancement of the red wing, as each time a photon packet reaches the base, its energy is resampled and k_BT < E_c. This process effectively removes photons from the blue wing and redistributes them to the left of the absorption line.

In Fig. A.1, we compare the angle-averaged spectra obtained using the two different boundary conditions for luminosity L = 0.6L^*. It is evident that when the thermalization scheme is implemented (blue line), the red wing dominates over the blue one, whereas for the reflection scheme (red line), the opposite is observed. However, the position of the absorption line is the same in both cases; the centroid energy of the CRSF is insensitive to the boundary condition.

Fig. A.1. Comparison of angle-averaged spectra with different boundary conditions for a typical source with luminosity L = 0.6L* and Ec = 25.55 keV (vertical gray dashed line). The solid black line represents the seed blackbody spectrum with kBT = 5.9 keV. The solid red (blue) line shows the emergent spectrum for the reflection (thermalization) case. The line wings exhibit differences; the red (blue) wing dominates in the thermalization (reflection) scheme. The centroid energy of the CRSF appears to be independent of the boundary condition.

In analogy with Fig. 4, Fig. A.2 shows the resulting spectra for 11 different luminosities L. The positive correlation between the position of the CRSF and luminosity is again evident. The thermalization scheme suppresses the blue wing and enhances the red wing, as described qualitatively above. Even at high luminosities, the blue wing remains relatively weak.

Fig. A.2. Same as Fig. 4, but for the thermalization scheme. The luminosity dependence of the CRSF properties is the same.

Figure A.3 is the same as Fig. 5, but also includes the results obtained using the thermalization scheme in triangles connected with dashed lines. From Fig. A.3a, it is evident that the centroid energy is not affected by the scheme used, regardless of the luminosity and/or viewing angle.

Fig. A.3. Same as Fig. 5, but including the results obtained using the thermalization scheme in triangles connected with dashed lines. While the energy of the cyclotron feature is independent of the boundary scheme, the line width shows deviations at high luminosities. We expect this, as the redistribution of photons into the line wings alters the CRSF shape, affecting its width.

In contrast, the line width is sensitive to the overall shape of the absorption feature, including its wings; consequently, we expect differences between the two schemes. This is apparent in Fig. A.3b, where the results for the two schemes diverge at luminosities L ≳ 0.5L^*. The discrepancy between the dashed (thermalization) and solid (reflection) lines increases with luminosity and is most pronounced at small values of cos θ. The primary cause of this divergence is the effective narrowing of the CRSF in the reflection scheme due to its prominent blue wing, which is absent in the thermalization case (see Figs. A.1 and A.2). We note, however, that the quantitative determination of the line width is subject to systematic uncertainties. The derived values depend heavily on the modeling of the underlying continuum and the choice of line profile (e.g., Gaussian or Lorentzian).

Appendix B: Absorption feature at low luminosity

This appendix aims at explaining the lack of a prominent absorption feature in the spectrum at low luminosities in our simulations. To this end, in Fig. B.1 we present the angle-resolved spectrum for a low luminosity value, L = 0.1L^*, calculated using the reflection scheme. All features seen in Fig. 3 are present here as well: the centroid energy is anti-correlated with cos θ, and the blue wing becomes increasingly prominent and shifts to higher energies as cos θ increases. However, in contrast to the L = 0.6L^* case, the superposition of lines emerging at different viewing angles in the low luminosity regime results in a shallow, almost flat CRSF, making it a challenging task to accurately identify the angle-averaged CRSF centroid. Nevertheless, the feature is readily distinguishable after subtraction of the underlying continuum, enabling a reliable determination of the centroid energy.

Fig. B.1. Same as Fig. 3, but for luminosity L = 0.1L*.

Appendix C: Calculation of line parameters

This appendix outlines the methodology used to measure the relevant CRSF parameters: the centroid energy (E_CRSF) and the line width (σ_CRSF). The CRSF shape is rather complicated, asymmetric, and deviates significantly from a Gaussian or a Lorentzian profile; thus, we adopted an approach that allows for a robust metric without any bias introduced by fitting specific analytical profiles. We defined the centroid energy as the flux-weighted average energy in the vicinity of the local minimum, and the line width as the full width at half minimum (FWHM) of the absorption feature’s depth.

Figure C.1 illustrates our approach, showing its implementation in a representative angle-averaged spectrum (cos θ ∈ (0, 1)) for a source with luminosity L = 0.6L^*, computed using the thermalization scheme. To estimate the line parameters, we first subtracted the background blackbody continuum (S_BB; black dashed line) from the total emergent spectrum (S_tot; solid gray line). Since the resulting residual spectrum (S_tot − S_BB), shown in solid red, does not saturate to zero away from the CRSF, we further fitted the remaining continuum with a sixth-order polynomial fit (C_poly; dotted purple curve), excluding the energy range (pink shaded rectangular), which is dominated by the absorption feature (i.e., 10 keV − 40 keV). We then subtracted it from the residual spectrum to obtain the isolated CRSF (S_tot − S_BB − C_poly; solid green curve). Finally, we measured the line width by computing the FWHM (blue line) and the centroid energy by calculating the flux-weighted average energy within the energy interval around the minimum where the absorption depth exceeds 90% of the maximum depth (indicated by the red segment in Fig. C.1).

Fig. C.1. Illustration of the procedure used to measure the line parameters for a run with L = 0.6L*. Stot is the total spectrum in gray, SBB is the seed blackbody continuum in dashed black, Stot − SBB is the residual spectrum in red, and Cpoly is the sixth-order polynomial fit of the residual continuum spectrum, excluding the interval around the CRSF (pink shaded region). We measured the line width from the isolated absorption feature (Stot − SBB − Cpoly); green curve) as the FWHM (blue segment). Likewise, we estimated the centroid energy as the flux-weighted average energy within the region around the minimum where the absorption depth is larger than 90% of the maximum depth (red segment).

We note that the choice of the excluded energy range in the polynomial fit can influence the final measurement; therefore, care must be taken in modeling the continuum to ensure stability. We varied the upper and lower limits of the exclusion zone by 10% to test the sensitivity of our inferences; the derived width and centroid energy change by less than 3%. We therefore consider the measurements robust, although a more sophisticated approach would further strengthen their validity.

All Tables

All Figures

Fig. 1. Illustration of the MC code’s geometry, coordinates, and calculation of the line-of-sight optical depth. A photon packet located at r0 is moving in direction $\widehat{n}$. We computed the accumulated optical depth by dividing the line of sight into a sequence of points r0 using an adaptive step-size algorithm. We then sampled a target optical depth, τrnd, and determined the position of the next scattering event via linear interpolation between the discretized points. We assumed a uniform magnetic field because the height of the line-forming region is much smaller than the NS radius.

In the text

Fig. 2. Representative emergent spectrum for an XRP with a fiducial luminosity L = 0.6L*. The solid black line represents the seed blackbody spectrum with temperature kBT = 5.9 keV. The solid red line shows the total normalized emergent spectrum. The dashed blue and green lines correspond to the angle-selected spectra observed at polar angles in the ranges (0, π/2) and (π/2, π), respectively. The vertical dashed gray line indicates the cyclotron energy at Ec = 25.55 keV. The spectrum features a redshifted, but prominent CRSF around E ∼ 20 keV followed by a bump. The angle-selected spectra are scaled relative to their contribution to the total photon flux.

In the text

Fig. 3. Angle-resolved spectra for a typical XRP source with luminosity L = 0.6L*. The solid black line corresponds to the angle-averaged spectrum in the range (0, π/2). The dashed colored lines show spectra extracted in five angular bins within this range. Each spectrum has been normalized to its angular flux contribution. The vertical gray line marks the local cyclotron energy. With increasing cos θ, the CRSF centroid energy shifts to lower values and the blue wing becomes more pronounced.

In the text

	Fig. 4. Angle-averaged emergent spectra for luminosities spanning the range L/L* ∈ (0, 1). The color gradient from dark blue to yellow corresponds to increasing luminosity. The inset highlights the spectral region around the CRSF. The vertical gray dashed line denotes the local cyclotron energy. The positive correlation between L and the CRSF energy ECRSF is evident.
In the text

Fig. 5. Dependence of the CRSF centroid energy ECRSF (panel a) and width σCRSF (panel b) on luminosity and viewing angle θ (measured with respect to the magnetic axis). Results are shown for five angular bins spanning the range θ ∈ (0, π/2), with Ec = 25.55 keV. Solid black lines denote quantities extracted from angle-averaged spectra. Both ECRSF and σCRSF exhibit a positive correlation with luminosity for all viewing angles. At a fixed luminosity, ECRSF (σCRSF) decreases (increases) toward smaller viewing angles, indicating substantial variability over the pulse cycle.

In the text

Fig. 6. Comparison of observed cyclotron line variability in GX 304-1 with our model. Panel a (b) shows the CRSF centroid energy (width). Red and black circles represent the observed ECRSF and σCRSF values of GX 304-1, obtained using the CYCLABS model (Lorentzian profile) and the GABS model (Gaussian profile), respectively, to fit the line feature (Mushtukov et al. 2015b). Green, blue, and purple lines correspond to the results derived from our MC simulations for three different viewing angle intervals. The model predictions, for specific viewing angles, agree sufficiently well with the data over an order of magnitude in luminosity.

In the text

Fig. A.1. Comparison of angle-averaged spectra with different boundary conditions for a typical source with luminosity L = 0.6L* and Ec = 25.55 keV (vertical gray dashed line). The solid black line represents the seed blackbody spectrum with kBT = 5.9 keV. The solid red (blue) line shows the emergent spectrum for the reflection (thermalization) case. The line wings exhibit differences; the red (blue) wing dominates in the thermalization (reflection) scheme. The centroid energy of the CRSF appears to be independent of the boundary condition.

In the text

	Fig. A.2. Same as Fig. 4, but for the thermalization scheme. The luminosity dependence of the CRSF properties is the same.
In the text

	Fig. A.3. Same as Fig. 5, but including the results obtained using the thermalization scheme in triangles connected with dashed lines. While the energy of the cyclotron feature is independent of the boundary scheme, the line width shows deviations at high luminosities. We expect this, as the redistribution of photons into the line wings alters the CRSF shape, affecting its width.
In the text

	Fig. B.1. Same as Fig. 3, but for luminosity L = 0.1L*.
In the text

Fig. C.1. Illustration of the procedure used to measure the line parameters for a run with L = 0.6L*. Stot is the total spectrum in gray, SBB is the seed blackbody continuum in dashed black, Stot − SBB is the residual spectrum in red, and Cpoly is the sixth-order polynomial fit of the residual continuum spectrum, excluding the interval around the CRSF (pink shaded region). We measured the line width from the isolated absorption feature (Stot − SBB − Cpoly); green curve) as the FWHM (blue segment). Likewise, we estimated the centroid energy as the flux-weighted average energy within the region around the minimum where the absorption depth is larger than 90% of the maximum depth (red segment).

In the text