# Spectrum formation in X-ray pulsars at very low mass accretion rate: Monte-
Carlo approach

Alexander A. Mushtukov,1,2 Valery F. Suleimanov,3,4,2 Sergey S. Tsygankov5,2
and Simon Portegies Zwart1
1 Leiden Observatory, Leiden University, NL-2300RA Leiden, The Netherlands
2 Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya
Str. 84/32, Moscow 117997, Russia
3 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1,
D-72076 Tübingen, Germany
4 Kazan Federal University, Kremlevskaya str. 18, Kazan 42008, Russia
5 Department of Physics and Astronomy, FI-20014 University of Turku, Finland
E-mail<EMAIL_ADDRESS>(AAM)

###### Abstract

It has been recently discovered that the transition of X-ray pulsars to the
low luminosity state ($L\lesssim 10^{35}\,{\rm erg\ \rm s^{-1}}$) is
accompanied by a dramatic spectral change. Namely, the typical power-law-like
spectrum with high energy cutoff transforms into a two-component structure
with a possible cyclotron absorption feature on top of it. It was proposed
that these spectral characteristics can be explained qualitatively by the
emission of cyclotron photons in the atmosphere of the neutron star caused by
collisional excitation of electrons to upper Landau levels and further
comptonization of the photons by electron gas. The latter is expected to be
overheated in a thin top layer of the atmosphere. In this paper, we perform
Monte Carlo simulations of the radiative transfer in the atmosphere of an
accreting neutron star while accounting for a resonant scattering of polarized
X-ray photons by thermally distributed electrons. The spectral shape is shown
to be strongly polarization-dependent in soft X-rays ($\lesssim 10\,{\rm
keV}$) and near the cyclotron scattering feature. The results of our numerical
simulations are tested against the observational data of X-ray pulsar A
0535+262 in the low luminosity state. We show that the spectral shape of the
pulsar can be reproduced by the proposed theoretical model. The applications
of the discovery to the observational studies of accreting neutron stars are
discussed.

###### keywords:

X-rays: binaries – stars: neutron – radiative transfer – scattering – stars:
magnetic field – polarization

††pubyear: 2021††pagerange: Spectrum formation in X-ray pulsars at very low
mass accretion rate: Monte-Carlo approach–B

##  1 Introduction

Classical X-ray pulsars (XRPs) are accreting strongly magnetized neutron stars
(NSs) orbiting around optical stellar companions (for a review see e.g.,
Walter et al., 2015). One of the richest families of XRPs are systems with Be
star as a companion (BeXRPs; Reig, 2011). These systems exhibit strong
transient activity that allows us to study radiation’s interaction with a
matter under conditions of extremely strong magnetic fields over a wide range
of mass accretion rates covering up to six orders of magnitude (Doroshenko,
2020).

From an observational point of view, spectra of XRPs at high luminosities
($>10^{36}$ erg s-1) have similar shapes which can be well fitted by a power-
law with an exponential cutoff at high energies (see, e.g., Nagase, 1989;
Filippova et al., 2005). However, it was recently discovered that decrease of
the observed luminosity below this value is accompanied by the dramatic
changes of the energy spectra in several XRPs (Tsygankov et al., 2019a;
Tsygankov et al., 2019b; Doroshenko et al., 2021; Lutovinov et al., 2021),
pointing to varied physical and geometrical properties of the emission region.

In particular, originally it was shown for the transient XRPs, GX 304$-$1 and
A 0535+26, that once their luminosity drops down to $~{}10^{34}-10^{35}$ erg
s-1, the "canonical" spectral shape of their emission undergoes a dramatic
transition into a two-component structure consisting of two humps peaking at
$\sim 5-7$ keV and $\sim 30-50$ keV (Tsygankov et al., 2019a; Tsygankov et
al., 2019b; Lutovinov et al., 2021). The only other source with a similar
spectral structure observed earlier is the persistent low-luminosity XRP X
Persei (Doroshenko et al., 2012), which, however, never showed spectral
transitions as a function of luminosity. In the case of A 0535+26, a cyclotron
absorption feature was also detected on top of the high energy component of
the spectrum. The appearance of the high-energy spectral component was
interpreted by the recombination of electrons collisionally excited to the
upper Landau levels in the heated layers of the NS atmosphere (see Fig. 1,
Tsygankov et al. 2019b). However, no quantitative description had been done.

Dramatic spectral changes in XRPs at very low mass accretion rates can shed
light on different physical processes responsible for the spectra formation at
different luminosities. Depending on the luminosity, one can distinguish
different regimes of accretion, i.e., the interaction of radiation with
accreting plasma:
(i) At high mass accretion rates, the luminosity and local X-ray energy flux
are high enough to stop accretion flow above NS surface in radiation-dominated
shock (Lyubarsky & Sunyaev, 1982). Below the shock region, the matter slowly
settles down to the stellar surface and forms the so-called accretion column
confined by a strong magnetic field (Inoue, 1975; Basko & Sunyaev, 1976; Wang
& Frank, 1981; Mushtukov et al., 2015b; Gornostaev, 2021). Because of magnetic
field confinement and opacity reduced in a strong magnetic field (Herold,
1979), the accretion column can survive under the extreme conditions of high
radiation pressure. The minimal accretion luminosity, which is sufficient to
stop the matter in radiation pressure dominated shock, is called the critical
luminosity $L_{\rm crit}$ (Basko & Sunyaev, 1976). The critical luminosity is
of the order of $10^{37}\,{\rm erg\ \rm s^{-1}}$, and depends on the geometry
of accretion flow and magnetic field strength at the NS surface (Becker et
al., 2012; Mushtukov et al., 2015a).
(ii) At low mass accretion rates (sub-critical regime, when $L<L_{\rm crit}$),
the radiation pressure is insufficient to stop accretion flow above NS
surface, and final braking of the flow happens in the atmosphere of a NS due
to collisions between particles. Sub-critical accretion results in hot spots
at the stellar surface located close to the NS’s magnetic poles. Even in the
sub-critical regime, the accretion flow dynamics can be affected by radiation
pressure, and X-ray spectra can be shaped by the interaction of photons with
the accreting material (see, e.g., Becker & Wolff 2007; Mushtukov et al.
2015c). Only in the case of very low mass accretion rates ($\dot{M}\lesssim
10^{15}\,{\rm g\,s^{-1}}$) the influence of interaction between accreting
matter and radiation above the NS surface is negligible. In this regime,
almost all energy of accretion flow is released in the geometrically thin NS
atmosphere, where the accreted protons transfer their kinetic energy to the
plasma heating. In this case, the geometry of the emitting region is
elementary, and the plane-parallel approximation for the particle heated
atmosphere is a reasonable approximation for a complete description of spectra
formation.

In this work, based on the Monte-Carlo simulations of radiative transfer in a
strong magnetic field, we propose the first simple two-slab model describing
the spectral formation in XRPs at very low mass accretion rates. In Section 2,
we introduce our model set up and the basic assumptions of radiative transfer.
The basic details of our numerical code are given in Section 3. Results of
numerical simulations and their comparison with the observational data are
given in Section 4. Summarizing our results in Section 5, we propose
applications to the observational studies and argue that the proposed concept
of spectra formation at extremely low state of accretion luminosity can be
used in detailed diagnostics of “propeller" state in accretion and be a base
of alternative estimations of NS magnetic field strength.

##  2 Model set up

Figure 1:  Schematic picture of the theoretical model. The X-ray energy
spectrum is originated from the atmosphere of a NS with the upper layer
overheated by low-level accretion. The accretion flow is stopped in the
atmosphere due to collisions. Collisions result in excitation of electrons to
upper Landau levels. The following radiative de-excitation of electrons
produces cyclotron photons. The cyclotron photons are partly reprocessed by
magnetic Compton scattering and partially absorbed in the atmosphere. The
reprocessed photons form a high-energy component of the spectrum, while the
absorbed energy is released in thermal emission and forms a low-energy part of
the spectrum.

###  2.1 Geometry of emitting region

Sub-critical mass accretion rates (see e.g. Basko & Sunyaev 1976; Mushtukov et
al. 2015a) result in the simplest geometry of emitting regions: the X-ray
photons are produced by hot spots at the NS surface. If the mass accretion
rate is well below the critical luminosity value $L\ll L_{\rm crit}$, the
influence of X-ray energy flux on the dynamics of the accretion process is
insignificant (Mushtukov et al., 2015c).

