© ESO, 2012

1. Introduction

In a series of papers (Hauschildt & Baron 2006, 2008, 2009, 2010; Baron & Hauschildt 2007; Baron et al. 2009b; Seelmann et al. 2010) we have developed a general framework for solving 3D radiative transfer problems in Cartesian, cylindrical, and spherical coordinates for both static and monotonic velocity fields in the comoving frame. We have also developed an Eulerian code for velocities ${\mathrm{~}}_{}\mathit{<}$1000 km s^-1 (Seelmann et al. 2010). The neglect of relativistic effects and resolution constraints limits the applicability of the Eulerian approach to v/c ≲ 0.01. Our affine method for solving the fully relativistic transfer equation is exact to all orders in v/c (Chen et al. 2007; Baron et al. 2009a). In a parallel series of papers we have introduced methods for solving the fully relativistic transfer equation in 1D spherical coordinates with arbitrary velocity fields (Baron & Hauschildt 2004; Knop et al. 2009a,b).

2. Basic formalism

Using the general formalism developed in Chen et al. (2007) we can derive the transfer equation in flat spacetime with arbitrary flows. We choose to work in spherical coordinates without loss of generality.

The photon’s four-momentum can be written ${\mathit{p}}^{\mathit{a}}\mathrm{\equiv }\frac{\mathrm{d}{\mathit{x}}^{\mathit{a}}}{\mathrm{d}\mathit{\xi }}\mathrm{=}\frac{\mathit{h}}{{\mathit{\lambda }}_{\mathrm{\infty }}}\mathrm{\left(}\mathrm{1}\mathit{,}\hat{{n}}\mathrm{\right)}\mathit{,}$(1)where h is Planck’s constant, ξ is the affine parameter, λ_∞ is the rest frame wavelength, and $\hat{{n}}$ is the 3D direction of the photon as seen by a distant, stationary, observer. The four-velocity of the comoving observer in an arbitrary flow can be written ${\mathit{u}}^{\mathit{a}}\mathrm{=}\mathit{\gamma }\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{[}\mathrm{1}\mathit{,}\hspace{0.17em}{\beta }\mathrm{(}{r}\mathit{,t}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(2)and the comoving wavelength λ can be obtained using $\frac{\mathit{h}}{\mathit{\lambda }}\mathrm{=}\mathrm{-}\mathrm{(}\mathit{u}\mathrm{·}\mathit{p}\mathrm{)}\mathit{.}$(3)Here, β = v/c, and γ = (1 − β²) ^− 1/2, are the usual quantities of special relativity. The 3D geodesic in the flat spacetime can be parametrized as ${r}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{{r}}_{\mathrm{0}}\mathrm{+}\hat{{n}}\hspace{0.17em}\mathit{s,}$(4)where r₀ is the starting point of the characteristics, and s is the rest frame physical distance related to the affine parameter ξ by $\mathit{s}\mathrm{\equiv }\frac{\mathit{h}}{{\mathit{\lambda }}_{\mathrm{\infty }}}\mathit{\xi }\mathit{.}$(5)The radiative transfer equation can be written in terms of the affine parameter ξ as (see Eq. (10) of Chen et al. 2007): ${\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \xi }}|}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{\xi }}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \lambda }}\mathrm{=}\mathrm{-}\left({\mathit{\chi }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{5}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{\xi }}\right){\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathit{,}$(6)where ${\mathit{I}}_{\mathit{\lambda }}\mathrm{\left(}{r}\mathit{,t}\mathrm{;}\hat{{n}}\mathrm{\right)}$ is the specific intensity measured by a comoving observer (note that I_λλ⁵ is observer-independent) at the (global) rest frame space-time point (r, t), toward the rest frame direction $\hat{{n}}$, and at the comoving wavelength λ. When expressing the 7D phase-space dependence of the comoving specific intensity, the only comoving variable we used was the wavelength λ, in particular, we did not use angles measured by a comoving observer to specify the direction of the photons. The advantages for this at first sight odd phase-space configuration have been explained in detail in Chen et al. (2007), and we applied this technique in Baron et al. (2009a).

We can rewrite Eq. (6) as $\begin{array}{ccc}{\frac{\mathrm{d}\mathrm{\left(}\mathit{ct}\mathrm{\right)}}{\mathrm{d}\mathit{\xi }}\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial t}}|}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{d}{r}}{\mathrm{d}\mathit{\xi }}\mathrm{·}\mathrm{\nabla }{\mathit{I}}_{\mathit{\lambda }}& \mathrm{+}& \frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{\xi }}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \lambda }}\mathrm{=}\\ & & \mathrm{-}\left({\mathit{\chi }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{5}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{\xi }}\right){\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathit{,}\end{array}$(7)or equivalently $\begin{array}{ccc}{\frac{\mathrm{d}\mathrm{\left(}\mathit{ct}\mathrm{\right)}}{\mathrm{d}\mathit{\xi }}\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial t}}|}_{\mathit{\lambda }}& \mathrm{+}& \frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{\xi }}\frac{\mathrm{d}{r}}{\mathrm{d}\mathit{s}}\mathrm{·}\mathrm{\nabla }{\mathit{I}}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{\xi }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{s}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \lambda }}\mathrm{=}\\ & & \mathrm{-}\left({\mathit{\chi }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{5}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{\xi }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{s}}\right){\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\frac{\mathit{h}}{\mathit{\lambda }}\mathrm{·}\end{array}$(8)Then using the definition of s from Eq. (5) and the fact that $\frac{\mathrm{d}\mathrm{\left(}\mathit{ct}\mathrm{\right)}}{\mathrm{d}\mathit{\xi }}\mathrm{=}\mathit{c}{\mathit{p}}^{\mathit{t}}\mathrm{=}\frac{\mathit{h}}{{\mathit{\lambda }}_{\mathrm{\infty }}}\mathit{,}$(9)we find $\begin{array}{ccc}{\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial t}}|}_{\mathit{\lambda }}\mathrm{+}{\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial s}}|}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{s}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \lambda }}& \mathrm{=}& \mathrm{-}\left({\mathit{\chi }}_{\mathit{\lambda }}\frac{{\mathit{\lambda }}_{\mathrm{\infty }}}{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{5}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{s}}\right){\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\frac{{\mathit{\lambda }}_{\mathrm{\infty }}}{\mathit{\lambda }}\mathrm{·}\end{array}$(10)Equation (10) can be put into our standard form: ${\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial t}}|}_{\mathit{\lambda }}\mathrm{+}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial s}}\mathrm{+}\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}\frac{\mathit{\partial }}{\mathit{\partial \lambda }}\mathrm{\left(}\mathit{\lambda }{\mathit{I}}_{\mathit{\lambda }}\mathrm{\right)}\mathrm{+}\mathrm{4}\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{I}}_{\mathit{\lambda }}\mathrm{=}\mathrm{-}{\mathit{\chi }}_{\mathit{\lambda }}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathit{,}$(11)where $\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{\equiv }\frac{{\mathit{\lambda }}_{\mathrm{\infty }}}{\mathit{\lambda }}\mathrm{=}\mathit{\gamma }\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{[}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }\mathrm{(}{r}\mathit{,t}\mathrm{\left)}\mathrm{\right]}$(12)is simply the Doppler factor, and $\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{1}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}\mathit{s}}\mathrm{·}$(13)Using Eqs. (4) and (12), a(s) is found to be $\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }}\left[\frac{\mathrm{d}}{\mathrm{d}\mathit{s}}\mathrm{(}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\frac{\mathrm{d}\mathit{\beta }}{\mathrm{d}\mathit{s}}\right]\mathit{,}$(14)where β is the magnitude of β, and $\frac{\mathrm{d}}{\mathrm{d}\mathit{s}}\mathrm{=}\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{+}\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{=}\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial s}}\mathrm{·}$(15)When we numerically integrate the radiation transfer equation, a(s) can be approximated as $\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{\approx }\frac{\mathit{\delta }\mathrm{(}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\hspace{0.17em}\mathit{\delta \beta }}{\mathit{\delta s}\mathrm{(}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}}\mathit{,}$(16)where δs is the differential step size (physical distance) along the characteristics, $\mathit{\delta }\mathrm{(}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}$ and δβ are the changes of $\hat{{n}}\mathrm{·}{\beta }$ and β, respectively, when we move one step forward, which includes the changes induced by both time and spatial advances, for instance $\mathit{\delta \beta }\mathrm{=}\mathit{\beta }\mathrm{\left(}{\mathit{s}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathit{,}{\mathit{t}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathrm{-}\mathit{\beta }\mathrm{\left(}{\mathit{s}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(17)Since few numerical schemes will be able to provide the fully implicit derivative, δβ, will often be obtained for example by the backward difference $\mathit{\delta \beta }\mathrm{=}\mathrm{\left(}\mathit{\beta }\mathrm{\right(}{\mathit{s}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathit{\beta }\mathrm{(}{\mathit{s}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{)}\mathrm{+}\mathrm{(}\mathit{\beta }\mathrm{\left(}{\mathit{s}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}\mathit{\beta }\mathrm{\left(}{\mathit{s}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathit{,}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathit{.}$(18)In the stationary case, both β and f(s) are independent of time and specializing Eqs. (4) and (12) to that case, a(s) becomes $\begin{array}{ccc}\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}& \mathrm{=}& \frac{\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\mathrm{(}\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathit{\beta }\mathrm{)}}{\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{·}{\beta }}\\ & \mathrm{=}& \end{array}$(19)where we have used the fact that along the characteristics, d/ds no longer contains the time derivative and is thus the directional derivative operator, that is, d/d$\mathit{s}\mathrm{=}\hat{{n}}\mathrm{·}\mathrm{\nabla }$. We recall that in the flat spacetime that we are considering, our characteristics are straight lines for all velocity flows.

In terms of its spherical components, β can be written ${\beta }\mathrm{=}{\mathit{\beta }}_{\mathit{r}}\hat{{e}}\mathit{r}\mathrm{+}{\mathit{\beta }}_{\mathit{\theta }}\hat{{e}}\mathit{\theta }\mathrm{+}{\mathit{\beta }}_{\mathit{\phi }}\hat{{e}}\mathit{\phi }\mathit{,}$(20)where $\hat{{e}}{}_{\mathit{r}}\mathit{,}\hat{{e}}{}_{\mathit{\theta }}\mathit{,}{\hat{{e}}}_{\mathit{\phi }}$ are the spherical orthonormal basis vectors at point r(r,θ,φ), i.e. $\begin{array}{ccc}\hat{{e}}\mathit{r}& \mathrm{=}& \mathrm{\left(}\sin \mathit{\theta }\cos \mathit{\phi ,}\sin \mathit{\theta }\sin \mathit{\phi ,}\cos \mathit{\theta }\mathrm{\right)}\mathit{,}\\ \hat{{e}}\mathit{\theta }& \mathrm{=}& \mathrm{(}\cos \mathit{\theta }\cos \mathit{\phi ,}\cos \mathit{\theta }\sin \mathit{\phi ,}\mathrm{-}\sin \mathit{\theta }\mathrm{)}\mathit{,}\\ {\hat{{e}}}_{\mathit{\phi }}& \mathrm{=}& \mathrm{(}\mathrm{-}\sin \mathit{\phi ,}\cos \mathit{\phi ,}\mathrm{0}\mathrm{)}\mathit{,}\end{array}$(21)and consequently the $\hat{{n}}\mathrm{·}{\beta }$ in Eq. (16) can be calculated using $\begin{array}{ccc}\hat{{n}}\mathrm{·}{\beta }& \mathrm{=}& {\mathit{\beta }}_{\mathit{r}}\hat{{n}}\mathrm{·}\hat{{e}}\mathit{r}\mathrm{+}{\mathit{\beta }}_{\mathit{\theta }}\hat{{n}}\mathrm{·}\hat{{e}}\mathit{\theta }\mathrm{+}{\mathit{\beta }}_{\mathit{\phi }}\hat{{n}}\mathrm{·}\hat{{e}}\mathit{\phi }\\ & \mathrm{\equiv }& {\mathit{\beta }}_{\mathit{r}}{\mathit{n}}_{\mathit{r}}\mathrm{+}{\mathit{\beta }}_{\mathit{\theta }}{\mathit{n}}_{\mathit{\theta }}\mathrm{+}{\mathit{\beta }}_{\mathit{\phi }}{\mathit{n}}_{\mathit{\phi }}\end{array}$(22)(note that along the characteristics, $\hat{{n}}$ has constant Cartesian components, n_x, n_y, n_z, but changing spherical components, n_r, n_θ, n_φ). Writing $\hat{{n}}\mathrm{=}\mathrm{\left(}{\mathit{n}}_{\mathit{x}}\mathit{,}{\mathit{n}}_{\mathit{y}}\mathit{,}{\mathit{n}}_{\mathit{z}}\mathrm{\right)}\mathit{,}$ the spherical components n_r, n_θ,n_φ can be easily computed using Eq. (21).

2.1. Comparison with Mihalas

At first glance, comparing Eq. (11) with Eq. (2.12) of Mihalas (1980) something seems amiss. Like Mihalas (1980), we work in the frame where spatial coordinates and clocks are measured by an observer at rest. However, Mihalas’ time derivative contains a Doppler factor, whereas ours does not. Also, our terms on the right-hand side contain Doppler factors, f(s), whereas those of Mihalas do not. The discrepency has been noted in passing by Chen et al. (2007) and arises because our s is a true distance measured in the observer’s frame, whereas that of Mihalas, s_M, contains an extra Doppler factor:

$\mathrm{d}{\mathit{s}}_{\mathrm{M}}\mathrm{=}\frac{{\mathit{\lambda }}_{\mathrm{\infty }}}{\mathit{\lambda }}\mathrm{d}\mathit{s}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{d}\mathit{s}$(note that f(s) is not a constant along the characteristics, and therefore s_M is not related to the physical distance s by a simple affine transformation). Thus, we can transform from s to s_M in Eq. (10) to find ${\frac{\mathit{\lambda }}{{\mathit{\lambda }}_{\mathrm{\infty }}}\frac{\mathrm{1}}{\mathit{c}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial t}}|}_{\mathit{\lambda }}\mathrm{+}{\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial }{\mathit{s}}_{\mathrm{M}}}|}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}{\mathit{s}}_{\mathrm{M}}}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\lambda }}}{\mathit{\partial \lambda }}\mathrm{=}\mathrm{-}\left({\mathit{\chi }}_{\mathit{\lambda }}\mathrm{+}\frac{\mathrm{5}}{\mathit{\lambda }}\frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}{\mathit{s}}_{\mathrm{M}}}\right){\mathit{I}}_{\mathit{\lambda }}\mathrm{+}{\mathit{\eta }}_{\mathit{\lambda }}\mathit{,}$(23)which is very similar to the equation of Mihalas, except that the coefficient multiplying the time derivative term is the inverse Doppler factor f(s)^-1, since we are working with I_λ instead of I_ν, as did Mihalas. Jack et al. (2009) accidentally forgot to convert the time derivative from I_ν to I_λ and thus their Eqs. (24), (25) have a coefficient of the time derivative with the Doppler factor in the numerator, rather than in the denominator. Thus, the time derivative terms in their Eqs. (24), (25) should be multiplied by f(s)^-2 to give the correct equations.

2.2. Nonmonotonic flows

We are now in a position to tie together the work of Baron & Hauschildt (2004), Knop et al. (2009a) and Baron et al. (2009b). The formal solution in the monotonic case is an initial value problem in wavelength, but in the arbitrary flow case both spatial coordinates and wavelengths are fully coupled. This poses a significant memory cost, since the matrix obtained by finite-differencing the equations now contains every wavelength and not just the spatial points along the characteristic. The computational cost is surprisingly low because the linear system can be solved using the semi-analytic method of Knop et al. (2009a). While the framework that was given in Baron & Hauschildt (2004) and Baron et al. (2009b) was formulated for just one wavelength discretization, we included here the fully implicit discretization developed in Hauschildt & Baron (2004). Furthermore, we used the new formal solution that avoids negative generalized opacities (Knop et al. 2009b).

The stationary equation of radiative transfer in its characteristic form for the specific intensity I along a path s reads $\frac{\mathrm{d}{\mathit{I}}_{\mathit{l}}}{\mathrm{d}\mathit{s}}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{\eta }}_{\mathit{l}}\mathrm{-}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{\chi }}_{\mathit{l}}{\mathit{I}}_{\mathit{l}}\mathrm{-}\mathrm{4}{\mathit{a}}_{\mathit{l}}{\mathit{I}}_{\mathit{l}}\mathrm{-}{\mathit{a}}_{\mathit{l}}\frac{\mathit{\partial }\mathrm{\left(}\mathit{\lambda }{\mathit{I}}_{\mathit{l}}\mathrm{\right)}}{\mathit{\partial \lambda }}\mathit{,}$(24)where η is the emissivity, χ the opacity, and the subscript l indicates dependence on wavelength. The $\frac{\mathit{\partial }}{\mathit{\partial \lambda }}$-term is the coupling term between the wavelengths and depends on the structure of the atmosphere and on the mechanism of the coupling (Mihalas 1980).

The wavelength derivative can be discretized in two ways as described in Hauschildt & Baron (2004). The different discretizations can be mixed via a Crank-Nicholson-like scheme with a mixing parameter ξ ∈ [ 0,1 ] . The wavelength discretized equation of radiative transfer can then be written as $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{I}}_{\mathit{l}}}{\mathrm{d}\mathit{s}}& \mathrm{=}& \mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{\eta }}_{\mathit{l}}\mathrm{-}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{\chi }}_{\mathit{l}}{\mathit{I}}_{\mathit{l}}\mathrm{-}{\mathit{a}}_{\mathit{l}}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\xi }{{\mathit{p}}_{\mathit{l}}^{\mathrm{|}}}^{\mathrm{)}}{\mathit{I}}_{\mathit{l}}\\ & & \mathrm{-}\mathit{\xi }{\mathit{a}}_{\mathit{l}}\mathrm{(}{\mathit{p}}_{\mathit{l}}^{\mathrm{-}}{\mathit{I}}_{\mathit{l}\mathrm{-}\mathrm{1}}\mathrm{+}{\mathit{p}}_{\mathit{l}}^{\mathrm{+}}{\mathit{I}}_{\mathit{l}\mathrm{+}\mathrm{1}}\mathrm{)}\\ & & \mathrm{-}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{]}{\mathit{a}}_{\mathit{l}}\mathrm{(}{\mathit{p}}_{\mathit{l}}^{\mathrm{-}}{\mathit{I}}_{\mathit{l}\mathrm{-}\mathrm{1}}\mathrm{+}{\mathit{p}}_{\mathit{l}}^{\mathrm{|}}{\mathit{I}}_{\mathit{l}}\mathrm{+}{\mathit{p}}_{\mathit{l}}^{\mathrm{+}}{\mathit{I}}_{\mathit{l}\mathrm{+}\mathrm{1}}\mathrm{)}\mathit{,}\end{array}$(25)where the ${\mathit{p}}_{\mathit{l}}^{\mathrm{•}}$ coefficients in an ordered wavelength grid λ_l − 1 < λ_l < λ_l + 1 are defined as $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}& \mathrm{for} & {\mathit{a}}_{{\mathit{\lambda }}_{\mathit{l}}}\mathrm{\ge }\mathrm{0}\\ \begin{array}{c}\\ & & \\ & & \\ & & \end{array}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}& \mathrm{for} & {\mathit{a}}_{{\mathit{\lambda }}_{\mathit{l}}}\mathit{<}\mathrm{0.}\end{array}$The dependence on the sign of a_λ is introduced to define local upwind schemes (see Baron & Hauschildt 2004).

