A&A 386, 757-762 (2002)
DOI: 10.1051/0004-6361:20020234

Comptonization and synchrotron emission in 2D accretion flows

I. A new numerical solver for the Kompaneets equation

A. Hujeirat1,2,3 - M. Camenzind2 - A. Burkert1


1 - Max-Planck-Institut für Astronomie, 69117 Heidelberg, Germany
2 - Landessternwarte-Koenigstuhl, 69117 Heidelberg, Germany
3 - Institute of Applied Mathematics, University of Heidelberg, 69120 Heidelberg, Germany

Received 11 October 2001 / Accepted 31 January 2002

Abstract
A numerical approach is presented for solving the hydro-magnetic radiative transfer equation in 4D, in which the action of the Kompaneets operator and synchrotron cooling are taken into account. The equation is solved using the pre-conditioned Black-White-Brown line Gauss-Seidel semi-iterative method. This alternating approach enables the use of vector and parallel machines without loss of efficiency. The method can be used to constrain and search for density-, temperature- and magnetic-field distributions with or without outflows. This is appropriate for displaying spectra of 2D accretion flows that are consistent with observations of low-mass X-ray binaries.

Key words: methods: numerical - accretion, accretion disks - radiative transfer - black hole physics


1 Introduction

The majority of the observed Low Mass X-ray Binaries (LMXBs) radiate a significant fraction of their total luminosity in X- and $\gamma-$ray bands (van der Klis 2000). It is thought that repeated scattering of low-frequency photons by hot electrons under certain conditions could be the main mechanism underlying the energy gain of the observed hard X-rays (Shapiro et al. 1976). However, it is not clear what the source of the seed photons is. There are at least three possibilities. They either emerge from the cool disk, in which the photon energies boost during their travel through the hot overlying corona (Haardt & Maraschai 1991; Poutanan & Svensson 1996), or from the cool matter at the surface of the star in the case of a neutron star, or from synchrotron radiation resulting from electron gyration around magnetic field lines threading the disk.

Compton up-scattering of soft photons is most efficient in unsaturated Comptonization regions where the Compton-Y parameter is of order unity. This parameter acquires large values in optically thick media, and small values in the corona, implying that the corona-disk interaction region and/or the innermost region of the disk are most appropriate for this process to operate efficiently. Consequently, Comptonization in accretion flows is intrinsically a 2D problem, and therefore requires a multi-dimensional treatment.

So far, Comptonization has been considered under strong assumptions that allow separation of variables leading to the separation of the Kompaneets operator from the radiative transfer equation. Here, the radiative intensity is assumed to be time-independent, isotropic and the plasma is isothermal. In this case, the generation and Comptonization of photons can be described by a second order differential equation in the frequency space (Illarionov & Sunyaev 1972; Falten & Rees 1972; Katz 1976; Shapiro et al. 1976).

The aim of this paper is to present, for the first time, a numerical tool for solving the photon density equation, self-consistently taking into account the effects of advection, diffusion, collision of matter with radiation and Comptonization in 3D axisymmetric accretion flows.

The present solver is complementary to the time-dependent radiative-magnetohydrodynamical implicit solver (IRMHD; see Hujeirat 2001 and the references therein). Using IRMHD, the spatial distributions of the physical variables obtained can then be used to calculate the corresponding spectra which can be compared with observations.

2 The governing equation

Different accretion models display different spectra. It is therefore essential to perform a diagnostic study to analyse their consistency with observations. This however requires solving the 7D radiation transfer equation self-consistently. The equation reads:

\begin{displaymath}\frac{1}{c} \frac{ \partial{I}}{\partial t} + n \cdot \nabla ... ...acute{\Omega} {\rm d}\acute{E}} + \varepsilon^{\rm mod}_{\nu}, \end{displaymath} (1)

