text:
# An efficient algorithm for two–dimensional radiative transfer in axisymmetric circumstellar envelopes and disks

## 1 Introduction

Molecular line and dust continuum observations are an important tool for studying the envelopes and disks around young stellar objects (YSO), post-AGB stars and AGN. One of the main difficulties in interpreting such observations is that optical depths effects play an important role in the emission of this radiation. It is, for instance, well known that self-absorption and non-LTE effects are the main processes at work in shaping the characteristic asymmetric double-peaked emission line profiles from collapsing protostellar cores (Zhou, Zhou (1992)). Radiative transfer computations, at an appropriate level of complexity, are therefore needed in order to reconstruct the density, velocity and temperature structure of an observed cloud.

If densities drop below the critical density, the system may deviate from local thermodynamic equilibrium (LTE). Line trapping and photon escape from the line wings can in some cases be treated in the large-velocity-gradient (LVG) limit, but this approach is only valid if systematic velocity fields are much greater than the local turbulent line width. In all other cases one must perform a full non-LTE line transfer computation.

For problems that can be formulated in 1-D slab or spherical geometry there exist many radiative transfer programs, many of which use sophisticated techniques such as Accelerated Lambda Iteration (ALI; see a review by Hubeny Hubeny (1989)) and Complete Linearization (CL; Auer & Mihalas Auer and Mihalas (1969)). But the requirement of 1-D geometry is often too restrictive. Distinctly non-spherical features are often observed from YSO, such as bipolar reflection nebulae (e.g. Lenzen Lenzen (1987)), bipolar outflows (see Bachiller Bachiller (1996)) and disks (e.g. McCaughrean & O’Dell McCaughrean and O’Dell (1996)). Even the progenitors of these YSOs, starless dense cloud cores, seem to appear as elongated structures at millimeter wavelengths (e.g. Myers et al. Myers et al. (1991)), indicating that even in the early stages of star formation spherical symmetry does not apply. The case of post-AGB stars and Planetary Nebulae is just as compelling, with the majority of these nebulae being bipolar. The Cygnus Egg Nebula (CRL 2688) and te Red Rectangle (HD 44179) are perhaps the most spectacular examples of such bipolarity.

To model such objects, clearly one must resort to multi-dimensional transfer computations. There is a vast literature on this topic. Methods roughly fall in one of three catagories: Monte Carlo methods, Discrete Ordinate methods and Angular Moment methods. Monte Carlo codes are very flexible and can be used for a large variety of problems in multidimensional geometries, such as UV continuum transfer (e.g. Spaans Spaans (1996)), optical and infrared continuum transfer (Wolf et al. Wolf et al. (1999)), molecular line transfer (Hogerheijde Hogerheijde (1998)), and Compton scattering (e.g. Pozdnyakov, Sobol & Sunyaev 1977; Haardt & Maraschi 1991). Such methods perform well at low to medium optical depths, but it is well known that at high optical depths they converge very slowly.

Angular Moment methods, on the other hand, are very well suited to treat the high optical depths regime, since they are related to (or variants of) the diffusion equation (see e.g. Yorke et al. Yorke et al. (1993), Sonnhalter et al. Sonnhalter et al. (1995), Murray et al. Murray et al. (1994)). However, it is not surprising that they fail at low optical depth, since the diffusion approximation was never meant for this regime.

In the Discrete Ordinate approach, not only space is discretized into cells, but also the photon propagation direction. Most multi-dimensional implementations of the Discrete Ordinate methods are based on the “Lambda Iteration” scheme (e.g. Collison & Fix Collison and Fix (1991), Efstathiou & Rowan-Robinson Efstathiou and Rowan-Robinson (1991), Philips Philips (1999)). The advantage of these methods over the Monte Carlo approach is that they do not involve random noise, and therefore provide ‘cleaner’ answers. But they suffer from the same convergence problems as Monte Carlo methods. However, for Lambda Iteration there are various ways to cure this disease. The most well known of these methods is Accelerated Lambda Iteration (ALI, e.g. Scharmer Scharmer (1981), Rybicki & Hummer Rybicki and Hummer (1991)).

In this paper we will focus on the Discrete Ordinate approach to radiative transfer because of its versatility, accuracy and the wide range of convergence acceleration techniques available. However, despite the relative efficiency of these methods, multidimensional calculations remain costly. Feasibility constraints can pose severe limits on the spatial and angular resolution, which could easily result in unacceptable numerical diffusion. Also, this limits the number of models one can reasonably make to fit observations, which could lead to dangerous undersampling of the parameter space.

