A&A 370, 707-714 (2001)
DOI: 10.1051/0004-6361:20010310

Analytical solution of the radiative transfer equation in the two-stream approximation

N. V. Kryzhevoi^1,2 - G. V. Efimov³ - R. Wehrse^1,2

1 - Institut f. Theoretische Astrophysik, Universität Heidelberg, Tiergartenstr. 15, 69121 Heidelberg, Germany
2 - Interdisziplinäres Zentrum f. Wissenschaftliches Rechnen, Universität Heidelberg, Im Neuenheimer Feld 368, 69120 Heidelberg, Germany
3 - Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, 141980 Dubna, Russia

Received 15 November 2000 / Accepted 27 February 2001

Abstract
The analytical solution of the plane-parallel radiative transfer equation is obtained in the two-stream approximation for a large class of continuous distributions of the de-excitation coefficient $\varepsilon$ (constant, linear, parabolic, with spikes etc.). We present also the method of the discrete space theory for obtaining solutions of the transfer equation in the media with strong density inhomogeneities. These sets of the analytical solutions can be used for the solution of the inverse problem. The deduction of the internal distribution of $\varepsilon$ from observational data is facilitated in the case of isothermal media, since the characteristic behavior of the solution refers to the certain behavior of $\varepsilon$. As an example, we find the corresponding parameters of the constant and linear distributions of $\varepsilon$ precisely.

Key words: radiative transfer - methods: analytical - methods: statistical

1 Introduction

From the time of the first derivation of the transfer equation about a century ago much work has been dedicated to its specific forms and their solution. Unfortunately, among the large number of published algorithms there are still only very few methods that provide exact solutions and they refer to rather particular situations. So, to our knowledge all analytical solutions have been obtained under assumption of a constant de-excitation coefficient $\varepsilon$ or for $\varepsilon$ approximated by a piecewise constant function. However, there are many cases, e.g. the transmission of light through the earth's atmosphere, the emission of radiation by nonisotropic high-temperature gas streams, where the properties of the medium may significantly vary with position. The difficulties associated with a complex behavior of $\varepsilon$ can be overcome by applying numerical methods. But in spite of the good success they give little insight into the general behavior of the solution of the transfer equation and sometimes they are extremely CPU-time and memory consuming. In addition the numerical methods are not efficient when one deals with media with many strong density inhomogeneities. Note that such media are ubiquitous as high resolution observations of the solar atmosphere, nova and supernova remnants, accretion disks etc. have shown. For the necessity to include the strong spatial variations in reliable models see e.g. Gu et al. (1997). Unfortunately, in contrast to the Navier-Stokes equations, for the radiative transfer equation there is not yet a general homogenization scheme available that can deal with such situations, and exististing algorithms (cf. Gierens et al. 1986; Lindsey & Jefferies 1990; Nikoghossian et al. 1997) have a very limited range of application.

In the present paper we have obtained the analytical solution of the radiative transfer equation for the different distributions of $\varepsilon$ in the plane-parallel media. The description of the problem, the basic equations and the simplifying assumptions are given in Sect. 2. In Sect. 3 we present the general solution which requires for the given $\varepsilon (\tau )$ only the knowledge of linearly independent solutions of a second order homogeneous differential Eq. (2). In Sect. 4 we give several examples of $\varepsilon (\tau )$ and corresponding linearly independent solutions. In the same section we offer also a method of the solution of the transfer equation which can be useful for $\varepsilon (\tau )$ with spikes. Having obtained the analytical solutions for the different behaviors of $\varepsilon (\tau )$ the solution of the inverse problem becomes possible. Section 5 is dedicated to the possibility of the diagnostic of $\varepsilon (\tau )$from observational data. In Sect. 6 the case of the medium with a stochastic distribution of $\varepsilon$ is considered. The solution is obtained by the method of the discrete space theory. Section 7 contains a summary.

   

Table 1: The different kinds of $\varepsilon (\tau )$ and corresponding linearly independent solutions of the homogeneous Eq. (2). The divergence at $\tau =0$ in the last three examples can be removed by an appropriate coordinate shift

