© ESO 2020

1. Introduction

Rotation-powered millisecond pulsars (RMPs), also known as radio millisecond pulsars, are rapidly rotating neutron stars (NSs), which are believed to have been spun-up due to the accretion from a companion star (Radhakrishnan & Srinivasan 1982; Alpar et al. 1982). After transiting from the accretion-powered to the rotation-powered phase, the NS polar caps are, as opposed to accretion, heated by the bombardment of electrons and positrons from a magnetospheric return current (see e.g. Zel’dovich & Shakura 1969; Alme & Wilson 1973; Ruderman & Sutherland 1975; Zampieri et al. 1995). The rotation of the resulting hotspots gives rise to X-ray pulsations.

Pulse profiles of millisecond pulsars carry information about the NS mass and radius (Poutanen & Gierliński 2003; Poutanen & Beloborodov 2006; Morsink et al. 2007; Miller & Lamb 2015). These pulses can be modelled using a general relativistic ray-tracing framework (e.g. Poutanen & Beloborodov 2006; Morsink et al. 2007; Nättilä & Pihajoki 2018; Salmi et al. 2018) to constrain the unknown equation of state (EOS) of cold dense matter (Lattimer 2012; Özel & Freire 2016; Watts et al. 2016, 2019). In accreting millisecond pulsars, the largest uncertainty in modelling the pulse profiles comes from the poorly known emission pattern of the hotspots (Poutanen & Gierliński 2003) as well as from the variability of the pulse profiles (e.g. Ibragimov & Poutanen 2009). In this respect, the RMPs are potentially better targets, as their pulse profiles are stable on the time scales of years and therefore can be accurately measured, for instance, with the Neutron star Interior Composition ExploreR (NICER) mission (Gendreau et al. 2016). Their emission pattern was also expected to be well-described by a hydrogen atmosphere model (Bogdanov 2013, 2016), which has also been used recently to model the NICER data (Bogdanov et al. 2019; Miller et al. 2019; Riley et al. 2019). However, the exact spectral energy distribution and the angular emission pattern of radiation emitted by the hotspots affect the pulse shape, and they are therefore important to be accurately modelled (Salmi et al. 2018). They depend on the details of the return-current particle distribution and the resulting atmosphere structure, which are currently unknown (Bauböck et al. 2019).

Atmospheres of the NSs have been previously modelled in many different studies. Using a plane-parallel atmosphere model for RMPs in radiative equilibrium, the angular and energy distribution of the radiation has been described, for example, by Bogdanov (2013; see also e.g. Zavlin et al. 1996; Heinke et al. 2006; Ho & Heinke 2009; Haakonsen et al. 2012). A similar model, but with an exact formulation for Compton scattering, is studied in Salmi et al. (2019). However, the assumption of the radiative equilibrium is only valid if the in-falling pair plasma heating the polar caps deposits its energy at very deep layers. It is therefore necessary to study the additional physical processes for atmospheres heated, in their different layers, by a beam of bombarding particles. Atmospheres heated from the top have previously been modelled by González-Caniulef et al. (2019) for strongly magnetised NSs (in the grey approximation) and for accreting NSs by Deufel et al. (2001) and Suleimanov et al. (2018). In the case of RMPs, the effects of stopping the magnetospheric particles have been studied by Bauböck et al. (2019), although using a simplified atmosphere model. In this work, we aim to study the return-current-heating effects using the full atmosphere model including the exact treatment of Compton scattering, which is expected to be important for the energy balance of the heated outer layers of the atmosphere and for the observed spectrum, in particular, its angular distribution.

The remainder of this paper is structured as follows. In Sect. 2, we explain the theoretical framework, including how we modelled the return-current-heated NS atmosphere, emphasising the differences to the existing models. In Sect. 3, we calculate NS atmospheres for a set of model parameters, especially for different energies or energy distributions of the bombarding particles, and compare the results with each other and to the non-heated model. We discuss the reasons for the differences and implications for the NICER parameter constraints in Sect. 4. We conclude in Sect. 5.

2. The model

2.1. Main equations

We consider steady-state plane-parallel atmospheres with additional heating in the surface layers caused by the magnetospheric return current. According to Arons (1981) and Harding & Muslimov (2002), we could expect that the particles penetrating NS atmosphere would primarily be positrons. However, we still consider the loss rates only for electrons, because they lose their energy in the same way to positrons and have a simpler form of energy-loss equations (Berger et al. 1984; Bauböck et al. 2019). We also assume the NS atmosphere to consist purely of hydrogen, which is justified by a short gravitational stratification timescale (Alcock & Illarionov 1980; Zavlin et al. 2002), leaving only the lightest elements to the layers which determine the properties of the escaping radiation.

The code that we use in atmosphere modelling is a modified version from Nättilä et al. (2015), which is based on Suleimanov et al. (2012) and Kurucz (1970). Our model has a set of equations very similar to those presented in Suleimanov et al. (2018). However, instead of considering the energy loss of accreted protons, we consider the case for relativistic return current electrons. The heating rate of the NS atmosphere is expressed in terms of relative luminosity $l=F/F_{\mathrm{Edd}}=T_{\mathrm{eff}}^4/T_{\mathrm{Edd}}^4$, where

$$\begin{array}{c}\hfill F_{\mathrm{Edd}}=\frac{\mathit{gc}}{{\kappa }_{\mathrm{e}}}={\sigma }_{\mathrm{SB}}{T}_{\mathrm{Edd}}^4\end{array}$$(1)

is the Eddington flux and T_Edd is the corresponding Eddington temperature, which is the maximum effective temperature on the NS surface for which local flux does not exceed the Eddington one for the given NS surface gravitational acceleration g (a free parameter of the model). Here σ_SB is the Stefan-Boltzmann constant, and the Eddington flux is evaluated using Thomson scattering opacity κ_e = σ_T/m_p ≈ 0.4 cm² g⁻¹, where σ_T is the Thomson cross-section and m_p is the proton mass.

The emergent bolometric flux is obtained as a sum of two components,

$$\begin{array}{c}\hfill F={\sigma }_{\mathrm{SB}}(T_{\mathrm{eff},\mathrm{i}}^4+T_{\mathrm{eff},\mathrm{h}}^4)={\sigma }_{\mathrm{SB}}{T}_{\mathrm{Edd}}^4(l_{\mathrm{i}}+l_{\mathrm{h}})=F_{\mathrm{i},0}+F_{\mathrm{h},0},\end{array}$$(2)

where the lower indices i and h refer to the intrinsic NS flux and that added by return-current-heating, respectively. In all of our models, we put l_i = l_h/100. The exact value is not expected to be important as long as l_h is much greater than l_i. The number flux (per unit area) of bombarding electrons is related to the heating rate as

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{e}}=\frac{F_{\mathrm{h},0}}{m_{\mathrm{e}}{c}^2(⟨{\gamma }_0⟩-1)}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1},\end{array}$$(3)

where ⟨γ₀⟩ is the average Lorentz factor of the in-falling particles and m_e is the electron (positron) mass.

The computation follows otherwise Suleimanov et al. (2012), but the absorption opacities and the electron number densities are obtained directly from the Kramers’ opacity law and the ideal gas law, since we assume a fully ionised hydrogen atmosphere with free-free opacity as the only source of absorption opacity. The contribution from return-current-heating is added in the hydrostatic equilibrium equation which reads as

