11institutetext:  Institute for Astronomy and Astrophysics, Section
Computational Physics, Eberhard Karls Universität Tübingen, Auf der
Morgenstelle 10, D-72076 Tübingen, Germany
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Dipartimento di Fisica, Università degli Studi di Torino, via
Pietro Giuria 1, 10125 Torino, Italy

# Radiation hydrodynamics integrated in the code PLUTO

Stefan M. Kolb 11 Matthias Stute 11 Wilhelm Kley 11 Andrea Mignone 22

(Received ; accepted 9. September 2013)

###### Abstract

Aims. The transport of energy through radiation is very important in many
astrophysical phenomena. In dynamical problems the time-dependent equations of
radiation hydrodynamics have to be solved. We present a newly developed
radiation-hydrodynamics module specifically designed for the versatile MHD
code PLUTO .

Methods. The solver is based on the flux-limited diffusion approximation in
the two-temperature approach. All equations are solved in the co-moving frame
in the frequency independent (grey) approximation. The hydrodynamics is solved
by the different Godunov schemes implemented in PLUTO , and for the radiation
transport we use a fully implicit scheme. The resulting system of linear
equations is solved either using the successive over-relaxation (SOR) method
(for testing purposes), or matrix solvers that are available in the PETSc
library. We state in detail the methodology and describe several test cases in
order to verify the correctness of our implementation. The solver works in
standard coordinate systems, such as Cartesian, cylindrical and spherical, and
also for non-equidistant grids.

Results. We have presented a new radiation-hydrodynamics solver coupled to the
MHD-code PLUTO that is a modern, versatile and efficient new module for
treating complex radiation hydrodynamical problems in astrophysics. As test
cases, either purely radiative situations, or full radiation-hydrodynamical
setups (including radiative shocks and convection in accretion discs) have
been studied successfully. The new module scales very well on parallel
computers using MPI. For problems in star or planet formation, we have added
the possibility of irradiation by a central source.

###### Key Words.:

radiation transport – irradiation – hydrodynamics – accretion disc

##  1 Introduction

Radiative effects play a very important role in nearly all astrophysical fluid
flows, ranging from planet and star formation to the largest structures in the
universe. Coupling the equations of radiation transport to those of
(magneto-)hydrodynamics (MHD) has been studied for decades, and comprehensive
treatments can be found for example in text books by Mihalas & Mihalas (1984)
or Pomraning (1973). The numerical implementation of two-temperature radiation
hydrodynamics (in the diffusion approximation) into multi-dimensional MHD/HD-
codes has been done already over twenty years ago in various implementations,
for example by Eggum et al. (1988), Kley (1989), in the ZEUS-code (Stone et
al., 1992), and later by Turner & Stone (2001).

In order to study, for example, the dynamics and characteristics of stellar
atmospheres together with convection, more accurate solvers for the radiation
transport based on the method of short characteristics have been developed,
see Davis et al. (2012) and Freytag et al. (2012) for the present status. This
can then be coupled to the hydrodynamics using the Variable Eddington Tensor
method (Jiang et al., 2012). Another approach is the $\mathrm{M}1$ closure
model where the radiative moment equations are closed at a higher level
(González et al., 2007; Aubert & Teyssier, 2008). Despite this progress it is
still useful and desirable to have a method at hand which solves the
interaction of matter and radiation primarily within the bulk part of the
matter which may be optically thick. In such type of applications, the method
of flux-limited diffusion (FLD, see Levermore & Pomraning (1981)) has its
clear merits and is still implemented into existing MHD-codes, for example in
NIRVANA (Kley et al., 2009) to study the planet formation process, in RAMSES
(Commerçon et al., 2011) for protostellar collapse simulations, and in
combination with a multi-frequency irradiation tool into PLUTO (Kuiper et al.,
2010) for massive star formation.

Since the 3D-MHD code PLUTO (Mignone et al., 2007) is becoming increasingly
popular within the computational astrophysics community, we added a publicly
available radiation module, which is based on the two-temperature FLD-
approximation, as described by Commerçon et al. (2011). PLUTO solves the
equations of hydrodynamics and magnetohydrodynamics including the non-ideal
effects of viscosity, thermal conduction and resistivity by means of shock-
capturing Godunov-type methods. Several Riemann solvers, several time stepping
methods and interpolation schemes can be chosen. Additionally, we added a ray-
tracing routine that allows for additional irradiation by a point source in
the center. Treating the irradiation in a ray-tracing approach, guarantees the
long-range character of the radiation better than FLD (Kuiper et al., 2012;
Kuiper & Klessen, 2013).

The paper is organized as follows. In section 2.1, we briefly introduce the
equations of hydrodynamics including radiation transport. Additionally we
describe the general idea behind the flux limited diffusion approximation. In
section 3, we present the discretization of the equations and the solver of
the resulting matrix equation, and present our numerical implementation of
irradiation. In section 4, we present six different test cases to show the
correctness of the implemented equations: four test cases with an analytical
solution (section 4.1 to 4.4) and two others in which our results are compared
with those from other codes (section 4.5 and 4.6). We end with a summary and
conclusions.

##  2 Radiation hydrodynamics

###  2.1 The equations

Even though the PLUTO-environment includes the full MHD-equations and non-
ideal effects such as viscosity, we restrict ourselves here to the Euler
equations of ideal hydrodynamics. Radiation effects are included in the two-
temperature approximation, which implies an additional equation for the
radiation energy. In order to follow the transport of radiation, we apply the
flux-limited diffusion approximation and treat the exchange of energy and
momentum between the gas and the radiation field with additional terms in the
gas momentum and energy equations. The system of equations then read:

$\displaystyle\frac{\partial}{\partial t}\rho+\nabla\cdot(\rho\@vec{v})$
$\displaystyle=0$ (1) $\displaystyle\frac{\partial}{\partial
t}\rho\@vec{v}+\nabla\cdot\left(\rho\@vec{v}\otimes\@vec{v}\right)+\nabla p$
$\displaystyle=+\rho(\@vec{a}_{\mathrm{ext}}+\@vec{a}_{\mathrm{rad}})$ (2)
$\displaystyle\frac{\partial}{\partial
t}e+\nabla\cdot\left[\left(e+p\right)\@vec{v}\right]$
$\displaystyle=+\rho\@vec{v}\cdot(\@vec{a}_{\mathrm{ext}}+\@vec{a}_{\mathrm{rad}})$
(3) $\displaystyle\phantom{=}\hskip 3.0pt-\kappa_{\mathrm{P}}\rho
c(a_{\mathrm{R}}T^{4}-E)$ $\displaystyle\frac{\partial}{\partial
t}E+\nabla\cdot\@vec{F}$ $\displaystyle=\kappa_{\mathrm{P}}\rho
c\left(a_{\mathrm{R}}T^{4}-E\right)$ (4)

The first three equations (1-3) describe the evolution of the gas motion,
where $\rho$ is the gas density, $p$ the thermal pressure, $\@vec{v}$ the
velocity, $e=\rho\,\epsilon+1/2\,\rho\,v^{2}$ the total energy density (i.e.,
the sum of internal and kinetic) of the gas without radiation, and
$\@vec{a}_{\mathrm{ext}}$ an acceleration caused by external forces (e.g.
gravity), not induced by the radiation field (see below). This system of
equations is closed by the ideal gas relation

$p=(\gamma-1)\,\rho\,\epsilon=\rho\frac{k_{\mathrm{B}}\,T}{\mu\,m_{\mathrm{H}}},$
(5)

where $\gamma$ is the ratio of specific heats, $T$ the gas temperature,
$k_{\mathrm{B}}$ the Boltzmann constant, $\mu$ the mean molecular weight, and
$m_{\mathrm{H}}$ the mass of hydrogen. The specific internal energy can be
written as $\epsilon=c_{\mathrm{V}}\,T$, with the specific heat capacity given
by

$c_{\mathrm{V}}=\frac{k_{\mathrm{B}}}{(\gamma-1)\mu m_{\mathrm{H}}}.$ (6)

Here, we assume constant $\gamma$ and $\mu$, which also implies a constant
$c_{\mathrm{V}}$.

The evolution of the radiation energy density $E$ is given by eq. (4), where
$\@vec{F}$ denotes the radiative flux, $\kappa_{\mathrm{P}}$ the Planck mean
opacity, $c$ the speed of light and $a_{\mathrm{R}}$ the radiation constant.
The fluid is influenced by the radiation in two different ways. First, the
radiation may be absorbed or emitted by the fluid leading to variation of its
energy density. This variation is given by the expression
$\kappa_{\mathrm{P}}\rho c\left(a_{\mathrm{R}}T^{4}-E\right)$, see right hand
side of equation (3) and (4). The second effect is that of radiation pressure.
We include this term as an additional acceleration to the momentum equation,
$\@vec{a}_{\mathrm{rad}}=\frac{\kappa_{\mathrm{R}}}{c}\@vec{F}$. The present
implementation does not include the advective transport terms for the
radiation energy and radiative pressure work in eqs. (3) and (4). For the
relatively low temperature protoplanetary disk application that we consider
here these terms are of minor importance. If required, these terms can be
treated in our implementation straightforwardly within PLUTO by adding
additional source terms.

###  2.2 The flux-limited diffusion approximation

