A&A 377, 955-963 (2001)
DOI: 10.1051/0004-6361:20011129

X-ray spectra from protons illuminating a neutron star

B. Deufel - C. P. Dullemond - H. C. Spruit

Max-Planck-Institut für Astrophysik, Karl-Schwarzschildstr. 1, 85740 Garching, Germany

Received 21 June 2001 / Accepted 8 August 2001

Abstract
We consider the interaction of a slowly rotating unmagnetized neutron star with a hot (ion supported, ADAF) accretion flow. The virialized protons of the ADAF penetrate into the neutron star atmosphere, heating a surface layer. Detailed calculations are presented of the equilibrium between heating by the protons, electron thermal conduction, bremsstrahlung and multiple Compton scattering in this layer. Its temperature is of the order 40-70 keV. Its optical depth increases with the incident proton energy flux, and is of the order unity for accretion at 10-2-10-1 of the Eddington rate. At these rates, the X-ray spectrum produced by the layer has a hard tail extending to 100 keV, and is similar to the observed spectra of accreting neutron stars in their hard states. The steep gradient at the base of the heated layer gives rise to an excess of photons at the soft end of the spectrum (compared to a blackbody) through an "inverse photosphere effect''. The differences with respect to previous studies of similar problems are discussed, they are due mostly to a more accurate treatment of the proton penetration process and the vertical structure of the heated layer.

Key words: accretion, accretion disks - radiative transfer - stars: neutron - X-rays: stars


   
1 Introduction

The spherical accretion of matter onto the surface of an unmagnetized neutron star (NS) has attracted much attention since Zel'dovich & Shakura (1969) [henceforce ZS69] first addressed this problem more than thirty years ago. As a model for the then newly discovered X-ray stars, these authors considered a neutron star whose surface was heated by radially infalling gas, and modeled this gas as consisting of freely falling ions. A hot X-ray emitting layer is formed, the temperature of which depends on the penetration depth of the protons into the atmosphere of the star. The model could reproduce the observed hard X-rays, but was eclipsed by accretion disk models once it was realized that accretion is not radial because of the angular momentum constraint. The disk model, however, only explains the sub-keV part of NS spectra, somewhat in contradiction with observed spectra. In addition to the optically thick disk, an additional source of hot gas thus had to be found to produce the hard X-rays. This was proposed in the form of the so-called two-temperature accretion flows (Shapiro et al. 1976; Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994, 1995). In this form of accretion, now called ADAF, the flow is geometrically thick, supported by virialized protons, while the electrons stay at a lower temperature around 100 keV due to their high radiative energy losses and the low rate at which Coulomb interactions transfer accretion energy from the ions to the electrons. While it is not entirely clear at the moment how such a flow would originate from the cool disk present at larger distances from the star's surface, the observed hard X-ray component is often taken as an indication of its existence. If accretion near the star takes indeed place in the form of such a two-temperature plasma, the star's surface is exposed to protons with energies around 50 MeV, and the physics is very much like the model proposed by ZS69. It is therefore interesting to pick up this line of theory and revisit proton-injection models in context of NS-ADAF accretion.

In ZS69 the importance of the penetration depth of the accreted protons for the outcoming spectrum was recognized. The penetration length in terms of the amount of material required to stop the protons was estimated and used as a model parameter. The injected energy of the infalling protons per unit time and mass in the neutron star atmosphere was then uniformly distributed in this stopping layer. A second parameter of their model was the accretion rate (or luminosity). Their computed spectra essentially show a blackbody spectrum with a high-energy tail due to Comptonization in the heated atmosphere.

A more detailed numerical approach was presented by Alme & Wilson (1973). They introduced more physics into their model and also realized that the proton deceleration depends on the atmospheric temperature. They solved the atmospheric structure and the beam deceleration simultaneously in a time dependent evolution of the model. But they did not account for the strong dependence of proton deceleration on the local proton velocity and used an estimate only for the amount of material needed to stop the infalling protons. Their resulting spectra again show a blackbody plus an additional high energy tail.

Bildsten et al. (1992) were interested in the fate of the accreted CNO elements in the neutron star atmosphere. They assumed an isothermal atmosphere at a temperature defined by the accretion rate and the proton stopping depth. Again a similar estimate for the penetration depth of the ions was used.

The model of spherical accretion onto a neutron star was revisited by Turolla et al. (1994). They confirmed the existence of the formerly known "cold'' solutions (as in Alme & Wilson 1973). In addition to this, they also found "hot'' solutions for the same luminosities. But their work was limited by the fact that the spectrum was just described in terms of a mean photon energy.

Zane et al. (1998) confirmed the existence of the hot solutions with a much more accurate treatment of the radiative processes including pair processes. The emergent spectra of the hot solution should be characterized by hard spectra peaking at $\approx$100 keV. Both models use the stopping depth as a free input parameter. The effect of thermal conduction was not included in these calculations.

In continuation of this series of papers Zampieri et al. (1995) presented solutions for neutron stars accreting at low rates, i.e. with luminosities $L_{\rm NS}=10^{-7}$- $ 10^{-3} L_{\rm Edd}$. They showed that the deviations from a blackbody spectrum with the effective temperature of the neutron star increase with decreasing accretion luminosities.

