© The Authors 2024

1. Introduction

Accretion-powered X-ray Pulsars (XRPs) are strongly magnetized neutron stars (NSs) that accrete matter from a donor companion star. Due to the NS’s enormous magnetic field (B ∼ 10¹² G), the inner region of the accretion disk is disrupted at a distance, which is called the Alfven radius, and the plasma becomes attached to the field lines. From there, the accreting matter is channeled onto the magnetic poles of the NS, where most of the gravitational energy is released, thereby resulting in the emergence of X-rays (Basko & Sunyaev 1976). For a recent review, we refer the reader to Mushtukov & Tsygankov (2022).

A quantum-mechanical phenomenon called magnetic resonant scattering gives the distinct opportunity to precisely measure the NS’s magnetic field strength on its surface. The motion of a charged particle (e.g., electron) in the presence of a strong magnetic field is quantized across its axis (Landau 1930; see also Landau & Lifshitz 1981), leading to discrete energy states that are commonly known as Landau levels. Hence, absorption or emission features are expected to appear in the spectra of highly magnetized NSs. These features are the so-called cyclotron lines, also known as cyclotron resonant scattering features (CRSFs). The difference in energy between two successive levels encodes information about the local magnetic field strength existing in the line-forming region, as it is directly proportional to its strength. The detection of such a line in the spectrum of a NS offers a direct and robust way of measuring the magnetic-field strength in NSs (Gnedin & Sunyaev 1974; Basko & Sunyaev 1975). On top of that, cyclotron lines are thought to be a solid diagnostic tool to probe the physics of accretion onto the magnetic poles of highly magnetized NSs (Staubert et al. 2019).

The first ever discovery of a cyclotron line in NSs occurred 46 years ago in the X-ray spectrum of the accreting binary pulsar Hercules X-1 (Trümper et al. 1977; Truemper et al. 1978), revealing this new powerful phenomenon to the X-ray community. Since then, more than 35 electron cyclotron lines have been detected in accreting magnetic NS X-ray spectra, sometimes with their harmonics (for a review see, Staubert et al. 2019) and, even more rarely, with anharmonics (see e.g., Fürst et al. 2018; Sharma et al. 2022; Yang et al. 2023), covering an energy range from ∼10 keV to ∼100 keV.

A large amount of observational data has been collected for most of these sources. Yet, the region where the CRSFs are formed and the mechanism that produces them remain a fundamental puzzle, preventing us from inferring the NS’s magnetic field strength accurately, as it varies noticeably on the NS surface and along the accretion column, as well.

In the standard picture, the formation of cyclotron lines in NSs is thought to take place in the lower part of the accretion column through resonant scattering (Gnedin & Sunyaev 1974; Basko & Sunyaev 1975). The accretion flow heats the magnetic polar cap resulting in a hot spot that emits soft X-ray photons. If the accretion rate is high enough, above a critical value, then the radiation pressure exerted on the in-falling plasma is efficient in decelerating the material well before it hits the surface of the NS, thereby giving rise to a radiative shock (for a detailed discussion, see Mushtukov et al. 2015a). While passing through the shock, the in-falling material slows down transferring most of its gravitational energy to photons via collisions of each electron with a multitude of photons and eventually settling down to the post-shock region, having relatively small flow velocity and mainly thermal one. Kylafis et al. (2014) demonstrated that this mechanism is capable of producing X-ray power-law spectra through a combination of bulk and thermal Comptonization of bremsstrahlung seed photons in the shock. Basko & Sunyaev (1976) suggested that a fraction of the radiation, which has energy comparable to the cyclotron energy, will undergo resonant scattering while criss-crossing the shock, eventually producing a CRSF in the emergent spectrum. Yet, no calculation has been performed so far for the formation of an X-ray power-law spectrum along with a cyclotron line in a radiative shock.

Over the last decades, several extensions to this scenario and other models have been proposed to address the problem of CRSF formation in NSs. They can be divided into two different classes; the first one postulates that the absorption feature emerges in the accretion column (e.g., Becker & Wolff 2007; Becker et al. 2012; Nishimura 2014; Mushtukov et al. 2015b; Becker & Wolff 2022), while the latter considers the reprocessing of X-rays emitted by the accretion column on the NS surface (see, e.g., Poutanen et al. 2013). Nevertheless, none of them have been confirmed to be quantitatively correct yet.

An observation that could be instrumental in constraining the currently viable models is the variation of the cyclotron line energy with the emergent X-ray luminosity. Several CRSFs, reported in XRP spectra, exhibit a correlation between the line energy E_c and the emergent X-ray luminosity L_X, indicating that the configuration of the line-forming region may depend on the accretion rate. It has been observed that XRPs of relatively low luminosity show a positive E_c − L_X correlation (e.g., Her X-1: Staubert et al. 2007; GX 304-1: Klochkov et al. 2012; Vela X-1: Fürst et al. 2014). In contrast, sources of high luminosity often exhibit a negative correlation (e.g., V0332+53: Tsygankov et al. 2006). This trend suggests that whether the correlation is positive or negative is dictated by a critical value of the X-ray luminosity (L_cr ∼ 10³⁷ erg s⁻¹) above which the correlation becomes negative. This phenomenon has already been observed in two XRPs (1A 0535+262: Kong et al. 2021; Shui et al. 2024 and V0332+53: Doroshenko et al. 2017; Vybornov et al. 2018; however, see also the GRO J1008-57: Chen et al. 2021).

Basko & Sunyaev (1976) were the first to show that such a critical value exists for the accretion luminosity above which the formation of a radiative shock takes place, stating that the wide range of observed luminosities in accreting XRPs is likely associated with different accretion regimes. This threshold naturally splits the accretion regimes into two distinct types: 1) accreting NSs of subcritical accretion luminosity with no shock at all (however, see Bykov & Krasil’Shchikov 2004) and 2) XRPs of supercritical luminosity that can support a radiation-dominated shock in the accretion column.

Subcritically accreting NSs have been studied extensively, and a solid model was proposed recently by Mushtukov et al. (2015b). According to this model, the deceleration of the in-falling material occurs on the NS surface, at a height comparable to the size of the hot-spot, via scattering of electrons by the hot-spot’s radiation. This scenario invokes the Doppler effect to explain the positive correlation between the emergent X-ray luminosity L_X and the cyclotron line centroid, E_c, due to the variation of the velocity profile with radiation pressure. It has been successful in describing the data of the XRP GX 304-1. Alternative scenarios include 1) the braking of the plasma solely by Coulomb collisions near the NS atmosphere proposed by Staubert et al. (2007); and 2) the emergence of a strong collisionless shock above the polar cap that dissipates the ion’s kinetic energy (Bisnovatyi-Kogan & Fridman 1970; Shapiro & Salpeter 1975; Langer & Rappaport 1982; Bykov & Krasil’Shchikov 2004; for an implementation of the model to observations see Rothschild et al. 2017; Vybornov et al. 2017). In the Coulomb collisions (collisionless shock) model, the larger the luminosity, the smaller the characteristic ion’s stopping length (shock height) and, therefore, the greater the local magnetic field strength (i.e., cyclotron energy) in the line-forming region. This naturally yields a positive E_c − L_X correlation.

On the other hand, accreting NSs that exhibit a super-critical luminosity are generally described by the model of Basko & Sunyaev (1976). In this scenario, the height of the shock scales approximately linearly with the X-ray luminosity, while the magnetic field strength drops with the increase of the shock’s height. Thus, it predicts an anticorrelation between the line energy and the X-ray luminosity, qualitatively in line with the observations of super-critical XRPs. However, it was pointed out by Poutanen et al. (2013) that in a dipole field, the predicted rate of change of E_c with L_X is much larger than the observed one. As a remedy, the authors put forward the idea that the observed cyclotron line is not produced in the shock, but from the reflection of radiation produced in the shock in the area surrounding the polar cap.

