A Fast and Accurate Method of Radiation Hydrodynamics Calculation in Spherical Symmetry

The monochromatic radiation transfer equation in three dimensions is

$$\begin{eqnarray}&&\displaystyle \frac{1}{c}\displaystyle \frac{{{dI}}_{\nu }}{{dt}}+{\boldsymbol{n}}\cdot {\rm{\nabla }}{I}_{\nu }={\kappa }_{\nu }\rho ({S}_{\nu }-{I}_{\nu }),\end{eqnarray} \tag{ 1 }$$

where I_ν, κ_ν, and S_ν are the intensity, opacity, and source function at frequency ν. In the following, we will leave out the index ν for convenience.

Equation (1) combines the effects of emission and absorption of radiation by the material. We shall neglect the explicitly time-dependent term dI_ν/dt throughout this paper. This is equivalent to assuming that the radiation field instantaneously adjusts to any changes in the hydrodynamic variables. This assumption is valid as long as the time for light to cross the system is much shorter than the time for the hydrodynamic variables to change significantly, which is the case in most astrophysical systems.

With the above assumption, the solution is

$$\begin{eqnarray}&&I({\boldsymbol{r}},{\boldsymbol{n}})\\ &&\quad =\,{\displaystyle \int }_{0}^{{D}_{{\rm{B}}}}\kappa ({\boldsymbol{r}}+{\boldsymbol{n}}D)\rho ({\boldsymbol{r}}+{\boldsymbol{n}}D)S({\boldsymbol{r}}+{\boldsymbol{n}}D){e}^{-\tau ({\boldsymbol{r}},{\boldsymbol{r}}+{\boldsymbol{n}}D)}{dD}\\ &&\quad \quad +\,{S}_{{\rm{O}}}{e}^{-{\tau }_{{\rm{B}}}},\end{eqnarray} \tag{ 2 }$$

where $\tau ({\boldsymbol{r}},{\boldsymbol{r}}+{\boldsymbol{n}}D)$ is the optical depth between positions ${\boldsymbol{r}}$ and ${\boldsymbol{r}}+{\boldsymbol{n}}D$, while D_B and τ_B refer to the distance and optical depth between ${\boldsymbol{r}}$ and the outer boundary of the system in the direction ${\boldsymbol{n}}$. The first term represents radiation emitted by other points within the system, while the second term represents radiation coming from outside. The latter's intensity is determined by S_O, the source function of the external radiation, which must be supplied as a boundary condition.

Equation (2) determines the direction-dependent intensity. The integral over all directions gives us the radiation energy density:

$$\begin{eqnarray}&&\epsilon ({\boldsymbol{r}})=\displaystyle \frac{1}{c}{\int }_{{\rm{\Omega }}}\left({\int }_{0}^{{D}_{{\rm{B}}}}\kappa ^{\prime} \rho ^{\prime} S^{\prime} {e}^{-\tau }{dD}+{S}_{{\rm{O}}}{e}^{-{\tau }_{{\rm{B}}}}\right)d{\rm{\Omega }}.\end{eqnarray} \tag{ 3 }$$

A prime on a position-dependent quantity indicates that it is to be evaluated at ${\boldsymbol{r}}+{\boldsymbol{n}}D$ rather than ${\boldsymbol{r}}$.

Let us now consider the total radiation term (heating plus cooling). Assuming that each fluid element emits isotropically,

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=\kappa (c\epsilon -4\pi S),\end{eqnarray} \tag{ 4 }$$

and inserting the previously derived expression for , we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=\kappa \left[{\int }_{{\rm{\Omega }}}\left({\int }_{0}^{{D}_{{\rm{B}}}}\kappa ^{\prime} \rho ^{\prime} S^{\prime} {e}^{-\tau }{dD}+{S}_{O}{e}^{-{\tau }_{{\rm{B}}}}\right)d{\rm{\Omega }}-4\pi S\right].\end{eqnarray} \tag{ 5 }$$

Assume a case where the temperature, and therefore the source function, is constant everywhere, $(S=S^{\prime} ={S}_{{\rm{O}}})$. In this case, $\tfrac{{dU}}{{dt}}$ must be zero, and it follows that

$$\begin{eqnarray}&&{\int }_{{\rm{\Omega }}}\left({\int }_{0}^{{D}_{{\rm{B}}}}\kappa ^{\prime} \rho ^{\prime} {e}^{-\tau }{dD}+{e}^{-{\tau }_{{\rm{B}}}}\right)d{\rm{\Omega }}=4\pi .\end{eqnarray} \tag{ 6 }$$

We can now replace the 4π in Equation (4) with this expression:

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=\kappa {\int }_{{\rm{\Omega }}}\left({\int }_{0}^{{D}_{{\rm{B}}}}\kappa ^{\prime} \rho ^{\prime} (S^{\prime} \mbox{--}S){e}^{-\tau }{dD}+({S}_{{\rm{O}}}-S){e}^{-{\tau }_{{\rm{B}}}}\right)d{\rm{\Omega }}.\end{eqnarray} \tag{ 7 }$$

Consider the physical meaning of the integration in the first term. For all directions, we integrate outward from a point ${\boldsymbol{r}}$ to the system's boundary. In other words, we are integrating over the volume of the physical system but using a coordinate system whose origin lies at ${\boldsymbol{r}}$, so the volume element is $d{\boldsymbol{r}}^{\prime} ={D}^{2}d{\rm{\Omega }}{dD}$, and we may write

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=\kappa \left({\int }_{V}\kappa ^{\prime} \rho ^{\prime} (S^{\prime} \,\mbox{--}\,S)\displaystyle \frac{{e}^{-\tau }}{{D}^{2}}d{\boldsymbol{r}}^{\prime} +{\int }_{{\rm{\Omega }}}({S}_{{\rm{O}}}-S){e}^{-{\tau }_{{\rm{B}}}}d{\rm{\Omega }}\right).\end{eqnarray} \tag{ 8 }$$