After introducing a generalized opacity (see Knop et al. 2009a,b) $\mathit{\chi ̂}\mathit{l}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{s}\mathrm{\right)}{\mathit{\chi }}_{\mathit{l}}\mathrm{+}\mathit{\xi }{\mathit{a}}_{\mathit{l}}{\mathit{p}}_{\mathit{l}}^{\mathrm{|}}\mathit{,}$(28)defining the source functions $\begin{array}{ccc}{\mathit{S}}_{\mathit{l}}& \mathrm{=}& \frac{{\mathit{\eta }}_{\mathit{l}}}{{\mathit{\chi }}_{\mathit{l}}}\\ \mathit{Ŝ}\mathit{l}& \mathrm{=}& \\ \begin{array}{}\mathrm{˜}\\ {\mathit{S}}_{\mathit{l}}\end{array}& \mathrm{=}& \end{array}$a formal solution of the radiative transfer problem can be formulated. We used a full characteristic method throughout the atmosphere. The spatial position on a characteristic is then discretized on a spatial grid. In the following a pair of subscript indices will mark the position in the spatial grid and in the wavelength grid. Commonly the spatial grid is mapped locally onto an optical depth grid via the relation . The formal solution of the radiative transfer Eq. (25) between two points s_i − 1 and s_i on a spatial grid along the photon path can be written in terms of the optical depth as follows: $\begin{array}{ccc}{\mathit{I}}_{\mathit{i,l}}& \mathrm{=}& {\mathit{I}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}}{\mathrm{e}}^{\mathrm{-}\mathrm{\Delta }\mathit{\tau }}\mathrm{+}\mathit{\delta }\mathit{Î}\mathit{i,l}\mathrm{+}\mathit{\delta }\begin{array}{}\mathrm{˜}\\ {\mathit{I}}_{\mathit{i,l}}\end{array}\\ \mathit{\delta }\mathit{Î}\mathit{i,l}& \mathrm{=}& \\ \mathit{\delta }\begin{array}{}\mathrm{˜}\\ {\mathit{I}}_{\mathit{i,l}}\end{array}& \mathrm{=}& \end{array}$with Δτ = τ_i + 1,l − τ_i,l and ${\mathit{\tau }}_{\mathit{i,l}}\mathrm{=}{{}^{\mathrm{\int }}\mathit{s}_{\mathit{i}}}^{}\mathit{\chi ̂}{}_{\mathit{l}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{d}\mathit{s}$. The α-β-γ coefficients are described in Olson & Kunasz (1987) and Hauschildt (1992). $\mathit{\delta }\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{I}}_{\mathit{l}}\end{array}$ in Eq. (34) is linearly interpolated and in general different from the coefficients in Eq. (33) and is therefore marked with a tilde.