In the recent work of Kylafis et al. (2021), we investigated the cyclotron-line formation scenario proposed by Poutanen et al. (2013). The reprocessing of photons on the NS surface has, in principle, the ability to produce a negative correlation between cyclotron-line-centroid and X-ray luminosity. As the luminosity increases, the shock height also increases, and thus more and more photons are scattered at lower and lower magnetic latitudes, where the magnetic field is lower. The magnetic field at the equator is half that of the pole. We employed a Monte Carlo (MC) code to derive the reflected spectrum, and we demonstrated that this mechanism cannot produce prominent CRSFs in the emergent X-ray spectrum similar to the ones observed in NS spectra. We also showed that the reflection model cannot explain the negative correlation between the cyclotron energy and X-ray luminosity reported in the emergent spectrum of V0332+53 source (Tsygankov et al. 2010). The reason is that scattering preserves the number of photons. Thus, photons that are missing in an absorption feature will appear as an emission feature -a “bump”- at slightly lower energy. Adding the spectra from all illuminated magnetic latitudes washes out the expected correlation.

In this work, we aim to revisit the model of Basko & Sunyaev (1976) and explore the possibility of cyclotron lines being formed in the radiative shock in the accretion column of a NS by means of a MC simulation. We considered an accreting magnetic NS in the supercritical regime with a radiative shock in a cylindrical accretion column. Our input photons were bremsstrahlung ones emitted below the shock. We let the photons criss-cross the shock as many times as they wanted before they escape sideways. We derived the emergent spectrum from the accretion shock taking into account both thermal- (post-shock) and bulk-Comptonization (pre-shock), as well as the resonant Compton scattering of photons by electrons, implementing the prescription developed in Loudas et al. (2021). This allowed for the emergence of a CRSF in the total spectrum. The importance of the Doppler effect between the shock’s reference frame and the bulk-motion’s frame in the CRSF’s energy centroid was also quantified. This model is an extension of the one proposed by Kylafis et al. (2014) to explain the hard X-ray power-law tails observed in the spectra of anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs).

The layout of this paper is as follows. In Sect. 2, we offer a qualitative presentation of our model, putting special emphasis on the physical processes that occur in the radiative shock and offering useful scaling relations among the parameters considered in this study. In Sect. 3, we describe the MC code implemented in this study. In Sect. 4, we exhibit our results, providing a qualitative comparison with observations. In Sect. 5, we comment on the results and draw our conclusions.

2. The model

To demonstrate our underlying ideas, we used a relatively simple model for describing accreting NSs of super-critical luminosity. Here, we considered a highly magnetized NS that is accreting matter at an accretion rate Ṁ. The strong dipole magnetic field channels the plasma to the magnetic poles of the NS. For simplicity, we assumed the accretion column to be a cylinder of cross-sectional area πa² (though, accretion physics dictates that the accretion column is not cylindrical, but bow-like around the magnetic pole; see, e.g., Basko & Sunyaev 1976; Romanova et al. 2004) and a radiative shock at height H above the magnetic pole. Above the shock, matter falls freely at a constant speed (pre-shock region), while below the shock (post-shock region) it has negligible flow velocity, electron temperature T_e, and density seven times the pre-shock one. The free-falling plasma is decelerated by the radiation field while passing through the shock, which is treated here as a mathematical discontinuity. To fit real data, more detailed models are required.

Bremsstrahlung photon beams from below the shock were injected isotropically upward, and each beam was followed by our MC code (see Sect. 3) until its intensity has dropped below 0.1% of the injected value. The crucial parameters in this calculation were the free-fall velocity, the accretion luminosity (i.e., accretion rate), the magnetic field strength at the shock, and the transverse Thomson optical depth, since the longitudinal optical depth in both the pre-shock and the post-shock regions is infinite by construction in our model.

2.1. Qualitative considerations

2.1.1. Accretional energy

Assuming for demonstration purposes that the accreting plasma is composed only of protons and electrons, due to charge and current neutrality, each free-falling proton is accompanied by an electron that travels at the same speed. Thus, all the accretional energy resides in the protons. Ignoring small thermal motion, we can treat all proton-electron pairs as if they are accreting at speed:

$$\begin{array}{c}\hfill v_{\mathrm{ff}}=\sqrt{2GM/r},\end{array}$$(1)

where G is the gravitational constant, M is the NS’s mass, r is the distance from its center, that is, r = R + H, where R is the NS’s radius. In a typical NS, the infalling matter can reach a mildly relativistic speed (e.g., v_ff ∼ 0.5 c). Each pair at r > R + H carries energy, which is equal to

$$\begin{array}{c}\hfill E_{\mathrm{ff}}=E_{\mathrm{p}}+E_{\mathrm{e}}\approx GM{m}_{\mathrm{p}}/r.\end{array}$$(2)

For a typical NS, this energy is on the order of 10² MeV, which is well above the X-ray band.

2.1.2. Deceleration mechanism and shock parameters

At first sight, it seems impossible for the accretion flow’s deceleration to occur in the radiative shock, as the material below the shock is thermal and emits bremsstrahlung radiation, which, in photon number, dominates in the soft X-ray band. For simplicity, we assumed that initially, each proton-electron pair has energy E_ff on the order of 100 MeV, while the X-ray photons have energies E_γ that lie in the (1 − 100) keV range. After a head-on collision of a photon with a relativistic electron, the photon-energy gain is at most ΔE_γ ≈ (a few) × E_γ, which is significantly smaller than the proton-electron pair’s accretional energy (see Eq. (2)). Hence, the pair loses only a tiny fraction of its accretional energy in a scattering. As we demonstrate below, however, the electron’s mean free path is orders of magnitude smaller than the photon’s mean free path. Thus, each electron is scattered by many photons while passing through the shock. This mechanism makes it possible for the accretion flow to transfer most of its energy to photons and slow down. The phenomenon is similar to the slowing down of a diver by the water molecules in the top layers of water in a pool.

We envision the deceleration process as follows: in-falling electrons undergo multiple scatterings with the X-ray photons that try to escape the accretion column upward; hence, they gradually lose energy, while radiation escapes sideways. However, charge neutrality implies that after each scattering event between a photon and an electron, a “companion” proton supplies a small fraction of its energy to the electron in order to counterbalance its loss due to scattering. As a result, the electron acts as a mediator to transfer the proton’s energy to photons while traversing the shock. This process takes place in the accretion shock and is terminated when the protons and the electrons acquire the same temperature behind the shock. Afterwards, the plasma sinks down in the post-shock region, having negligible flow velocity -but only a thermal one- and losing energy via bremsstrahlung. A fraction of the bremsstrahlung photons go upward and create the shock, and in the process the upscattered photons that escape constitute the largest fraction of the observed luminosity, the called accretion luminosity, which equals

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=G\frac{M\dot{M}}{r}=\frac{1}{2}\dot{M}{\beta }_{\mathrm{ff}}^2{c}^2,\end{array}$$(3)

where β_ff is the free-fall plasma velocity v_ff normalized to the speed of light c.

The photons’ energy density u_γ in the shock region is estimated by

$$\begin{array}{c}\hfill u_{\gamma }=\frac{L_{\mathrm{acc}}}{\pi a^{2}c}.\end{array}$$(4)

It is convenient to express u_γ in terms of the mass-accretion rate, as it can then be estimated directly by invoking the continuity equation. Substituting Eqs. (3) into (4), we obtain

$$\begin{array}{c}\hfill u_{\gamma }=\frac{\dot{M}}{\pi a^2}\frac{v_{\mathrm{ff}}^2}{2c}.\end{array}$$(5)

Conservation of mass yields the following approximate formula for the plasma density in the shock region:

$$\begin{array}{c}\hfill \rho =\frac{\dot{M}}{\pi a^2{v}_{\mathrm{ff}}},\end{array}$$(6)

which, after taking into account charge neutrality (implying n_e = n_p) and that m_e ≪ m_p, becomes

$$\begin{array}{c}\hfill n_{\mathrm{e}}=\frac{\dot{M}/m_{\mathrm{p}}}{\pi a^2{v}_{\mathrm{ff}}}.\end{array}$$(7)

Using Eq. (5) along with Eq. (7), we find