In the case of magnetic field dominated by dipole component, the area of the
spots at the NS surface can be roughly estimated as

$\displaystyle S\approx 3\times
10^{9}\,\Lambda^{-7/8}\,m^{-13/20}\,R_{6}^{19/10}\,B_{12}^{-1/2}\,L_{37}^{2/5}\quad\mbox{cm}^{2},$
(1)

where $\Lambda$ is a constant which depends on accretion flow geometry:
$\Lambda<1$ for the case of accretion through the disc with $\Lambda=0.5$
being a commonly used value (see e.g. Ghosh & Lamb 1978; Lai 2014), $m$ is the
NS mass measured in units of solar masses $M_{\odot}$, $R_{6}$ is the NS
radius measured in $10^{6}\,{\rm cm}$, $B_{12}$ is the magnetic field strength
at the pole of a NS measured in units of $10^{12}\,{\rm G}$, and $L_{37}$ is
the accretion luminosity of XRP measured in units of $10^{37}\,{\rm erg\ \rm
s^{-1}}$. The corresponding effective temperature of a hot spots can be
estimated as

$\displaystyle T_{\rm eff}=\frac{L}{2\sigma_{\rm
SB}S}=6.6\,\Lambda^{7/32}\,m^{13/80}\,R_{6}^{-19/40}\,B_{12}^{1/8}\,L_{37}^{3/20}\quad\mbox{keV},$
(2)

where $\sigma_{\rm SB}$ is the Stefan–Boltzmann constant.

In the case of accretion on weakly magnetized NSs, the protons of accretion
flow brake in a stellar atmosphere due to Coulomb interactions with
atmospheric electrons (see Chapter 3 in Frank et al. 2002). In the case of
highly magnetized atmospheres, the physical picture is similar, but a
significant fraction of proton kinetic energy turns into excitation of
electrons to upper Landau levels (Nelson et al., 1995). Further de-excitation
of electrons results in production of cyclotron photons at energy $E_{\rm
cyc}\approx 11.6\,B_{12}\,{\rm keV}$ (it is a case if $B\lesssim 10^{13}\,{\rm
G}$, see e.g. Harding & Preece 1987). It is expected that the sources of
cyclotron photons decay exponentially with the optical depth $\propto
e^{-\tau/\tau_{\rm br}}$, where $\tau_{\rm br}$ is a typical optical depth due
to Thomson scattering, where the accretion flow is braked due to collisions.
The total number of cyclotron photons emitted in the atmosphere per ${\rm
cm}^{2}$ in a ${\rm sec}$ due to accretion flow braking can be estimated as

$\displaystyle\frac{1}{S}\frac{{\rm d}N_{\rm cyc}}{{\rm d}t}\approx 6.2\times
10^{35}\frac{L_{37}}{E_{\rm cyc,keV}S_{10}}\quad{\rm cm^{-2}s^{-1}},$ (3)

where $S_{10}$ is the area of a hot spot at the NS surface in units of
$10^{10}\,{\rm cm}$.

A temperature structure of the accretion heated weakly magnetized atmosphere
could be divided into two qualitatively different layers. The most of energy
is released in deep optically thick layers (if $\tau_{\rm br}>1$) where the
plasma mass density is high enough and cooling due to free-free emission can
compensate the heating due to braking of the accretion flow. The deep heated
layer is close to being isothermal, with the temperature close to the
effective one (2). However, the upper rarefied layers of the atmosphere have
too low density to be able to cool with free-free emission. As a result, the
upper layers are cooled by Compton scattering and overheated up to tens and
even a hundred keV (see, e.g., Suleimanov et al. 2018). We assume that the
highly magnetized accretion heated atmospheres have a similar qualitatively
structure. Thus, we consider a plane-parallel magnetized semi-infinite
atmosphere with a thin overheated upper layer and an isothermal deeper layer
(see Fig. 1). The temperatures of both parts are free parameters of our model,
as well as the optical thickness of the overheated slab. The magnetic field is
taken to be perpendicular to the NS surface, which is a good approximation for
the regions located close to the magnetic poles of a NS. The resulting
spectrum is computed by solving the radiation transfer equation using the
Monte-Carlo method.

###  2.2 Radiative transfer

Both opacity and refractive index are dependent on photon polarization in a
strongly magnetized plasma and vacuum. As a result, the solution of the
radiative transfer problem has to account for the polarization of photons. In
the general case, the polarization of X-ray photons can be described in terms
of four Stokes parameters. At the same time, the radiative transfer equation
turns into a set of four equations: one for each Stokes parameter. Because the
plasma in a strong magnetic field is anisotropic and birefringent (Ginzburg,
1970) the polarization of X-ray photon can vary along its trajectory. It makes
the description of polarised radiative transfer by the Stokes parameters even
more complicated. However, any X-ray photon can be represented as a linear
composition of two orthogonal normal modes, which conserve their polarisation
state along their trajectories (see, e.g., Zheleznyakov 1977). The linear
composition of the normal modes changes its polarisation due to the difference
in the phase velocity of the modes. If the difference between the phase
velocities of normal modes is large enough, the radiative transfer problem
reduces and can be effectively solved for two normal modes, i.e., the
radiation can be described by specific intensities in two normal polarization
modes $I_{E}^{(s)}$, where $s=1\,\,(s=2)$ corresponds to X-mode (O-mode).

The ellipticity of normal modes is determined by plasma conditions
(temperature, mass density and chemical composition, see e.g., Ginzburg 1970;
Kirk 1980) and vacuum polarization (see Harding & Lai 2006 for review). In
this paper, we neglect the effects of vacuum polarization (see Appendix B),
assume that the dielectric tensor in the atmosphere coincides with the
dielectric tensor of cold plasma, the normal modes are transverse and use the
approximated ellipticity of the modes affected by plasma effects only:

$\displaystyle K_{s}(E,\theta)$ $\displaystyle\equiv$
$\displaystyle-i\left(\frac{E_{x}}{E_{y}}\right)$ $\displaystyle=$
$\displaystyle\frac{2\cos\theta}{\frac{E_{\rm
cyc}}{E}\sin^{2}\theta-(-1)^{s}\sqrt{\frac{E^{2}_{\rm
cyc}}{E^{2}}\sin^{4}\theta+4\cos^{2}\theta}},$

where $E_{x}$ ($E_{y}$) is the photon’s electric field component along
(perpendicular to) the $\mathbf{k}-\mathbf{B}$ plane, where $\mathbf{k}$ is
the vector of photon momentum and $\mathbf{B}$ is the vector of external
magnetic field, and $\theta$ is the angle between the photon momentum and
$B$-field direction (Gnedin & Syunyaev, 1973).

The main processes of interaction between radiation and matter in our model
are

  * •

magnetic Compton scattering,

  * •

cyclotron absorption,

  * •

bremsstrahlung affected by a strong magnetic field.

The radiative transfer equation can be written as

$\displaystyle\cos\theta\frac{{\rm d}I_{E}^{(s)}(\Omega)}{{\rm d}\tau_{\rm
T}}$ $\displaystyle=$ $\displaystyle\frac{\alpha_{\rm
abs}^{(s)}(\Omega)}{\alpha_{\rm T}}I_{E}^{(s)}(\Omega)$ $\displaystyle-$
$\displaystyle\sum\limits_{j=1}^{2}\int\limits_{0}^{\infty}{\rm
d}E^{\prime}\int\limits_{(4\pi)}{\rm d}\Omega^{\prime}$
$\displaystyle\times\left[R_{sj}(E,\Omega|E^{\prime},\Omega^{\prime})I_{E^{\prime}}^{(j)}(\Omega^{\prime})\right.$
$\displaystyle-\left.R_{js}(E^{\prime},\Omega^{\prime}|E,\Omega)I_{E}^{(s)}(\Omega)\right]$
$\displaystyle-$ $\displaystyle\frac{\alpha_{\rm
abs}^{(s)}(\Omega)}{\alpha_{\rm T}}\frac{B_{E}}{2}-S_{\rm
ini}^{(s)}(E,\Omega),$