As long as the integral is over the whole volume, we may use the usual coordinate system and volume element. From this point onward, we shall abbreviate the product κρ as χ (the extinction/emission coefficient). The physical interpretation of this equation is as follows.

The volume integral represents the exchange of radiative energy between different positions within the system. Thus, the term with S (cooling term) is the radiation emitted at ${\boldsymbol{r}}$ and absorbed at some other point ${\boldsymbol{r}}^{\prime}$, while the term with S' (heating term) is radiation emitted at ${\boldsymbol{r}}^{\prime}$ and absorbed at ${\boldsymbol{r}}$.

Likewise, the S_O term in the solid angle integral represents radiation coming from outside the system that is absorbed at ${\boldsymbol{r}}$, while the S term represents radiation emitted at ${\boldsymbol{r}}$ that reaches the boundary and escapes from the system.

The factors $\kappa \chi ^{\prime} {e}^{-\tau }/{D}^{2}$ and $\kappa {\int }_{{\rm{\Omega }}}{e}^{-{\tau }_{{\rm{B}}}}d{\rm{\Omega }}$ can be regarded as exchange coefficients, quantifying the efficiency of the radiative exchange between ${\boldsymbol{r}}$ and ${\boldsymbol{r}}^{\prime}$ and ${\boldsymbol{r}}$ and the outside radiation field, respectively. Furthermore, Equation (6) tells us that the volume integral over all of these exchange coefficients, plus the exchange coefficient for the boundary, must be equal to 4πκ. This is simply another way of saying that the radiative cooling is equal to −4πκ S.

In principle, Equation (8) can be integrated numerically to obtain the radiation term, but in realistic problems, the computational cost is prohibitively high.

Many hydrodynamical problems in astrophysics can be reasonably well approximated as spherically symmetric, so that all variables depend only on radius. While this reduces the computational cost, it is still quite slow because the equations include the distance and optical depth between different points. Even in spherical symmetry, these depend on all three coordinates. In this section, we shall manipulate the equation in such a way as to reduce this computational load for the case of a spherically symmetric system.

Consider the radiative exchange between a point ${\boldsymbol{r}}$ and an infinitesimally thin, spherical shell at radius r' (Figure 1). Calculating the radiation term at point ${\boldsymbol{r}}$ is sufficient, because spherical symmetry implies that it is the same for every point with radius r.

Zoom In Zoom Out Reset image size

If we choose the coordinate system such that ${\boldsymbol{r}}$ lies on the Cartesian y axis, there is no longer any dependency on the azimuthal angle, essentially reducing the system to 2D. Equation (8) tells us that the contribution of the (infinitesimal) shell r' to the radiation term at ${\boldsymbol{r}}$ is

$$\begin{eqnarray}{R}_{r\leftarrow r^{\prime} } & = & \kappa \chi ^{\prime} (S^{\prime} \mbox{--}S){\displaystyle \int }_{\theta ^{\prime} }{\displaystyle \int }_{\phi ^{\prime} }\displaystyle \frac{{e}^{-\tau }}{{D}^{2}}r{{\prime} }^{2}\sin (\theta ^{\prime} )d\theta ^{\prime} d\phi ^{\prime} {dr}^{\prime} \\ & = & 2\pi \kappa \chi ^{\prime} (S^{\prime} \mbox{--}S){\displaystyle \int }_{\theta ^{\prime} }\displaystyle \frac{{e}^{-\tau }}{{D}^{2}}r{{\prime} }^{2}\sin (\theta ^{\prime} )d\theta ^{\prime} {dr}^{\prime} \\ & = & 2\pi \kappa \chi ^{\prime} (S^{\prime} \mbox{--}S){A}_{r\leftarrow r^{\prime} }.\end{eqnarray} \tag{ 9 }$$

Numerical integration over r' is, of course, unavoidable. The main challenge is the calculation of the angle integral:

$$\begin{eqnarray}&&{A}_{r\leftarrow r^{\prime} }=r{{\prime} }^{2}{\int }_{0}^{\pi }\displaystyle \frac{{e}^{-\tau (\theta ^{\prime} )}}{{D}^{2}(\theta ^{\prime} )}\sin (\theta ^{\prime} )d\theta ^{\prime} {dr}^{\prime} .\end{eqnarray} \tag{ 10 }$$

Figure 1 illustrates the physical meaning of the variables in this equation. The next section deals with methods of calculating ${A}_{r\leftarrow r^{\prime} }$.

In order to perform numerical simulations, we subdivide our system into a finite number of spherically symmetric shells. Within each shell, density, opacity, and temperature are assumed to be constant. The contribution of shell j to the radiation term at radius r is determined by

$$\begin{eqnarray}&&{A}_{r\leftarrow j}={\int }_{{r}_{\mathrm{jL}}}^{{r}_{\mathrm{jR}}}r{{\prime} }^{2}{\int }_{0}^{\pi }\displaystyle \frac{{e}^{-\tau }}{{D}^{2}}\sin (\theta ^{\prime} )d\theta ^{\prime} {dr}^{\prime} .\end{eqnarray} \tag{ 11 }$$

Here r_jL and r_jR are the inner and outer boundaries of shell j, respectively. In order to obtain the exchange coefficient between shells i and j, the above expression would have to be integrated over shell i, i.e.,

$$\begin{eqnarray}&&{A}_{i\leftarrow j}=\displaystyle \frac{{\int }_{{m}_{\mathrm{iL}}}^{{m}_{\mathrm{iR}}}{A}_{r\leftarrow j}(m){dm}}{{M}_{i}}.\end{eqnarray} \tag{ 12 }$$

Here m is the mass coordinate, and M_i is the mass contained within shell i.