The system of equations shown cannot be solved without further assumptions for
the radiative flux $\@vec{F}$. Here we use the flux limited diffusion
approximation (FLD) where the radiation flux is given by a diffusion
approximation

$\@vec{F}=-\lambda\,\frac{c}{\kappa_{\mathrm{R}}\,\rho}\,\nabla E\,,$ (7)

with the Rosseland mean opacity $\kappa_{\rm R}$. The flux-limiter $\lambda$
describes approximately the transition from very optically thick regions with
$\lambda=1/3$ to optically thin regimes, where $\@vec{F}\rightarrow-
cE\frac{\nabla E}{|\nabla E|}$. This leads to the formal definition of the
flux-limiter which is a function of the dimensionless quantity

$R=\frac{|\nabla E|}{\kappa_{\mathrm{R}}\rho E}\,,$ (8)

with the following behaviour:

$\lambda(R)=\left\\{\begin{array}[]{ll}\frac{1}{3},&R\to 0\\\
\frac{1}{R},&R\to\infty\end{array}\right.$ (9)

Physically sensible flux-limiters thus have to fulfil the equation (9) in the
given limits and describe the behaviour between the limits approximately. We
have implemented three different flux-limiters:

$\displaystyle\lambda(R)$ $\displaystyle=\frac{1}{R}\left(\coth
R-\frac{1}{R}\right)$ (10) $\displaystyle\lambda(R)$
$\displaystyle=\left\\{\begin{array}[]{l l}\frac{2}{3+\sqrt{9+12R^{2}}}&\quad
0\leq R\leq\frac{3}{2}\\\
\frac{1}{1+R+\sqrt{1+2R}}&\quad\frac{3}{2}<R\leq\infty\\\ \end{array}\right.$
(13) $\displaystyle\lambda(R)$ $\displaystyle=\left\\{\begin{array}[]{l
l}\frac{2}{3+\sqrt{9+10R^{2}}}&\quad 0\leq R\leq 2\\\
\frac{10}{10R+9+\sqrt{180R+81}}&\quad 2<R\leq\infty\\\ \end{array}\right.$
(16)

from Levermore & Pomraning (1981), Minerbo (1978), and Kley (1989),
respectively. A comparison of them is presented in Kley (1989).

In general it is necessary to solve the equations for each frequency which
appears in the physical problem. However, here we use the grey approximation
in which all radiative quantities including the opacities are integrated over
all frequencies. In our treatment scattering is not accounted for directly,
but it is included in the effective isotropic absorption and emission
coefficients.

##  3 Solving the radiation part

###  3.1 Reformulation of the equations

Instead of solving system of equations (1-4) directly as a whole, the problem
is split into two steps. In the first step, PLUTO is used to solve the
equations of fluid dynamics with the additional force caused by the radiation.
This corresponds to the equations (1) to (3) with the additional acceleration,
$\@vec{a}_{\mathrm{rad}}$, but without the interaction term between the matter
and radiation (last term in eq. 3). By using PLUTO for solving the non-
radiative part of the equations, we are not limited to the Euler equations,
but are able to use the full capabilities of PLUTO for solving the equations
of hydrodynamics or magnetohydrodynamics, including the effects of viscosity
and magnetic resistivity.

In a second, additional step we solve the radiation energy equation (4) and
for the corresponding heating-cooling term in the internal energy of the
fluid:

$\left.\begin{array}[]{lll}\displaystyle{\frac{\partial}{\partial
t}E-\nabla\cdot\left(\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}\nabla
E\right)}&=&\phantom{-}\kappa_{\mathrm{P}}\rho
c\left(a_{\mathrm{R}}T^{4}-E\right)\\\ \displaystyle{\frac{\partial}{\partial
t}\rho\epsilon}&=&-\kappa_{\mathrm{P}}\rho
c\left(a_{\mathrm{R}}T^{4}-E\right)\end{array}\quad\right\\}$ (17)

In order to obtain the radiation energy density, we solve the system of
coupled equations (17). Within one time step PLUTO advances the hydrodynamical
quantities, i.e. the density $\rho$, the velocity $\@vec{v}$ and a temporary
pressure $p$ from time $t^{n}$ to the time $t^{n+1}$, where the time step,
$\Delta t=t^{n+1}-t^{n}$, is determined by PLUTO using the CFL conditions,
presently without including the radiation pressure. These depend on the used
time stepping method in PLUTO , for more information see Mignone et al. (2007)
and the userguide of PLUTO .

The physical process of radiation transport takes place on time scales much
shorter than the one in hydrodynamics. In order to use the same time step for
hydrodynamics and the radiation transport, we apply an implicit scheme to
handle the radiation diffusion and the coupling between matter and radiation
described by equation (17). Because of the coupling of the equations, the
method will update $T$ and $E$ simultaneously, which leads formally to a
nonlinear set of coupled equations. As outlined below, the system is solved
for the radiation energy density $E$. From the new values for $E$, we compute
the new fluid temperature (see eq. 21 below) and update the fluid pressure by
using the ideal gas relation from equation (5). This is then used within PLUTO
to calculate a new total gas energy $e$.

###  3.2 Discretization

In order to discretize the equations (17), we apply a finite volume method.
For that purpose we integrate over the volume of a grid cell and transform the
divergence into a surface integral. Furthermore, we replace the gradient of
$E$ by finite differences, and apply an implicit scheme. The discretization
scheme has been implemented in 3D for Cartesian, cylindrical and spherical
polar coordinates including all the necessary geometry terms for the
divergence and gradient. Since the density has been updated already in the
hydrodynamical part of the solver, we can replace
$\frac{\partial\rho\epsilon}{\partial t}$ with
$\rho\,c_{\mathrm{V}}\,\frac{\partial T}{\partial t}$, which is valid for a
constant heat capacity. Then the resulting discretized equations for the
radiative part can be written as

$\displaystyle\frac{E_{i,j,k}^{n+1}-E_{i,j,k}^{n}}{\Delta t}$ (18)
$\displaystyle=$ $\displaystyle
G_{x1}^{r}K_{i+\frac{1}{2},j,k}^{n}\frac{E_{i+1,j,k}^{n+1}-E_{i,j,k}^{n+1}}{\Delta{x_{1}}_{i+\frac{1}{2}}}-G_{x1}^{l}K_{i-\frac{1}{2},j,k}^{n}\frac{E_{i,j,k}^{n+1}-E_{i-1,j,k}^{n+1}}{\Delta{x_{1}}_{i-\frac{1}{2}}}$
$\displaystyle+$ $\displaystyle
G_{x2}^{r}K_{i,j+\frac{1}{2},k}^{n}\frac{E_{i,j+1,k}^{n+1}-E_{i,j,k}^{n+1}}{\Delta{x_{2}}_{j+\frac{1}{2}}}-G_{x2}^{l}K_{i,j-\frac{1}{2},k}^{n}\frac{E_{i,j,k}^{n+1}-E_{i,j-1,k}^{n+1}}{\Delta{x_{2}}_{j-\frac{1}{2}}}$
$\displaystyle+$ $\displaystyle
G_{x3}^{r}K_{i,j,k+\frac{1}{2}}^{n}\frac{E_{i,j,k+1}^{n+1}-E_{i,j,k}^{n+1}}{\Delta{x_{3}}_{k+\frac{1}{2}}}-G_{x3}^{l}K_{i,j,k-\frac{1}{2}}^{n}\frac{E_{i,j,k}^{n+1}-E_{i,j,k-1}^{n+1}}{\Delta{x_{3}}_{k-\frac{1}{2}}}$
$\displaystyle+$
$\displaystyle{\kappa_{\mathrm{P}}^{n}}_{i,j,k}\,\rho_{i,j,k}^{n}c\left(a_{\mathrm{R}}(T_{i,j,k}^{n+1})^{4}-E_{i,j,k}^{n+1}\right)\,,$

and for the thermal energy (or temperature, respectively)

$\frac{T_{i,j,k}^{n+1}-T_{i,j,k}^{n}}{\Delta
t}=-\frac{{\kappa_{\mathrm{P}}}_{i,j,k}^{n}\,c}{c_{\mathrm{V}}}\,\left(a_{\mathrm{R}}\left(T_{i,j,k}^{n+1}\right)^{4}-E_{i,j,k}^{n+1}\right)\,.$
(19)

Here, the superscript $n$ refers to the values of all variables after the most
recent update from the hydrodynamical step. In order to simplify the notation
for the separate radiation module, we assume the update takes place from time
$n$ to $n+1$. The subscripts $i,j,k$ refer to the 3 spatial directions of the
computational grid, where all variables are located at the cell centers. Half-
integer indices refer to cell interfaces. The physical sizes (proper length)
of each cell in the 3 spatial directions $m$ ($m=1,2,3$) are given by
$\Delta{x_{m}}$, where we additionally allow for non-equidistant grids. The
effective radiative diffusion coefficient (defined at cell centers) is given
by

$K_{i,j,k}^{n}=\frac{c\lambda(R_{i,j,k})}{{\kappa_{\mathrm{R}}}_{i,j,k}^{n}\,\rho_{i,j,k}^{n}}\,,$

where $R_{i,j,k}$ is calculated from eq. (8) by central differencing. Values
at cell interfaces are obtained by linear interpolation. The factors
$G^{l,r}_{xm}$ are geometrical terms defined, respectively, as the left and
right surface areas divided by the cell volume in the direction given by
$m=1,2,3$. In the recent work by Bitsch et al. (2013b) the difference
equations have been written out in more detail for Cartesian, equidistant
grids. The required opacities are evaluated using the values of $\rho$ and $T$
after the hydrodynamical update at time $t^{n}$.

As mentioned before, equations (17) constitute a set of coupled nonlinear
equations. The non-linear term $(T_{i,j,k}^{n+1})^{4}$ that appears in
equation (19) is linearised using the method outlined in Commerçon et al.
(2011)

$(T_{i,j,k}^{n+1})^{4}=(T_{i,j,k}^{n})^{4}\left(1+\frac{T_{i,j,k}^{n+1}-T_{i,j,k}^{n}}{T_{i,j,k}^{n}}\right)^{4}\approx
4(T_{i,j,k}^{n})^{3}T_{i,j,k}^{n+1}-3(T_{i,j,k}^{n})^{4}\,.$ (20)

Using this approximation, we obtain an equation for computing the new
temperature in terms of the new radiation energy density, $E_{i,j,k}^{n+1}$,
and the old temperatures, $T_{i,j,k}^{n}$

$T_{i,j,k}^{n+1}=\frac{{\kappa_{\mathrm{P}}}_{i,j,k}^{n}c\left(3a_{\mathrm{R}}(T_{i,j,k}^{n})^{4}+E_{i,j,k}^{n+1}\right)\Delta
t+c_{\mathrm{V}}T_{i,j,k}^{n}}{c_{\mathrm{V}}+4{\kappa_{\mathrm{P}}}_{i,j,k}^{n}ca_{\mathrm{R}}(T_{i,j,k}^{n})^{3}\Delta
t}\,.$ (21)

The expression can be substituted into eq. (18) to obtain a linear system of
equations for the new radiation energies $E_{i,j,k}^{n+1}$, that can be solved
using standard matrix solvers, see section 3.4. The new temperature can then
be calculated from eq. (21). We implemented several boundary conditions for
the radiation energy density including periodic, symmetric and fixed value.

###  3.3 Irradiation

In order to couple possible irradiation to the radiation transport equations,
a new source term, $S$, has to be added to the right hand side of the thermal
energy equation in system (17)

$\frac{\partial\rho\epsilon}{\partial t}=-\kappa_{\mathrm{P}}\rho
c(a_{\mathrm{R}}T^{4}-E)+S\,.$ (22)

This results in an additional term, $S_{i,j,k}/(\rho_{i,j,k}c_{\rm V})$, in
Eq. 19, correspondingly in equation (21), and in a modification of the right
hand side of the resulting matrix equation for $E_{i,j,k}^{n+1}$.

For the present implementation, we assume that the irradiating source is
located at the centre of a spherical coordinate system. Therefore it is
straightforward to compute the optical depth $\tau_{i,j,k}$ even for
simulations using parallel computers. Assuming that a ray of light travels
along the radial direction from the origin to the grid cell $i,j,k$ under
consideration, the optical depth from the inner radius $r_{0}$ to the $i$th
grid cell with radius $r_{i}$ can be simply expressed as the integral along
the radial coordinate,

$\tau_{i,j,k}=\int_{r_{0}}^{r_{i}}\kappa_{\star}\rho(r)\,dr\approx\sum_{n=0}^{i}\kappa_{\star}{{}_{n,j,k}}\,\rho_{n,j,k}\Delta
r_{n}$ (23)

where $\Delta r_{n}$ is the radial length of the $n$th grid cell, and
$\kappa_{\star}$ the opacity used for irradiation. For the sake of
readability, we write $\tau_{i}$ instead of $\tau_{i,j,k}$ in the following.
We use $\kappa_{\star}=\kappa_{\rm P}$ in the test case with irradiation
presented in section 4.3. Additionally $\kappa_{\star}$ can be defined by the
user as well as the other opacities. Re-emission of the photons which were
absorbed in the cell volume is handled in our treatment by the heating-cooling
term see equation (17).

The luminosity of the source is given by

$L_{\star}=4\pi R_{\star}^{2}\sigma T_{\star}^{4}\,,$ (24)

where $\sigma$ denotes the Stefan-Boltzmann constant, $T_{\star}$ is the
temperature of the star and $R_{\star}$ its radius. In order to compute the
amount of irradiated energy which is absorbed by a specific grid cell we have
to know the surface area $A$ of a grid cell oriented perpendicular to the
radiation from the star and the flux $f$ at the radius $r$. This surface area
$A$ is given by the expression

$A_{i,j,k}=\int\limits_{\theta_{j}}^{\theta_{j+1}}\int\limits_{\phi_{k}}^{\phi_{k+1}}\,dA=r_{i}^{2}(\phi_{k+1}-\phi_{k})(\cos\theta_{j}-\cos\theta_{j+1})\,,$
(25)

where $\theta$ is the azimuthal and $\phi$ the polar angle in the spherical
coordinate system. Without absorption the flux $f$ is given by the expression

$f=\frac{L_{\star}}{4\pi r^{2}}=\sigma
T_{\star}^{4}\left(\frac{R_{\star}}{r}\right)^{2}\,.$ (26)

The amount of energy per time which arrives at the surface of the grid cell
$(i,j,k)$ is

$H_{i,j,k}=A_{i,j,k}f=(\phi_{k+1}-\phi_{k})(\cos\theta_{j}-\cos\theta_{j+1})\sigma
T_{\star}^{4}R_{\star}^{2}\,,$ (27)

again without absorption. If the irradiated energy is partly absorbed, the
remaining amount of energy per time is then $H_{i,j,k}e^{-\tau_{i,j,k}}$.
Using these results, we can compute the energy density per time, $S$, which is
absorbed by one grid cell ($i,j,k$)

$\displaystyle S_{i,j,k}$
$\displaystyle=\frac{H_{i,j,k}e^{-\tau_{i}}-H_{i,j,k}e^{-\tau_{i+1}}}{V_{i,j,k}}=\frac{H_{i,j,k}\left(e^{-\tau_{i}}-e^{-\tau_{i+1}}\right)}{V_{i,j,k}}$
$\displaystyle=\frac{3\sigma
T_{\star}^{4}R^{2}_{\star}\left(e^{-\tau_{i}}-e^{-\tau_{i+1}}\right)}{(r_{i+1}^{3}-r_{i}^{3})}\,,$
(28)

with the volume of a grid cell

$\displaystyle V_{i,j,k}$
$\displaystyle=\int\limits_{r_{i}}^{r_{i+1}}\int\limits_{\theta_{j}}^{\theta_{j+1}}\int\limits_{\phi_{k}}^{\phi_{k+1}}r^{2}sin\theta\,dr\,d\theta\,d\phi$
$\displaystyle=\frac{1}{3}(r_{i+1}^{3}-r_{i}^{3})(\cos\theta_{j}-\cos\theta_{j+1})(\phi_{k+1}-\phi_{k})\,.$
(29)

The absorbed energy density per time, $S_{i,j,k}$, is computed for each grid
cell before solving the matrix equation. A similar treatment of irradiation
has been described recently by Bitsch et al. (2013b), for a multi-frequency
implementation see Kuiper et al. (2010).

###  3.4 The matrix solver

We implemented two different solvers for the matrix equation. The first one
uses the method of successive over-relaxation (SOR), and as a faster and more
flexible solver we use the PETSc111For more information visit the website
http://www.mcs.anl.gov/petsc or have a look at Balay et al. (2012). library.
From the PETSc library we use the Krylov subspace iterative method and a
preconditioner to solve the matrix equation. For all test cases described we
used gmres (Generalized Minimal Residual) as iterative methode and bjacobi
(Block Jacobi) as preconditioner. Beside others the convergence of the SOR
algorithm and the PETSc library can be estimated using the following criteria

$\left\|\@vec{r}^{(k)}\right\|<\max(\epsilon_{\mathrm{r}}\cdot\left\|\@vec{b}\right\|,\epsilon_{\mathrm{a}})$
(30)

where $\@vec{b}$ is the right hand side of the matrix equation
$A\@vec{x}=\@vec{b}$, $\@vec{r}^{(k)}=\@vec{b}-A\@vec{x}^{(k)}$ is the
residual vector for the k-th iteration of the solver and $\@vec{x}$ is the
solution vector (here the radiation energy density). As norm we used here the
$L_{2}$ norm. The quantities $\epsilon_{\mathrm{r}}$ and
$\epsilon_{\mathrm{a}}$ are the relative and absolute tolerance, respectively,
and are problem dependent, with a common value of $10^{-50}$ for
$\epsilon_{\mathrm{a}}$. For the test cases in section 4 we use relative
tolerances $\epsilon_{\mathrm{r}}$ between $10^{-5}$ and $10^{-8}$. The
criterion (30) is the default one used by the PETSc library. For more
information about the convergence test in PETSc the reader should refer to
section 4.3.2. of Balay et al. (2012). The solver performance in a parallel
environment is described in section 4.6.4.

##  4 Test cases

In order to verify the implemented method, we simulated several test problems
and compared the results with either corresponding analytical solutions or
calculations done with different numerical codes. Most of the tests correspond
to one-dimensional problems. In order to model those, we have used quasi one-
dimensional domains, with a very long cuboid that has the height $h$, width
$w$ and a length $l$. The length $l$ is much larger than the width or height,
and for simplicity we use $w=h$. We performed some of the tests in all three
implemented coordinate systems (Cartesian, cylindrical and spherical) and in
three different alignments of the cuboid along each coordinate direction. This
is done to check whether the geometry factors are correct. In the case of a
non-Cartesian coordinate system we placed the cuboid at large distances $r$
from the origin such that the domain approximately describes a Cartesian
setup.

We use for all test cases the solver based on the PETSc library with the
default iterative solver gmres and the pre-conditioner bjacobi.

###  4.1 Linear diffusion test

Figure 1: Linear diffusion test at the time $t=4.2\cdot 10^{-12}\,\mathrm{s}$.
The simulated (read dots) and the analytical (black line) solution is plotted.
We also plot the solution with the absolute value of the relative error (blue
dashed line) which belongs to the axis on the right. Figure 2:  Time evolution
for the linear diffusion test from time $t=0\,\mathrm{s}$ to $4.2\cdot
10^{-12}\,\mathrm{s}$ at three different positions at $x=0\,\mathrm{cm}$
(black lines), at $x=0.5\,\mathrm{cm}$ (blue lines) and at
$x=1.0\,\mathrm{cm}$ (red lines). The dotted lines for each position belong to
the simulated solution, the solid lines to the analytical solution and the
dashed lines show the relative error which belong to the axis on the right.

The following test is adapted from Commerçon et al. (2011). The initial
profile of the radiation energy density is set to a delta function which is
then evolved in time and compared to the analytical one dimensional solution.
We perform this test in all implemented coordinate systems (Cartesian,
cylindrical and spherical coordinates) as described above, which results in
nine different simulations. The used domain is quasi one-dimensional and the
equations of hydrodynamics are not solved in this test. Only the radiation
diffusion equation

$\frac{\partial}{\partial
t}E=\nabla\cdot\left(\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}\nabla E\right)$
(31)