where $\tau_{\rm T}$ is the optical depth due to non-magnetic Thomson
scattering, $\alpha_{\rm abs}^{(s)}$ is the coefficient of true absorption
(due to free-free and cyclotron mechanisms), $\alpha_{\rm T}$ is the
absorption coefficient due to non-magnetic Thomson scattering, $B_{E}(T)$ is
the Planck function. The redistribution function
$R_{s\,s^{\prime}}(E,\Omega|E^{\prime},\Omega^{\prime})$ describes
redistribution of X-ray photons over energy, momentum and polarization state,
and is related to the probability of photon transition from energy
$E^{\prime}$, direction given by solid angle $\Omega^{\prime}$ and
polarisation state $s^{\prime}$ to energy $E$, direction $\Omega$ and
polarization state $s$ due to magnetic Compton scattering. The radiative
transfer equation (2.2) neglects non-linear effects like induced scattering
because they are expected to be negligible in low-luminosity states.

The first term in the right-hand side of equation (2.2) describes absorption
due to bremsstrahlung and cyclotron mechanism. With the second term, we
account for the photon redistribution due to Compton scattering. The third
term introduces the thermal emission of photons, which is calculated in the
assumption of local thermodynamic equilibrium (LTE). The last term gives the
primary sources of cyclotron emission, which are distributed exponentially in
our particular case $S_{\rm ini}\propto e^{-\tau/\tau_{\rm br}}$.

Figure 2:  The scheme of Monte Carlo simulations performed in the paper.
Figure 3:  $E\,F_{E}$ spectrum (in the arbitrary units) of X-ray radiation
leaving the atmosphere of a NS at very low mass accretion rates. The fiducial
case is given by red solid line and corresponds to the case of atmosphere with
describing by the following set of parameters: $E_{\rm cyc}=50\,{\rm keV}$,
$T_{\rm bot}=3\,{\rm keV}$, $T_{\rm up}=80\,{\rm keV}$, $\tau_{\rm up}=0.1$,
$\tau_{\rm br}=1$, $Z=1$. Different panels show the influence of different
parameters on the final X-ray energy spectrum. (a) The effect of optical
thickness $\tau_{\rm up}$ of the overheated upper layer. (b) The influence of
effective depth $\tau_{\rm br}$, where the accretion flow is stopped by
collisions. (c) The influence of the temperature $T_{\rm up}$ of the
overheated upper layer. (d) The influence of the temperature $T_{\rm bot}$ of
the atmosphere below the overheated upper layer. (d) The influence of the
chemical composition of the atmosphere. (f) The influence of the angle
$\theta_{\rm B}$ between the magnetic field lines and normal to the NS
surface. The run results from $2\times 10^{6}$ photons, leaving the
atmosphere.

####  2.2.1 Compton scattering

The redistribution function due to Compton scattering in (2.2) is given by

$\displaystyle R_{s_{\rm f},s_{\rm i}}(E_{\rm f},\Omega_{\rm f}|E_{\rm
i},\Omega_{\rm i})=\frac{1}{\sigma_{\rm T}}\frac{{\rm d}\sigma_{s_{\rm
f},s_{\rm i}}}{{\rm d}\Omega_{\rm f}{\rm d}E_{\rm f}}([f_{\rm
e}(p_{b})],E_{\rm i},\Omega_{\rm i},\Omega_{\rm f}),$ (6)

where the double differential cross section of scattering by the electron gas
with the electron distribution over the momentum $f_{\rm e}(p_{b})$ is

$\displaystyle\frac{{\rm d}\sigma_{s_{\rm f},s_{\rm i}}}{{\rm d}\Omega_{\rm
f}{\rm d}E_{\rm f}}([f_{\rm e}(p_{b})],E_{\rm i},\Omega_{\rm i},E_{\rm
f},\Omega_{\rm f})$ $\displaystyle=$ $\displaystyle\frac{{\rm d}\sigma_{s_{\rm
f},s_{\rm i}}}{{\rm d}\Omega_{\rm f}}(p_{b}^{*},E_{\rm i},\Omega_{\rm
i},\Omega_{\rm f})$ (7) $\displaystyle\times\left(\frac{{\rm d}p_{b}^{*}}{{\rm
d}E_{\rm f}}\right)f_{\rm e}(p_{b}^{*}),$

the solid angles are related to the spherical coordinates describing the
direction of photon motion $(\theta_{{\rm i},{\rm f}},\varphi_{{\rm i},{\rm
f}})$ as ${\rm d}\Omega_{{\rm i},{\rm f}}=\sin\theta_{{\rm i},{\rm f}}{\rm
d}\theta_{{\rm i},{\rm f}}{\rm d}\varphi_{{\rm i},{\rm f}}$, $p^{*}_{b}$ is
the electron momentum required for the photon transition $\\{E_{\rm
i},\Omega_{\rm i}\\}\longrightarrow\\{E_{\rm f},\Omega_{\rm f}\\}$,

$\displaystyle\frac{{\rm d}\sigma_{s_{\rm f},s_{\rm i}}}{{\rm d}\Omega_{\rm
f}}(p_{b},E_{\rm i},\Omega_{\rm i},\Omega_{\rm f})$ $\displaystyle=$
$\displaystyle\frac{{\rm d}\sigma_{s_{\rm f},s_{\rm i}}}{{\rm d}\Omega_{\rm
f}}(p_{b}=0,E^{*}_{\rm i},\theta^{*}_{\rm i},\varphi_{\rm i},\theta^{*}_{\rm
f},\varphi_{\rm f})$ (8)
$\displaystyle\times\frac{(1-\beta^{2})}{(1-\beta\cos\theta_{\rm f})^{2}},$

is the ordinary scattering cross section by the electron of momentum $p_{b}$
along magnetic field lines, and

$\displaystyle E^{*}_{\rm i}$ $\displaystyle=$ $\displaystyle E_{\rm
i}\gamma(1-\beta\cos\theta_{\rm i}),$ (9) $\displaystyle\cos\theta^{*}_{\rm
i,f}$ $\displaystyle=$ $\displaystyle\frac{\cos\theta_{\rm
i,f}+\beta}{1+\beta\cos\theta_{\rm i,f}}$ (10) $\displaystyle\gamma$
$\displaystyle=$ $\displaystyle[1+p_{b}^{2}/(m_{\rm
e}c)^{2}]^{1/2}=(1-\beta^{2})^{-1/2}$ (11)

account for the Doppler effect and aberration due to transition between
different reference frames.

Compton scattering is affected by a strong magnetic field in the atmosphere of
a NS (Daugherty & Harding, 1986). The scattering cross-section strongly
depends on photon energy $E$, momentum, and polarization state. The scattering
of photons which energy is close to the cyclotron energy $E_{\rm cyc}\approx
11.6\,(B/10^{12}\,{\rm G})\,{\rm keV}$ becomes resonant. In this paper,
Compton scattering is considered in a non-relativistic manner (see e.g.,
Herold 1979). The differential cross sections of Compton scattering by
electron at rest are expressed by the complex amplitudes of scattering
$a^{\rm(p)}_{s_{\rm f}s_{\rm i}}$:

$\displaystyle\frac{{\rm d}\sigma_{s_{\rm f}s_{\rm i}}}{{\rm d}\Omega_{\rm
f}}(p_{b}=0,E_{\rm i},\theta_{\rm i},\varphi_{\rm i},\theta_{\rm
f},\varphi_{\rm f})=\frac{3}{32\pi}\sigma_{\rm T}|a^{\rm(p)}_{s_{\rm f}s_{\rm
i}}|^{2},$ (12)

where $E_{\rm i}$ is the initial photon energy, $s_{\rm i}$ and $s_{\rm f}$
denote the photon polarization state before and after the scattering
respectively. The complex amplitudes $a^{\rm(p)}_{s_{\rm f}s_{\rm i}}$ depend
on exact expression for the polarization modes. The amplitudes can be combined
in a matrix

$\displaystyle\widehat{a}^{\rm(p)}=\left(\begin{array}[]{cc}a^{\rm(p)}_{\rm
11}&a^{\rm(p)}_{\rm 12}\\\ a^{\rm(p)}_{\rm 21}&a^{\rm(p)}_{\rm
22}\end{array}\right),$ (15)

which is given by

$\displaystyle\widehat{a}^{\rm(p)}=\widehat{M}_{\rm pv}(E_{\rm f},\theta_{\rm
f})\,\widehat{a}^{\rm(v)}\,\widehat{M}^{-1}_{\rm pv}(E_{\rm i},\theta_{\rm
i}),$ (16)

where the unitary matrix