$$\begin{array}{c}\hfill \frac{u_{\gamma }}{n_{\mathrm{e}}}={\beta }_{\mathrm{ff}}\frac{1}{2}{m}_{\mathrm{p}}{v}_{\mathrm{ff}}^2.\end{array}$$(8)

Next, writing u_γ as the photons’ number density n_γ times a characteristic photon’s energy E_γ, the ratio n_γ/n_e reads

$$\begin{array}{c}\hfill \frac{n_{\gamma }}{n_{\mathrm{e}}}={\beta }_{\mathrm{ff}}\frac{m_{\mathrm{p}}{v}_{\mathrm{ff}}^2}{2E_{\gamma }}={\beta }_{\mathrm{ff}}\frac{E_{\mathrm{ff}}}{E_{\gamma }},\end{array}$$(9)

where E_ff is given by Eq. (2).

The characteristic mean free path for propagation of electrons in the “bath” of outgoing photons is

$$\begin{array}{c}\hfill {\overline{\lambda }}_{\mathrm{e}}=\frac{1}{n_{\gamma }{\sigma }_{\tau }},\end{array}$$(10)

while the corresponding photon’s mean free path is given by

$$\begin{array}{c}\hfill {\overline{\lambda }}_{\gamma }=\frac{1}{n_{\mathrm{e}}{\sigma }_{\tau }},\end{array}$$(11)

where, for simplicity in these back-of-the-envelope calculations, we employ the classical Thomson cross-section σ_τ for the e − γ interaction. We note that in our MC calculation, we made use of the accurate magnetic resonant Compton cross-section, which accounts for the photon’s direction of propagation, energy, and polarization; and depends on the magnetic field (see Sect. 3.4).

Combining Eqs. (9)–(11), we find the ratio ${\overline{\lambda }}_{\gamma }/{\overline{\lambda }}_{\mathrm{e}}$ to be

$$\begin{array}{c}\hfill \frac{{\overline{\lambda }}_{\gamma }}{{\overline{\lambda }}_{\mathrm{e}}}={\beta }_{\mathrm{ff}}\frac{E_{\mathrm{ff}}}{E_{\gamma }}\lesssim E_{\mathrm{ff}}/E_{\gamma }\sim \mathcal{O}({10}^3-{10}^4),\end{array}$$(12)

where we took into account that β_ff ≲ 1. Equation (12) indicates that the electron’s mean free path is significantly smaller than that of the photon, as the number density of photons is way greater than that of the electrons (see Eq. (9)). Given that a characteristic length scale for the shock’s width is the mean free path of the photons (see Eq. (15)), this result implies that each electron is scattered by many photons in the shock. Since the ratio of the photon’s mean free path to that of the electron is approximately inversely proportional to their corresponding energy ratio, the effective deceleration of the free-falling plasma in the shock arises naturally. This is in line with the statement we made above.

By construction in our model, the accretion column is effectively infinite along the magnetic field’s axis; thus, the longitudinal optical depth for the photon’s propagation is infinite, implying that all photons escape sideways. So, an important quantity for our calculations is the transverse Thomson optical depth τ_⊥, which is given by

$$\begin{array}{c}\hfill {\tau }_{\perp }=n_{\mathrm{e}}{\sigma }_{\tau }a=a/{\overline{\lambda }}_{\gamma }.\end{array}$$(13)

Substituting Eqs. (7) into (13), we obtain a relation between the transverse optical depth and the mass accretion rate, namely

$$\begin{array}{c}\hfill {\tau }_{\perp }=\frac{\dot{M}}{m_{\mathrm{p}}\pi a{v}_{\mathrm{ff}}}{\sigma }_{\tau }.\end{array}$$(14)

Finally, we can infer a characteristic value for the longitudinal width of the radiative shock; that is, the shock is roughly defined as the region in which the in-falling plasma loses most of its energy, transferring it to photons through multiple scatterings. As long as E_ff ≫ E_γ, after each scattering event the proton-electron pair’s energy is reduced by an amount of ΔE ∼ 3E_γ (if β_ff = 0.5), on average. Thus, N ∼ E_ff/E_γ scattering events are required for an electron-proton pair to be effectively decelerated, supplying its energy to X-ray photons. We demonstrated earlier on that the ratio E_ff/E_γ is roughly equal to n_γ/n_e (see Eq. (9)). Hence, if each photon has undergone one scattering by an electron, on average, then each electron will have been scattered ∼n_γ/n_e times, on average; this is hence how the the accretional energy is converted into X-ray radiation for the reasons mentioned above. This simple consideration leads us to define the characteristic shock’s width, ζ, as the longitudinal mean free path of the photons ${\overline{\lambda }}_{\gamma }$ given in Eq. (11), since it represents the statistically average distance that photons travel between two sequential scattering events. Therefore,

$$\begin{array}{c}\hfill \zeta \simeq {\overline{\lambda }}_{\gamma }=\frac{1}{n_{\mathrm{e}}{\sigma }_{\tau }}.\end{array}$$(15)

For the transverse size of the shock across the column, we adopted the work of Kylafis et al. (2014). They offer the following simple formula for the radius of the accretion column at the shock

$$\begin{array}{c}\hfill a\approx 270{\left(\frac{\dot{M}}{{10}^{16}\mathrm{g}{s}^{-1}}\right)}^{1/5}{\left(\frac{B_0}{{10}^{12}\mathrm{G}}\right)}^{-1/4}\mathrm{m},\end{array}$$(16)

where B₀ is the magnetic field strength on the NS surface.

2.2. Convenient expressions

In Sect. 2.1, we derived simple formulae for every parameter that is considered important in studying the radiation that emerges from the radiative shock in the accretion column. It is helpful to derive characteristic values for these parameters. For a typical accreting magnetized NS, we have M = 1.4 M_⊙, r = 10 km (assuming H ≪ r), B₀ = 10¹² G, and Ṁ = 5 × 10¹⁶ g s⁻¹. Thus, Eqs. (1), (3), (14), and (15) are for convenience written as

$$\begin{array}{c}\hfill {\beta }_{\mathrm{ff}}\simeq 0.65{\left(\frac{M}{1.4M_{\odot }}\right)}^{1/2}{\left(\frac{r}{10\mathrm{km}}\right)}^{-1/2},\end{array}$$(17)

$$\begin{array}{c}\hfill L_{\mathrm{acc}}\simeq 2\times {10}^{36}\left(\frac{\dot{M}}{{10}^{16}\mathrm{g}{s}^{-1}}\right)\left(\frac{M}{1.4M_{\odot }}\right){\left(\frac{r}{10\mathrm{km}}\right)}^{-1}\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$(18)

$$\begin{array}{c}\hfill {\tau }_{\perp }\simeq 2.5{\left(\frac{\dot{M}}{{10}^{16}\mathrm{g}{s}^{-1}}\right)}^{4/5}{\left(\frac{M}{1.4M_{\odot }}\right)}^{-1/2}{\left(\frac{r}{10\mathrm{km}}\right)}^{1/2}{\left(\frac{B_0}{{10}^{12}\mathrm{G}}\right)}^{1/4},\end{array}$$(19)

$$\begin{array}{cc}& \zeta \simeq 107{\left(\frac{\dot{M}}{{10}^{16}\mathrm{g}{s}^{-1}}\right)}^{-3/5}{\left(\frac{M}{1.4M_{\odot }}\right)}^{1/2}\hfill \\ \hfill & \times {\left(\frac{r}{10\mathrm{km}}\right)}^{-1/2}{\left(\frac{B_0}{{10}^{12}\mathrm{G}}\right)}^{-1/2}\mathrm{m}.\hfill \end{array}$$(20)

We note that for high-luminosity sources (ṅ ∼ 10¹⁷ g s⁻¹), we can safely consider that ζ ≪ H.

3. The MC code

The MC code that was employed in this study is an extension of the one used in Kylafis et al. (2014), with the use of resonant cross sections (see Sect. 3.4 below). The inner workings of the code are based on the works of Cashwell & Everett (1959) and Pozdnyakov et al. (1983), and the radiative transfer method is the so-called forced first collision, as it allows us to properly treat the abrupt change in the mean free path from the pre-shock region to the post-shock one and results in high accuracy calculations, even for optically thin media (i.e., τ_⊥ ≪ 1). For a recent review on MC radiative transfer schemes, see Noebauer & Sim (2019).