is solved which we obtain from equations (17) by setting
$\kappa_{\mathrm{P}}=0$. An analytical solution to equation (31) can be
calculated in the one dimensional case with a constant flux-limiter
$\lambda=\frac{1}{3}$ and a constant product of the Rosseland opacity and
density, here we set $\kappa_{\mathrm{R}}\rho=1\,\mathrm{cm^{-1}}$. The
equation to solve is then given by

$\frac{\partial}{\partial t}E(x,t)=\frac{c}{3}\frac{\partial^{2}}{\partial
x^{2}}E(x,t)$ (32)

with solution

$E(x,t)=\frac{\tilde{E}_{0}}{\sqrt{\frac{4}{3}c\pi
t}}e^{-\frac{3x^{2}}{4ct}}\,,$ (33)

where $\tilde{E}_{0}$ is the integral over the initial profile of the energy
density, $E(x,t=0)$. Note that in the quasi one-dimensional case (using a
stretched 3D domain) $\tilde{E}_{0}$ has the units $\mathrm{erg\,cm^{-2}}$.

####  4.1.1 Setup

The domain is a cuboid with a length of $4\,\mathrm{cm}$ and a width and
height of $0.04\,\mathrm{cm}$. We used here $301\times 3\times 3$ grid cells.
The initial profile of the radiation energy density in the quasi one-
dimensional case is set by