$\displaystyle\widehat{M}_{\rm pv}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{1+|K_{2}|^{2}}}\left(\begin{array}[]{cc}-iK_{2}&1\\\
i&K_{2}\end{array}\right)$ (19) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{1+|K_{1}|^{2}}}\left(\begin{array}[]{cc}-i&K_{1}\\\
iK_{1}&1\end{array}\right),$ (22)

and the matrix $\widehat{a}^{\rm(v)}$ is composed of the scattering amplitudes
calculated for the linearly polarised vacuum modes Herold 1979:

$\displaystyle a^{\rm(v)}_{\rm 11}=\frac{E_{\rm i}}{E_{\rm i}+E_{\rm
cyc}}e^{i(\varphi_{\rm i}-\varphi_{\rm f})}+\frac{E_{\rm i}}{E_{\rm i}-E_{\rm
cyc}}e^{-i(\varphi_{\rm i}-\varphi_{\rm f})}$ (23) $\displaystyle
a^{\rm(v)}_{\rm 22}=2\sin\theta_{\rm i}\sin\theta_{\rm f}$
$\displaystyle+\cos\theta_{\rm i}\cos\theta_{\rm f}\left(\frac{E_{\rm
i}}{E_{\rm i}+E_{\rm cyc}}e^{i(\varphi_{\rm i}-\varphi_{\rm f})}+\frac{E_{\rm
i}}{E_{\rm i}-E_{\rm cyc}}e^{-i(\varphi_{\rm i}-\varphi_{\rm f})}\right)$
$\displaystyle\left(\begin{array}[]{cc}a^{\rm(v)}_{\rm 12}\\\ a^{\rm(v)}_{\rm
21}\end{array}\right)=\left(\begin{array}[]{cc}-i\cos\theta_{\rm i}\\\
i\cos\theta_{\rm f}\end{array}\right)$ (28)
$\displaystyle\quad\quad\quad\quad\quad\times\left(\frac{E_{\rm i}}{E_{\rm
i}+E_{\rm cyc}}e^{i(\varphi_{\rm i}-\varphi_{\rm f})}-\frac{E_{\rm i}}{E_{\rm
i}-E_{\rm cyc}}e^{-i(\varphi_{\rm i}-\varphi_{\rm f})}\right).$

The thermal motion of electrons significantly affects both scatterings below
the cyclotron resonance and at the cyclotron resonance. 111The thermal
broadening of the resonance is calculated according to one-dimensional
electron distribution and natural width of Landau levels approximately
calculated according to Pavlov et al. 1991. In our Monte-Carlo code, we use
pre-calculated tables, which describe the redistribution of photons over the
energy, momentum, and polarization states.

####  2.2.2 Cyclotron absorption

Cyclotron absorption at the fundamental depends on the photon polarization
state and is calculated according to Zheleznyakov 1977. In the case of $k_{\rm
B}T_{\rm e}\ll m_{\rm e}c^{2}$, where $k_{\rm B}$ is the Boltzmann constant,
$T_{\rm e}$ is the electron temperature, and $m_{\rm e}$ is the electron mass,
the adopted absorption cross section for extraordinary photons is given by

$\displaystyle\sigma_{\rm s}^{(1)}\simeq 5.67\times 10^{3}\frac{\sigma_{\rm
T}}{B_{12}}\frac{(1+\cos^{2}\theta)}{|\cos\theta|}\left(\frac{m_{\rm
e}c^{2}}{k_{\rm B}T_{\rm e}}\right)^{1/2}$ (29)
$\displaystyle\times\exp\left[-\frac{m_{\rm e}c^{2}}{2k_{\rm B}T_{\rm
e}\cos^{2}\theta}\left(\frac{E-E_{\rm cyc}}{E_{\rm cyc}}\right)^{2}\right].$

The cross section of cyclotron absorption for O-mode photons is smaller:
$\sigma_{\rm s}^{(2)}\approx\sigma_{\rm s}^{(1)}(k_{\rm B}T_{\rm e}/m_{\rm
e}c^{2})$. In the case of $k_{\rm B}T_{\rm e}\sim m_{\rm e}c^{2}$ the
absorption cross section is obtained numerically under the assumption that the
distribution of electrons over the momentum $p_{b}$ on the ground Landau level
is given by

$\displaystyle f_{\rm e}(p_{b})=\frac{1}{2}\frac{e^{-y\gamma}}{K_{1}(y)},$
(30)

where $\gamma=[1+(p_{b}/(m_{\rm e}c))^{2}]^{1/2}$ is the Lorentz factor,
$y=m_{\rm e}c^{2}/(kT_{\rm e})$ is the inverse dimensionless temperature and
$K_{1}$ is the modified Bessel function of the second kind (Mushtukov et al.,
in prep.). Cyclotron absorption at the fundamental results in excitation of an
electron from the ground Landau level to the first excited level. Because the
de-excitation rate of electrons due to emission of cyclotron photons is much
larger than the de-excitation rate due to collisions between particles $r_{\rm
coll}$ (Bonazzola et al., 1979), the majority of cyclotron absorption events
are followed by almost immediate photon emission and can be considered as an
event of Compton scattering (see e.g., Harding & Lai 2006). The probability
$P_{\rm abs,true}$ that the cyclotron absorption event will end up with a true
absorption can be estimated from cyclotron $r_{\rm cyc}$ and collision $r_{\rm
coll}$ de-excitation rates:

$\displaystyle P_{\rm abs,true}\simeq\frac{r_{\rm coll}}{r_{\rm cyc}}\sim
1.7\times 10^{-7}n_{\rm e,21}B_{12}^{-7/2},$ (31)

where $n_{\rm e,21}=n_{\rm e}/10^{21}\,{\rm cm^{-3}}$ is local number density
of electrons, which is taken in our calculations to be dependent on optical
depth in the atmosphere.

####  2.2.3 Free-free absorption

Free-free (bremsstrahlung) opacity is polarization, direction and density
dependent. We calculate the opacity according to Lai & Ho 2003 (see Appendix
A) and under assumption that the ellipticity of normal modes is given by
(2.2). The mass density $\rho$ in the atmosphere is calculated under
assumption of hydrostatic equilibrium: ${\rm d}P/{\rm d}\xi=-\rho g$, where
$P$ is the gas pressure, $\xi$ \- vertical coordinate in the atmosphere, and
$g$ is the gravitational acceleration. In the case of atmosphere composed of a
few layers of fixed temperature, the mass density at optical depth $\tau$ is
given by 222 We use the optical depth due to non-magnetic Thomson scattering,
assuming that the opacity is $\kappa_{\rm T}=0.34\,{\rm cm^{2}\,g^{-1}}$.

$\displaystyle\rho(\tau)=\rho(\tau_{i})+40.6\,\frac{m}{T^{(i)}_{\rm
keV}R_{6}^{2}}[\tau-\tau_{i}]\quad{\rm g\,cm^{-3}},$ (32)

where $\tau_{i}$ are optical depth at the top border of $i$ layer,
$T^{(i)}_{\rm keV}$ is the temperature in the layer.

The initial cyclotron photons are emitted at $\sim E_{\rm cyc}$ with the
thermal broadening, which depends on the direction in respect to magnetic
field lines. The angular distribution of the cyclotron photons is taken to be
isotropic.

##  3 Monte Carlo code

We perform Monte Carlo simulations tracing X-ray photons in plane parallel
atmosphere of magnetized NS. Magnetic field strength and its direction are
fixed and assumed to be constant in the atmosphere, which is reasonable
assumption because of small geometrical thickness of the atmosphere and small
size of hot spots at the NS surface. Non-linear effects of radiative transfer
are not important at low luminosity states and neglected in our simulations.
Tracing X-ray photons in the atmosphere we get angular-dependent energy
spectra of polarized radiation leaving the atmosphere (see Fig. 2).