Equation (32) can be written in matrix notation for any given characteristic: $\mathbf{I}\mathrm{=}{A}\mathrm{·}{I}\mathrm{+}\mathbf{\Delta }\mathit{I}\mathit{.}$(35)Here I is a vector with all intensities, A is a square matrix that describes the influence of the different intensities upon each other, and ΔI is a vector with the thermal emission and scattering contribution of the source function. For a characteristic with n_i spatial points and n_l points in the wavelength grid, the intensity vector I has n_i × n_l entries. In the following a superscript of k will label the characteristic at hand. The components of the matrix A from Eq. (35) at the spatial point i and the wavelength point l are given by $\begin{array}{ccc}{\mathit{A}}_{\mathit{i,l}}^{\mathrm{-}\mathit{,k}}& \mathrm{=}& \\ {\mathit{B}}_{\mathit{i,l}}^{\mathrm{-}\mathit{,k}}& \mathrm{=}& \mathrm{-}\mathrm{(}\mathit{\xi }{\mathit{\beta }}_{\mathit{i,l}}^{\mathit{k}}\mathrm{+}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{]}\hspace{0.17em}\stackrel{}{{{\mathit{\beta }}_{\mathit{i,l}}^{\mathit{k}}}_{\mathrm{˜}}}\mathrm{)}\frac{{\mathit{a}}_{\mathit{i,l}}^{\mathit{k}}}{\mathit{\chi ̂}\begin{array}{}\mathit{k}\\ \mathit{i,l}\end{array}}{\mathit{p}}_{\mathit{i,l}}^{\mathrm{-}\mathit{,k}}\\ {\mathit{C}}_{\mathit{i,l}}^{\mathrm{-}\mathit{,k}}& \mathrm{=}& \mathrm{-}\mathit{\xi }{\mathit{\gamma }}_{\mathit{i,l}}^{\mathit{k}}\frac{{\mathit{a}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,l}}^{\mathit{k}}}{\mathit{\chi ̂}\begin{array}{}\mathit{k}\\ \mathit{i}\mathrm{+}\mathrm{1}\mathit{,l}\end{array}}{\mathit{p}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,l}}^{\mathrm{-}\mathit{,k}}\\ {\mathit{A}}_{\mathit{i,l}}^{\setminus \mathit{,k}}& \mathrm{=}& \exp \mathrm{(}\mathrm{-}\mathrm{\Delta }{\mathit{\tau }}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}}^{\mathit{k}}\mathrm{)}\mathrm{-}\stackrel{}{{{\mathit{\alpha }}_{\mathit{i,l}}^{\mathit{k}}}_{\mathrm{˜}}}\frac{{\mathit{a}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}}^{\mathit{k}}}{\mathit{\chi ̂}\begin{array}{}\mathit{k}\\ \mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}\end{array}}\mathrm{[}\mathrm{4}\mathrm{+}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{]}{{\mathit{p}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}}^{\mathrm{|}\mathit{,k}}}^{\mathrm{]}}\\ {\mathit{B}}_{\mathit{i,l}}^{\setminus \mathit{,k}}& \mathrm{=}& \mathrm{-}\stackrel{}{{{\mathit{\beta }}_{\mathit{i,l}}^{\mathit{k}}}_{\mathrm{˜}}}\frac{{\mathit{a}}_{\mathit{i,l}}^{\mathit{k}}}{\mathit{\chi ̂}\begin{array}{}\mathit{k}\\ \mathit{i,l}\end{array}}\mathrm{[}\mathrm{4}\mathrm{+}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\xi }\mathrm{]}\hspace{0.17em}{{\mathit{p}}_{\mathit{i,l}}^{\mathrm{|}\mathit{,k}}}^{\mathrm{]}}\\ {\mathit{C}}_{\mathit{i,l}}^{\setminus \mathit{,k}}& \mathrm{=}& \mathrm{0}\\ {\mathit{A}}_{\mathit{i,l}}^{\mathrm{+}\mathit{,k}}& \mathrm{=}& \\ {\mathit{B}}_{\mathit{i,l}}^{\mathrm{+}\mathit{,k}}& \mathrm{=}& \\ {\mathit{C}}_{\mathit{i,l}}^{\mathrm{+}\mathit{,k}}& \mathrm{=}& \end{array}$Following Knop et al. (2009a), the naming scheme of the quantities in Eqs. (36)–(44) indicates with which specific intensity element they are associated. For an index pair i and l a • ⁻  superscript refers to an intensity at wavelength l − 1, a •^∖ superscript to the same wavelength, and • ⁺  to the next wavelength point l + 1. The A,B,C terms refer to the spatial points i − 1,i,i + 1, respectively. Equations (25)–(44) are nearly identical to those of Knop et al. (2009a) except for the explicit Doppler factor f(s), which arises because the photon direction is measured by a distant stationary observer here rather than by a comoving observer, as was the case in Knop et al. (2009a). We also clarified a problem with ${\mathit{C}}_{\mathit{i,l}}^{\setminus \mathit{,k}}$ that was confusing in Knop et al. (2009a). The formal solution matrix is therefore identical in form to that shown in Fig. 1 of Knop et al. (2009a).