$E_{i}=\begin{cases}1\,\mathrm{\frac{erg}{cm^{3}}},&i=1,2,\dots,N\text{ with
}i\neq\frac{N}{2}\\\\[5.0pt] \hskip 10.0pt\frac{\tilde{E}_{0}}{\Delta
x},&i=\frac{N}{2}\\\ \end{cases}$ (34)

where $\Delta x$ is the length of a grid cell. For numerical reasons, we have
set $E_{i}$ for $i\neq\frac{N}{2}$ to the value $1\,\mathrm{erg\,cm^{-3}}$
instead of $0\,\mathrm{erg\,cm^{-3}}$. This choice is not problematic, since
$\tilde{E}_{0}/\Delta x\gg 1\,\mathrm{erg\,cm^{-3}}$ for our chosen value of
$\tilde{E}_{0}=10^{5}\,\mathrm{erg\,cm^{-2}}$. The initial values for pressure
and density are $p=1\,\mathrm{g\,cm^{-1}\,s^{-2}}$ and
$\rho=1\,\mathrm{g\,cm^{-3}}$. Furthermore we use
$\kappa_{\mathrm{R}}=1\,\mathrm{cm^{2}\,g^{-1}}$ for the Rosseland opacity.
All boundary conditions are set to periodic except for the boundary conditions
at the beginning and end of the quasi one-dimensional domain, which are set to
outflow. For the matrix solver we used a relative tolerance of
$\epsilon_{r}=10^{-8}$. The simulation starts at $t=0\,\mathrm{s}$ with an
constant time step of $\Delta t=1\cdot 10^{-14}\,\mathrm{s}$ and stops at
$t=4.2\cdot 10^{-12}\,\mathrm{s}$.

####  4.1.2 Results

The numerical solution $E_{\mathrm{n}}$ and the analytical solution
$E_{\mathrm{a}}$ from equation (33) are plotted in the figures 1 and 2
together with the absolute value of the relative error. In figure 1 the
radiation energy density is plotted against the position at the time
$t=4.2\cdot 10^{-12}\,\mathrm{s}$. The relative error in the relevant range
from $-1\,\mathrm{cm}$ to $1\,\mathrm{cm}$ is always below one percent. In
figure 2 the time evolution from $t=0\,\mathrm{s}$ to $4.2\cdot
10^{-12}\,\mathrm{s}$ is shown for the positions
$x=\\{0,0.5,1.0\\}\,\mathrm{cm}$ coded in the colors black, blue and red,
respectively. The results shown in this figure depend strongly on the
position. For the position $x=0\,\mathrm{c}m$ the error is, for all times
later than $t=4\cdot 10^{-13}\,\mathrm{s}$, below one percent and decreases
with time. For the other positions, the behaviour is different. The relative
error rises and after a while it decreases. This behaviour can be explained by
looking at figure 1. The error is higher at the diffusion front. This region
moves with time and causes the effect for the other positions. The test shows
that the time evolution of the radiation energy density is reproduced
correctly. As described, this test was performed in different coordinate
systems and orientations, with the same results.

###  4.2 Coupling test

Figure 3: Coupling test from $t=10^{-20}\,\mathrm{s}$ to
$t=10^{-4}\,\mathrm{s}$ with three different initial gas energy densities. The
reference solution (black lines) and the simulated results for the initial
energy density $e_{0}=10^{10}\,\mathrm{erg\,cm^{-3}}$ (red dots),
$e_{0}=10^{6}\,\mathrm{erg\,cm^{-3}}$ (blue dots) and
$e_{0}=10^{2}\,\mathrm{erg\,cm^{-3}}$ (green dots) are plotted.

The purpose of this test from Turner & Stone (2001) is to check the coupling
between radiation and the fluid. For this purpose we simulate a stationary
fluid which is initially out of thermal equilibrium. In this simulation the
radiation energy density is the dominant energy which is constant over the
whole simulation. The system of equations (17) decouples in this case and, in
addition, it is not necessary to solve the matrix equation for $E$. By setting
$\sigma_{\mathrm{P}}=\kappa_{\mathrm{P}}\rho$ and $T=\frac{p}{\rho}\frac{\mu
m_{\mathrm{H}}}{k_{\mathrm{B}}}$ from eq. (5) with the assumption that
$\sigma_{\mathrm{P}}$ and $\rho$ are constant, we can rewrite the thermal
energy equation of the system (17) as

$\frac{\,\mathrm{d}e}{\,\mathrm{d}t}=\underbrace{c\sigma_{\mathrm{P}}E}_{C_{\mathrm{1}}}-\underbrace{c\sigma_{\mathrm{P}}a_{\mathrm{R}}\left(\frac{\gamma-1}{\rho}\frac{\mu
m_{\mathrm{H}}}{k_{\mathrm{B}}}\right)^{4}}_{C_{\mathrm{2}}}e^{4}\,.$ (35)

With the used approximations, the coefficients $C_{\mathrm{1}}$ and
$C_{\mathrm{2}}$ are constant. The solution to Eq. (35) can be calculated
analytically in terms of an algebraic equation which would have to be solved
iteratively. Hence, we integrate Eq. (35) numerically using a Runge Kutta
solver of $4$-th order scheme with adaptive step size. In the following we
refer to this solution as the reference solution.

Information on the expected behaviour of the solution can be obtained directly
from the differential equation. It is clear that in the final equilibrium
state (with $\frac{\,\mathrm{d}e}{\,\mathrm{d}t}=0$) the gas temperature has
to be equal to the radiation temperature
$T=\sqrt[4]{\frac{E}{a_{\mathrm{R}}}}$, thus the final gas energy density will
be

$e_{\mathrm{final}}=\left(\frac{C_{\mathrm{1}}}{C_{\mathrm{2}}}\right)^{\frac{1}{4}}\,.$
(36)

If the initial gas energy density $e_{0}$ is much lower than
$e_{\mathrm{final}}$, we can neglect the second term in eq. (35) at the
beginning, thus $e(t)=C_{\mathrm{1}}\,t+e_{0}$. The corresponding coupling
time can be estimated to

$\tau=\frac{e_{\mathrm{final}}-e_{0}}{C_{\mathrm{1}}}\,.$ (37)