Unfortunately, the integration over shell i can only be done numerically, so this method is slow. In addition, if we numerically integrate over j as well (as in Section 2.3.3), there are large inaccuracies for neighboring shells near the singularity that occurs at r = r' for θ' = 0. For these reasons, we do not integrate over shell i, but instead we set ${A}_{i\leftarrow j}={A}_{{r}_{i}\leftarrow j}$, where r_i is the radius of shell i's center. This definition allows us to accurately calculate the radiation term at radius r_i, but note that energy conservation between shells requires that ${R}_{i\leftarrow j}=-{R}_{j\leftarrow i}$. Since the radiation term is power per mass, this implies that ${A}_{i\leftarrow j}\ {M}_{i}={A}_{j\leftarrow i}\ {M}_{j}$. With this definition, this is generally not fulfilled, since the radiation term at the shell center is not equal to the shell average. In other words, when updating shell i, we calculate the effect of shell j on radius r_i, but when updating shell j, we calculate the effect of shell i on radius r_j. This is not symmetric between i and j, which ultimately causes a violation of energy conservation.

Alternatively, it is possible to use the following definition:

$$\begin{eqnarray}&&{A}_{i\leftarrow j}=\displaystyle \frac{{A}_{{r}_{i}\leftarrow j}{M}_{i}+{A}_{{{\rm{r}}}_{j}\leftarrow i}{M}_{j}}{2{M}_{i}}.\end{eqnarray} \tag{ 13 }$$

The radiation term is now defined as an average between the two directions of transfer. This means the radiation term at r_i is no longer calculated as accurately, but this expression is symmetric to exchanging i and j and therefore conserves energy.

The method we develop in the following sections is able to calculate the radiation term at a specific radius exactly in the limit of infinite resolution. Therefore, for the test calculations in this paper and to investigate the effect of different resolutions, we use the nonsymmetrized definition, and from now on, we will refer to ${A}_{{r}_{i}\leftarrow j}$ simply as ${A}_{i\leftarrow j}$. In a real application, the symmetrized definition may be preferred, since it guarantees energy conservation.

The fastest way of calculating ${A}_{i\leftarrow j}$ is to substitute $\tau ={\overline{\chi }}_{{ij}}D$. This is an approximation: it assumes that ${\overline{\chi }}_{{ij}}$, the proportionality constant between distance and optical depth, is a constant independent of θ'. But unless χ is actually constant throughout the whole system, this is clearly not true. Consider again Figure 1. If χ increases toward the center, as is usually the case, then rays going through the inner region will have larger values of χ_ij than those that move only outward.

It is ultimately up to the user which value to choose for ${\overline{\chi }}_{{ij}};$ the simplest method, which we adopt here, is to take the value for θ' = 0. If χ increases toward the center, we expect this to cause an overestimation of ${A}_{i\leftarrow j}$, especially for the case i < j, and a somewhat smaller overestimation for i > j. Despite this inaccuracy, this method can be useful, because it allows us to solve the angle integral analytically and is therefore very fast. We start by expressing the distance D as a function of r_i, r_j, and θ':

$$\begin{eqnarray}{D}^{2} & = & {{r}_{i}}^{2}+{{r}_{j}}^{2}-2{r}_{i}{r}_{j}\cos (\theta ^{\prime} ),\\ & \Rightarrow & 2D\displaystyle \frac{{dD}}{d\theta ^{\prime} }=2{r}_{i}{r}_{j}\sin (\theta ^{\prime} ),\\ & \Rightarrow & \sin (\theta ^{\prime} )d\theta ^{\prime} =\displaystyle \frac{D}{{r}_{i}{r}_{j}}{dD}.\end{eqnarray} \tag{ 14 }$$

We insert this into the definition of ${A}_{i\leftarrow j}$ (Equation (11)) but ignore the integration over shell j for the moment and set r' = r_j. Using $\tau ={\overline{\chi }}_{{ij}}D$, we obtain

$$\begin{eqnarray}{A}_{i\leftarrow j} & = & \displaystyle \frac{{r}_{j}}{{r}_{i}}{\displaystyle \int }_{{D}_{\min }}^{{D}_{\max }}\displaystyle \frac{{e}^{-{\overline{\chi }}_{{ij}}D}}{D}{dD}{\rm{\Delta }}{r}_{j},\\ & = & \displaystyle \frac{{r}_{j}}{{r}_{i}}{\displaystyle \int }_{{\tau }_{\min }}^{{\tau }_{\max }}\displaystyle \frac{{e}^{-\tau }}{\tau }d\tau {\rm{\Delta }}{r}_{j}.\end{eqnarray} \tag{ 15 }$$

The indices min and max refer to the minimum and maximum distances/optical depths between the radii r_i and r_j (see Figure 1). The integral over τ can be reformulated using exponential integrals. The exponential integral of order n is defined as

$$\begin{eqnarray}&&{E}_{n}(x)={\int }_{x}^{\infty }\displaystyle \frac{{e}^{-}t}{{t}^{n}}{dt}.\end{eqnarray} \tag{ 16 }$$

Using this, Equation (15) becomes

$$\begin{eqnarray}&&{A}_{i\leftarrow j}=\displaystyle \frac{{r}_{j}}{{r}_{i}}({E}_{1}({\tau }_{\min })-{E}_{1}({\tau }_{\max })){\rm{\Delta }}{r}_{j}.\end{eqnarray} \tag{ 17 }$$

We have avoided the numerical integration over θ' and instead have to solve two exponential integrals. This is much faster, since analytic formulae are available for the latter. The τ_max term is negligible in optically thick systems, since in that case ${E}_{1}({\tau }_{\max })\ll {E}_{1}({\tau }_{\min })$.

We now consider the integration over shell j, so we replace the constant r_j with the variable r'. Based on Equation (17), the dependency on r' is of the form $r^{\prime} {E}_{1}(\tau (r^{\prime} ))$. This can be integrated analytically using the relation between an exponential integral and its derivative, $\int {E}_{{\rm{n}}}(x){dx}=-{E}_{{\rm{n}}+1}(x)+C$.