where $I=I(t,r,\theta,\varphi,\vartheta,\phi,\nu)$ is the radiative intensity which depends on time t, the spherical coordinates $(r,\theta,\varphi)$, two ordinates $(\vartheta,\phi)$that determine the direction of the photons on the unit sphere, and on the frequency $\nu$. $\kappa_{\nu}$ and $\sigma$ are the absorption and scattering coefficients. $S_{\nu}$ is a source function. $ I_{\rm int} = \int{CI~{\rm d}\acute{\Omega} {\rm d}\acute{E}}$ describes the scattering of photons through electrons, and C is the scattering kernel. $\varepsilon^{\rm mod}_{\nu}$ is the modified synchrotron emission.

To make the problem tractable, the following approximations have been employed:

Using the last approximation, $ I_{\rm int}$ can be expanded up to second order in $\epsilon$ which reduces it to the so-called Kompaneets operator (Payne 1980):

\begin{displaymath}I_{\rm int} \Longrightarrow{\cal K}_{\nu} = -\frac{\nu}{m_{\r... ...T\nu}{m_{\rm e}c^2}\frac{\partial^2}{\partial{\nu}^2} (\nu I). \end{displaymath} (3)

In this case, the radiative transfer equation with respect to a rest frame of reference reads:

 \begin{displaymath}\frac{1}{c}\left[\frac{ \partial{E_{\nu}}}{\partial t} + \n... ...{\nu}-E_{\nu}) + {\cal K}_{\nu} + \varepsilon^{\rm mod}_{\nu}, \end{displaymath} (4)

where $E_{\nu}=\frac{4\pi}{c}J= E(t,r,\theta,\nu)$, $\chi_{\nu} = \rho(\kappa_{\nu}+\sigma)$, $\lambda_{\nu}$ is the flux limited diffusion coefficient (Levermore & Pomraning 1981), which forces the radiative flux to adopt the correct form in optically thin and thick regions, i.e.,

\begin{displaymath}{\nabla \cdot \left[\frac{\lambda_{\nu}}{\chi_{\nu}}\nabla E_... ... {\sf if} \hspace*{0.75cm} \tau \ll 1. \end{array} \right. } \end{displaymath} (5)

It provides also a smooth matching in transition regions.

The above two different behaviour of the operator can be combined as follows:

\begin{displaymath}\nabla \cdot \frac{\lambda_{\nu}}{\chi}\nabla_{\nu} E_{\nu} \doteq \nabla \cdot \eta_{\rm r}\nabla E_{\nu}, \end{displaymath} (6)

where $\eta_{\rm r} = (1-\alpha) \frac{\nabla E_{\nu}}{ \vert\nabla E_{\nu}\vert} + \alpha \frac{1}{3\chi} ,$ $\alpha = {\rm e}^{-R_{\rm FLD}}$, and $R_{\rm FLD} = \nabla E/\rho(\kappa_{\nu}S_{\nu}+\sigma E_{\nu})$ (Hujeirat & Papaloizou 1998).

$\varepsilon_{\nu}^{\rm mod}$ in Eq. (4) corresponds to the modified synchrotron emission of photons by relativistic electrons gyrating around magnetic field lines, which reads:

\begin{displaymath}\varepsilon_{\nu}^{\rm mod} = \xi~\varepsilon_{\nu} + (1-\xi)\varepsilon^{BB}_{\nu}, \end{displaymath}

where $\xi (\doteq e^{-({{\nu}_{\rm c}}/{\nu})^2})$ is a switch on/off operator which bridges optically thin and thick media to synchrotron radiation, and ${\nu}_{\rm c}$is a critical frequency (see below).

An appropriate approximation for $\varepsilon_{\nu}$ in optically thin medium reads (Mahadevan et al. 1996):

$\displaystyle \varepsilon_{\nu} = 2.73\times 10^{-5} \frac{\rho \nu }{K_{2}(1/\... ...m e})}{\tilde{\cal I}}(\nu,B,\Theta) ~{\rm ergs~cm^{-3}~s^{-1}~ {\rm Hz}^{-1}},$     (7)

where K2 is the Bessel function of the second kind and $\tilde{\cal I} = \frac{4.05}{\zeta^{1/6}}(1 + \frac{0.4}{\zeta^{1/4}} + \frac{0.53}{\zeta^{1/2}}) {\rm e}^{-1.89 \zeta^{1/3}}.$
Here $\zeta = 2.38\times 10^{-7} ({\nu}/{B{\theta}^{2}_{\rm e}})$ and $\theta_{\rm e} = {k T_{\rm e}}/{m_{\rm e} c^2}.$
Below a certain critical frequency ${\nu}_{\rm c}$, the media become self-absorbing to synchrotron emission. In this case, $ \varepsilon_{\nu}^{\rm mod} \approx \varepsilon^{\rm BB}_{\nu} ={2\pi}\frac{\nu^2}{c^2}kT. $To find ${\nu}_{\rm c}$, we use the local non-linear Newton iteration procedure applied to the equation

\begin{displaymath}\int_{V}{\varepsilon_{\nu}{\rm d}V} = \int_{S}{\varepsilon^{\rm BB}_{\nu}} {\rm d}S. \end{displaymath} (8)

Having obtained ${\nu}_{\rm c}$, the switch on/off operator $\xi$ can then be constructed.

3 Method of solution

Let ${\cal L}E=0$ be the equivalent operator form of Eq. (4) in the continuous space $\Omega_{\rm C}$. $\cal{L} E$ consists of several terms, each of which requires a careful and different representation in the finite spherical discretization space $\Omega_{\rm h}$ which is defined as $[t_{1},t_{2},...,t_{N}]\otimes[r_{1},r_{2},...,r_{J}] \otimes[\theta_{1},\theta_{2},...,\theta_{\rm K}]\otimes[\nu_{1},\nu_{2},...,\nu_{\rm M}]$. ${[t],~[r],~[\theta]}$ and $[\nu]$ correspond to time, radius (spherical), latitude, and frequency intervals, respectively.

In most astrophysical problems the radiative effects occur on relatively short time scales compared to the hydro- or magneto-hydrodynamical ones. Therefore it is essential to treat these effects using unconditionally implicit numerical solvers. This requires however that all terms of Eq. (4) should be evaluated on the new time-level. Further, the spatial discretization used should guarantee that the resulting coefficient-matrix $A_{\rm r\theta\nu} = \partial{\cal L} E/\partial E$ in $\Omega_{\rm h}$ is strictly diagonally dominant. To fulfill this requirement, we perform the following procedures.
1) The term $\nabla \cdot V E_{\nu} $ which is an advection term in $\Omega_{\rm C}$ is transformed to $\Omega_{\rm h}$ using a second order up-winding discretization scheme.
2) $\nabla \cdot [\frac{\lambda_{\nu}}{\chi}\nabla E_{\nu}]$ is a second order diffusion term $\mapsto$ second order central-difference operator in $\Omega_{\rm h}$.
3) $\lambda_{\nu} (\nabla \cdot V) E_{\nu}$, $ \kappa \rho (B_{\nu}-E_{\nu})$ and $\varepsilon_{\nu}$ are local functions $\mapsto$ locally-evaluated functions at cell centers.
4) ${\cal K}_{\nu}$ constitutes of an advection and diffusion terms in the frequency space $\mapsto$ second order central-difference operators in $\Omega_{\rm h}$.
Combining the contributions from all the terms of Eq. (4), the local equation at each point (j,k,m) reads:

\begin{displaymath}{ \underline{S}^{\rm r} E^{\rm new}_{j-1,k,m} + \overline{S}^... ...w}_{j+1,k,m} + \underline{S}^{\theta} E^{\rm new}_{j,k-1,m} }\end{displaymath}


\begin{displaymath}{+ \overline{S}^{\theta} E^{\rm new}_{j,k+1,m} } {+ \underli... ...\rm new}_{j,k,m-1} + \overline{S}^{\nu} E^{\rm new}_{j,k,m+1}} \end{displaymath}


\begin{displaymath}{ + (D^{\rm r} + D^{\theta} + D^{\nu}) E^{\rm new}_{j,k,m} = RHS}, \end{displaymath} (9)

where $\underline{ S}^{\rm r} = \partial{\cal L} E_{\nu}/\partial E_{j-1,k,m}$, $\overline{ S}^{\rm r} = \partial{\cal L} E_{\nu}/\partial E_{j+1,k,m}$, ${\rm D^{\rm r} + D^\theta + D^\nu} = \partial{\cal L} E_{\nu}/\partial E_{j,k,m}$, $\underline{ S}^{\theta} = \partial{\cal L} E_{\nu}/\partial E_{j,k-1,m}$, $\overline{S}^{\theta} = \partial{\cal L} E_{\nu}/\partial E_{j,k+1,m}$, $\overline{ S}^{\nu} = \partial{\cal L} E_{\nu}/\partial E_{j,k,m+1}$ and $\underline{ S}^{\nu} = \partial{\cal L} E_{\nu}/\partial E_{j,k,m-1}$. Note that $\underline{ S}^{\rm r}, D^{\rm r},$ and $ \overline{ S}^{\rm r}$ respectively correspond to the sub-diagonal, diagonal and super-diagonal entries of the coefficient-matrix $A_{\rm r\theta\nu}$ in the radial direction. A similar description applies to the $\theta-$ and $\nu-$directions also. Also, note that the contribution from operators in other directions to the diagonal elements of the matrix is $ ( D^{\rm r} + D^{\theta} + D^{\nu}) \delta E_{j,k,m} $. These pronounce the diagonal-dominance of the matrix and might enable a stable inversion procedure. Thus, the equivalent form of Eq. (4) in $\Omega_{\rm h}$ at all grid points is the matrix equation: $ A_{\rm r\theta\nu} E^{\rm new} = E^{\rm old}$, or more generally Aq=b.