$$\begin{array}{c}\hfill \frac{1}{\rho }\frac{\mathrm{d}P}{\mathrm{d}r}=-g-g_{\mathrm{ram}}+g_{\mathrm{rad}},\end{array}$$(4)

where P is the gas pressure, g_rad is the radiative acceleration (defined in Suleimanov et al. 2012), and g_ram is the additional ram pressure acceleration (Suleimanov et al. 2018). For monoenergetic in-falling particles the ram acceleration is

$$\begin{array}{c}\hfill g_{\mathrm{ram}}=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}\frac{\mathrm{d}}{\mathrm{d}m}(\gamma v)=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}c\frac{\gamma }{p}\frac{\mathrm{d}\gamma }{\mathrm{d}m},\end{array}$$(5)

where dm = −ρdr is the column density, ρ is the plasma density of the atmosphere, v = βc is the velocity of the particles, γ is the corresponding Lorentz factor, and $p=\sqrt{{\gamma }^2-1}$ is the dimensionless momentum at the given position m.

The energy balance equation is written as

$$\begin{array}{c}\hfill 2\pi {\displaystyle {\int }_0^{\infty }\mathrm{d}x{\int }_{-1}^{+1}[\sigma (x,\mu )+k(x)][I(x,\mu )-S(x,\mu )]\mathrm{d}\mu =-Q^{+},}\end{array}$$(6)

where k(x) is the ‘true’ absorption opacity, and σ(x, μ) is the electron scattering opacity. The general form of the local energy dissipation rate Q⁺ is (Suleimanov et al. 2018):

$$\begin{array}{c}\hfill Q^{+}={\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}{q}^{+}=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2\frac{\mathrm{d}}{\mathrm{d}m}(\gamma -1)=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2\frac{\mathrm{d}\gamma }{\mathrm{d}m},\end{array}$$(7)

where q⁺ is the specific particle deceleration.

The dependence of particle deceleration on depth dγ/dm fully determines the atmosphere model. To solve it, we can compute the energy loss rate along the radial direction r (which we assume to be the same as the displacement direction because the particles penetrate very closely to the normal direction) as

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t}=-\beta c\frac{\mathrm{d}E}{\mathrm{d}r}=-\beta c\frac{\mathrm{d}(m_{\mathrm{e}}{c}^2(\gamma -1))}{\mathrm{d}r}=\rho m_{\mathrm{e}}{c}^3\beta \frac{\mathrm{d}\gamma }{\mathrm{d}m}·\end{array}$$(8)

Another rather general expression for the energy losses can be written as

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t}=-n_{\mathrm{e}}{\sigma }_{\mathrm{T}}{m}_{\mathrm{e}}{c}^3\beta \mathrm{\Phi }(\gamma ),\end{array}$$(9)

where Φ is an effective cross-section in units of the Thomson cross-section, and n_e is the number density of the electrons in the atmosphere. These two equations then transform to the following differential equation for the Lorentz factor

$$\begin{array}{c}\hfill \frac{\mathrm{d}\gamma }{\mathrm{d}\tau }=-\mathrm{\Phi }(\gamma ),\end{array}$$(10)

where dτ = κ_edm is the Thomson optical depth. This differential equation can be solved numerically to obtain γ(m) with a given T(m), ρ(m), and n_e(m). The solution is then used to determine the local energy dissipation rate and additional acceleration due to ram pressure. We note that the solution γ(m) is largely independent of the atmospheric structure, except for a residual dependence of the cross-section Φ(γ) on density (Berger et al. 1984). The dependence of the cross-section on the particle Lorentz factor is considered in the next section.

2.2. Stopping of fast electrons

In a highly or moderately relativistic case, the energy loss of the electrons takes place via the bremsstrahlung radiation (dominating above 500 MeV or γ ≈ 1000), due to the Coulomb collisions with particles in the plasma, and due to excitation of collective modes of plasma (Langmuir waves). However, it could be that electrons with energies much larger than the thermal energy of the atmosphere electrons cannot excite Langmuir waves effectively. The total energy loss cross-section Φ(γ) (in Eq. (10)) is obtained by summing the cross-sections for all the different mechanisms.

An accurate approximation at all energies for the effective cross-section of electron-proton bremsstrahlung derived from the analytic solutions of Heitler (1954) using Born approximation in the non-relativistic and extreme relativistic limits is given by Haug (2004):

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{rad}}^{{\mathrm{e}}^{-}\mathrm{p}}(\gamma )& =\frac{3\alpha }{\pi }\frac{{\gamma }^3}{p^2+{\gamma }^2}\{\frac{ln(\gamma +p)}{\beta }-\frac{1}{3}\hfill \\ \hfill & +\frac{p^2}{{\gamma }^6}[\frac{2}{9}{\gamma }^2-\frac{19}{675}\gamma p^2-0.06\frac{p^4}{\gamma }]\},\hfill \end{array}$$(11)

where α is the fine structure constant. For the electron-electron bremsstrahlung cross-section, we use different analytic approximations at different electron energies, that are also presented in Haug (2004). For example, in case of ultra-relativistic energies (γ >  1000) we have

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{rad}}^{{\mathrm{e}}^{-}{\mathrm{e}}^{-}}(\gamma )\approx \alpha \frac{3}{2\pi }(\gamma -1)[ln(2\gamma )-\frac{1}{3}].\end{array}$$(12)

The energy loss of the electron due to the combined effect of Coulomb collisions with electrons (assuming that the electron velocity v is much larger than the thermal velocity v_th of the plasma electrons at that depth) and Langmuir waves, is given by Solodov & Betti (2008; see also Gould 1972):

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{C}+\mathrm{L}}(\gamma )& =\frac{3}{4{\beta }^2}[ln\left(\frac{(\gamma -1){\beta }^2}{2{ϵ}_{\mathrm{p}}^2}\right)+1\hfill \\ \hfill & +\frac{1}{8}{\left(\frac{\gamma -1}{\gamma }\right)}^2-\left(\frac{2\gamma -1}{{\gamma }^2}\right)ln2],\hfill \end{array}$$(13)

where ϵ_p = ℏω_p/m_ec², ${\omega }_{\mathrm{p}}=\sqrt{4\pi n_{\mathrm{e}}{e}^2/m_{\mathrm{e}}}$ is the plasma frequency. The effect of electron-proton Coulomb collisions can be neglected because the energy loss is inversely proportional to the mass of the target.

2.3. Return current with energy distribution

In case the in-falling electrons have some energy distribution, the total energy deposition rate is obtained by integrating the product of local dissipation rate and the distribution function of the particles over the initial Lorentz factor γ₀ as

$$\begin{array}{c}\hfill Q^{+}(m)={\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}f({\gamma }_0)q^{+}({\gamma }_0,m)\mathrm{d}{\gamma }_0,}\end{array}$$(14)

where we set the maximum to 10⁵ MeV (γ_max ≈ 2 × 10⁵), and the minimum γ_min is our free parameter. We use a power-law distribution of the return-current particle number over the Lorentz factor as suggested by observations of gamma-ray pulsars and simulations of pulsar magnetospheres (Harding & Muslimov 2001; Cerutti et al. 2016; Brambilla et al. 2018):

