Issue 
A&A
Volume 619, November 2018



Article Number  A114  
Number of page(s)  14  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201833581  
Published online  15 November 2018 
Accretion heated atmospheres of Xray bursting neutron stars
^{1}
Institut für Astronomie und Astrophysik, Kepler Center for Astro and Particle Physics, Universität Tübingen, Sand 1, 72076
Tübingen, Germany
email: suleimanov@astro.unituebingen.de
^{2}
Astronomy Department, Kazan (Volga region) Federal University, Kremlyovskaya str. 18, 420008 Kazan, Russia
^{3}
Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya str. 84/32, 117997 Moscow, Russia
^{4}
Tuorla Observatory, Department of Physics and Astronomy, 20014 University of Turku, Finland
^{5}
Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
Received:
6
June
2018
Accepted:
29
August
2018
Some thermonuclear (type I) Xray bursts at the neutron star surfaces in lowmass Xray binaries take place during hard persistent states of the systems. Spectral evolution of these bursts is well described by the atmosphere model of a passively cooling neutron star when the burst luminosity is high enough. The observed spectral evolution deviates from the model predictions when the burst luminosity drops below a critical value of 20–70% of the maximum luminosity. The amplitude of the deviations and the critical luminosity correlate with the persistent luminosity, which leads us to suggest that these deviations are induced by the additional heating of the accreted particles. We present a method for computation of the neutron star atmosphere models heated by accreted particles assuming that their energy is released via Coulomb interactions with electrons. We computed the temperature structures and the emergent spectra of the atmospheres of various chemical compositions and investigate the dependence of the results on the velocity of accreted particles, their temperature and the penetration angle. We show that the heated atmosphere develops two different regions. The upper one is the hot (20–100 keV) coronalike surface layer cooled by Compton scattering, and the deeper, almost isothermal optically thick region with a temperature of a few keV. The emergent spectra correspondingly have two components: a blackbody with the temperature close to that of the isothermal region and a hard Comptonized component (a power law with an exponential decay). Their relative contribution depends on the ratio of the energy dissipation rate of the accreted particles to the intrinsic flux from the neutron star surface. These spectra deviate strongly from those of undisturbed, passively cooling neutron star atmospheres, with the main differences being the presence of a highenergy tail and a strong excess in the lowenergy part of the spectrum. They also lack the iron absorption edge, which is visible in the spectra of undisturbed lowluminosity atmospheres with solar chemical composition. Using the computed spectra, we obtained the dependences of the dilution and colorcorrection factors as functions of relative luminosities for pure helium and solar abundance atmospheres. We show that the helium model atmosphere heated by accretion corresponding to 5% of the Eddington luminosity describes well the late stages of the Xray bursts in 4U 1820−30.
Key words: accretion / accretion disks / stars: neutron / stars: atmospheres / methods: numerical / Xrays: binaries / Xrays: bursts
© ESO 2018
1. Introduction
Type I thermonuclear Xray bursts have been used to obtain neutron star (NS) parameters. In particular, bursts happening during the hard persistent states of the systems can be used for such analysis (see discussions in Suleimanov et al. 2011a, 2016; Kajava et al. 2014; Poutanen et al. 2014)^{1}. In order to use full information on the variations of the burst spectrum temperature and normalization, the so called cooling tail method or its modification (Suleimanov et al. 2011a, 2017b; Poutanen et al. 2014) have been used. More accurate results can be obtained by directly fitting the burst spectra at different flux levels with NS atmosphere models (Nättilä et al. 2017). However, these hardstate bursts show a spectral evolution that deviates from the theoretically predicted behavior for passively cooling NSs when the burst luminosity drops below a certain level. Furthermore, the higher the persistent flux, the more significant are the deviations. For instance, deviations start to be visible at burst luminosities below 20% of the Eddington value L_{Edd} for the bursts taking place at persistent luminosities of about 0.01 L_{Edd} (Nättilä et al. 2016). On the other hand, the Xray bursts of the helium accreting NS in 4U 1820−30 happen at the relatively high persistent luminosity of about 0.07 L_{Edd}, and such deviations begin at significantly higher burst luminosity of ≈0.7 L_{Edd} (Suleimanov et al. 2017a). This leads us to the conclusion that additional heating by the accreted gas is the cause of the deviations.
Therefore we need to develop NS atmosphere models with an additional energy dissipation in the surface layers. Such models would be applicable to the hardstate bursts. On the other hand, in the soft state the accreted matter likely spreads over the NS surface forming a spreading/boundary layer (Inogamov & Sunyaev 1999; Suleimanov & Poutanen 2006). Its interaction with the NS atmosphere cannot be described by a formalism involving interaction of single particles, but a hydrodynamical treatment would be needed. Thus the modeling of bursts happening in the soft state will not be considered here.
Xray spectra of accreting compact objects such as NSs and black holes in their hard persistent states form in a hot rarefied medium. In black holes, the most likely source of this radiation is the inner geometrically thick and optically thin hot accretion flow (see e.g., reviews by Done et al. 2007; Poutanen & Veledina 2014; Yuan & Narayan 2014). In accreting NSs this radiation maybe produced also at the NS surface (Deufel et al. 2001), which even can dominate the total power. In this paper we consider old weakly magnetized NSs in lowmass Xray binaries (LMXBs).
Zel’dovich & Shakura (1969) were among the first to discuss the emission from a NS surface heated by accreted matter. They considered two physical processes that lead to deceleration of protons in NS atmospheres. The simplest one is Coulomb interaction with electrons. In this case, the main part of the proton kinetic energy is released deep in the optically thick atmospheric layers resulting in the rather thermal emergent spectra. Another possibility is excitation of collective plasma processes by the protons. Here the proton energy is dissipated within an optically thin atmospheric layer and the emergent spectrum may have a shape far different from the blackbody. Alme & Wilson (1973) confirmed these conclusions using first numerical models of accretion heated NS atmospheres. They considered hydrogen NS atmospheres heated by protons falling radially with freefall velocity. The proton deceleration was considered in a selfconsistent way. Two input parameters of the model were considered, namely the NS mass and the accretion luminosity. The most important properties of the accretion heated atmospheres, such as an overheated Comptoncooled outer layers and spectral hardening with the increase of the accretion luminosity were found in this work.
Later the problem was considered in detail by Bildsten et al. (1992). They provided useful relations for describing the proton velocity profiles and the vertical distribution of the accretion heating, which were used by other authors for modeling accretion heated atmospheres (Turolla et al. 1994; Zampieri et al. 1995; Zane et al. 1998). They also considered hydrogen atmospheres heated by accreting protons. They confirmed the basic properties of the heated model atmospheres found by Alme & Wilson (1973), and additionally declared existing of the so called hot solutions, where the temperature can reach almost 10^{12} K (Turolla et al. 1994). The properties of such “hot” solutions were investigated by Zane et al. (1998), who, in particular, demonstrated the importance of electronpositron pair creation. Zampieri et al. (1995) concentrated on the lowluminosity accretion heated atmospheres. They assumed the energy release mainly in the optically thick layers of the atmosphere. Therefore, the temperature structures of these models at the spectrum formation depths were close to that of the undisturbed model atmospheres with the same luminosities. As a result, they obtained spectra harder than the blackbody and very similar to those of undisturbed NS atmospheres (see, e.g., Romani 1987; Zavlin et al. 1996). The main reason is the dependence of the freefree opacity on the photon energy k_{ff} ∼ E^{−3}. Photons prefer to escape from an atmosphere in the most transparent energy bands. Zampieri et al. (1995) also found a clear division between the outer heated atmospheric layers cooled by Compton scattering from the almost isothermal inner part cooled by freefree processes.
Deufel et al. (2001) computed the accretion heated NS atmospheres and proton deceleration selfconsistently and used the results for interpretation of the hard persistent spectra of some LMXBs. For the first time they showed the importance of the temperature of accreted protons on the atmosphere properties. They also considered the bulk velocity as a free parameter.
The results obtained by Zampieri et al. (1995) and Deufel et al. (2001) were widely used for interpretation of the NS spectra in LMXBs during low massaccretionrate states (see e.g., Homan et al. 2014; Wijnands et al. 2015). Also deceleration of the protons in magnetized NS atmospheres was considered before (see e.g., Nelson et al. 1993), and the results were used for the modeling of such atmospheres by Zane et al. (2000).
In this paper we develop a method to compute models of NS atmosphere heated by accreted matter. It is based on the approach developed by Deufel et al. (2001), which we slightly modified. In contrast to previous papers devoted to pure hydrogen atmospheres heated by protons, here we consider arbitrary chemical composition of the atmospheres and an arbitrary mix of protons and αparticles for the accreted particles. We present the basic properties of the accretionheated NS atmospheres and their dependence on the input parameters such as the velocity of accreted particles and their direction, the accretion and the intrinsic NS luminosities. We also compare the developed models to the observed spectral evolution of hardstate Xray bursts.
2. The model
2.1. Main equations
We consider steadystate hot NS atmospheres in planeparallel approximation with additional heating in the surface layers caused by suprathermal accreted particles. The basic principles of the modeling of hot, undisturbed NS atmospheres were presented in our earlier works (Suleimanov et al. 2011b, 2012). This work is based on the formalism described in the latter paper.
We consider a uniform nonrotating NS with mass M, radius R and intrinsic luminosity L. The input parameters for the NS atmosphere model are the surface gravity g, the effective temperature T_{eff}, and the chemical composition. The effective temperature here is a parameter describing the intrinsic bolometric surface flux
where σ_{SB} is the Stefan–Boltzmann constant. Here the argument of F_{0, ν}(0) means that the flux is measured at the surface, while the lower index 0 corresponds to the intrinsic flux. Another parameter, the relative intrinsic NS luminosity ℓ = L/L_{Edd}, could be used instead of T_{eff}. Here and below a relative luminosity means some luminosity normalized to the Eddington luminosity. Both the surface gravity g and the Eddington luminosity L_{Edd} are computed for Schwarzschild metric (see details in Suleimanov et al. 2016; Degenaar & Suleimanov 2018). The Eddington luminosity depends also on the hydrogen mass fraction X.
In the previous works (Suleimanov et al. 2011b, 2012; Nättilä et al. 2015), we computed models of atmospheres consisting of pure hydrogen, pure helium, solar hydrogen/helium mixture enriched by various fractions of heavy elements, from the solar abundance down to an onehundredth of that, as well as heavy metals. Here we consider pure helium and solar composition NS atmospheres only. We express the additional luminosity arising due to external heating by fast particles through the relative accretion luminosity ℓ_{a} = L_{a}/L_{Edd}. We assume that the kinetic energy of the fast particles is thermalized inside the atmosphere and is radiated away through the surface. As a result, the emergent bolometric flux is higher by this relative fraction:
Let us consider for simplicity identical particles with mass m_{i} moving along the normal to the NS surface with the same velocity υ_{0} before penetrating into the atmosphere (more general cases will be described later). The total accretion luminosity is
where
is the Lorentz factor. Here ṁ is the local mass accretion rate (per unit area), which can be related to ℓ_{a} and the NS surface gravity g as
The value of ṁ is conserved in the atmosphere, but the number density n of the accreted particles and their velocity v can change according to the mass conservation law
The flux of the accreting particles creates the ram pressure P_{ram} = ṁ γυ, which manifests itself as an additional acceleration pressing the atmosphere down
where ρ is the plasma density of the atmosphere and m is a Lagrangian independent variable–column density, defined as dm = −ρdr. We note that the velocity υ is positive in the direction of increasing m. The ram pressure force has to be included into the hydrostatic equilibrium equation, the first equation determining the structure of the atmosphere
where P is the gas pressure and g_{rad} is the radiative acceleration (as defined in Suleimanov et al. 2012 for nonisotropic Compton scattering).
The accreted particles do not directly affect the equation of radiative transfer for the specific intensity I(x, μ), and it is the same as for the undisturbed atmosphere (see details in Suleimanov et al. 2012). Here μ = cos θ is the cosine of the angle between the surface normal and the direction of radiation propagation, x = hν/m_{e}c^{2} is the photon energy in units of electron rest mass. We note that the fully relativistic, angular dependent redistribution function is used to describe Compton scattering (see Nagirner & Poutanen 1993; Poutanen & Svensson 1996).
The accretion of fast particles heats the NS surface producing additional radiation flux ṁc^{2}(γ−1). The kinetic energy of fast particles is released in the NS atmosphere gradually and the local energy dissipation rate has to be proportional to the deceleration rate of the particles considered in the next section. The general form of the local energy dissipation rate, which is important for modeling of the heated atmospheres, takes the following form:
We assume that all the local dissipated energy is transformed to radiation, therefore, we can present the change in the radiation flux as
The energy balance equation now takes the form
and has to be fulfilled at the every depth in the atmosphere (for the undisturbed atmosphere Q^{+} = 0). Here k(x) is the “true” absorption opacity, and σ(x, μ) is the electron scattering opacity (see definitions in Suleimanov et al. 2012).
Thus the hydrostatic equilibrium and the energy balance equations describing the NS atmosphere, can be modified easily to account for accretion of fast particles with equal mass and velocity along the normal. It is clear that the dependence of the particle deceleration dv/dm on depth fully determines the atmosphere model (in addition to the standard parameters defining the undisturbed atmosphere). Below we derive the deceleration function and generalize the treatment to a distribution of particles obliquely accreting to the NS surface.
2.2. Stopping of fast particles in the atmosphere
Let us consider a fast particle having velocity υ directed into the atmosphere, mass m_{i} and charge Ze moving through the plasma with temperature T and electron number density n_{e}.
Here and further e is the charge of electron and m_{e} is the mass of electron. The particle loses energy via Coulomb interactions with the plasma electrons. The basic physics of deceleration of such particles is well established (see e.g., Alme & Wilson 1973; Deufel et al. 2001, and chapter 3 in Frank et al. 2002). The main parameter describing deceleration rate is the slowingdown timescale
which can be solved for the deceleration
The corresponding energy loss is
Here ln Λ is the Coulomb logarithm:
for the case υ > α c ≈ c/137 (Larkin 1960) (α is the finestructure constant). At low particle velocities (υ ≪ α c), Eq. (15) gives too small values and we use instead the standard Coulomb logarithm for a thermal plasma
where k is Stefan–Boltzmann constant. In computations we choose the larger value of ln Λ and ln Λ_{T}.
The function f(x_{e}) takes into account the reduction of the deceleration force at low particle velocities and is expressed as
where ϕ(x) is the error function, and x_{e} = (m_{e}υ^{2}/kT)^{1/2} is the ratio of the particle velocity to the averaged thermal electron velocity. For the computations we use the approximation
suggested by Alme & Wilson (1973).
2.3. Heating by thermally distributed fast particles
The Xray spectra of LMXBs in their hard persistent states are well described with the spectrum of radiation Comptonized in a hot electron slab with optical depth τ_{e} ∼ 1 and temperature kT_{e} ∼ 15 − 50 keV (Barret et al. 2000; Burke et al. 2017). This radiation may be produced in the optically thin and geometrically thick accretion flow which could be described with the advectiondominated accretion flow (ADAF) model (see e.g., Yuan & Narayan 2014). The ADAF model differs from the standard accretion disk model (Shakura & Sunyaev 1973) at least in two respects. The gravity of the central object is balanced not only by the rotation, but the radial gas pressure gradient as well. It is also important that the ion temperature in ADAF model is close to the virial one at a given radius r
and can reach 10^{12} K in the vicinity of the NS. Such high ion temperature provides a significant accretion flow thickness h, comparable with the radial distance, h/r ∼ 1, and a high radial velocity of the matter, comparable with the freefall velocity . These properties imply a quite large impact angle of the accreting ions to the NS surface normal, and we consider this angle Ψ as a parameter. Another parameter is the ratio of the bulk velocity of the particles to the freefall velocity at the NS surface η = V_{0}/V_{ff}(R) (we follow here Deufel et al. 2001). Because the mean ion velocity at the virial temperature is close to V_{ff}, the distribution of the ions over velocities has to be taken into account (Deufel et al. 2001). The ADAF could be so hot and rarefied that the Coulomb interactions are not efficient enough to establish the Maxwellian distribution (Mahadevan & Quataert 1997). Thus potentially the ion velocity distribution may significantly deviate from the Maxwellian. In this work, however, we assume that ions follow the relativistic Maxwellian distribution and we describe the ion temperature through its ratio to the virial temperature, χ = T_{i}/T_{vir}(R).
It is necessary to note that in addition to the bulk kinetic energy, thermally distributed accreted particles contribute their thermal energy (enthalpy) to the atmosphere. Therefore, the accretion luminosity is larger than that given by Eq. (3):
where
is the bulk Lorentz factor. The contribution of thermal energy is important at χ > 0.1. We note that the input parameter in our model is a relative accretion luminosity ℓ_{a}, and therefore in our numerical computations we correct the local mass accretion rate obtained using Eq. (5) to keep ℓ_{a} fixed. The second important thing is that the parameters η and χ are not completely free. The total gravitational potential energy of the accreted matter related to the freefall Lorentz factor Γ_{ff} as ṁc^{2}(Γ_{ff}−1) is released by three means. It is partially radiated away by the hot accretion flow, with the corresponding flux F_{flow}. The rest is split up between the enthalpy of the hot matter and the energy of the bulk motion. In nonrelativistic approximation, the latter can be presented by the equation . Thus we have an additional constraint on the parameters η^{2} + χ < 1.
Furthermore, for the isotropic particle distribution in the comoving frame there is a significant number of particles penetrating into the atmosphere and contributing their energy even when the bulk velocity is parallel to the surface, when the contribution of the bulk kinetic energy formally is absent. This effect is especially significant for the ram pressure. We account for this effect in our calculations.
The relativistic generalization of the Maxwellian distribution is the MaxwellJüttner distribution in the comoving frame
where γ′ is the particle Lorentz factor
υ′ is the particle velocity, p′=γ′υ′/c is the dimensionless momentum, p′=p′, Θ is the normalized Gibbs parameter Θ = m_{a}c^{2}/kT_{i}, and K_{2}(Θ) is the Macdonald function (the modified Bessel function of the second kind). The distribution is normalized to unity:
When the gas moves with some bulk velocity defined by the dimensionless average momentum P and corresponding Lorentz factor , the MaxwellJüttner distribution takes the form (Belyaev & Budker 1956; Nagirner & Poutanen 1994)
Let us introduce two Cartesian coordinate systems with z axis directed along the bulk velocity of the fast particles (see Fig. 1). The first one moves with the bulk velocity (the comoving frame), and the second one (laboratory) is connected with the atmosphere. In the comoving frame, the particle distribution follows the MaxwellianJüttner distribution Eq. (22). In the laboratory frame, two components of the momentum, p_{x} = p′_{x} and p_{y} = p′_{y}, retain their values, but the third one is transformed as
Fig. 1. Coordinate systems used in the paper. We also show the momentum of a single particle and its components with the corresponding angles. 
We note that the normalization of the distribution Eq. (25) remains unchanged
We note that the Lorentz factor Γ in the denominator of Eq. (25) appears due to the Lorentz contraction of the elementary volume (dz = dz′/Γ).
To estimate the energy deposition rate, we consider the impact of the particles distributed over the momentum components according to Eq. (25), and assume that they move independently inside the atmosphere. Let us consider a particle with the momentum p = (p_{x}, p_{y}, p_{z}). It has the initial velocity
and moves along the line that makes angle ψ with the surface normal (directed toward the surface):
where p_{r} = −(p_{z} cosΨ + p_{x} sinΨ), see Fig. 1. For the steady state flux of accreting particles we can replace the total time derivative in Eq. (14) with the space derivative along the displacement direction s, and compute the deceleration
Here we used a relativistic representation for the kinetic energy of the particle γm_{i}c^{2} and replaced the derivative from Euler (r) to Lagrangian (m) coordinates. Finally, using the energy dissipation rate due to Coulomb interaction given by Eq. (14) we get
For a given atmosphere model (with given T(m), ρ(m) and n_{e}(m)), we solve this equation to find υ(m). The solution is then used to determine the specific energy deposition rate per unit mass (see also Eq. (9)):
Integrating this equation over the possible initial momenta (with p_{r} < 0), we get the total energy deposition rate by accretion:
The total energy released in the atmosphere is found by the integration over the depth
Similarly, we compute the contribution per unit mass to the ram pressure force along r (see Eq. (7)):
We note that only the normal component of the force is important. Thus the total ram pressure acceleration is obtained by integrating over all momenta (with p_{r} < 0)
Here and before we determine the contributions to the energy dissipation rate and ram pressure of fast particles per unit mass in order to separate the mass accretion rate ṁ in the final Eqs. (33) and (36).
2.4. Thermal conduction
Previous works (Alme & Wilson 1973; Deufel et al. 2001) demonstrated that the temperature gradient in the accretionheated atmospheres can be significant, and the energy transport via thermal conduction may be not negligible. The thermal conduction flux is
where
is the heat transfer coefficient. The heating/cooling rate of the matter due to thermal conduction is given by
The total heating rate in Eq. (9) then becomes
2.5. Computation of atmosphere structure
The general method of computation of the atmosphere structure is the iterative temperature correction approach. It was described in our previous paper (see Sect. 2 of Suleimanov et al. 2012). Here we report only the differences caused by the inclusion of additional energy dissipation in the atmosphere. As in Suleimanov et al. (2012), we start iterations with the gray atmosphere model defined on a grid of Rosseland optical depths. Other physical quantities such as gas pressure, density, number densities of the ions and atoms included into the model, opacities, and the column density grid m are computed for the given temperature distribution on the Rosseland optical depth grid. Further computations are performed on the column density grid m. Using the current model atmosphere we compute the accretion and the conductivity heating using the method described above. We also derive the vertical radiation flux distribution through the atmosphere. As we have not accounted for the contribution of the thermal energy of the accretion flow (see Eq. (20)) when we evaluated the local mass accretion rate (Eq. (5)), the obtained additional emergent flux F_{a} differs from the value determined by the input parameter ℓ_{a}, namely ℓ_{a}F_{0}/ℓ. Therefore, we correct the mass accretion rate using the correction factor C determined by
The ram pressure acceleration Eq. (36) as well as the energy generation rate Eq. (33) are corrected using ṁ_{c} instead of ṁ.
Then we solve the radiation transfer equation with Compton scattering taken into account (see details in Suleimanov et al. 2012) for the current model atmosphere. The equation is solved on the chosen photon energy grid covering all the range where the atmosphere radiates. In a selfconsistent atmosphere the integral radiation flux at every depth
has to be equal to the values determined before known from the physics of the atmosphere. For instance, the integral flux has to be constant over the depth in an undisturbed atmosphere. In our case the integral total flux (a sum of the radiation flux and the flux due to thermal conductivity (F(m)+F_{c}(m)) at every depth is determined by two terms: the intrinsic radiation flux and the flux generated by accretion at depths deeper than the considered depth. Therefore, the final radiation flux distribution over the depth is determined by the following equation
Here we also take into account the correction to the mass accretion rate. In fact, this flux distribution is the integral form of the energy balance equation. It has to be fulfilled together with the differential form of the energy balance Eq. (11).
Both forms of the energy balance equation are not satisfied after the first and many following iterations. The imbalance values for both forms at every column density are used for computation of the correction to the temperature at every depth (see details in Suleimanov et al. 2012). All the computations are repeated for the new temperature distribution, and such iterations are repeated to satisfy the energy balance at every atmospheric point with the relative flux accuracy better than 1%.
2.6. Comparison with previous works
The main difference between the model presented above and the models used by other authors is the inner boundary condition. We considered atmospheres with a given intrinsic radiative flux F_{0}, while Deufel et al. (2001) and Turolla et al. (1994) assumed that it is absent, F_{0} = 0. Formally, Zampieri et al. (1995) introduced the intrinsic luminosity as a parameter, but in fact, they put the inner boundary slightly above the penetration depth of fast particles and did not give the value for the intrinsic luminosity. We believe it means that the intrinsic luminosity is small in comparison with the accretion luminosity. Therefore, we can compare the previously computed models only with our models at low intrinsic luminosity.
3. Results
The physical approach presented in Sect. 2 was incorporated into the existing stellar atmosphere code described in Suleimanov et al. (2012). Using this code we computed a set of models for nonmagnetized NS atmospheres heated by accretion. The NS mass and radius were fixed at M = 1.663 M_{⊙} and R = 12 km, corresponding to the gravitational redshift 1 + z = 1.27, the surface gravity log g = 14.3, the free fall velocity V_{ff} = (2GM/R)^{1/2} = 1.92 × 10^{10} cm s^{−1} = 0.64 c and the virial temperature T_{vir} ≈ Ā × 10^{11} K ≈ Ā 67 Me V at the NS surface. Here Ā is the average ion mass in units of the proton mass; Ā = 1 for pure hydrogen matter, Ā = 4 for pure helium matter, and Ā ≈ 1.3 for a solar H/He mix. The Eddington luminosity also depends on the chemical composition and can be determined once the hydrogen mass fraction X is fixed, L_{Edd} ≈ 3.8 × 10^{38} (1 + X)^{−1} erg s^{−1}; X = 0.7374 for the solar H/He mixture. The varied parameters are the relative internal luminosity ℓ = L/L_{Edd}, the relative accretion luminosity ℓ_{a} = L_{a}/L_{Edd}, the angle between the bulk velocity vector of the accreted ions and the surface normal Ψ, the bulk velocity of the accreted ions with respect to the free fall velocity η = V_{0}/V_{ff}, and the relative temperature of the accreted ions χ = T_{i}/T_{vir}.
3.1. Atmospheres of low intrinsic luminosity
Let us first consider a model atmosphere with the low intrinsic luminosity, ℓ = 0.001, and with a relatively low accretion luminosity, ℓ_{a} = 0.05, which is, however, fifty times larger than the intrinsic luminosity. For simplicity, we consider a pure helium atmosphere heated by accreted αparticles. As the fiducial set of parameters, we chose the temperature of the inflowing ions, χ = 0.2, and their bulk velocity η = 0.75. As the first example, we consider models with different angles of the incoming flow Ψ = 0°, 30°, 60°, and 85°. The emergent spectra and the temperature structures for these four models are shown in Fig. 2. Almost all the energy is released at the column density of a few g cm^{−2} (Fig. 3, bottom panel) and the total radiative flux grows rapidly at those depths (Fig. 3, top panel). The column density of the maximum energy release depends on the angle Ψ: the smaller the angle the deeper the energy release maximum. The atmosphere heated by the fast particles can be divided into an optically thick inner part (m ≥ 1 − 10 g cm^{2}) and an optically thin outer part. The plasma temperature in the optically thick, almost isothermal part is slightly higher in comparison with that of the undisturbed atmosphere and the radiation flux is controlled by the temperature gradient. We can see that this part of the atmosphere is hotter for smaller angles Ψ. At small column densities, the atmosphere becomes optically thin and the energy generation here cannot be balanced with thermal radiation by plasma with kT ≈ kT_{eff} and these layers are overheated. As a result their energy balance is determined by Compton scattering.
Fig. 2. Emergent energy spectra (top panel) and temperature structures (bottom panel) for four lowluminosity (ℓ = 0.001) helium model atmospheres heated by αparticles with the accretion luminosity ℓ_{a} = 0.05 and parameters η = 0.75 and χ = 0.2. The models correspond to different incoming angles Ψ = 0°, 30°, 60°, and 85°. The emergent spectrum of the most inclined (Ψ = 85°) particle flow is well fitted at photon energies higher than 5 keV with an exponential cutoff powerlaw model with photon index Γ = 2.45 and cutoff energy E_{cut} = 85 keV (dashed curve). The blackbody spectrum of temperature kT = 1.1 keV is also shown by the dotted curve. 
Fig. 3. Bolometric flux (top panel) and heating rate (bottom panel) as functions of the column density for the models presented in Fig. 2. 
The heating rate of the upper atmospheric layers by the accreted particles increases with increasing angle Ψ. Therefore, models with the larger angle Ψ have higher temperatures of the upper layers and they transit to overheated states at larger m. The upper layers of the atmosphere with Ψ = 85° are heated up to 10^{9} K, and the emergent spectrum is close to the Comptonized spectrum of a hot electron slab with Thomson optical depth of about unity. The emergent flux can be fitted by a power law with exponential cutoff, F_{E} ∝ E^{−(Γ − 1)}exp(−E/E_{cut}), at photon energies larger than 5 keV (Fig. 2, top panel). The photon energy cutoff E_{cut} = 85 keV is approximately equal to the surface temperatures kT. On the other hand, the blackbody component with a temperature close to the temperature of the heated optically thick part of the atmosphere dominates the spectra of the models with small Ψ, and the optically thin layers add some high energy tails to the spectra. The spectrum of the model with Ψ = 60° has some intermediate shape with a clear blackbody component and a significant hard tail. Such division of the emergent spectra into two components is obvious in the lowluminosity models computed by Alme & Wilson (1973) and Deufel et al. (2001). The temperature structure of these models as well as of those computed by Zampieri et al. (1995) show an almost isothermal inner heated slab and the overheated Comptoncooled upper layers with surface temperatures between 10^{8} and 10^{9} K.
We note also, that the emergent spectra have significant excess at low photon energies in comparison with the blackbody. The atmospheres are very opaque at low energies due to the freefree opacity and the emergent radiation at these energies therefore forms in the overheated upper layers with the escaping monochromatic fluxes close to the blackbody flux with temperatures of those layer. This effect was mentioned by Deufel et al. (2001) as the “inverse photosphere effect”.
The angular distribution of the emergent spectra is also very unusual for stellar spectra (see Fig. 4). The normal limb darkening at the energies below 10 keV alternates with the limb brightening at higher energies, as it was suggested earlier by Poutanen & Gierliński (2003) for the accreting millisecond pulsars. If the bulk velocity vector of the inflowing particles is highly inclined to the atmosphere normal, their contribution to the ram pressure gradient and to the ram pressure itself is small in comparison with that for the particle inflowing along the normal (see Fig. 5).
Fig. 4. Emergent specific intensity spectra of lowluminosity (ℓ = 0.001) helium atmosphere model heated by αparticles with accretion luminosity ℓ_{a} = 0.05 at three zenith angles of 27°.5 (dotted curve), 60° (solid curve), and 83°.5 (dashed curve). The parameters of the model are Ψ = 85°, η = 0.75 and χ = 0.2. 
Fig. 5. Relative ram pressure acceleration (top panel) and the relative ram pressure (bottom panel) as functions of the column density for the models presented in Fig. 2. 
We took into account the thermal conductivity flux in the models, but it is not significant. The relative heat flux does not exceed a few tenth of percent in comparison with the intrinsic radiative flux even for the atmosphere with the hottest upper layers (see Fig. 6, top panel). The contribution of the heat conduction flux derivative into the differential form of the energy conservation law (see Eqs. (11), (39), (40)) is also insignificant (see bottom panel of Fig. 6). However, the energy generation rate by accretion in the upper layers is small and comparable with the corresponding rate provided by conduction (compare bottom panel of Fig. 6 to that of Fig. 3). In spite of the fact that the atmosphere is nearly isothermal, the heat flux by conduction becomes here significant, because of the strong temperature dependence (see Eqs. (37) and (38)). In this case, small numerical fluctuations in temperature of random signs lead to large derivatives of the conduction flux, comparable to the accretion heating. This leads to numerical and, possibly, physical instability of the upper layers where m < 10^{−4} g cm^{−2}. We suppressed this instability in our computations, but it has to be studied in future investigations.
Fig. 6. Relative heat conduction flux (top panel) and heat conduction flux derivative (bottom panel) as functions of the column density for the models presented in Fig. 2. 
Despite the high temperatures of the upper layers, the computed atmospheres are geometrically thin, as it was mentioned by other authors before (see, e.g., Turolla et al. 1994), with a thickness of a few tens of meters at most (see Fig. 7). It is not surprising because the highest model temperature (≈10^{9} K) is significantly lower than the virial temperature for the protons. However, the physical picture can potentially change if the electronpositron pair production is taken into account (Zane et al. 1998). The virial temperature for the electronpositron pairs is close to 10^{9} K, therefore, we can expect an outflow of pairs at some model parameters. The effects of pairs on the atmosphere models will be considered in a separate paper.
3.2. Accretion heated atmospheres with various input parameters
We have considered above how the physical properties of the heated atmospheres depend on the incoming angle of the accretion flow. We now concentrate on the influence of other parameters on the emergent spectrum and the temperature structure. We consider pure helium models with two intrinsic luminosities ℓ = 0.001 and 0.5, having accretion luminosity ℓ_{a} = 0.05, and other parameters Ψ = 60°, η = 0.75 and χ = 0.2 as in the fiducial model. We study then how variations in every input parameter affect the emergent spectrum and the temperature structure.
The obtained results for the intrinsic luminosity ℓ = 0.001 are shown in Figs. 8–10. Varying the temperature and the bulk velocity affects the temperature structures and the emergent spectra of the models in the qualitatively similar way as the inclination angle of the bulk velocity vector made (Figs. 8 and 9). Such a similarity arises because fast particles penetrate into different atmospheric depths when we vary their temperature or the bulk velocity. And, therefore, the ratio between the energy, which is released in the optically thick, heated slab and in the surface overheated layer also changes. Indeed, for higher temperature, the fraction of the particles with high individual velocities increases as well. They penetrate into deeper layers, and the relative contribution of the optically thick heated layer also increases. It becomes hotter, whereas the temperature of the upper overheated layers decreases. The relative contributions of the blackbody and the Comptonized components also change accordingly (Fig. 8). In a similar way the particles with the larger bulk velocities penetrate deeper (Fig. 9) and heat the optically thick part of the atmosphere to higher temperatures. On the other hand, the accretion heating rate of the very surface layers is larger for the smaller bulk velocities.
Fig. 8. Emergent energy spectra (top panel) and temperature structures (bottom panel) of four low luminosity (ℓ = 0.001) helium model atmospheres heated by αparticles with luminosity ℓ_{a} = 0.05 and parameters η = 0.25 and Ψ = 60°. The presented models have different temperature parameter χ = 0.1, 0.2, 0.5, and 0.7. The fit to the spectrum above 5 keV with an exponential cutoff power law with Γ = 2.15, E_{cut} = 85 keV is shown for χ = 0.1. The dotted curve shows the blackbody with temperature kT = 1 keV. 
Fig. 9. Emergent spectrum (top panel) and temperature structure (bottom panel) of the fiducial model (with η = 0.75) in comparison with those of the models with smaller bulk velocities of the accreting particles (η = 0.50 and 0.25). The fit to the spectrum above 5 keV with an exponential cutoff power law with Γ = 2.25, E_{cut} = 80 keV is shown for η = 0.25. The dotted curve shows the blackbody with temperature kT = 1.1 keV. 
Fig. 10. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with different accretion luminosities (ℓ_{a} = 0.01, 0.03, 0.05, 0.1, 0.3, and 0.5). The fit to the spectrum above 15 keV with an exponential cutoff power law with Γ = 2.25, E_{cut} = 42 keV is shown for ℓ_{a} = 0.5. 
We note that the dependence of the accretion heated model atmospheres on the temperature of the fast particles was investigated before by Deufel et al. (2001). They showed that fast protons of lower temperature heat the overheated upper layers to the higher temperatures, whereas the temperature of the isothermal inner slab is lower than for the corresponding case of hightemperature protons (see their Fig. 3). As a result the emergent spectra of the models heated by lowtemperature protons are harder. Our models (see Fig. 8) confirm these results.
The models with varying mass accretion rate have identical dependence of the heating rate over the depth, with the normalizations directly proportional to ṁ. The depth of the transition to the overheated parts of the atmosphere and its temperature are determined by the accretion luminosity ℓ_{a}. Higher accretion luminosity leads to higher temperatures of the heated layers, both of the overheated upper layers and of the isothermal region, and to a smoother transition from the isothermal region to the overheated layers. The strength of the hard Comptonized component also depends on ℓ_{a}. At small accretion luminosities, the blackbody component dominates and the hard tail in the emergent spectra is relatively weak. At high accretion luminosities, L_{a} ≫ L, the emergent spectra are dominated by the hard Comptonized component. These results confirm those obtained before by Alme & Wilson (1973) and Deufel et al. (2001). They demonstrated that a rise of the accretion luminosity leads to a harder emergent spectrum and a higher temperature of the overheated layers and the isothermal slabs. This fact is especially clear from Fig. 4 in Deufel et al. (2001).
The influence of the accretion heating decreases when the intrinsic stellar luminosity increases (Fig. 11). The spectrum of the model with the lowest intrinsic luminosity ℓ = 0.001 is completely determined by the heated part of the atmosphere whereas the spectrum of the model with the highest intrinsic luminosity ℓ = 0.5 differs from the spectrum of the undisturbed model very little.
Fig. 11. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with different intrinsic luminosities (ℓ = 0.5, 0.1, and 0.001), but the same accretion luminosity ℓ_{a} = 0.05. The corresponding models for the undisturbed atmosphere are shown by the dashed curves. 
We also investigated the influence of the parameters for the case with the high intrinsic luminosity ℓ = 0.5. The qualitative dependences for the temperature structures are the same as they were described for the models with the low intrinsic luminosity ℓ = 0.001. However, the spectra are less sensitive to the variations of the parameters (see Fig. 12).
Fig. 12. Same as Fig. 9 but for higher intrinsic luminosity ℓ = 0.5. At the bottom panel, the curves are labeled with the corresponding values of parameter η. The spectrum and the temperature structure of the undisturbed model are also shown by the dashed curves. 
3.3. Atmospheres with the solar composition
The accreting matter in LMXBs has helium rich composition only in the ultracompact systems (for instance, in 4U 1820−30). Most of the systems have solar chemical composition of the accreting matter. Therefore, we have to consider how a change in the chemical composition affects the results. As can be seen from Eq. (31), protons and αparticles decelerate the same way (Bildsten et al. 1992), because the velocity derivative depends on the combination Z^{2}/A only. Here we consider only these two types of fast particles and ignore the nuclei of heavy elements. The solar abundance model atmospheres heated by the solar mix of the protons and αparticles are similar to the pure helium models (Fig. 13).
Fig. 13. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with solar chemical abundance and three intrinsic luminosities, ℓ = 0.5, 0.1 and 0.001, but with the same parameters of the accreting particles (ℓ_{a} = 0.05, η = 0.75, χ = 0.2, Ψ = 60°). The spectrum and the temperature structure of the undisturbed models are also shown by the dashed curves. 
However, there is some qualitative difference in the emergent spectra. The spectra of the lowluminosity undisturbed model atmospheres show absorption edges. The photoionisation edge of hydrogenlike iron is the most prominent among them (see model spectra with ℓ = 0.001 and 0.1 in Fig. 13). The opacity at photon energies above the photoionisation threshold is larger than that below the threshold. Therefore photons above the threshold escape from smaller column densities, where the temperature is lower. This is the reason why absorption edges exist. In the heated atmospheres, iron is more ionized and the jump in the opacity across the photoionisation threshold is smaller. Furthermore, the atmosphere is nearly isothermal in the region where the photons around the edge are formed (see bottom panel of Fig. 14). As a result, the spectra of the accretionheated atmospheres are almost featureless (see top panel of Fig. 14).
Fig. 14. Magnification of Fig. 13 but for one intrinsic luminosity ℓ = 0.1 only. The formation depths of the photons with energies of 8 and 10 keV are indicated by thin vertical lines for the heated and the undisturbed atmospheres. 
3.4. Application to Xray bursting neutron stars
Often the persistent spectra of Xray bursting NSs are subtracted from the spectra detected during Xray bursts for obtaining “pure” burst spectra (f_{E′}) (see, however, Worpel et al. 2013; Worpel et al. 2015; Degenaar et al. 2016; Kajava et al. 2017). These spectra are well fitted (Galloway et al. 2008) with a blackbody
where the normalization K and the color temperature T_{BB} are the fitting parameters. On the other hand, the observed flux is determined by the model atmosphere flux F_{E}(T_{eff})
where D is the distance to the source, and z is the gravitational redshift,
Here R_{S} = 2GM/c^{2} is the Schwarzschild radius. We note also that E = E′(1 + z). The model spectra of hot NS atmospheres are close to diluted blackbodies because of the effective interaction between plasma and radiation by Compton scattering (London et al. 1986; Lapidus et al. 1986). They can be fitted as
with two fitting parameters, the color correction factor f_{c} and the dilution factor . Both fitting parameters depend on the relative luminosity of the NS ℓ = L/L_{Edd}. Therefore, the observed spectrum will be defined as follows
or
Now we can determine the parameters K and T_{BB} from comparison of Eqs. (44) and (49):
and
It is clear that the observed normalization K can change with the observed bolometric flux at the cooling burst phases only because of the changes in the dilution factor w. Therefore, K has to change in a similar fashion as the dilution factor w varies with the relative luminosity ℓ. If we account for deviation of the bolometric flux from the observed blackbody flux F_{BB}, the observed relation K − F_{BB} should be fitted with the theoretical curve . Here is the bolometric correction, as was shown by Suleimanov et al. (2017b). From such fitting two parameters, the observed Eddington flux F_{Edd, ∞} and the NS observed solid angle (R(1 + z)/D)^{2} can be found, and NS mass and radius can be limited. This approach is known as the cooling tail method, see details in Suleimanov et al. (2011b, 2017b) and Poutanen et al. (2014). This method would be correct if a passively cooling NS were a correct model for the evolution of the burst spectra in the cooling phase. Indeed, the normalization of the spectra of bursts taking place during the hard persistent states evolve according to this model at a relatively high luminosity ℓ > 0.2 − 0.7. This allowed us to obtain constraints on the NS masses and radii (Suleimanov et al. 2011a, 2017b; Poutanen et al. 2014; Nättilä et al. 2016).
However, the observed dependences K − F_{BB} deviate from the model curves at low ℓ. The amplitude of the deviations depends on the relative persistent luminosity of the system before the Xray burst. Compare, for instance, small deviations of the model curves from the observed ones for the systems with low persistent emission, L_{per} ≈ 0.01 L_{Edd} (Nättilä et al. 2016), and the large deviations for the system 4U 1820−30 (Suleimanov et al. 2017b) with a relatively high persistent luminosity, L_{per} ≈ 0.04 − 0.07 L_{Edd}. We suggest that these deviations are determined by additional heating of the NS atmospheres by the accretion flow. Let us evaluate the effect of such a heating on the theoretical dependences to study this hypothesis.
To obtain the model spectra we use the same approach as used for the analysis of the observed spectra of Xray bursting NSs (see, e.g., Galloway et al. 2008). We assume that the spectrum of the heated model atmosphere with the lowest intrinsic luminosity ℓ = 0.001 represents the spectrum of the persistent emission from an accreting NS. We then subtract this spectrum from the spectra corresponding to the higher intrinsic luminosities to mimic pure burst model spectra. The examples of such residual spectra are presented in Fig. 15. The residual spectra can be fitted by diluted blackbodies even better than the spectra of undisturbed model atmospheres and their color temperatures are higher than the corresponding color temperatures of the spectra of undisturbed model atmospheres. We also note, that the absorption edge clearly seen in the model spectrum of the solarabundance undisturbed atmosphere completely disappears in the residual spectrum (Fig. 15, bottom panel). In the first approximation the reduced spectra are slightly dependent on the parameters of the accretion flow (Fig. 16).
Fig. 15. Residual spectra (solid curves) of the heated model atmospheres with parameters ℓ = 0.1, ℓ_{a} = 0.05, χ = 0.2, η = 0.75, and Ψ = 60° for two chemical compositions, pure helium (top panel) and solar abundance (bottom panel). The spectra of undisturbed atmospheres are shown by dashed curves. The bestfit diluted blackbody spectra to the residual spectra are shown by the dotted curves, and the blackbodies of the effective temperatures are shown by the pink dashdotted curves. 
Fig. 16. Residual spectra (solid curves) of the heated helium atmosphere model with parameters ℓ = 0.5, ℓ_{a} = 0.05, χ = 0.2, and Ψ = 60° for three different values of η = 0.75, 0.50, and 0.25 (see Fig. 12). The spectrum of an undisturbed atmosphere is shown by the blue dashed curve and the blackbody of the effective temperature is shown by the pink dashdotted curve. 
Using the grid of residual model spectra in the intrinsic luminosity range ℓ from 0.001 to 0.5 we computed the color correction and dilution factors for pure He and solar abundance atmospheres (see Figs. 17 and 18). We do not consider luminosities above ℓ = 0.5, because at those luminosities the radiation pressure should affect significantly the dynamics of the accreted fast particles. It is also possible that during Xray bursts the accretion luminosity smoothly increases when the burst luminosity decreases. We leave the models with high intrinsic luminosities for future investigations.
Fig. 17. Color correction (top panel) and dilution (bottom panel) factors as functions of the corrected relative luminosity. They are computed using residual spectra of the heated helium model atmospheres for three accretion rates, ℓ_{a} = 0.05 (dark blue), 0.03 (red), and 0.01 (light blue). Other parameters are χ = 0.2, η = 0.75, and Ψ = 60°. The corresponding factors computed using spectra of undisturbed atmospheres are also shown by black curves. 
Fig. 18. Same as in Fig. 17, but for atmospheres of solar composition. Here the accretion flow temperature parameter is χ = 0.3. 
We considered three values of the relative accretion luminosities, ℓ_{a} = 0.05, 0.03 and 0.01. We took the accretion flow parameters to be η = 0.75, Ψ = 60°, and χ = 0.2 for pure helium models, and χ = 0.3 for the solar abundance models. The residual spectra were fitted with the diluted blackbody spectra in the blueshifted observational energy range of the RXTE observatory, (3 − 20)(1 + z) keV. Here the gravitational redshift is 1 + z = 1.27 for the accepted NS parameters and logg = 14.3, see Sect. 3. According to expectations, the color correction factors are larger for the heated model atmospheres, and their values depend on the mass accretion rate. We note that the prominent dip at ℓ ≈ 0.1 in the theoretical curve computed for the undisturbed solarabundance atmospheres becomes weaker for slightly heated model atmospheres (ℓ_{a} = 0.01) and completely disappears in the curves computed for higher accretion luminosities.
In the previous work (Suleimanov et al. 2017a) we considered Xray bursts occurring in a hard state of the ultracompact system 4U 1820−30 which likely accretes helium from a white dwarf companion. We applied the direct cooling tail method of Suleimanov et al. (2017b) to obtain constraints on the NS mass and radius. The persistent luminosities before the considered bursts were high enough (0.047–0.063 L_{Edd}) and the observed dependences K − F_{BB} deviated from the model curves at ℓ < 0.6 (see Fig. 19). Therefore, we analyzed only the highluminosity part of the observed data. We suggested that deviations of the theoretical curve from the data are likely caused by accretion during the later stages of the cooling tail. Taking the accretion luminosity of 5% of the Eddington value (ℓ_{a} = 0.05), the model curves for the heated pure helium atmospheres shift down sufficiently to pass through the data points at ℓ < 0.5. At the relative luminosities higher than ℓ ≈ 0.3, i.e. F_{BB} ≈ 0.2 × 10^{−7} erg s^{−1} cm^{−2}, the data lie slightly above the model curve. We suggest that when the burst luminosity drops from ℓ ≈ 0.6 to ≈0.3, the accretion rate increases causing a smooth shift of the blackbody normalization K from a theoretical curve with no accretion ℓ_{a} = 0 to a curve corresponding to ℓ_{a} = 0.05. There are a few possible reasons why the accretion restarts again after some time. First, the inner parts of the accretion flow can be destroyed during the photospheric radius expansion phase. The typical time of the photosphere to settle down to the NS surface can be much shorter than the viscous time for the inner accretion flow. Therefore, the accretion restarts again gradually with some time gap. Second, is the influence of radiation pressure. The radiation field has two effects on the particle dynamics near the NS (Miller & Lamb 1993). At high luminosity, the direct radial radiation force may reduce the radial velocity decreasing mass accretion rate, while at low luminosity the radial velocity may increase as a result of the angular momentum loss due to the PoyntingRobertson drag. Finally, the burst radiation may cool the accretion flow making it thinner. As a result, only a relatively small part of the NS surface is affected by accretion. It may well be possible that all three reasons play a role.
Fig. 19. Comparison of the observed dependences K − F_{BB} obtained for selected Xray bursts of 4U 1820−30 (see Suleimanov et al. 2017a) with the model curves computed for undisturbed (blue curves) and heated (red curves) pure helium atmospheres. The axes of the fluxes are shown in linear (top panel) and logarithmic (bottom panel) scales. The fitting parameters F_{Edd, ∞} = 0.6 × 10^{−7} erg s^{−1} cm^{−2} and Ω = 500 (km/10 kpc)^{−2} are the same for both model curves. 
4. Summary
In this paper we presented a computational method to construct NS model atmospheres heated by accreted particles. Our results confirm previous findings (e.g., Alme & Wilson 1973; Zampieri et al. 1995; Deufel et al. 2001) that fast heavy particles (protons, αparticles, or their mix) penetrate down to some depth in the atmosphere and release most of their kinetic energy in a relatively narrow layer due to Coulomb interaction with electrons. The heating by accretion is balanced by the cooling due to mainly freefree emission and Comptonization. The relatively dense optically thick region, where the cooling is dominated by freefree emission, is only slightly hotter than the corresponding layer of an undisturbed atmosphere. It forms a quasiisothermal transition region between the inner atmosphere unaffected by accretion heating and the hot (∼10^{8} − 10^{9} K) rarefied upper layer. The upper hot layers cool mainly by Compton scattering. The width of the transition region and the temperature of the upper layer depend mainly on the accretion luminosity L_{a} and the intrinsic NS luminosity L.
The emergent spectrum of an accretionheated NS atmosphere is wider than the spectrum of the corresponding undisturbed atmosphere and can be roughly represented as a sum of two components: a blackbodylike radiation from the transition region and Comptonized (cutoff powerlawlike) spectrum of the hot upper layers. The relative contribution of the components is determined by the ratio of the intrinsic luminosity L to the accretion luminosity L_{a}. If L_{a} ≫ L, the total spectrum is dominated by the Comptonized spectrum of the upper, hot rarefied layers. This kind of spectra are similar to the observed spectra of LMXBs in their low hard spectral states, as mentioned by Deufel et al. (2001), and their slope depends on the input parameters of the accreted particles: velocity, temperature, and the penetration angle. We note, however, that all computed spectra have photon indexes Γ > 2. The observed spectra of LMXBs in their low hard spectral states often have harder spectra with Γ < 2. It means that the contribution of the accretion flow emission to the total observed spectra may not be negligible. In the opposite case, when L_{a} ≪ L, the total spectra are very similar to the spectra of the undisturbed models with harder Wien tails and increased flux at low photon energies. In that case, the final spectra depend very little on the properties of accreted particles.
The atmospheres heated by material of solar abundance have an important qualitative feature in comparison to pure helium atmospheres heated by αparticles only. The spectra of the relatively lowluminosity (L ≪ L_{Edd}) undisturbed model atmospheres show significant absorption edges due to photoionisation of hydrogenlike iron. These edges disappear in the spectra of the accretionheated atmospheres because of a significant change in the temperature structure.
We also investigated the influence of the accretion heating on the model curves and , which are used in the (direct) cooling tail method. As it was expected the colorcorrection factors f_{c} are larger for the heated atmospheres in comparison with the undisturbed atmospheres, and the dilution factors w are smaller. The model curve computed for heated helium atmospheres (with ℓ_{a} = 0.05) is well fitted to the lowluminosity part of the observed data K − F_{BB} obtained for Xray bursts taken place during the hard state of the heliumaccreting system 4U 1820−30 (see Suleimanov et al. 2017a).
The model curves computed for the heated solarabundance atmospheres do not have a dip at ℓ ≈ 0.1, which is clearly seen in the model curves computed for the undisturbed solarabundance atmospheres. This fact may have important implications for interpretation of the Xray data on bursting NSs that accrete gas of solar abundance, for example, the Clocked Burster GS 1826−24 (see discussion in Zamfir et al. 2012). We plan to apply the method described here for interpretation of the spectral evolution of that and other bursters in a followup paper.
We note here that the bursts occurring in the soft, highaccretionrate state never show spectral evolution consistent with theoretical prediction (Kajava et al. 2014). Therefore, this kind of burst cannot be used to determine NS parameters at all.
Acknowledgments
The authors acknowledge A. M. Beloborodov for useful discussions and an anonymous referee for the very insightful comments. This research has been supported by the grant 14.W03.31.0021 of the Ministry of Education and Science of the Russian Federation. V. F. S. also thanks Magnus Ehrnrooth foundation and Deutsche Forschungsgemeinschaft (DFG) for financial support (grant WE 1312/511). We thank the German Academic Exchange Service (DAAD, project 57405000) and the Academy of Finland (project 317552) for travel grants. This work benefited from discussions at the BERN18 Workshop supported by the National Science Foundation under Grant No. PHY1430152 (JINA Center for the Evolution of the Elements).
References
 Alme, M. L., & Wilson, J. R. 1973, ApJ, 186, 1015 [NASA ADS] [CrossRef] [Google Scholar]
 Barret, D., Olive, J. F., Boirin, L., et al. 2000, ApJ, 533, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, S. T., & Budker, G. I. 1956, Dokl. Adac. Nauk SSSR, 107, 807 [Google Scholar]
 Bildsten, L., Salpeter, E. E., & Wasserman, I. 1992, ApJ, 384, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Burke, M. J., Gilfanov, M., & Sunyaev, R. 2017, MNRAS, 466, 194 [NASA ADS] [CrossRef] [Google Scholar]
 Degenaar, N., & Suleimanov, V. 2018, ArXiv eprints [arXiv:1806.02833] [Google Scholar]
 Degenaar, N., Koljonen, K. I. I., Chakrabarty, D., et al. 2016, MNRAS, 456, 4256 [NASA ADS] [CrossRef] [Google Scholar]
 Deufel, B., Dullemond, C. P., & Spruit, H. C. 2001, A&A, 377, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Done, C., Gierliński, M., & Kubota, A. 2007, A&ARv, 15, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge: Cambridge University Press) [Google Scholar]
 Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Homan, J., Fridriksson, J. K., Wijnands, R., et al. 2014, ApJ, 795, 131 [Google Scholar]
 Inogamov, N. A., & Sunyaev, R. A. 1999, Astron. Lett., 25, 269 [NASA ADS] [Google Scholar]
 Kajava, J. J. E., Nättilä, J., Latvala, O.M., et al. 2014, MNRAS, 445, 4218 [NASA ADS] [CrossRef] [Google Scholar]
 Kajava, J. J. E., SánchezFernández, C., Kuulkers, E., & Poutanen, J. 2017, A&A, 599, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lapidus, I. I., Sunyaev, R. A., & Titarchuk, L. G. 1986, Sov. Astron. Lett., 12, 383 [NASA ADS] [Google Scholar]
 Larkin, A. I. 1960, Sov. Phys. JETP, 37, 186 [Google Scholar]
 London, R. A., Taam, R. E., & Howard, W. M. 1986, ApJ, 306, 170 [Google Scholar]
 Mahadevan, R., & Quataert, E. 1997, ApJ, 490, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 1993, ApJ, 413, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Nagirner, D. I., & Poutanen, J. 1993, A&A, 275, 325 [NASA ADS] [Google Scholar]
 Nagirner, D. I., & Poutanen, J. 1994, Astrophys. Space Phys. Rev., 9, 1 [NASA ADS] [Google Scholar]
 Nättilä, J., Suleimanov, V. F., Kajava, J. J. E., & Poutanen, J. 2015, A&A, 581, A83 [NASA ADS] [EDP Sciences] [Google Scholar]
 Nättilä, J., Steiner, A. W., Kajava, J. J. E., Suleimanov, V. F., & Poutanen, J. 2016, A&A, 591, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nättilä, J., Miller, M. C., Steiner, A. W., et al. 2017, A&A, 608, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelson, R. W., Salpeter, E. E., & Wasserman, I. 1993, ApJ, 418, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., & Veledina, A. 2014, Space Sci. Rev., 183, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J., Nättilä, J., Kajava, J. J. E., et al. 2014, MNRAS, 442, 3777 [NASA ADS] [CrossRef] [Google Scholar]
 Romani, R. W. 1987, ApJ, 313, 718 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., Revnivtsev, M., & Werner, K. 2011a, ApJ, 742, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2011b, A&A, 527, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V. F., Poutanen, J., Klochkov, D., & Werner, K. 2016, Eur. Phys. J. A, 52, 20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suleimanov, V. F., Kajava, J. J. E., Molkov, S. V., et al. 2017a, MNRAS, 472, 3905 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V. F., Poutanen, J., Nättilä, J., et al. 2017b, MNRAS, 466, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Turolla, R., Zampieri, L., Colpi, M., & Treves, A. 1994, ApJ, 426, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Wijnands, R., Degenaar, N., Armas Padilla, M., et al. 2015, MNRAS, 454, 1371 [NASA ADS] [CrossRef] [Google Scholar]
 Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Worpel, H., Galloway, D. K., & Price, D. J. 2015, ApJ, 801, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Zamfir, M., Cumming, A., & Galloway, D. K. 2012, ApJ, 749, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Zampieri, L., Turolla, R., Zane, S., & Treves, A. 1995, ApJ, 439, 849 [Google Scholar]
 Zane, S., Turolla, R., & Treves, A. 1998, ApJ, 501, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Zane, S., Turolla, R., & Treves, A. 2000, ApJ, 537, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Zavlin, V. E., Pavlov, G. G., & Shibanov, Yu A 1996, A&A, 315, 141 [NASA ADS] [Google Scholar]
 Zel’dovich, Y. B., & Shakura, N. I. 1969, Sov. Astron., 13, 175 [NASA ADS] [Google Scholar]