The bottleneck lies in the integration of the formal transfer equation. The most straightforward way of performing these integrals is by the method of “Long Characteristics”, which is accurate, but rather costly in CPU time. A more efficient algorithm for doing this in two dimensions is the method of “Short Characteritics” (SC; Mihalas et al. Mihalas et al. (1978), Kunasz & Auer Kunasz and Auer (1988), Auer & Paletou Auer and Paletou (1994), Stone et al. (Stone et al. (1992)). These algorithms are designed specifically with cartesian or cylindrical coordinates in mind, and are not straightforward to generalize to other coordinate systems. For circumstellar envelopes, however, there are several arguments favoring the use of spherical (polar) coordinates, as opposed to cylindrical coordinates. Most circumstellar nebulae have density and temperature profiles that are peaked towards the center. This means that the radiation field is dominated by photons emitted in the central regions, which are subsequently reprocessed in the outer parts of the nebula. The numerical scheme must therefore be able to resolve both the very concentrated central regions and the extended outer regions simultaneously. Also, it must guarantee that all radiation emitted at small radii will eventually emerge at large radii, which amounts to saying that flux must be conserved over a large range of radii. Using spherical coordinates and a logarithmic radial grid is the most natural way to cover such a large dynamic range and guarantee flux conservation.

The goal of this paper is to describe, test and demonstrate an algorithm that generalizes the Short Characteristics method to spherical coordinates111Just prior to submission we became aware of a paper by Busche & Hillier (Busche and Hillier (2000)), who describe a method of short characteristics in spherical coordinates that is quite similar to ours.. It is an algorithm specifically suited for axisymmetric circumstellar nebulae and disks. It has been implemented in a multi-purpose radiative transfer code named RADICAL, which is designed to perform 2-D computations in dust continuum emission/absorption, multiple scattering off dust grains, non-LTE line transfer, and (for application to X-ray binaries and AGN) Comptonization. In this paper we describe the method of short characteristics in spherical coordinates, and focus our attention to the cases of simple isotropic scattering off dust grains and non-LTE line transfer for multi-level molecules. For a more extensive discussion of the algorithm and its applications, see Dullemond (Dullemond (1999)).

The structure of this paper is as follows. In Section 2 we will present the equations of transfer we wish to solve. In Section 3 we shall review the method of short characteristics as it is often presented in the literature. In Section 4 we will show how this method can be generalized to spherical coordinates. Then we will put the algorithm to the test in Section 5. Finally we will present two example applications in Sections 6 and 7.

## 2 Equations of radiative transfer

The method that will be described in this paper is an Accelerated Lambda Iteration method. In such an algorithm the integration of the formal transfer equation is performed using a “Lambda Operator”. In this section we will present the equations that are to be solved, and we define the Lambda Operator. The numerical details of the Lambda Operator will be given in the next sections.

The formal transfer equation is

$$\frac{dI_\nu }{ds}=\alpha _\nu (S_\nu I_\nu ),$$
(1)
with $`I_\nu `$ is the intensity, $`S_\nu `$ the source function, $`\alpha _\nu `$ the opacity, and $`s`$ the path length. This equation must hold along every straight line through the medium. Its integral form along a ray through a point $`P`$ reads:

$$I_\nu (P)=e^{\tau _\nu }I_\nu (0)+_0^{\tau _\nu }e^{\tau _\nu ^{}}S(\tau _\nu ^{})𝑑\tau _\nu ^{},$$
(2)
where $`\tau _\nu `$ is the optical depth along the ray, between point $`P`$ and the edge of the medium. After evaluating this integral for all angles $`𝝎`$, one can compute the mean intensity $`J_\nu `$

$$J_\nu =\frac{1}{4\pi }I_\nu (𝝎)𝑑𝝎.$$
(3)
The entire operation of computing $`J_\nu (P)`$ at every point $`P`$, for a given source function $`S_\nu `$, can be written as the action of a linear Lambda Operator $`\mathrm{\Lambda }`$:

$$J_\nu =\mathrm{\Lambda }[S_\nu ],$$
(4)
Using this Lambda Operator we can write down the complete transfer equation for a simple problem of thermal emission and isotropic (dust) scattering

$$S_\nu =ϵB_\nu (T)+(1ϵ)\mathrm{\Lambda }[S_\nu ],$$
(5)
where $`ϵ\alpha _\nu ^{\text{abs}}/\alpha _\nu `$ is the thermalization coefficient (with $`\alpha _\nu ^{\text{abs}}`$ the thermal absorption opacity), and $`B_\nu (T)`$ is the Planck function. Solving the transfer problem for isotropic scattering and thermal emission amounts to solving Eq.(5) for $`S_\nu `$. The Lambda Iteration procedure amounts to iteratively applying the Lambda Operator and computing the new $`S_\nu `$ until convergence is reached. The Accelerated Lambda Iteration procedure, which converges much faster, is a variant of this procedure, involving an approximate operator $`\mathrm{\Lambda }^{}`$. For details we refer to Hubeny (Hubeny (1989)) and Rutten (Rutten (1999)).

For multi-level line transfer, we follow the treatment of Rybicki & Hummer (Rybicki and Hummer (1991)). Consider an atom or molecule having $`N`$ levels, with spontaneous radiative downward transition rates $`A_{ij}`$, Einstein coefficients $`B_{ij}`$ and collision rates $`C_{ij}`$ between levels $`i`$ and $`j`$. The line profile function $`\stackrel{~}{\phi }_{ij}(\nu )`$ determines at which frequencies the line emits and absorbs. When no systematic fluid velocities are present, the line profile function is isotropic, and is normalized to unity. For the application to circumstellar envelopes, the dominant broadening mechanisms are turbulent and thermal broadening. These two mechanisms produce a Gaussian profile:

$$\stackrel{~}{\phi }_{ij}(\nu )=\frac{c}{a_{\text{tot}}\nu _{ij}\sqrt{\pi }}\mathrm{exp}\left(\frac{c^2(\nu \nu _{ij})^2}{a_{\text{tot}}^2\nu _{ij}^2}\right).$$
(6)
Here $`c`$ is the speed of light, $`\nu _{ij}`$ the line-center frequency of the transition between levels $`i`$ and $`j`$, and $`a_{\text{tot}}`$ is the line width,

$$a_{\text{tot}}=a_{\text{turb}}+\sqrt{\frac{2kT_{\text{kin}}}{m_{\text{mol}}}},$$
(7)
where $`T_{\text{kin}}`$ is the (kinetic) temperature of the gas, $`m_{\text{mol}}`$ the mass of the molecule, and $`a_{\text{turb}}`$ is the turbulent line width. A systematic fluid velocity can cause the line profile function to be angle-dependent in the lab frame as a result of Doppler shift,

$$\phi _{ij}(𝝎,\nu )=\stackrel{~}{\phi }_{ij}\left(\nu (1𝝎𝐯/c)\nu _{ij}\right),$$
(8)
The opacity in the line associated with this line profile is:

$$\alpha _{ij}(𝝎,\nu )=\frac{h\nu }{4\pi }N(n_jB_{ji}n_iB_{ij})\phi _{ij}(𝝎,\nu ),$$
(9)
where $`n_i`$ are the fractional level populations, and $`N`$ the number density of molecules. We assume complete redistribution for the lines. The source function is then independent of frequency and angle:

$$S_{ij}=\frac{n_iA_{ij}}{n_jB_{ji}n_iB_{ij}}.$$
(10)
The transfer equation for this source function is then

$$\frac{dI_{ij}(𝝎,\nu )}{ds}=\alpha _{ij}(𝝎,\nu )\left[S_{ij}I_{ij}(𝝎,\nu )\right],$$
(11)
where we assume non-overlapping lines.

The source term $`S_{ij}`$ is known once the fractional level populations $`n_i`$ are known. They are a solution of the statistical equilibrium equation. Using the definition of the line-integrated Lambda Operator $`\overline{\mathrm{\Lambda }}_{ij}[S_{ij}]\overline{J}_{ij}`$, with

$$\overline{J}_{ij}=\frac{1}{4\pi }I(𝝎,\nu )\phi _{ij}(𝝎,\nu )𝑑𝝎𝑑\nu ,$$
(12)
the statistical equilibrium equations become:

$$\begin{array}{cc}\hfill \underset{j>l}{}& \left[n_jA_{jl}+(n_jB_{jl}n_lB_{lj})\overline{\mathrm{\Lambda }}_{jl}[S_{jl}]\right]\hfill \\ \hfill & \underset{j<l}{}\left[n_lA_{lj}+(n_lB_{lj}n_jB_{jl})\overline{\mathrm{\Lambda }}_{lj}[S_{lj}]\right]\hfill \\ \hfill +& \underset{j}{}\left[n_jC_{jl}n_lC_{lj}\right]=0.\hfill \end{array}$$
(13)
The non-locality of radiative transfer is now hidden in the $`\overline{\mathrm{\Lambda }}_{jl}`$ operator, so that Eq.(13) now represents the complete (non-linear) set of equations for line transfer. Lambda Iteration now proceeds by iteratively applying the $`\overline{\mathrm{\Lambda }}_{jl}`$ operator and solving the matrix equation represented by Eq.(13). Accelerated Lambda Iteration proceeds according to the MALI scheme of Rybicki & Hummer (Rybicki and Hummer (1991)).

## 3 Short characteristics in cartesian coordinates

To carry out the Lambda Iteration or Accelerated Lambda Iteration procedure, we need a numerical implementation of the Lambda Operator $`\mathrm{\Lambda }`$. In Cartesian coordinates, the formal transfer equation Eq.(1) becomes

$$\frac{dI_\nu }{ds}\omega _x\frac{I_\nu }{x}+\omega _y\frac{I_\nu }{y}=\alpha _\nu (S_\nu I_\nu )$$
(14)
where translational symmetry in the $`z`$–direction was assumed.

The numerical implementation of the Lambda Operator amounts to integrating Eq.(14) for given $`S_\nu `$ and $`\alpha _\nu `$. This must be done on a 2–dimensional spatial grid $`𝐱=(x_i,y_j)`$, for a discrete set of directions $`𝝎=\{𝝎_i\}`$ and frequencies $`\nu =\{\nu _i\}`$. This will provide the specific intensity $`I(x_i,y_j;𝝎{}_{k}{}^{},\nu _l)`$ for all $`i,j,k,l`$. Let us focus on a given point $`P=(x_i,y_j)`$ and on a single direction and frequency, $`𝝎=𝝎{}_{k}{}^{},\nu =\nu _l`$. The integral of Eq.(14) can be performed numerically along the entire characteristic starting at the upstream boundary, heading in the downstream direction (i.e. the direction where the radiation comes from) and ending at point $`P`$ (see Fig. 1).

This direct approach is called the method of Long Characteristics (LC). Provided the discretization in angle $`𝝎`$ is appropriate, this method is quite accurate and reliable. But it has a computational redundancy, and hence it is overly time–consuming. Consider, for instance, a spatial grid $`N\times N`$, a set of $`N_𝝎`$ directions and of $`N_\nu `$ frequencies. The long characteristics integral of Eq.(14) typically requires in the order of $`N`$ integration steps. This means that while the dimension of the grid is $`N^2\times N_𝝎\times N_\nu `$, the total computational time scales as

$$t_{CPU}N^3\times N_𝝎\times N_\nu .$$
(15)
The Short Characteristics method of integration (SC; see Fig. 2) does not have this redundancy.

Instead of performing the integral along the entire ray (the long characteristic), we perform the integral only along that portion of the ray (the short characteristic) which connects a point $`U`$ on the grid upstream of $`P`$ to the closest intersection downstream of $`P`$ itself. The intensity at $`P`$ is given by

$`I_\nu (P;𝝎)`$ $`=`$ $`e^{\tau _\nu }I_\nu (U;𝝎)+`$ (16)
$`{\displaystyle _0^{\tau _\nu }}e^{\tau _\nu ^{}}S_\nu (𝐱(\tau _\nu ^{});𝝎)𝑑\tau _\nu ^{}`$
where $`\tau _\nu `$ is the optical depth between points $`U`$ and $`P`$. The upstream intensity $`I_\nu (U;𝝎)`$ can be found from the intensities at $`A`$, $`B`$ and $`C`$ by 3–point quadratic interpolation

$$I_\nu (U;𝝎)=aI_\nu (A;𝝎)+bI_\nu (B;𝝎)+cI_\nu (C;𝝎)$$
(17)
where $`a`$, $`b`$ and $`c`$ are the usual Lagrange coefficients for polynomial interpolation. Quadratic or higher order interpolation is necessary in order to reproduce the diffusion limit for high optical depth, which is governed by a second order partial differential equation.

The integral from $`U`$ to $`P`$ can be computed with second order accuracy by interpolating the source function $`S_\nu (𝐱(\tau _\nu ^{});𝝎)`$ between the points $`D`$, $`P`$ and $`U`$. Following Olson & Kunasz (Olson and Kunasz (1987)), one finds

$`I_\nu (P;𝝎)`$ $`=`$ $`e^{\tau _\nu }I_\nu (U;𝝎)+u_\nu S_\nu (U;𝝎)`$ (18)
$`+p_\nu S_\nu (P;𝝎)+d_\nu S_\nu (D;𝝎),`$
with

$`u_\nu `$ $`=`$ $`e_0+[e_2(2\tau _\nu +\overline{\tau }_\nu )e_1]/[\tau _\nu (\tau _\nu +\overline{\tau }_\nu )]`$ (19)
$`p_\nu `$ $`=`$ $`[(\tau _\nu +\overline{\tau }_\nu )e_1e_2]/[\tau _\nu \overline{\tau }_\nu ]`$ (20)
$`d_\nu `$ $`=`$ $`[e_2\tau _\nu e_1]/[\overline{\tau }_\nu (\tau _\nu +\overline{\tau }_\nu )]`$ (21)
$`e_0`$ $`=`$ $`1e^{\tau _\nu }`$ (22)
$`e_1`$ $`=`$ $`\tau _\nu e_0`$ (23)
$`e_2`$ $`=`$ $`\tau _\nu ^22e_1`$ (24)
where $`\tau _\nu `$ and $`\overline{\tau }_\nu `$ are the depths at $`U`$ and $`D`$, respectively. It should be noted that this quadrature formula may have pathological behaviour if the source function and/or the opacity varies strongly between the points $`U`$, $`P`$ and $`D`$. This problem can be solved by limiting the resulting integrals between zero and $`\mathrm{max}(j_\nu (P),j_\nu (U))\mathrm{\Delta }s`$, where $`\mathrm{\Delta }s`$ is the path length along the short characteristic, and $`j_\nu =\alpha _\nu S_\nu `$.

By systematically performing the integrals over all the short characteristics, one can find an approximate formal solution of the transfer equation (Kunasz & Auer Kunasz and Auer (1988), Auer & Paletou Auer and Paletou (1994), Auer et al. Auer et al. (1994)). A key ingredient for the SC method to work is that the integrals should be performed in the right order, so that the upstream intensities $`I_\nu (A;𝝎)`$, $`I_\nu (B;𝝎)`$ and $`I_\nu (C;𝝎)`$ are known before the integral is performed and Eq.(18) evaluated. In order to do so, the grid must be swept from the two upstream boundaries towards the two downstream boundaries.

The method of Short Characteristics is computationally less time consuming than the method of Long Characteristics, because now the transfer integral is performed over a much shorter path. For the same discretization introduced earlier in this section, the computational time scales as

$$t_{CPU}N^2\times N_𝝎\times N_\nu $$
(25)
which is typically an factor $`N`$ shorter than in the case of Long Characteristics.

## 4 Short characteristics in spherical coordinates

We now wish to generalize the Short Characteristics to spherical coordinates. In the following we refer to a standard spherical coordinate system $`(R,\mathrm{\Theta },\mathrm{\Phi })`$ where $`\mathrm{\Theta }`$ is the latitude and $`\mathrm{\Phi }`$ the azimuth. By assuming axial symmetry, any dependence on $`\mathrm{\Phi }`$ is suppressed, although radiation is still allowed to travel along $`/\mathrm{\Phi }`$, as well as in the radial and meridional directions.

In order to describe the radiation field at each spatial point $`P=(R,\mathrm{\Theta })`$ we need to set up a local coordinate system to characterize the photon direction at $`P`$. We introduce two independent angles on the sky of the local observer: $`\theta `$ and $`\varphi `$. The north pole of this local sky-map is chosen to coincide with the outward–pointing radial direction. The $`\varphi `$ angle is gauged in such a way that $`\varphi =0`$ points parallel to the equator of the global coordinate system (see Fig. 3).

As is customary in transfer theory, we use $`\mu \mathrm{cos}\theta `$ instead of $`\theta `$ itself, so the specific intensity depends upon the two spatial variables $`R,\mathrm{\Theta }`$, the ray direction $`\mu ,\varphi `$, and the frequency $`\nu `$, $`I=I_\nu (R,\mathrm{\Theta };\mu ,\varphi )`$. The transfer equation, Eq.(1), in spherical coordinates reads

$$\begin{array}{cc}\hfill \frac{dI_\nu }{ds}\mu \frac{I_\nu }{R}& \frac{\sqrt{1\mu ^2}}{R}\mathrm{sin}\varphi \frac{I_\nu }{\mathrm{\Theta }}+\frac{1\mu ^2}{R}\frac{I_\nu }{\mu }\hfill \\ & \frac{\mathrm{cos}\varphi }{\mathrm{tan}\mathrm{\Theta }}\frac{\sqrt{1\mu ^2}}{R}\frac{I_\nu }{\varphi }=\alpha _\nu (S_\nu I_\nu ).\hfill \end{array}$$
(26)
An important consequence of the use of spherical coordinates is that, contrary to what happens for cartesian coordinates, the photon angles $`\mu `$ and $`\varphi `$ are no longer constant along the rays. The variation of $`R`$, $`\mathrm{\Theta }`$, $`\mu `$ and $`\varphi `$ along the path are

$`{\displaystyle \frac{dR}{ds}}`$ $`=\mu ,`$ $`{\displaystyle \frac{d\mathrm{\Theta }}{ds}}`$ $`={\displaystyle \frac{\sqrt{1\mu ^2}}{R}}\mathrm{sin}\varphi ,`$ (27)
$`{\displaystyle \frac{d\mu }{ds}}`$ $`={\displaystyle \frac{1\mu ^2}{R}},`$ $`{\displaystyle \frac{d\varphi }{ds}}`$ $`={\displaystyle \frac{\sqrt{1\mu ^2}}{R}}{\displaystyle \frac{\mathrm{cos}\varphi }{\mathrm{tan}\mathrm{\Theta }}},`$ (28)
where $`s`$ is the path length. Solving these equations yields

$`R^2`$ $`=b^2+s^2,`$ (29)
$`\mathrm{cos}\mathrm{\Theta }`$ $`={\displaystyle \frac{z_0+s\mathrm{cos}\mathrm{\Theta }_{\mathrm{}}}{\sqrt{b^2+s^2}}},`$ (30)
$`\mu `$ $`={\displaystyle \frac{s}{\sqrt{b^2+s^2}}},`$ (31)
$`\mathrm{sin}\varphi `$ $`={\displaystyle \frac{b^2\mathrm{cos}\mathrm{\Theta }_{\mathrm{}}z_0s}{b\sqrt{b^2+s^2(z_0+s\mathrm{cos}\mathrm{\Theta }_{\mathrm{}})^2}}},`$ (32)
where $`b`$ is the impact parameter of the ray with respect to the origin, $`z_0`$ is the height above the midplane of closest approach to the symmetry-axis, and $`\mathrm{\Theta }_0`$ is the inclination at infinity. When projected into the subspace spanned by $`R,\mathrm{\Theta }`$ the trajectory becomes a hyperbola, as is shown in Fig. 4. We stress that this shape is caused by eliminating the dependence on the $`\mathrm{\Phi }`$–angle, and is purely a projection effect.

For the numerical implementation of the short characteristics scheme, we are interested in those characteristics that pass through a grid point $`P=(R_k,\mathrm{\Theta }_l)`$ and are tangent to one of the local discrete ordinates $`(\mu _i,\varphi _j)`$. Clearly, once $`R_k,\mathrm{\Theta }_l,\mu _i,\varphi _j`$ are fixed, such a characteristic is unique and the values of its parameters are

$`b^2`$ $`=R_k^2(1\mu _i^2)`$ (33)
$`\mathrm{cos}\mathrm{\Theta }_{\mathrm{}}`$ $`=\mu _i\mathrm{cos}\mathrm{\Theta }_l+\sqrt{1\mu _i^2}\mathrm{sin}\mathrm{\Theta }_l\mathrm{sin}\varphi _j`$ (34)
$`z_0`$ $`=R_k[(1\mu _i^2)\mathrm{cos}\mathrm{\Theta }_l`$
$`\mu _i\sqrt{1\mu _i^2}\mathrm{sin}\mathrm{\Theta }_l\mathrm{sin}\varphi _j].`$ (35)
The short characteristic passing through $`(R_k,\mathrm{\Theta }_l;\mu _i,\varphi _j)`$ is defined as the section of this curve that starts at the closest intersection with the grid lines upstream of $`P`$ (point $`U`$), passes through $`P`$ and ends at the closest intersection with the grid lines downstream of $`P`$ (point $`D`$). The location of the points $`U`$ and $`D`$ is specified by the corresponding values of parameter $`s`$ along the ray, $`s_U`$ and $`s_D`$, which are found solving equations (29)–(30) with $`R=R_K`$ and $`\mathrm{\Theta }=\mathrm{\Theta }_L`$, where $`K=k1,k,k+1`$ and $`L=l1,l,l+1`$. Both $`R=R_k`$ and $`\mathrm{\Theta }=\mathrm{\Theta }_l`$ need to be included because the characteristic may intersect the same $`\mathrm{\Theta }`$ or $`R`$ grid line twice. In principle, each equation has two solutions for a given value of $`K`$ and $`L`$, giving 12 possible roots

$`s_{1\mathrm{}6}`$ $`=`$ $`\pm \sqrt{R_K^2b^2}`$ (36)
$`s_{7\mathrm{}12}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{cos}^2\mathrm{\Theta }_{\mathrm{}}\mathrm{cos}^2\mathrm{\Theta }_L}}\{z_0\mathrm{cos}\mathrm{\Theta }_{\mathrm{}}\pm `$ (37)
$`\mathrm{cos}\mathrm{\Theta }_L\sqrt{b^2(\mathrm{cos}^2\mathrm{\Theta }_{\mathrm{}}\mathrm{cos}^2\mathrm{\Theta }_L)+z_0^2}\}.`$
However, two of these solutions always give $`s=s_P`$, i.e. $`P=(R_k,\mathrm{\Theta }_l)`$ itself, and are of no interest. Of the remaining 10 roots, some are complex and must be rejected. Between the real solutions, the one representing point $`D`$ ($`U`$) is selected asking that $`s>s_P`$ ($`s<s_P`$) and that $`|ss_P|`$ is minimum. For convenience, in the following we will denote with $`\stackrel{~}{R}_1`$, $`\stackrel{~}{\mathrm{\Theta }}_1`$, $`\stackrel{~}{\mu }_1`$ and $`\stackrel{~}{\varphi }_1`$ the values of the independent variables along the ray at point $`U`$.

Although, as we have just shown, short characteristics can be easily defined in spherical coordinates, two major problems have to be solved before they can be of any use in building a transfer algorithm. The first point concerns the fact that, as it was mentioned earlier, $`\mu `$ and $`\varphi `$ change along the ray. This means that in addition to spatial interpolation (see Eq. 17), we are forced to interpolate in $`\mu `$ and $`\varphi `$ as well in order to evaluate $`I_\nu (U,𝝎)`$. This is because the intensities at points $`A`$, $`B`$ and $`C`$ are known only for a discrete set of directions which are different, in general, from ($`\stackrel{~}{\mu }_1,\stackrel{~}{\varphi }_1)`$.

The second, more fundamental difficulty arises because in spherical coordinates the concept of upstream and downstream boundaries is different from the Cartesian case. Radial infinity is both the upstream and the downstream boundary, while in $`\mathrm{\Theta }`$ there is no obvious upstream or downstream boundary. If the grid is swept from $`\mathrm{\Theta }=0`$ to $`\mathrm{\Theta }=\pi `$ or vice versa, one will encounter situations in which the intensity at one of the points $`A`$, $`B`$, $`C`$ is not known before the evaluation the transfer integral along the short characteristic (Eq. 16) is performed. An example of such a situation is shown in Fig. 5. Interpolation makes use of the points $`A`$, $`B`$ and $`C`$, but since $`A`$ coincides with $`P`$, the intensity at the point $`AP`$, has not been computed yet.

### 4.1 Extended Short Characteristics

The problem of unknown intensities can be solved by modifying the definition of short characteristics to be the part of the ray that connects $`P`$, not with just the nearest gridline intersection, but with the nearest $`R=R_k`$ gridline intersection, i.e. the nearest radial shell. Such an “extended short characteristic” (ESC) is illustrated in Fig. 6.

The starting point $`U`$ of such an ESC will be located either at $`R=R_{k1}`$, $`R=R_{k+1}`$ or back at $`R=R_k`$. This means that in between $`U`$ and $`P`$ the ESC may intersect one or more $`\mathrm{\Theta }`$ grid lines. The point $`D`$ on the downstream side remains the same as for standard SCs.

By using ESC instead of SC the “problem of unknown upstream intensities” can be eliminated. In fact, if a proper sweeping scheme is chosen (see Subsection 4.2), the problem of unknown intensities only occurs in those situations when a short characteristic curves back onto the same $`\mathrm{\Theta }`$-gridline from which it originates, as is illustrated in Fig. 5. By extending only those short characteristics, and leaving the rest truly short, one can also avoid unknown intensities in the sweeping scheme. Just for notation we call this scheme the Minimally Extended Short Characteristics scheme (MESC). MESC is almost as accurate as ESC, but ESC is more closely similar to its one-dimensional spherical analogues, and is slightly less numerically diffusive.

In the following we denote with $`D`$ (or $`U_1`$) the single downstream intersection with a grid line, and with $`U_i(i=1,\mathrm{},m)`$ the multiple intersections upstream of $`P`$. The point $`U_m`$ is therefore the true upstream starting point of the ESC, where the intensity must be found by interpolation. In Fig. 6 this is the point $`U_2`$ and the ESC consists of two segments in this case.

### 4.2 The sweeping scheme

Using the (minimally) extended short characteristics defined above, we can systematically sweep the grid without encountering unknown intensities. We start at the outer boundary $`R_{\mathrm{}}`$ and integrate inwards only those ESCs for which $`\mu _i0`$. The sweeping order in $`\mathrm{\Theta }`$ is from $`\mathrm{\Theta }=0`$ to $`\mathrm{\Theta }=\pi `$ and then back.

The intensity at each $`P=(R_k,\mathrm{\Theta }_l)`$ is found for all $`\mu _i0`$ and $`\varphi _j`$ by tracing the ESCs back to their upstream starting point $`U_m`$. At point $`U_m`$ the values of $`R`$, $`\mathrm{\Theta }`$, $`\mu `$ and $`\varphi `$ are different from those at $`P`$, and will be denoted with $`\stackrel{~}{R}_m`$, $`\stackrel{~}{\mathrm{\Theta }}_m`$, $`\stackrel{~}{\mu }_m`$ and $`\stackrel{~}{\varphi }_m`$. The find the intensity at $`(\stackrel{~}{R}_m,\stackrel{~}{\mathrm{\Theta }}_m,\stackrel{~}{\mu }_m,\stackrel{~}{\varphi }_m))`$ by interpolation. We then integrate the formal transfer equation along each segment of the ESC connecting $`U_m`$ to $`P`$, according to Eq.(18). This gives

$$\begin{array}{cc}\hfill I_\nu (P)=& \mathrm{exp}(\tau _m)I_\nu (U_m)+\underset{i=1,m}{}\mathrm{exp}(\tau _{i1})\times \hfill \\ & \left(u_{\nu ,i}S_{\nu ,i}+p_{\nu ,i1}S_{\nu ,i1}+d_{\nu ,i2}S_{\nu ,i2}\right)\hfill \end{array}$$
(38)
where $`\tau _i`$ is the depth from $`P`$ to $`U_i`$ and the index $`i`$ denotes the quantity evaluated at $`U_i`$ (e.g. $`S_iS(U_i)`$; $`i=0`$ refers to $`P`$ and $`i=1`$ to $`D`$). The $`u`$’s, $`p`$’s and $`d`$’s are defined as in equations (19), (20) and (21), but with $`\tau `$ replaced by $`(\tau _i\tau _{i1})`$ and $`\overline{\tau }`$ by $`(\tau _{i1}\tau _{i2})`$.

The integration is then repeated moving towards smaller radii, until the inner boundary at $`R=R_1`$ is reached. Here we can include the contribution of a central source or any other boundary condition.

Then we start integrating back towards larger radii, until we reach the outer edge. By now the radiation field on the grid $`I_\nu (R_k,\mathrm{\Theta }_l;\mu _i,\varphi _j)`$ is known.

### 4.3 Tangent-ray discretization of photon direction

Equations (31) and (32) show that $`\mu `$ and $`\varphi `$ change along the ESC. If we follow the ESC upstream towards a point $`U_m`$ (see Fig. 6), then the values of these two angles at $`U_m`$ are generally not exactly at the discrete values $`\{\mu _i\}`$ and $`\{\varphi _j\}`$ of the sample of directions. This means that we must interpolate not only in space (between $`A`$, $`B`$ and $`C`$), but also in direction $`\mu `$, $`\varphi `$. Although, this is not, in principle, a fundamental problem, the use of interpolations should be reduced to a minimum to avoid unnecessary numerical diffusion.

Fortunately one can eliminate the interpolation in $`\mu `$ by means of a suitable choice of the $`\mu `$–grid so that all ESCs always start and end at one of the $`\{\mu _i\}`$ points. We let the $`\mu `$–discretization depend on $`R_k`$,

$$\{\mu _{k,j};j=m_k,\mathrm{},m_k\},$$
(39)
and choose the $`\mu _{k,j}`$ in such a way that for each $`j=m_k,\mathrm{},m_k`$ there is a $`j^{}=m_{k1},\mathrm{},m_{k1}`$ such that

$$(1\mu _{k,j}^2)R_k^2=(1\mu _{k1,j^{}}^2)R_{k1}^2.$$
(40)
This choice is based on the fact that the values of $`\mu `$ and $`R`$ along an ESC depend only on each other and on $`b`$, as can be seen by combining equations (29) and (31):

$$(1\mu ^2)R^2=b{}_{}{}^{2}.$$
(41)
By choosing the $`\mu _{k,j}`$ according to Eq.(40), a ray which originates at $`R_k`$ with, say, $`\mu =\mu _{k,j}`$ and arbitrary $`\varphi `$, reaches any other radius along the path with a value of $`\mu `$ which coincides with one of the points of the local $`\mu `$ grid there, thus eliminating the need for interpolation.

A $`\mu `$-grid that is consistent with Eq.(40) is:

$$\mu _{k,\pm i}=\pm \sqrt{1\frac{R_{ki}^2}{R_k^2}},$$
(42)
in agreement with the angular spacing induced in spherical symmetry by the “tangent ray method” (see e.g. Mihalas Mihalas et al. (1975); Zane et al. Zane et al. (1996)). Actually, it can be easily shown that in 1-D spherical symmetry the ESC method, with $`\mu _{k,j}`$ given by Eq.(42), is fully equivalent to the tangent ray method. This is an important feature of the algorithm since it then exactly recognizes spherical symmetry. And, even in the absence of spherical symmetry, it transports radiation outward without any numerical diffusion in $`R`$ or $`\mu `$.

However, the $`\mu `$ spacing implied by Eq.(42) has the tendency to give a poor sampling around $`\mu =0`$. This problem can be easily solved by introducing some (typically one or two) extra points around $`\mu 0`$ to enhance the angular resolution there. Obviously this violates the original prescription and therefore requires the use of interpolation for these extra $`\mu `$–points, producing a small amount of angular diffusion for $`\mu 0`$. Generally this diffusion is small.

Unfortunately the interpolation in $`\varphi `$ can never be avoided. The $`\varphi `$ angle depends in a complicated way on $`s`$ (see Eq. 32) and it can change rapidly even within one element of a ESC. Both first and higher order interpolation in $`\varphi `$ have been tested in our numerical code. We have found that in most cases the $`\varphi `$ diffusion is not very large and, in general, influences the solution less than the spatial ($`\mathrm{\Theta }`$) diffusion.

### 4.4 Special treatment of radiation near $`\mu 1`$

The tangent-ray discretization of $`\mu `$ allows the algorithm to accurately conserve radial flux. However, such a choice of $`\mu `$–angles requires a large number of $`\mu `$ points at larger radii, typically $`m_kk`$. One cannot make do with a smaller number of $`\mu `$ points without facing the risk of loosing flux. This is illustrated in the following argument. If radiation is emitted at a radius $`R`$, an observer at $`R_kR`$ can see the radiation from the emitting region even if its eyes cannot resolve the source. This is because the observer’s eyes measure the flux and not the intensity. The ESC algorithm, on the other hand, deals with intensity, and intensity is converted into flux by performing an integral over $`d\mu d\varphi `$. For this integral to be reasonably accurate, the emitting region must be resolved in $`\varphi `$ and $`\mu `$, leading to the requirement $`m_kk`$.

Unfortunately this means that the computational cost scales as $`N_R^2`$ if one wishes to extend the span of the radial domain. Since the ability to deal with many orders of magnitude in radius is crucial to solving transfer problems in circumstellar envelopes, this scaling is undesirable.

An easy way to solve this scaling problem is related to the simple observation that all photons with $`\mu 1`$ follow roughly a radially outgoing trajectory and they tend to travel more radially (i.e. with $`\mu `$ closer to unity), the further they propagate outwards. In the “radial streaming” limit ($`1\mu 1`$), Eq.(40) becomes approximately

$$\frac{1\mu _{k,j}}{1\mu _{kl,jl}}=\frac{\mathrm{\Omega }_{k,j}}{\mathrm{\Omega }_{kl,jl}}\frac{R_{kl}^2}{R_k^2}$$
(43)
where $`\mathrm{\Omega }_{k,j}`$ is the solid angle (bounded by $`\mu _{k,j}`$) at $`R_k`$. Eq.(43) is just a restatement of the $`1/R^2`$ law which is exactly obeyed by a point source and by any radiation in the radial streaming limit.

This property of radially outward radiation makes it possible to bundle all $`\mu `$–points with sufficiently large $`\mu `$ into a single collective flux-like bin. The intensity of that bin will be treated as the average intensity within that collective bin. The idea is to divide the $`\mu `$–range $`[1,1]`$ into three parts

$`[1,\mu _{rs}]`$ inward intensity bin: $`\mu `$ $`1`$
$`[\mu _{rs},\mu _{rs}]`$ intensity samples: $`\mu `$ $`=\mu _{k,j}`$
$`[\mu _{rs},1]`$ outward averaged bin: $`\mu `$ $`1.`$
This way the number of $`\mu `$–gridpoints $`m_k`$ at each radius can be limited, depending on how close $`\mu _{rs}`$ is to unity. We choose a global value for $`\mu _{rs}`$, and do not allow this to differ from one radius to another.

Because the radial outward bin represents an integrated intensity, i.e. an average of the true intensity over a solid angle $`\mathrm{\Omega }_{rs}=2\pi (1\mu _{rs})`$, it requires a special treatment. Let us denote the average intensity in this bin as $`I_\nu ^+(R_k,\mathrm{\Theta }_l)`$. The integration formula, Eq.(38), for $`I_\nu ^+(R_k,\mathrm{\Theta }_l)`$ becomes

$$\begin{array}{cc}\hfill I_\nu ^+(& R_k,\mathrm{\Theta }_l)=\mathrm{exp}(\tau )[\frac{R_{k1}^2}{R_k^2}I_\nu ^+(R_{k1},\mathrm{\Theta }_l)+\hfill \\ & (1\frac{R_{k1}^2}{R_k^2})\frac{1}{2\pi }_0^{2\pi }I_\nu (R_{k1},\mathrm{\Theta };\mu _{m_k1},\varphi )d\varphi ]\hfill \\ & +u_\nu S_\nu (R_{k1},\mathrm{\Theta }_l)+p_\nu S_\nu (R_k,\mathrm{\Theta }_l)\hfill \\ & +d_\nu S_\nu (R_{k+1},\mathrm{\Theta }_l).\hfill \end{array}$$
(44)
This formula simulates the correct $`1/R^2`$ behavior of the flux for optically thin media provided $`\mu _{rs}`$ is sufficiently close to unity. It reduces to the standard expression when the medium is optically thick. An estimate of the error introduced by the assumption of radial streaming can be made by comparing equations (40) and (43), and is of order $`(1\mu _{rs})/2`$. The radiation contained in the $`\mu =1`$ bin will not accumulate any numerical interpolation errors because it moves strictly along the radial grid lines.

The inward collective bin will always behave as a real intensity, so that the $`1/R^2`$ behavior does not need to be simulated.

### 4.5 Spectra and images

Once the iterative part of the transfer has been completed and the source function is known, the next step is to produce images and spectra. An image is produced by formal integration of the source function along long characteristics through the medium (ray tracing). Each ray represents one pixel of the image. One can produce spectra by making images at a range of frequencies, and integrating these images over the “detector” aperture.

Here, as in the Lambda iteration, we face resolution problems if the source under consideration spans a large range in $`\mathrm{log}(R)`$. The central parts of the image are often much brighter than the rest, but cover a much smaller fraction of the image. The spectrum may therefore contain significant contributions of flux from both the central parts and the outer regions of the image. Unless the image resolve all spatial scales of the object, the spectra produced in such a way are unreliable.

If a rectangular arrangement of pixels is used, one must make sure to use a variable spacing in both $`x`$ and $`y`$, in such a way that the small scales around the star are sufficiently resolved. If one is mainly interested in the images themselves, this seems the most reliable and straightforward way to go.

For the production of spectra we propose a different approach. Rather than arranging the pixels over a rectangle, we arrange them in concentric rings. The impact parameters of the circles are related to the radial grid points of the transfer calculation. For a reliable evaluation of the spectra it is generally enough to have one circle for each $`R_i`$, plus some more, about 5, to resolve the central region. The number of circles $`N_b`$ in each image is therefore roughly the same as the number of radial grid points: $`N_b=N_R+5`$. The number of pixels in each circle $`N_\phi `$ is slightly less straightforward to choose, but for reliable spectra it is generally sufficient to take $`N_\phi =2N_\mathrm{\Theta }`$, where $`N_\mathrm{\Theta }`$ is the number of $`\mathrm{\Theta }`$ grid points, counted from pole to pole. Using this method, the images automatically resolve all relevant scales, while using only a fairly limited number of pixels.

## 5 Testing the ESC Lambda Operator

The Extended Short Characteristic implementation of the Lambda Operator is not exact, as opposed to the one based on Long Characteristics. The interpolations used in the ESC algorithm introduce numerical diffusion, even in the optically thin regime, and this constitutes a potential threat to the reliability of the method. In order to test the accuracy of the ESC Lambda Operator we have performed a series of runs for a number of simple setups, comparing the results of the ESC calculation with those obtained by means of an exact LC Lambda Operator. Here we present the analysis for three such tests.

### 5.1 Optically thick annulus

The first test problem concerns the determination of the radiation field produced by an optically thick, isothermal, sharply–edged annuls, bounded by $`R_0<R<R_M`$ and $`\mathrm{\Theta }_M<\mathrm{\Theta }<\pi \mathrm{\Theta }_M`$. In the actual calculation we have taken $`R_0=2.86`$, $`R_M=4.83`$, $`\mathrm{\Theta }_M=1.26`$ and the absorbtion $`\alpha _a`$ is given by

$$\alpha _a(R,\mathrm{\Theta })=\{\begin{array}{cc}10^3& \text{inside the annulus}\\ 0& \text{elsewhere}.\end{array}$$
(45)
For the sake of simplicity all variables are dimensionless and the temperature has been taken such that $`B_\nu (T)=1`$. We use a spatial grid with $`20`$ radial points, logarithmically spaced such that

$$\frac{R_{i+1}R_i}{R_i}=1.1402,$$

and the angular grid in $`\mathrm{\Theta }`$ consists of 20 equally spaced points from pole to pole. Fig. 7

shows the configuration for the test problem. Mirror symmetry in the equator reduces the number of actual points to 10 in the range $`0<\mathrm{\Theta }<\pi /2`$. The mesh in photon momentum space consists of 32 equally spaced points in $`\varphi `$, covering the range 0–$`2\pi `$, and 41 points in $`\mu `$, 38 chosen according to Eq.(40), plus $`\mu =0`$ and $`\mu =\pm 0.24`$ to ensure sufficient resolution at small $`\mu `$. Because of the symmetry, the transfer needs to be solved only for the 16 points in the ranges $`0<\varphi \pi /2`$ and $`3\pi /2<\varphi 2\pi `$.

The radiation field emitted by such a source can be semi–analytically determined. As seen by an observer at some point $`P`$, it is simply the projection of the object on the sky of the observer. Since the object is sharply–edged and highly optically thick ($`\tau 100`$), its projection will be sharply–edged as well. The intensity is simply given by

$$I_\nu (P;\mu ,\varphi )=\{\begin{array}{cc}0& \text{for rays missing the object}\\ B_\nu (T)& \text{for rays hitting the object}\end{array},$$
(46)
so it is easy to compute the projections on the sky of the observer at various positions in space by using some independent ray–tracing algorithm or a semi–analytical computation of the image. This image can then be compared to that produced by the ESC and the LC transfer algorithms to evaluate the accuracy of the ESC Lambda Operator.

Let $`I_{ESC}(\mu ,\varphi )`$ and $`I_{LC}(\mu ,\varphi )`$ be the intensities, at the observer location, computed using the ESC and the LC algorithms, for the same discretization in $`\mu `$ and $`\varphi `$. Let, moreover, $`I(\mu ,\varphi )`$ be the true intensity, which may be found by tracing individual rays with very high resolution in $`\mu `$ and $`\varphi `$. We define the standard error of the ESC algorithm as

$$\sigma _{ESC}^2=\frac{(I_{ESC}I)^2𝑑\mu 𝑑\varphi }{I^2𝑑\mu 𝑑\varphi }$$
(47)
and the error in the mean intensity $`J`$ as

$$ϵ_{ESC}=\frac{(I_{ESC}I)𝑑\mu 𝑑\varphi }{I𝑑\mu 𝑑\varphi }\frac{J_{ESC}J}{J}.$$
(48)
Similar definitions apply to the LC errors.

Radiative transfer for this setup has been performed using the ESC method. The results are shown in Fig. 8 for the 9 gridpoints labeled $`AI`$ in Fig. 7.

The contours of the real images are overplotted. The same calculation has been repeated using the LC method. The errors of both the ESC and the LC calculations are listed in Table 1. These figures show that the errors of the ESC method are not very much greater than those of the LC method (which result from the discretization alone) and strenghten the reliability of the ESC algorithm.

### 5.2 Optically thin annulus

Contrary to the optically thick case of the previous subsection, the radiation field from an optically thin annulus cannot be determined by a simple analytic formula such as Eq.(46). At large radii, however, the mean intensity should follow the $`1/R^2`$ law and be independent of $`\mathrm{\Theta }`$. We can verify if the solution produced by the ESC algorithm indeed has this expected behavior. All details are the same as in the previous test, with the only difference that now the (constant) absorption is taken $`10^6`$.

We focus on the behavior at large radii. In radial streaming $`H(R_i,\mathrm{\Theta }_j)J(R_i,\mathrm{\Theta }_j)`$ which, for an optically thin source, is simply proportional to the volume integral of the emissivity $`j(R,\mathrm{\Theta })`$ and is independent of $`\mathrm{\Theta }`$

$$J(R)=\frac{1}{4\pi R^2}j(R^{},\mathrm{\Theta }^{})(R^{})^2𝑑R^{}\mathrm{sin}\mathrm{\Theta }^{}d\mathrm{\Theta }^{}.$$
(49)
Although close to the source the mean intensity is still very dependent on $`\mathrm{\Theta }`$, at $`RR_M`$ the code should be able to recover the correct behaviour $`JR^2`$ at large radii. The mean intensity $`J(R_i,\mathrm{\Theta }_J)`$ resulting from the ESC transfer calculation is shown in Fig. 9, where it is multiplied with $`R^2`$. It can be clearly seen that the while $`J`$ depends on $`\mathrm{\Theta }`$ for small radii, it becomes almost independent on $`\mathrm{\Theta }`$ as the radius increases and that the inverse square law is very well reproduced. The dependence of $`J`$ on $`\mathrm{\Theta }`$ at the largest radius is shown in Fig. 10. Typical errors are $`5\%`$, which lies within the error expected from the coarseness of the grid.

### 5.3 Spherically symmetric test problem

Although the above tests show that the ESC algorithm performs well on its own, they are not enough to prove that it will produce accurate and reliable results when applied to, for instance, a non–LTE line transfer computation. Unfortunately it is not easy to test this, because to our knowledge there exists no useful benchmark test case yet for 2–D axisymmetric radiative transfer in circumstellar clouds.

The least we can do is to test our 2–D algorithm on a 1–D spherically symmetric test case, and check the output against that produced by an independent 1–D transfer calculation. This way we can at least test two of the special features of the ESC algorithm: the additional $`\mu `$-points close to $`\mu 0`$, and the special treatment of the intensity near $`\mu 1`$. These are features that are not particularly related to 2-D, and can therefore also be tested in 1-D just as well.

Our test cloud is a spherically symmetric power law model with hydrogen density specified by

$$N_{\text{H}_2}(R)=N_{\text{H}_2}^0\left(\frac{R_0}{R}\right)^2\text{cm}^3,$$
(50)
where $`R`$ is the radius in cm, and $`N_{\text{H}_2}^0=2.0\times 10^7\text{cm}^3`$ is the number density at $`R=R_01.0\times 10^{15}\text{cm}`$. We take a constant kinetic temperature $`T_{\text{kin}}(R)=20\text{K}`$. The abundance of our molecule is also a constant, $`X_{\text{mol}}(R)N_{\text{mol}}(R)/N_{\text{H}_2}(R)=1.0\times 10^6`$. The systematic velocity is taken zero. The models are computed in spherical coordinates, in the domain $`R_{\text{in}}RR_{\text{out}}`$. We take $`R_{\text{in}}=1.0\times 10^{15}\text{cm}`$ and $`R_{\text{out}}=7.8\times 10^{18}\text{cm}`$. At the inner boundary we put a reflective boundary condition. The incoming radiation at the outer boundary is the $`T=2.728\text{K}`$ microwave background radiation.

We choose a fictive 2-level molecule which is specified by

$`E_2E_1`$ $`=6.0\text{cm}^1=8.63244\text{K}`$ (51)
$`g_2/g_1`$ $`=3.0`$ (52)
$`A_{21}`$ $`=1.0\times 10^4\text{s}^1`$ (53)
$`K_{21}`$ $`=2.0\times 10^{10}\text{cm}^3\text{s}^1`$ (54)
from which the downward collision rate follows: $`C_{21}=N_{\text{H}_2}K_{21}\text{s}^1`$. The total (thermal+turbulent) line width $`a_{\text{tot}}`$ is $`a_{\text{tot}}=0.150\text{km}\text{s}^1`$ (see Eq. 6).

The test problem presented here has high optical depth and a very sub-critical density. It is therefore well suited to test whether non-LTE effects are properly computed.

The line transfer is computed in a passband of 40 frequency points equally spaced between $`0.40\text{km}\text{s}^1`$ and $`+0.40\text{km}\text{s}^1`$. The $`\mu `$ angle is discretized using $`43`$ points, arranged according to the tangent-ray prescription of Eq.(42) with $`3`$ additional $`\mu `$-points around $`\mu =0`$ on each side. Our convergence criterion is simply:

$$\mathrm{max}(\delta n_i/n_i)<1\times 10^4$$
(55)
at all radii. For the radius we use an equally spaced logarithmic grid with $`(R_{i+1}R_i)/R_i\mathrm{\Delta }R/R=0.1720`$. We perform 4 runs: Lambda Iteration and Accelerated Lambda Iteration with and witout Ng acceleration.

The results for the upper level population is plotted in Fig. 11 and compared to the results obtained independently with SIMLINE written by V. Ossenkopf (Ossenkopf (1999)). The convergence plots for four different methods are shown in Fig. 12.

More test problems, and a more extensive discussion of them was presented by Dullemond (Dullemond (1999)).

## 6 A simple model for the Egg Nebula

Now that the algorithm has been tested, we demonstrate here how it can be used in practice. Our first example is a simple model of the optical appearance of the Cygnus Egg Nebula (CRL 2688). This object has been extensively studied ever since its discovery by Ney et al. (Ney et al. (1975)). It is a bipolar reflection nebula surrounding an F5 supergiant of $`T_{\text{eff}}=6500\text{K}`$ (Crampton Crampton et al. (1975)). At optical wavelength it appears as a diffuse bi–lobed nebula with two sharply edged “searchlight beams” emerging from each of the poles (Sahai et al. Sahai et al. (1998)). The lobes are separated by a dark equatorial lane which completely obscures the central star.

The optical emission from this nebula can be understood as reflected starlight escaping from the nebula through polar cavities (Latter et al.Latter et al. (1993), Morris Morris (1981)). It is clear that the “searchlight beams” are due to single scattering of direct starlight by dust grains. The lobes are, however, more likely to be the result of multiple–scattering.

We will model this multiple–scattering process in 2–D with the MESC algorithm, in an attempt to reproduce the complex optical appearance of the nebula. Our setup consists of an almost spherical wind with a cavity at both poles. The density in the cavities is small, but it is still high enough to reflect sufficient amounts of starlight. To reproduce the twin–beams at both poles, we place a small blob of matter at the polar axis in the cavity, causing a shadow. A star is placed at the center of the coordinate system to illuminate the nebula from within.

We model only a single frequency in the optical, at 600 nm. For this reason we refrain from taking actual realistic dust opacities, and specify total optical depth and scattering albedo instead. The dust density is shown as contour plots in Fig. 13. The total optical depth at the equator is about 60. The ratio of absorption over scattering is $`0.3`$. The dust scattering is assumed to be isotropic, which suffices for the present simplified example. We do not specify the dust temperature for this setup since thermal emission at 600 nm is negligible.

The simulation was performed by RADICAL, using the MESC algorithm, and applying Accelerated Lambda Iteration and Ng acceleration. The image subsequently produced by formal integration is shown in Fig. 14. The model reproduces the searchlight beams and the diffuse glow. It also naturally reproduces the intensity difference between the north and south lobe. This is a result of the slight inclination at which the object is seen. For light emerging from the south lobe, the path length through the outer regions of the nebula is larger than for the north lobe.

Although the model resemblances the HST image of Sahai et al., it should be noted that the density structure that we have used may not be consistent with observations at other wavelengths. For example, HCN observations by Bieging & Nguyen-Q-Rieu (Bieging and Q.-Rieu (1996)) seem to rule out the presence of a cavity in the wind. Also the rather high albedo may be difficult to reconcile with the fact CRL 2688 is carbon–rich.

## 7 A model of line transfer in a collapsing protostellar cloud

In this section it will be demonstrated how the 2-D transfer algorithm can be used in the observational study of low-mass star formation in dense molecular cloud cores.

Low mass star formation takes place in dense molecular cloud cores. According to the spherically symmetric model of Shu (Shu (1977)), such a core develops a cusp with density close to a $`1/R^2`$ powerlaw. Once a gravitational instability is trigged, the centeral part of the cloud collapses, and forms a star. An expansion wave propagates into the cloud towards larger radii, allowing more and more matter to fall supersonically down the potential well, and add to the protostar’s mass.

Both observational evidence and theoretical arguments, however, indicate that purely spherical collapse is rare. Any slight amount of angular momentum in the primordial cloud will cause deviations from sphericity as centrifugal forces tend to dominate over radial infall deep down the potential well. And even before the collapse stage these primordial clouds often appear to be non-spherical (Myers et al. Myers et al. (1991)). Theoretical models of non-spherical protostellar collapse include, among others, Ulrich (Ulrich (1976)), Cassen & Moosman (Cassen and Moosman (1981)), Terebey et al. (Terebey et al. (1984), Galli & Shu Galli and Shu (1993).

The models of Cassen & Moosman and Ulrich (hereafter CMU) focus on the inner free-falling part of the collapsing cloud, and assume that the material originates from an originally spherical cloud with some angular momentum. Their model is almost spherical at large radii, but flattens off closer towards the center, and forms a disk near the centrifugal radius. This model was later extended by Hartmann et al. (Hartmann et al. (1996), hereafter HCB) to include flattening of the parent cloud. These models show that the inner free-fall part of an initially flattened cloud naturally tends to form a bipolar cavity, which is often observed in YSO. These models are distinctly non-spherical at all radii, despite the fact that centrifugal forces only dominate at small radii. It is therefore evident that fitting such models to molecular line observations requires 2-D axi-symmetric radiative transfer computations.

In this section we perform such a calculation, using the algorithms of this paper. We solve the non-LTE level populations for the first 7 rotational levels of HCO+, and compute the predicted single–dish spectra.

### 7.1 Description of the model

In our HCB models we assume that the radius of the expansion wave $`R_{\text{acc}}`$ is outside our domain, so we shall confine our study to the free-fall inner region of the collapsing sheet-like molecular cloud. We assume that matter in the parent cloud had a small amount of rotation in the plane of the sheet before it collapsed. According to the HCB model, the velocity field of the gas is given by the formulae of Ulrich (Ulrich (1976)):

$`v_R`$ $`=\left({\displaystyle \frac{GM}{R}}\right)^{1/2}\left(1+{\displaystyle \frac{\mu }{\mu _0}}\right)^{1/2}`$ (56)
$`v_\mathrm{\Theta }`$ $`=\left({\displaystyle \frac{GM}{R}}\right)^{1/2}\left({\displaystyle \frac{\mu _0\mu }{\sqrt{1\mu ^2}}}\right)\left(1+{\displaystyle \frac{\mu }{\mu _0}}\right)^{1/2}`$ (57)
$`v_\varphi `$ $`=\left({\displaystyle \frac{GM}{R}}\right)^{1/2}\left({\displaystyle \frac{1\mu _0^2}{1\mu ^2}}\right)^{1/2}\left(1{\displaystyle \frac{\mu }{\mu _0}}\right)^{1/2},`$ (58)
where $`\mu \mathrm{cos}(\mathrm{\Theta })`$ and $`\mu _0\mathrm{cos}(\mathrm{\Theta }_0)`$. The angle $`\mathrm{\Theta }_0`$ is the $`\mathrm{\Theta }`$-coordinate that the gas parcel had when it started its free-fall at large radius. For a given $`(R,\mathrm{\Theta })`$ the value of $`\mathrm{\Theta }_0`$ can be found from

$$\frac{R}{R_\text{c}}\left(1\frac{\mu }{\mu _0}\right)=(1\mu _0^2),$$
(59)
where $`R_\text{c}`$ is the centrifugal radius, i.e. the radius at which centrifugal forces equal gravity. This is the outer radius of the disk that is formed as a result of the rotation.

The density of the gas for the HCB model is given by

$$\begin{array}{cc}\hfill \rho _{HCB}(R,\mathrm{\Theta })=& \frac{\dot{M}}{4\pi (GMR^3)^{1/2}}\frac{\eta \text{sech}^2(\eta \mu _0)}{\mathrm{tanh}(\eta )}\hfill \\ & \times \left(1+\frac{\mu }{\mu _0}\right)^{1/2}\left(\frac{\mu }{\mu _0}+\frac{2\mu _0^2R_\text{c}}{R}\right)^1\hfill \end{array}$$
(60)
where $`\eta `$ is a dimensionless flattening parameter, which is roughly equal to the ratio of the accretion radius $`R_{\text{acc}}`$ to the sheet thickness $`H`$. HCB argue that this value must be somewhere in between $`\eta =0`$ and $`\eta =\eta _{\text{max}}4`$. For $`\eta =0`$ the CMU models are reproduced. Density contours of this free-falling envelope, for different flattening parameters, are shown in Fig. 15.

The centrifugal radius of the infalling envelope is at $`R_\text{c}=100\text{AU}`$. We place a thin disk with a radius of $`R_{\text{disk}}=R_\text{c}=100\text{AU}`$ at the equator. A zoom-in of the density distribution, down to the scale of the disk, is shown in Fig. 16. We will ignore any emission from the disk, and merely treat it as a light-blocking boundary condition at the equator.

### 7.2 Non-LTE line transfer

We compute the line transfer problem for four HCB models, with flattening parameter $`\eta =0,1,2,3`$ for models 1,2,3,4 respectively. The adopted valued for the accretion rate and the turbulent line width are $`\dot{M}=2.4\times 10^6M_{}\text{yr}^1`$ and $`a_{\text{turb}}=0.25\text{km}\text{s}^1`$.

For these models we compute the non-LTE line transfer problem of the lowest few rotational levels of HCO+, including the effects of the moving medium. We assume a cosmic background (CMB) continuum as incident radiation at the outer edge of the computational domain. Dust emission and opacity are neglected in the line transfer, which is justified for the lower-lying HCO+ lines because the nebula is optically thin to dust in the millimeter and sub-millimeter, and radiative pumping by dust continuum is not important for HCO+. Also, we need not include the dust emissivity in the computation of the emerging spectra, since we shall show only the spectra with the dust- and CMB-continuum removed. The radiative transfer is computed within a range of $`R[70,10^4]\text{AU}`$. As an inner boundary we have a vacuum. The cross sections for H2-HCO+ collisional transitions were taken from Monteiro (Monteiro (1985)) and Green (Green (1975)). We adopt an HCO+ abundance of $`2\times 10^9`$ The gas temperature is taken to be $`T=20\text{K}`$ throughout the cloud.

We perform the non-LTE line transfer for all four models. The resulting non-LTE level populations for model 3 are shown in Fig. 17, and the corresponding excitation temperatures are shown in Fig. 18. One can see that at the equator ($`\mathrm{\Theta }=\pi /2`$) the levels are almost thermalized, except at large radii. This is due to the much larger density at the equator than at the pole. The drop in excitation temperature at large radii is a result of the decoupling of radiation and matter. At large radii the level populations will be strongly influenced by the cosmic background radiation. Another interesting phenomenon occurs at small radii near the pole: the excitation temperature exceeds the gas temperature. This effect was discussed by Leung & Liszt (Leung and Liszt (1976)) for the CO 1-0 transition. It can be understood as resulting from an overpopulation of the $`J=1`$ level due to the large ratio of radiative rates ($`A_{21}/A_{10}10`$).

Once the level populations have been computed, the line spectra are produced. The spectra are centered on the origin of the object. First the circular images are produces in a range of frequencies. This circular rendering of the images ensures that no details at large or small radii are missed, and thus that no flux is accidently lost. The antenna temperatures are then computed by integrating the images, after they have been multiplied by the beam pattern centered on the origin of the object. We use an “Airy” beam, with a beam size corresponding to a single dish of 15 m diameter. The object is placed at 140 parsec distance. The spectra of the four models, computed for the first four radiative transitions at three different inclinations, are shown in Fig. 19.

From the spectra one can clearly see the effect of flattening of the HCB cloud, in particular for models 3 and 4. At near pole-on inclination (5) hardly any self-absorption is seen in these models, because one looks straight into the “cavity”. At near edge-on inclination (85) the “torus” blocks the central regions from view near line center, which results in the clear self-absorption features seen in the line shapes. A similar manifestation of the non-spherical symmetry of the cloud has been discussed recently by Van der Tak et al. van der Tak et al. (1999). An interesting feature of the line spectra of models 3 and 4 is that the edge-on line profiles are wider than the pole-on profiles. This can be attributed to the fact that the density and the excitation temperature is lower at the pole than at the equator. At pole-on inclination the high density equatorial matter will emit near line-center instead of in the line wings, thus making the line profile narrower.

The asymmetry between the red-shifted and the blue-shifted peaks are typical for protostellar collapse. The rotation is hardly seen in these spectra. This is because the rotational velocity is everywhere much smaller than the free-fall velocity, except at very small radii where the emission barely contributes to the single-dish spectra shown here.

## 8 Conclusions

Numerical radiative transfer modeling on desktop workstations is extremely cheap. Not doing so in cases where this is possible would mean an enormous waste of valuable information that lies encoded in observed data. However the success of such modeling depends on the algorithms that are available. We have developed a robust and accurate method, called the “extended short characteristics” (ESC) method, by which complicated 2–D axi–symmetric multi–frequency radiative transfer calculations can be performed. By using spherical coordinates, this method can accurately treat circumstellar envelopes and disks from the stellar surface all the way up to parsec scale, without the need of grid refinement. By making a special choice of discrete photon angles and bundling ‘almost–radially–moving’ rays into a single bin, the conservation of radial flux can be guaranteed even over many orders of magnitude in radius, and without excessive computational cost.

The ESC method, and a slight variation called MESC, forms the core of a multi–purpose 2-D radiative transfer code called RADICAL. We have tested the ESC/MESC algorithm on a simple test problem which we described in this paper. We have also verified that the 2-D algorithm, when applied to a 1-D spherically symmetric problem, indeed reproduces what an independent 1-D algorithm would produce for the same problem. The errors remained within a few percent, the exact value of which depends on the grid resolution.

The ESC/MESC algorithm is designed for a variety of applications. We have demonstrated in this paper how the method can be used for the problem of dust scattering in a bipolar proto–planetary nebula and the problem of non-LTE line transfer in a collapsing cloud. But the method can also easily be applied to other radiative processes, such as dust continuum emission with radiative equilibrium for the dust grains, thermal Bremsstrahlung, electron scattering and even Comptonization in hot plasmas. The ESC/MESC algorithms are therefore a useful tool for a wide range of astrophysical problems.

## Acknowledgments

We wish to thank Vincent Icke for making us aware of the applicability of our algorithm (which was originally intended only for Comptonization) to non–LTE line processes and dust absoption/emission in stellar nebulae and winds, and for suggesting the Egg Nebula application. We are grateful to Volker Ossenkopf for sending us his solution to the test problem. CPD thanks Gerd–Jan van Zadelhof for useful discussions on line transfer, Ewine van Dishoeck, Rens Waters and Alex de Koter for discussions on the applicability of the algorithm to envelopes around young and old stars. CPD thanks Björn Heijligers and Vincent de Heij for their help with IDL. Special thanks to Jeremy Yates, Floris v.d. Tak and Volker Ossenkopf for their collaboration in comparing the results of RADICAL with those of their codes for molecular line transfer.