$$\begin{array}{c}\hfill f({\gamma }_0)=N{\gamma }_0^{-\delta },\end{array}$$(15)

where $N=1/{\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}{\gamma }_0^{-\delta }\mathrm{d}{\gamma }_0$ is the normalisation constant, and the exponent δ is a free parameter in our model. The average Lorentz factor for the incoming return-current beam is given as (if δ ≠ 2)

$$\begin{array}{c}\hfill ⟨{\gamma }_0⟩={\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}N{\gamma }_0^{1-\delta }\mathrm{d}{\gamma }_0=\frac{N({\gamma }_{\max }^{2-\delta }-{\gamma }_{\min }^{2-\delta })}{2-\delta }·}\end{array}$$(16)

If δ = 2, we have

$$\begin{array}{c}\hfill ⟨{\gamma }_0⟩=Nln\left(\frac{{\gamma }_{\max }}{{\gamma }_{\min }}\right).\end{array}$$(17)

For the energy deposition rate, we get now

$$\begin{array}{c}\hfill Q^{+}(m)=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}f({\gamma }_0)\frac{\mathrm{d}\gamma }{\mathrm{d}m}\mathrm{d}{\gamma }_0,}\end{array}$$(18)

where dγ/dm is a function of γ₀ and m. Integration over the column density should give us the total energy deposition rate per unit area:

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\infty }{Q}^{+}(m)\mathrm{d}m=F_{\mathrm{h},0}={\dot{N}}_{\mathrm{e}}(⟨{\gamma }_0⟩-1)m_{\mathrm{e}}{c}^2,}\end{array}$$(19)

which determines the number flux of electrons Ṅ_e appearing in Eq. (18).

The expression for the additional acceleration due to ram pressure created by the penetrating particles (used in the hydrostatic balance equation) becomes (see again Suleimanov et al. 2018)

$$\begin{array}{c}\hfill g_{\mathrm{ram}}=-{\dot{N}}_{\mathrm{e}}{m}_{\mathrm{e}}c{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}f({\gamma }_0)\frac{\mathrm{d}(\gamma \beta )}{\mathrm{d}m}\mathrm{d}{\gamma }_0.}\end{array}$$(20)

2.4. Computation of the atmosphere structure

We computed the model atmospheres using similar setup in energy and optical depth grids as in Suleimanov et al. (2012). The formal solution for radiation transfer equation was found using the short-characteristic method (Olson & Kunasz 1987) first in three angles, but finally with eleven angles in the end of the temperature iterations. The parabolic approximation of the solution was replaced with the linear approximation to avoid negative intensities. The full solution was found using the accelerated Λ-iteration method described in Suleimanov et al. (2012).

Compared to Suleimanov et al. (2012), the temperature corrections were modified so that the relative flux error (used in Avrett-Krook flux correction) was

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{F}}(m)=1-\frac{F_{\mathrm{i},0}+F_{\mathrm{h},0}-{\int }_0^m{Q}^{+}(m\prime )\mathrm{d}m\prime }{{\int }_0^{\infty }{F}_{\mathrm{x}}(m)\mathrm{d}x},\end{array}$$(21)

where x is the photon energy. In addition, the energy balance error, used in the temperature corrections for upper atmospheric layers, became

$$\begin{array}{cc}& {ϵ}_{\lambda }(m)=\frac{w}{2}{\displaystyle {\int }_0^{\infty }\mathrm{d}x{\int }_{-1}^{+1}[\sigma (x,\mu )+k(x)][I(x,\mu )-S(x,\mu )]\mathrm{d}\mu }\hfill \\ \hfill & +w{Q}^{+}(m),\hfill \end{array}$$(22)

where w is a weight used to adjust the correction at the beginning of the iterations.

The energy dissipation rate Q⁺(m) was calculated for each iteration using Eq. (18) for a distribution of particles or Eq. (7) for mono-energetic particles. When solving Eq. (10) to find γ(m) we assumed that the bombarding electrons have lost all their energy at the depth where their velocity drops below γ = 1.1 (as mentioned above, the equations given in Sect. 2.2 are only valid for fast electrons).

We performed our model computations by starting from a model atmosphere with equal intrinsic and return-current-heated temperatures, T_eff, i = T_eff, h. During the first tens of iterations, we applied only the flux correction in order to have the innermost part of the atmosphere converged. After that, we steadily increased T_eff, h up to its final value, and started linearly increasing the weight, w, on the energy balance error (and decreasing the maximum allowed temperature correction from flux error) when making temperature corrections. Typically, during the last hundreds of iterations, only the outer layers of the NS were heated until reaching the energy balance.

3. Results

3.1. Comparison of energy loss mechanisms

Let us first discuss the details of the particle stopping and energy loss mechanisms. We began by computing the deceleration of the return current particles as a function of the Thomson optical depth using a simplified atmosphere model (grey atmosphere) to confirm that the results are similar to those of Bauböck et al. (2019). We also studied, when the inclusion of the bremsstrahlung energy losses (neglected in all the previous works) becomes important. The comparison of electron deceleration is shown in Fig. 1, in case of a grey atmosphere, with log g = 14.0, T_eff = 4.64 MK (0.4 keV), and assuming mono-energetic return current particles with γ₀ = 2, 10, 50, or 500.

Fig. 1. Electron deceleration as a function of Thomson optical depth for different values of incoming beam Lorentz factor γ0 for grey atmosphere. Blue, orange, green, and red curves correspond to energy losses due to Coulomb collisions and wave excitation for γ0 = 2, 10, 50, and 500, respectively. The dashed lines show the same for energy losses where also bremsstrahlung is taken into account. The dotted line shows the result for γ0 = 500 when we use a return-current-heated atmosphere instead of grey atmosphere (for model parameters given in Sect. 3.1). Inclusion of bremsstrahlung losses is found to be important for particles with γ0 ≳ 100.

In addition, for γ₀ = 500 we present the result using a return-current-heated atmosphere from Sect. 3.4 (with δ = 3, γ_min = 10, T_eff, h = 5 MK, T_eff, i = 1.6 MK, and log g = 14.3). This demonstrates that dγ/dτ depends only slightly on the atmosphere model. It only depends on n_e through the plasma frequency in Eq. (13), and therefore the heated atmosphere with much more rarefied upper layers (structure discussed more later in Sects. 3 and 4) causes slightly faster deceleration. Our calculations also indicate that the bremsstrahlung losses are important for stopping electrons that have initial Lorentz factors higher than γ₀ ≈ 100.

Similarly to Bauböck et al. (2019), we also detect a sharp peak at the effective stopping depth of the particle. This is caused by the rapid increase of the energy deposition rate at relatively low particle velocities. The height of the peak is determined by our choice of the velocity where the particles are considered to have lost all of their energy (see Sect. 2.4). The equations presented in Sect. 2.2 are not valid for velocities comparable to thermal velocities of the plasma electrons because the cross-sections diverge at β = 0. In any case, the peak is not important to our calculation because its contribution to the energy loss is minor; the peak is also not resolved when interpolating the calculated energy losses to the more coarse optical depth grid, where the radiative transfer and temperature corrections were calculated. Similar selection of the grid parameters was also used in Suleimanov et al. (2012). This coarse grid is used in all the subsequent figures after Fig. 2.