All Figures
Fig. 1. Coordinate systems used in the paper. We also show the momentum of a single particle and its components with the corresponding angles. 

In the text 
Fig. 2. Emergent energy spectra (top panel) and temperature structures (bottom panel) for four lowluminosity (ℓ = 0.001) helium model atmospheres heated by αparticles with the accretion luminosity ℓ_{a} = 0.05 and parameters η = 0.75 and χ = 0.2. The models correspond to different incoming angles Ψ = 0°, 30°, 60°, and 85°. The emergent spectrum of the most inclined (Ψ = 85°) particle flow is well fitted at photon energies higher than 5 keV with an exponential cutoff powerlaw model with photon index Γ = 2.45 and cutoff energy E_{cut} = 85 keV (dashed curve). The blackbody spectrum of temperature kT = 1.1 keV is also shown by the dotted curve. 

In the text 
Fig. 3. Bolometric flux (top panel) and heating rate (bottom panel) as functions of the column density for the models presented in Fig. 2. 

In the text 
Fig. 4. Emergent specific intensity spectra of lowluminosity (ℓ = 0.001) helium atmosphere model heated by αparticles with accretion luminosity ℓ_{a} = 0.05 at three zenith angles of 27°.5 (dotted curve), 60° (solid curve), and 83°.5 (dashed curve). The parameters of the model are Ψ = 85°, η = 0.75 and χ = 0.2. 