The important question at this stage is how to invert A efficiently.

Since the matrix $\,A\,$ is sparse (i.e., the majority of its entries are zeros), it can be easily be approximated by a matrix $\,{\tilde A}\,$ which, (a) is easy to construct and to invert, and (b) does not differ significantly from $\,A\,$ in the sense that  $\,{\tilde A}\,$ and $\,A\,$ share the essential spectral properties ("spectral equivalence'', see Hackbusch 1985 for further details). In the mathematical literature  $\,{\tilde A}\,$ is usually called a "pre-conditioner''. In our case $\tilde{A}$ is obtained by using a first order upwinding scheme for the advection terms. Thus, starting with some initial guess $\,q^0\,$, consistency with the original problem is maintained by adopting the so called "defect-correction'' iteration procedure (Hujeirat 2001). A defect correction iteration for computing $\,q^i(=q^{\rm new})\,$ from $\,q^{i-1}(=q^{\rm previous})\,$ is organized as follows:

1.
Compute the defect $\,d^{i} = b\!-\!A q^{i}\,$;
2.
Solve the defect equation $\,{\tilde A}\mu^{i} = d^{i}\,$;
3.
Correct $\,q^{i} = q^{i-1}\!+\!\omega_{i}\mu^i\,$ ( $\,\omega_{i}$ a damping factor);
4.
If $\,\Vert\mu^{i}\Vert/\Vert q^{i}\Vert \le \epsilon\,$, then STOP, otherwise GOTO 1.
Common examples for pre-conditioners $\,{\tilde A}\,$ are: More sophisticated pre-conditioners can be obtained by using "projection methods'', "Krylov-subspace iterative methods -KSIM'', "multigrid methods'', or "domain decomposition methods'' (see Hujeirat & Rannacher 2001 for further details and applications).

Given the properties of A, an optimal choice would be to use a direct solver. Unfortunately, almost all known direct solvers rely on Gauss elimination which undermines the advantageous of the sparse-property of A (they fill-in the zero-entries), requires an extremely large memory and prohibitive amount of algebraic operations ($\sim $N3, where N is the number of grid points). For example, using Gauss elimination to invert a matrix that corresponds to $200\times 100$ grid points in $r-\theta$ plane and 200 frequencies requires approximately 1019 algebraic operations. Therefore we can safely rule out direct methods from our present considerations.

On the other hand, iterative methods are most suited for sparse-matrices and their convergence rates increase the more the matrix A is diagonally dominant. The essential feature of efficient iterative methods and in particular KISMs is that the system matrix $\,A\,$ is not explicitly involved in the process but rather its multiplication with a vector. This is particularly useful for sparse matrices for which storage is no longer a critical problem as in direct methods (Saad 1996). However, the convergence rate of most iterative methods has been found to depend strongly on the proper choice of the pre-conditioner, and hence indirectly on on the mesh size. Thus, the much finer the mesh cells are, the slower is the convergence rate.

A promising method frequently used in CFD whose convergence rate is independent of mesh size is the so called multi-grid method (MGM).

The main idea underlying MGM is the fast reduction of high-frequency error components by "cheap'' relaxation ("smoothing'') on the fine meshes and the reduction of the remaining low-frequency error components by defect correction on coarser meshes ("coarse-grid correction''). Since most of the algebraic operations are actually done on coarser grids, the combined method requires significantly less work than iterating only on a fixed fine grid. For an introduction to MGM see Hackbusch (1985). Since the convergence rate of the MGM is independent of the grid size, this method is particularly suited for systems which result from very fine meshes. Usually, the performance is satisfactory when the equations to be solved are of elliptic or parabolic character. This is because here error propagation is extremely fast compared to the flow velocities. On the other hand, MGM has problems when the flow is advection-dominated. Here, errors that appear to slow the convergence rate on the fine grid may easily be advected by the flow and there may be no gain from the coarse-grid correction. Unfortunately, MGMs are not the optimal solver for Eq. (4), mainly because all second order operators degenerate into first order advective ones in optically thin regions where the radiative flux is being advected with the speed of light, as well as in regions where the values of the Compton-Y parameter is effectively small compared to unity.

Other possible approaches for sparse systems are the Alternating Direction Implicit (ADI) and the Approximate Factorization Methods (AFM). However, ADI is not appropriate for searching steady solution in three or more dimensions, as it is numerically unstable (Fletcher 1988). Alternatively, we have tried the AFM as a pre-conditioner. However, it turns out that the AFM converges slower than our favorite iterative method: "Black-White-Brown'' line Gauss-Seidel method (henceforth BWB-LGS, see Fig. 1 for illustration and Hirsch 1990 for further details). The latter method pronounces the diagonal dominance of $\tilde{A}$, and hence converges faster than AFM. It should be noted that the line Gauss-Seidel method in its classical form is not appropriate for vector and parallel machines, mainly because the vector-length is proportional to the number of unknowns in one direction. One possibility to overcome this barrier is to solve for all unknowns located on even-numbered grid points, followed by solving for those on odd-numbered grid points. The resulting vector-length in this case is proportional to the number of unknowns in the plane under consideration, and therefore enabling to perform the calculations on vector and parallel machines with high efficiency.

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig_1new.eps}\end{figure} Figure 1: The Black-White-Brown line Gauss-Seidel relaxation method. For each frequency $\nu _{\rm m}$ we solve for the unknowns in the plane $(R,\theta )$ in two sweeps: unknowns along odd-numbered lines are solved and followed by those on even-numbered lines. Similarly, this procedure is applied to the $(\theta ,\nu )$ and $(R,\nu )$ planes also.
Open with DEXTER