Fig. 2. Comparison of average electron deceleration (see the definition in Sect. 3.1) for multi-energetic particles using grey atmosphere with different power-law indices δ and minimum cutoff energies γmin (cf. Fig. 5 in Bauböck et al. 2019). The green curve is for δ = 2 and γmin = 10, the blue curve for δ = 2 and γmin = 2, the red curve for δ = 3 and γmin = 10, and the orange curve for δ = 3 and γmin = 2; in all cases the upper limit for the energy distribution is γmax = 2 × 105. The dashed curves present the results where the bremsstrahlung losses are taken into account.

The average electron deceleration rates for multi-energetic particles (see Sect. 2.3 for more details) and for grey atmosphere are shown in Fig. 2, where $\langle \text{d}\gamma /\text{d}\tau \rangle ={\int }_{{\gamma }_{\text{min}}}^{{\gamma }_{\text{max}}}f({\gamma }_0)(\text{d}\gamma /\text{d}\tau )\text{d}{\gamma }_0$. In this case, the electrons lose their energy and stop at significantly lower depths if the bremsstrahlung energy loss is taken into account. We also note that we have a significantly higher deceleration rate (and thus energy loss) in the deeper layers than in Bauböck et al. (2019; where the energy deposition rates are cut off around τ ≈ 10). This could be an outcome of a different choice for γ_max when integrating over the energy distribution. In the electron deceleration shown in Fig. 2, we detect a cutoff at high depths only if the bremsstrahlung losses are taken into account. We also note that we made the stopping of the particle less abrupt by using an exponential function to decrease the energy loss rate when the velocity becomes smaller. This decreases the numerical noise in the results when integrating over different particle Lorentz factors, γ₀ (due to the sharp peaks seen in Fig. 1) and simultaneously avoids the need to have an extremely fine grid resolution for γ₀.

3.2. Atmospheres heated by mono-energetic particles

Let us now discuss the emergent radiation from atmospheres heated with mono-energetic particles. In this and in the following sections, we calculate atmosphere models with T_eff, h = 5.0 MK, T_eff, i = 1.6 MK, and log g = 14.3, unless otherwise stated. These values are feasible for millisecond pulsars, which are our primary targets. The effective temperature was chosen roughly correspond to the hottest observed millisecond pulsars (Guillot et al. 2019) in order to emphasise the effects of external heating (see also Sect. 3.5 for different temperatures). For simplicity, the bremsstrahlung energy losses were not taken into account in the case of mono-energetic particles, because the effect is very small for the chosen particle energies. The results of taking the effect into account are only shown in Sect. 3.3, where the particle distributions extend to higher energies.

We began by comparing the spectra and temperature structures from atmosphere models where return-current-heating dominates to that which includes only the intrinsic heat from NS (with T_eff, i = 5 MK). We also compared the results for mono-energetic particles of different initial energies. The results are shown in Fig. 3. The calculated atmospheres attain a two-layer structure: the beam of bombarding particles flows first through a hot and rarefied over-heated outer layer and then stops in a cooler and denser inner layer. The over-heated layer has such high temperature because the energy, which is released in this rarefied plasma, cannot be emitted by thermal radiation, and only Compton down-scattering can cool the plasma. The last process is effective at high temperatures only.

Fig. 3. Atmosphere structure (left) and the spectrum of the escaping radiation (right) for mono-energetic return current. Left panels: dependence of the electron number density ne, energy loss of the return current mQ+, and the temperature T on the column density. The red curves correspond to a model with Lorentz factor γ0 = 500, the blue curves are for γ0 = 50, and the black curves are for γ0 = 10. The parameters are Teff, h = 5 MK, Teff, i = 1.6 MK, and log g = 14.3. The corresponding structure and the spectrum for model without external heating (Teff, i = 5 MK, Teff, h = 0) are shown with the magenta dashed curve.

For high energy of the in-falling particles, the energy deposition happens in very deep layers and the resulting atmosphere structure resembles closely the model with deep-heated atmosphere in a radiative equilibrium. The difference in temperature at high depths is closely related to the fact that the intrinsic NS flux in the two models is very different (the effective temperature is three times smaller in heated atmosphere models) and T ∝ T_eff, iτ^1/4. Close to the surface layers, we see a temperature inversion where the temperature jumps to T ≈ (6 − 7) × 10⁷ K (see red curves in left panels of Fig. 3). The discontinuity in the electron number density is formed at the same depth.

The penetrating particles that have the lowest energies, deposit most of their energy in the upper layers, therefore the models with γ₀ = 10 produce a very hot skin, where T ≈ 10⁹ K, and the depth where the temperature inversion occurs increases (see black curves in Fig. 3) compared to the case of high-energy incoming particles. This model has largest deviation from a non-heated atmosphere.

The spectra of the escaping radiation, also presented in Fig. 3, show that the models with higher contribution of low-energy particles have their spectral peak at lower energies (especially with γ₀ = 10). They also have both a high-energy tail and increased emission at low photon energies, compared to the almost black-body like spectra in case of high-energy penetrating particles. The spectrum with γ₀ = 500 is very similar to the deep-heating model, whereas the spectrum with γ₀ = 10 deviates highly from them. The model with γ₀ = 50 shows similar discrepancies when comparing to the deep-heating model, but with a slightly smaller magnitude.

3.3. Atmospheres heated by distribution of particles

Let us now consider a more realistic return-current energy distribution. More specifically, we used a power-law distribution of bombarding particle energies, expected for realistic magnetosphere return currents (Cerutti et al. 2016). We considered different slope δ of the distribution and studied its influence on the atmosphere structure and the emitted spectra (the fiducial model parameters are shown in Table 1). The results for δ = 1, 2 and 3 are shown in Fig. 4 (assuming that γ_min = 10) and compared with corresponding quantities of a non-heated atmosphere. In addition, the effect of including the bremsstrahlung energy losses is also shown. We see that the atmosphere structure and the spectra deviate significantly from those of the non-heated atmosphere if there is contribution of low-energy particles.

Fig. 4. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles. Left panels: energy loss of the return current mQ+ and the temperature T dependence on the column density. The red curves correspond to δ = 1, the blue curves to δ = 2, and the black curves to δ = 3. Other parameters are given in Table 1. The dotted curves show the corresponding results when also the bremsstrahlung energy losses are taken into account. The dashed magenta curves correspond to the temperature structure and escaping spectrum of the atmosphere without external heating (almost completely overlapping the red solid curve).

For a steep slope with δ = 3 (corresponding average Lorentz factor ⟨γ₀⟩≈20) there are many low-energy particles that heat the surface layers resulting in the 100 keV skin and a significant high-energy photon tail above a few keV. On the other hand, the case with δ = 1 correspond to a high average Lorentz factor ⟨γ₀⟩≈10⁴ and the results differ very little from the non-heated case. In the case with δ = 2 (⟨γ₀⟩≈100), the outer layers of the atmosphere are already largely heated, but the resulting spectra shows major deviations from the non-heated case only at very low and high energies. All the cases are also well expected based on the comparison of their average Lorentz factors to the mono-energetic results of Sect. 3.2. In addition, the colour temperature again decreases (as spectral peak shifts) when having lower energy return current particles. Implications of this (and other deviations from the non-heated model) to NS parameter fitting are discussed in Sect. 4.1.