On the other hand, if $e_{0}\gg e_{\mathrm{final}}$, we can neglect the first
term in eq. (35) and derive

$e(t)\propto\left(C_{\mathrm{2}}\,t\right)^{-\frac{1}{3}}\quad\mathrm{and}\quad\tau=\frac{1}{e_{\mathrm{final}}^{3}\,C_{\mathrm{2}}}\,.$
(38)

####  4.2.1 Setup

The computational domain is identical to that of the linear diffusion test in
section 4.1. For the grid we use a resolution of $25\times 3\times 3$ grid
cells. As before we do not solve the equations of hydrodynamics and the
boundary conditions are quite simple. All boundaries are set to periodic
boundary conditions. The constants we used are set to: radiation energy
density $E=10^{12}\,\mathrm{erg\,cm^{-3}}$, density
$\rho=10^{-7}\,\mathrm{g\,cm^{-3}}$, opacity $\sigma_{\mathrm{P}}=4\cdot
10^{-8}\,\mathrm{cm^{-1}}$, mean molecular weight $\mu=0.6$ and the ratio of
specific heats $\gamma=5/3$. The simulations starts at $t=0\,\mathrm{s}$ with
an initial time step of $\Delta t=10^{-20}\,\mathrm{s}$ and evolves until
$t=10^{-4}\,\mathrm{s}$. After each step the time step is increased by 1% in
order to speed-up the computation. The simulation is done with three different
initial gas energy densities, $e_{0}=10^{10}\,\mathrm{erg\,cm^{-3}}$,
$e_{0}=10^{6}\,\mathrm{erg\,cm^{-3}}$ and
$e_{0}=10^{2}\,\mathrm{erg\,cm^{-3}}$.

####  4.2.2 Results

Figure 3 shows the numerical gas energy density and the reference solution
plotted against time for the three different initial values of $e$. The
agreement of both results is excellent for all initial values. From the figure
we see that in the limit of small and large initial $e_{0}$, we find exactly
the behaviour as predicted by the estimates for eq. (35). The analytic
estimates for the coupling time $\tau$ from equation (37) agrees very well
with our results. The estimate for $e_{0}=10^{2}\,\mathrm{erg\,cm^{-3}}$ is
$\tau=5.88\cdot 10^{-8}\,\mathrm{s}$ and for
$e_{0}=10^{6}\,\mathrm{erg\,cm^{-3}}$ we calculated $\tau=5.78\cdot
10^{-8}\,\mathrm{s}$. In the case that $e_{\mathrm{final}}<e_{0}$ the estimate
in eq. (38) is approximate $\tau=5.88\cdot 10^{-8}\,\mathrm{s}$. We have to
mention here that this test verifies primarily the correctness of equation
(21). As in the linear diffusion test, this test was performed in three
different coordinate systems in different orientations, with the same results.

###  4.3 Coupling test with irradiation

Figure 4: Coupling test with enabled irradiation from $t=10^{-7}\,\mathrm{s}$
to $t=10^{3}\,\mathrm{s}$ at three different distances $d$ from the inner
boundary of the domain. The reference solution (black lines) and the simulated
results for the energy density $e=10^{2}\,\mathrm{erg\,cm^{-3}}$ are plotted
at distance $d=0\,\mathrm{cm}$ (red dots), $d=3\cdot 10^{4}\,\mathrm{cm}$
(blue dots) and $d=3\cdot 10^{5}\,\mathrm{cm}$ (green dots). Figure 5: Radial
dependency of the gas energy density for the coupling test with enabled
irradiation. The gas energy density is plotted at different times,
$t=1.12\cdot 10^{-4}\,\mathrm{s}$ (red dots), $t=1.62\cdot
10^{-2}\,\mathrm{s}$ (blue dots), $t=1.43\,\mathrm{s}$ (green dots),
$t=1.72\cdot 10^{1}\,\mathrm{s}$ (yellow dots) and $t=9.18\cdot
10^{2}\,\mathrm{s}$ (magenta dots). The dots represent the numerical solution
and the solid black lines the reference solution. The position $d$ is again
measured relative to the inner boundary of the quasi one-dimensional domain.

This test is in its basic setup the same as that from section 4.2, but with
irradiation enabled, i.e. equation (22) is solved instead of the second
equation in (17). As described in section 3.3, irradiation is limited to
spherical coordinates, which we use for this test. With the same assumptions
as in section 4.2 , i.e., that $\sigma_{\mathrm{P}}$ and $\rho$ are constant
and with the definitions for $\sigma_{\mathrm{P}}$, $e$, $p$ as well as for
$T$, it is possible to rewrite $S$ from equation (28) to

$S(r)=\frac{3\sigma
T_{\star}^{4}R^{2}_{\star}e^{-\sigma_{\mathrm{P}}\left(r-r_{\mathrm{0}}\right)}\left(1-e^{-\sigma_{\mathrm{P}}\Delta
r}\right)}{\left(r+\Delta r\right)^{3}-r^{3}}\,,$ (39)

and obtain for equation (22)

$\frac{\,\mathrm{d}e}{\,\mathrm{d}t}=\underbrace{S(r)+c\sigma_{\mathrm{P}}E}_{C_{\mathrm{1}}(r)}-\underbrace{c\sigma_{\mathrm{P}}a_{\mathrm{R}}\left(\frac{\gamma-1}{\rho}\frac{\mu
m_{\mathrm{H}}}{k_{\mathrm{B}}}\right)^{4}}_{C_{\mathrm{2}}}e^{4}\,.$ (40)

The reference solution is computed in the same way as before although it now
depends on the distance $r$ from the star. The quasi one-dimensional domain
starts at $r=9000\cdot 10^{5}\,\mathrm{cm}$ and ends at $r=9003\cdot
10^{5}\,\mathrm{cm}$ and we use $300\times 3\times 3$ grid cells. The domain
size in $\theta$ and $\phi$ direction was chosen in a way such that the grid
cells are nearly quadratic. For the simulation we use a constant radiation
energy density of $E=10^{-2}\,\mathrm{erg\,cm^{-3}}$, a density of
$\rho=10^{-5}\,\mathrm{g\,cm^{-3}}$ , a Rosseland opacity of $\kappa_{\rm
R}=10\,\mathrm{cm^{2}\,g^{-1}}$ and a Planck opacity of
$\kappa_{\mathrm{P}}=\kappa_{\mathrm{R}}$ which corresponds to
$\sigma_{\mathrm{P}}=10^{-4}\,\mathrm{cm^{-1}}$. The opacity for the
irradiation $\kappa_{\star}$ is set to $\kappa_{\mathrm{P}}$. For the star,
the temperature was set to $T_{\star}=6000\,\mathrm{K}$ and the radius to
$R_{\star}=8.1\cdot 10^{8}\,\mathrm{cm}$. Additionally we make the assumption
that there is no absorption in the region between the surface of the star and
the inner boundary of the computation domain. Figure 4 shows the gas energy
density plotted against time with an initial gas energy of
$e=10^{2}\,\mathrm{erg\,cm^{-3}}$ at three different positions
$d=0\,\mathrm{cm}$, $d=3\cdot 10^{4}\,\mathrm{cm}$ and $d=3\cdot
10^{5}\,\mathrm{cm}$ where $d$ is measured relative to the inner boundary of
the quasi one-dimensional domain. In figure 5 the radial dependency of the gas
energy density is plotted for the same simulation at five different times. As
expected, the results show that the gas energy density at a time later than
$t=10^{2}\,\mathrm{s}$ becomes constant and depends on the distance from the
star. The simulated and reference solution show a excellent agreement.

###  4.4 A steady state test

Figure 6: Comparison between the numerical (red dots) and the analytical
(black line) solution of the steady state test after $t=1200\,\mathrm{s}$. In
addition the absolute value of the relative error (blue dashed line) is
plotted related to the axis on the right side. Note this axis is logarithmic.

The original version of this test was published in Flaig (2011). We consider a
one-dimensional stationary setup with a given density stratification. In the
steady state, the time derivatives in the equations (17) vanish and the system
is reduced to the following equation for the radiation energy density

$0=\nabla\cdot\,\left(\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}\nabla
E\right)\,.$ (41)

A further reduction is obtained when we rewrite this equation in one dimension
along the z-axis in Cartesian coordinates. The equation is then much simpler
and can be written as

$\frac{\,\mathrm{d}}{\,\mathrm{d}z}\left(\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}\frac{\,\mathrm{d}}{\,\mathrm{d}z}E\right)=0\,.$
(42)

In general the expression $\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}$ is not
known analytically for realistic opacities. In order to circumvent this
problem, we define the effective optical depth
$\tau_{\mathrm{eff}}=\int\,\mathrm{d}\tau_{\mathrm{eff}}=\int_{z_{b}}^{z_{a}}\kappa_{\mathrm{eff}}\rho\,\mathrm{d}z$
where $z_{a}$ and $z_{b}$ are the lower and upper boundaries of the quasi one-
dimensional domain, respectively, and $\kappa_{\mathrm{eff}}$ is the effective
opacity given by
$\kappa_{\mathrm{eff}}=\frac{1}{3}\frac{\kappa_{\mathrm{R}}}{\lambda}$. By
using
$\,\mathrm{d}\tau_{\mathrm{eff}}=\kappa_{\mathrm{eff}}\rho\,\mathrm{d}z$,
equation (42) can be rewritten as:

$\frac{\,\mathrm{d}}{\,\mathrm{d}z}\frac{\,\mathrm{d}}{\,\mathrm{d}\tau_{\mathrm{eff}}}E=0\,.$
(43)

The solution of this equation is then given by

$E=\left(E(\tau_{\mathrm{eff}}=1)-E(\tau_{\mathrm{eff}}=0)\right)\tau_{\mathrm{eff}}+E(\tau_{\mathrm{eff}}=0)$
(44)

where $E(\tau_{\mathrm{eff}}=0)$ and $E(\tau_{\mathrm{eff}}=1)$ are the
radiation energy density at the position where the effective optical depth has
the values zero or one, respectively. Thus, in the static case the radiation
energy has a linear dependence on the optical depth $\tau_{\mathrm{eff}}$ for
all opacity laws.

####  4.4.1 Initial setup

The domain was chosen to have an arbitrary length of $300\,\mathrm{cm}$ and a
width and height of $3\,\mathrm{cm}$ and $300\times 3\times 3$ grid cells were
used. This test is performed without solving the hydrodynamical equations,
instead we solved equations (17) for a fixed density and opacity law, and
evolved the solution, until a stationary state has been reached. For the
radiation boundary conditions, we used boundary conditions with fixed values
of $E$ at the lower and upper boundary of the domain. At the lower boundary we
have chosen $E=a_{\mathrm{R}}T^{4}$ with a temperature of
$T=2000\,\mathrm{K}$.

Because the stratification is optically thin at the upper boundary we want to
allow the radiation to escape freely from the domain. For this reason we
simply set the temperature to a very small value at the upper boundary, here
$T=10\,\mathrm{K}$.

All other boundary conditions have been set to periodic. The density
stratification is given by

$\rho(z)=\rho_{0}e^{\frac{1}{2}\left(\frac{z-z_{a}}{0.46\cdot\left(z_{b}-z_{a}\right)}\right)^{2}}\,.$
(45)

The initial temperature profile can be chosen randomly in principle, but in
order to speed up the computation we used a linear temperature profile
starting at $z_{a}$ with $T=2000\,\mathrm{K}$ and ending at $z_{b}$ with
$T=10\,\mathrm{K}$. From this temperature profile we assigned pressure values
using equation (5). The radiation energy density $E$ inside the domain is also
set using the gas temperature profile and $E=a_{\mathrm{R}}T^{4}$. The ratio
of specific heats and the mean molecular weight are set to $\gamma=1.43$ and
$\mu=0.6$, respectively. As flux-limiter we have chosen equation (13), for the
Rosseland mean opacity $\kappa_{\rm R}$ we use data from Lin & Papaloizou
(1985), and the Planck mean opacity is set to $\kappa_{\rm P}=\kappa_{\rm R}$.
The initial time step is $\Delta t=0.3\,\mathrm{s}$ and it is increased
slightly with time in order to speed-up the computation and to keep the number
of iterations done by the matrix solver nearly constant. This simulation was
preformed with a relative tolerance of $\epsilon_{r}=10^{-6}$ for the matrix
solver.

####  4.4.2 Results

A steady state is reached approximately after $t=1200\,\mathrm{s}$. In figure
6, we plot the radiation energy density against the effective optical depth
$\tau_{\text{eff}}$ from our numerical solution (red dots) together with the
analytical solution from equation (44). The parameters
$E(\tau_{\mathrm{eff}}=1)-E(\tau_{\mathrm{eff}}=0)$ and
$E(\tau_{\mathrm{eff}}=0)$ have been obtained by fitting equation (44) to the
numerical solution. We have to note here that $E(\tau_{\mathrm{eff}}=0)$ is
determined by interpolation between ghost cells and active cells near the
upper boundary $z_{\mathrm{b}}$. Hence, the radiation temperature in the
active region can be much larger than $10\,\mathrm{K}$, a value was
specifically chosen to be very small. We also plot the absolute value of the
relative error $\left|(E_{n}-E_{a})/E_{a}\right|$. The results from the
simulation agree very well with the analytical prediction. As we can see from
figure 6, the largest deviation from the analytical solution is at small
values of $\tau_{\mathrm{eff}}$ with an relative error around one percent. As
the linear diffusion and the coupling test, this test was performed in all
three coordinate systems and in different orientations, with the same results.

###  4.5 Radiation shock

(a) Subcritical shock

(b) Supercritical shock

Figure 7: Sub- and supercritical shock test. In both cases we plot the
radiation temperature (blue line) and the gas temperature (red line) against
$s=z-v\cdot t$ where $z$ is the position along the quasi one-dimensional
domain and $v$ the piston velocity. The subcritical shock 7(a) is shown at
time $t=3.8\cdot 10^{4}\,\mathrm{s}$ and the supercritical shock 7(b) at
$t=7.5\cdot 10^{3}\,\mathrm{s}$.

In this section we extend the previous tests and solve now the full equations
of hydrodynamics and radiation transport simultaneously, testing the complete
new module within PLUTO environment.

####  4.5.1 Initial setup

Following a set-up from Ensman (1994), a shock is generated in a quasi one-
dimensional domain. This test case is more complex than the ones before, and
it is not possible to derive an analytical solution. Instead we compare our
results with the simulations of Commerçon et al. (2011). The computational
domain is chosen to have a length of $7\cdot 10^{10}\,\mathrm{cm}$ and a width
and height of $3.418\cdot 10^{7}\,\mathrm{cm}$ with $2048\times 4\times 4$
grid cells. The initial density and temperature are set to $\rho=7.78\cdot
10^{-10}\,\mathrm{g\,cm^{-3}}$ and $T=10\,\mathrm{K}$. The initial radiation
energy density is set by the equation $E=a_{\mathrm{R}}T^{4}$. For the flux-
limiter we employ the Minerbo-formulation according to eq. (13), and for the
opacity we use
$\kappa_{\mathrm{R}}\cdot\rho=\kappa_{\mathrm{P}}\cdot\rho=3.1\cdot
10^{-10}\,\mathrm{cm^{-1}}$. Furthermore the ratio of specific heats is set to
$\gamma=7/5$ and the mean molecular weight to $\mu=1$, in analogy to Commerçon
et al. (2011). The time step is computed through the CFL condition of PLUTO
for wich we assume a value of $0.4$. For the solver we took the not so
accurate but robust tvdlf which uses a simple Lax-Friedrichs scheme. For
generating the radiative shock, the following boundary conditions are used: in
the direction of the shock propagation, we employ a reflective boundary
condition at the lower boundary and a zero-gradient at the upper boundary of
the domain. The remaining boundaries are set to periodic. For the relative
tolerance used by the matrix solver we have chosen a value of
$\epsilon_{r}=10^{-5}$. The shock is generated by applying an initial velocity
$v$ to the gas. The velocity is directed towards the reflecting boundary
condition which acts as a wall. The shock propagates then from the wall back
into the domain. Depending on the velocity, the shock is sub- or
supercritical, i.e., the temperature behind the shock front is larger or equal
than the temperature upstream (in front of the shock front), respectively. In
this test we simulate both cases: the subcritical shock with a velocity of
$v=6\cdot 10^{5}\,\mathrm{cm\,s^{-1}}$ and the supercritical shock with
$v=20\cdot 10^{5}\,\mathrm{cm\,s^{-1}}$.

####  4.5.2 Results

For a better comparison with the results of simulations, where the material is
at rest and a moving piston causes the shock, we introduce the quantity $s$.
This quantity is given by the relation $s=z-v\cdot t$ where $z$ is the
position along the quasi one-dimensional domain. Note that this quantity is
called $z$ in Commerçon et al. (2011). Fig. 7 shows the radiation temperature
(blue line) and the gas temperature (red line) against the previously defined
quantity $s$ for both the subcritical (at $t=3.8\cdot 10^{4}\,\mathrm{s}$) and
supercritical case (at $t=7.5\cdot 10^{3}\,\mathrm{s}$). In the supercritical
case the pre- and post-shock gas temperature are equal, as expected. In the
subcritical case these temperatures can be estimated analytically (Ensman,
1994; Mihalas & Mihalas, 1984; Commerçon et al., 2011). In table 1, the
analytical estimates and the numerical values from our simulations and the
results from Commerçon et al. (2011) are shown together. Here $T_{2}$ is the
post-shock temperature, $T_{-}$ the pre-shock temperature and $T_{+}$ the
spike temperature. In the equations, $R_{\rm G}=\frac{k_{\mathrm{B}}}{\mu
m_{\mathrm{H}}}$ is the perfect gas constant,
$\sigma_{\mathrm{SB}}=\frac{ca_{\mathrm{R}}}{4}$ the Stefan-Boltzmann
constant, and $u$ is the velocity of the shock relative to the upstream
material (or vice versa) in our case $u=7.19\cdot
10^{5}\,\mathrm{cm\,s^{-1}}$.

