A&A 386, 757-762 (2002)
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
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
The majority of the observed Low Mass X-ray Binaries (LMXBs) radiate a significant fraction of their total luminosity in X- and 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.
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:
To make the problem tractable, the following approximations have been employed:
The above two different behaviour of the operator can be combined as follows:
in Eq. (4) corresponds to the modified synchrotron
emission of photons by relativistic electrons gyrating around magnetic field lines, which
An appropriate approximation for
in optically thin medium reads
(Mahadevan et al. 1996):
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
is strictly diagonally dominant. To fulfill this requirement, we perform the
1) The term which is an advection term in is transformed to using a second order up-winding discretization scheme.
2) is a second order diffusion term second order central-difference operator in .
3) , and are local functions locally-evaluated functions at cell centers.
4) constitutes of an advection and diffusion terms in the frequency space second order central-difference operators in .
Combining the contributions from all the terms of Eq. (4), the local equation at each point (j,k,m) reads:
The important question at this stage is how to invert A efficiently.
Since the matrix is sparse (i.e., the majority of its entries are zeros), it can be easily be approximated by a matrix which, (a) is easy to construct and to invert, and (b) does not differ significantly from in the sense that and share the essential spectral properties ("spectral equivalence'', see Hackbusch 1985 for further details). In the mathematical literature is usually called a "pre-conditioner''. In our case is obtained by using a first order upwinding scheme for the advection terms. Thus, starting with some initial guess , consistency with the original problem is maintained by adopting the so called "defect-correction'' iteration procedure (Hujeirat 2001). A defect correction iteration for computing from is organized as follows:
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 (N3, where N is the number of grid points). For example, using Gauss elimination to invert a matrix that corresponds to grid points in 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 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 ,
and hence converges faster
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.
|Figure 1: The Black-White-Brown line Gauss-Seidel relaxation method. For each frequency we solve for the unknowns in the plane 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 and planes also.|
|Open with DEXTER|
More specifically, in each plane we perform two sweeps: in the first sweep we consider
the unknowns in the
plane, i.e., we solve the system of equations:
As a reference configuration we consider a 2D accretion disk around a black hole. The disk is threaded by an equipartition large scale magnetic field and embedded in a hot, tenuous and static corona, where . The thickness of the disk is taken to be H = 0.1 r, and are set to adopt radial profiles similar to those in standard accretion disks, i.e., . 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 , where and and is the Schwarzchild radius. and . The domain of calculation is sub-divided respectively into strongly stretched cells in the r, and 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 , 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 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
This is a consequence of using the
modified source function
which, unlike the classical Planck function, does
not vanish in the corona-disk transition region but induces a considerable smoothing
of the radiative flux.
|Figure 2: Uniformaly distributed isolines of the density (solid lines) and magnetic field strength (dashed lines) threading an accretion disk with around a black hole. The disk extends from to , where R is given in unites of . This configuration is used to calculate disk spectra in the following models.|
|Open with DEXTER|
|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|
|Figure 4: The normalized radiative flux of optically thick disks versus frequencies viewed from three different angles: and . In this special case, we switch off Comptonization and synchrotron emission and allow a constant accretion rate 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|
|Open with DEXTER|
This feature becomes more pronounced when truncating the disk at , 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 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 ). 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
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
|Figure 5: The normalized radiative flux of optically thin disks versus frequencies. The profiles here correspond to a truncated disk at 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 and corresponds to data in which synchrotron emission has been included.|
|Open with DEXTER|
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.
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.