There are two sources of seed photons in the model: the cyclotron photons due
to radiative de-excitation of electrons excited to upper Landau levels by
collisions of accretion flow with the atmosphere, and thermal photons. The
sources of cyclotron photons are distributed exponentially in the atmosphere.
The cyclotron photons are emitted close to the cyclotron energy within the
thermally broadened line. The distribution of initial thermal photons is
determined by the absorption coefficients in the atmosphere and local
temperature (see the third term in radiative transfer equation 2.2). The
transfer of X-ray photons originated from different initial sources is
considered separately. For each simulation we trace history of $N_{\rm
cyc}^{({\rm i})}$ cyclotron photons and $N_{\rm th}^{({\rm i})}$ thermal
photons. Because we account for free-free and cyclotron absorption in the
atmosphere, the number of cyclotron and thermal photons leaving the atmosphere
- $N_{\rm cyc}^{({\rm f})}$ and $N_{\rm th}^{({\rm f})}$ \- is smaller than
the number of seed photons: $N_{\rm cyc}^{({\rm f})}\leq N_{\rm cyc}^{({\rm
i})}$ and $N_{\rm th}^{({\rm f})}\leq N_{\rm th}^{({\rm i})}$. Simulating
radiative transfer we aim to reach a certain number ($N_{\rm th}^{({\rm
f})},N_{\rm cyc}^{({\rm f})}\sim 2\times 10^{6}$) of photons that leave the
atmosphere. The results of the separately calculated radiative transfer
problems are combined in a final angular and polarization dependent spectra:

$\displaystyle\frac{{\rm d}N_{\rm cyc}^{({\rm f})}(E,\theta,s)}{{\rm
d}E},\quad\quad\frac{{\rm d}N_{\rm th}^{({\rm f})}(E,\theta,s)}{{\rm d}E}.$
(33)

In order to construct the final spectrum accounting for both sources of seed
photons, we normalize the contribution of both sources:

$\displaystyle\frac{{\rm d}N^{({\rm f})}(E,\theta,s)}{{\rm
d}E}=A_{1}\,\frac{{\rm d}N_{\rm cyc}^{({\rm f})}(E,\theta,s)}{{\rm
d}E}+A_{2}\,\frac{{\rm d}N_{\rm th}^{({\rm f})}(E,\theta,s)}{{\rm d}E},$ (34)

where $A_{1}$ and $A_{2}$ are constants. The normalization is performed on the
base of the energy conservation law in the atmosphere. We start with the
simulation of radiative transfer of the cyclotron photons. The accretion
luminosity is proportional to the total energy of seed cyclotron photons in
the simulation:

$\displaystyle L_{\rm acc}\propto\sum_{j=1}^{N_{\rm cyc}^{\rm(i)}}E_{j,{\rm
cyc}}^{\rm(i)},$ (35)

while the part of luminosity due to thermal emission of the atmosphere is
determined by the difference between total energy of seed cyclotron photons
and total energy of reprocessed cyclotron photons leaving the atmosphere:

$\displaystyle\sum_{j=1}^{N_{\rm th}^{\rm(f)}}E_{j,{\rm
th}}^{\rm(f)}=\left[\sum_{j=1}^{N_{\rm cyc}^{\rm(i)}}E_{j,{\rm
cyc}}^{\rm(i)}-\sum_{j=1}^{N_{\rm cyc}^{\rm(f)}}E_{j,{\rm
cyc}}^{\rm(f)}\right].$ (36)

Using conditions (35) and (36) we get constants $A_{1}$ and $A_{2}$ in (34).

In order to perform Monte Carlo simulations and track the photons we use a set
of pre-calculated tables describing photon redistribution due to magnetic
Compton scatterings:

Tables A of total scattering cross-section, where the cross-section is given
as a function of photon energy, polarisation state before and after the
scattering event, the angle between the $B$-field direction and photon
momentum, temperature and bulk velocity of the electron gas. For each
combination of the initial and final polarization state of a photon, the
tables are pre-calculated in a fixed grid in photon energy and angle
$\theta_{\rm i}$, and for a fixed temperature and bulk velocity of the gas. To
get a scattering cross-section for given photon energy and momentum, we use
quadratic interpolation in the photon energy grid and further quadratic
interpolation in the angle grid.

Tables B of photon probabilities to be scattered into a certain segment of the
solid angle $(\theta_{\rm f}+\Delta\theta_{\rm f},\varphi_{\rm
f}+\Delta\varphi_{\rm f})$. The tables are pre-calculated on a grid of photon
initial parameters (energy, polarization state, momentum) and for both
possible final polarization states.

Steps of tracing of photon history are

  1. 1.

We make a choice of seed photon origin: thermal emission of the atmosphere or
cyclotron emission. If the photon is due to the cyclotron emission, we get the
optical depth where the photon is emitted:

$\tau_{0}=-\tau_{\rm br}\ln X_{1},$

where $X_{1}\in(0;1)$ is a random number, and the photon energy assuming that
the photon is emitted within a thermally broadened cyclotron line. If the
photon is due to the thermal emission, we get the optical depth of its
emission in the atmosphere out of pre-calculated cumulative distribution
functions of photon emission accounting for the assumed temperature structure
in the atmosphere and free-free absorption coefficient.

  2. 2.

We get the free path of the photon accounting for scattering and absorption
opacity.

  3. 3.

Using the initial coordinate of the photon and free path length, we get a new
coordinate of a photon, where it is scattered or absorbed. If the new
coordinate is located out of the atmosphere, the photon contributes to the
spectra of X-ray radiation leaving the atmosphere, and we come back to the
step 1. If the photon is still in the atmosphere, we move on to step 4.

  4. 4.

Comparing the cross section of Compton scattering, free-free and cyclotron
absorption, we specify the elementary process at a new photon coordinate. If
the photon is absorbed, we stop trace history of the photon and start tracing
a new photon, i.e. come back to step 1. If the photon is scattered by
electrons, we get its new energy, momentum direction and polarization state on
the base of pre-calculated tables, and come back to step 2.

##  4 Results of numerical simulations

###  4.1 Influence of NS atmosphere conditions on the X-ray spectra

For the case of fixed local magnetic field strength, there are five main
parameters of the performed numerical simulations and resulting spectra of
X-ray photons. These are (i) the optical thickness of the overheated upper
layer due to Thomson scattering $\tau_{\rm up}$, (ii) the typical length of
accretion flow braking in the atmosphere measured in optical depth due to
Thomson scattering $\tau_{\rm br}$, (iii) the temperature of the atmosphere
under the overheated upper layer $T_{\rm bot}$, (iv) the temperature of the
overheated upper layer $T_{\rm up}$, and (v) the chemical composition of the
atmosphere given by the atomic number $Z$. Here we investigate how these
parameters affect the spectrum. To separate effects caused by different
reasons, we compare the results of numerical simulations with the results
based on the following set of fiducial parameters: $T_{\rm bot}=3\,{\rm keV}$,
$T_{\rm up}=80\,{\rm keV}$, $\tau_{\rm up}=0.1$, $\tau_{\rm br}=1$, $Z=1$,
$\theta_{\rm B}=0$ (red line in Fig. 3 a-f).

In most cases, we see the X-ray energy spectrum consisting of two components.
The low energy component is a result of black-body radiation comptonized by
electrons in the atmosphere. The high energy component is a result of the
initial emission of cyclotron photons and their further comptonization by
electrons. The multiple scatterings of cyclotron photons are strongly affected
by the resonance at the cyclotron energy broadened by the thermal motion.
X-ray photons hardly escape the atmosphere at the energies close to the
cyclotron energy, and because of that, the photons tent to escape in red and
blue wings of a cyclotron line. Note that thermal emission also contributes to
the initial photons at cyclotron energy because free-free absorption is
resonant at $E_{\rm cyc}$ in X-mode.