3.1. Accretion column: Geometry and plasma

We assumed a radiative shock with infinitesimal width at height H, above the NS surface, in a cylindrical accretion column of radius a (see Eq. (16)). The shock was of Rankine-Hugoniot form. Above the shock (pre-shock region; z > H), matter falls freely at a constant speed (see Eq. (17)), while in the post-shock region, below the shock (z < H), the density of the matter is increased by a factor of seven (appropriate for a radiative shock), hot, and has negligible flow velocity. It should be noted that the downstream sinking velocity would be v_d = v_ff/7 ≲ 0.1 c. However, the rms in the velocity due to thermal motion is $v_{\mathrm{rms}}=\sqrt{k{T}_{\mathrm{e}}/m_{\mathrm{e}}{c}^2}c\approx 0.1{(k{T}_{\mathrm{e}}/5\mathrm{keV})}^{1/2}c$. Thus, for temperatures considered in this study, the rms velocity is greater than the sinking flow velocity in the post-shock region and thus electrons motion is mainly thermal. Even if downstream electrons were much colder, the effect of this nonrelativistic flow would be minimal, as it corresponds to a Lorentz factor of γ ≈ 1.005. In essence, the fractional shift in the cyclotron energy due to Doppler would be ∼0.5%, which is definitely negligible for the scope of this study. Therefore, we could safely ignore the sinking motion in the post-shock region.

The temperature of electrons in the post-shock region is not well determined from theory, though it is expected to be proportional to the kinetic accretional energy of electrons (see Kylafis et al. 2014). Thus, we considered various values of kT_e in this study (see Sect. 4.3). The energy (or better the Lorentz factor) distribution of the thermal electrons in the post-shock region was taken to be the relativistic Maxwell-Jüttner one, namely

$$\begin{array}{c}\hfill f_{\mathrm{e}}(\gamma ,T_{\mathrm{e}})\mathrm{d}\gamma =\frac{\gamma \sqrt{{\gamma }^2-1}}{\vartheta K_2(1/\vartheta )}exp(-\gamma /\vartheta )\mathrm{d}\gamma ,\end{array}$$(21)

where $\gamma =1/\sqrt{1-{(v/c)}^2}$ is the Lorentz factor, v is the electron velocity, K₂(x) is the modified Bessel function of the second kind and second order, and ϑ = kT_e/m_ec².

3.2. Seed photons

Seed photon beams from below the shock were injected upward, and each of them was followed by the MC code until its intensity drops by at least a factor of 1000, that is, 99.9% of the injected beam escapes. The modeling of their injection was as follows: photons are emitted isotropically upward from a disk at z = H. That means that we did not perform radiative transfer on the isotropically emitted seed photons in the post-shock region, but instead we assumed an upward flux of such photons.

The input photon energy distribution was considered to be nonmagnetic bremsstrahlung (Becker & Wolff 2007), namely

$$\begin{array}{c}\hfill E\frac{\mathrm{d}N}{\mathrm{d}E}\propto G(E/2k{T}_{\mathrm{e}})exp(-E/k{T}_{\mathrm{e}}),\end{array}$$(22)

where

$$\begin{array}{c}\hfill G(x)=e^x{K}_0(x)\end{array}$$(23)

is the Gaunt factor in the Born approximation for the cross-section (Greene 1959), E is the seed photon energy, and K₀(x) is the modified Bessel function of the second kind, zeroth order.

3.3. Mean free path

A proper treatment of the mean free path is required in order for accurate radiative transfer calculations to be performed. Since we have a two-zone model, the mean free path does not only depend on the direction of propagation and on the incident photon energy, but it varies significantly from the pre-shock region to the post-shock one.

The mean free path for every position in the accretion column and in any direction is given by (Pozdnyakov et al. 1983; Schwarm et al. 2017)

$$\begin{array}{c}\hfill {\langle \lambda (E,\theta )\rangle }_{f_{\mathrm{e}}}=\frac{1}{n_{\mathrm{e}}{\langle \sigma (E,\theta )\rangle }_{f_{\mathrm{e}}}},\end{array}$$(24)

where E,  θ are the incident photon energy and angle with respect to the magnetic field axis, respectively, as measured in the lab (shock) reference frame, n_e is the local number density of electrons, and 〈σ(E,θ)〉_fe is the averaged (over electron speeds) magnetic resonant Compton cross-section, which is given by

$$\begin{array}{c}\hfill {\langle \sigma (E,\theta )\rangle }_{f_{\mathrm{e}}}=\frac{{\int }_{-1}^1{f}_{\mathrm{e}}(\beta )(1-\beta cos\theta )\sigma (E_{\mathrm{rf}},{\theta }_{\mathrm{rf}})\mathrm{d}\beta }{{\int }_{-1}^1{f}_{\mathrm{e}}(\beta )\mathrm{d}\beta },\end{array}$$(25)

where f_e(β) refers to the electron’s velocity distribution function (see Eq. (21)). The subscript rf stands for quantities that are measured in the electron’s rest frame and can be easily obtained through Lorentz transformations from the lab frame to that of the electron, and σ(E, θ) is the scattering cross-section (see Sect. 3.4). Following Harding & Daugherty (1991; see also Schwarm et al. 2017), we averaged the cross-section over only the longitudinal electron temperature; that is, we presumed the motion of electrons to be effectively along the magnetic field axis. Thus, the integration in Eq. (25) was over the electron velocity distribution along the field axis.

In the pre-shock region, the electron velocity distribution function is essentially a Dirac δ-function (i.e., f_e(β) = δ(β − β_ff)),  and therefore the calculation of the mean free path (Eq. (24)) is trivial. In the post-shock region, the electron velocity function is related to Eq. (21), and thus the integral was carried out numerically (see Appendix A).

3.4. Resonant Compton scattering prescription

The radiative transfer calculation requires accurate expressions for the description of the radiation-matter interaction. The complete, quantum electrodynamical Compton cross-sections are thoroughly discussed in the literature (e.g., Harding & Daugherty 1991; Sina 1996; Nobili et al. 2008; Gonthier et al. 2014; Mushtukov et al. 2016), but they are rather cumbersome and their implementation to MC codes turns out to be a computationally expensive task. In a previous paper (Loudas et al. 2021), we offered a simple but accurate prescription for resonant Compton scattering calculations along with approximate expressions for the polarization-dependent differential cross-sections, which are valid for photons with an energy level of E ≪ m_ec² and for cyclotron line energies of E_c ≪ m_ec², that is, for magnetic field strengths of $B\ll B_{\text{cr}}\equiv m_{\text{e}}^2{c}^2/e\hslash =4.41\times {10}^{13}\text{G}$. We employed this prescription in our MC code, as most of the cyclotron lines that have been observed are at cyclotron energies ≪m_ec² (Staubert et al. 2019), and therefore our work is applicable to most of the observed accreting NSs that exhibit a CRSF in their spectra. For CRSFs at energies of ∼m_ec², the implementation of the complete cross-section expressions is inevitable.

The approximate polarization-dependent resonant Compton cross-sections read (Loudas et al. 2021)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\sigma }_{11}}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\approx \frac{3\pi r_{0}c}{8}({cos}^2{\theta }_{\mathrm{rf}}{cos}^2{\theta }_{\mathrm{rf}}^{\prime }{L}_{-}+\frac{1}{2}\frac{B}{B_{\mathrm{cr}}}{L}_{+}),\end{array}$$(26)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\sigma }_{12}}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\approx \frac{3\pi r_{0}c}{8}({cos}^2{\theta }_{\mathrm{rf}}{L}_{-}+\frac{1}{2}\frac{B}{B_{\mathrm{cr}}}{cos}^2{\theta }_{\mathrm{rf}}^{\prime }{L}_{+}),\end{array}$$(27)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\sigma }_{21}}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\approx \frac{3\pi r_{0}c}{8}({cos}^2{\theta }_{\mathrm{rf}}^{\prime }{L}_{-}+\frac{1}{2}\frac{B}{B_{\mathrm{cr}}}{cos}^2{\theta }_{\mathrm{rf}}{L}_{+}),\end{array}$$(28)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\sigma }_{22}}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\approx \frac{3\pi r_{0}c}{8}(L_{-}+\frac{1}{2}\frac{B}{B_{\mathrm{cr}}}{cos}^2{\theta }_{\mathrm{rf}}{cos}^2{\theta }_{\mathrm{rf}}^{\prime }{L}_{+}),\end{array}$$(29)