$\varepsilon (\tau )$	$Y_1(\tau)$	$Y_2(\tau)$
b2	$\cosh(b\tau)$	$\frac1b \sinh(b\tau)$
$a\tau+b$	${\rm Ai}\left(\frac{a\tau+b}{\vert a\vert^{2/3}}\right)$	$\pi\frac{\vert a\vert^{2/3}}{a}{\rm Bi}\left(\frac{a\tau+b}{\vert a\vert^{2/3}}\right)$
$a^2\tau^2+b$	$a^{1/4}{\rm e}^{-\frac{a\tau^2}{2}} F(\frac{a+b}{4a},\frac12; a\tau^2)$	$a^{-1/4}{\rm e}^{-\frac{a\tau^2}{2}}\tau F(\frac{3a+b}{4a},\frac32; a\tau^2)$
$\frac{b^2m^2-\frac14}{\tau^2}+\frac{b^2a^2\tau^{2b}}{\tau^2}$	$\sqrt{\tau}\,K_m(a\tau^b)$	$\frac1b\sqrt{\tau}\,I_m(a\tau^b)$
$\frac14-\frac{b-2a}{2\tau}-\frac{b(2-b)}{4\tau^2}$	$\tau^{\frac b2}{\rm e}^{-\frac{\tau}2}\, F(a,b; \tau)$	$\frac1{1-b}\tau^{\frac {1-b}2}{\rm e}^{-\frac{\tau}2}\, F(a-b+1,2-b; \tau)$
$\frac{b^2m^2-\frac14}{\tau^2}-\frac{kab^2\tau^b}{\tau^2}+ \frac{b^2a^2\tau^{2b}}{4\tau^2}$	$\tau^{\frac{1-b}2} M_{k,-m}(a\tau^b)$	$\frac1{2mab}\tau^{\frac{1-b}2} M_{k,m}(a\tau^b)$

${\rm Ai}(\tau)$, ${\rm Bi}(\tau)$ - the Airy functions,

$I_m(\tau)$, $K_m(\tau)$ - the modified Bessel functions,

$F(a,b; \tau)$ - the Kummer confluent hypergeometric function,

$M_{k,m}(\tau)$ - the Whittaker function, $2m\neq\pm1,\pm2,\pm3...$

2 The radiative transfer equation

The radiative transfer equation for the specific intensity $I(\tau)$ in plane-parallel media reads

 

$$\mu{{\rm d}\over {\rm d}\tau}I(\tau) -I(\tau)+(1-\varepsilon(\tau)) J(\tau)+\varepsilon (\tau)B(\tau),$$ = (1)

where $\varepsilon (\tau )$ is the de-excitation coefficient, $B(\tau)$ is the Planck function, $J(\tau)$ is the mean intensity, $\mu=\cos\theta$, $\theta$ is an angle between the normal and a ray of radiation.