An element of the source function vector ΔI is given by (Knop et al. 2009a) $\mathrm{\Delta }{\mathit{I}}_{\mathit{i,l}}^{\mathit{k}}\mathrm{=}{\mathit{\alpha }}_{\mathit{i,l}}^{\mathit{k}}\mathit{Ŝ}\begin{array}{}\mathit{k}\\ \mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}\end{array}\mathrm{+}{\mathit{\beta }}_{\mathit{i,l}}^{\mathit{k}}\mathit{Ŝ}\begin{array}{}\mathit{k}\\ \mathit{i,l}\end{array}\mathrm{+}{\mathit{\gamma }}_{\mathit{i,l}}^{\mathit{k}}\mathit{Ŝ}\begin{array}{}\mathit{k}\\ \mathit{i}\mathrm{+}\mathrm{1}\mathit{,l}\end{array}\mathit{.}$(45)Note that all Doppler factors are explicitlyhandled by including them in the opacity as χ_l ∗ f(s).

From Eq. (35) the solution for the specific intensity for a given spatial point and wavelength reads $\begin{array}{ccc}{\mathit{I}}_{\mathit{i,l}}^{\mathit{k}}& \mathrm{=}& {\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{B}}_{\mathit{i,l}}^{\mathrm{diag}\mathit{,k}}}^{\mathrm{)}}}^{-1}\mathrm{·}(\mathrm{\Delta }{\mathit{I}}_{\mathit{i,l}}^{\mathit{k}}\mathrm{+}{\mathit{B}}_{\mathit{i,l}}^{\mathrm{sub}\mathit{,k}}{\mathit{I}}_{\mathit{i,l}\mathrm{-}\mathrm{1}}^{\mathit{k}}\mathrm{+}{\mathit{B}}_{\mathit{i,l}}^{\mathrm{super}\mathit{,k}}{\mathit{I}}_{\mathit{i,l}\mathrm{+}\mathrm{1}}^{\mathit{k}}\\ & & \mathrm{+}{\mathit{A}}_{\mathit{i,l}}^{\mathrm{sub}\mathit{,k}}{\mathit{I}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}\mathrm{-}\mathrm{1}}^{\mathit{k}}\mathrm{+}{\mathit{A}}_{\mathit{i,l}}^{\mathrm{diag}\mathit{,k}}{\mathit{I}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}}^{\mathit{k}}\mathrm{+}{\mathit{A}}_{\mathit{i,l}}^{\mathrm{super}\mathit{,k}}{\mathit{I}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,l}\mathrm{+}\mathrm{1}}^{\mathit{k}}\\ & & \end{array}$(46)