More specifically, in each plane we perform two sweeps: in the first sweep we consider the unknowns in the $r-\theta$ plane, i.e., we solve the system of equations:

\begin{displaymath}\underline{S}^{\rm r} \delta E^{\rm new}_{j-1,k,m} + (D^{\r... ...,m} + \overline{S}^{\rm r} \delta E^{\rm new}_{j+1,k,m} = RHS, \end{displaymath}

where $j=1\rightarrow J$ and k runs over odd numbers. In the second sweep, we solve:

\begin{displaymath}\underline{S}^{\rm r} \delta E^{\rm new}_{j-1,k,m} + (D^{\r... ...}_{j,k,m} + \overline{S}^{\rm r} \delta E^{\rm new}_{j+1,k,m}\end{displaymath} \begin{displaymath}= RHS + \underline{S}^{\theta} \delta E^{\rm new}_{j,k-1,m} + \overline{S}^{\theta} \delta E^{\rm new}_{j,k+1,m}, \end{displaymath}

where $j=1\rightarrow J$ and k here runs over even numbers. Therefore, we actually perform 6-inversion procedures per each time step. Here the 3-dimensional problem is being replaced by three one-dimensional problems that are solved iteratively to recover the solution of the original problem. The method is relatively efficient, as the overall number of arithmetic operations scales linearly with the number of grid points ($\sim $ $6\times 9\times N$).

4 Applications

As a reference configuration we consider a 2D accretion disk around a $10~{\cal M}_{\odot}$ black hole. The disk is threaded by an equipartition large scale magnetic field and embedded in a hot, tenuous and static corona, where $(\rho,{T}) \approx (10^{-6}\rho_{\rm disk},T_{\rm virial})$. The thickness of the disk is taken to be H = 0.1 r, and $(\rho,T)$ are set to adopt radial profiles similar to those in standard accretion disks, i.e., $\sim(r^{-15/8},r^{-3/4})$. Inside the disk-region the matter is isothermal in the vertical direction, yielding therefore an exponential decrease of the density with latitude (Frank et al. 1991).

The domain of calculation is restricted to $[r_{\rm in},r_{\rm out}]\times[0, \pi/2]\times[\nu_1,\nu_2]$, where $r_{\rm in} = 2.75 R_{\rm Sch}$ and $ r_{\rm out}= 50 R_{\rm in},$ and $R_{\rm Sch}$ is the Schwarzchild radius[*]. $\nu_1 = 10^{-2}~kT(r_{\rm out})/h$ and $\nu_2 = 2.41\times 10^{19}~ {\rm Hz}$. The domain of calculation is sub-divided respectively into $200\times 50 \times 200$ strongly stretched cells in the r, $\theta$ and $\nu$ directions. Grid stretching in the present case is essential in order to capture the fine structure of the disk, as well as to cover a considerable dynamical range in the radial direction and in the frequency space. For example, using 200 uniformally distributed grid points to cover the interval $[1\le r\le 100]$, yields a grid spacing that is five orders of magnitude larger than the pressure scale height of a typical cold disk. Hence, a uniform grid distribution is not appropriate for cold accretion flows. The same applies in the frequency space, if soft and hard X-rays are to be considered. Therefore, we use a grid spacing in the innermost region that is five orders of magnitude smaller than in the outermost region. Such an extremely stretched mesh is not useful for time-explicit schemes, as they enhance their numerical instability even more.

We note that although the time derivative of the radiative density (first term of Eq. (4)) should be neglected if time-independent solutions are of interest only, considering this term while gradually increasing the time-step size with each iteration, has been verified to enhance the global convergence to steady state solutions. This pseudo time-stepping procedure has been verified to be useful in the 2D case, especially if ADI is used as a direct solver for the Possion equation when modeling self-gravitating flows or for the radiation transfer equation under optically thick conditions (Yorke et al. 1993). In Fig. 3, the evolution of the maximum time-independent residual as function of the iteration number is shown (see Fig. 2 for the initial distributions of the density and the magnetic field). The numerical scheme presented here appears, as it is obvious from Fig. 3 and from other model test calculations (e.g., using different density and temperature distributions), to be stable and capable to provide steady state solutions for a large variety of parameters. In particular, it can be used to constrain and find the $\rho-,T-$ and B-distributions appropriate for displaying spectra that are consistent with observations.

When applying this scheme to accretion disks around black holes, the method is capable to reproduce the Rayleigh-Jeans branch of the spectra in the low frequency regime (see Fig. 4). As the disk here is embedded in a hot and tenuous medium, the radiative flux adopts a flatter profile than ${ F_{\nu} \sim {\nu}^{1/3}}.$ This is a consequence of using the modified source function $S_{\nu}$ which, unlike the classical Planck function, does not vanish in the corona-disk transition region but induces a considerable smoothing of the radiative flux.

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig_2.ps}\end{figure} Figure 2: Uniformaly distributed isolines of the density (solid lines) and magnetic field strength (dashed lines) threading an accretion disk with ${\,\dot{\cal M}}=0.1 {\,\dot{\cal M}}_{\rm Edd}$ around a $ 10{\,{\cal M}_{\odot }}$ black hole. The disk extends from $R_{\rm in}= 1 $ to $R_{\rm out} = 50$, where R is given in unites of $2.75~R_{\rm Sch}$. This configuration is used to calculate disk spectra in the following models.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig_3.ps}\end{figure} Figure 3: The evolution of the residual in the maximum norm versus the iteration number. The adopted density and temperature profiles correspond to standard accretion disks (see Fig. 2).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig_4.ps}\end{figure} Figure 4: The normalized radiative flux of optically thick disks versus frequencies viewed from three different angles: $\Theta_1=10^0,~~\Theta_2=30^0$ and $\Theta _3=80^0$. In this special case, we switch off Comptonization and synchrotron emission and allow a constant accretion rate ${\,\dot{\cal M}}=0.1 {\,\dot{\cal M}}_{\rm Edd}$ to enter the domain through the outer boundary. In these models the disks are set to be embedded in a hot and tenuous medium The solid lines in this plot correspond to standard accretion disks, while the dotted lines correspond to truncated disks at $r=11~R_{\rm Sch}.$
Open with DEXTER