Let us suppose that the radiation field can be characterized by a discrete number of directed streams ("discrete ordinates'') to mimic the true variation of intensity with angle. To simplify the problem we consider only two rays in the opposite directions $\mu=\pm1$ that gives already reasonably accurate results. The Feautrier technique (Mihalas 1978) allows us to transform Eq. (1) into a second order differential equation

 

$$\left(-{{\rm d}^2\over {\rm d}\tau^2}+\varepsilon(\tau)\right) J(\tau)=\varepsilon(\tau) B(\tau),$$ (2)

where the mean intensity $J(\tau)$ and the flux $F(\tau)$ are expressed through the specific intensities in the positive and negative directions by the following way

$$J(\tau)=\frac12(I^+(\tau)+I^-(\tau)),\quad F(\tau)=\frac12(I^+(\tau)-I^-(\tau)).$$

We look for the solution in a slab whose mid plane is a symmetry plane. The optical depth is measured away from the symmetry plane. It equals $\Delta$ and $-\Delta$ at the upper and lower boundary, respectively. Due to the symmetry it is sufficient to obtain the solution e.g. for the upper part $0\leq\tau\leq\Delta$. We assume no incident radiation at the boundary $I^-(\Delta)=0$. At the symmetry plane the radiation field satisfies the reflection condition I⁺(0)=I^-(0). Since in the two-stream approximation $F(\tau)=-J'(\tau)$, the boundary conditions can be written as

 

$$J(\Delta)+J'(\Delta)=0,$$ (3)

where differentiation with respect to $\tau$ is denoted by prime.

	Figure 1: The runs of the mean intensity in the optically thin ( left) and optically thick ( right) isothermal, $B(\tau )=1$, slabs with the different distributions of $\varepsilon (\tau )$ (see insert)
Open with DEXTER

	Figure 2: The mean intensity as the function of the optical depth and $\varepsilon (\tau )$ whose internal distribution can be approximated by a resonance curve. The dependence on the shape and the position of the resonance is shown
Open with DEXTER

3 The general solution

Let $\{Y_1,(\tau), Y_2(\tau)\}$ be linearly independent solutions of the homogeneous Eq. (2) satisfying the normalization condition of the Wronskian

$$\begin{eqnarray*}W&=&Y_1(\tau)Y_2'(\tau)-Y_2(\tau)Y_1'(\tau)=1. \end{eqnarray*}$$

Then the formal solution of Eq. (2) can be written as the following

  

$$J(\tau) C_1\,Y_1(\tau)+C_2\,Y_2(\tau)+\int_{0}^{\tau}\!\!\left(Y_1(\tau)Y_2(\tau')-Y_1(\tau')Y_2(\tau)\right) \,{\cal B(\tau')}\,{\rm d}\tau',$$ = (4)

$$J'(\tau) C_1\,Y_1'(\tau)+C_2\,Y_2'(\tau)+\int_{0}^{\tau}\!\!\left(Y_1'(\tau)Y_2(\tau')-Y_1(\tau')Y_2'(\tau)\right) \,{\cal B(\tau')}\,{\rm d}\tau',$$ = (5)

where ${\cal B}(\tau)=\varepsilon(\tau)\,B(\tau)$. C₁ and C₂ are arbitrary constants whose values are defined by the boundary conditions (3).

After some transformations (see Appendix A) one can get

 

$$J(\tau) \frac{({\vec Z}(\Delta)\sigma{\vec Y}(\tau))} {({\vec Z}(\Delta)\... ...elta}\! ({\vec Z}(\Delta)\sigma{\vec Y}(\tau'))\,{\cal B(\tau')}\,{\rm d}\tau',$$ = (6)

where we use vectors

$$\begin{eqnarray*}{\vec Y}(\tau)=\left(\begin{array}{c} Y_1(\tau)\\ Y_2(\tau) \e... ...\prime}(\tau)\\ Y_2(\tau)+Y_2^{\prime}(\tau) \end{array}\right) \end{eqnarray*}$$

and matrix $\sigma= \left(\begin{array}{cc} 0&1\\ -1&0 \end{array}\right)$. Since $({\vec Z}(\Delta)\sigma{\vec Y}(\Delta))=-1$, the mean intensity at the boundary $\tau=\Delta$ becomes

 

$$J(\Delta) \int_{0}^{\Delta}\!\frac {({\vec Y}(\tau')\sigma{\vec Y}'(0))}{({\vec Z}(\Delta)\sigma{\vec Y}'(0))} \,{\cal B(\tau')}\,{\rm d}\tau'.$$ = (7)

As one can see the exact solution of the radiative transfer equation for the given run of $\varepsilon (\tau )$ requires only the knowledge of the linearly independent solutions of the homogeneous Eq. (2).

4 Examples

Some examples of the continuous distributions of $\varepsilon (\tau )$ and corresponding solutions $Y_1(\tau)$ and $Y_2(\tau)$ (taken from Abramowitz & Stegun 1972; Kamke 1965) have been collected in Table 1. The free parameters must be chosen in such a way to satisfy the condition of the location of $\varepsilon (\tau )$ in the interval between 0 to 1. In spite of the small variation range of $\varepsilon (\tau )$, the solutions obtained for different $\varepsilon (\tau )$ may have significant difference, especially in optically thick media. So, in Fig. 1 the solutions of Eq. (2) with constant and linear $\varepsilon (\tau )$ are shown. In optically thin isothermal media (left part) the difference does not exceed 10%. However, it becomes larger with the increasing of the total optical thickness and in the optically thick media can reach at some points 50% (right part).

Although the functions presented in Table 1 are suitable for the approximation of a large variety of internal distributions of $\varepsilon (\tau )$, they cannot be applied for the description of media with strong density condensations. Furthermore, the solution of the homogeneous Eq. (2) can hardly be found directly with $\varepsilon (\tau )$ approximated by a function with spikes. To avoid these difficulties we suggest the following procedure: if $\varepsilon (\tau )$ can be represented as

 

$$\varepsilon(\tau) -A(\tau){{\rm d}\over {\rm d}\tau} \left({A'(\tau)\over A^2(\tau)}\right)+ \varepsilon^{\rm T}(\tau)A^4(\tau),$$ = (8)

then the corresponding solutions of the homogeneous Eq. (2) can be expressed through already known solutions in the following way

$$Y_i(\tau) {1\over A(\tau)}Y_i^{\rm T}(\phi(\tau)), ~~~~~~~~(i=1,2)$$ = (9)

where functions $\varepsilon^{\rm T}(\tau)$, $Y_1^{\rm T}(\tau)$ and $Y_2^{\rm T}(\tau)$ are taken from Table 1, $A(\tau)$ is known function and $\phi(\tau)$ can be found from differential equation $\phi'(\tau)=A^2(\tau)$.

For example, the choice of

$$\begin{eqnarray*}A(\tau)&=&0.9+\frac{0.3}{0.3+(\tau-8)^2},\quad \varepsilon^{\rm T} =0.047 \end{eqnarray*}$$

leads to the dotted curves in the right part of Fig. 2. The altering of the parameters results in the other curves.

	Figure 3: The mean intensity at the boundary as a function of the optical thickness $\Delta$. The curves correspond to the different runs of $\varepsilon (\tau )$ (see insert): $\varepsilon (\tau )=\tau /\Delta$ ( dashed), $\varepsilon (\tau )=0.6(\tau /\Delta )^2+0.3$ ( dashed-dotted), $\varepsilon (\tau )=-\tau /\Delta +1$ ( dotted) and $\varepsilon (\tau )=0.5$ ( solid)
Open with DEXTER

5 Diagnostics of $\vec{\varepsilon(\tau)}$

The prediction of the internal structure of a medium from observational data is one of the most important task in astrophysics. The only observable quantity of the problem is the emergent intensity that is a function of wavelength $\lambda$. The total optical thickness of a layer also depends on $\lambda$. The knowledge of these functions allows us to plot $J(\lambda(\Delta))=J(\Delta)$ and therefore makes the prediction of $\varepsilon (\tau )$ possible. In the general case when the solution depends both on the temperature and on $\varepsilon (\tau )$ the diagnostic of $\varepsilon (\tau )$ is hardly possible. However, in the isothermal media the features of the solutions associate only with the definite behavior of $\varepsilon (\tau )$ and therefore the derivation of the corresponding parameters of such behavior seems not to be so hopeless. In order to confirm that we consider a slab with B=1, constant and linear $\varepsilon (\tau )$. In these cases the integration in (7) gives

	Figure 4: The mean intensity at the boundary of a medium with a strong density condensation
Open with DEXTER

 

$$J(\Delta) \frac{\sqrt{\varepsilon}\tanh(\sqrt{\varepsilon}\Delta)} {1+\sqrt... ...}\tanh(\sqrt{\varepsilon}\Delta)} ~~ {\rm for } ~\varepsilon(\tau)=\varepsilon,$$ = (10)

$$J(\Delta) \frac{{\cal A}}{1+{\cal A}} ~~~~~~~~~ {\rm for } ~ \varepsilon(\tau)=\left(\frac{a}{\Delta}\right)\tau+b,$$ = (11)

where

$$\begin{eqnarray*}{\cal A}&=&\frac{a\,\Delta^{1/3}}{\vert a\vert^{2/3}} \left\{{\... ...left(\frac{a+b}{(\vert a\vert/\Delta)^{2/3}}\right)\right\}\cdot \end{eqnarray*}$$

These expressions can be used now for the fitting of observational data and in the case of a good fit the parameters of $\varepsilon (\tau )$ can be derived without much effort.

The derivation of the corresponding parameters can be done much easier if we take into account the behavior of these curves at the limit of large and small $\Delta$. At the limit of small $\Delta$ these functions are proportional to $\Delta$ whereas at large $\Delta$ they saturate (Fig. 3) in accordance with the following

 

$${J(\Delta)= \displaystyle\frac{\sqrt{\varepsilon}}{1+\sqrt{\varepsilon}} \atop J(\Delta)=\displaystyle \frac{\sqrt{a+b}}{1+\sqrt{a+b}}} {\rm for} \quad \Delta\gg1,$$ (12)

$${J(\Delta)=\varepsilon\Delta-\varepsilon^2\Delta^2 \atop J(\Delta)=\bar{\varepsilon}\Delta-\bar{\varepsilon}^2\Delta^2} {\rm for} \quad \Delta\ll1,$$ (13)

where $\bar{\varepsilon}=(a+2b)/2$.
We want to stress, however, that some features such as spikes which occur in intermediate points may be missed during the determination of the internal behavior of $\varepsilon (\tau )$ by means of this method (Fig. 4).

As mentioned above the presence of some peculiarities can point to the definite behavior of $\varepsilon (\tau )$ and, thus, simplify its diagnostic. For example, the maximum of the function $J(\Delta )$ may indicate to the linearly decreasing $\varepsilon (\tau )$ (Fig. 3), although other distributions of the de-excitation coefficient may also result in such feature.

6 Stochastic $\mathsf{\varepsilon}$

In order to solve the radiative transfer equation in the media with many strong density inhomogeneities we use the method of the discrete space theory. For this purpose we divide the slab into N layers

$$[0,\Delta]=\bigcup\limits_{j=1}^N[(j-1)\delta,j\delta], \quad \delta={\Delta\over N}\cdot$$ (14)

The Planck function and the de-excitation coefficient are assumed to be constant in each layers but their values differ from one layer to an other

$${\varepsilon(\tau)=\varepsilon_j={\rm const} \atop B(\tau)=B_j={\rm const}} {\rm for} \quad \tau\in[(j-1)\delta,j\delta].$$

In addition, $\varepsilon_j$ is a function of a random number r_j whose values may be independent as well as obey correlations from layer to layer.

We do not use the formalism of Peraiah (1984) or that of Schmidt & Wehrse (1987) based on the interaction principle which relates the incident and emergent intensities in a layer. Instead of these we propose another method which relates the mean intensity and the flux at one boundary of a layer with the same functions at another one. So, a vector

$${\vec J}_j(\tau) \left(\begin{array}{c} J_j(\tau) \\ J'_j(\tau) \end{array}\right)$$ = (15)

representing corresponding values in the jth cell can be represented in terms of that in sell (j-1) by the following equation (see Appendix B)

 

$${\vec J}_j {\vec U}_j{\vec J}_{j-1}-{\vec K}_jB_j,$$ = (16)

where

$$\begin{eqnarray*}{\vec U}_j&=&\left(\begin{array}{cc} \cosh(\omega_j\delta) & {... ...\delta) \end{array}\right), \quad \omega_j=\sqrt{\varepsilon_j}. \end{eqnarray*}$$

B_j and $\varepsilon_j$ denote the value of the Planck function and the de-excitation coefficient at the upper boundary of each layer.

A successive application of Eq. (16) - with the corresponding boundary conditions - allows us to study the evolution of the mean intensity in the medium (see Appendix B)

$$J_j={({\vec e}{\vec B}(N,j+1)) ({\vec f}{\vec U}(j,1){\vec f... ... ({\vec f}{\vec W}(j,1))\over({\vec e}{\vec U}(N,1){\vec f})},$$ (17)

as well as to compute the value of J at the boundary

$$J_N={({\vec f}{\vec W}(N,1))\over({\vec e}{\vec U}(N,1) {\vec f})},$$ (18)

where

$$\begin{eqnarray*}{\vec U}(j,i)&=&{\vec U}_j{\vec U}_{j-1}...{\vec U}_i, ~~~~~~(i... ...p(1,1){\vec W}_2)+...+B_{j-1}\,({\vec U}^\top(j-1,1){\vec W}_j), \end{eqnarray*}$$

with

 

$${\vec W}_j {\vec U}_j^\top\sigma{\vec K}_j =\left(\begin{array}{c} \omega_j\sinh(\omega_j\delta) \\ \cosh(\omega_j\delta)-1 \end{array}\right)$$ =

and

$$\begin{eqnarray*}{\vec f}&=&\left(\begin{array}{r} 1 \\ 0 \end{array}\right),~... ...~~ {\vec g}=\left(\begin{array}{r} 1 \\ -1 \end{array}\right). \end{eqnarray*}$$

The representations given for the matrix ${\vec U}_j$ and the vectors ${\vec K}_j$and ${\vec W}_j$ are unfortunately not well suited for numerical calculations, since the hyperbolic functions involved lead to machine overflows for large $\delta$. In order to overcome this problem we extract factors $X=\cosh(\omega_j\delta)$ from the expressions. It can be shown that then the X terms in the numerator and denominator of expressions (17)-(18) cancel, i.e we can use the following formulae

$${\vec U}_j \left(\begin{array}{cc} 1 &{1\over\omega_j}\tanh(\omega_j\delta) \\ \omega_j\tanh(\omega_j\delta) &1 \end{array}\right),$$ = (19)

$${\vec K}_j \left(\begin{array}{c} 1-{\rm sech}(\omega_j\delta) \\ \omega_j\... ...ega_j\tanh(\omega_j\delta) \\ 1-{\rm sech}(\omega_j\delta) \end{array}\right).$$ =

One is left with equations that involve tanh and sech functions only. If it is necessary, the ${\rm sech}(x)$ for large arguments can be approximated by $2{\rm e}^{-x}$. In this way we have obtained convenient expressions that are well suited for all optical depths.

In our example we divide the slab into 200 layers. We require that values of $\varepsilon_j$ lie in interval [0,1] and probability of appearance of small $\varepsilon_j$ be higher. As an example of such $\varepsilon_j$ we take the following

$$\varepsilon_j=10^{-4r_j},$$ (20)

where r_j are random numbers from the interval [0,1].

One realization of $\varepsilon$ is shown in Fig. 5. In Fig. 6 one can see the set of solutions for 50 realizations of $\varepsilon$ as well as the solution for $\bar{\varepsilon}$ in a medium with B=1. Their statistical distributions at different $\Delta$ are shown in Fig. 7.

	Figure 5: An example of the stochastic distribution of $\varepsilon$
Open with DEXTER

	Figure 6: Solutions obtained for 50 different realizations of $\varepsilon$. The bold curve was obtained for $\bar{\varepsilon}$
Open with DEXTER

Figure 7: Statistical distribution of $J(\Delta )$ from Fig. 6. Top left: The mean value is $\langle J\rangle=1.04~10^{-5}$, the standard deviation $\sigma=0.14~10^{-5}$, the value of the bold curve (see Fig. 6) at this point is $J_{\bar{\varepsilon}}=1.09~10^{-5}$. Top right: $\langle J\rangle=1.03~10^{-2}$, $\sigma=0.14~10^{-2}$, $J_{\bar{\varepsilon}}=1.08~10^{-5}$. Bottom left: $\langle J\rangle=0.24$, $\sigma=1.65~10^{-2}$. $J_{\bar{\varepsilon}}=0.247$. Bottom right: $\langle J\rangle=0.22$, $\sigma =0.052$, $J_{\bar{\varepsilon}}=0.247$

Open with DEXTER

Since the number of layers does not change the width of each layer becomes larger with the increasing of $\Delta$ and $\varepsilon$takes, thus, a block structure. This seems to be a reason of large scattering at large $\Delta$. In particular, the dependence of the variation range on the number of layers shown in Fig. 8 confirms this assumption.

7 Summary

The analytical solutions of the plane-parallel radiative transfer equation have been obtained for the large variety of the internal distributions of $\varepsilon$. We proposed also a method for obtaining the solution for $\varepsilon$ with spikes. This allows us to deal with media with a small number of the density inhomogeneities.

However, in very inhomogeneous media whose properties can only be treated statistically this method no longer works. Therefore we had to use another technique, namely, the method of the discrete space theory. The solution obtained is written as a sequence of the products of $2\times2$-matrices and two-components vectors that is very easy to be implemented. It is no longer necessary to solve the system of 2(N-1) linear equations (N - number of layers) (Wehrse 1981) or to use the method of the forward elimination and back substitution (Peraiah 1984) for the determination of the internal distribution of J. The problems related to the finding of the corresponding inverse matrix and the keeping of many coefficients are thus avoided, so that the numerical calculations are speeded up.

The presence of the analytical solutions enables us to solve the inverse problem. The more accurate diagnostic of $\varepsilon$ can be done in the isothermal media, since there the features of the solutions refer to the definite behavior of $\varepsilon$ only. By using the characteristic behaviors of the solutions in the limit of the large and small $\Delta$the exact derivation of the corresponding parameters is possible.

As mentioned in the Introduction, there are many classes of objects for which a spectral analysis requires the consideration of complicated depth dependencies of $\varepsilon$ (including stochastic ones) and their consequences for the system parameters. We plan to apply the algorithm developed here to the radiation fields of accretion disks.

	Figure 8: The variation range vs. the number of layers taken at $\Delta =50$
Open with DEXTER

Acknowledgements

We gratefully acknowledge support by the Deutsche Forschungsgemeinschaft (Sonderforschungsbereich 359/C2) and the Graduiertenkolleg at the Interdisciplinary Center for Scientific Computing at the University of Heidelberg.

Appendix A:

The solutions of Eq. (2) are given by (4). Using the boundary conditions (3) we get the integration constants C₁ and C₂ in following form

$$-Y_2'(0)\int_{0}^{\Delta}\! \frac{({\vec Z}(\Delta)\sigma{\vec Y}(\tau'))}{({\vec Z}(\Delta)\sigma{\vec Y}'(0))} \,{\cal B(\tau')}\,{\rm d}\tau',$$ C₁ = (A.1)

$$Y_1'(0)\int_{0}^{\Delta}\! \frac{({\vec Z}(\Delta)\sigma{\vec Y}(\tau'))}{({\vec Z}(\Delta)\sigma{\vec Y}'(0))} \,{\cal B(\tau')}\,{\rm d}\tau.$$ C₂ = (A.2)

Substituting them into (4) we have

$$J(\tau) -\frac{Y_1(\tau)Y_2'(0)}{({\vec Z}(\Delta)\sigma{\vec Y}'(0))} \i... ...Delta}\!({\vec Z}(\Delta)\sigma{\vec Y}(\tau')) \,{\cal B(\tau')}\,{\rm d}\tau'$$ =