2.2.1. The operator splitting method

Now that we have the formal solution, the full scattering problem can be solved by using operator splitting. The mean intensity J_λ is obtained from the source function S_λ by a formal solution of the radiative transfer equation, which is symbolically written using the Λ-operator Λ_λ as ${\mathit{J}}_{\mathit{\lambda }}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{\lambda }}{\mathit{S}}_{\mathit{\lambda }}\mathit{.}$(47)For the transition of a two-level atom, we have $\mathit{J̅}\mathrm{=}\mathrm{\Lambda ̅}\mathit{S,}$(48)where , with the normalized line profile φ(λ).

The Λ-iteration method, i.e. solving Eq. (48) by a fixed-point iteration scheme of the form $\mathit{J̅}\mathrm{new}\mathrm{=}\mathrm{\Lambda ̅}{\mathit{S}}_{\mathrm{old}}\mathit{,} {\mathit{S}}_{\mathrm{new}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{J̅}\mathrm{new}\mathrm{+}\mathit{ϵB,}$(49)fails in the case of high optical depths and small ϵ. This is because the largest eigenvalue of the amplification matrix (for Doppler-profiles) is approximately λ_max ≈ (1 − ϵ)(1 − T^-1), where T is the optical thickness of the medium (Mihalas et al. 1975). For small ϵ and high T, this is very close to unity and, therefore, the convergence rate of the Λ-iteration is very poor. A physical description of this effect can be found in Mihalas (1980).

Fig. 1 Comparison of monotonic flow case with spherical symmetry solved using the full arbitrary velocity field method (red) to the well-tested 1D solution (black). The comoving mean intensity is plotted at each point in the computational volume. The agreement is at the 1% level.

The idea of the ALI or operator splitting method is to reduce the eigenvalues of the amplification matrix in the iteration scheme (Cannon 1973) by introducing an approximate -operator (ALO) Λ^∗ and to split according to $\mathrm{\Lambda ̅}\mathrm{=}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{+}\mathrm{(}\mathrm{\Lambda ̅}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{)}$(50)and rewrite Eq. (48) as $\mathit{J̅}\mathrm{new}\mathrm{=}{\mathrm{\Lambda }}^{\mathrm{\ast }}{\mathit{S}}_{\mathrm{new}}\mathrm{+}\mathrm{(}\mathrm{\Lambda ̅}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{)}{\mathit{S}}_{\mathrm{old}}\mathit{.}$(51)This relation can be written as (Hamann 1987) $\left[\mathrm{1}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\right]\mathit{J̅}\mathrm{new}\mathrm{=}\mathit{J̅}\mathrm{fs}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{J̅}\mathrm{old}\mathit{,}$(52)where . Equation (52) is solved to obtain the new values of , which are then used to compute the new source function for the next iteration cycle.

Mathematically, the ALI method belongs to the same family of iterative methods as the Jacobi or the Gauss-Seidel methods (Golub & Van Loan 1989). These methods have the general form $\mathit{M}{\mathit{x}}^{\mathit{k}\mathrm{+}\mathrm{1}}\mathrm{=}\mathit{N}{\mathit{x}}^{\mathit{k}}\mathrm{+}\mathit{b}$(53)for the iterative solution of a linear system Ax = b where the system matrix A is split according to A = M − N. For the ALI method we have M = 1 − Λ^∗(1 − ϵ) and, accordingly, for the system matrix . The convergence of the iterations depends on the spectral radius, ρ(G), of the iteration matrix G = M^-1N. For convergence the condition ρ(G) < 1 must be fulfilled, this puts a restriction on the choice of Λ^∗. In general, the iterations will converge faster for a smaller spectral radius. To achieve a significant improvement compared to the Λ-iteration, the operator Λ^∗ is constructed so that the eigenvalues of the iteration matrix G are much less than unity, resulting in swift convergence. Using parts of the exact matrix (e.g., its diagonal or a tri-diagonal form) will optimally reduce the eigenvalues of the G. The calculation and the structure of Λ^∗ should be simple to make the construction of the linear system in Eq. (52) fast. For example, the choice  is best in view of the convergence rate (it is equivalent to a direct solution by matrix inversion) but the explicit construction of is more time-consuming than the construction of a simpler Λ^∗.

In the following discussion we use the notation of Hauschildt (1992) and Hauschildt & Baron (2006). The basic framework and the methods used for the formal solution and the solution of the scattering problem via operator splitting are discussed in detail in Hauschildt & Baron (2006) and will therefore not be repeated here. We have extended the framework to solve line transfer problems with a background continuum. The basic approach is similar to that of Hauschildt (1993). In the simple case of a two-level atom with background continuum that we consider here as a test case, we use a wavelength grid that covers the profile of the line including the surrounding continuum. We then use the wavelength-dependent mean intensities J_λ and approximate Λ operators Λ^∗ to compute the profile-integrated line mean intensities and $\overline{{\mathrm{\Lambda }}^{\mathrm{\ast }}}$ via

$\mathit{J̅}\mathrm{=}\mathrm{\int }\mathit{\phi }\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}{\mathit{J}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{d}\mathit{\lambda }$and

$\overline{{\mathrm{\Lambda }}^{\mathrm{\ast }}}\mathrm{=}\mathrm{\int }\mathit{\phi }\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}{\mathrm{\Lambda }}^{\mathrm{\ast }}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\mathit{.}$ and $\overline{{\mathrm{\Lambda }}^{\mathrm{\ast }}}$ are then used to compute an updated value for and the line source function

$\mathit{S}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{J̅}\mathrm{+}\mathit{ϵB,}$where ϵ is the line thermalization parameter (0 for a purely absorptive line, 1 for a purely scattering line). B is the Planck function, B_λ, profile averaged over the line