The energy loss due to bremsstrahlung radiation seems not to affect significantly the calculated spectrum and atmosphere structure, except for the temperature profile in the case with δ = 1. In that case, the inclusion of bremsstrahlung losses makes the highest-energy particles lose their energy significantly faster and at lower depths, which is enough to create a much hotter outer layer. However, the differences in the spectra are rather small, since the hot skin is not very extended and most of the photons escape from deeper layers.

We also studied how the results depend on the minimum Lorentz factor γ_min (assuming δ = 2). They are shown in Fig. 5 (with γ_min = 10, γ_min = 50, and γ_min = 200, and with no bremsstrahlung losses included). The lowest energy case (black line) is the same as the intermediate energy case (blue line) of the previous comparison. The two other cases correspond to average Lorentz factors ⟨γ₀⟩≈400 and ⟨γ₀⟩≈1400.

Fig. 5. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 for different γmin = 10 (black curves), γmin = 50 (blue), and γmin = 200 (red). Other NS atmosphere parameters are given in Table 1. The temperature structure and the spectrum without heating are shown with magenta dashed curves.

The results show again, that a larger contribution of low-energy particles (smaller γ_min) leads to hotter upper layers and to a spectrum that deviates more from the non-heated model. We also note that using the chosen distribution, the choice of the lower limit γ_min has a much larger impact on the results than the choice of the upper limit γ_max. For example, we tested for γ_min = 10, that γ_max = 10⁴ would give closely the same structure and spectrum as γ_max = 2 × 10⁵ that was used throughout the paper. This could also be seen from the results shown in Fig. 4, where inclusion of bremsstrahlung energy losses acted in the same way as having significantly smaller γ_max.

3.4. Atmospheres with different surface gravities

We have also studied the dependency of the results on the surface gravity in the case of distribution of particles with δ = 2 and γ_min = 10. The results for log g = 13.7, 14.0, 14.3, and 14.6 are shown in Fig. 6 (with no bremsstrahlung losses included). We see that the surface gravity has a much smaller effect on the structure and the spectra of the atmosphere compared to the effects of the energy distribution parameters δ and γ_min of the bombarding particles, discussed in the previous section. However, the effect is still significant, and similarly as growing δ or decreasing γ_min, decreasing the value of g leads to higher amount of hard energy photons. In this case, the return current energy loss pattern does not change, as seen from the upper left panel of Fig. 6, meaning that the changes in the spectrum and temperature profile are caused by other means (see discussion in Sect. 4). It seems that the lowest surface gravity allows largest inward expansion of the hot skin, resulting in the hardest spectrum.

Fig. 6. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 and γmin = 10 for different NS surface gravities: log g = 13.7 (green curves), 14.0 (black), 14.3 (blue), and 14.6 (red). Other NS atmosphere parameters are given in Table 1.

3.5. Atmospheres with different effective temperatures

Additionally, we have considered the dependency of the results on the effective temperature T_eff, h (keeping the ratio T_eff, h/T_eff, i ≈ 3 same as in the other computations). The fiducial 5 MK temperature (used in the models of other sections) is relatively high compared to the temperature estimates of recently observed RMPs (Guillot et al. 2019). Therefore, we have calculated the models (again with δ = 2 and γ_min = 10) for T_eff, h = 2, 3, 4, and 5 MK. The results are shown in Fig. 7 (solid lines) and compared to the non-heated models of the same effective temperature (dashed lines). The results show that lower temperature is related to less-pronounced hard tails in the spectra, because for colder atmospheres the temperature inversion happens at lower depths. This is caused by the cooling rate being more extensively dominated by free-free radiation instead of Compton scattering (see discussion in Sect. 4.). However, the increased emission due to heated upper layers can still be detected at high energies even for T = 2 MK.

Fig. 7. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 and γmin = 10 for different NS effective temperatures: Teff, h = 2 (green solid curves), 3 (black), 4 (blue), and 5 MK (red). Other NS atmosphere parameters are given in Table 1. The dashed curves correspond to the temperature structure and escaping spectrum of the atmosphere without external heating (for Teff, i = 2, 3, 4, and 5 MK).

3.6. Beaming patterns

Let us now discuss the angular distribution of the emergent radiation. The emission pattern of radiation escaping from the atmosphere of log g = 14.3 and T_eff, h = 5 MK for a power-law distribution of in-falling particles with γ_min = 10 and three slopes δ = 1, 2 and 3 is shown in both in polar (left) and Cartesian (right) coordinates in Fig. 8. The specific intensities at different cosines of zenith angles (denoted as μ) were obtained directly using the radiative transfer equation, that was solved finally in eleven angles (instead of three angles used during the major part of temperature iterations). The intensities for the intermediate angles were obtained using a third order spline interpolation from the calculated points.

Fig. 8. Emission pattern of the specific intensity in polar (left) and Cartesian (right) coordinates for the NS atmosphere models presented in Fig. 4. The dashed curves show the patterns for non-heated atmosphere model, and the solid curves are for a power-law distribution of return current particles. Top, middle and bottom panels correspond to δ = 1, 2 and 3, respectively. The black, blue, green, orange, and red colours correspond to intensities at 0.1, 0.3, 1, 3, and 10 keV, respectively.

We see that in the case of low-energy particle heating (i.e. δ = 3) the beaming of the radiation at low and high energies is clearly different, unlike in the case of non-heated atmosphere. The energy, where the limb darkening at low energies changes to the limb brightening at higher energies, depends on the return-current electron energy distribution. In case of δ = 1, the angular distribution of radiation resembles very closely the distribution of a non-heated atmosphere, and no limb brightening is seen. When δ = 2, we observe a rapid transition from limb darkening to limb brightening at the highest zenith angles and highest energies (10 keV). Even in the thermal part of the spectrum, we see large deviations of the angular distribution from that predicted by standard non-heated model. With energy distribution having even higher small-γ contribution (i.e. δ = 3), the energies above 3 keV show even stronger brightening, for all angles, towards the surface tangent. The deviations between heated and non-heated models become slightly larger also at the smaller energies, when a significant fraction of the energy of the beam is dissipated high up in the atmosphere (δ ≥ 2).

The emergent spectra for different zenith angles are shown in Fig. 9. From this we see that the high-energy tail of the spectrum is larger for high zenith angles, as expected. On the other hand, the peak of the spectrum is higher for small angles. At very low energies again the radiation intensity becomes slightly larger at high zenith angles.

Fig. 9. Emergent specific intensity spectra for power-law distribution of the particle beam for the NS atmosphere model shown in Fig. 4 in case of δ = 2 and γmin = 10. The results are shown for three emission angles (out of the eleven that were computed): μ = 0.99 (red solid curve), μ = 0.50 (blue dashed), and μ = 0.01 (black dotted).

We also found similar discrepancies between the beaming patterns of heated and non-heated models in case of a colder atmosphere (T_eff, h = 2 MK or without heating T_eff, i = 2 MK). The results are shown in Fig. 10. We see again rapid transition from limb darkening to limb brightening at 10 keV, although the energies should not be directly compared to those computed with different temperature (since the thermal peak is shifted). Nevertheless, the differences in the angular distributions remain at 10–30% level even at smaller energies. This demonstrates that the effects of return-current-heating effects should be relevant also for RMPs with relatively low temperatures studied recently (Miller et al. 2019; Riley et al. 2019).