In the text 
Fig. 5. Relative ram pressure acceleration (top panel) and the relative ram pressure (bottom panel) as functions of the column density for the models presented in Fig. 2. 

In the text 
Fig. 6. Relative heat conduction flux (top panel) and heat conduction flux derivative (bottom panel) as functions of the column density for the models presented in Fig. 2. 

In the text 
Fig. 7. Atmosphere density dependence on the geometrical height for the models presented in Fig. 2. 

In the text 
Fig. 8. Emergent energy spectra (top panel) and temperature structures (bottom panel) of four low luminosity (ℓ = 0.001) helium model atmospheres heated by αparticles with luminosity ℓ_{a} = 0.05 and parameters η = 0.25 and Ψ = 60°. The presented models have different temperature parameter χ = 0.1, 0.2, 0.5, and 0.7. The fit to the spectrum above 5 keV with an exponential cutoff power law with Γ = 2.15, E_{cut} = 85 keV is shown for χ = 0.1. The dotted curve shows the blackbody with temperature kT = 1 keV. 

In the text 
Fig. 9. Emergent spectrum (top panel) and temperature structure (bottom panel) of the fiducial model (with η = 0.75) in comparison with those of the models with smaller bulk velocities of the accreting particles (η = 0.50 and 0.25). The fit to the spectrum above 5 keV with an exponential cutoff power law with Γ = 2.25, E_{cut} = 80 keV is shown for η = 0.25. The dotted curve shows the blackbody with temperature kT = 1.1 keV. 