|  analytical estimate |  numerical solution |  Commerçon et al.
---|---|---|---
$T_{2}\approx\frac{2(\gamma-1)u^{2}}{R_{\rm G}(\gamma+1)^{2}}$ |  $\sim\phantom{1}865\,\mathrm{K}$ |  $\phantom{1}816.6\,\mathrm{K}$ |  $\phantom{1}825\,\mathrm{K}$
$T_{-}\approx\frac{\gamma-1}{\rho uR_{\rm G}}\frac{2\sigma_{\mathrm{SB}}T_{2}^{4}}{\sqrt{3}}$ |  $\sim\phantom{1}315\,\mathrm{K}$ |  $\phantom{1}331.9\,\mathrm{K}$ |  $\phantom{1}275\,\mathrm{K}$
$T_{+}\approx T_{2}+\frac{3-\gamma}{\gamma+1}T_{-}$ |  $\sim 1075\,\mathrm{K}$ |  $1147.1\,\mathrm{K}$ |  $1038\,\mathrm{K}$
Table 1: Comparison of the results from the radiation shock test with
analytical estimates and the results from Commerçon et al. (2011) for the pre-
shock $T_{-}$ and post-shock $T_{2}$ gas temperature as well as the spike
temperature $T_{+}$.

The results agree in general with the analytical estimates and the results
from Commerçon et al. (2011). The analytical estimate for the post-shock
temperature is higher than the numerical results with both codes. We have to
note here that the analytical estimate depends on $u$ and differs therefore
from the values given in Commerçon et al. (2011). The pre-shock and spike
temperatures agree reasonably well with the analytical estimates in our
simulations but are higher than the results from Commerçon et al. (2011). The
differences of our numerical solution to the analytical estimates might due to
the fact that we ignored the advective terms in the radiation energy density
in eq. (4) that may play a role in this dynamic situation. Additionally, it is
noteworthy that the position of the shock front is very well reproduced. This
test was performed in Cartesian coordinates.

###  4.6 Accretion disc

The goal of this last test is to compare the results of different codes on a
more complex two-dimensional physical problem that involves the onset of
convective motions. For this purpose we model a section of an internally
heated, viscous accretion disc in spherical coordinates $(r,\theta,\phi)$
where $r$ is the distance to the centre of the coordinate system, $\theta$ the
polar angle measured from the $z$-axis in cylindrical coordinates and $\phi$
the azimuth angle. The setup follows the standard disc model used in Kley et
al. (2009). The tests proceed in two steps. In a first setup we reduce the
complexity of the problem and consider a static problem, i.e., without solving
the equations of hydrodynamics. This will demonstrate that the equilibrium
between viscous heating and radiative cooling is treated correctly in our
implementation. In the second setup we consider the full hydrodynamic problem
and study the onset of convection in discs.

####  4.6.1 The initial setup

For both, the static and the dynamical case we use the same initial setup. The
radial extent ranges from $r_{\text{min}}=0.4$ to $r_{\text{max}}=2.5$, where
all lengths are given in units of the semi-major axis of Jupiter
$a_{\text{jup}}=5.2\,\mathrm{AU}$. In the vertical direction the domain
extends from $\theta_{\text{min}}=83\,\mathrm{{}^{\circ}}$ to
$\theta_{\text{max}}=90\,\mathrm{{}^{\circ}}$ and in $\phi$ direction from
$\phi_{\text{min}}=0\,\mathrm{{}^{\circ}}$ to
$\phi_{\text{max}}=360\,\mathrm{{}^{\circ}}$. In the three coordinate
directions $(r,\theta,\phi)$ we use $256\times 32\times 4$ grid cells. The
disc aspect ratio $h$ is set to $h=\frac{H}{s}=0.05\,\mathrm{}$ where
$s=r\sin\theta$ describes the (radial) distance from the $z$-axis in
cylindrical coordinates, and $H$ is the disc’s vertical scale height. The
viscosity $\nu$ is set to a value of $\nu=10^{15}\,\mathrm{cm^{2}\,s^{-1}}$,
and the mean molecular weight to $\mu=2.3\,\mathrm{}$. For the ratio of
specific heats we have used different values, as specified below. The density
stratification can be obtained from vertical hydrostatic equilibrium, assuming
a temperature that is constant on cylinders, $T=T(s)$. It follows (Masset et
al., 2006)

$\rho(r,\theta)=\rho_{0}\cdot
s^{-1.5}\exp\left(\frac{\sin{\theta}-1}{h^{2}}\right)$ (46)

where the quantity $\rho_{0}$ was chosen such that the total mass of the disc
is $M_{\text{disc}}=0.01\cdot M_{\star}$, where $M_{\star}$ is the mass of the
central star of the system which is set to the mass of the sun,
$M_{\star}=M_{\odot}$. The mass within the computational domain is then
$1/2M_{\text{disc}}$ because we only compute the upper half of the disc. The
radial variation leads to a surface density profile of $\Sigma\propto
r^{-1/2}$, which is the equilibrium profile for constant viscosity, and
vanishing mass flux through the disc. The pressure $p$ is set by the
isothermal relation $p=\rho c_{\mathrm{s}}^{2}$, with the speed of sound
$c_{\mathrm{s}}=H\Omega_{\mathrm{K}}$ and the Keplerian angular velocity

$\Omega_{\mathrm{K}}=\sqrt{\frac{GM_{\star}}{s^{3}}}\,,$

with the gravitational constant $G$. The temperature can be computed through
equation (5) and results in $T=\frac{\mu
m_{\mathrm{H}}}{k_{\mathrm{B}}}\frac{p}{\rho}$. The initial velocities are set
to zero except for the angular velocity $v_{\phi}$ which is set to

$v_{\phi}=\sqrt{\frac{\left(1-2h^{2}\right)GM_{\star}}{s}}\,.$

For the Rosseland mean opacity $\kappa_{\rm R}$ we use data from Lin &
Papaloizou (1985), and the Planck mean opacity is set to $\kappa_{\rm
P}=\kappa_{\mathrm{R}}$. The displayed simulations have been performed in the
rotating frame in which the coordinate system rotates with the constant
angular velocity of $\Omega_{\mathrm{K}}$ at $a_{\text{jup}}$, but for non-
rotating systems identical results are obtained. As before the radiation
energy density is initialised to $E=a_{\mathrm{R}}T^{4}$.

For density, pressure and radial velocity we apply reflective radial boundary
conditions and the angular velocity is set to the Keplerian values. In the
azimuthal direction periodic boundary conditions are used for all variables.
In the vertical direction we apply an equatorial symmetry and reflective
boundary condition for $\theta_{min}$. The radiation boundary conditions are
set to reflective for the $r$ direction (both lower and upper), in
$\theta$-direction we use a fixed value of $E=a_{\mathrm{R}}T^{4}$ with
$T=5\,\mathrm{K}$ at $\theta_{min}$ (which denotes the disc surface), and a
symmetric boundary condition holds at the disc’s midplane $\theta_{max}$. For
the $\phi$-direction we use periodic boundary conditions.

In both cases we used for the matrix solver a relative tolerance of
$\epsilon_{r}=10^{-8}$. In the simulation with hydrodynamics we use the
Riemann-solver hllc222Harten, Lax, Van Leer approximate Riemann Solver with
the contact discontinuity.

(a) $t=10\text{ orbits}$

(b) $t=100\text{ orbits}$

Figure 8: Radial mid-plane temperature profile in the simulations with PLUTO
(red dots) and with the code RH2D (black line) after $t=10\,\mathrm{orbits}$
8(a) and $t=100\,\mathrm{orbits}$ 8(b), together with the absolute value of
the relative error (blue dashed line) which belongs to the log axis on the
right.

####  4.6.2 The static case

In this test case only the radiative equations are solved without the
hydrodynamics. In order to account for the viscous heating in this case, we
add an additional dissipation contribution, $D$, to the right hand side of the
internal energy equation in (17). We consider standard viscous heating, and
include only the main contribution due to the approximately Keplerian shear
flow. At the individual grid points the dissipation is then given by

$D_{i,j,k}=r_{i}^{2}\rho_{i,j,k}\nu\left(\frac{\partial\Omega_{i,j,k}}{\partial
r_{i}}\right)^{2}\,,$ (47)

where $\nu$ is the constant viscosity and $\Omega_{i,j,k}$ the angular
velocity at the individual grid points. In summary we solve the same equations
as in the case with irradiation, when we substitute $S_{i,j,k}$ with
$D_{i,j,k}$.

(a) $\gamma=5/3$

(b) $\gamma=1.1$

Figure 9: Radial mid-plane temperature profile in the simulation with PLUTO
(red line), RH2D (black line) and NIRVANA (blue line) in the quasi-equilibrium
state after 100 orbits in the case with $\gamma=5/3$ without convection 8(a)
and in the strongly convective case with $\gamma=1.1$ 8(b). Additionally we
added the results of a simulation performed with PLUTO where we use a
logarithmic grid in $r$-direction (green line).

In the steady state, the time derivatives in the equations (17) vanish and the
system is reduced to the following equation for the radiation energy density

$\nabla\cdot\,\left(\frac{c\lambda}{\kappa_{\mathrm{R}}\rho}\nabla
E\right)=D.$ (48)

In optically thick regions, $E=a_{\rm R}T^{4}$ and eq. (48) determines the
temperature stratification within the disc.

Figure 10: Vertical slice of the disc temperature at $t=100\,\mathrm{orbits}$
in the dynamical case with $\gamma=1.1$ showing convection cells. Also plotted
in the inset is the enlarged region from $r=0.4\,a_{\text{jup}}$ to
$0.6\,a_{\text{jup}}$ with the velocity field in the $r-\theta$ plane (black
arrows).