Fig. 10. Emission pattern of the specific intensity for the NS atmosphere models in a low temperature case with Teff = 2 MK (with δ = 2) presented in Fig. 7. The dashed curves show the patterns for non-heated atmosphere model, and the solid curves are for a power-law distribution of return current particles. The black, blue, green, orange, and red colours correspond to intensities at 0.1, 0.3, 1, 3, and 10 keV, respectively.

The differences between the beaming patterns can also be coarsely characterised by inspecting the values of intensity emitted towards tangential direction μ = 0. A summary of all the computed models (excluding those where bremsstrahlung losses and mono-energetic particles were considered) and the values obtained for the normalised intensity at μ = 0 (defined as a = I(μ = 0)/I(μ = 1)) for four different energies (0.3, 1, 3, and 10 times the colour temperature kT_c of each model) are shown in Table 2. The colour temperatures (and the colour correction factors f_c) were obtained by fitting a diluted blackbody function to the emergent spectrum. The fitting procedure was similar to the first method described in Suleimanov et al. (2011), although the energy band was adjusted to (0.1 − 6) × (1 + z) keV (for RMPs detected by NICER), where z is the gravitational redshift, obtained from log g by assuming a NS mass of 1.4 M_⊙.

We can see that the angular dependency of the intensity is different from non-heated models also when making the comparison at the same energies relative to the spectral peak (see Table 2). The location of the peak varies between the models, so that more over-heated upper layers, or higher effective temperature, produce smaller f_c. From Table 2 we also see that the beaming at large angles is most insensitive (although not entirely) to the used model at close to the peak of the spectrum (a₃ at 3kT_c). At these energies the angular dependency can be approximated by a linear function as I(μ)/I(1) = a + (1 − a)μ. For energies above and below the spectral maximum, the deviations in the beaming become larger, and the approximately linear dependency breaks down.

4. Discussion

4.1. Difference from deep-heating models and implications for NICER results

The results presented in Sect. 3 demonstrate that the spectra and beaming patterns for the return-current-heated atmosphere models can significantly differ from those assuming that heat is released at the bottom of the atmosphere (non-heated models). As we have shown, the results depend on the unknown energy distribution of the return-current particles. The difference is very large for the smallest energies of the bombarding particles, but insignificant for the highest energies that we have considered. Even when the spectral differences are small (e.g. for the cases with γ₀ = 500, or δ = 1, or γ_min = 200), the two-layer structure with an overheated upper layers can still be produced. The heated part of the atmosphere consists of a hot optically thin skin and a cooler optically thick inner part. The increased temperature, due to the return current, forces the outer layers to expand, as the pressure does not change, and creates a discontinuity in the particle concentration at the depth where the temperature inversion occurs. This allows for a possibility to represent the RMP atmosphere with a simplified two-layer model consisting of a hot Comptonizing slab above a usual deep-heated atmosphere.

We note that the depth of the temperature inversion is likely related to the change in the dominant opacity mechanism. The cooling rates due to free-free radiation and due to Compton scattering are almost equal always at the point of the temperature inversion. In case of a higher surface gravity, the free-free opacity dominates at larger range of depths in the atmosphere (because of higher gas pressure and electron number densities for the same temperature) and therefore the temperature inversion occurs also at a lower depth and the spectrum is closer to the blackbody (as seen in Fig. 6). In case of a fixed surface gravity, the amount of energy deposited in the upper layers determines the inversion depth, so that the inversion occurs at higher depth and outer layers become hotter when the amount of the deposited energy in the upper layers is increasing. Nevertheless, most of the observed radiation still escapes from the layers deeper than the inversion depth in all of our models. For models with less energetic bombarding particles, the layers producing thermal radiation are less efficiently heated, and that can explain why the peak of the spectrum is shifted towards the lower energies (see Figs. 3–5). As a consequence, this leads to a substantial bias in the estimate of the atmosphere effective temperature (and also in the size of the emitting spot) when constraining NS parameters using phase-resolved spectra, if incorrect assumptions about the atmosphere model are made.

As we have shown above, the radiation spectrum and the angular distribution of escaping radiation may deviate significantly from those predicted by standard hydrogen atmosphere models without external heating. What is important is that the angular distribution even of the thermal part of the spectrum (see blue and green curves in the middle panel of Fig. 8 and in Fig. 10), when the high-energy tail is completely negligible (see Fig. 7), shows 10–50% difference from the corresponding angular distributions of non-heated models. Obviously, a different emission pattern of the hotspots on the surface of RMPs would produce a significantly different pulse profile (see e.g. Poutanen & Beloborodov 2006), with characteristic difference at the tens of per cent level. This is orders of magnitude larger than the accuracy of 0.1% the NICER team aims at in their pulse profile modelling (Bogdanov et al. 2019). Because in both NICER papers on NS parameters constraints from PSR J0030+0451 (Miller et al. 2019; Riley et al. 2019) the spectral shape and the angular distribution were taken from standard non-heated models, the obtained constraints on the NS radius and the structure of magnetosphere may be heavily biased. Detailed analysis and improved constraints on the NS parameters and magnetospheric geometry obtained from the pulse profile modelling of RMPs using return-current-heated atmosphere models will be considered in a forthcoming publication.

Conversely, investigation of the thermal spectra of RMPs opens a possibility to study the energy distribution of the relativistic particles by using the atmosphere model to constrain the particle distribution parameters. This would allow us to also indirectly probe the exact physics of pulsar magnetospheres via the return-current structure. Such studies can provide an interesting new observational method for actually validating, for example, the different proposed pulsar magnetosphere gap models or pair-production sites in the magnetosphere since both of these physical details alter the incoming return current properties (see e.g. Chen & Beloborodov 2014; Cerutti & Beloborodov 2017; Beskin 2018; Philippov & Spitkovsky 2018).

4.2. Comparison with previous works

Compared to the results of Bauböck et al. (2019), who also studied the return-current-heated atmospheres, our calculations show similar features in the spectra and beaming patterns, although having a significantly higher temperature of the outer atmosphere layers, using the same parameters for the energy distribution of the return current. Thus, also our spectrum deviates more from a blackbody, typical for a non-heated atmosphere. In our case, spectra are harder and the resulting beaming pattern slightly more limb brightened at high energies for low return-current electron energies.

The main difference of our atmosphere model compared to that of Bauböck et al. (2019) is in the modelling of the radiative transfer. We account exactly for the effects of electron scattering and perform calculations at a grid of photon energies, while Bauböck et al. (2019) consider grey atmosphere. This is likely the reason why the atmosphere structure and the resulting spectra are significantly different. As a part of accurate radiation transfer modelling, we also use energy-dependent opacities instead of averaged ones. We checked that in the overheated layers the Planck-weighted opacity differs significantly from the opacity that is weighted using correct intensities.