We express r' as a function of τ, assume $\tfrac{{dr}^{\prime} }{d\tau }=$ const., and perform integration by parts:

$$\begin{eqnarray}&&\displaystyle \int r^{\prime} (\tau ){E}_{1}(\tau )\displaystyle \frac{{dr}^{\prime} }{d\tau }d\tau \\ &&\quad =\displaystyle \frac{{dr}^{\prime} }{d\tau }\left(-r^{\prime} {E}_{2}+\displaystyle \frac{{dr}^{\prime} }{d\tau }\displaystyle \int {E}_{2}d\tau \right)+C\\ &&\quad =\displaystyle \frac{{dr}^{\prime} }{d\tau }\left(-r^{\prime} {E}_{2}-\displaystyle \frac{{dr}^{\prime} }{d\tau }{E}_{3}\right)+C.\end{eqnarray} \tag{ 18 }$$

The function $r^{\prime} (\tau )$ differs depending on whether we are looking at τ_min or τ_max and whether shell i is inside or outside of shell j:

$$\begin{eqnarray}&&r^{\prime} ={r}_{i}\pm \displaystyle \frac{{\tau }_{\min }}{{\overline{\chi }}_{{ij}}}={r}_{i}+\displaystyle \frac{{\tau }_{\max }-2{\overline{\chi }}_{{ij}}{r}_{i}}{{\overline{\chi }}_{{ij}}}.\end{eqnarray} \tag{ 19 }$$

Here ± means + in the case of r_i < r' and – in the case of r_i > r'. In every case, $\tfrac{{dr}^{\prime} }{d\tau }$ is either ${\overline{\chi }}_{{ij}}^{-1}$ or $-{\overline{\chi }}_{{ij}}^{-1}$, validating the above assumption $\tfrac{{dr}^{\prime} }{d\tau }\,=$ const.

To summarize, we have three different relations between r' and τ: ${r}_{-}^{{\prime} }({\tau }_{\min })$, ${r}_{+}^{{\prime} }({\tau }_{\min })$, and $r^{\prime} ({\tau }_{\max })$.