$$\frac{Y_2(\tau)Y_1'(0)}{({\vec Z}(\Delta)\sigma{\vec Y}'(0))} \in... ...^{\tau}\!({\vec Y}(\tau)\sigma{\vec Y}(\tau')) \,{\cal B(\tau')}\,{\rm d}\tau',$$ +

or

 

$$J(\tau)=\frac{({\vec Y}'(0)\sigma{\vec Y}(\tau))}{({\vec Z}(\Delt... ...^{\tau}\!({\vec Y}(\tau)\sigma{\vec Y}(\tau')) \,{\cal B(\tau')}\,{\rm d}\tau'.$$ (A.3)

For arbitrary vectors ${\vec U}$, ${\vec V}$, ${\vec X}$, ${\vec Y}$ the following identity is valid

 

$$({\vec U}(u)\sigma{\vec V}(v)) ({\vec X}(x)\sigma{\vec Y}(y))=({\... ...igma{\vec V}(v))- ({\vec U}(u)\sigma{\vec X}(x))({\vec Y}(y)\sigma{\vec V}(v)).$$ (A.4)

Breaking the interval of integration in the first integral in (A.3) and using (A.4) we get (6)

Appendix B:

The Eqs. (4)-(5) can be written in a matrix form

 

$${\vec J}(\tau) {\vec R}(\tau){\vec C}+ {\vec R}(\tau)\sigma\int_0^\tau {\vec Y}(\tau')\,{\cal B}(\tau')\,{\rm d}\tau',$$ = (B.1)

with

$$\begin{eqnarray*}{\vec R}(\tau)=\left(\begin{array}{cc} Y_1(\tau) & Y_2(\tau) \\... ... {\vec C}=\left(\begin{array}{c} C_1 \\ C_2 \end{array}\right). \end{eqnarray*}$$

Taking into account the boundary conditions

  $$\begin{eqnarray*} {\vec J}(0)&=&\left(\begin{array}{c} J(0) \\ 0 \end{array}\... ...J(\Delta) \\ -J(\Delta) \end{array}\right)= J(\Delta){\vec g}, \end{eqnarray*}$$

we get the formal solution as

 

$${\vec J}(\tau) {\vec U}(\tau,0){\vec J}(0)+ {\vec R}(\tau)\sigma\int_0^\tau {\vec Y}(\tau'){\cal B}(\tau')\,{\rm d}\tau',$$ = (B.2)

where

  $$\begin{eqnarray*} {\vec U}(\tau_1,\tau_2)&=&{\vec R}(\tau_1) {\vec R}^{-1}(\tau_2),~~~~~~(\tau_1\geq\tau_2). \end{eqnarray*}$$

In the case of constant $\varepsilon (\tau )$ and $B(\tau)$ we obtain (16).

Let us now introduce new variables $\Phi$ and $\Psi$

 

$${\vec J}_j \Phi_j-\Psi_j,~~~~~~~(j=1,...,N)$$ = (B.3)

which obey the following recurrent equations

$$\begin{eqnarray*}\Phi_j&=&{\vec U}_j\Phi_{j-1},~~~~~~~~~~~~~~~\Phi_0={\vec J}_0,\\ \Psi_j&=&{\vec U}_j\Psi_{j-1}+{\vec K}_jB_{j-1}, ~~~\Psi_0=0. \end{eqnarray*}$$

Using notations

$$\begin{eqnarray*}{\vec U}(j,k)&=&{\vec U}(j,i+1){\vec U}(i,k),~~~~~ (k\leq i\leq j),\\ {\vec U}(j,j)&=&{\vec U}_j \end{eqnarray*}$$

we get

$$\begin{eqnarray*}\Phi_j&=&{\vec U}(j,1)\Phi_0,\\ \Psi_j&=&{\vec B}(j,1)={\vec U... ...B_1+...\\ &+&{\vec U}_j{\vec K}_{j-1}B_{j-2}+{\vec K}_jB_{j-1}. \end{eqnarray*}$$

The solution for jth interval thus becomes

 

$${\vec J}_j {\vec U}(j,1){\vec J}_0-{\vec B}(j,1).$$ = (B.4)

The value of ${\vec J}$ at the boundary reads

$$\begin{eqnarray*}{\vec J}_N&=&{J}_N{\vec g}={\vec U}(N,1) {\vec f}\,J_0-{\vec B}(N,1). \end{eqnarray*}$$

Multiplying both parts of the equation by the vector ${\vec e}$ we get

$$\begin{eqnarray*}J_0&=&{({\vec e}{\vec B}(N,1))\over ({\vec e}{\vec U}(N,1){\vec f})}\cdot \end{eqnarray*}$$

The substitution of J₀ into (B.4) gives

 

$$J_j\! \!=\! \!({\vec f}{\vec J}_j)\!=\! ({\vec f}{\vec U}(j,1){\vec f})\cdot\... ...\vec B}(N,1))\over({\vec e}{\vec U}(N,1){\vec f})} \!-\!({\vec f}{\vec B}(j,1))$$

$${A(N,j,1)\over({\vec e}{\vec U}(N,1){\vec f})}\cdot$$ = (B.5)

Taking into account the following identity

$$\begin{eqnarray*}{\vec B}(N,1)&=&{\vec U}(N,j+1){\vec B}(j,1)+{\vec B}(N,j+1) \end{eqnarray*}$$

with

$$\begin{eqnarray*}{\vec U}(N,N+1)=1,~~~~{\vec B}(N,N+1)=0 \end{eqnarray*}$$

and

$$\begin{eqnarray*}{\vec f}\cdot{\vec f}^\top+{\vec h}\cdot{\vec h}^\top=I \end{eqnarray*}$$

we get the expression of A(N,j,1)

$$\begin{eqnarray*}A(N,j,1)&=&({\vec f}{\vec U}(j,1){\vec f})({\vec e}{\vec B}(N,1... ...N,j+1))+({\vec e}{\vec U}(N,j+1){\vec h})({\vec f}{\vec W}(j,1)) \end{eqnarray*}$$

and thus Eq. (17).

$$\begin{eqnarray*}({\vec f}{\vec W}(j,1))&=&({\vec f}{\vec U}(j,1){\vec f})({\vec... ... B}(j,1))\\ &=&({\vec f}{\vec U}^\top(j,1)\sigma{\vec B}(j,1)), \end{eqnarray*}$$

where we used an identity

$$\begin{eqnarray*}{\vec f}\cdot{\vec h}^\top-{\vec h}\cdot{\vec f}^\top &=&\sigma. \end{eqnarray*}$$

Taking into account the following property

$$\begin{eqnarray*}{\vec U}^\top(j,i)\sigma{\vec U}(j,i)=\sigma, \end{eqnarray*}$$

we have

$$\begin{eqnarray*}({\vec f}{\vec W}(j,1))&=&({\vec f}{\vec U}^\top(j,1)\sigma[ {\... ..._2)\, B_1+...+({\vec f}{\vec U}^\top(j-1,1){\vec W}_j)\,B_{j-1}, \end{eqnarray*}$$

where ${\vec W}_j$ are defined by (19).

References

 

• Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (Dover Publication Inc, New York) In the text  
• Gierens, K. M., Traving, G., & Wehrse, R. 1986, JQRST, 37, 361 In the text  
• Gu, Y., Jefferies, J. T., Lindsey, C., & Avrett, E. H. 1997, ApJ, 484, 960 In the text NASA ADS  
• Kamke, E. 1965, Differentialgleichungen, I. Gewöhnliche Differentialgleichungen (Akadem. Verlagsges., Leipzig) In the text  
• Lindsey, C., & Jefferies, J. T. 1990, ApJ, 349, 286 In the text  
• Mihalas, D. 1978, Stellar Atmospheres, 2nd Ed. (Freeman, San Francisco) In the text  
• Nikoghossian, A. G., Pojoga, S., & Mouradian, Z. 1997, A&A, 325, 813 In the text  
• Peraiah, A. 1984, in Methods in Radiative Transfer, ed. W. Kalkofen (University Press, Cambridge), 281 In the text  
• Schmidt, M., & Wehrse, R. 1987, in Numerical Radiative Transfer, ed. W. Kalkofen (University Press, Cambridge), 341 In the text  
• Wehrse, R. 1981, MNRAS, 195, 553 In the text NASA ADS