3 Heating and cooling processes of the model

   
3.1 Proton illumination of the accretion disk

3.1.1 Coulomb interactions in an ionized plasma

At the typical energies of the protons incident on the cool disk, the energy loss is mostly by long-range Coulomb interactions with the electrons in the disk (small-angle scattering on the large number of electrons in a Debye sphere). This is opposite to the case of protons with a temperature near that of the plasma in which they move. In the latter case, the equilibration among the protons is faster than between electrons and protons, by a factor of order $(m_{\rm p}/m_{\rm e})^{1/2}$. To see how this apparent contradiction is resolved, consider the basic result for the energy loss of a charged particle moving in a fully ionized, charge-neutral plasma. This was derived by Spitzer (1962) (making use of Chandrasekhar's (1942) earlier results on dynamical friction). Introduce as a measure of distance in the plasma the Thomson optical depth $\tau$, i.e. ${\rm d}\tau=\sigma_{\rm T}n {\rm d}l$, where $\sigma_{\rm T}=8\pi e^4/(3m_{\rm e}^2c^4)$ is the Thomson cross section, n the electron density and l the distance. The rate of change of energy $E=m_{\rm p}v^2/2$ of a proton moving with velocity v in a field of particles with charge e and mass $m_{\rm f}$ (the "field particles'' in Spitzer's nomenclature) is then given by Spitzer's Eq. (5)-(15). In our notation, this can be written in terms of the energy loss length $\tau_{\rm f}$ for interaction with the field particles $\rm f$,