$\mathit{B}\mathrm{=}\mathrm{\int }\mathit{\phi }\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}{\mathit{B}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{d}\mathit{\lambda }$via the standard iteration method $\left[\mathrm{1}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\right]\mathit{J̅}\mathrm{new}\mathrm{=}\mathit{J̅}\mathrm{fs}\mathrm{-}{\mathrm{\Lambda }}^{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{ϵ}\mathrm{)}\mathit{J̅}\mathrm{old}\mathit{,}$(54)where . This equation is solved directly to obtain the new values of , which are then used to compute the new source function for the next iteration cycle.

We construct the line directly from the wavelength-dependent Λ^∗ generated by the solution of the continuum transfer problems.

Given the form of Eq. (46) for the formal solution, the construction of the Λ^∗-operator can proceed exactly as described in Baron & Hauschildt (2004). However, to conserve memory, we have implemented the Λ^∗-operator, retaining the spatial off-diagonal terms, but neglecting the off-diagonal terms in wavelength. We still find good convergence at considerable savings in memory (see below).

3. Test calculations

In this section we present the results of test calculations we performed to test the new algorithm in terms of accuracy by regression testing. We compare these to results of the homologous case (Baron et al. 2009b), and to the spherical nonmonotonic cases (Baron & Hauschildt 2004).

The test calculations were performed on Opteron CPUs running Linux (Franklin and Hopper at NERSC), on Intel CPUs (Carver at NERSC, ICE2 at HLRN, and our own local Xserve clusters), and on IBM CPUs (PWR-4 and PWR-5). The code was compiled with Gfortran/gcc/g++, ifort/icc/icpc (versions 11 and 12), and xlf/xlc/xlC and with NAG f95/gcc/g++. Using the varied compiler suites and CPUs allowed us to find numerous errors.

3.1. Regression with monotonic case

Figure 1 shows the profile of the mean intensity J for homologous flow and spherical symmetry. The test problem is similar to that of Baron et al. (2009b). These are the basic model parameters:

• 1.

  An inner radius R_in = 10¹¹ cm and an outer radius R_out = 1.01 × 10¹³ cm.

• 2.

  A minimum optical depth in the continuum ${\mathit{\tau }}_{\mathrm{std}}^{\min }\mathrm{=}{\mathrm{10}}^{-4}$ and a maximum optical depth in the continuum ${\mathit{\tau }}_{\mathrm{std}}^{\max }\mathrm{=}{\mathrm{10}}^{\mathrm{4}}$.

• 3.

  A gray temperature structure with K. That is the temperature solution of the spherical gray atmosphere problem with effective temperature (Mihalas 1978).

• 4.

  An outer boundary condition ${\mathit{I}}_{\mathrm{bc}}^{\mathrm{-}}\mathrm{\equiv }\mathrm{0}$ and an inner boundary condition I_λ = B_λ for all wavelengths.

• 5.

  For the initial wavelength the boundary condition is taken from that given by the 3D homologous calculation for homologous tests and set equal to the Planck function B_λinit for non-homologous tests.

• 6.

  A continuum extinction χ_c = C/r², with the constant C fixed by the radius and optical depth grids.

• 7.

  A parametrized coherent and isotropic continuum scattering given by

  ${\mathit{\chi }}_{\mathrm{c}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{c}}{\mathit{\kappa }}_{\mathrm{c}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{ϵ}}_{\mathrm{c}}\mathrm{)}{\mathit{\sigma }}_{\mathrm{c}}$(55)

  with 0 ≤ ϵ_c ≤ 1. κ_c and σ_c are the continuum absorption and scattering coefficients. In this work we have neglected scattering in the continuum.

The line of the simple two-level model atom is parameterized by the ratio of the profile-averaged line opacity χ_l to the continuum opacity χ_c and the line thermalization parameter ϵ_l. For the test cases presented below, we used ϵ_c = 1 and the line strength is given by Γ ≡ χ_l/χ_c = 10² to simulate a strong line, with varying ϵ_l (see below).

The test model is just an optically thick sphere put into the 3D grid. The velocity at the outer radius was set to be relativistic, v_max = 8 × 10⁴ km s^-1. The calculation was performed on a spherical grid with 33³ spatial points, 33² solid angle directions, and 22 wavelength points. This and all calculations presented below were parallelized over characteristics and were run on parallel clusters. Ng (Ng 1974; Auer 1987) acceleration was used to significantly speed up the operator splitting iterations for cases with scattering. Figure 2 shows the line profile at the surface, here we compare the 3D monotonic calculation (Baron et al. 2007) to the same calculation using the 3D arbitrary velocity algorithm. Because of the spatial resolution and the way that characteristics end at different points in a voxel (as opposed to always ending on a radial grid point in the 1D case), it is better to compare 3D cases to 3D cases.

Fig. 2 Comparison of monotonic flow case with spherical symmetry solved using the full arbitrary velocity field method (green) to the 3D monotonic flow solution (red). The comoving mean intensity is plotted at the surface of the sphere. The test cases are identical, vmax = 8 × 104 km s-1, 129 × 17 × 17 spatial voxels, 114 wavelength points, and 256 × 256 directions. The thermalization parameter in the line is ϵ = 10-3 and the line strength is Γ = 100. The lower panel shows the relative difference, δJλ/Jλ, as a function of wavelength.

Figure 3 shows the convergence rate of a monotonic flow case, with ϵ = 10^-4, with and without Ng accleration. Not only does Ng acceleration clearly improve the convergence rate, it gives almost the exact same errors as the same test setup using both the 3D and 1D homologous algorithms.

Fig. 3 Comparison of the convergence rate of a monotonic flow case with spherical symmetry with the thermalization parameter ϵ = 10-4 with and without Ng acceleration. The Ng acceleration is identical to that produced by the pure homologous 3D module.

Fig. 4 Velocity profile used for the nonmonotonic velocity test with a damped sine wave.

3.2. Sine velocity flow

Here we again consider a spherically symmetric case with the same physical parameters as in Sect. 3.1, but the velocity, while still radial is now given modulated by a damped sine wave. We have also introduced a line into the opacity with thermalization parameter ϵ_l = 1. The thermalization parameter in the continuum is ϵ_c = 1. This case was calculated in 1D in Baron & Hauschildt (2004). Figure 4 shows the velocity as a function of radial optical depth in the continuum. The spatial grid was 129 × 65 × 17 and the solid angle resolution was set to 512². The maximum velocity is v ≈ 10 000 km s^-1 and the total number of wavelength points is 226. Figure 5 shows the comoving mean intensity J_λ from each surface voxel as a function of wavelength λ. The 1D result is plotted (the green curve) and the agreement is good to the sub-percent level. The variation in the velocity leads to an asymmetric line profile even in the comoving frame. But since this test has such a high spatial resolution, it requires a significant amount of memory per process. We explored the effects of lower spatial resolution, 129 × 17 × 17, while keeping the high solid angle resolution 512 × 512. Figure 6 shows the mean intensity as a function of wavelength for the outmost part of the sphere. For computational expediency we set the line thermalization parameter ϵ_l = 1. The spread in the results of 1.5% is indicative of the accuracy of these calculations, whereas the deviation of the envelope of the solution is indicative of the low spatial resolution of this calculation. While modern nodes may have 24 CPUs, but only about 1 GB of RAM per CPU, it is unfortunate that the evelope of the low spatial resolution calculation is offset from the 1D result. This shows that not only solid angle resolution is important, but spatial resolution is quite important as well. Figure 7 shows the same calculation with the spatial grid enhanced to 129 × 65 × 65 and the solid angle grid reduced to 32 × 24 (Hauschildt & Baron 2010). As can be seen in Fig. 7, the reduced solid angle resolution increases the spread by roughly a factor of two, but now the 3D solution envelopes the 1D nonmonotonic result. With this higher spatial resolution, the spread can be further reduced by increasing the solid angle resolution, which shows near-perfect weak scaling, and thus does not increase the wallclock time required for these calculations (although it does require more CPUs).