where index 1 (2) stands for the ordinary (extraordinary) mode; θ_rf and ${\theta }_{\text{rf}}^{\prime }$ are the incident and scattered angles, respectively, with respect to the local magnetic field, B, as measured in the electron’s rest frame; r₀ is the classical electron radius; and c is the speed of light. The quantities L₋ and L₊ are the normalized Lorentz profiles, namely

$$\begin{array}{c}\hfill L_{\pm }=\frac{{\mathrm{\Gamma }}_{\pm }/2\pi }{{({\omega }_{\mathrm{rf}}-{\omega }_{\mathrm{rf}}^r)}^2+{({\mathrm{\Gamma }}_{\pm }/2)}^2},\end{array}$$(30)

where ω_rf and ${\omega }_{\text{rf}}^r$ are the photon frequency and the resonant frequency measured in the electron’s rest frame, respectively. The cyclotron line widths (or transition rates) Γ_± are given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{-}=\frac{4}{3}\frac{m_{\mathrm{e}}{c}^2}{ħ}\alpha \frac{B^2}{B_{\mathrm{cr}}^2},\end{array}$$(31)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{+}=\frac{2}{3}\frac{m_{\mathrm{e}}{c}^2}{ħ}\alpha \frac{B^3}{B_{\mathrm{cr}}^3}\end{array}$$(32)

(Herold et al. 1982), where α is the fine structure constant.

In this paper, we present polarization-independent calculations, and thus we only made use of the polarization-averaged differential cross-section (Loudas et al. 2021), namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\approx \frac{3\pi r_{0}c}{16}(1+{cos}^2{\theta }_{\mathrm{rf}})(1+{cos}^2{\theta }_{\mathrm{rf}}^{\prime })(L_{-}+\frac{1}{2}\frac{B}{B_{\mathrm{cr}}}{L}_{+}).\end{array}$$(33)

The resonant energy $E_{\text{rf}}^r=\hslash {\omega }_{\text{rf}}^r$, expressed in the electron rest frame, is given by (e.g., Nobili et al. 2008)

$$\begin{array}{c}\hfill E_{\mathrm{rf}}^r=\frac{2E_{\mathrm{c}}}{1+\sqrt{1+2(B/B_{\mathrm{cr}}){sin}^2{\theta }_{\mathrm{rf}}}},\end{array}$$(34)

where E_c = ℏω_c is the classical cyclotron energy and ω_c = eB/m_ec. For each scattering event in the code, the new photon energy in the rest frame of the electron after scattering was computed by Eq. (8) of Loudas et al. (2021).

3.5. Formation of the power law and the scattering feature

Some of the photons that are injected isotropically upward escape unscattered, while the rest of them are scattered once or multiple times. It is the latter category of photons, that may criss-cross the shock many times, increasing their energy considerably and eventually producing the high-energy power-law spectrum (see Kylafis et al. 2014). However, a fraction of those photons have energy comparable to the energy difference between the fundamental and the first Landau level (i.e., the cyclotron energy), which undergo resonant scattering with electrons and thus become trapped in the matter until they lose or gain enough energy to escape in the wings. The combination of these processes inevitably leads to the formation of a hard X-ray power-law continuum, including a prominent CRSF in the overall spectrum.

4. Results

In this section, we present the results of our MC polarization-independent calculations in the form of dN/dE, in arbitrary units, as a function of the photon energy E in units of keV. Every spectrum is normalized so that the integral of dN/dE over all energies in the plot is unity. The main parameters in our calculations were the following:

$$\begin{array}{cc}{\tau }_{\perp }\hfill & \text{Transverse Thomson optical depth.}\hfill \\ {\beta }_{\mathrm{ff}}\hfill & \text{Free-fall velocity over the speed of light.}\hfill \\ B\hfill & \text{Magnetic-field strength at the shock"s altitude.}\hfill \\ T_{\mathrm{e}}\hfill & \text{Thermal electrons" post-shock temperature.}\hfill \end{array}$$

In all cases, the seed photon spectrum was bremsstrahlung (see Sect. 3.2) of temperature equal to the post-shock temperature T_e. We studied the emergent spectrum under diverse conditions, that is, under different sets of the aforementioned parameters of our model. To achieve high-resolution spectra, we used a minimum of 10⁸ photons in every MC run.

4.1. Benchmark model

To assess the performance and the validity of our MC code and in order to see the formation of a CRSF in the emergent spectrum clearly, we started with a benchmark model. We injected photon beams in the shock, the energy distribution of which is bremsstrahlung and of a typical post-shock temperature of kT_e = 10 keV. We set the local magnetic-field strength to B = 0.05B_cr, which corresponds to a cyclotron line energy of E_c = 25.55 keV (typical for accreting NSs), but we considered an unrealistic, nonrelativistic bulk velocity of v_ff = 0.05c to suppress the relativistic effects that play a leading role in the spectral formation, as we highlight in Sect. 4.2. We ran the code for three different values of the parameter τ_⊥: 0.01,  0.1,  1. The low optical depths are non-physical, but we used them just for this benchmark calculation. Also, we avoided using optical depths greater than one, as we only aimed to test the effects of resonant scattering at this stage.

In Fig. 1, we show the results. The dashed black line represents the input bremsstrahlung spectrum, while the rest of them correspond to the resulting angle- and polarization-averaged spectra of our MC simulations for different values of the transverse Thomson optical depth: τ_⊥ = 0.01 (blue), τ_⊥ = 0.1 (golden), τ_⊥ = 1 (green). For visualization purposes, we added a vertical, light blue dotted line that indicates the value of the classical cyclotron energy used in these runs. Although the continuum is not altered noticeably, a prominent absorption feature arises around the classical cyclotron energy and becomes broader and deeper as the transverse optical depth increases. Even for τ_⊥ = 0.1 (golden line), which refers to an optically thin medium for the continuum, the absorption feature is clearly formed in the spectrum. In addition, bump features appear in the wings, with the one in the right wing (i.e., for E ≳ E_c) being more prominent.

Fig. 1. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for various transverse Thomson optical depths. The bulk velocity was set to vff = 0.05c. The black dashed line represents the seed Bremsstrahlung photon spectrum of temperature kTe = 10 keV. The solid blue line stands for τ⊥ = 0.01, the solid gold line is for τ⊥ = 0.1, and the solid green line refers to τ⊥ = 1. The vertical blue dotted line indicates the classical cyclotron line energy Ec = 25.55 keV. All spectra are normalized to unity.

The emergent spectra, shown in Fig. 1, exhibit the anticipated characteristics. In particular, a prominent absorption feature appears near the cyclotron energy, which reflects the fact that the mean free path for photons of energy comparable to the resonant energy is very small (see Sect. 3.4), though the medium is optically thin with regard to Thomson scattering. Essentially, a considerable fraction of these photons undergo resonant scattering and may criss-cross the shock (up and down) multiple times until their energies become significantly higher or lower than the resonant one and they acquire a large mean free path and escape the column sideways. These photons appear in the absorption feature’s left or right wing as bump features. The emergence of bumps in the neighborhood of the absorption features is a natural consequence of the shock’s medium being pure scattering, which dictates that the total number of photons must be conserved.