Comparing the results of different numerical simulations, we can make some
conclusions:
(a) The larger optical thickness of the overheated upper layer $\tau_{\rm up}$
does not affect much the low-energy part of X-ray spectra, but influence the
high energy component (see Fig. 3a) affecting photon redistribution around the
cyclotron line.
(b) Smaller optical depth of accretion flow braking $\tau_{\rm br}$ leads to
stronger high-energy component of X-ray spectrum (see Fig. 3b). It is natural
because at smaller $\tau_{\rm br}$ it becomes easier for cyclotron photons to
leave the atmosphere, starting their diffusion from a smaller optical depth.
Additionally, the smaller optical depth of the accretion flow braking results
in a smaller fraction of absorbed cyclotron photons, contributing to the
thermal low-energy part of X-ray spectra.
(c) The overheated upper layer of the atmosphere affects the high-energy end
of X-ray spectra. The lower temperature of the upper layer $T_{\rm up}$ makes
the blue wing of a cyclotron line weaker (see Fig. 3c). If the temperature of
the upper layer is much smaller than the cyclotron energy $E_{\rm cyc}$, the
most of cyclotron photons are scattered into a red wing of the line and leave
the atmosphere at $E\lesssim E_{\rm cyc}$.
(d) The temperature $T_{\rm bot}$ of the atmosphere below the overheated upper
layer shapes the thermal radiation of the atmosphere and affects the low
energy part of X-ray spectra. Increase of the temperature of the bottom
atmosphere results in a corresponding shift of low-energy component (see Fig.
3d).
(e) The chemical composition in the atmosphere affects the cross section of
free-free absorption: the larger the atomic number, the larger the cross
section of free-free absorption. Because the low-energy component is dominated
by the thermal emission of the atmosphere, this component is affected stronger
by the variations in free-free absorption coefficient. Specifically, we see
that the energy spectra tend to be slightly suppressed at low energies for
larger effective atomic numbers $Z$ (see Fig. 3e).
(f) The direction of magnetic field in respect to the NS surface does not
affect much the final spectra integrated over the solid angle (see Fig. 3f).
In our simulations, we see only a slight decrease of the width of the
absorption feature at $E\sim E_{\rm cyc}$. This decrease is likely due to the
fact that the thermal broadening of the cyclotron resonances in Compton
scattering cross section is weaker for the photons propagating across the
field direction.

Figure 4:  Specific intensity $I_{E}$ at the stellar surface at different
directions given by an angle $\theta$ between local $B$-field direction and
photon momentum. Different lines are given for different angles $\theta=0$
(solid), $0.125\pi$ (dashed), $0.25\pi$ (dotted), $0.375\pi$ (dashed-dotted).
Note, that the intensity in red and blue wings of a cyclotron line are
strongly variable. Parameters for simulated spectrum: $E_{\rm cyc}=50\,{\rm
keV}$, $T_{\rm bot}=3\,{\rm keV}$, $T_{\rm up}=80\,{\rm keV}$, $\tau_{\rm
up}=0.1$, $\tau_{\rm br}=1$, $Z=1$.

The specific intensity of X-ray radiation leaving the atmosphere is angular
dependent (see Fig. 4). In particular, a strong angular dependence of the
specific intensity at red and blue wings of a cyclotron line is expected. The
intensity integrated over the energies composes a pencil beamed diagram, which
is slightly suppressed in the direction orthogonal to the stellar surface.
This phenomenon is natural for the case of the atmosphere with the overheated
upper layer: the contribution of the overheated upper layer into the intensity
is smaller in the direction perpendicular to the stellar surface.

The X-ray energy flux leaving the atmosphere is polarization-dependent (see
Fig. 5). The polarization dependence is particularly strong near the cyclotron
energy, which is natural because the strength of cyclotron resonance and even
its existence depends on the polarization state of a photon. At low energies
$E\ll E_{\rm cyc}$ the flux is dominated by X-mode photons because the
scattering cross-section is smaller for this polarization state (Herold, 1979;
Daugherty & Harding, 1986; Mushtukov et al., 2016). However, the difference in
X-ray energy flux at different polarization states is not dramatic because of
the inverse temperature profile in the atmosphere: the upper layer is assumed
to be much hotter than the underling atmosphere. Note, however, that the exact
predictions for polarization require detailed analyses of effects arising from
vacuum polarization and complicated behavior of plasma dielectric tensor under
conditions of high temperatures. The detailed analysis of these effects is
beyond the scope of this paper and will be discussed in a separate
publication.

Figure 5:  X-ray energy spectra at the NS surface in X-mode (blue dashed line)
and O-mode (red dashed-dotted line) and complete spectra (black solid line).
Parameters for simulated spectrum: $E_{\rm cyc}=50\,{\rm keV}$, $T_{\rm
bot}=3\,{\rm keV}$, $T_{\rm up}=80\,{\rm keV}$, $\tau_{\rm up}=0.1$,
$\tau_{\rm br}=1$, $Z=1$.

###  4.2 Comparison to the observational data

In order to verify our model, we compared results of the simulations with the
data obtained during a low-luminosity state of transient XRP A 0535+262 with
$L_{\rm X}=7\times 10^{34}$ erg s-1. The data were adopted from Tsygankov et
al. (2019b) and cover broad-energy band from 0.3 to 79 keV using Swift/XRT and
NuSTAR instruments.

In Fig. 6 we represent the model calculated for the following set of
parameters (see red line):

$\displaystyle E_{\rm cyc}=48\,{\rm keV},\quad\quad\quad
Z=1,\quad\quad\quad\theta_{\rm B}=0$ $\displaystyle T_{\rm bot}=3.5\,{\rm
keV},\quad\quad\quad T_{\rm up}=100\,{\rm keV},$ $\displaystyle\tau_{\rm
up}=0.1,\quad\quad\quad\tau_{\rm br}=0.5.$

The theoretical model represents the spectrum integrated over the solid angle
at the NS surface and is based on Monte Carlo simulation with $2\times 10^{6}$
photons, leaving the atmosphere. In order to compare the theoretical model
with the observed X-ray spectra, we have accounted for spectral changes due to
the gravitational red shift. We assume that the photon energy detected by a
distant observer $E^{\infty}$ is related to the photon energy at the NS
surface $E$ as

$\displaystyle E_{\infty}=E(1-u)^{1/2},$ (37)

where $u=R_{\rm S}/R$, and Schwarzschild radius $R_{\rm S}\approx 3m\,{\rm
km}$. Mass and radius of a NS were taken to be $M=1.4\,M_{\odot}$ and
$R=12\,{\rm km}$. As can be seen, our theoretical predictions are able to
describe complex spectral shape of the source including all observed features.

We note, however, that our theoretical model shows a lack of X-ray photons in
a red wing of a cyclotron line at energy $\sim 15-20\,{\rm keV}$. It is hard
to eliminate this discrepancy in the energy spectrum integrated over the solid
angle. We suppose that this problem can be solved if one accounts for the
exact geometry of NS rotation in the observer’s reference frame and the
precise process of pulse profile formation. We also note that the radiative
transfer at the high energy part of X-ray spectra can be affected by the
effect of vacuum polarisation (see Appendix B), which was not taken into
account in our simulations. This analysis is beyond the scope of the paper and
a matter of further investigations, which will be discussed in a separate
publication.

Figure 6:  Observed $E\,F_{E}$ spectrum of A 0535+262 at accretion luminosity
of $~{}7\times 10^{34}\,{\rm erg\ \rm s^{-1}}$ is given by black circles
(NuSTAR FPMA and FPMB data) and squares (Swift/XRT). Simulated spectrum is
represented by red line. Parameters for simulated spectrum: $E^{\infty}_{\rm
cyc}=39\,{\rm keV}$, $T_{\rm bot}=2.8\,{\rm keV}$, $T_{\rm up}=100\,{\rm
keV}$, $\tau_{\rm up}=0.1$, $\tau_{\rm br}=0.5$, $Z=1$. The run results from
$5\times 10^{5}$ photons leaving the atmosphere.

##  5 Summary and Discussion

###  5.1 Spectra formation at low mass accretion rates

We performed numerical simulations for spectra formation in XRPs at very low-
luminosity state, when the interaction between radiation and accretion flow
above NS surface does not affect X-ray spectra and dynamics of the accretion
flow. Our simulations are coherent with a physical model (see Fig. 1) where
the accretion flow is braked in upper layers of NS atmosphere due to
collisions between particles, and most of kinetic energy is released initially
in the form of cyclotron photons. The spectra leaving the atmosphere of a NS
is a matter of radiative transfer, strongly affected by magnetic Compton
scattering. The essential ingredient of the model is an overheated upper layer
of the NS atmosphere proposed earlier by Suleimanov et al. 2018 for the case
of low-level accretion onto weakly magnetized NSs. Simulated radiative
transfer in the atmosphere was performed under the assumption of LTE and
accounts for Compton scattering of X-ray photons by thermally distributed
electrons, cyclotron, and free-free absorption. Two components in the spectrum
correspond to comptonized thermal radiation (low-energy hump), and comptonized
cyclotron photons (high-energy hump) originated from collisions of accreting
particles with the electrons in the NS atmosphere and further radiative
transition of electrons to the ground Landau level. The absorption feature on
top of the high-energy hump is due to the resonant scattering of X-ray photons
at cyclotron energy, which forces cyclotron photons to leave the atmosphere in
the wings of a cyclotron line.