In the text 
Fig. 10. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with different accretion luminosities (ℓ_{a} = 0.01, 0.03, 0.05, 0.1, 0.3, and 0.5). The fit to the spectrum above 15 keV with an exponential cutoff power law with Γ = 2.25, E_{cut} = 42 keV is shown for ℓ_{a} = 0.5. 

In the text 
Fig. 11. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with different intrinsic luminosities (ℓ = 0.5, 0.1, and 0.001), but the same accretion luminosity ℓ_{a} = 0.05. The corresponding models for the undisturbed atmosphere are shown by the dashed curves. 

In the text 
Fig. 12. Same as Fig. 9 but for higher intrinsic luminosity ℓ = 0.5. At the bottom panel, the curves are labeled with the corresponding values of parameter η. The spectrum and the temperature structure of the undisturbed model are also shown by the dashed curves. 

In the text 
Fig. 13. Emergent spectra (top panel) and temperature structures (bottom panel) of the models with solar chemical abundance and three intrinsic luminosities, ℓ = 0.5, 0.1 and 0.001, but with the same parameters of the accreting particles (ℓ_{a} = 0.05, η = 0.75, χ = 0.2, Ψ = 60°). The spectrum and the temperature structure of the undisturbed models are also shown by the dashed curves. 

In the text 
Fig. 14. Magnification of Fig. 13 but for one intrinsic luminosity ℓ = 0.1 only. The formation depths of the photons with energies of 8 and 10 keV are indicated by thin vertical lines for the heated and the undisturbed atmospheres. 