Regarding the shape of the absorption feature, we notice that as the transverse Thomson optical depth increases, more and more photons of energies around the resonant one are resonantly scattered, resulting in a deeper and broader absorption feature. This trend is clearly depicted in Fig. 1. On the other hand, whether most of the resonantly scattered photons escape in the left wing or in the right one is determined mainly by the post-shock temperature. In fact, every time a photon enters the post-shock region, it is more likely to experience resonant scattering, since the cross-section for propagation of photons in a thermal plasma is broadened due to the thermal motion of the electrons (see Appendix A); hence, photons are scattered by the thermal plasma until they gain or lose enough energy to escape. For low-temperature plasma (cold), photons lose energy on average, and they escape primarily in the left wing; while for high-temperature (hot) plasma, photons gain energy on average and escape primarily in the right wing (see also Sect. 4.3).

The results of the benchmark calculations shown in Fig. 1 enable us to understand the main physical processes, but they are not appropriate for describing shocks in accretion columns of NSs, because we have considered a nonrelativistic bulk velocity above the shock. For a typical accreting NS, the free-fall velocity would be something like v_ff = 0.6c.

4.2. Effects of the bulk motion

As we discuss in Sect. 2, the effect of the bulk-motion Comptonization is of unique importance in shaping the resulting spectrum, as the accretion flow is decelerated in the shock supplying its accretional energy to photons, thereby altering their energy distribution. In order to investigate the effect of the relativistic bulk motion of the accreting plasma (pre-shock region), we ran the MC code for τ_⊥ = 1, 3, and 5 with a bulk velocity, which is now chosen to be v_ff = 0.6c. This is appropriate for an NS of mass M = 1.4 M_⊙ and radius R = 10 km.

In Fig. 2, we present the results of our MC simulations for different values of the transverse Thomson optical depth, kT_e = 10 keV, and E_c = 25.55 keV. The input spectrum is displayed with a dashed black line, while the emergent spectra are displayed with solid colored lines: τ_⊥ = 1 (green), τ_⊥ = 3 (orange), and τ_⊥ = 5 (purple). The vertical, light blue dotted line indicates the classical cyclotron energy. As expected, a prominent absorption feature arises in the spectrum. However, in this case, the spectra differ significantly from those presented in Sect. 4.1. First, the centroid of the absorption feature is shifted to lower energies by ∼30% with respect to the classical cyclotron line, and there is a hint that the shift slightly decreases with increasing τ_⊥, though higher levels of precision and fitting the whole spectrum would be required to achieve reliable evidence of this trend. The depth of the absorption feature increases with increasing τ_⊥, in line with Sect. 4.1. However, the absorption feature’s width is significantly broader than the one shown in the results of the benchmark model discussed in Sect. 4.1. Third, a bump appears in the right wing of the line, followed by a hard X-ray tail, similar to the reported spectra of accreting XRPs (e.g., V 0332+53: Doroshenko et al. 2017; 4U 1538-522: Sharma et al. 2023). In addition, a power-law-like spectrum arises at low energies and extends up to the CRSF. The power-law index decreases; that is, the spectrum becomes harder as the transverse optical depth increases, which is in line with Kylafis et al. (2014).

Fig. 2. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for various transverse Thomson optical depths. The bulk velocity was set to vff = 0.6c. The black dashed line represents the seed Bremsstrahlung photon spectrum of temperature kTe = 10 keV. The solid green line refers to τ⊥ = 1, the solid orange line accounts for τ⊥ = 3, and the solid purple line is for τ⊥ = 5. The vertical blue dotted line indicates the classical cyclotron line energy, Ec = 25.55 keV. All spectra are normalized to unity.

The physical explanation of the characteristics of the emergent spectrum follows. As we discuss in Sect. 4.1, the resonant cross-section implies an extremely small mean free path for the propagation of photons of energy comparable to the resonant energy. However, one should take into account that photons emitted just below the shock (shock/lab frame) appear blueshifted in the in-falling plasma (bulk frame) as a consequence of relativity (we recall that v_ff = 0.6c). As a result, photons that have energies lower than the classical cyclotron energy -depending on their direction of propagation- will undergo resonant scattering with the free-falling electrons and subsequently will appear in the blue (right) wing of the CRSF, on average, due to relativistic beaming (i.e., an inverse Lorentz transformation from the bulk-frame back to the shock-frame). Therefore, photons that are resonantly scattered from the in-falling electrons are Doppler boosted twice, thus gaining a considerable amount of energy, on average. This discussion explains both the shifting of the CRSF to lower energies than the classical cyclotron one, due to the Doppler boosting effect, and the creation of a bump in the right wing of the absorption feature with a hard X-ray tail, because scattering preserves the total number of photons. We highlight that the resonantly upscattered photons acquire energies significantly larger than bremsstrahlung allows. Since photons are emitted isotropically upward from below the shock, photons of different directions of propagation (with respect to the magnetic field axis) will experience a different Doppler boosting from the shock frame to the bulk frame; thus, for resonant scattering between photons and free-falling electrons, a range of different photon energies is sampled in view of multiple directions of photon propagation. This leads to the formation of a broad absorption feature, though shifted to lower energies due to the Doppler effect, as can be seen in Fig. 2.

As the transverse optical depth increases, a larger fraction of photons undergo resonant scattering and criss-cross the shock; thus, the CRSF becomes more prominent. However, in the post-shock region, the electrons are thermal with negligible bulk velocity, and, as a result, the absorption feature that is formed below the shock by the photons that enter this region is almost centered at the classical cyclotron energy, though it is thermally broadened due to the thermally averaged cross-section (Appendix A). It turns out that the combination of the Doppler boosting (pre-shock region) and the thermal Doppler broadening (post-shock region) result in a wide, prominent CRSF that is shifted toward lower energies than the classical cyclotron energy. Finally, the continuum on the left of the CRSF becomes harder as the τ_⊥ increases and exhibits a power-law-like shape as a natural consequence of the first-order Fermi energization (for a thorough discussion, see Kylafis et al. 2014).

4.3. Effects of the post-shock temperature

As we argue in Sect. 2, the post-shock temperature is unknown, although it is related to the accretional energy of the electrons. For this reason, we ran our MC code for a wide range of kT_e values to study its effect on the emergent spectrum.

Figure 3 is similar to Fig. 2, but for kT_e = 40 keV. The characteristics of the spectra are the same as the ones shown in Fig. 2. Again, a prominently shifted absorption line is formed, accompanied by a bump in the right wing and a hard X-ray tail. The shift of the line centroid again seems to slightly decrease with increasing τ_⊥. This trend, although not well established, could likely be a consequence of the interplay between the bulk-motion resonant scattering in the pre-shock region and the thermal resonant scattering in the post-shock one. In fact, for larger τ_⊥, a greater fraction of photons criss-cross the shock, which implies that the number of photons that are resonantly scattered both in the pre-shock region and in the post-shock one increases, leading, in this way, to the overlapping of the two different resonant scattering processes. This, in turn, could reduce the effect of the bulk motion resonant scattering, especially when the post-shock temperature is high, and eventually the CRSF may be only partially shifted to the left. However, the width and the depth of the CRSF are somewhat greater than those in Fig. 2 for the same value of τ_⊥. The broadening and the deepening of the absorption feature with the post-shock temperature are consequences of two factors: 1) the seed photon spectrum around the resonant energy is flatter for greater temperatures, thereby allowing for a larger number of photons to be resonantly scattered, and 2) the higher the kT_e, the broader the thermally averaged resonant cross-section (see Appendix A), and thus the wider the absorption feature generated in the post-shock region. The combination of these reasons and the preservation of the total number of photons explains the fact that the peak of the bump (in the right wing) appears at relatively higher energies and is stronger and broader for larger post-shock temperatures with respect to the ones in Fig. 2 for the same τ_⊥. Figures 4 and 5 are similar to Fig. 2, but for kT_e = 5 keV and kT_e = 2 keV, respectively. The conclusions are the same.

Fig. 3. Same as in Fig. 2, but for kTe = 40 keV.

Fig. 4. Same as in Fig. 2, but for kTe = 5 keV.

Fig. 5. Same as in Fig. 2, but for kTe = 2 keV.