Using the constructed model, it was possible to reproduce the observed
spectrum of X-ray pulsar A 0535+262 at a very low luminosity state (see Fig.
6, Tsygankov et al. 2019b). Qualitative agreement between simulated and
observed X-ray spectra confirms assumptions of the underlying physical model.
We argue that two-component spectra should be typical for low-level accretion
onto strongly magnetized NSs.

###  5.2 Applications to the observational studies

####  5.2.1 Investigation of the “propeller" state

Decreasing of the mass accretion rates in XRPs down to the very low values
results in the transition of the source either to (i) the “propeller" state,
when the accretion flow cannot penetrate through the centrifugal barrier set
up by the rotating magnetosphere of a NS (see, e.g., Illarionov & Sunyaev
1975; Lipunov 1987; Ustyugova et al. 2006), or (ii) to the regime of stable
accretion from a cold disc (see Tsygankov et al. 2017). “Propeller" effect was
recently detected in a few XRPs (see e.g., Tsygankov et al., 2016; Lutovinov
et al., 2017). Moreover, in some sources, transitions into the “propeller"
state were discovered to be accompanied by dramatic changes of X-ray energy
spectra (Tsygankov et al., 2016). Specifically, in the energy range below
$\sim 10$ keV, the spectra were shown to become significantly softer with
shape changed from the power-law to a black body with a typical temperature
around 0.5 keV. However, it is still unknown if the centrifugal barrier blocks
the accretion process entirely, and the detected soft X-ray spectra are
observational evidence of cooling NS surface, or leakage of matter through the
barrier is still possible and responsible for some fraction of the observed
emission (Wijnands & Degenaar, 2016; Rouco Escorial et al., 2017).

Our theoretical model of spectra formation provides a natural way to
distinguish low-level accretion from the cooling NS surface. The hard
component of X-ray spectra is a direct result and a specific feature of the
accretion process. As a result, low-level accretion in the case of the leakage
of the centrifugal barrier should result in two-component X-ray spectra, while
the “propeller" state without penetration of material through the barrier
should result in single hump soft spectra.

####  5.2.2 Measurements of magnetic field strength

Two-component X-ray energy spectra at low mass accretion rate provide a way to
estimate magnetic field strength at the NS surface. Indeed, the hard component
of X-ray spectra is formed around cyclotron energy, which is directly related
to the field strength: $E_{\rm cyc}\approx 11.6\,B_{12}\,$keV. Thus, the
detection of hard-energy hump provides a way to estimate the cyclotron energy
in the case when cyclotron absorption feature is not seen in the source
spectrum.

For instance, the spectra of X Persei (Di Salvo et al., 1998), according to
our model, implies the cyclotron energy $E_{\rm cyc}\gtrsim 100\,{\rm keV}$
and corresponding magnetic field strength $B\gtrsim 10^{13}\,{\rm G}$. Note
that this estimation is coherent with the results based on timing analyses and
torque models applied to this particular source (Doroshenko et al., 2012).

## Acknowledgements

This work was supported by the Netherlands Organization for Scientific
Research Veni Fellowship (AAM), the Väisälä Foundation (SST) and the Academy
of Finland travel grant 324550. VFS thanks Deutsche Forschungsgemeinschaft
(DFG) for financial support (grant WE 1312/51-1). The authors also thank the
Russian Science Foundation (grant 19-12-00423) for financial support. We are
grateful to an anonymous referee for a number of useful comments and
suggestions which helped us improve the paper.

## Data availability

The calculations presented in this paper were performed using a private code
developed and owned by the corresponding author, please contact him for any
request/question about. Data appearing in the figures are available upon
request. The observational data used in the manuscript are adopted from those
reported in Tsygankov et al. 2019b.

## References

  * Basko & Sunyaev (1976) Basko M. M., Sunyaev R. A., 1976, MNRAS, 175, 395
  * Becker & Wolff (2007) Becker P. A., Wolff M. T., 2007, ApJ, 654, 435
  * Becker et al. (2012) Becker P. A., et al., 2012, A&A, 544, A123
  * Bonazzola et al. (1979) Bonazzola S., Heyvaerts J., Puget J. L., 1979, A&A, 78, 53
  * Daugherty & Harding (1986) Daugherty J. K., Harding A. K., 1986, ApJ, 309, 362
  * Di Salvo et al. (1998) Di Salvo T., Burderi L., Robba N. R., Guainazzi M., 1998, ApJ, 509, 897
  * Doroshenko (2020) Doroshenko V., 2020, MNRAS, 491, 1857
  * Doroshenko et al. (2012) Doroshenko V., Santangelo A., Kreykenbohm I., Doroshenko R., 2012, A&A, 540, L1
  * Doroshenko et al. (2021) Doroshenko V., Santangelo A., Tsygankov S., Ji L., 2021, arXiv e-prints,
  * Filippova et al. (2005) Filippova E. V., Tsygankov S. S., Lutovinov A. A., Sunyaev R. A., 2005, Astronomy Letters, 31, 729
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: Third Edition
  * Ghosh & Lamb (1978) Ghosh P., Lamb F. K., 1978, ApJ, 223, L83
  * Ginzburg (1970) Ginzburg V. L., 1970, The propagation of electromagnetic waves in plasmas
  * Gnedin & Syunyaev (1973) Gnedin Y. N., Syunyaev R. A., 1973, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 65, 102
  * Gornostaev (2021) Gornostaev M. I., 2021, MNRAS, 501, 564
  * Harding & Lai (2006) Harding A. K., Lai D., 2006, Reports on Progress in Physics, 69, 2631
  * Harding & Preece (1987) Harding A. K., Preece R., 1987, ApJ, 319, 939
  * Herold (1979) Herold H., 1979, Phys. Rev. D, 19, 2868
  * Illarionov & Sunyaev (1975) Illarionov A. F., Sunyaev R. A., 1975, A&A, 39, 185
  * Inoue (1975) Inoue H., 1975, PASJ, 27, 311
  * Kirk (1980) Kirk J. G., 1980, Plasma Physics, 22, 639
  * Lai (2014) Lai D., 2014, in European Physical Journal Web of Conferences. p. 01001 (arXiv:1402.1903), doi:10.1051/epjconf/20136401001
  * Lai & Ho (2003) Lai D., Ho W. C. G., 2003, ApJ, 588, 962
  * Lipunov (1987) Lipunov V. M., 1987, Ap&SS, 132, 1
  * Lutovinov et al. (2017) Lutovinov A. A., Tsygankov S. S., Krivonos R. A., Molkov S. V., Poutanen J., 2017, ApJ, 834, 209
  * Lutovinov et al. (2021) Lutovinov A., et al., 2021, arXiv e-prints, p. arXiv:2103.05728
  * Lyubarsky & Sunyaev (1982) Lyubarsky Y. E., Sunyaev R. A., 1982, Technical report, Comptonization in radiation-dominated shocks and spectra of X-ray pulsars
  * Mushtukov et al. (2015a) Mushtukov A. A., Suleimanov V. F., Tsygankov S. S., Poutanen J., 2015a, MNRAS, 447, 1847
  * Mushtukov et al. (2015b) Mushtukov A. A., Suleimanov V. F., Tsygankov S. S., Poutanen J., 2015b, MNRAS, 454, 2539
  * Mushtukov et al. (2015c) Mushtukov A. A., Tsygankov S. S., Serber A. e. V., Suleimanov V. F., Poutanen J., 2015c, MNRAS, 454, 2714
  * Mushtukov et al. (2016) Mushtukov A. A., Nagirner D. I., Poutanen J., 2016, Phys. Rev. D, 93, 105003
  * Nagase (1989) Nagase F., 1989, PASJ, 41, 1
  * Nelson et al. (1995) Nelson R. W., Wang J. C. L., Salpeter E. E., Wasserman I., 1995, ApJ, 438, L99
  * Pavlov et al. (1991) Pavlov G. G., Bezchastnov V. G., Meszaros P., Alexand er S. G., 1991, ApJ, 380, 541
  * Reig (2011) Reig P., 2011, Ap&SS, 332, 1
  * Rouco Escorial et al. (2017) Rouco Escorial A., Bak Nielsen A. S., Wijnands R., Cavecchi Y., Degenaar N., Patruno A., 2017, MNRAS, 472, 1802
  * Suleimanov et al. (2010) Suleimanov V. F., Pavlov G. G., Werner K., 2010, ApJ, 714, 630
  * Suleimanov et al. (2018) Suleimanov V. F., Poutanen J., Werner K., 2018, A&A, 619, A114
  * Tsygankov et al. (2016) Tsygankov S. S., Lutovinov A. A., Doroshenko V., Mushtukov A. A., Suleimanov V., Poutanen J., 2016, A&A, 593, A16
  * Tsygankov et al. (2017) Tsygankov S. S., Mushtukov A. A., Suleimanov V. F., Doroshenko V., Abolmasov P. K., Lutovinov A. A., Poutanen J., 2017, A&A, 608, A17
  * Tsygankov et al. (2019a) Tsygankov S. S., Rouco Escorial A., Suleimanov V. F., Mushtukov A. A., Doroshenko V., Lutovinov A. A., Wijnand s R., Poutanen J., 2019a, MNRAS, 483, L144
  * Tsygankov et al. (2019b) Tsygankov S. S., Doroshenko V., Mushtukov A. e. A., Suleimanov V. F., Lutovinov A. A., Poutanen J., 2019b, MNRAS, 487, L30
  * Ustyugova et al. (2006) Ustyugova G. V., Koldoba A. V., Romanova M. M., Lovelace R. V. E., 2006, ApJ, 646, 304
  * Walter et al. (2015) Walter R., Lutovinov A. A., Bozzo E., Tsygankov S. S., 2015, A&ARv, 23, 2
  * Wang & Frank (1981) Wang Y. M., Frank J., 1981, A&A, 93, 255
  * Wijnands & Degenaar (2016) Wijnands R., Degenaar N., 2016, MNRAS, 463, L46
  * Zheleznyakov (1977) Zheleznyakov V. V., 1977, Electromagnetic waves in cosmic plasma. Generation and propagation.