This situation arises because the minimum optical depth (case θ' = 0) either increases or decreases with r' depending on whether r' > r_i or r' < r_i, but the maximum optical depth (case θ' = π) always increases with r'.

We now abandon the assumption of a constant ${\overline{\chi }}_{{ij}}$, which makes it impossible to solve the angle integral analytically. For numerical integration, one may use either variant (7) or (8) of the basic equation. We choose the latter for consistency with the previous section, but the method discussed here is essentially independent of the formulation used.

Recall that our definition of ${A}_{i\leftarrow j}$ is

$$\begin{eqnarray}&&{A}_{i\leftarrow j}={\int }_{{r}_{\mathrm{jL}}}^{{r}_{\mathrm{jR}}}r{{\prime} }^{2}{\int }_{0}^{\pi }\displaystyle \frac{{e}^{-\tau }}{{D}^{2}}\sin (\theta ^{\prime} )d\theta ^{\prime} {dr}^{\prime} .\end{eqnarray} \tag{ 20 }$$

We consider rays emitted from point ${\boldsymbol{r}}$ in directions between $\tilde{\theta }=0$ and π (see again Figure 1 for a visualization of the situation). It is important to note that, unless r = 0, the direction angle $\tilde{\theta }$ is different from the polar angle θ'. Here τ depends on D and $\tilde{\theta }$, so if we have a total of n shells and the number of direction angles considered is a, for each individual frequency, we create 2D arrays of size n ∗ a in which to store the values of τ, D, and θ' for each combination of j and $\tilde{\theta }$. Once that is done, the integral can be calculated numerically. In addition, implementing the numerical integration over the radius of shell j is straightforward: instead of stopping only at each shell center and writing τ, D, and θ' there, we do so multiple times within each shell.

The rest of this section describes how to obtain these arrays. For each direction angle $\tilde{\theta }$, we start from position ${\boldsymbol{r}}$ and then move along the ray from shell to shell until we reach the outer boundary. We will now summarize the algorithm for an individual direction.

First, we define variables representing the current position (initially ${\boldsymbol{r}}$), distance, and optical depth (both initially zero). Then, we move along the ray from one shell to the next, all the while updating distance and optical depth and writing the results in our arrays. We can easily calculate θ' from the position. This continues until reaching the outer boundary. This scheme is similar to the "impact parameter" method introduced by Hummer & Rybicki (1971), in the sense that we move from shell to shell along a ray. In that method, parallel rays at different distances (impact parameters) from the center are considered. The calculation then jumps from one intersection of such a ray with a discrete spherical shell to the next in order to track the changing intensity. This can then be used to calculate the moments of radiation and the Eddington factor. The main difference in our method is the fact that we do not follow the intensity but instead calculate the distance and optical depth between any two shells in a certain direction, in order to obtain our coefficients for the radiative exchange between those two.

To move along the ray, at every step we must define a "target radius" to reach with the next step. For outward-moving rays ($\tilde{\theta }\lt \tfrac{\pi }{2}$), the radial coordinate is always increasing, so the target radius is always larger than the current radius. For inward-moving rays, such as the one depicted in Figure 2, the situation is more complicated, since the radius decreases down to a minimum before increasing again. In addition, the ray crosses smaller shells not once, but either twice or never (the blue and magenta shells in Figure 2).

Zoom In Zoom Out Reset image size

If our current position is $({x}_{c},{y}_{c})$ at radius r_c, we wish to move along the direction of the ray by a distance ΔD until we reach the target radius r_t:

$$\begin{eqnarray}&&{\rm{\Delta }}y=\displaystyle \frac{{\rm{\Delta }}x}{\tan (\tilde{\theta })},\Rightarrow \,{r}_{{\rm{t}}}^{2}={({x}_{c}+{\rm{\Delta }}x)}^{2}+{\left({y}_{c}+\displaystyle \frac{{\rm{\Delta }}x}{\tan (\tilde{\theta })}\right)}^{2}.\end{eqnarray} \tag{ 21 }$$

Defining $g=1/\tan (\tilde{\theta })$, after some manipulation, we obtain the following quadratic equation for Δx:

$$\begin{eqnarray}&&{\rm{\Delta }}{x}^{2}+\displaystyle \frac{2{x}_{c}+2{y}_{c}g}{1+{g}^{2}}{\rm{\Delta }}x+\displaystyle \frac{{r}_{c}^{2}-{r}_{{\rm{t}}}^{2}}{1+{g}^{2}}=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\Rightarrow \,{\rm{\Delta }}{x}_{\mathrm{1,2}}=-\displaystyle \frac{{x}_{c}+{y}_{c}g}{1+{g}^{2}}\pm \sqrt{{\left(\displaystyle \frac{{x}_{c}+{y}_{c}g}{1+{g}^{2}}\right)}^{2}-\displaystyle \frac{{r}_{c}^{2}-{r}_{{\rm{t}}}^{2}}{1+{g}^{2}}}.\end{eqnarray} \tag{ 23 }$$

This equation tells us how far in the x direction we need to move along the ray to reach the target radius. There are four possible cases, whose geometric meaning is depicted in Figure 2.

• 1.  
  r_t > r_c: There are two solutions, one positive and one negative. Since $0\leqslant \tilde{\theta }\leqslant \pi$, all the rays we consider move in the positive x direction, so we choose the positive solution.
• 2.  
  r_t < r_c, and the value under the root is positive: There are two positive solutions. We choose the smaller of the two, since otherwise we would go through the whole inner region in a single step.
• 3.  
  r_t < r_c, and the value under the root is negative: There is no real solution. This means that the target radius is not reached by this ray. We set r_t = r_c and recalculate.
• 4.  
  r_t = r_c: One solution is zero, and the other one is $-(2{x}_{c}+2{y}_{c}g)/(1+{g}^{2})$. We choose the latter.

Note that for outward-moving rays ($\tilde{\theta }\lt \pi /2$), only case 1 occurs, since the radius is always increasing.

We have not yet considered the second term in the basic Equation (8). This "boundary term" is given by

$$\begin{eqnarray} & & \kappa {\displaystyle \int }_{{\rm{\Omega }}}({S}_{{\rm{O}}}-S){e}^{-{\tau }_{{\rm{B}}}}d{\rm{\Omega }}\\ & = & 2\pi \kappa {\displaystyle \int }_{\tilde{\theta }}({S}_{{\rm{O}}}-S){e}^{-{\tau }_{{\rm{B}}}}\sin (\tilde{\theta })d\tilde{\theta }.\end{eqnarray} \tag{ 24 }$$

The index O here means the value at the outer boundary, while the index B means "from radius r to the outer boundary." We shall not consider the possibility of a ϕ-dependent S_O here, so that we can eliminate the azimuthal angle and treat the system as 2D. Including the boundary term in the numerical angle integration (Section 2.3.3) is straightforward: we simply add an additional element to our array for τ at the radius of the outer boundary and integrate numerically.

For the constant-χ method, we must numerically integrate as well, since, even if S_O is isotropic, there is no analytic solution. Here τ_B can be expressed as ${\overline{\chi }}_{{\rm{B}}}{D}_{{\rm{B}}}$, and

$$\begin{eqnarray}&&{D}_{{\rm{B}}}=\sqrt{{r}^{2}+{r}_{{\rm{O}}}^{2}-2{{rr}}_{{\rm{O}}}\cos (\theta ^{\prime} )}.\end{eqnarray} \tag{ 25 }$$

Here the polar angle θ' appears again, so we require a relation between θ' and $\tilde{\theta }$.

Here $\tilde{\theta }$ is simply the polar angle of a coordinate system centered at ${\boldsymbol{r}}$, and D is equivalent to the radial coordinate in that system. Since ${\boldsymbol{r}}$ lies on the y axis, for any point (x', y'),

$$\begin{eqnarray}&&x^{\prime} =D\sin (\tilde{\theta }),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&y^{\prime} =D\cos (\tilde{\theta })+r,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\Rightarrow \,\tilde{\theta }=\mathrm{atan}\left(\displaystyle \frac{x^{\prime} }{y^{\prime} -r}\right).\end{eqnarray} \tag{ 28 }$$

Numerical integration is performed over θ'. For each θ', calculate D_B, x', and y', which leads to $\tilde{\theta }$ and, ultimately, ${\rm{\Delta }}\tilde{\theta }(\theta ^{\prime} )$, allowing numerical integration.

While the lack of an analytic solution makes this comparatively slow, it hardly affects the total calculation time because the boundary term only needs to be calculated once per shell, not n times per shell.

Nevertheless, it is possible to find an analytic solution by approximating $\tilde{\theta }=\theta ^{\prime}$. This approximation only holds if r ≪ r_O, but this can be reasonable if, for instance, one is interested in the behavior of a small object inside a large, optically thin envelope.

If S_O does not depend on $\tilde{\theta }$, we can use Equation (14) and $\tau ={\overline{\chi }}_{{\rm{B}}}{D}_{{\rm{B}}}$ to write

$$\begin{eqnarray}&&{\displaystyle \int }_{\tilde{\theta }}{e}^{-{\tau }_{{\rm{B}}}}\sin (\theta ^{\prime} )d\theta ^{\prime} =\displaystyle \frac{1}{{{rr}}_{{\rm{B}}}{\overline{\chi }}_{{\rm{B}}}^{2}}{\displaystyle \int }_{{\tau }_{\mathrm{Bmin}}}^{{\tau }_{\mathrm{Bmax}}}{e}^{-{\tau }_{{\rm{B}}}}d{\tau }_{{\rm{B}}}\\ &&\quad =\displaystyle \frac{1}{{{rr}}_{{\rm{O}}}{\overline{\chi }}_{{\rm{B}}}^{2}}((1+{\tau }_{\mathrm{Bmin}}){e}^{-{\tau }_{\mathrm{Bmin}}}-(1+{\tau }_{\mathrm{Bmax}}){e}^{-{\tau }_{\mathrm{Bmax}}}).\end{eqnarray} \tag{ 29 }$$

We begin with the simplest possible case: a static, fully homogeneous medium with constant temperature. Of course, the radiation term is zero everywhere in this case, but nevertheless, we can assess the accuracy of our method by looking at the radiative cooling only. Even in more complicated cases, the radiative cooling is known to be 4πκ S, giving us an analytic result to compare our method to. Unfortunately, an analytic solution for the heating and therefore the total radiation term is generally not available, except in this trivial case and the thermal relaxation test (Section 3.1.2).

The accuracy of our method depends on the optical depth structure of the system. In practice, this is different for every frequency, but for our test calculations, it is sufficient to look at a single optical depth array. We therefore employ the gray approximation, using the Planck opacity and expressing S as σ T⁴/π.

Instead of the radiation term, we calculate the total radiative cooling by simply setting S' and S_O to zero in all the equations. The optical depth per shell and the number of shells determine which of the three terms in our method dominate. We consider optically thick systems with three different values of τ_S (the optical depth per shell in the radial direction).

• 1.  
  τ_S = 0.1. The external term dominates most of the system, but the boundary term is important in the outer region.
• 2.  
  τ_S = 1. The internal and external terms are comparable. The boundary term is negligible except in the outermost shell.
• 3.  
  τ_S = 10. The internal term dominates everywhere.

We discuss the relative error in the cooling rates (Figure 3). As expected, the constant-χ method is nearly exact in this case, since the system's density and opacity are homogeneous. The only inaccuracy comes from the limited angular resolution of the boundary term (recall from Section 2.3.4 that the boundary term is integrated numerically over all angles θ', even for the fast method). This resolution should be chosen for each shell individually, since for shells very close to the outer boundary, the dependency of distance and optical depth on θ' is much stronger than for shells further inside. We do not investigate the effect of this resolution in detail, since its impact on the overall calculation time is very small.

Zoom In Zoom Out Reset image size

Let us now turn to the numerical integration. We investigate the dependency of the error on two resolutions, namely, the angular resolution (number of angles: a) and the radial resolution in shell j (number of points considered in shell j: x).

Consider first the upper two plots. Both have a total optical depth of 20 and a radial size of 1 per shell, but the optical depth per shell and the number of shells differ.

For τ_S = 1, in the lowest-resolution case (a = 20, x = 2), the numerical integration overestimates the result by about 6%. Increasing both a and x helps to reduce this error. In contrast, for τ_S = 0.1, the angular resolution alone determines the accuracy. This is expected, since for optically thin shells, the exponential term e^−τ hardly changes between the inner and outer boundary of the shell, so once individual (sub)shells are optically thin, a further increase of x has no effect.

We further notice that in the case of optically thin shells, the error is not constant but increases toward the system center up to a maximum and then decreases again. In addition, higher angular resolution is required to achieve a similar accuracy as the τ_S = 1 case. To understand this behavior, refer to case 3 in Figure 2. When calculating the exchange between shell i and a smaller shell j, only a subset of all direction angles $\tilde{\theta }$ actually interacts with shell j. This subset becomes smaller as the ratio of r_j to r_i becomes smaller, effectively decreasing the angular resolution, until it eventually becomes zero. In other words, the exchange with shells smaller than i becomes first inaccurate and then neglected, causing an underestimation of ${A}_{i\leftarrow j}$. This also explains why the error decreases again near the center, since shells there mainly interact with shells larger than themselves. In the case of τ_S = 1, this error does not appear, since only the first few neighboring shells contribute, and their radius is very similar to r_i.

The magnitude of this error is limited by the fact that smaller shells also have a smaller surface area, decreasing their overall contribution, but care must be taken: just because some shell j is negligible for the cooling (i.e., it only absorbs a very small amount of the radiation emitted by shell i) does not necessarily mean that it is negligible for the heating as well, since the shell may be very hot. In such a case, it must be ensured that the angular resolution is high enough to capture all shells that are relevant for the heating.

Finally, the case of τ_S = 10 shows only a tiny error, even with very low angular and radial resolution. This is understandable, since the internal term, which is the only term that matters in this case, depends on neither of these. It does not matter how accurately the external and boundary terms are calculated because they are negligible anyway.

A small temperature perturbation in an otherwise homogeneous system will be smoothed out over time by radiation ("thermal relaxation"), and in this special case, an analytic solution can be found. This was first done by Spiegel (1957) for a plane wave, and Masunaga et al. (1998) showed that the thermal relaxation rate for a spherical wave is exactly the same, provided that the initial perturbation has the form of the zeroth-order spherical Bessel function ${j}_{0}(r)=\sin ({kr})/({kr})$ (see Figure 4 for the concept). This was also confirmed by Stamatellos et al. (2007), in which they applied their radiative transfer scheme for SPH simulations to a spherical thermal relaxation test.

Zoom In Zoom Out Reset image size

Similarly, we test our method by superimposing small perturbations on a base temperature of 10 K. We then measure the decay rate for various values of k/χ. Figure 5 shows very good agreement with the analytic result, confirming that our method works in the case of a homogeneous system.

Zoom In Zoom Out Reset image size

Consider first the results for the constant-χ method (Figure 7). Recall that for this method, we took the average value of χ between i and j in the radial direction and assumed that this same value could be used for every direction. Since in reality, χ near the center is now much larger than that further outside, we expect this method to significantly overestimate the radiative exchange in the region where the external exchange is important. The plot confirms this; while the result is nearly exact close to the outer boundary, the overestimation increases toward the inside and reaches a maximum near the photosphere. After that, the error decreases again for several reasons. In the optically thick region, directions that are much different from the radial direction play a small role; near the center, χ does actually approach direction independency due to symmetry; and, for the very optically thick cases, the external term becomes negligible compared to the internal term.

Zoom In Zoom Out Reset image size

The maximum error is about 20% in the optically intermediate (τ = 1) case but increases to 35% when the core is fully optically thick. At this point, the core's optical thickness is essentially infinite, so further increases have no more effect.

The error in the cooling rate at the same timesteps as in Section 3.2.1 is shown in Figure 8.

Zoom In Zoom Out Reset image size

As expected, there is no error in the far outer and deep inner regions. Near the photosphere, we see two different kinds of erroneous behavior. For all optical depths, there is an underestimation of the result reminiscent of Figure 3's τ_S = 0.1 case. The magnitude of this error depends on the angular resolution. In addition, for the very optically thick cases, there is also an overestimation behind the photosphere. The magnitude of this error depends on the subshell resolution.

The latter error naturally increases with the optical depth per shell but reaches a maximum at τ_S ≈ 3. At even larger values, the internal term starts to dominate the external term, so the error in the latter becomes irrelevant. The internal term does not suffer from this error, since it integrates over the radius analytically.

Except for the thermal relaxation test (Section 3.1.2), up to this point, we have tried to evaluate our method by comparing the total radiative cooling to the analytic value. As was mentioned in Section 3.2, unlike the radiative cooling, the error in the heating depends not only on the optical depth but also on the temperature profile. This makes it difficult to predict how a given error in the cooling will translate to an error in the total.

Figure 9 compares the total radiation terms for different resolutions and the constant-χ method at various stages in the collapse. We also include a high-resolution calculation for the symmetric version of the method discussed in Section 2.3.

Zoom In Zoom Out Reset image size

We concentrate on the photospheric region, since the differences are due to the external term, which is negligible in the deep interior and far outside. We notice several tendencies.

• 1.  
  The net effect of radiation is a radially outward transfer of energy, as would be expected (net cooling inside, net heating outside).
• 2.  
  The constant-χ method consistently shows a smaller value than the numerical integration. This is understandable, since, as was mentioned briefly in Section 2.3.2, the overestimation of ${A}_{i\leftarrow j}$ is smaller for i > j than for i < j. This leads to the relative overestimation in the heating being smaller than in the cooling, since the former comes mostly from the inside.
• 3.  
  As the total optical depth increases, so does the jump in χ at the boundary of the first core, which becomes more and more sharply defined. Because of this, the number of shells where the external term is relevant continues to decrease. This is advantageous for the computational speed; however, if a higher resolution of the (post-)photospheric region is desired, some kind of mesh refinement method is needed.
• 4.  
  The differences are most noticeable in the τ ≈ 100 case, when there is an optically thick region where individual shells are still optically thin. The error in the fast method is more significant here. Because of this, the pure constant-χ method may not be a good choice for a realistic calculation.
• 5.  
  Concerning the resolution of the numerical integration: Similar to the results in Figure 8, we see that low angular resolution leads to an overall underestimation, while low subshell resolution leads to an overestimation in the optically thick region but has almost no effect further outside. The behavior of the fast method essentially corresponds to a calculation with minimal angular and infinite subshell resolution. The differences between a = 100, x = 4 and a = 500, x = 10 are small, showing that the method has converged.
• 6.  
  The use of the symmetric definition of A_ij has only a small effect for optically thin shells, since there is little difference then between fixing the position in shell i and integrating over shell j or doing it the other way around. It also does not matter for very optically thick shells, where the internal term dominates. However, for the intermediate case (interior region of the middle panel in Figure 9), the result is significantly different from the nonsymmetric case, despite using the same angular and subshell resolutions. Note that the differences here do not necessarily mean that the symmetric version is less accurate; it is rather that the two versions attempt to do different things. The nonsymmetric version is designed to calculate the radiation term at the shell center as accurately as possible, while the symmetric version is intended to calculate a shell average that guarantees energy conservation.
• 7.  
  As would be expected, the flux-limited diffusion approximation causes significant errors in the photospheric region where individual shells are not yet optically thick.

These results are specific for this particular optical depth and temperature profile, so when applying the method to a different system, the results may differ as well. However, the general situation of a hot, optically thick interior surrounded by a cold, thin envelope is a rather common state in astrophysics.

The differences between the symmetric and nonsymmetric versions of the method are expected to disappear in the limit of optically thin shells, since in that case, the variation of the radiation term within a shell should vanish. In principle, we could confirm this by performing a protostellar collapse calculation with sufficiently high resolution to keep all shells optically thin throughout the simulation. In practice, the computational cost makes this impossible: the system optical depth after first core formation is already larger than 10³, so the required number of shells would be 10⁴ at the very least and 10⁵ for truly optically thin shells. The presence of optically thick shells is not an issue, since they are dominated by the internal term, which is the same for both variants of the methods. However, there will always be a region in which shells are neither fully thin nor fully thick. This is the problematic region in which the external term is important. Increasing the resolution will not make this region disappear but only move it closer to the system center.

Therefore, we test the convergence of the two variants by using a simplified system in which the density is constant, the velocity is zero, and all hydrodynamic calculations are switched off. We introduce a temperature gradient (temperature decreasing outward from the center) to create a situation similar to the realistic simulation, and the (gray) opacity increases quadratically with temperature. The latter is important to introduce some asymmetry between the different shells. If the density, opacity, and radial size of all shells are the same, then the asymmetry that causes the differences between the two variants of the method (Section 2.3.1) disappears, and both variants give identical results irrespective of the temperature profile.

Figure 10 clearly shows that the two variants, whose results differ substantially in the low-resolution case (optically intermediate shells), are converged for sufficiently high resolution (optically thin shells). Interestingly, the high-resolution results are closer to the low-resolution results of the symmetric than the nonsymmetric method, which indicates that, at least for this setup, the former may be preferable. The flux-limited diffusion results are also included, but as would be expected in this region of optically thin/intermediate shells, they are quite different from the results of our method.

Zoom In Zoom Out Reset image size

In Section 3.2.3, we found that the constant-χ method is too inaccurate to be used for the whole calculation, at least during the time when there is an optically thick region with optically thin shells. The numerical angle integration is much more accurate but also much slower. In this section, we will attempt to combine the two in order to arrive at a compromise between speed and accuracy.

This can be done in the following way. At timestep t, both the slow and fast methods are used to calculate the ${A}_{i\leftarrow j}$. Then, the ratio of the two is written into an array of correction factors ${C}_{i\leftarrow j}={A}_{i\leftarrow j}(\mathrm{fast})/{A}_{i\leftarrow j}(\mathrm{slow})$. For the following timesteps, only the fast method is used, and the resulting ${A}_{i\leftarrow j}$ are divided by ${C}_{i\leftarrow j}$. After n timesteps, both methods are used and the correction factors recalibrated. From Figure 11, we can see that after five steps, the results have already become quite inaccurate, so this combined method could only reduce the computational cost by a factor of a few at most. However, this is not true over the whole calculation. Figure 11 represents the worst case, as it is during the transition phase when the optical depth structure is changing quickly. In earlier and later stages, acceptable values of n can be much larger. Therefore, instead of being constant, n should be changed dynamically during runtime. We suggest the following method.

Zoom In Zoom Out Reset image size

First, we choose a number m (for example, m = 0 initially) and set n = 2^m. Every n timestep, we perform both slow and fast calculations as described above to obtain the radiation terms dU/dt_fast and dU/dt_slow. We then calculate the relative error in the internal energy that would result after n timesteps:

$$\begin{eqnarray}&&\delta U=\displaystyle \frac{| {dU}/{{dt}}_{\mathrm{fast}}-{dU}/{{dt}}_{\mathrm{slow}}| n{\rm{\Delta }}t}{U}.\end{eqnarray} \tag{ 37 }$$

We evaluate this for every shell, and if it is larger than some critical value in any shell, we reduce m by 1. If it is smaller than some other critical value in every shell, we increase m by 1. In this way, the frequency of slow and accurate radiative calculations changes dynamically during the calculation. This method still risks having large errors in some timesteps when m is too large before it can adapt. We can avoid this by including a rollback function as a safety. Every time we perform a slow calculation, we also store all the basic variable arrays, such as density, temperature, etc. If the error in the next slow timestep (n timesteps later) is too large, in addition to reducing m, we also jump back to the stored values. This guarantees that the maximum error is always smaller than some critical value defined by the user. In practice, this method succeeds in significantly increasing the computational speed compared to the pure slow method, since numerical calculations are only very rarely needed after the first core has established itself and the optical depth structure no longer changes quickly.

Table 1 shows the wall-clock time for calculations using the pure fast or pure slow methods. In our calculation, there are three different constraints on the timestep.

• 1.  
  The basic hydrodynamical timestep, determined by the condition that neighboring shells must not overtake each other, i.e., Δt = xΔr/(Δv+c_s), where Δr and Δv are the differences in radius and velocity of two neighboring shells, and x is a number between 0 and 1. Here we choose x = 0.1.
• 2.  
  The radiation timestep. Since the radiation term is very sensitive to temperature, it must be recalculated after the temperature in any shell has changed by more than a certain percentage. Here we choose 1%. Note that the exchange coefficients A_ij are not updated at this point.
• 3.  
  The opacity timestep. The A_ij depend on the optical depth structure of the system, which in turn is determined mostly by the density structure. We update these coefficients whenever the density in any shell has changed by more than 1% since the last update. Whenever this is done, the radiation term is also recalculated.

The purely hydrodynamical timestep requires much less computational effort than the ones related to the radiative transfer, as evidenced by the fact that the calculation without radiative transfer (RT) takes less than a minute. The radiation timestep and especially the opacity timestep are far more time-consuming, the latter even more so if numerical angle integration is used. The calculations shown in Table 1 required roughly 7500–10,000 hydrodynamical timesteps (depending on the number of shells), 3000 radiation timesteps without A_ij updates, and 1500–2000 opacity timesteps. Note that the ratios of these numbers may be vastly different in other systems, depending on the timescales at play. The flux-limited diffusion method is comparable in speed to the fast method at 200 shells but scales better with higher resolutions, since it is $O({fn})$. Its speed depends on how many steps are allowed for the radiation field to converge, however.

As mentioned before, the external term is by far the slowest of the three terms in the radiative calculation. If the number of shells is n and the number of frequencies considered is f, the fast method's external term scales as $O({{fn}}^{2})$, while the numerical integration over a angles and with x subshells scales as $O({{axfn}}^{2})$, which for reasonable resolutions is two to three orders of magnitude slower. This is why optimization requires minimizing the number of numerical integrations, as we did in the previous section.

The separation of the computationally expensive external term from the other two, and the fact that it can be ignored in many regions, is of primary importance. In order to optimize speed, it is imperative to neglect the external term wherever possible, especially in the slow method. In our collapse simulations, we may ignore the external term in the far outer regions and, after first core formation, also in the interior, because individual shells there are optically thick. In addition, there may be regions where the external term is not negligible but high accuracy is not required. In this case, instead of completely neglecting it, one may calculate it using only the fast method. These optimizations can reduce the actual time required by an order of magnitude or more compared to the values seen in Table 1. In any case, the optimization strategy must be adapted according to the problem.