For completeness, and for direct comparison, we plot in Fig. 6 the emergent spectrum for B = 0.05 B_cr, v_ff = 0.6c, and τ_⊥ = 1 for four different values of the post-shock temperature: kT_e = 2 keV (blue), kT_e = 5 keV (golden), kT_e = 10 keV (green), kT_e = 40 keV (orange). Figure 7 is the same as Fig. 6, but for τ_⊥ = 5. The two figures indicate that the post-shock temperature and the transverse optical depth affect the shape of the CRSF. The shift of the CRSF is greater for larger temperatures and relatively small optical depths (τ_⊥ ∼ 1). This is because the spectrum is flatter for higher temperatures, and a larger fraction of the seed photons undergo resonant scattering with the in-falling electrons, leading to stronger (shifted-) absorption and bump features (see Fig. 6). On the other hand, the effect of thermal resonant scattering (post-shock region) is cancelled out by the strong bump feature produced via bulk-resonant scattering, and hence it becomes less efficient when it comes to reducing the shift of the line produced by the bulk resonant scattering as long as the transverse optical depth is small (see Fig. 6). However, as τ_⊥ increases, more and more photons are scattered downward by the accretion flow, crossing the shock and entering the post-shock region. Therefore, a greater fraction of photons are thermal-resonantly scattered below the shock, leading to a stronger, non-shifted absorption feature, which impacts the overall shape of the emergent absorption feature and the right wing to a greater extent. This is clearly inferred by comparing Fig. 7 to Fig. 6.

Fig. 6. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for τ⊥ = 1, vff = 0.6c, and Ec = 25.55 keV. Solid lines correspond to different post-shock temperatures kTe: 2 keV (blue); 5 keV (golden); 10 keV (green); and 40 keV (orange). Vertical blue dotted line indicates the classical cyclotron line energy. All spectra are normalized to unity.

Fig. 7. Same as in Fig. 6, but for τ⊥ = 5.

4.4. Effects of the transverse optical depth

As we discuss in the previous sections, the transverse Thomson optical depth has a dominant role in forming the overall spectrum. Its effects are summarized below. First, the depth of the absorption feature increases remarkably, that is, it becomes more prominent for higher values of τ_⊥. Second, the larger the optical depth, the smaller the shift of the line centroid to lower energies seems to be, although resolution limits and noise prevent us from obtaining conclusive evidence of this point. Third, the bump in the right wing becomes stronger with increasing τ_⊥. Fourth, the continuum becomes harder and the power law at energies lower than the cyclotron energy becomes flatter (i.e., the power-law index decreases) for higher values of τ_⊥. Fifth, the hard X-ray tail that emerges on the right of the CRSF extends to higher energies as the optical depth increases and reaches energies significantly higher than the corresponding temperature dictates.

4.5. Effects of the magnetic-field strength

As we discuss in Sect. 3.4, the resonant energy is approximately proportional to the magnetic-field strength. Thus, the emergent CRSF scales almost linearly with B, and it appears near the classical cyclotron energy, though shifted to lower energies, as we demonstrate in Sect. 4.2. Although the magnetic-field strength on the surface of an accreting magnetic NS spans from ∼10¹² G to ∼10¹³ G, we chose to present results for a typical magnetic field strength that corresponds to E_c = 25.55 keVin this work (i.e., B/B_cr = 0.05). For demonstration purposes, we ran our MC code for the set of parameters used in Fig. 4, except for B/B_cr = 0.1, which corresponds to E_c = 51 keV. We present the results in Fig. 8. The solid black line is the input spectrum, while the other lines refer to the emergent spectra for different values of the transverse optical depth: τ_⊥ = 1 (green), τ_⊥ = 3 (orange), and τ_⊥ = 5 (purple). The conclusions are identical to those made in the previous sections. As expected, the CRSF appears near the resonant energy, though redshifted.

Fig. 8. Same as in Fig. 4, but for B = 0.1 Bcr (or Ec = 51 keV).

We note that the magnetic-field strength varies significantly along the column under the dipolar magnetic-field assumption. Therefore, the emergent CRSF’s centroid is expected to vary with the shock’s height. In fact, as the shock’s height H increases, the B-field strength drops as (1 + H/R)⁻³, and hence one might consider that the rate of change of the observed cyclotron line energy with the shock’s height is dominated by this drop of the magnetic-field strength along the column. However, as we demonstrated in the previous sections, the bulk-motion speed also plays a crucial role in the line-forming process. In particular, it shifts the line centroid to lower energies and the amount of shifting depends on the free-fall speed at height H. Given that the free-fall velocity drops as the shock’s height increases, the shift due to the bulk resonant scattering decreases and thus affects the correlation between the line centroid and the shock’s height (i.e., the X-ray luminosity). Therefore, the redshifting effect must be considered when one tries to explain the variation of the cyclotron line energy with the X-ray luminosity. However, a further discussion on this topic is beyond the scope of this work and will be addressed thoroughly in a follow-up study.

5. Discussion and conclusions

In this study, we explored the possibility of cyclotron-line formation in the radiative shock in the accretion column of an accreting magnetic NS of super-critical X-ray luminosity, which was first proposed by Basko & Sunyaev (1976). We developed an MC code, implementing a fully relativistic numerical scheme to carry out radiative-transfer calculations in the shock and compute the emergent spectrum. We used a model in which the shock is a pure mathematical discontinuity in a cylindrical column, assuming that above the shock, matter falls freely, while below it, the plasma is thermal and has negligible flow velocity. We assumed a bremsstrahlung seed photon distribution. The main parameters of our model are the transverse Thomson optical depth, the local magnetic-field strength at the shock’s altitude, the post-shock electron’s temperature, and the free-fall velocity at the shock’s altitude. We properly treated both the bulk-resonant scattering by free-falling electrons in the pre-shock region and the thermal-resonant scattering by thermal electrons in the post-shock one, employing an accurate prescription for resonant scattering in accreting NSs (Loudas et al. 2021).

We provide insights into whether cyclotron line formation in the radiative shock is possible and show that prominent CRSFs are generated in the radiative shock and appear in the overall spectrum. We find that the characteristics of the emergent spectrum are 1) a power law at energies below the energy of the CRSF; 2) a deep, broad absorption feature shifted to lower energies; 3) a bump in the right wing of the CRSF; and 4) a hard X-ray tail that extends above ∼100 keV. We remark that these features are generally observed in accreting NSs (e.g., V0332+53; Doroshenko et al. 2017), though fitting data are still premature. We would like to make it clear that the part of the spectrum next to the funtamental CRSF (blue wing) could be altered by the redistribution of photons undergoing resonant scattering due to harmonics, leading to a sequence of approximately evenly spaced cyclotron lines (corresponding to harmonics) imprinted in the hard X-ray tail of the spectrum. This has already been observed in V 0332+53 (Tsygankov et al. 2006), but such a calculation is beyond the scope of this study.

The power-law spectrum that is formed at low energies is a natural consequence of the combination of bulk and thermal Comptonization, and is due to the first-order Fermi energization. The spectral power-law index depends both on the transverse Thomson optical depth and the post-shock temperature. The larger the optical depth, the flatter the spectrum. Similarly, the higher the post-shock temperature, the lower the power-law index. The power-law-like feature extends up to the CRSF. Thus, its extent depends on the magnetic field strength at the shock.

We demonstrate that the absorption feature’s depth, shift, and width are mainly affected by the free-fall velocity, the transverse optical depth, and the post-shock temperature. The absorption feature becomes deeper and broader as the optical depth or the post-shock temperature increases. The line centroid is shifted to lower energies due to the bulk resonant scattering, in line with the findings of Nishimura (2014). The amount of shift is greater for higher values of the free-fall velocity. For v_ff = 0.6c, typical for accreting NSs, the CRSF is shifted by ∼(20 − 30)%,  and the shift seems to slightly decrease as the optical depth increases (e.g., Fig. 3), although this hint needs to be tested by increasing the resolution in the simulations and fitting accurately the line profile. This result implies that for precise measurements of the magnetic field strength through cyclotron line observations, this effect must be quantified and taken into account.