Fig. 5 Line profile for the damped sine-wave velocity case with high spatial nr = 129, nθ = 65, and nφ = 17, and angular resolution nΩ = 512 × 512. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 10-4 and the line strength is Γ = 100.

Fig. 6 Line profile for the damped sine-wave velocity case with lower spatial nr = 129, nθ = 17, and nφ = 17, and higher angular resolution nΩ = 512 × 512 than in Fig. 5. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 1 for computational expediency, and the line strength is Γ = 100. With the very low spatial resolution, the 3D result no longer envelopes the 1D result (which by assuming spatial spherical symmetry has essentially infinite resolution in nθ, and nφ).

Fig. 7 Line profile for the damped sine-wave velocity case with higher spatial resolution, but lower angular resolution than in Fig. 5. The spatial resolution is nr = 129, nθ = 65, and nφ = 65, and the angular resolution is reduced to nΩ = 32 × 24. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 1 for computational expediency, and the line strength is Γ = 100.

3.3. Example of radial nonhomologous flows

The previous tests were all totally spherically symmetric, with radial variations in the velocity field. We now assume the flow to be radial and azimuthally symmetric, i.e. ${\beta }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathit{p}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{q}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hat{{e}}\mathit{r}\mathit{.}$(56)For this case, we have $\begin{array}{ccc}\hat{{n}}\mathrm{·}{\beta }& \mathrm{=}& \\ \mathrm{(}\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{)}\mathit{\beta }& \mathrm{=}& \end{array}$and $\hat{{n}}\mathrm{·}\mathrm{\nabla }\mathrm{(}\hat{{n}}\mathrm{·}{\beta }\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}\mathrm{[}{\mathit{p}}^{\mathrm{\prime }}\mathit{q}{\mathit{n}}_{\mathit{r}}{\mathit{n}}_{\mathit{\theta }}\mathrm{+}\mathit{p}\left({\mathit{q}}^{\mathrm{\prime }}\mathit{r}\right){\mathit{n}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathit{pq}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{n}}_{\mathit{r}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{]}}\mathit{,}$(59)where we made use of $\mathit{ṙ}\mathrm{=}{\mathit{n}}_{\mathit{r}}\mathit{,}\mathit{\theta ̇}\mathrm{=}\frac{{\mathit{n}}_{\mathit{\theta }}}{\mathit{r}}\mathit{,}\mathit{\phi ̇}\mathrm{=}\frac{{\mathit{n}}_{\mathit{\phi }}}{\mathit{r}\sin \mathit{\theta }}\mathit{,}$(60)and $\mathit{ṅ}\mathit{r}\mathrm{=}{\mathit{n}}_{\mathit{\theta }}\mathit{\theta ̇}\mathrm{+}\sin \mathit{\theta }{\mathit{n}}_{\mathit{\phi }}\mathit{\phi ̇}\mathit{.}$(61)Here an over dot implies $\frac{\mathrm{d}}{\mathrm{d}\mathit{s}}\mathit{,}$ e.g., $\mathit{ṙ}\mathrm{\equiv }\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{s}}\mathit{,}$ etc., and a prime denotes differentiation with respect to μ = cosθ.

Inserting Eqs. (57)–(59) into Eq. (19), we find an analytic expression for a(s): $\begin{array}{ccc}\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{p}}^{\mathrm{\prime }}\mathit{q}{\mathit{n}}_{\mathit{r}}{\mathit{n}}_{\mathit{\theta }}\mathrm{+}\mathit{p}\mathrm{\left(}{\mathit{q}}^{\mathrm{\prime }}\mathit{r}\mathrm{\right)}{\mathit{n}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathit{pq}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{n}}_{\mathit{r}}^{\mathrm{2}}\mathrm{)}}{\mathit{r}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{pq}{\mathit{n}}_{\mathit{r}}\mathrm{)}}\\ & & \mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathit{pq}\left[\mathit{p}{\mathit{n}}_{\mathit{r}}{\mathit{q}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{q}{\mathit{n}}_{\mathit{\theta }}\right]\mathit{,}\end{array}$(62)where the spherical components of the unit vector $\hat{{n}}\mathrm{=}\mathrm{\left(}\sin {\mathit{\theta }}_{\mathit{n}}\cos {\mathit{\phi }}_{\mathit{n}}\mathit{,}\sin {\mathit{\theta }}_{\mathit{n}}\sin {\mathit{\phi }}_{\mathit{n}}\mathit{,}\cos {\mathit{\theta }}_{\mathit{n}}\mathrm{\right)}$ are $\begin{array}{ccc}{\mathit{n}}_{\mathit{r}}& \mathrm{=}& \sin \mathit{\theta }\sin {\mathit{\theta }}_{\mathit{n}}\cos \mathrm{(}\mathit{\phi }\mathrm{-}{\mathit{\phi }}_{\mathit{n}}\mathrm{)}\mathrm{+}\cos \mathit{\theta }\cos {\mathit{\theta }}_{\mathit{n}}\mathit{,}\\ {\mathit{n}}_{\mathit{\theta }}& \mathrm{=}& \cos \mathit{\theta }\sin {\mathit{\theta }}_{\mathit{n}}\cos \mathrm{(}\mathit{\phi }\mathrm{-}{\mathit{\phi }}_{\mathit{n}}\mathrm{)}\mathrm{-}\sin \mathit{\theta }\cos {\mathit{\theta }}_{\mathit{n}}\mathit{,}\\ {\mathit{n}}_{\mathit{\phi }}& \mathrm{=}& \cos {\mathit{\theta }}_{\mathit{n}}\mathit{.}\end{array}$(63)For p(θ), we use an expansion in terms of Legendre polynomials (which form a complete and orthonormal basis for axially symmetric spherical functions) $\mathit{p}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathit{k}}{\mathit{c}}_{\mathit{n}}{\mathrm{𝒫}}_{\mathit{n}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(64)with and ${\mathrm{𝒫}}_{\mathrm{2}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{3}{\mathit{\mu }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}$ etc., where μ = cosθ. We consequently obtain ${\mathit{p}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{k}}{\mathit{c}}_{\mathit{n}}{\mathrm{𝒫}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$(65)with ${\mathrm{𝒫}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\sin \mathit{\theta ,}$ and ${\mathrm{𝒫}}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\sin \mathrm{2}\mathit{\theta ,}$ etc. (note that the case ${\mathit{c}}_{\mathit{n}}\mathrm{=}{\mathit{c}}_{\mathrm{0}}{\mathit{\delta }}_{\mathit{n}}^{\mathrm{0}}$ degenerates to homologous flow). We can use this finite expansion in terms of the Legendre polynomial to construct azimuthally symmetric jets.

As a simple example of nonspherically symmetric flow, we run a test case where ${\beta }\mathrm{=}{\mathit{c}}_{\mathrm{0}}\mathit{rp}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hat{{e}}\mathit{r}\mathit{,}$(66)i.e. we take q(r) = c₀r (c₀ is a constant), which simplifies a(s) to be $\mathit{a}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{0}}\left[\frac{\mathit{p}\mathrm{+}{\mathit{p}}^{\mathrm{\prime }}{\mathit{n}}_{\mathit{r}}{\mathit{n}}_{\mathit{\theta }}}{\mathrm{1}\mathrm{-}\mathit{r}{\mathit{n}}_{\mathit{r}}\mathit{p}}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathit{\beta }\mathrm{(}\mathit{p}{\mathit{n}}_{\mathit{r}}\mathrm{+}{\mathit{p}}^{\mathrm{\prime }}{\mathit{n}}_{\mathit{\theta }}\mathrm{)}\right]\mathit{.}$(67)Furthermore, we assume $\mathit{p}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{𝒫}}_{\mathrm{2}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{3}}{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\right)\mathit{,}$(68)i.e., we perturb the homologous flow by adding a second-order Legendre polynomial with perturbation coefficient 0.5. We show the resulting mean intensity plot J(θ,φ) at the boundary R = R_out in Fig. 8. We obtained perfect azimuthal symmetry as expected from the form of β in Eq. (66), although it was not explicitly enforced. Furthermore, we also recovered the symmetry with respect to the north/south pole (i.e., symmetry under a reflection with respect to the equatorial plane, θ → π − θ), which is characteristic for Legendre polynomials of even order, see Eq. (68). Figure 9 shows a longitudinal slice of the comoving intensity, which shows that the comoving mean intensity varies by roughly a factor of ten from the pole to the equator. This is just due to the effect that the continuum level varies with the maximum velocity (Baron et al. 2009b).

Fig. 8 Comoving mean intensity J(θ,φ) plotted at the surface R = Rout for case where ${\beta }\mathrm{=}{\mathit{c}}_{\mathrm{0}}\mathit{rp}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hat{{e}}{}_{\mathit{r}}\mathit{,}$ and $\mathit{p}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{𝒫}}_{\mathrm{2}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}$$\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{)}$. The velocity flow is radial, but not spherically symmetric to approximate a jet-like flow in the ± z direction.