Qualitatively, our results are similar to those of Suleimanov et al. (2018), who instead of RMPs, studied the accretion-heated atmospheres of NS. The angular distribution of the emergent radiation is comparable to ours, as we see the limb darkening at lowest energies switching to the limb brightening at higher energies. In addition, they have the pronounced two-layer structure of the atmosphere, and emergent spectra showing a similar excess at lower photon energies in comparison with the blackbody. This is a consequence of the free-free opacity being high at low energies causing the observed low energy photons to be produced in the overheated upper layers. This ‘reverse photosphere effect’ was already discussed by Deufel et al. (2001), who also obtained similar structures and spectra using a pure hydrogen atmosphere heated by accreted protons. It could also be possible to detect the low energy excess as an increased emission at optical wavelengths. However, the observed optical excess, usually associated with highly magnetised neutron stars (van Kerkwijk & Kaplan 2007) may be explained by magnetic atmosphere models instead (Ho et al. 2007; González-Caniulef et al. 2019).

4.3. Connection to pulsar physics

In pulsar physics, instead of the effective temperature T_eff, h, the parameter of interest is the number density of precipitating particles in terms of the Goldreich-Julian number density defined as (Goldreich & Julian 1969)

$$\begin{array}{c}\hfill n_{\mathrm{GJ}}\approx \frac{B}{\mathit{ePc}},\end{array}$$(23)

where B is the magnetic field strength, P is the period of the pulsar, and e is the elemental charge. It defines a characteristic density in the pulsar magnetosphere that is needed to screen the longitudinal electric field near the neutron star surface (Beskin 2018). The pulsar magnetosphere has a total number density n_tot = ℳn_GJ, where ℳ is the pair multiplicity parameter.

The number density of the inward penetrating particles suggested by our atmosphere models (see Eq. (19)), on the other hand, is given as

$$\begin{array}{c}\hfill n_{\mathrm{in}}\approx \frac{{\dot{N}}_{\mathrm{e}}}{c}=\frac{{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff},\mathrm{h}}^4}{(⟨{\gamma }_0⟩-1)m_{\mathrm{e}}{c}^3}=\zeta \mathcal{M}{n}_{\mathrm{GJ}},\end{array}$$(24)

where ζ = n_in/n_tot is the ratio of in-going to the total number of pairs. Both ζ and ℳ are currently unknown and depend on the exact details of pulsar magnetosphere structure and pair cascade physics.

Multiplicities corresponding to our models are shown in Table 2, when assuming typical values for millisecond pulsars (B = 10⁹ G and P = 1 ms). For simplicity, we have also taken ζ = 1. The multiplicity depends on our model parameters T_eff, h, δ and γ_min (the two latter through ⟨γ₀⟩) as seen from Eq. (24). The presented values range from ℳ ∼ 10³ to 10⁶ which are reasonable given the theoretical uncertainties in pulsar magnetosphere physics (Timokhin & Arons 2013).

We note that the already large uncertainties of ζ and M are increased when considering millisecond pulsars. The details of the pair-production mechanism in millisecond pulsars must differ from slower rotating pulsars that can enable efficient pair production via a strong magnetic field; in millisecond pulsars the process could be seeded, for example, by high-energy inverse Compton photons instead. However, to have at least a modest view on the return current physics and particle distribution, we briefly consider next our models in the context of recent developments in simulations of rotation-powered pulsars.

One-dimensional particle-in-cell (PIC) simulations of the polar cap pair cascades, created near the magnetic poles, have shown that a large fraction of the energy may be carried by particles having very high Lorentz factor of γ ∼ 10⁷ (Timokhin 2010; Timokhin & Arons 2013). However, in these studies also a significant population of low-energy particles flow towards the NS surface heating the upper layers of the atmosphere. Global multi-dimensional PIC simulations of pulsar magnetospheres has been performed as well, but using scaled energies of the return current particles (Chen & Beloborodov 2014, 2017; Philippov & Spitkovsky 2018). Therefore, the typical Lorentz factor of secondary pairs cannot be directly inferred. However, the particle energy distribution depends on the details of the model and the assumed type of the pulsar. These uncertainties translate to a large uncertainty in the ζ parameter, which could be, for example, around 0.1 for strong-field pulsars (Timokhin & Arons 2013).

One-dimensional PIC simulations have also indicated that pair creation must be time-dependent and exhibit a quasi-periodic behaviour in accelerators where the pair formation is not suppressed (Beloborodov & Thompson 2007; Beloborodov 2008; Timokhin & Arons 2013). The non-stationary nature of the process means that the physically relevant parameters to consider are, in reality, time-averaged quantities $\overline{\zeta }$ and $\overline{M}$, if the characteristic relaxation time scale of the atmosphere is longer than time scale of the fluctuations. The latter may be from fractions of microseconds to tens of microseconds depending on the altitude where the cascade starts (Timokhin & Arons 2013). The thermal relaxation time scale of the atmosphere can also be about the same order of magnitude, depending however on the model. Therefore, a demand for time-dependent atmospheres models (which are not considered here) cannot be ruled out. However, we note that if the amplitude of the fluctuations is small enough, as it may be in multi-dimensional simulations, the atmosphere is not affected.

4.4. Caveats

There are a few uncertainties in our calculated model atmospheres. One is our assumption that all the dissipated energy (from the return current particles) translates to the outward directed flux at every depth point, as can be seen in our definition of the theoretical flux (the nominator in Eq. (21)). The same assumption is implicitly used also in our energy balance correction, as the sign of Q⁺ has been chosen to correspond to a negative flux derivative implying an outward flux. The assumption might not be physical, and inward flux could occur, for instance, just below the stopping depth of the penetrating particles, if the initial heat of the NS is significantly smaller than heat due to the stopped particles. However, considering the internal heating of the NS, caused by the return current, is not in the scope of this paper. Also, if all the incoming energy is lost at optical depths smaller than 1, as seems to be in our models, most of the energy is expected to be radiated outwards anyway.

In our analysis we neglected the effects of thermal conduction, affecting the heating and cooling rate of the matter. However, based on the results of Suleimanov et al. (2018) this effect is expected to be rather small. Also, the electron-positron pair production effects were neglected. For some of the model parameters, where we obtain outer layers as hot as T = 10⁹ K, we could have an outflow of pairs (Zane et al. 1998). Detailed study of this effect is left for a future work. Furthermore, as we mentioned above, the bombarding particles are actually positrons, and therefore they are annihilated with the background electrons once their energy is low enough, producing an annihilation line at 511 keV. Half of those photons may reach the observer directly and another half are Compton back-scattered, making an extended tail containing more than 1/γ₀ fraction of the total energy.

Another simplification in our model is related to the bremsstrahlung radiation as a source of heating. We assume that all the bremsstrahlung radiation of the stopping particles is converted to heat of the surrounding gas at the same optical depth. This assumption should be valid in the deepest layers (where generated high-energy photons rapidly lose their energy due to Compton scattering) but some deviations in the optically thin upper layers could be expected (photons would rather heat the layers below them). However, as seen from those results where bremsstrahlung effects were included, the resulting spectra and temperature structure do not have major changes because of this energy loss mechanism. In addition, the bremsstrahlung radiation should still be taken into account as an additional source of photons in the radiation transfer equation. The inward moving high-energy bremsstrahlung photons produced at relatively low depths are Compton back-scattered causing a bump at energies above the thermal peak. This effect is most important when bremsstrahlung energy losses are highest and the spectrum would otherwise look thermal, like in non-heated models.

A major uncertainty lies in the energy spectrum of the magnetospheric particles, especially in the low-γ regime as discussed also by Bauböck et al. (2019). The spectra for the particle beam could deviate from the power-law. There could be, for example, different contribution of low-γ, either because of different shape of the distribution, or because of using different lower γ_min and the upper limits γ_max of the distribution. However, this uncertainty was considered in this study, as we examined the results using different values for the distribution parameters.