The bump always appears in the right wing and is followed by a hard X-ray tail. The main parameters affecting these features’ shapes are the free-fall velocity, the optical depth, and the post-shock temperature. The bump is mainly formed by the upscattering of photons that undergo resonant scattering with the in-falling electrons and is reprocessed by the thermal electrons in the post-shock region. Its strength and width increase, and the hard X-ray tail extends to higher energies for greater optical depth or higher post-shock temperature values. Its peak energy depends primarily on the free-fall velocity through the Doppler effect, but it also increases significantly with the optical depth and the post-shock temperature.

It should be noted that the derived spectra, although similar, are generally harder than observed, presumably indicating less bulk Comptonization of the emergent photons, but the overall conclusion also depends on the contribution of the emission of the accretion column below the shock. If one takes into account in the emergent spectrum the radiation emitted in the accretion column far below the shock, provided the emission profile along the column is known, and assumes that kT_e < E_c, then the total emergent spectrum from the accretion column may be softer than the one produced solely at the radiative shock.

One source of uncertainty in our calculations arises from the fact that we considered a rather simplified picture for the accretion column/shock structure. We assumed the shock to be 1) stationary, 2) plane parallel, and 3) of infinitesimal width (i.e., a mathematical discontinuity), satisfying the standard Rankine-Hugoniot jump conditions for a radiation-dominated shock of extremely high Mach number far upstream (appropriate for accreting XRPs in the supercritical luminosity regime). Such assumptions result in a velocity drop of ( ≈ 1/7) across the shock (Basko & Sunyaev 1976). Recent and more advanced studies on the topic have revealed a more complex picture, however (e.g., Postnov et al. 2015). Zhang et al. (2022) show that the shock front exhibits high-frequency oscillations, has a paraboloidal-like shape, and is of finite size (i.e., continuous) rather than a mathematical discontinuity. How these findings would alter our conclusions is difficult to determine without performing more detailed MC radiative transfer simulations. One scenario would be that the Fermi energization process of soft photons at the shock may not be that efficient due to the larger spatial scale of the shock transition (i.e., a smoother density and velocity gradient across the shock), resulting in an overall weaker energization. Apparently, this calculation, although important, far exceeds the scope of this paper.

Another caveat in our consideration is related to the choice of the seed photon energy spectrum. Although the plasma in the post-shock temperature is strongly magnetized, we chose to employ a field-free bremsstrahlung emission spectrum below the shock. A more rigorous choice would be to consider magnetic bremsstrahlung photons that have undergone thermal Comptonization in the post-shock region. However, in this study, we focused on studying the first-order Fermi energization of soft photons in the shock, accounting for the effects of bulk and thermal Comptonization and resonant scattering; yet, we ignored the detailed shape of the seed photon distribution as minor changes in the seed photon spectrum would not alter the results, provided the shock resides in the optically thick regime.

We note that magnetic bremsstrahlung emission (e.g., Riffert et al. 1999) induces anisotropy and yields an angle-dependent seed photon emission spectrum, while cyclotron resonant scattering leads to an angle-dependent emission line in the magnetic bremsstrahlung spectrum. These differences may affect the subtle properties of the line features obtained in this work, but would not alter the conclusions, as the thermally and angle-averaged emission spectrum of magnetic Bremsstrahlung agrees fairly well with the field-free one, so long as the photon energy is lower than the cyclotron energy (see Sect. 7.3 in Becker & Wolff 2007). Given that in most cases in our study the temperature was lower than the cyclotron energy (except when kT_e = 40 keV > E_c = 25.5 keV), the choice between those seed photon distributions would not affect the results significantly.

In summary, we present, for the first time, spectral formation in a radiative shock in the accretion column of an accreting magnetic NS, putting special emphasis on the formation of the CRSF. We offer a qualitative description of the spectral formation, provide a quantitative demonstration of the characteristics of the emergent spectra, and describe how we studied their dependence on the conditions encountered in the radiative shock. Radiation magneto-hydrodynamic (RMHD) simulations (e.g., Zhang et al. 2022) in accretion-column-like structures or hollow funnel geometries, allowing for resonant scattering that enhances matter-radiation interaction, could be extremely beneficial in advancing cyclotron line formation calculations in the future.

Acknowledgments

We would like to thank the anonymous referee for providing a constructive review with insightful comments and important suggestions that helped us improve this manuscript and offered new ideas to explore in the future. NL acknowledges support by the European Research Council (ERC) under the HORIZON ERC Grants 2021 programme under grant agreement No. 101040021. NL acknowledges the “Summer School for Astrostatistics in Crete” for providing training on the statistical methods adopted in this work.

References

Appendix A: Thermally averaged cross-section

This appendix offers multiple plots of the thermally and polarization-averaged scattering cross-section used in our MC code for the scattering of photons by thermal electrons. To compute the thermally averaged cross-section, we used the polarization-averaged cross-section shown in Eq. (33). We substituted it into Eq. (25) and averaged over the electron’s velocity along the magnetic field axis, employing a relativistic Maxwellian probability distribution function (see Eq. 21). In Fig. A.1, we plot the thermally averaged cross-section for a wide range of magnetic field strengths B, electron temperatures kT_e, and angles of incidence θ. As can be inferred, the width of the cross-section increases remarkably for higher temperatures, which is a natural consequence of the thermal electrons’ motion. This implies that magnetic resonant Compton scattering of photons with thermal electrons in the post-shock region can occur for photons that have energy that lies in a larger energy range, thereby leading to a wider absorption feature (thermal broadening). At the same time, though, the peak value decreases since the scattering cross-section is averaged over a normalized probability distribution function. The peak value also decreases for higher magnetic field strengths as the cyclotron line widths Γ_± scale with the magnetic field (see Eqs. 31, 32). For a more detailed analysis of the thermally averaged cross-section, we refer the reader to the recent work of Schwarm et al. (2017).

Fig. A.1. Thermally and polarization-averaged scattering cross-section as a function of the incident photon energy for various magnetic field strengths, electron temperatures, and incident angles. Each subplot refers to a specific (B, kTe) set. In each subplot, different line colors correspond to different incident photon angles. The vertical gray line indicates the classical cyclotron energy.

All Figures

Fig. 1. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for various transverse Thomson optical depths. The bulk velocity was set to vff = 0.05c. The black dashed line represents the seed Bremsstrahlung photon spectrum of temperature kTe = 10 keV. The solid blue line stands for τ⊥ = 0.01, the solid gold line is for τ⊥ = 0.1, and the solid green line refers to τ⊥ = 1. The vertical blue dotted line indicates the classical cyclotron line energy Ec = 25.55 keV. All spectra are normalized to unity.

In the text

Fig. 2. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for various transverse Thomson optical depths. The bulk velocity was set to vff = 0.6c. The black dashed line represents the seed Bremsstrahlung photon spectrum of temperature kTe = 10 keV. The solid green line refers to τ⊥ = 1, the solid orange line accounts for τ⊥ = 3, and the solid purple line is for τ⊥ = 5. The vertical blue dotted line indicates the classical cyclotron line energy, Ec = 25.55 keV. All spectra are normalized to unity.

In the text

	Fig. 3. Same as in Fig. 2, but for kTe = 40 keV.
In the text

	Fig. 4. Same as in Fig. 2, but for kTe = 5 keV.
In the text

	Fig. 5. Same as in Fig. 2, but for kTe = 2 keV.
In the text

	Fig. 6. Angle- and polarization-averaged emergent spectrum dN/dE as a function of energy E for τ⊥ = 1, vff = 0.6c, and Ec = 25.55 keV. Solid lines correspond to different post-shock temperatures kTe: 2 keV (blue); 5 keV (golden); 10 keV (green); and 40 keV (orange). Vertical blue dotted line indicates the classical cyclotron line energy. All spectra are normalized to unity.
In the text

	Fig. 7. Same as in Fig. 6, but for τ⊥ = 5.
In the text

	Fig. 8. Same as in Fig. 4, but for B = 0.1 Bcr (or Ec = 51 keV).
In the text

	Fig. A.1. Thermally and polarization-averaged scattering cross-section as a function of the incident photon energy for various magnetic field strengths, electron temperatures, and incident angles. Each subplot refers to a specific (B, kTe) set. In each subplot, different line colors correspond to different incident photon angles. The vertical gray line indicates the classical cyclotron energy.
In the text