##  Appendix A Free-free absorption in a strong magnetic field

The opacity due to free-free absorption is calculated in our simulations
according to Lai & Ho 2003. The opacity is dependent on magnetic field
strength, polarization of X-ray photons, their energy and direction of the
momentum in respect to local direction of $B$-field. The opacity for mode $j$
of ellipticity $K_{j}$ can be written as

$\displaystyle\kappa_{j}^{\rm
ff}=\kappa_{+}|e_{+}^{j}|^{2}+\kappa_{-}|e_{-}^{j}|^{2}+\kappa_{0}|e_{0}^{j}|^{2},$
(38)

where $e_{0}$, $e_{\pm}$ are the spherical components of the photon’s unit
polarization vector with

$\displaystyle|e_{\pm}^{j}|^{2}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{\sqrt{2}}(e_{x}^{j}\pm ie_{y}^{j})\right|^{2}$
$\displaystyle=$
$\displaystyle\frac{[1\pm(K_{j}\cos\theta+K_{z,j}\sin\theta)]^{2}}{2(1+K_{j}^{2}+K_{z,j}^{2})},$
$\displaystyle|e_{z}^{j}|^{2}$ $\displaystyle=$
$\displaystyle\frac{(K_{j}\sin\theta-
K_{z,j}\cos\theta)^{2}}{1+K_{j}^{2}+K_{z,j}^{2}},$ (40)

where the ellipticities $K_{j}$ are given by (2.2) and $K_{z,j}$ is taken to
be $0$. The coefficients in (38) are given by

$\displaystyle\kappa_{\pm}=\frac{\omega}{c\rho}v_{e}\Lambda_{\pm}\gamma_{ei}^{\perp},\quad\quad\kappa_{0}=\frac{\omega}{c\rho}v_{e}\gamma_{ie}^{\parallel},$
(41)

where

$\displaystyle\Lambda_{\pm}=[(1\pm u_{e}^{1/2})^{2}(1\mp
u_{i}^{1/2})^{2}+\gamma_{\pm}^{2}]^{-1},$ (42)
$\displaystyle\gamma_{\pm}=\gamma_{ei}^{\perp}+(1\pm
u_{e}^{1/2})\gamma_{ri}+(1\mp u_{i}^{1/2})\gamma_{re}.$ (43)

Dimensionless quantities:

$\displaystyle u_{e}=\left(\frac{E_{\rm cyc}}{E}\right)^{2},\quad
u_{i}=\left(\frac{E_{Bi}}{E}\right)^{2},\quad
v_{e}=\left(\frac{E_{pe}}{E}\right)^{2},$ (44)

where the energies corresponding to electron cyclotron frequency, ion
cyclotron frequency, and the electron plasma frequency are given by

$\displaystyle E_{\rm cyc}$ $\displaystyle=$ $\displaystyle
11.58\,B_{12}\,{\rm keV},$ (45) $\displaystyle E_{Bi}$ $\displaystyle=$
$\displaystyle 6.305\times 10^{-3}\,B_{12}\left(\frac{Z}{A}\right)\,{\rm
keV},$ (46) $\displaystyle E_{pe}$ $\displaystyle=$ $\displaystyle 2.871\times
10^{-2}\left(\frac{Z}{A}\right)^{1/2}\rho^{1/2}\,{\rm keV}.$ (47)

The dimensionless dumping rates due to electron-ion collisions and photon
emission by electrons and ions are given by

$\displaystyle\gamma_{ei}^{\perp,\parallel}$ $\displaystyle=$ $\displaystyle
4.55\times 10^{-8}\,\frac{Z^{2}n_{i,21}}{T_{\rm keV}^{1/2}E_{\rm
keV}^{2}}\left(1-e^{-E/k_{\rm B}T}\right)g^{\rm ff}_{\perp,\parallel},$
$\displaystyle\gamma_{re}$ $\displaystyle=$ $\displaystyle 9.52\times
10^{-6}\,E_{\rm keV},$ (48) $\displaystyle\gamma_{ri}$ $\displaystyle=$
$\displaystyle 7.1\times 10^{-7}\,\frac{Z^{2}}{A}\,E_{\rm keV},$

where $g^{\rm ff}_{\perp,\parallel}$ are magnetic Gaunt factors, which are
calculated according to Suleimanov et al. 2010 in our simulations.

##  Appendix B On the influence of vacuum polarization

Assuming constant temperature $T$ in the atmosphere one can get the dependence
of the mass density on the vertical coordinate $\xi$:

$\displaystyle\rho(\xi)=\rho(\xi_{0})\exp\left[\frac{gm_{\rm p}}{k_{\rm
B}T}(\xi-\xi_{0})\right],$ (49)

where $g\simeq 1.3\times 10^{14}\,mR_{6}^{-2}\,{\rm cm\,s^{-2}}$ is the
acceleration due to gravity at the NS surface. If the opacity is given by
Thomson opacity $\kappa_{\rm e}=0.34\,{\rm cm^{2}\,g^{-1}}$, the local mass
density is related to the optical depth in the atmosphere as

$\displaystyle\rho=40.6\,\tau_{\rm T}\frac{m}{T_{\rm keV}R_{6}^{2}}\quad{\rm
g\,cm^{-3}}.$ (50)

The critical mass density, when the contribution of vacuum polarisation
becomes comparable to the contribution of plasma to the dielectric tensor, can
be estimated as (Lai & Ho, 2003)

$\displaystyle\rho_{\rm V}=9.64\times 10^{-5}\,Y_{\rm e}^{-1}B_{12}^{2}E_{\rm
keV}^{2}f^{-2}\quad{\rm g\,cm^{-3}},$ (51)

where $Y_{\rm e}=Z/A$ and $f\sim 1$. At $\rho\ll\rho_{\rm V}$ the dielectric
tensor is dominated by vacuum effects, while at $\rho\gg\rho_{\rm V}$, the
dielectric tensor is dominated by plasma effects. Using (50) and (51) we get
the optical depth due to Thomson scattering corresponding to $\rho_{\rm V}$:

$\displaystyle\tau_{\rm T,V}=2.4\times 10^{-6}\,Y_{\rm e}^{-1}B_{12}^{2}T_{\rm
keV}E_{\rm keV}^{2}f^{-2}m^{-1}R_{6}^{2}.$ (52)

Therefore, for the case of physical conditions expected in XRP A 0535+262 we
get

$\tau_{\rm T,V}\sim 3\times 10^{-5}T_{\rm keV}E_{\rm keV}^{2},$

which means that the polarization of X-ray photons in low-energy hump is well
described by (2.2), while polarization of X-ray photons in high-energy hump
might be affected by the effects of vacuum polarization, which will be
investigated in a separate publication.