Proton illumination is not only restricted to the accretion of matter onto the surface of a neutron star. Spruit (1997), Spruit & Haardt (2000) and Deufel & Spruit (2000) have shown that this process might also be of importance for accretion disks embedded in a hot corona. With an improved version of the code used in Deufel & Spruit (2000) we reconsider the accretion of matter onto the surface of a neutron star (i.e. ion illumination of a neutron star surface). The code now includes a better treatment of the radiative processes by including thermal emission due to bremsstrahlung. We also allow for energy redistribution due to electron thermal conduction within the atmosphere. The interaction between the accreted protons and the neutron star atmosphere and the radiative processes are computed time dependently in a one-dimensional, plane-parallel approximation. The density distribution through the atmosphere is found from hydrostatic equilibrium including the pressure from the penetrating ions.

In the following we give a complete description of our model. First the heating of the electrons via Coulomb interactions is presented. We show that the use of a stopping depth treated as a free parameter is not sufficiently accurate for the calculation of the proton stopping. The radiative processes are treated by solving the radiative transfer equation instead of the Monte Carlo method used before. This has advantages in terms of computation speed and the ability to treat layers of arbitrary absorption optical depth.

   
2 The neutron star illumination model

2.1 Proton stopping by Coulomb collisions

We consider protons in the vicinity of a neutron star of mass $M_{\rm *}$ and radius $R_{\rm *}$. In ADAF accretion the typical energy of such protons at the surface of a neutron star is the virial temperature,

 
$\displaystyle E_{\rm p} = k\,T_{\rm vir} = \frac{G M_{\rm *} m_{\rm p}}{3 R_{\r... ...frac{M_{\rm *}}{M_{\odot}}\right ) \left( \frac{10~{\rm km}}{R_{\rm *}}\right),$     (1)

where G is the gravitational constant, $m_{\rm p}$ the proton mass and k the Boltzmann constant.

If protons with such energies encounter a much a cooler and dense medium, they will be stopped very efficiently by Coulomb interactions with the cold electrons inside their Debye sphere. The stopping of a fast particle in an ionized plasma was quantitatively discussed by Spitzer (1962). In a time ${\rm d} t$, such a proton with energy $E_{\rm p}$ will loose an amount of energy ${\rm d} E$,

 \begin{displaymath} \frac{{\rm d} E}{{\rm d} t}= -\varepsilon(E_{\rm p})\, n_{\rm e}\, v_{\rm p}, \end{displaymath} (2)

where $n_{\rm e}$ is the electron density of the neutron star atmosphere and $v_{\rm p}^2 \approx 3 k T / m_{\rm p}$ is the proton velocity.  $\varepsilon(E_{\rm p})$ is called the stopping power of a plasma. Ryter et al. (1970) show that it is given through Spitzers's formalism by

 \begin{displaymath}\varepsilon(E)=\frac{4 \pi e^4}{m_{\rm e} v_{\rm p}^2}\,\ln \Lambda\,[\psi(x)-x\psi'(x)]. \end{displaymath} (3)

Here $\ln \Lambda=\ln[(3/2e^3)(k^3 T_{\rm e}^3/\pi n_{\rm e})^{1/2}]$ is the Coulomb logarithm, $\sigma_{\rm T}$ is the Thomson cross-section, $\psi(x)$ and $\psi'(x)$ are the error function and its derivative, and $x=(m_{\rm e} v_{\rm p}^2 / 2 k T_{\rm e})^{1/2} $ is the ratio of the incident proton velocity to the thermal velocity of the electrons. This formula holds for non relativistic electron and proton temperatures.

Now we can express the energy loss of a proton moving at an angle $\theta$ with respect to the vertical coordinate z, per unit of vertical Thomson depth, ${\rm d}\tau = n_{\rm e} \sigma_{\rm T} {\rm d} z$, by

 \begin{displaymath} \frac{{\rm d} E}{{\rm d}\tau_{\rm T}}= -\frac{1}{\sigma_{\r... ...rm e} v_{\rm p}^2\cos\theta}\,\ln\Lambda\,[\psi(x)-x\psi'(x)]. \end{displaymath} (4)

In a time-dependent evolution of the model we follow the protons numerically through the neutron star atmosphere with temperature $T(\tau)$ and density $n_{\rm e}(\tau)$. We record the loss of kinetic energy of the protons as a function of optical depth. This yields the local time dependent heating rate $\Lambda_{\rm p}(\tau)$. In this way the stopping depth of the protons is obtained self consistently.

As mentioned in Sect. 1 the penetration depth significantly influences the outcoming spectrum of the proton illumination model. Usually a global estimate is used to calculate the stopping depth of the protons in an atmosphere. The estimate for the stopping length $\lambda_{\rm s}$ is obtained following Spitzer (1962),


 \begin{displaymath} \lambda_{\rm s} = v_{\rm p} \,t_{\rm s}\;, \end{displaymath} (5)

where $t_{\rm s}$ is the slowing down time scale given by Spitzer (1962), Eqs. (5)-(28). If we express the stopping length in terms of Thomson optical depth, we get


 \begin{displaymath} \tau_{\rm s}=\frac{m_{\rm p}}{m_{\rm e}}\frac{2}{3 \ln \Lambda} \frac{(v_z/c)^4}{\psi - x\psi'}, \end{displaymath} (6)

where $v_{\rm z}$ is the vertical component of the proton thermal velocity and c is the speed of light. This expression is inaccurate because of several factors. Most importantly the proton velocity enters this formula as the fourth power. As the proton velocity decreases, this factor drops rather fast and shortens the penetration depth considerably. Counteracting is the function $\psi-x\psi'$. For x>2 this function is unity. But for 0<x<2 it is small, increasing the stopping length but also reducing the energy loss [cf. Eq. (4)]. Even the Coulomb logarithm should not be taken as constant (as it is often done). In the wide range of electron temperatures and densities of interest (as solutions from e.g. 1994 suggest) the Coulomb logarithm can easily vary by a factor of two.
 

 
Table 1: Estimated stopping depth (est) and calculated stopping depth (num) [according to Eq. (4)] in units of Thomson optical depth for a isothermal atmosphere with temperature $T_{\rm e}$ and uniform electron density $n_{\rm e} = 10^{22}$ cm-3 for various proton energies.
  Temperature $T_{\rm e}$[keV]
$T_{\rm p}\, [{\rm MeV}]$ 1 10 100
  num est num est num est
23 0.2 0.8 0.4 0.8 6.8 7.8
46 0.8 3.1 0.9 2.4 10.2 12.4
92 3.1 12.4 2.6 8.9 15.9 21.6


Table 1 gives some values for the estimated stopping depth [from Eq. (6)] and the calculated stopping depth [according to Eq. (4)] by following a penetrating proton through an isothermal electron layer until it has lost all of its energy in excess with respect to the ambient electrons. The estimated stopping depth is generally higher than the numerical result due to the reasons mentioned above. For low electron temperatures these discrepancies are quite conspicuous.

The stopped protons accumulate at their stopping depth in the neutron star atmosphere. To maintain charge neutrality, an equal number of electrons has to move to the same location. The ADAF electrons entering the disk are stopped over a much short distance, compared with the protons. To make them move to the location where the protons come to rest in the disk, an electric field has to develop in the neutron star atmosphere, such that the resulting electron current has just this property. We assume that such an equilibrium actually develops, on the grounds that any charge imbalance will quickly lead to the buildup of a strong restoring electric field.

Equation (4) is only valid for non-relativistic conditions, whose validity needs to be checked for the high proton temperatures in an ADAF. Analytic expressions for relativistic temperatures exist only in special cases. Stepney (1983) has derived such an expression for the rate of transfer of energy between populations with relativistic Maxwellian distributions in terms of an integral over the scattering cross-section. In the case of hot protons heating cooler electrons the Rutherford cross-section is the relevant cross-section. An expression in closed form for the heating rate (in ergs cm-3s-1) is given by Stepney & Guilbert (1983),

   \begin{displaymath}% q_{\rm Stp} = -{\displaystyle\frac{3 m_{\rm e}}{2 m_{\rm p... ...})}{K_2(1/\theta_{\rm e})K_2(1/\theta_{\rm p})}}\; \ln\Lambda \end{displaymath}


\begin{displaymath}\times \left[{\displaystyle\frac{2(\theta_{\rm e}+\theta_{\rm... ...theta_{\rm p}} {\theta_{\rm e}\theta_{\rm p}}}\right)\right], \end{displaymath} (7)

where $\theta_{\rm e}=k T_{\rm e}/m_{\rm e} c^2$, $\,\theta_{\rm p}=k T_{\rm p}/m_{\rm p} c^2$, and K0, K1, K2 are Bessel functions.

The proton-electron heating rate according to the stopping formula derived from Spitzer's theory is given by

 \begin{displaymath} q(v_{\rm p}) = n_{\rm e}\,n_{\rm p}\,v_{\rm p}\,\varepsilon(E_{\rm p}). \end{displaymath} (8)


  \begin{figure} \par\includegraphics[width=8.3cm,clip]{h2960_f1.ps}\end{figure} Figure 1: Energy loss rates of a hot thermal distribution of protons in a cooler ionized hydrogen plasma (with temperature indicated along the curves). Solid line: fully relativistic expression of Stepney & Guilbert (1983), dash-dotted line: the non relativistic expression from Spitzer's (1962) theory. Dashed line: Spitzer's formula, but averaged over a relativistic instead of a classical Maxwellian. Spitzer's treatment is accurate to better than 5% for proton temperatures <100 MeV.
Open with DEXTER

By averaging Eq. (8) over a thermal velocity distribution we obtain a rate that can be compared directly with Eq. (15) from Stepney & Guilbert. We do this averaging both with a classical Maxwell distribution and a relativistic Maxwell distribution. The results of this comparison are shown in Fig. 1. The heating rates are almost indistinguishable up to proton temperatures of roughly 100 MeV. At higher proton temperatures the deviations increase steadily. But for the proton temperatures used in our simulation ( $k T_{\rm p}<100$ MeV) Spitzer's formalism is accurate to better than 5% compared to Stepney's result. Therefore we are confident that the usage of the classical stopping formula following Spitzer (1962) together with a classical Maxwell distribution for the protons will yield the correct results for to energy transfer from virialized protons to electrons near a neutron star surface.

   
2.2 The model atmosphere

We compute the density distribution through our plane-parallel one dimensional model atmosphere from the equation of hydrostatic equilibrium,


 \begin{displaymath} \frac{{\rm d} P}{{\rm d} \tau} = \frac{g}{\kappa_{\rm es}} + \frac{\partial P_{\rm p}(\tau)}{\partial\tau}\cdot \end{displaymath} (9)

Here $g = \frac{G M_{*} }{ R_{*}^2}$ is the gravitational acceleration, which we treat as constant as the atmosphere is thin compared to the neutron star radius, and $\kappa_{\rm es} = 0.40$cm2 g-1 is the electron scattering opacity. $\frac{\partial P_{\rm p}}{\partial\tau}(\tau)$ denotes the force exerted by the deceleration of the protons as a function of optical depth. The change of momentum (exerted on the atmosphere) of an incoming proton with respect to the vertical can be expressed by


 \begin{displaymath} \Delta p_{\perp} = -m_{\rm p}\,\Delta v \cos \theta . \end{displaymath} (10)

The force $\partial P_{\rm p}/\partial\tau$ is evaluated by recording the change of velocity $\Delta v$ of the incoming protons as a function of optical depth. Since we only study cases with luminosities below the Eddington limit we neglect the contribution of the radiative pressure.

We solve our pressure profile by starting with $\tau_0=0$ at the top of our atmosphere and integrate Eq. (10) to the maximum optical depth of our model atmosphere $\tau_{\rm max}$. We usually set $\tau_{\rm max}=10$. At this optical depth the energy flux of the protons is already negligible.

2.3 Accretion from the ADAF

The protons from the ADAF have a temperature of the order of the virial temperature [Eq. (1)]. As mentioned above we use a Maxwell distribution for the proton velocities. Additionally, the ADAF has a rotation rate, which is somewhat smaller than the Keplerian rate $\Omega_{\rm K}$ (e.g. Narayan & Yi 1995). When computing the penetration of the protons we take the Maxwellian velocity component of the protons. Further we take into account the component of velocity tangential to the neutron star surface due to the ADAF rotation. Instead of using a detailed model of an ADAF in which these velocity components can in principle be determined quantitatively, we model it with a temperature $T=\xi T_{\rm vir}$ and tangential velocity $v=\eta R\Omega_{\rm K}$, i.e. we use the parameters $\xi$and $\eta$ to scale these quantities. Both the tangential and vertical velocity components add to the energy deposited, but are not equivalent since the thickness of the heated layer is smaller the more the protons enter tangentially.

   
2.4 The radiative transfer

The radiative transfer equation appropriate for our model is

 \begin{displaymath} \mu\frac{{\rm d}I_{\mu,\nu}}{{\rm d}z} = j^\mathit{ff}_{\nu}... ...u} - (\alpha^\mathit{ff}_{\nu}+\alpha^{c}_{\nu}) I_{\mu,\nu}, \end{displaymath} (11)

where $I_{\mu,\nu}$ is the intensity as a function of frequency $\nu$and photon angle $\mu=\cos(\theta)$. The coordinate z is measured from the top of the atmosphere where the density is virtually zero (z=0) downwards into the neutron star atmosphere until $z=z_{{\rm max}}$ (corresponding to $\tau_{\rm max}$).

The bremsstrahlung absorption coefficient $\alpha^\mathit{ff}_{\nu}$ is given by

 \begin{displaymath} \alpha^\mathit{ff}_\nu = 3.7\times 10^8\;n_{\rm e} n_{\rm i}... ...c{1-{\rm e}^{-h\nu/kT}} {\nu^3\,\sqrt{T}} \,\bar g_\mathit{ff} \end{displaymath} (12)

in units of cm-1. The symbol $\bar g_\mathit{ff}$ represents the velocity averaged Gaunt factor for bremsstrahlung [see e.g. Rybicki & Lightman 1979]. The Compton extinction coefficient $\alpha^{c}_{\nu}$ is simply

 \begin{displaymath} \alpha^{c}_\nu = n_e \sigma^c_\nu, \end{displaymath} (13)

where $\sigma^c_\nu$ is the electron cross section with Klein-Nishina corrections. For the non-relativistic limit this equals the Thompson cross section $\sigma_{\rm T}$.

The bremsstrahlung emission $j^\mathit{ff}_{\nu}$ is given by

 \begin{displaymath} j^\mathit{ff}_{\nu} = 5.4\times 10^{-39} \;n_{\rm\rm e} n_{\rm i} \, \frac{{\rm e}^{-h\nu/kT}}{\sqrt{T}}\,\bar g_\mathit{ff} \end{displaymath} (14)

in units of ergcm-3 s-1 ster-1. The Compton "emissivity'' (i.e. radiation scattered to $\mu,\nu$ from some other $\mu',\nu'$) is

 \begin{displaymath} j^{c}_{\mu,\nu} = n_{\rm\rm e} \int \left(\frac{\nu}{\nu'}\r... ...igma_c\,(\mu',\nu'\rightarrow\mu,\nu)\,{\rm d}\mu'{\rm d}\nu'. \end{displaymath} (15)

Here $\sigma_{\rm c}(\mu',\nu'\rightarrow\mu,\nu)$ is the Compton Scattering Kernel (CSK). An efficient algorithm for run-time computation of the CSK is described by Kershaw et al. (1986) and Kershaw (1987). Both $\sigma^c_\nu$ and $\sigma_{\rm c}(\mu',\nu'\rightarrow\mu,\nu)$depend on the electron temperature.

Whenever the temperatures of our solutions are low enough, $kT\ll m_{\rm e}c^2$, we can simplify the Compton emissivity Eq. (15) without compromising accuracy. We introduce the angle-averaged Compton cross section $\sigma_{\rm c}(\nu'\rightarrow\nu)$ defined as

\begin{displaymath}\sigma_{\rm c}(\nu'\rightarrow\nu) \equiv \frac{1}{2} \int \sigma_{\rm c}(\mu',\nu'\rightarrow\mu,\nu) {\rm d}\mu'{\rm d}\mu. \end{displaymath} (16)

A semi-analytic expression for $\sigma_c(\nu'\rightarrow\nu)$ is given by (Poutanen & Svensson 1996). The angle-averaged intensity $J_\nu$ is defined as

\begin{displaymath}J_\nu = \frac{1}{2}\int I_{\mu,\nu} {\rm d}\mu. \end{displaymath} (17)

By replacing $\sigma_{\rm c}(\mu',\nu'\rightarrow\mu,\nu)$ with its average $\sigma_{\rm c}(\nu'\rightarrow\nu)$ in Eq. (15) one obtains an approximative expression for the Compton emissivity

 \begin{displaymath} j^{c}_{\nu} = n_{\rm e} \int \left(\frac{\nu}{\nu'}\right) J_{\nu'} \sigma_{\rm c}(\nu'\rightarrow\nu){\rm d}\nu', \end{displaymath} (18)

which is independent of $\mu$, and which involves only a single integral instead of the double integral of Eq. (15).

Equations (11)-(15), (18) constitute the complete set of equations for the radiative transfer. Since Eq. (11) and Eqs. (15), (18) are mutually dependent, the solution cannot be found by direct evaluation. To solve the system, we use a standard Lambda Iteration scheme [see e.g. Rutten 1999 and references therein; see also Zane et al. 1996; Poutanen & Svensson 1996].

We discretize frequency $\nu$ with equal spacing in $\log\nu$, and photon angle $\mu=\cos(\theta)$ according to the roots of the Legendre polynomials of order n. At the start of the procedure we choose $j^c_{\mu,\nu}=0$, or any other initial guess that might be appropriate. For each $\nu_i$ and $\mu_j$ we now integrate the transfer equation (Eq. (11)) from z=0 to $z=z_{{\rm max}}$ for $\mu_j>0$, or from $z=z_{{\rm max}}$ to z=0 for $\mu_j<0$. At z=0 we impose zero inward flux as boundary condition, while at $z=z_{{\rm max}}$ we choose a Planck function at temperature $T_{{\rm BC}}$ as the starting value for the integration. The temperature $T_{{\rm BC}}$ of this Planck function is chosen such that:

\begin{displaymath}\sigma T^4_{{\rm BC}} = 2\pi\int_0^1 \mu\, {\rm d}\mu \int_0^\infty {\rm d}\nu\; I_{\mu,\nu}(z=z_{{\rm max}}) \end{displaymath} (19)

so that at the bottom of the atmosphere the net flux is zero.

After performing this operation for every $\nu_i$ and $\mu_j$, we can evaluate the new $j^c_{\mu,\nu}$ by employing Eqs. (15), (18). This new $j^c_{\mu,\nu}$ is inserted into Eq. (11) and the integration of Eq. (11) is repeated along the lines sketched above. After each iteration the relative difference of the solution with the previous iteration step is computed. If this relative difference drops below 10-3, we assume that the solution has converged. Since in our solutions the scattering optical depths remains always of the order of a few, we need not worry about the well-known convergence problems of the lambda iteration procedure sketched above.

Once the radiative transfer solution is obtained, the radiative heating and cooling rates can be evaluated by flux differences. The flux at each position x is defined as

\begin{displaymath}F(z) = 2\pi\int_{-1}^{+1} \mu\,{\rm d}\mu \int_0^{\infty} {\rm d}\nu \;I_{\mu,\nu}(x). \end{displaymath} (20)

The net cooling rate is then

 \begin{displaymath} \Lambda_{\rm rad}^{-}(z) = \frac{{\rm d}F(z)}{{\rm d}z}\cdot \end{displaymath} (21)

The radiative transfer is solved on a grid set up in terms of Thomson optical depth when solving the pressure profile (Sect. 2.2). This grid is then transferred to the vertical zcoordinates for the radiative transfer code by the usual

\begin{displaymath}z(\tau) = \int_0^{\tau_{\rm max}}(n_{\rm e}(\tau) \sigma_{\rm T})^{-1} {\rm d} \tau. \end{displaymath} (22)

Thus we immediately obtain the net cooling rates $\Lambda_{\rm rad}^{-}(\tau)$ as a function of Thomson optical depth.

We also checked our solution for the importance of pair processes, calculating the pair number densities according to 1998. We find that in our model the electron temperatures are too low and the electron densities too high for pair production to become important.

  \begin{figure} \par\includegraphics[width=12cm]{h2960_f2.eps} \end{figure} Figure 2:  Upper left: emergent spectrum (solid line), thermalized downward directed spectrum (dashed line) and blackbody input spectrum at the base (dotted - dashed line); upper right panel: combined heating rates from proton heating and electron thermal conductivity (solid line) and radiative cooling rates due to Comptonization and bremsstrahlung (squares) - the dotted line shows the rates from electron conductivity alone; lower panels show from left to right and from top to bottom the electron temperature $T_{\rm e}$, pressure P, electron density $n_{\rm e}$ and the geometrical depth of the layer z, for a solution with $F_{\rm p} = 1.6\times 10^{24}$ ergcm-2 s-1, $T_{\rm p} = \frac{1}{2} T_{\rm vir}$ and $\Omega = \frac{1}{2}\Omega_{\rm K}$.
Open with DEXTER


  \begin{figure} \par\mbox{ \includegraphics[keepaspectratio=false,width=\columnwi... ...hics[keepaspectratio=false,width=\columnwidth]{h2960_f3b.ps} }\par\end{figure} Figure 3: Emergent model spectra and temperature profiles at a fixed proton energy flux $F_{\rm p} = 1.6\times 10^{24}$ ergcm-2 s-1 for various initial proton temperatures, $\xi \,T_{\rm vir}$, with $\xi =0.25$,0.5,0.75, and $\eta =0.5$. Higher proton temperatures produce lower electron temperatures in the hot part of the atmosphere and softer emergent spectra. The dotted dashed line shows the corresponding blackbody spectrum with the effective temperature of the layer, $T_{\rm eff} = (F_{\rm p}/\sigma )^{1/4}$.
Open with DEXTER


  \begin{figure} \par\mbox{ \includegraphics[keepaspectratio=false,width=\columnwi... ...egraphics[keepaspectratio=false,width=\columnwidth]{h2960_f4b.ps} } \end{figure} Figure 4: Emergent spectra and temperature profiles for various values of the proton energy flux (in ergcm-2 s-1, denoted by the numbers of the lines) at constant proton temperature, $\xi =0.5$. With decreasing energy flux the optical depth of the heated layer shrinks.
Open with DEXTER


  \begin{figure} \par\mbox{ \includegraphics[keepaspectratio=false,width=\columnwi... ...egraphics[keepaspectratio=false,width=\columnwidth]{h2960_f5b.ps} } \end{figure} Figure 5: Emergent spectra and temperature profiles at low proton energy flux (in ergcm-2 s-1, denoted by the numbers of the lines) and constant proton temperature, $\xi =0.5$. The hot part is extremely thin in terms of Thomson optical depths. The dotted dashed lines show the corresponding blackbody spectra with the effective temperature of the layer.
Open with DEXTER

2.5 Energy balance from heating and cooling

We start our calculations with an isothermal atmosphere in hydrostatic equilibrium according to Sect. 2.2. For the initial electron temperature we set $T_{\rm e}$ = 1 keV throughout the layer. After each time step we obtain the heating rates $\Lambda_{\rm p}^+(\tau)$ from the Coulomb interactions and the cooling rates $\Lambda_{\rm rad}^{-}(\tau)$ due to the radiative processes bremsstrahlung and multiple Compton scattering as a function of optical depth. The time step of every cycle in our simulation is adjusted to the shortest energy exchange time scale occurring in the calculation.

Additionally we include the energy redistribution due to electron thermal conductivity. Generally the flow of heat per unit area, Q, in the presence of a temperature gradient, is given by (see e.g. Spitzer 1962)


 \begin{displaymath} Q \approx -\kappa\,\nabla T. \end{displaymath} (23)

The heat conductivity $\kappa$ for a fully ionized gas is given by

 \begin{displaymath} \kappa \approx 1.85 \times 10^{-5}\,\frac{T^{5/2}}{\ln\Lambda} \;\; {\rm erg\;s}^{-1}~{\rm K}^{-1}~{\rm cm}^{-1}. \end{displaymath} (24)

Thus the change of enthalpy $\partial w/\partial t$ due to electron conductivity can be expressed by the divergence of the heat flux Q,

 \begin{displaymath} \frac{\partial w}{\partial t} = \rho\, c_{\rm p} \,\frac{\p... ... \nabla \cdot (\kappa \, \nabla T) =: \Lambda_{\rm cond}(\tau) \end{displaymath} (25)

where $c_{\rm p}$ is the specific heat at constant pressure.

Now we can calculate the total change of enthalpy per time step,


 \begin{displaymath} \frac{\Delta w(\tau)}{\Delta t} = \rho c_{\rm p} \frac{\Delt... ...(\tau) + \Lambda_{\rm rad}^-(\tau) + \Lambda_{\rm cond}(\tau). \end{displaymath} (26)

With Eq. (26) we are able to calculate the change of temperature and thus the new temperature profile after each time-step. With the new temperature profile we can update the hydrostatic structure according to Sect. 2.2. We follow the simulation until the Coulomb heating is balanced by the radiative cooling and energy redistribution by conductivity and the temperature and density structure as well as the spectra have reached a stationary state.

   
3 Results of the model computations

For our model calculations we use fixed values for the mass $M_{\rm *} = 1~M_{\odot}$ and the radius $R_{\rm *} = 10^6$ cm for the neutron star. For the rotational angular velocity of the accreting material we use $\eta =0.5$. The results depend on the proton energy flux per unit area of the neutron star surface. In order to make the results interpretable in terms of a total luminosity, we assume that 60% of the neutron star surface is involved in the accretion process. This takes into account the approximate vertical extent of an ADAF (Popham & Sunyaev 2001).

3.1 Dependence on proton energy flux and temperature

Figure 2 shows an equilibrium solution for a proton energy flux $F_{\rm p} = 1.6\times 10^{24}$ ergcm-2s-1, $\xi=\eta=0.5$, yielding a luminosity of $L_{\rm *}=0.1~L_{\rm Edd}$. The neutron star atmosphere is clearly divided into two parts. A hot surface layer is separated by a sharp temperature front from the much cooler bottom part. The maximum temperature in the hot part of the atmosphere is $T_{\rm e}\approx50$keV. This is considerably hotter than previous calculations showed (e.g. Alme & Wilson 1973; Turolla et al. 1994) but also significantly cooler than the "hot'' solutions from Turolla et al. (1994) and Zane et al. (1998). We did not find any comparable hot solutions.

Figure 2 also shows that energy redistribution through electron conductivity gets important at the transition zone from the hot to the cool part. The emergent model spectrum in Fig. 2 illustrates that the hot part acts as an effective Comptonization layer with an optical thickness of order unity. The downward directed photons are almost completely thermalized in the cold part. This justifies our assumption of Sect. 2.4 that the downward flux of photons is converted into a black body flux of upward soft photons of temperature $T_{{\rm BC}}$ at the base of our model since almost all the thermalization has already taken place.

Figure 3 shows the emergent spectra and temperature profiles through the atmosphere at a fixed proton energy flux level $F_{\rm p} = 1.6\times 10^{24}$ ergcm-2 for various proton temperatures, $\xi T_{\rm vir}$. The higher the initial proton temperature the lower is the electron temperature in the hot part of the atmosphere and vice versa for the cold part. The optical depth of the hot part is of order unity in every case. As a consequence harder spectra are produced with lower initial proton temperatures. This can be explained with the decreasing penetration depth of the protons at smaller velocities: the heating rates in the upper part of the atmosphere are increased and thus the temperature is higher there.

Figures 4 and 5 show the dependence of the solution on the proton energy flux at constant proton temperature, $\xi =0.5$. At low proton energy flux (Fig. 5) the optical depth of the hot part is small ( $\tau_{\rm hot} < 10^{-2}$). The solutions for these cases are comparable to those found by Zampieri et al. (1995). But the temperature of the hot part in our solutions turns out to be much higher than in their solutions. This hot part causes an (energetically unimportant) bump in the emergent spectra at high frequencies due to optically thin bremsstrahlung emission in this layer. As in Zampieri et al. (1995) the emergent spectra are significantly harder than the blackbody at the layer effective temperature, $T_{\rm eff} = (F_{\rm p}/\sigma )^{1/4}$. We also see a harder spectrum with respect to the blackbody of the layer effective temperature at the lowest luminosity.

By further increasing the proton energy flux the temperature jump moves inward to greater optical depth reaching order unity at the highest flux levels. Thereby the temperature in the hot part is roughly kept at a constant value in contrast to the temperature of the cold part, which increases continously. At $F_{\rm p} = 1.6\times 10^{24}$ ergcm-2 s-1 the hot part acts as an effective Comptonization layer for the cool thermal photons from the interior. In addition the electron thermal conduction smoothes the temperature profile and the shock is not as sharp as in the low energy flux cases.

3.2 A reverse photosphere effect

Evident in our model spectra is an excess of soft photons (see e.g. Figs. 2, 3, 5) with respect to a blackbody with the effective temperature of the corresponding energy flux. Though energetically unimportant it is an interesting feature of our model spectra. This soft photon excess is caused by the reverse temperature profile of the atmosphere, i.e. the high temperatures are above the cooler parts (in contrast to e.g. the solar photosphere). Figure 6 shows the absorption optical depth as a function of the photon energy and electron temperature for an equilibrium solution with $\xi =0.5$ and $\eta =0.5$. The line for $\tau _{\rm a}=1$ is emphasized. For low photon energies $\tau _{\rm a}=1$ is found high in the neutron star atmosphere, where the temperature is high. For high frequency photons the atmosphere gets more and more transparent and $\tau _{\rm a}=1$ occurs at lower temperatures. Thus the high frequency photons decouple at lower temperatures and propagate outwards. As a result there is an excess of low energy photons, which decouple from the atmosphere at higher temperatures.

The same arguments explain the harder spectra at low proton energy fluxes. Here the temperature profile in the atmosphere is dominated by the increase of temperature with optical depth, as has already been explained by Zampieri et al. (1995). Thereby the high energy photons decouple in the deeper, hotter layers and propagate out freely. This contribution to the high energy part of the spectrum causes the harder spectrum compared to the blackbody of the layer effective temperature.

  \begin{figure} \par\includegraphics[width=8.3cm,clip]{h2960_f6.ps}\end{figure} Figure 6: Absorption optical depth as a function of photon energy and atmospheric temperature for an equilibrium solution with $\xi =0.5$ and $\eta =0.5$. The line depicting $\tau _{\rm a}=1$ is emphasized.
Open with DEXTER

   
4 Discussion and conclusions

The main conclusion from our work is that a neutron star surface embedded in a hot ADAF-type accretion flow does not act like a simple blackbody thermalizer. Instead the interaction with the hot ions produces spectra with pronounced high energy, comptonized tails.

We have considered the accretion of protons onto the surface of a neutron star for a wide range of proton energy fluxes. In our model the stopping depth of the protons is obtained self consistently. We allow for the scaling of the proton temperature $T_{\rm p} = \xi T_{\rm vir}$, which in our model controls the vertical velocity component of the protons. Further we consider accretion from an ADAF, where the protons have a tangential velocity component $v=\eta \Omega_{\rm K}$. By variation of the proton temperature and the proton energy flux a range of X-ray spectra can be produced. Our model computations show that the local proton energy flux determines the optical depth of the hot part whereas the temperatures through the atmosphere are determined by the initial proton temperatures.

Previous models considered spherical accretion onto the neutron star. The free infall velocity of the protons was used in such models. If accretion takes place from an accretion disk or an ADAF, the vertical velocity component is considerably smaller than than the free infall velocity. This has important consequences on the emergent spectra as our solutions suggest: at the same local energy flux a decrease in the vertical proton velocity increases the temperature of the hot part of the atmosphere, and X-ray spectra similar to those of neutron stars in their hard states can be produced. Figure 7 shows a comparison of a model spectrum with observed spectra of X-ray bursters in their low/hard states (taken from Gilfanov et al. 1998).

A question left unanswered by our studies is how the spectrum from the surface is modified by the overlying accretion flow. E.g. an overlying optically thick boundary layer might significantly thermalize hard photons from the neutron star surface again. But we can conclude that the input spectrum at the bottom of a boundary layer is not Planckian at high proton energy fluxes. If the boundary layer is optically thin the contribution to the observed hard spectral component at high proton energy fluxes might be quite important. A discussion whether there is an optically thin boundary layer can be found in e.g. Barret et al. (2000). For low accretion rates the boundary layer is indeed optically thin (e.g. King & Lasota 1987).

Detailed studies about the accretion disk boundary layer around a NS have been performed by Popham & Sunyaev (2001). They show that the accretion disk near the NS is radially and vertically extended and that the angular velocity $\Omega$ is reduced with respect to the Keplerian value. Such a boundary layer naturally produces a proton illumination scenario as described in our model. Their boundary layer is optically thin to absorption. For high luminosities this region has a radial extent larger than one stellar radius, i.e. the hard radiation from the NS surface, produced by proton illumination, can easily propagate outwards without being thermalized. The boundary layer even further enhances comptonization as the gas there is hot ($\gtrsim $108 K).

Our results emphasize the importance of the neutron star surface for the contribution to the hard spectral component in low mass X-ray binaries (LMXBs). A prerequisite for the occurrence of a hard spectral component is a high local proton energy flux. The local proton energy flux onto the neutron star surface depends on both the mass accretion rate and the size of the accretion belt on the neutron star. Our model computations have local validity, i.e. if the involved accretion belt is small, high proton energy fluxes on the neutron star surface can even occur at accretion rates well beneath the Eddington limit. Therefore our model does not exclude the possibility of hard spectral components from the neutron star surface at low luminosities, as is observed.

  \begin{figure} \par\includegraphics[width=7.8cm,clip]{h2960_f7.ps}\end{figure} Figure 7: Comparison of an emergent model spectrum with $\xi =0.25$ and $\eta =0.5$ with observed neutron star spectra in their hard states. The spectra are arbitrarily shifted in the vertical direction for clarity.
Open with DEXTER

Acknowledgements
We thank Roberto Turolla for the usage of the CSK subroutines. We are grateful to Marat Gilfanov for the neutron star spectra. This research has made use of data obtained through the HEASARC Online Service, provided by the NASA/GSFC. This work was done in the research network "Accretion onto black holes, compact stars and proto stars'' funded by the European Commission under contract number ERBFMRX-CT98-0195.

References

 
Copyright ESO 2001