This feature becomes more pronounced when truncating the disk at ${r=11 R_{\rm Sch}}$, indicating therefore the importance of the corona-disk interaction region in the vicinity of the black hole in supplying energy to the outer regions.

In order to study the effect of Comptonization, we have carried out three model calculations assuming an accretion rate ${\,\dot{\cal M}}=0.025 {\,\dot{\cal M}}_{\rm Edd}$ across the outer boundary. The truncation radius and the 2D distribution of density and temperature are taken to be similar to our previous disk-calculations (Hujeirat & Camenzind 2000).

In the first case, only soft-photons from the cold disk are considered. The corresponding spectrum is shown in Fig. 5 (the profile with $\beta=0$). Apparently, a small fraction of the soft photons originating from the disk are Compton up-scattered by the hot electrons surrounding the disk. This process becomes more efficient when the magnetic field strength is gradually increased with respect to the gas pressure (see the profiles with $\beta=0.1 $ and $\beta =0.5$ in Fig. 5)[*]. In addition to Comptonization of soft photons from the cold disk, Compton up-scattering of synchrotron radiation has a significant effect on the photon energy-enhancement. It also pronounces the power-law dependence of the flux on the frequencies.

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{fig_5.ps}\end{figure} Figure 5: The normalized radiative flux of optically thin disks versus frequencies. The profiles here correspond to a truncated disk at $r=11~R_{\rm Sch}$ embedded in a hot and tenuous medium. The solid line corresponds to data obtained with Comptonization of soft photons from the cold disk without magnetic fields. The profiles with $\beta =P_{\rm mag}/P_{\rm gas}= 0.1$ and $\beta =0.5$ corresponds to data in which synchrotron emission has been included.
Open with DEXTER