Fig. 9 Comoving mean intensity J(θ,φ0) plotted at the surface R = Rout for a fixed (arbitrarily chosen) value of φ = φ0. The variation of nearly an order of magnitude from the pole to the equator is due to the fact the the continuum level is a function of the velocity.

3.4. Checkpointing

We implemented a checkpointing scheme that allows for restarts and new starts with higher resolution in solid angle. We simply write out (using stream I/O) the Λ^∗ operator once it has been calculated, and the value of J at each voxel after each ALO iteration. Since these calculations require significant numbers of processors, which may go down, or calculations may run out of time, this allows us to perform perfect restarts. We checked that restarting from the checkpoint files works perfectly and that the calculations are converged in a single iteration when restarting from a converged checkpoint file. We write Λ^∗ and J to separate files, since that allows us to read just J, and then recalculate Λ^∗ if we restart with a higher resolution. We checked this, and it allows us to use a low-resolution calculation (performed on fewer processors, for example) and then converge the higher resolution calculation with roughly half as many iterations for each resolution increase. Table 1 shows the rate of convergence for a test with homologous velocity fields (chosen for computational expedience), v_max = 10⁵ km s^-1 and ϵ = 10^-5. While the restart result still requires several iterations, this scheme allows for some speedup. The small ϵ and high maximum velocity make this test particularly demanding. The calculations were all restarted such that Ng acceleration cannot begin until the fourth iteration, and thus the restart calculations converge more slowly than they could. The scheme could indeed possibly be made more efficient by keeping the previous three values of , so that the restart iteration could immediately use Ng accleration.

4. Conclusion

We presented algorithm strategies and details for the solution of the radiative transfer problem in 3D atmospheres with arbitrary wavelength couplings.

Future possible applications are the velocity profiles of cool stellar winds, treating partial redistribution, calculating radiative transfer in shock fronts like in accretion shocks, calculating quasar jets, and general relativistic situations such as a rotating black hole.

Acknowledgments

This work was supported in part by SFB 676, GRK 1354 from the DFG, NSF grant AST-0707704, and US DOE Grant DE-FG02-07ER41517. Support for Program number HST-GO-12298.05-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. B.C. thanks the University of Okalahoma Foundation for a fellowship. This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231; and the Höchstleistungs Rechenzentrum Nord (HLRN). We thank both these institutions for a generous allocation of computer time.

References

All Tables

All Figures

	Fig. 1 Comparison of monotonic flow case with spherical symmetry solved using the full arbitrary velocity field method (red) to the well-tested 1D solution (black). The comoving mean intensity is plotted at each point in the computational volume. The agreement is at the 1% level.
In the text

Fig. 2 Comparison of monotonic flow case with spherical symmetry solved using the full arbitrary velocity field method (green) to the 3D monotonic flow solution (red). The comoving mean intensity is plotted at the surface of the sphere. The test cases are identical, vmax = 8 × 104 km s-1, 129 × 17 × 17 spatial voxels, 114 wavelength points, and 256 × 256 directions. The thermalization parameter in the line is ϵ = 10-3 and the line strength is Γ = 100. The lower panel shows the relative difference, δJλ/Jλ, as a function of wavelength.

In the text

	Fig. 3 Comparison of the convergence rate of a monotonic flow case with spherical symmetry with the thermalization parameter ϵ = 10-4 with and without Ng acceleration. The Ng acceleration is identical to that produced by the pure homologous 3D module.
In the text

	Fig. 4 Velocity profile used for the nonmonotonic velocity test with a damped sine wave.
In the text

Fig. 5 Line profile for the damped sine-wave velocity case with high spatial nr = 129, nθ = 65, and nφ = 17, and angular resolution nΩ = 512 × 512. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 10-4 and the line strength is Γ = 100.

In the text

Fig. 6 Line profile for the damped sine-wave velocity case with lower spatial nr = 129, nθ = 17, and nφ = 17, and higher angular resolution nΩ = 512 × 512 than in Fig. 5. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 1 for computational expediency, and the line strength is Γ = 100. With the very low spatial resolution, the 3D result no longer envelopes the 1D result (which by assuming spatial spherical symmetry has essentially infinite resolution in nθ, and nφ).

In the text

Fig. 7 Line profile for the damped sine-wave velocity case with higher spatial resolution, but lower angular resolution than in Fig. 5. The spatial resolution is nr = 129, nθ = 65, and nφ = 65, and the angular resolution is reduced to nΩ = 32 × 24. The comoving mean intensity Jλ for each voxel on the surface (there are 65 × 17 = 1105 of them) is plotted as a function of λ (red lines). The green curve is the 1D result using the methods of Baron & Hauschildt (2004) and Knop et al. (2009b,a). Here the scattering fraction is ϵ = 1 for computational expediency, and the line strength is Γ = 100.

In the text

Fig. 8 Comoving mean intensity J(θ,φ) plotted at the surface R = Rout for case where ${\beta }\mathrm{=}{\mathit{c}}_{\mathrm{0}}\mathit{rp}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hat{{e}}{}_{\mathit{r}}\mathit{,}$ and $\mathit{p}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{𝒫}}_{\mathrm{2}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}$$\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{)}$. The velocity flow is radial, but not spherically symmetric to approximate a jet-like flow in the ± z direction.

In the text

	Fig. 9 Comoving mean intensity J(θ,φ0) plotted at the surface R = Rout for a fixed (arbitrarily chosen) value of φ = φ0. The variation of nearly an order of magnitude from the pole to the equator is due to the fact the the continuum level is a function of the velocity.
In the text