The simulation starts at $t=0\,\mathrm{orbits}$ and is evolved until
$t=100\,\mathrm{orbits}$ are reached, where one orbit corresponds to the
Keplerian orbital period at the distance of $a_{\text{jup}}$ which is given
here by $3.732\cdot 10^{8}\,\mathrm{s}$. The initial and overall time step was
chosen as $\Delta t=10^{-3}\,\mathrm{orbits}=3.732\cdot 10^{5}\,\mathrm{s}$.
The results for the static case are shown in figure 8 using here a value of
$\gamma=7/5$ for the adiabatic index. The plots show the radial temperature
profile of the accretion disc in the mid-plane for the simulations after 10
orbits (top panel) and after 100 orbits (bottom panel). We display results of
two different simulations, one done with the code PLUTO (red dots) using the
described methods, and the second (black lines) run with the code RH2D (Kley,
1989). The result from both codes are nearly identical. Even after
$100\,\mathrm{orbits}$ the absolute value of the relative error is always less
than $2\%$. The test shows that the time-scale of the radiative evolution, as
well as the equilibrium state is captured correctly. We note that the code
RH2D uses the one-temperature approach of radiation transport in this case.

####  4.6.3 The dynamical case

The final equilibrium of the described static case does not depend on the
magnitude of $\gamma$, because the viscous heating is independent of it, see
eq. (48). The situation is different, however, for the dynamical cases, where
the hydrodynamical evolution of the flow is taken into account. Since the time
scale of the radiative transport depends on $\gamma$ (through eq. 6), one
might expect the possibility of convective instability, see for example the
recent work by Bitsch et al. (2013a). This is indeed the case for small enough
values of $\gamma$. In order to demonstrate the correctness of our
implementation also for the full dynamical problem, we modelled two discs, one
with $\gamma=5/3$ which clearly shows no convection, and the other with
$\gamma=1.1$ which shows strong convection. The initial setup was identical to
that described before, but now we solve the equations of viscous hydrodynamics
with radiation transport, but without irradiation and explicit dissipation.
Please note that for viscous flows the energy generation due to viscous
dissipation is automatically included in the total energy equation. The
equations (1) to (3) are solved by PLUTO , and the system of equations (17)
are solved as described in section 3. Since this setup is very dynamical and
requires a more complex interplay of hydrodynamics and radiative transport, we
use an additional third code, NIRVANA, for comparison. The NIRVANA code has
been used in Kley et al. (2009) and Bitsch et al. (2013a) on very similar
setups. The results of the two cases are shown in Fig. 9. In the top panel (a)
we display the result for the $\gamma=5/3$ case which is not convective. Here,
the agreement between the codes is excellent with the maximum deviation in the
percentage range. In the lower panel (b) we display the results for the
$\gamma=1.1$ case. Here the radiative transport time-scale is enhanced which
leads to a strongly convective situation, which can be seen in the raggedness
of the curves. In this simulation we doubled the spatial resolution, compared
with the $\gamma=5/3$ case, such that the convection cells are reasonably well
resolved, see figure 10. The agreement between the three different codes is
very good, despite of the very different solution methods for the
hydrodynamics equations: PLUTO uses the total energy equation with a Riemann-
solver while RH2D and NIRVANA use a second-order upwind scheme and the thermal
energy equation. Additionally, the latter two codes use the full dissipation
function and the one-temperature approach.

####  4.6.4 Parallel scaling

Figure 11: Parallel scaling benchmark results for the static accretion disc
test case. We plot here the number of processor cores against
$\frac{t_{64}}{t_{\mathrm{N}}}$ where $t_{N}$ is the runtime used on $N$
processors accordingly for $t_{64}$. The used run-times with full
hydrodynamics and radiation transport for $64$, $128$, $256$, $512$ and $1024$
cpu cores (red crosses) are shown together with the ideal case (black dashed
line).

In order to test the parallel scaling of our new implementation, we used the
same setup as in section 4.6.3 and increased the number of grid cells to
$1024\times 64\times 256$. The computations were only run until
$t=5\,\mathrm{orbits}$, and we used the solver PETSc. So we were able to run
the test on $64$ up to $1024$ processor cores within a reasonable time. The
simulations were run on clusters of the BWGrid which are equipped with Intel
Xeon E5440 cpus and have a low latency InfiniBand network. In figure 11 we
show the results of the simulations performed with full hydrodynamics and
radiation transport. The run-time increases nearly by a factor of two when
doubling the number of cores. With this setup, solving the hydrodynamics
equations needs between $40\%$ and $50\%$ of the computation time and the
radiation transport the remaining $60\%$ to $50\%$, however, these numbers are
strongly problem-dependent. Therefore even up to 1024 cores, we see good
agreement with ideal scaling. According to Amdahl’s law the full code,
including the original part of PLUTO and our implementation of the radiation
transport, is well parallelised.

##  5 Summary and conclusions

We described the implementation of a new radiation module to the PLUTO code.
The module solves for the flux-limited diffusion approximation in the two-
temperature approach. For discretisation the finite volume method is used, and
the resulting difference equations couple the updates of the temperature and
radiation energy density. Due to possibly severe time step limitations, the
set of equations is solved implicitly. For treating the non-linearity of the
temperature in the matter-radiation coupling term, we utilize the method of
Commerçon et al. (2011).

The accuracy of the implementation has been verified using different physical
and numerical setups. The first set of tests deals with purely radiative
problems that include the purely diffusive evolution towards an equilibrium,
and special setups to test the coupling terms between radiative and thermal
energy. A newly developed setup checks for the correct inclusion of the
irradiation from a central source in a spherical coordinate system.

In the second test suite we study the full simultaneous evolution of
hydrodynamics and radiation. First, sub- and super-critical radiative shock
simulations are performed and their outcomes agree very well with published
results of identical setups. Finally, we study the onset of convection in
internally heated viscous discs, and find very good agreement between 3
different, independent hydrodynamical codes. This last test also allowed us to
test the correct implementation in a spherical coordinate system and a non-
equidistant logarithmic grid. Our numerical performance tests indicate
excellent parallel scaling, up to at least 1024 processors.

The current version of the radiation module comes with routines for the
Rosseland mean opacity from Lin & Papaloizou (1985) and Bell & Lin (1994).
Additionally it is possible to use the Rosseland and Planck mean opacities
from Semenov et al. (2003).

The described radiation module can be easily used within the PLUTO
-environment. It can be found on the webpage 333http://www.tat.physik.uni-
tuebingen.de/~pluto/pluto_radiation/ as a patch for the version 4.0 of PLUTO .

###### Acknowledgements.

We gratefully thank the bwGRiD project444bwGRiD (http://www.bw-grid.de),
member of the German D-Grid initiative, funded by the Ministry for Education
and Research (Bundesministerium für Bildung und Forschung) and the Ministry
for Science, Research and Arts Baden-Württemberg (Ministerium für
Wissenschaft, Forschung und Kunst Baden-Württemberg). for the computational
resources. We gratefully acknowledges support through the German Research
Foundation (DFG) through grant KL 650/11 within the Collaborative Research
Group FOR 759: The formation of Planets: The Critical First Growth Phase. We
thank Rolf Kuiper for many stimulating discussions, either physical or
technical.

## References

  * Aubert & Teyssier (2008) Aubert, D. & Teyssier, R. 2008, MNRAS, 387, 295
  * Balay et al. (2012) Balay, S., Brown, J., , et al. 2012, PETSc Users Manual, Tech. Rep. ANL-95/11 - Revision 3.3, Argonne National Laboratory
  * Bell & Lin (1994) Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987
  * Bitsch et al. (2013a) Bitsch, B., Boley, A., & Kley, W. 2013a, A&A, 550, A52
  * Bitsch et al. (2013b) Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013b, A&A, 549, A124
  * Commerçon et al. (2011) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35
  * Davis et al. (2012) Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9
  * Eggum et al. (1988) Eggum, G. E., Coroniti, F. V., & Katz, J. I. 1988, ApJ, 330, 142
  * Ensman (1994) Ensman, L. 1994, ApJ, 424, 275
  * Flaig (2011) Flaig, M. 2011, PhD thesis, Universität Tübingen
  * Freytag et al. (2012) Freytag, B., Steffen, M., Ludwig, H.-G., et al. 2012, Journal of Computational Physics, 231, 919
  * González et al. (2007) González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429
  * Jiang et al. (2012) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14
  * Kley (1989) Kley, W. 1989, A&A, 208, 98
  * Kley et al. (2009) Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971
  * Kuiper et al. (2012) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2012, A&A, 537, A122
  * Kuiper et al. (2010) Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81
  * Kuiper & Klessen (2013) Kuiper, R. & Klessen, R. S. 2013, A&A, 555, A7
  * Levermore & Pomraning (1981) Levermore, C. D. & Pomraning, G. C. 1981, ApJ, 248, 321
  * Lin & Papaloizou (1985) Lin, D. N. C. & Papaloizou, J. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 981–1072
  * Masset et al. (2006) Masset, F. S., D’Angelo, G., & Kley, W. 2006, ApJ, 652, 730
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mihalas & Mihalas (1984) Mihalas, D. & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics
  * Minerbo (1978) Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541
  * Pomraning (1973) Pomraning, G. C. 1973, The equations of radiation hydrodynamics
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
  * Stone et al. (1992) Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819
  * Turner & Stone (2001) Turner, N. J. & Stone, J. M. 2001, ApJS, 135, 95