5. Conclusions

We have presented a model for return-current-heated atmospheres of RMPs, using the full radiative transfer calculation and accounting exactly for Compton scattering. We have assumed that the magnetospheric return current that heats the polar caps of NSs consists of pair plasma and considered different energy distributions. We compared various energy loss mechanisms of the return current, concluding that the effect of the bremsstrahlung energy loss is small. Finally, we computed the temperature structures, emergent spectra and angular distribution of radiation for different model parameters.

We found that the usual deep-heating approximation deviates significantly from the return-current-heated model if the bombarding particles have a significant contribution at Lorentz factors smaller than about 100. The models including high contribution of low-energy particles resulted in a very hot skin, reaching the temperatures around 10⁹ K, that led to a strong high-energy tail in the spectrum and a switch of the beaming pattern from limb darkening at lower energies to a strong limb brightening at highest energies. In the opposite case of high contribution of high-energy particles, a rather hot skin could still be produced, although at a lower temperature, but the spectrum and the beaming pattern were very similar to those of the non-heated atmospheres. We also found that the surface gravity of the NS affects the strength of the high-energy tail of the spectrum.

Large change in the emergent radiation spectrum and its beaming properties can have a significant impact on the inferred NS and magnetosphere parameters such as those obtained by the NICER team from RMPs using pulse profile modelling (Miller et al. 2019; Riley et al. 2019). The deviations of the angular distribution of specific intensity from that predicted by standard deep-heating models may exceed tens of per cent even in the thermal part of the spectrum, when the high-energy tail is nearly invisible. These deviations produce potentially orders of magnitude larger changes to the pulse profiles than the systematic uncertainty assumed in NICER studies (Bogdanov et al. 2019). The results of our work may be used to improve NS radius constraints using these data. Finally, because the exact properties of the emergent radiation are strongly dependent on the magnetospheric return-current energy distribution, the presented new atmosphere models can, in the future, be used to probe the still unknown structure and pair-production physics of the pulsar magnetospheres.

Acknowledgments

This research was supported by the University of Turku Graduate School in Physical and Chemical Sciences (TS), by the Ministry of Science and Higher Education of the Russian Federation grant 14.W03.31.0021 (JP, VFS), the Deutsche Forschungsgemeinschaft (DFG) grant WE 1312/51-1 (VFS), the German Academic Exchange Service (DAAD) travel grants 57405000 and 57525212 (VFS), and the Academy of Finland grants 317552, 331951 and 333112 (TS, JP). TS thanks Cole Miller and Ilia Kosenkov for useful discussions. The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Grid and Cloud Infrastructure project are acknowledged.

References

All Tables

All Figures

Fig. 1. Electron deceleration as a function of Thomson optical depth for different values of incoming beam Lorentz factor γ0 for grey atmosphere. Blue, orange, green, and red curves correspond to energy losses due to Coulomb collisions and wave excitation for γ0 = 2, 10, 50, and 500, respectively. The dashed lines show the same for energy losses where also bremsstrahlung is taken into account. The dotted line shows the result for γ0 = 500 when we use a return-current-heated atmosphere instead of grey atmosphere (for model parameters given in Sect. 3.1). Inclusion of bremsstrahlung losses is found to be important for particles with γ0 ≳ 100.

In the text

Fig. 2. Comparison of average electron deceleration (see the definition in Sect. 3.1) for multi-energetic particles using grey atmosphere with different power-law indices δ and minimum cutoff energies γmin (cf. Fig. 5 in Bauböck et al. 2019). The green curve is for δ = 2 and γmin = 10, the blue curve for δ = 2 and γmin = 2, the red curve for δ = 3 and γmin = 10, and the orange curve for δ = 3 and γmin = 2; in all cases the upper limit for the energy distribution is γmax = 2 × 105. The dashed curves present the results where the bremsstrahlung losses are taken into account.

In the text

Fig. 3. Atmosphere structure (left) and the spectrum of the escaping radiation (right) for mono-energetic return current. Left panels: dependence of the electron number density ne, energy loss of the return current mQ+, and the temperature T on the column density. The red curves correspond to a model with Lorentz factor γ0 = 500, the blue curves are for γ0 = 50, and the black curves are for γ0 = 10. The parameters are Teff, h = 5 MK, Teff, i = 1.6 MK, and log g = 14.3. The corresponding structure and the spectrum for model without external heating (Teff, i = 5 MK, Teff, h = 0) are shown with the magenta dashed curve.

In the text

Fig. 4. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles. Left panels: energy loss of the return current mQ+ and the temperature T dependence on the column density. The red curves correspond to δ = 1, the blue curves to δ = 2, and the black curves to δ = 3. Other parameters are given in Table 1. The dotted curves show the corresponding results when also the bremsstrahlung energy losses are taken into account. The dashed magenta curves correspond to the temperature structure and escaping spectrum of the atmosphere without external heating (almost completely overlapping the red solid curve).

In the text

	Fig. 5. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 for different γmin = 10 (black curves), γmin = 50 (blue), and γmin = 200 (red). Other NS atmosphere parameters are given in Table 1. The temperature structure and the spectrum without heating are shown with magenta dashed curves.
In the text

	Fig. 6. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 and γmin = 10 for different NS surface gravities: log g = 13.7 (green curves), 14.0 (black), 14.3 (blue), and 14.6 (red). Other NS atmosphere parameters are given in Table 1.
In the text

Fig. 7. Atmosphere structure (left) and the spectrum of escaping radiation (right) for a power-law distribution of bombarding particles with δ = 2 and γmin = 10 for different NS effective temperatures: Teff, h = 2 (green solid curves), 3 (black), 4 (blue), and 5 MK (red). Other NS atmosphere parameters are given in Table 1. The dashed curves correspond to the temperature structure and escaping spectrum of the atmosphere without external heating (for Teff, i = 2, 3, 4, and 5 MK).

In the text

Fig. 8. Emission pattern of the specific intensity in polar (left) and Cartesian (right) coordinates for the NS atmosphere models presented in Fig. 4. The dashed curves show the patterns for non-heated atmosphere model, and the solid curves are for a power-law distribution of return current particles. Top, middle and bottom panels correspond to δ = 1, 2 and 3, respectively. The black, blue, green, orange, and red colours correspond to intensities at 0.1, 0.3, 1, 3, and 10 keV, respectively.

In the text

	Fig. 9. Emergent specific intensity spectra for power-law distribution of the particle beam for the NS atmosphere model shown in Fig. 4 in case of δ = 2 and γmin = 10. The results are shown for three emission angles (out of the eleven that were computed): μ = 0.99 (red solid curve), μ = 0.50 (blue dashed), and μ = 0.01 (black dotted).
In the text

Fig. 10. Emission pattern of the specific intensity for the NS atmosphere models in a low temperature case with Teff = 2 MK (with δ = 2) presented in Fig. 7. The dashed curves show the patterns for non-heated atmosphere model, and the solid curves are for a power-law distribution of return current particles. The black, blue, green, orange, and red colours correspond to intensities at 0.1, 0.3, 1, 3, and 10 keV, respectively.

In the text