5 Summary

In this paper we have presented a numerical tool for solving the angle-averaged time-dependent radiation transfer equation, taking into account the Kompaneets operator to model up-scattering of soft photons by hot electrons. The advection of the radiative density has been considered to enable coupling with our available 3D-axisymmetric implicit MHD numerical solver. The radiative solver presented here is capable to handle Comptonization of soft photons originating either from a cold accretion disk or from synchrotron radiation of charged particles in a magnetized relativistic plasma.

As a preliminary test, we have applied the numerical solver to model the spectra of an accretion disk surrounding a Schwarzchild BH. The results obtained indicate that Comptonization of synchrotron radiation in the corona and in the corona-disk interaction regions is likely to be the main mechanism underlying the emission of the hard X-rays observed in LMXBs. This is in agreement with the recent astronomical data of the X-ray nova XTE J1118+480 (Esin et al. 2000).

More consistent calculations should allow a two-temperature description of the flow, and solving the internal energy equation fully coupled to Eq. (4). This together with more specific applications is the subject of the forthcoming paper.

Acknowledgements
AH's research has been supported by the Deutsche Forschungsgemeinschaft (DFG) through the Schwerpunktprogramm Physik der Sternenentstehung. Computations were performed on facilities provided by Interdisciplinary Center for Scientific Computing (IWR) at the University of Heidelberg.

References

 

Copyright ESO 2002