$$\tau_{\rm f}^{-1}\!=\!{1\over E}\!\left({{\rm d}E\over {\rm d... ...1+{m_{\rm p}\over m_{\rm f}}\right)\!F\left(x_{\rm f}\right),$$ (3)

where $\ln\Lambda$ is the Coulomb logarithm (which is determined by the size of the Debye sphere). Here

$$F(x)=\psi(x)-x\psi^\prime (x),$$ (4)

where $\psi$ is the error function and $\psi^\prime$ its derivative, and $x_{\rm f}= v[{m_{\rm f}}/(2kT)]^{1/2}$ is (up to a numerical factor) the ratio of the incident proton's velocity to the thermal velocity of the field particles. The limiting forms of F are

$$F(x)\rightarrow x^3/3\quad (x\rightarrow 0),\qquad F\rightarrow 1\quad (x\rightarrow\infty).$$ (5)

We can evaluate (3) under the assumption that the field particles that are most relevant for the energy loss are the protons or the electrons, respectively, and compare the loss lengths. If the incident proton has velocity comparable with the thermal velocity of the field protons, we have $x_{\rm p} \sim 1$, and $x_{\rm e}=(m_{\rm e}/m_{\rm p})^{1/2}x_{\rm p}\ll 1$. $F(x_{\rm p})$ is then of order unity and $F(x_{\rm e})\approx x^3_{\rm e}/3$. Setting ${\rm f}={\rm e}$ respectively ${\rm f=p}$ in (3) and taking the ratio, we have

$${\tau_{\rm e}\over\tau_{\rm p}}\approx \left({m_{\rm p}\over m_{\rm e}}\right)^{1/2}\cdot$$ (6)

The loss length for interaction with the electrons is thus much longer than for interaction with the protons, and the interaction with electrons can be neglected. This well known result is the relevant limit for the relaxation of a proton distribution in a plasma that is not too far from its thermal equilibrium.

For incoming protons of high energy, however, the result is different because $x_{\rm e}$ is not sufficiently small any more. In the high-v limit, $F(x_{\rm e})= F(x_{\rm p})=1$, and one has

$${\tau_{\rm e}\over\tau_{\rm p}}\approx 2{m_{\rm e}\over m_{\rm p}}=10^{-3}.$$ (7)

In this limit, the energy loss is thus predominantly to the electrons. A related case is that of the ionization losses of fast particles in neutral matter (for references see Ryter et al. 1970). The case of an ionized plasma is simpler, since the electrons are not bound in atoms. The change from proton-dominated loss to electron-dominated loss takes place at an intermediate velocity $v_{\rm c}$, at which $F(x_{\rm p})\approx 1$ but $x_{\rm e}$ still small, so that $F(x_{\rm e})\approx x^3_{\rm e}/3$. Equating $\tau_{\rm e}$ and  $\tau_{\rm p}$ then yields

$$x_{\rm e,c}\approx \left({m_{\rm e}\over m_{\rm p}}\right)^{1/3},$$ (8)

or

$${E_{\rm c}\over kT}\approx \left({m_{\rm p}\over m_{\rm e}}\right)^{1/3}\approx 12.$$ (9)

For the electron temperatures we encounter in our models, $T\sim 100$ keV, energy loss to the field protons can thus be neglected for incoming protons with energy $E\ga 1$ MeV. This is the case in all calculations presented here.

3.1.2 Corrections at high and low energies

Spitzer's treatment is non-relativistic, while virialized ADAF protons near the hole can reach sub-relativistic temperatures. A fully relativistic treatment of the Coulomb interactions in a plasma has been given by Stepney & Guilbert (1983). We have compared the classical treatment according to Spitzer's theory with this relativistic result in Deufel et al. (2001), and found it to be accurate to better than 5% for proton temperatures <100 MeV. The classical approximation in Spitzer's analysis therefore does not introduce a significant error for the problem considered here.

For high energies, the Coulomb energy loss becomes so small that loss by direct nuclear collisions becomes competitive. This happens (cf. Stepney & Guilbert 1983) at $E\ga 300$ MeV, an energy that is not reached by virialized protons except in the tail of their distribution. We ignore these direct nuclear collisions. Note, however, that a gradual nuclear processing by such collisions can be important (Aharonian & Sunyaev 1984), in particular for the production of the lithium. The lithium overabundances seen in the companions of LMXB (Martín et al. 1994a), may in fact be a characteristic signature of the interaction of an ADAF and a disk described here (Martín et al. 1994b; Spruit 1997).

As the protons slow down, they eventually equilibrate with the field protons. This last part of the process is not accurately described by the energy loss formula (3). In addition to the simple energy loss of a particle moving on a straight path through the plasma, one has to take into account the random drift in direction and energy resulting from the interaction with the fluctuating electric field in the plasma. This drift can be ignored to first order (end of Sect. 5.2 in Spitzer 1962), but takes over in the final process of equilibration with the plasma. This last phase involves negligible energy transfer compared with the initial energy of the protons in the present calculations, and can be ignored here.

3.1.3 Charge balance

The protons penetrating into the disk imply a current that has to be balanced by a "return current''. As in all such situations, this return current results from the electric field that builds up due to the proton current. This field drives a flow of electrons from the ADAF to the disk which maintains the charge balance. Since the electron density in the disk is high, the return current does not involve a high field strength.

3.1.4 Calculations

The numerical method used to calculate the electron heating rates of the accretion disk due to proton illumination is described in detail in Paper I; here we give a short overview.

We follow the evolution of an initially Maxwellian distribution of protons, which we place above the cool accretion disk, as the protons penetrate through the disk atmosphere. For the temperature of the incident protons in our model we take the local virial temperature

$$T_{\rm p}=T_{\rm vir}=\frac{G M m_{\rm p}}{3 k_{\rm B} R}\approx \frac{156}{r}\hspace{2mm}\mbox{MeV},$$ (10)

where G is the gravitational constant, M the mass of the black hole, $m_{\rm p}$ the proton mass, R the radius from the central object and $k_{\rm B}$ is the Boltzmann constant. r is the dimensionless radius, $r=R/R_{\rm S}$, where $R_{\rm S} = 2GM/c^2$ is the Schwarzschild radius.

The energy loss of the incoming protons is computed from (3), with ${\rm f}={\rm e}$. Denote by $\tau$ the Thomson optical depth measured vertically into the layer, from its top. The rate of change of energy with depth for a proton at depth $\tau$ moving at an angle $\theta$ with respect to the vertical is then

$${1\over E}{{\rm d}E\over {\rm d}\tau}={3\ln\Lambda\over\cos\t... ...m_{\rm p}}\left({c\over v}\right)^4[\psi(x)-x\psi^\prime (x)],$$ (11)

with

$$x=x_{\rm e}=v\left({m_{\rm e}\over 2kT}\right)^{1/2},$$ (12)

while $T(\tau)$ is the depth-dependent temperature of the layer, and v the depth and angle-dependent velocity of the incident proton. The angular distribution of the protons is discretized in 50 equidistant points. The Debye-length in the plasma, measured in Thomson optical depths, is quite small so that temperature variations over a Debye-length can be ignored and the local energy loss rate at any depth is adequately described by (11).

The energy flux of the incident flux $q_{\rm p}$ is proportional to the density in the ADAF, which depends on its accretion rate and viscosity. Instead of using detailed models we parameterize these dependences by scaling the energy flux of the incident protons with the local energy dissipation rate in the ADAF. Using the thin-disk expression for this dissipation rate, we set

$$q_{\rm p}=fQ(R) = f\cdot\frac{3GM\dot{M}}{8\pi R^3}\cdot J(R),$$ (13)

where f is an adjustable numerical factor of order unity and $J(R)=1-(R_{\rm i}/R)^{1/2}$. This scaling of the proton flux is also useful if one has a coronal model in mind instead of an ADAF. The factor f then is interpreted as the fraction of the viscous energy release that is dissipated in the corona.

We follow the protons from the corona numerically through the disk atmosphere and record the energy loss as a function of optical depth $\tau$ according to Eq. (11). This yields the local time dependent heating rate $\Lambda_{\rm P}(\tau)$.

   
3.2 Radiative transfer

We treat the radiative transfer by solving the radiative transfer equation

$$\mu\frac{{\rm d}I_{\mu,\nu}}{{\rm d}z} = j^{\it ff}_{\nu} + ... ...u} - (\alpha^{\it ff}_{\nu}+\alpha^{\rm c}_{\nu}) I_{\mu,\nu}$$ (14)

where $I_{\mu,\nu}$ is the intensity as a function of frequency $\nu$and photon angle $\mu=\cos(\theta)$, $\alpha^{\it ff}_{\nu}$ is the bremsstrahlung absorption coefficient, $\alpha^{\rm c}_{\nu}$ the Compton extinction coefficient, $j^{\it ff}_{\nu}$ the bremsstrahlung emissivity and $j^{\rm c}_{\mu,\nu}$ the Compton emissivity. The method for solving the radiative transfer equation is described in Deufel et al. (2001). From the solution of the radiative transfer equation we obtain the radiative cooling rates $\Lambda_{\rm rad}(\tau)$.

3.3 Pair processes

As discussed below, the cases of moderate optical depth can become rather hot so that pair processes have to be included. We include photon-photon pair production ( $\gamma\gamma\rightarrow {\rm e}^{+}{\rm e}^{-}$) as well as pair production due to ee collisions ( ${\rm ee}\rightarrow {\rm e\, ee}^{+} {\rm e}^{-}$) in a steady state, i.e. the pair production rate equals the pair annihilation rate.

The photon-photon pair production rate can be expressed as a quadruple integral over the dimensionless photon energies $x_1=h\nu_1/m_{\rm e}c^2$, $x_2=h\nu_2/m_{\rm e}c^2$ and photon angles $\mu_1$, $\mu_2$, respectively (cf. Zane et al. 1995):

$$R_{\gamma\gamma} = \frac{4\pi}{h^2c} \int {\rm d}x_1 {\rm d}... ...ac{I_{\mu_1,\nu_1}I_{\mu_2,\nu_2}}{x_1^2x_2^2} F(x_{+},x_{-}),$$ (15)

where F(x₊,x_-) is defined as

$$F(x_{+},x_{-}) = \int_{x_{-}}^{x_{+}} \frac{\sigma(x)~x^3{\rm d}x}{\sqrt{(x_{+}^2-x^2)(x^2-x_{-}^2)}},$$ (16)

with $x_{\pm}^2\equiv x_1x_2[1-\cos(\theta_1\pm\theta_2)]/2$(Stepney & Guilbert 1983). The pair cross-section $\sigma(x)$ is (Jauch & Rohrlich 1976)

$$\sigma(x)\!= \!\frac{1}{2}(1-\beta^2)\!\left[(3-\beta^4)\log\... ...ft( \frac{1+\beta}{1-\beta}\right)\!-2\beta(2-\beta^2)\right],$$ (17)

with $\beta\equiv\sqrt{x^2-1}/x$.

In our models we have simplified the above integral by replacing the intensities by their corresponding mean-intensities: $I_{\mu_1,\nu_1}\rightarrow J_{\nu_1}$ and $I_{\mu_2,\nu_2}\rightarrow J_{\nu_2}$. For any given value of $y\equiv x_1x_2$ we can now evaluate the function $G(y)\equiv\int F(x_{+},x_{-}) {\rm d}\mu_1 {\rm d}\mu_2$, which is independent of the radiation field, and which can be tabulated and stored beforehand. The pair production rate now reduces to

$$R_{\gamma\gamma} \approx \int G(x_1x_2) J_{\nu_1}J_{\nu_2} {\rm d}x_1 {\rm d}x_2\;.$$ (18)

Thus the quadruple integral is reduced to a double integral, which is much faster to evaluate. An accurate evaluation of this integral can be made by writing ${\rm d}x_1{\rm d}x_2=x_1x_2 {\rm d}\log(\sqrt{x_1x_2}) {\rm d}\log(x_1/x_2)$, since the function G(x₁x₂) is independent of x₁/x₂. The pair production rate due to electron-electron collisions is e.g. given in Svensson (1982) as

$$\dot{n}_{+}^{\rm ee}=(n_{+}+n_{-})^2 c ~r_{\rm e}^2 ~P_{\rm ee}(\theta)~$$ (19)

where we have introduced $\theta = k T_{\rm e}/ m_{\rm e}c^2$ and

$$P_{\rm ee} = \frac{1}{2}~\alpha_f^2~\frac{28}{27\pi}(2 \ln 2 \theta)^3.$$ (20)

Here $\alpha_f$ is the fine structure constant. This production rate is only included for $\theta>0.6$.

The pair annihilation rate is also given by Svensson (1982) as

$$\dot{n}_{+}^{\rm an} = n_{+}~n_{-}~c~r_{\rm e}^2~A(\theta)~,$$ (21)

where

$$A(\theta)=\frac{\pi}{1+2\theta^2/\ln(2\eta\theta + 1.3)}~,$$ (22)

$\eta\simeq 0.56146$. Now we can set up the pair balance equation in the stationary case,

$$R_{\gamma\gamma} + n^2~c~r_{\rm e}^2 \; [ (2z+1)^2 P(\theta) - z(1+z) A(\theta)] = 0.$$ (23)

We have used the charge neutrality condition n_- = n₊ + n and z = n₊/n is the ratio of the pair number density to the proton number density n. Equation (23) is a quadratic equation for the pair fraction z. The positive root of this equation is

$$z(\tau)=\frac{1}{2}\left[-1 + \frac{\sqrt{A(\theta)~c~n^2~r_{... ...~r_{\rm e}~\sqrt{c}\sqrt{A(\theta) - 4P(\theta)}} \right]\cdot$$ (24)

The additional electrons from pair production serve as further scattering partners for the Coulomb collisions and the Compton scattering. Effects due to electron-positron bremsstrahlung are not included nor is a treatment of the annihilation line.

3.4 Energy balance from heating and cooling

We begin our computation with an isothermal temperature profile in hydrostatic equilibrium according to Sect. 2. The initial temperature throughout the layer is $T_{\rm d}$ (see Paper I) for the thick disks and we set $T_{\rm e}=1$ keV for disks with moderate optical depths. The solutions do not depend on the initial temperature profile. First we calculate 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. Additionally the energy redistribution due to electron thermal conductivity is included using Spitzer's classical value, as in Deufel et al. (2001). This process adds the contribution $\Lambda_{\rm cond}(\tau)$ to the energy balance. The validity of classical electron conduction was checked by evaluating the electron mean free path $\lambda$, which turns out to be of the order 10² cm or less in the black hole candidate cases shown in Fig. 1. This is much smaller than the temperature gradient length (more than 10⁴ cm), so the condition for validity of the electron formula used is amply satisfied. The same holds for the AGN cases.

The equilibrium is computed following the time evolution of the layer until a balance between the heating and cooling processes is obtained. In the optically thick models, the time scales for approaching equilibrium turn out to be a sharp function of depth in the model. To deal with this, an adaptive time stepping process is used in which the time step depends on both time and depth in the model. Stability of this process was obtained by scaling this step with the shortest of the energy exchange time scales associated with the contributing heating and cooling processes. With this procedure the time evolution of the model is not realistic, but the final equilibrium obtained is. For the cases with moderate optical depth we use a depth-independent time step. In these cases, the time evolution of the model is also realistic.

The total change of enthalpy per time step as function of optical depth can then be expressed by

$$\frac{\Delta w(\tau)}{\Delta t} = \rho c_{\rm p} \frac{\Delta... ...(\tau) + \Lambda_{\rm rad}^-(\tau) + \Lambda_{\rm cond}(\tau),$$ (25)

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

Now the change of temperature as a function of optical depth can be calculated. With the new temperature profile the hydrostatic structure is updated according to Sect. 2. We follow the simulation until the Coulomb heating is balanced by the radiative cooling and energy redistribution due to electron conductivity. The energy balance of all our computations is better than 10%.