In the text 
Fig. 15. Residual spectra (solid curves) of the heated model atmospheres with parameters ℓ = 0.1, ℓ_{a} = 0.05, χ = 0.2, η = 0.75, and Ψ = 60° for two chemical compositions, pure helium (top panel) and solar abundance (bottom panel). The spectra of undisturbed atmospheres are shown by dashed curves. The bestfit diluted blackbody spectra to the residual spectra are shown by the dotted curves, and the blackbodies of the effective temperatures are shown by the pink dashdotted curves. 

In the text 
Fig. 16. Residual spectra (solid curves) of the heated helium atmosphere model with parameters ℓ = 0.5, ℓ_{a} = 0.05, χ = 0.2, and Ψ = 60° for three different values of η = 0.75, 0.50, and 0.25 (see Fig. 12). The spectrum of an undisturbed atmosphere is shown by the blue dashed curve and the blackbody of the effective temperature is shown by the pink dashdotted curve. 

In the text 
Fig. 17. Color correction (top panel) and dilution (bottom panel) factors as functions of the corrected relative luminosity. They are computed using residual spectra of the heated helium model atmospheres for three accretion rates, ℓ_{a} = 0.05 (dark blue), 0.03 (red), and 0.01 (light blue). Other parameters are χ = 0.2, η = 0.75, and Ψ = 60°. The corresponding factors computed using spectra of undisturbed atmospheres are also shown by black curves. 

In the text 
Fig. 18. Same as in Fig. 17, but for atmospheres of solar composition. Here the accretion flow temperature parameter is χ = 0.3. 

In the text 
Fig. 19. Comparison of the observed dependences K − F_{BB} obtained for selected Xray bursts of 4U 1820−30 (see Suleimanov et al. 2017a) with the model curves computed for undisturbed (blue curves) and heated (red curves) pure helium atmospheres. The axes of the fluxes are shown in linear (top panel) and logarithmic (bottom panel) scales. The fitting parameters F_{Edd, ∞} = 0.6 × 10^{−7} erg s^{−1} cm^{−2} and Ω = 500 (km/10 kpc)^{−2} are the same for both model curves. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.