# Radiation transport methods in star formation simulations

Richard Wünsch 1,∗

###### Abstract

## 1

Radiation transport plays a crucial role in star formation models, as certain
questions within this field cannot be accurately addressed without taking it
into account. Given the high complexity of the interstellar medium from which
stars form, numerical simulations are frequently employed to model the star
formation process. This study reviews recent methods for incorporating
radiation transport into star formation simulations, discussing them in terms
of the used algorithms, treatment of radiation frequency dependence, the
interaction of radiation with the gas, and the parallelization of methods for
deployment on supercomputers. Broadly, the algorithms fall into two
categories: (i) moment-based methods, encompassing the flux-limited diffusion
approximation, M1 closure, and variable Eddington tensor methods, and (ii)
methods directly solving the radiation transport equation, including forward
and reverse ray tracing, characteristics-based methods, and Monte Carlo
techniques. Beyond discussing advantages and disadvantages of these methods,
the review also lists recent radiation hydrodynamic codes implemented the
described methods.

##  2 Keywords:

star formation, interstellar medium, numerical techniques, radiation
transport, astrophysical fluid dynamics

##  3 Introduction

Electromagnetic radiation plays a fundamental role in the process of star
formation as it regulates the dynamics, thermodynamics, and chemistry of star-
forming regions. The interstellar radiation field (ISFR) heats the
interstellar dust, dissociates molecules and affects the ionization balance,
which is critical for the chemistry of the interstellar medium and the
formation of molecular clouds. During the collapse of molecular cores, the
dust absorbs the cooling radiation and once the medium becomes optically
thick, the first proto-stellar core is formed. When the first stars are
formed, their radiation exerts a profound influence on the surrounding gas and
dust shaping the morphology of molecular clouds. This set of processes is
called the radiative feedback as it regulates the rate and efficiency of star
formation. Extreme ultraviolet photons from young, massive stars can ionize
the surrounding gas, creating HII regions. The ionizing radiation also
influences the chemistry of the gas, affecting the abundance of molecules
necessary for the cooling and fragmentation of the cloud. Furthermore,
radiation pressure from newly formed stars can counteract gravitational
collapse, slow down the gas inflow and drive the winds and outflows.

The inclusion of radiation transport in star formation simulations is
therefore desirable. However, it is also very hard due to both the high
computational costs and the complexity of the processes that have to be taken
into account. The radiation transport equation (RTE) in its general form reads

$\frac{1}{c}\frac{\partial I_{\nu}}{\partial t}+\mathbf{n}\cdot\nabla
I_{\nu}=-(\kappa_{\nu,a}+\kappa_{\nu,s})\rho
I_{\nu}+j_{\nu}\rho+\frac{1}{4\pi}\kappa_{\nu,s}\rho\int_{\Omega}I_{\nu}d\Omega$
(1)

where the radiation field is represented by the specific intensity $I_{\nu}$.
It is defined through the element of the radiation energy $dE_{\nu}$ flowing
across an area element $da$ located at position $\mathbf{x}$ in time $dt$ in
the solid angle $d\Omega$ about the direction $\mathbf{n}$ in the frequency
interval $\nu+d\nu$ as

$dE_{\nu}=I_{\nu}(\mathbf{x},\mathbf{n},t)\cos(\theta)\,d\nu\,da\,d\Omega\,dt$
(2)

where $\theta$ is the angle that $\mathbf{n}$ makes with a normal to the area
element $da$. The first term on the r.h.s. of Equation 1 represents the
radiation losses due to absorption (described by the absorption opacity
$\kappa_{\nu,a}$) and scattering (described by the scattering opacity
$\kappa_{\nu,}$); $\rho$ is the mass density of the medium through which the
radiation propagates. The second term on the r.h.s. of Equation 1 represents
the emission of radiation describe by the emission coefficient $j_{\nu}$, the
last r.h.s term represents contribution of the radiation scattered from other
directions. Note that in the above equations, the specific intensity is
expressed in the frame of coordinate system, however, emission coefficient and
opacities are defined in the comoving frame.

The structure of Equations 1 and 2 reveals that the computational costs
associated with solving them directly in their complete generality are
exceedingly high. Specifically, in the case of three-dimensional problems,
$I_{\nu}$ becomes a function of three components of $\mathbf{x}$, two
components of $\mathbf{n}$, time $t$, and the frequency of the radiation $\nu$
(written as a subscript, as it is often grouped into several ranges known as
wavebands) – amounting to seven independent variables. This issue becomes
particularly pronounced when implementing the radiation transport algorithm
within a star formation simulation that demands high resolution and,
consequently, a substantial number of computational points along the
aforementioned variables. To address this, various computational techniques
have been developed that solve only a portion of Equation 1 (e.g.,
incorporating only a subset of the processes it describes) or utilize diverse
approximations, including those transforming it into a different (albeit
related) equation that can be solved more computationally efficiently.

This paper reviews the most common techniques to solve the radiation transport
equation within (magneto-)hydrodynamic simulations of star formation. The
focus of this work is primarily on continuum radiation transport, rather than
line transport (atomic or molecular). There is an extensive amount of work on
solving radiation transport on static grids or as tools for synthetic
observations, however, covering it is beyond the scope of this work.
Additional information can be found e.g. in classical works on radiation
hydrodynamics Mihalas and Mihalas (1984); Castor (2004), in a review of dust
radiative transfer by Steinacker et al. (2013) and in a review of numerical
methods in simulations of star formation by Teyssier and Commerçon (2019). The
paper is organised as follows. Section §4 describes different way how
radiation transport codes can be classified. Section §5 reviews the most
common radiation transport algorithms used in star formation, consisting of
§5.1 discussing the moment-based methods and §5.2 describing the methods
solving Equation 1 directly. Finally, Section §6 summarizes this review.

##  4 Classifications of the radiation transport treatments

There are several criteria based on which radiation transport codes can be
classified.

Firstly, the method can either explicitly include the first term of Equation 1
(partial time derivative) or assume the infinite speed of light and solve the
RTE equation implicitly. The former approach, also known as the evolutionary
method, is typically employed by moment-based methods. In some cases, the
speed of light is artificially reduced to avoid excessively small time steps
mandated by the stability conditions of an explicit solver – a technique
demonstrated to be useful in other fields as well (e.g. Cohen et al., 1999).
The implicit approach, alternatively termed the instantaneous method, is often
utilized by codes based on the method of characteristics, ray tracing, and
Monte Carlo methods.

Secondly and related to the above, the method can either directly solve
Equation 1 (albeit with some approximations) or transform it into another type
of equation. The former approach includes ray tracing, long and short
characteristic methods, and Monte Carlo methods. While ray tracing and long
characteristic methods are among the most accurate, they are also
computationally expensive. The computational cost of ray tracing typically
scales with the number of discrete radiation sources, making it particularly
expensive for simulations with a high number of sources. In the latter
approach, Equation 1 is commonly transformed into a hyperbolic system (M1
closure methods) or a parabolic system (flux-limited diffusion method).

Thirdly, various methods concentrate on distinct wavelengths of radiation,
incorporating different physical processes related to the interaction between
radiation and matter. In the simplest approach, termed monochromatic, all
photons are assumed to possess a single (representative) frequency. A slightly
more complex method is the grey approximation, where a specific constant
spectrum – often the black body spectrum – is assumed for all radiation. In
certain scenarios, it becomes necessary to categorize radiation based on its
frequency into multiple groups that exhibit qualitative differences (e.g.,
where only some frequencies ionize a specific species). This approach is
referred to as multi-group frequency treatment. Lastly, certain codes have the
capability to encompass full frequency dependence by sampling the frequency
with a designated number of bins.

An essential consideration in implicit radiation transport methods is whether
the opacity at a specific location is depends on the radiation intensity at
that point. A notable instance of this behavior is the treatment of ionizing
(EUV) radiation using the on-the-spot approximation. In this approach,
ionizing photons are effectively destroyed by case B recombinations, as they
do not result in the replication of ionizing photons. The number of case B
recombinations occurring in a given gas parcel needs to be distributed among
all the rays passing through it, particularly when using ray tracing or the
long characteristic method. Consequently, such methods necessitate iteration,
leading to an increase in computational cost.

Finally, different methods pose varying challenges in terms of
parallelization, with particular emphasis on the parallelization using domain
decomposition suitable for distributed memory machines. Many hydrodynamic
codes opt for this type of parallelization, as distributed memory machines
offer extensive memory and processor core capabilities. Unfortunately, this
parallelization scheme is also the most complex from a coding perspective.
Radiation transport, being primarily a long-distance interaction problem,
presents significant difficulties in parallelization. This is especially true
for ray tracing, the long characteristics method, and Monte Carlo, where
numerous calculations over considerable distances must be executed within a
limited timeframe.

##  5 Radiation transport algorithms

###  5.1 Moment-based methods

Instead of the specific intensity $I_{\nu}$, the radiation field can be
described by its moments obtained by integrating $I_{\nu}$ over all
directions. The first three moments are defined as (see e.g. Shu, 1991):

$\left(\begin{array}[]{c}cE_{\nu}\\\ \mathbf{F}_{\nu}\\\
c\mathbb{P}_{\nu}\end{array}\right)=\oint\left(\begin{array}[]{c}1\\\
\mathbf{n}\\\ \mathbf{n\otimes n}\end{array}\right)I_{\nu}d\Omega$ (3)

where $\mathbf{n}$ is a unit vector. The mean radiation energy density
$E_{\nu}$, the radiation flux $\mathbf{F}_{\nu}$, and the radiation pressure
tensor $\mathbb{P}_{\nu}$ are the zeroth, first and second moments of the
radiation field, respectively. The radiation pressure is often expressed as
$\mathbb{P}_{\nu}=\mathbb{D}E_{\nu}$ where $\mathbb{D}$ is the so called
Eddington tensor.

Moments of Equation 1 can be obtained by multiplying it with a certain power
of $\mathbf{n}$ and integrating it over all directions. The scattering is
assumed to be isotropic. Then, the first two moment equations are

$\frac{\partial E_{\nu}}{\partial t}+\nabla\cdot\mathbf{F}_{\nu}=4\pi\rho
j_{\nu}-\rho\kappa_{\nu,a}cE_{\nu}$ (4)

and

$\frac{1}{c}\frac{\partial\mathbf{F}_{\nu}}{\partial
t}+c\nabla\cdot\mathbb{P}_{\nu}=-\rho\kappa_{\nu,a}\mathbf{F}_{\nu}\ .$ (5)

In the above equations, the opacity and the emission coefficient are expressed
in the frame of reference co-moving with the fluid (where they are typically
isotropic), however, $E_{\nu}$, $\mathbf{F}_{\nu}$ and $\mathbb{P}_{\nu}$ are
expressed in the coordinate system of the simulation (the laboratory frame).
This problem is usually solved by either transforming Equations 4 and 5 into
the co-moving frame (CMF) or by transforming $\kappa_{\nu,a}$ and $j_{\nu}$
into the laboratory frame using the first order expansion (mixed frame
formulation). The radiation moment equations in the co-moving frame have a
form (see Mihalas and Mihalas, 1984)

$\frac{\partial E_{\nu}}{\partial
t}+\nabla\cdot(\mathbf{u}E_{\nu})+\nabla\cdot\mathbf{F}_{\nu}+\mathbb{P}:\mathbf{u}=4\pi\rho
j_{\nu}-\rho\kappa_{\nu,a}cE_{\nu}$ (6)

and

$\frac{1}{c}\frac{\partial\mathbf{F}_{\nu}}{\partial
t}+\frac{1}{c}\nabla\cdot(\mathbf{uF}_{\nu})+c\nabla\cdot\mathbb{P}_{\nu}=-\rho\kappa_{\nu,a}\mathbf{F}_{\nu}$
(7)

where $\mathbf{u}$ is the velocity of the gas. For details on the mixed frame
formulation see e.g. Mihalas and Klein (1982). The advantage of the mixed
frame approach is that the radiation transfer equations are significantly
simpler and that it conserves the total energy, however, it is not suitable
for computing for instance line opacities in the supersonic flow due to their
rapid changes (Moens et al., 2022).

####  5.1.1 Flux limited diffusion (FLD)

The simplest moment-based method was derived as an extension of the radiation
transport in a diffusion limit (Minerbo, 1978; Levermore and Pomraning, 1981).
If the photon free mean path is small and the optical depth is high, the
radiation is isotropic and the Eddington tensor has form
$\mathbb{D}=\frac{1}{3}\mathbb{I}$ where $\mathbb{I}$ is the identity tensor.
Assuming the radiation field is stationary, i.e. neglecting the first two
terms in Equation 7, and considering
$\nabla\cdot\mathbb{P}_{\nu}=\frac{1}{3}\nabla E_{\nu}$, the radiation flux is
related to the radiation energy density as

$\mathbf{F}_{\nu}=-\frac{c}{3\kappa_{\nu}\rho}\nabla E_{\nu}$ (8)

and the radiation transport problem reduces to the Fick’s diffusion law.

Relation 8 is not valid if the optical depth is low, because
$\mathbf{F}_{\nu}$ can exceed the physical limit
$\mathbf{F}_{\mathrm{max}}=cE_{\nu}$ there. Therefore, the FLD method
introduces the so called flux limiter $\lambda$ and modifies Equation 8 in the
following way

$\mathbf{F}_{\nu}=-\frac{\lambda(R)c}{\kappa_{\nu}\rho}\nabla E_{\nu}\ .$ (9)

The flux limiter $\lambda(R)$ is defined as a function of the normalized
radiation energy density gradient

$R=\frac{|\nabla E_{\nu}|}{\kappa_{\nu}\rho E_{\nu}}\ ,$ (10)

in such a way that it ensures Equation 9 reduces to the Fick’s law in the
optically thick regime (i.e. $\lambda=1/3$) and to the maximum allowed flux,
$\mathbf{F}_{\mathrm{max}}$, in the optically thin regime (i.e.
$\lambda=1/R$). Various flux limiters fulfilling the above requirements have
been suggested in the literature (see e.g. Minerbo, 1978; Levermore and
Pomraning, 1981; Levermore, 1984; Turner and Stone, 2001). As an example we
give the flux limiter of Levermore and Pomraning (1981) in the form

$\lambda(R)=\frac{1}{R}\left(\mathrm{coth}(R)-\frac{1}{R}\right)\simeq\frac{2+r}{6+3R+R^{2}}$
(11)

and the corresponding approximate Eddington tensor

$\mathbb{D}=\frac{1}{2}(1-f)\mathbb{I}+\frac{1}{2}(3f-1)(\mathbf{n}\otimes\mathbf{n})\qquad\mathrm{where}\qquad
f=\lambda+\lambda^{2}R^{2}\ .$ (12)

The FLD method then solves only the $0$-th moment equation in the form

$\frac{\partial E_{\nu}}{\partial
t}+\nabla\cdot(\mathbf{u}E_{\nu})+\nabla\cdot\frac{c\lambda(R)}{\kappa_{\nu}\rho}\nabla
E_{\nu}+\mathbb{P}:\mathbf{u}=4\pi\rho j_{\nu}-\rho\kappa_{\nu,a}cE_{\nu}$
(13)

which is, from a mathematical point of view, a parabolic equation, and hence
the standard parabolic solvers can be used.

The above equations are monochromatic, i.e. they are written for a single
wavelength of the radiation. However, it is often practical to use the FLD
method (and other methods as well) for radiation with a certain range of
wavelengths. Then we speak about the grey radiation transport if there is a
single range covering the broad wavelength range important for a given
problem, or about the multi-group / multi-band radiation transport if there
are several wavelength ranges. Additionally, the grey methods often use the
Planck black body law, $B_{\nu}(T)$, as the source function, and then the
$0$-th moment equation has form

$\frac{\partial E_{r}}{\partial
t}+\nabla\cdot(\mathbf{u}E_{r})+\nabla\cdot\frac{c\lambda}{\kappa_{\mathrm{R}}}\nabla
E_{r}+\mathbb{P}:\mathbf{u}=\kappa_{\mathrm{P}}(a_{r}T^{4}-cE_{r})$ (14)

where $E_{r}\equiv\int E_{\nu}d\nu$ is the bolometric radiation energy
density, $T$ is the gas/dust temperature, $a_{r}=4\sigma_{\mathrm{B}}/c$ is
the radiation constant, $\sigma_{\mathrm{B}}$ is the Stefan-Boltzmann
constant, and $\kappa_{\mathrm{P}}$ and $\kappa_{\mathrm{R}}$ are Planck and
Rosseland mean opacities defined as

$\kappa_{\mathrm{P}}\equiv\frac{\int\kappa_{\nu}B_{\nu}(T)d\nu}{\int
B_{\nu}(T)d\nu}$ (15)

and

$\frac{1}{\kappa_{\mathrm{R}}}\equiv\frac{\int\frac{1}{\kappa_{\nu}}\frac{\partial
B_{\nu}(T)}{\partial T}d\nu}{\int\frac{\partial B_{\nu}(T)}{\partial T}d\nu}\
.$ (16)

RHD codes using the FLD method are discussed below and summarized in Table 1
giving a reference to the code paper, code name, whether it uses co-moving or
mixed frame of reference, whether it is parallel with the domain decomposition
and some additional information. The first usage of the FLD method is probably
by Alme and Wilson (1973) for the problem of the X-ray emission resulting from
the accretion of material onto a neutron star. In the field related to star
formation, the grey FLD method was used by Kley (1989); Bodenheimer et al.
(1990); Klahr et al. (1999); Boss (2001) for simulations of accretion discs
around young low mass stars. It yielded reasonably correct temperatures in the
optically thick and optically thing regions, however, it showed deviations
from the correct solution in the transition regions (Boley et al., 2007). The
treatment of radiation in such models has been significantly improved by
Kuiper et al. (2010) who developed a hybrid scheme combining the grey FLD
method with the frequency dependent ray tracing, and implemented it to the MHD
code Pluto. They at first calculate the stellar radiation up to the first
absorption by a simple ray tracing algorithm and subsequently, they use FLD to
calculate the radiation re-emitted by the dust. This approach turned out to be
very successful and many other authors implemented this or similar schemes
into their codes. Bitsch et al. (2013) used it with the Nirvana code to study
the influence of opacity and stellar irradiation on the disc structure and on
the migration of planets. Flock et al. (2013) and Kolb et al. (2013) created
two independent implementations of this hybrid method for the Pluto code.
Klassen et al. (2014) generalized this method for multiple sources and
arbitrary geometry wrote it as a module for the MHD code Flash. Similarly,
Ramsey and Dullemond (2015) implemented a generalized version of this method
for the AZEuS code – a Zeus family AMR code with a fully staggered mesh.

The simplicity of the FLD method allowed it to be implemented into many
general purpose (magneto-)hydrodynamic codes used in astrophysics. The widely
used Zeus2d code was originally released with the radiation module based on
the tensor variable Eddington factor, however, the FLD method was written for
it by Turner and Stone (2001). The Orion code used for many simulations of
star formation, including the feedback from massive stars, was developed with
FLD method being part of it (Krumholz et al., 2007). Gittings et al. (2008)
developed the AMR code Rage including the grey FLD, for which they introduce a
new technique of solving the diffusion and material energy equations, allowing
larger time steps and more robust behavior. Commerçon et al. (2011)
implemented FLD into the AMR code Ramses, and Commerçon et al. (2014) added
adaptive time stepping and the implicit time integration. The block AMR code
Crash (van der Holst et al., 2011) was developed with grey FLD where electrons
and ions are allowed to have different temperatures. The grey FLD method in
its mixed-frame formulation is also part of the AMR code Castro (Zhang et al.,
2011), one of the first 3D RHD codes using the unsplit PPM solver. Trhd (Sijoy
and Chaturvedi, 2015) is two-dimensional unstructured mesh Lagrangian code
using the grey FLD with three temperature approach (ions, electrons and
radiation are allowed to have different temperatures). The Flash code (Fryxell
et al., 2000), a computational framework consisting of modules implementing
various physical processes including hydrodynamics and radiation transport,
was provided by the FLD module implemented by Chatzopoulos and Weide (2019)
and tested on 1D simulations of supernova Ia explosions. The general fluid
dynamics code Gizmo (Hopkins and Grudić, 2019) includes several radiation
transport solvers including the one based on FLD. Recently, the grey FLD
method was implemented into the finite-volume MHD code Mpi-Amrvac (Moens et
al., 2022).

The FLD method was implemented also into the Smoothed Particle Hydrodynamics
(SPH) codes. Basic methods for calculating radiative heat diffusion were
described in seminal works on SPH by Lucy (1977) and Brookshaw (1985). The
full grey FLD method was implemented into the SPH by Whitehouse and Bate
(2004) and Whitehouse et al. (2005). Later, this code originating from code
SphNG (Benz, 1990) was used by Price and Bate (2009); Bate and Keto (2015);
Bate (2019) and many other works for simulating star forming clouds, where the
FLD method can be efficiently used to follow collapsing pre-stellar cores with
high optical depths. Mayer et al. (2007) used implemented FLD into SPH code
Gasoline and used it to study the gravitational instability in protoplanetary
discs. Fryer et al. (2006) developed code Snsph including the FLD radiation
transport and used it for simulating supernova explosions, where FLD can also
be used efficiently due to high optical depths. Another implementations of FLD
into SPH codes are by Viau et al. (2006) and recently by Bassett et al. (2021)
(code Spheral).

Reference | Name/code | frame | DD | Note
---|---|---|---|---
Kley (1989) | - | CMF | N | 2D axisym, protoplanetary discs
Bodenheimer et al. (1990) | - | CMF | N |  2D HD code by Rozyczka (1985)
Klahr et al. (1999) | Tramp | CMF | N |
Boss (2001) | - | CMF | N | 3D disc forming giant planet
Turner and Stone (2001) | Zeus2d | CMF | N |
Fryer et al. (2006) | Snsph | CMF | Y | SPH code, neutrino diffusion
Viau et al. (2006) | - | CMF | N | SPH code
Krumholz et al. (2007) | Orion | mixed | Y |
Mayer et al. (2007) | Gasoline | CMF | Y |
Gittings et al. (2008) | Rage | CMF | Y |  rad-matter coupling $\rightarrow$ larger dt
Kuiper et al. (2010) | Pluto | negl. | Y |  \+ $\nu$-dep ray tr.
Commerçon et al. (2011) | Ramses | CMF | Y |
van der Holst et al. (2011) | Crash | CMF | Y | multiple frequency groups
Zhang et al. (2011) | Castro | mixed | Y | unsplit PPM solver
Bitsch et al. (2013) | Nirvana | negl. | N |  \+ $\nu$-dep ray tr. in rad. dir.
Flock et al. (2013) | Pluto | negl. | Y |  \+ $\nu$-dep ray tr.
Kolb et al. (2013) | Pluto | negl. | Y |  \+ $\nu$-dep ray tr.
Klassen et al. (2014) | Flash | CMF | Y |  \+ $\nu$-dep ray tr., multiple src.
Bate and Keto (2015) | sphNG | CMF | Y |  \+ TreeCol (Clark et al., 2012)
Sijoy and Chaturvedi (2015) | Trhd | CMF | N | 3-temp, 2D unstructured-mesh
Ramsey and Dullemond (2015) | AZEuS | CMF | N |  \+ $\nu$-dep ray tr.
Chatzopoulos and Weide (2019) | Flash | mixed | Y | opacities from OPAL
Hopkins and Grudić (2019) | Gizmo | mixed/CMF | Y | multiple frequency groups
Bassett et al. (2021) | Spheral | CMF | Y | SPH code
Moens et al. (2022) | Mpi-Amrvac | CMF | Y | finite-volume MHD code
Table 1: Radiation transport codes using the flux limited diffusion (FLD)
approximation. The first column gives the reference to the work describing the
code or algorithm; the second column give the name of the radiation transport
code or RHD code (if they exist); the third column specifies the frame of
reference in which the radiation energy equation is formulated (mixed, co-
moving – CMF, or radiation energy advection terms are neglected); the fourth
column tells whether the code is parallel using the domain decomposition; and
the last column gives eventually an additional information.

####  5.1.2 M1 closure

A more complex and accurate two-moments method is the M1 closure (Dubroca and
Feugeas, 1999). It assumes that the Eddington tensor is axially symmetric
around the flux vector and hence can be expressed as

$\mathbb{D}=\frac{1-\chi}{2}\mathbb{I}+\frac{3\chi-1}{2}\mathbf{n}\otimes\mathbf{n}$
(17)

where $\chi$ is a scalar function called the Eddington factor. The M1 closure
method uses one of its simplest forms discussed in Levermore (1984),
corresponding to the Lorentz boosted isotropic distribution of the radiation,

$\chi=\frac{3+4|\mathbf{f}|^{2}}{5+2\sqrt{4-3|\mathbf{f}|^{2}}}$ (18)

where $\mathbf{f}=|\mathbf{F}_{\nu}|/cE_{\nu}$ is the reduced flux. Then, the
method solves Equations (6) and (7) or their equivalent in the mixed frame
approach.

The advantage of the M1 closure is that it captures correctly the two extremes
of the radiation field: (i) the free streaming (or optically thin) limit with
$|\mathbf{f}|=1$ and $\chi=1$, and (ii) the diffusion (or optically thick)
limit with $|\mathbf{f}|=0$ and $\chi=1/3$. The ability to maintain the
original direction of the free streaming radiation is a pronounced improvement
in comparison to the FLD method leading to the artificial diffusion of the
free streaming radiation. On the other hand, the assumed simple form of the
radiation pressure tensor cannot correctly describe more complex radiation
fields as for instance the radiation from multiple sources in which case the
nonphysical interaction would occur.

A specific feature of the M1 closure making is suitable for the implementation
into (magneto-)hydrodynamic codes is that the system of Equations 6 and 7 is
hyperbolic and hence it can be solved with the same algorithms as the
equations of the fluid dynamics. Table 2 gives a list of star formation and
related fields codes implementing the M1 closure method. The first such code
was Heracles (González et al., 2007) designed to study radiative shocks (and
compare the results with laboratory experiments), jets of young stars,
formation of pre-stellar cores and protoplanetary disks. This method was
implemented into the cosmic reionisation code Aton (Aubert and Teyssier,
2008), and later ported to the general RHD code Ramses (Rosdahl et al., 2013;
Rosdahl and Teyssier, 2015). Sądowski et al. (2013) included the M1 closure
method into the general relativistic code Koral formulating the RHD equations
in the covariant form. Skinner and Ostriker (2013) describe the M1 closure
module for the Athena code, in which they use the reduced speed of light
approximation to improve the performance. The general fluid dynamics code
Gizmo includes the M1 closure as one of its radiation transport solvers.
Kannan et al. (2019) describe the radiation transport module for the
unstructured moving mesh code Arepo, based on a new higher order
implementation of the M1 closure method. This method is also included into
code Fornax (Skinner et al., 2019) primarily intended for core collapse
supernova simulations, general fuild dynamics code Pluto (Melon Fuksman and
Mignone, 2019), and hybrid OpenMP/MPI/GPU parallel code Ark-rt (Bloch et al.,
2021). Chan et al. (2021) present module Sph-m1rt for task-based parallel SPH
code Swift.

Reference | Name/code | frame | DD | Note
---|---|---|---|---
González et al. (2007) | Heracles | CMF | Y |
Rosdahl et al. (2013) | Ramses-rt | CMF | Y |  from ATON (Aubert and Teyssier, 2008)
Sądowski et al. (2013) | Koral | covar. | N | general relativistic MHD code
Skinner and Ostriker (2013) | Athena | mixed | Y | reduced speed of light approx.
Hopkins and Grudić (2019) | Gizmo | CMF/mixed | Y | multiple frequency groups
Kannan et al. (2019) | Arepo | CMF | Y | higher order scheme of M1
Melon Fuksman and Mignone (2019) | Pluto | mixed | Y |
Skinner et al. (2019) | Fornax | CMF | Y | photon and neutrino rad. fields
Bloch et al. (2021) | Ark-rt | CMF | Y | hybrid OpenMP/MPI/GPU parallel
Chan et al. (2021) | Sph-m1rt | CMF | Y | in SPH code SWIFT
Table 2: Radiation transport codes using the M1 closure method. The meaning of
columns is the same as in Table. 1.

####  5.1.3 Variable Eddington Tensor methods

Another approach to solve moment Equations 4 and 5, as opposed to the closure
assuming an approximate form of the Eddington tensor, is to calculate the
Eddington tensor self consistently by directly solving the radiation transport
Equation 1. Methods using it are referred to as Variable Eddington Tensor
(VET) methods, and the Eddington tensor, $\mathbb{D}$, is calculated by one of
the methods described in §5.2 (ray tracing, long or short characteristics, or
Monte Carlo). Subsequently, $\mathbb{D}$ is combined with the radiation moment
equations to calculate the radiation energies and fluxes.

The VET method was developed for 1D calculations of plane-parallel stellar
atmospheres by Auer and Mihalas (1970), however, the idea of iterating the
ratio of second and zeroth radiation moments dates back to Eddington (1926).
The first multi-dimensional astrophysical code using this method was Zeus2d
(Stone et al., 1992), where $\mathbb{D}$ was calculated by the short
characteristics. The algorithm was later updated in code Zeus-mp (Hayes and
Norman, 2003) by adding more methods to calculate $\mathbb{D}$ and by
separating the radiation and gas energy densities. The 1D RHD code Titan
(Gehmeyr and Mihalas, 1994) uses VET and calculates $\mathbb{D}$ by ray
tracing. Gnedin and Abel (2001) developed an algorithm based on the Optically
Thin Variable Eddington Tensor (OTVET) approximation where $\mathbb{D}$ is
calculated in the optically thin regime. This is done by integrating
contributions of radiation sources with attenuation factor $1/r^{2}$ over the
computational domain taking advantage of the solver used for the self-gravity.
They implemented this algorithm into the Soften Lagrangian Hydrodynamic code
(SLH) following all physical quantities on a moving deformed mesh. This
algorithm was later implemented by Petkova and Springel (2009) into the SPH
code Gadget. Sekora and Stone (2010) developed a hybrid Godunov method based
on VET with user defined $\mathbb{D}$ and implemented it into 1D RHD code
Nike. This method was later implemented into code Athena by Jiang et al.
(2012) where $\mathbb{D}$ is calculated by the short characteristics method
described by Davis et al. (2012). The general relativistic radiation
magnetohydrodynamic (GR-RMHD) code Inazuma (Asahina et al., 2020) also uses
the VET based method to treat the radiation, and it calculates $\mathbb{D}$
integrating the local radiation intensity over the angular directions.
Recently, the RHD code Quokka (Wibking and Krumholz, 2022) implements the
general two-moment radiation transport method. The code uses the block
structured adaptive mesh refinement handled by the AmREX library and it is
parallelized for both MPI and usage on graphic processing units (GPUs) using
the CUDA library. Its radiation solver is written to use a general Eddington
tensor, the default option is the M1 closure.

Reference | Name/code | frame | DD | Note
---|---|---|---|---
Stone et al. (1992) | Zeus2d | CMF | N |  $\mathbb{D}$ cmp. by short chars.
Gehmeyr and Mihalas (1994) | Titan | CMF | N |  1D code; $\mathbb{D}$ cmp. by ray tracing
Gnedin and Abel (2001) | Otvet | Lag | N |  $\mathbb{D}$ cmp. in optically thin approx
Sekora and Stone (2010) | Nike | Lag | N |  hybrid Godunov RHD, $\mathbb{D}$ user defined
Jiang et al. (2012) | Athena | mixed | Y |  $\mathbb{D}$ cmp. by short chars.
Asahina et al. (2020) | Inazuma | CMF | N |  $\mathbb{D}$ cmp. by integrating intensity
Wibking and Krumholz (2022) | Quokka | mixed | Y |  GPU, $\mathbb{D}$ general, M1 by default
Table 3: Radiation transport codes using the Variable Eddington Tensor (VET)
approach. The meaning of columns is the same as in Table. 1

###  5.2 Methods solving RTE directly

Table 4 lists codes that employ the direct solution of the RTE: forward and
reverse ray tracing, methods of characteristics, and Monte Carlo methods.

####  5.2.1 Ray tracing

The most straightforward method to solve the radiation transport Equation 1 is
to follow the radiation from its sources along the lines of its propagation, a
technique called ray tracing. It is widely used in many other fields, in
particular in the computer image generation. The original idea of ray tracing
dates back to 16th century to Albrecht Dürer who described multiple techniques
for projecting 3D scenes onto an image plane (Hofmann, 1990).

As previously mentioned, the primary challenge in directly solving Equation 1
is its high dimensionality, resulting in high computational costs. Therefore,
a key characteristic of a ray tracing algorithm is a set of approximations
designed to reduce the complexity of the problem. A crucial aspect is the
method by which the radiation field is discretized in both space and its
propagation directions. The simplest approach, employed by many early and
recent codes, involves calculating the radiation on the same grid as the gas
dynamics. This confines the direction of radiation propagation to the
coordinate axes. Examples of this approach include Wolfire and Cassinelli
(1986) in a 1D spherical grid, Murray et al. (1994); Garcia-Segura and Franco
(1996) in a 2D spherical grid, and numerous other works, the review of which
exceeds the scope of this study. This method is particularly well-suited for
treating radiation emanating from a single source (e.g., a central star)
located at the center of the coordinate system. It can be effectively combined
with other methods that describe the diffuse radiation component (e.g., FLD,
as mentioned in Section 5.1.1 and Table 1).

In many cases, confining the radiation propagation solely along the hydrocode
grid axes is overly restrictive. Therefore, general radiation transport codes
typically employ an independent method to define rays, making ray tracing
algorithms particularly suitable for implementation into unstructured grid or
smooth particle hydrodynamic codes. One approach is to cast rays from a source
in spherical coordinates, associating each ray with a spherical segment of the
same size (Abel et al., 1999; Razoumov and Scott, 1999). This method was later
refined in the HEAlPix library (Górski et al., 2005), which has become a de
facto standard. HEALPix facilitates the tessellation of the sphere surface
into equal area regions in a hierarchical manner, a feature utilized by
adaptive ray tracing algorithms that split rays in regions requiring finer
angular resolution (Abel and Wandelt, 2002). This approach has proven highly
successful and has been employed by various authors, including Pawlik and
Schaye (2008) (Traphic code linked to SPH code Gadget2), Bisbas et al. (2009)
(SPH code Seren), Wise and Abel (2011) (Moray code linked to AMR code Enzo),
Baczynski et al. (2015) (Fervent code within the AMR code Flash), and others.

Depending on the number of cast rays, ray tracing algorithms can achieve high
accuracy, although with correspondingly higher computational costs that
typically scale with the number of radiation sources. As a result, they are
occasionally employed in conjunction with other methods where ray tracing
specifically handles radiation from a relatively low number of discrete
sources, such as massive stars. An example is the Harm2 code (Rosen et al.,
2017), which calculates ionizing stellar radiation using adaptive ray tracing
and the diffuse component using the FLD method in the ORION code. Similarly,
the Torus-3dpdr code (Bisbas et al., 2015a) utilizes adaptive ray tracing for
FUV radiation and the Monte Carlo method for the ionizing and diffuse
components.

A slightly different approach to spatially adaptive ray tracing was used by
Dale et al. (2007) who integrated the RTE on the line-of-sight between the
source and the target particle using evaluation points defined by SPH
particles close to the line-of-sight selected using method by Kessel-Deynet
and Burkert (2000).

Ritzerveld and Icke (2006) developed algorithm SimpleX to solve Equation 1 on
an unstructured grid calculated from the properties of the medium in which the
radiation propagates, using the Delaunay triangulation (Delaunay, 1934). The
computational costs of this method do not depend on the number of sources,
similarly to moment-based methods or some reverse ray tracing schemes. It is
typically less diffusive that moment based methods and it is more suitable for
parallelisation with domain decomposition than most of the ray tracing
methods. It is particularly suitable for being used with SPH or unstructured
grid codes, Clementel et al. (2014) implemented it into an SPH code. The
SimpleX algorithm was improved and parallelized by Paardekooper et al. (2010)
and implemented as code Sprai integrated into code Arepo by Jaura et al.
(2018).

####  5.2.2 Characteristics-based methods

A specific way of ray tracing are methods of characteristics, which solve the
Radiative Transfer Equation (RTE) along predefined rays (characteristics)
covering the entire computational domain rather than casting rays from
discrete sources. In this manner, characteristics-based methods can describe
the diffuse radiation component. They can be categorized into long, short and
hybrid characteristics based on the end points of the rays. Long
characteristics extend across the entire computational domain, with their
separation comparable to the size of the computational elements and defined in
specific directions. While long characteristics are among the most accurate
methods, they are also the most expensive, depending on the number of rays.
Another disadvantage, applicable to some extent to all ray tracing methods, is
that they are difficult to get parallelized within codes that use domain
decomposition on distributed memory machines. This challenge arises because
rays crossing multiple domains typically represent a substantial amount of
information that needs to be communicated over the (relatively) slow network.
An example of a code utilizing long characteristics is Lampray (Frostholm et
al., 2018), implemented into the Ramses code. Lampray calculates both diffuse
and discrete sources radiation on an adaptive oct-tree mesh, allowing the
definition of several frequency bins for both ionizing and non-ionizing
radiation. It is coupled with the non-equilibrium chemistry code Krome (Grassi
et al., 2014), providing the chemical networks to evolve species abundances.

The short characteristics methods create rays only among the neighbour cells
(or generally computational elements) and interpolate the radiation intensity
at the cell borders. This makes this class of methods more diffusive, but also
computationally cheaper and much easier to parallelize with domain
decomposition codes. Short characteristics are used e.g. by code C2ray
(Mellema et al., 2006), code Athena (Davis et al., 2012), or code Pion (Mackey
et al., 2021). They have been also used together with the moment-based VET
method to calculate the local radiation pressure tensor $\mathbb{D}$ (see
§5.1.3).

Rijkhorst et al. (2006) developed the hybrid characteristics method, combining
the advantages of both long and short characteristics, particularly in
reducing the necessary parallel communication. This method is implemented
within the block-based daptive AMR code Flash. It divides the radiation field
into two components: (i) the local component within the blocks, calculated
using the long characteristics method, and (ii) the global component,
transporting radiation among blocks (and parallel sub-domains) using
interpolation similar to the short characteristics method. This algorithm was
further enhanced by Peters et al. (2010) and later re-implemented by
Buntemeyer et al. (2016) to also handle diffuse radiation and operate on GPU-
accelerated architectures.

Jiang (2021) developed an implicit radiation transport algorithm for MHD code
Athena++. It represents the radiation by its specific intensities along
discrete rays, and evolves them using a conservative finite volume method. It
avoids using a closure as in moment-based methods and the related artificial
diffusion. On the other hand, it solves time-dependent radiation transport
equation, contrary to majority of ray or characteristics based methods. The
original implementation was grey, however, Jiang (2022) extended it to include
the frequency dependence using the multi-group approach.

####  5.2.3 Reverse ray tracing

An alternative approach is the so-called reverse ray tracing. In contrast to
the "standard" forward ray tracing, which casts rays from the source outward,
reverse ray tracing schemes cast rays from the target point where the
radiation interacts with the gas (e.g., grid cell or SPH particle) and
transport the radiation inward to the target point. This approach offers
several advantages, including the highest density of rays (i.e., the highest
spatial resolution) near the location where radiation interacts with the gas.
Another advantage is its suitability for handling external diffuse radiation.
Hasegawa and Umemura (2010) implemented this algorithm into the SPH code
Start, incorporating tree-based acceleration to group distant sources into
oct-tree nodes. A grid-based version of this algorithm was developed by
Okamoto et al. (2012) for the code Argot. Altay and Theuns (2013) created the
reverse ray tracing code Urchin and integrated it into an SPH code focused on
cosmological simulations. Another implementation of this approach is Trevr for
the SPH code Gasoline by Grond et al. (2019), who introduced adaptive
refinement along rays and directionally dependent absorption coefficients for
tree nodes. Wünsch et al. (2021) developed the code Treeray, implementing
tree-accelerated reverse ray tracing for the Adaptive Mesh Refinement (AMR)
code Flash. It utilizes a fast parallel tree solver described in Wünsch et al.
(2018) and is coupled with the chemical network developed by Glover and Mac
Low (2007). Later, it was extended by Klepitko et al. (2023), who added
diffuse grey infrared radiation and related radiation pressure, and by Gaches
et al. (2023), who added frequency-dependent X-ray radiation. Treeray has been
used in the SILCC project (Walch et al., 2015, and follow-up papers) using
RMHD simulations to model the full cycle of molecular clouds in a part of a
galactic disc.

####  5.2.4 Monte Carlo methods

Monte Carlo Radiation Transport (MCRT) algorithms belong to a broad category
of statistical sampling methods that share this name. They simulate the
movement of individual photons through a medium by stochastically sampling
random events, including interactions with matter. Given that tracking
individual photons is prohibitively expensive, MCRT algorithms introduce
"photon packets" that carry information about the position, direction,
frequency distribution, etc., of the group of photons they represent. This
approach closely mirrors the reality of propagating radiation, making it
highly versatile, and facilitating the inclusion of various complex processes
such as scattering.

However, the statistical nature of this method introduces Poisson noise, and
maintaining it below a required level often necessitates a very large number
of photon packets. Consequently, MCRT codes have historically mostly been used
for static problems (e.g., the code Mocassin by Ercolano et al. 2003, the code
Hyperion by Robitaille 2011, or the code Radmc3D by Dullemond et al. 2012).
Due to the potential long distances photon packets can travel, Monte Carlo
methods are relatively challenging to parallelize using domain decomposition.
On the other hand, thread-based parallelization on shared memory architectures
is typically straightforward.

As this method does not rely on any grid, it is natural to use it with SPH or
unstructured grid codes. An algorithm to incorporate MCRT into SPH codes was
proposed by Lucy (1999). A Monte Carlo-based ray-tracing algorithm was
developed by Altay et al. (2008) and implemented into the publicly available
code Sphray, primarily focused on cosmological simulations. Another
implementation of MCRT into an SPH code was done by Nayakshin et al. (2009),
who used it to calculate the radiation pressure force. Roth and Kasen (2015)
developed a 1D MCRT code and coupled it with both Eulerian and Lagrangian
hydrodynamic codes. The first star formation simulations with MCRT were
performed by Lomax and Whitworth (2016), who developed the code Spamcart
coupled with SPH hydrodynamics. They used it to calculate synthetic intensity
maps and spectra of an embedded protostellar multiple system. Vandenbroucke
and Wood (2018) developed the radiation hydrodynamic code CMacIonize, which
works on an unstructured moving grid and is based on the Monte Carlo method
calculating $\nu$-dependent transport of ionizing photons. The code is
parallelized using a task-based scheme and has an option to work in a
distributed memory configuration. Code Torus, developed by Harries et al.
(2019), is an MCRT that includes native AMR grid-based hydrodynamics and can
also be coupled with the SPH code Sphng. It implements numerous microphysics
modules, including atomic and molecular line transport in moving media, dust
radiative equilibrium, photoionization equilibrium, and time-dependent
radiative transfer, and has been used e.g. for modelling massive stars
feedback in turbulent clouds (Ali et al., 2018). It is also designed for
postprocessing the output of several other (magneto-)hydrodynamic codes.
Recently, MCRT has also been implemented into the unstructured moving mesh
code Arepo by Smith et al. (2020).

Reference | Name/code |  $t$-dep | DD | MF | Note
---|---|---|---|---|---
Abel and Wandelt (2002) | - | N | N | ion. | gen., adaptive RT with HealPix
Mellema et al. (2006) | C2ray | N | N | ion. | grid, short chars.
Rijkhorst et al. (2006) | Flash | N | Y | grey | grid, hybrid characteristics
Dale et al. (2007) | Sphng | N | Y | ion. | SPH, adaptive RT
Altay et al. (2008) | Sphray | N | N |  $\nu$-dep | SPH, Monte Carlo
Pawlik and Schaye (2008) | Traphic/Gadget2 | Y | Y |  $\nu$-dep | SPH, photon packets, cones
Bisbas et al. (2009) | Seren | N | N | ion. | SPH, adaptive RT
Nayakshin et al. (2009) | - | Y | N | grey | SPH, rad. pressure
Hasegawa and Umemura (2010) | Start | N | N | ion. | SPH, reverse RT, tree accel.
Wise and Abel (2011) | Moray/Enzo | Y | Y | MG | grid, adaptive RT
Okamoto et al. (2012) | Argot | N | Y | ion. | grid-based version of START
Clark et al. (2012) | Treecol | N |  |  | reverse RT / column densities
Davis et al. (2012) | Athena | N | Y | grey | grid, short chars.
Altay and Theuns (2013) | Urchin | N | N |  $\nu$-dep | SPH, reverse RT
Baczynski et al. (2015) | Fervent/Flash | N | Y | MG | grid, forward RT
Bisbas et al. (2015a) | Torus-3Dpdr | N | Y |  $\nu$-dep | SPH, adp. RT (FUV) + MC (ion.&diff.)
Roth and Kasen (2015) | - | Y | N |  $\nu$-dep | grid, 1D Eulerian and Lagrangian
Buntemeyer et al. (2016) | Flash | N | Y | grey | grid, hybrid chars on adaptive mesh, GPU
Lomax and Whitworth (2016) | Spamcart | Y | N |  $\nu$-dep | SPH, first MCRT SF sims.
Rosen et al. (2017) | Harm2/Orion | N+Y | Y |  $\nu$-dep |  grid, forward RT ($\nu$-dep) + grey FLD
Frostholm et al. (2018) | Lampray/Ramses | N | Y | MG | grid, long chars., coupled to KROME
Jaura et al. (2018) | Sprai/Arepo | N | Y | ion. | unstruct. grid, SimpleX2 alg.
Vandenbroucke and Wood (2018) | CMacIonize | N | Y |  $\nu$-dep | unstruct. grid, Monte Carlo
Grond et al. (2019) | Trevr/Gasoline | N | Y | ion. | SPH, reverse RT, tree
Harries et al. (2019) | Torus | Y | Y |  $\nu$-dep | grid, Monte Carlo
Smith et al. (2020) | Arepo-MCRT | Y | Y | grey | unstruct. grid, Monte Carlo
Jiang (2021) | Athena++ | Y | Y | MG | grid, rays + FVM, implicit
Mackey et al. (2021) | Pion | N | Y | MG | grid, short chars.
Wünsch et al. (2021) | TreeRay/Flash | N | Y | MG | grid, reverse RT, tree
Table 4: Radiation transport codes using ray tracing and methods of
characteristics. The first column gives the reference to the work describing
the code or algorithm; the second column give the name of the radiation
transport code or RHD code (if they exist); the third column whether the code
solves time-dependent RTE; the fourth column tells whether the code is
parallel using the domain decomposition; the fifth column gives how the code
deals with the frequency dependence (grey, ion. - a single band of ionising
photons, MG - multi-group approach, $\nu$-dep - full frequency dependence);
and the last column gives eventually an additional information.

###  5.3 Other approximations

Many studies conducting hydrodynamic simulations of star formation further
simplify the radiation transport problem. For example, Vázquez-Semadeni et al.
(2017) include stellar ionizing radiation feedback by calculating the
Strömgren sphere radius from the EUV photon production rate of a star and the
characteristic density obtained as the geometric mean of the density at the
star and at the target cell. Subsequently, they set the gas temperature within
the sphere to $10^{4}$ K.

A specific case involves non-ionizing diffuse ambient interstellar radiation
fields, crucial for the physics of molecular clouds as they provide heating
and affect molecular gas chemistry. In molecular clouds, this radiation
intensity strongly depends on shielding, i.e., the optical depth between a
given location and the edge of the cloud. This consideration led Clark et al.
(2012) to develop the Treecol algorithm, implemented into the code Gadget2. It
efficiently calculates optical depths in all directions (discretized by
HealPix) for each SPH particle using a tree, enabling the determination of
heating and dissociation rates. This algorithm was later implemented by
Valdivia and Hennebelle (2014) into Ramses and by Wünsch et al. (2018) into
Flash.

A different approach to solving a very similar problem was suggested by
Stamatellos et al. (2007), who estimate the mean optical depth from local
density, temperature, and gravitational potential using a set of precalculated
models of collapsing clouds. This method provides reasonably good results for
almost zero computational costs.

The simplest treatment of radiation is to modify the gas equation of state to
calculate the equilibrium gas temperature from local gas density, thus
accounting for the most probable result of the heating-cooling balance. A
popular choice is the barotropic equation of state (e.g., Bonnell, 1994;
Whitworth et al., 1995; Bate, 1998), a simple analytical formula designed to
mimic the thermodynamics of spherically symmetric collapse of a single,
isolated protostar.

##  6 Summary

This work reviews methods to solve radiation transport equation (RTE) in
simulations of star formation and related fields. According to the algorithm,
they can be roughly divided into two groups: moment-based methods and methods
solving the RTE directly.

Moment-based methods approximate the radiation specific intensity by
integrating it over all directions and expressing it in terms of several
moments. The order of the moments in the sequence closure determines the level
of approximation, leading to artificial diffusion of the radiation, unphysical
self-interaction, and other artificial effects. These methods typically
include the time-dependent term of the RTE and sometimes artificially decrease
the speed of light to reduce computational costs. Their computational expenses
are independent of the number of sources, making them well-suited for
calculating diffuse radiation. Frequency-wise, moment-based methods often use
the grey approximation or, occasionally, the multi-group frequency approach.
They are sometimes combined with ray tracing, which calculates radiation from
one or several discrete sources, while the moment-based method addresses the
diffuse radiation component. The most commonly used methods of this type
include flux-limited diffusion, the M1 closure, and the variable Eddington
tensor method.

Methods that directly solve the RTE integrate it along specific lines,
commonly referred to as rays. These methods include forward and reverse ray
tracing, long, hybrid, and short characteristics, as well as Monte Carlo
methods. Ray tracing and characteristics-based methods typically do not
include the time-dependent RTE term, utilizing the infinite speed of light
approximation; however, there are exceptions. Monte Carlo methods are usually
time-dependent. Depending on angular discretization, ray tracing and
long/hybrid characteristics can achieve high accuracy but at the expense of
high computational costs. Short characteristics methods tend to be more
diffusive. While Monte Carlo methods can be highly accurate, achieving this
requires a substantial number of photon packets, leading to increased
computational costs. Monte Carlo methods are also versatile, suitable e.g. for
treating radiation scattering. Computational costs for most of these methods
scale with the number of sources, with exceptions like reverse ray tracing
methods that group distant sources and methods combining ray tracing with
advection schemes. Various treatments of radiation frequency dependence are
employed; generally, methods designed for a small number of sources can
incorporate detailed frequency dependence. Forward ray tracing and long
characteristics methods are typically challenging to parallelize using domain
decomposition, whereas short characteristics are relatively straightforward.
Hybrid characteristics, reverse ray tracing, and Monte Carlo methods fall in
the middle ground regarding parallelization difficulty.

Radiation transport has become a well-established component of star formation
simulations, with many methods and codes available, continuously advancing in
both maturity and complexity. Consequently, it becomes imperative to not only
validate these codes against standard tests but also to benchmark them against
each other and, ideally, against realistic problems inherent in the realm of
star formation. In this regard, initiatives such as the code comparison
projects, such as those conducted by Iliev et al. (2006, 2009) for
cosmological radiation transport codes and the STARBENCH project led by Bisbas
et al. (2015b) that compares radiation transport codes for star formation,
prove to be very valuable. Moreover, there is a growing need for more
comparison projects to enhance our understanding and confidence in the
capabilities of these codes.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Funding

This study has been supported by the institutional project RVO:67985815.

## Acknowledgments

I thank the anonymous referees for their constructive comments that have
helped to improve this work. This study has been supported by the
institutional project RVO:67985815.

## References

  * Abel et al. (1999) Abel, T., Norman, M. L., and Madau, P. (1999).  Photon-conserving Radiative Transfer around Point Sources in Multidimensional Numerical Cosmology.  _ApJ_ 523, 66–71.  10.1086/307739
  * Abel and Wandelt (2002) Abel, T. and Wandelt, B. D. (2002).  Adaptive ray tracing for radiative transfer around point sources.  _MNRAS_ 330, L53–L56.  10.1046/j.1365-8711.2002.05206.x
  * Ali et al. (2018) Ali, A., Harries, T. J., and Douglas, T. A. (2018).  Modelling massive star feedback with Monte Carlo radiation hydrodynamics: photoionization and radiation pressure in a turbulent cloud.  _MNRAS_ 477, 5422–5436.  10.1093/mnras/sty1001
  * Alme and Wilson (1973) Alme, M. L. and Wilson, J. R. (1973).  X-Ray Emission from a Neutron Star Accreting Material.  _ApJ_ 186, 1015–1026.  10.1086/152566
  * Altay et al. (2008) Altay, G., Croft, R. A. C., and Pelupessy, I. (2008).  SPHRAY: a smoothed particle hydrodynamics ray tracer for radiative transfer.  _MNRAS_ 386, 1931–1946.  10.1111/j.1365-2966.2008.13212.x
  * Altay and Theuns (2013) Altay, G. and Theuns, T. (2013).  URCHIN: a reverse ray tracer for astrophysical applications.  _MNRAS_ 434, 748–764.  10.1093/mnras/stt1067
  * Asahina et al. (2020) Asahina, Y., Takahashi, H. R., and Ohsuga, K. (2020).  A Numerical Scheme for General Relativistic Radiation Magnetohydrodynamics Based on Solving a Grid-based Boltzmann Equation.  _ApJ_ 901, 96.  10.3847/1538-4357/abaf51
  * Aubert and Teyssier (2008) Aubert, D. and Teyssier, R. (2008).  A radiative transfer scheme for cosmological reionization based on a local Eddington tensor.  _MNRAS_ 387, 295–307.  10.1111/j.1365-2966.2008.13223.x
  * Auer and Mihalas (1970) Auer, L. H. and Mihalas, D. (1970).  On the use of variable Eddington factors in non-LTE stellar atmospheres computations.  _MNRAS_ 149, 65.  10.1093/mnras/149.1.65
  * Baczynski et al. (2015) Baczynski, C., Glover, S. C. O., and Klessen, R. S. (2015).  <monospace>Fervent</monospace>: chemistry-coupled, ionizing and non-ionizing radiative feedback in hydrodynamical simulations.  _MNRAS_ 454, 380–411.  10.1093/mnras/stv1906
  * Bassett et al. (2021) Bassett, B. R., Owen, J. M., and Brunner, T. A. (2021).  Efficient smoothed particle radiation hydrodynamics I: Thermal radiative transfer.  _Journal of Computational Physics_ 429, 109996.  10.1016/j.jcp.2020.109996
  * Bate (1998) Bate, M. R. (1998).  Collapse of a Molecular Cloud Core to Stellar Densities: The First Three-dimensional Calculations.  _ApJ Let_ 508, L95–L98.  10.1086/311719
  * Bate (2019) Bate, M. R. (2019).  The statistical properties of stars and their dependence on metallicity.  _MNRAS_ 484, 2341–2361.  10.1093/mnras/stz103
  * Bate and Keto (2015) Bate, M. R. and Keto, E. R. (2015).  Combining radiative transfer and diffuse interstellar medium physics to model star formation.  _MNRAS_ 449, 2643–2667.  10.1093/mnras/stv451
  * Benz (1990) Benz, W. (1990).  Smooth Particle Hydrodynamics - a Review.  In _Numerical Modelling of Nonlinear Stellar Pulsations Problems and Prospects_ , ed. J. R. Buchler. 269
  * Bisbas et al. (2015a) Bisbas, T. G., Haworth, T. J., Barlow, M. J., Viti, S., Harries, T. J., Bell, T., et al. (2015a).  TORUS-3DPDR: a self-consistent code treating three-dimensional photoionization and photodissociation regions.  _MNRAS_ 454, 2828–2843.  10.1093/mnras/stv2156
  * Bisbas et al. (2015b) Bisbas, T. G., Haworth, T. J., Williams, R. J. R., Mackey, J., Tremblin, P., Raga, A. C., et al. (2015b).  STARBENCH: the D-type expansion of an H II region.  _MNRAS_ 453, 1324–1343.  10.1093/mnras/stv1659
  * Bisbas et al. (2009) Bisbas, T. G., Wünsch, R., Whitworth, A. P., and Hubber, D. A. (2009).  Smoothed particle hydrodynamics simulations of expanding H II regions. I. Numerical method and applications.  _A &A_ 497, 649–659.  10.1051/0004-6361/200811522
  * Bitsch et al. (2013) Bitsch, B., Crida, A., Morbidelli, A., Kley, W., and Dobbs-Dixon, I. (2013).  Stellar irradiated discs and implications on migration of embedded planets. I. Equilibrium discs.  _A &A_ 549, A124.  10.1051/0004-6361/201220159
  * Bloch et al. (2021) Bloch, H., Tremblin, P., González, M., Padioleau, T., and Audit, E. (2021).  A high-performance and portable asymptotic preserving radiation hydrodynamics code with the M1 model.  _A &A_ 646, A123.  10.1051/0004-6361/202038579
  * Bodenheimer et al. (1990) Bodenheimer, P., Yorke, H. W., Rozyczka, M., and Tohline, J. E. (1990).  The Formation Phase of the Solar Nebula.  _ApJ_ 355, 651.  10.1086/168798
  * Boley et al. (2007) Boley, A. C., Durisen, R. H., Nordlund, Å., and Lord, J. (2007).  Three-Dimensional Radiative Hydrodynamics for Disk Stability Simulations: A Proposed Testing Standard and New Results.  _ApJ_ 665, 1254–1267.  10.1086/519767
  * Bonnell (1994) Bonnell, I. A. (1994).  A New Binary Formation Mechanism.  _MNRAS_ 269, 837–848.  10.1093/mnras/269.3.837
  * Boss (2001) Boss, A. P. (2001).  Gas Giant Protoplanet Formation: Disk Instability Models with Thermodynamics and Radiative Transfer.  _ApJ_ 563, 367–373.  10.1086/323694
  * Brookshaw (1985) Brookshaw, L. (1985).  A method of calculating radiative heat diffusion in particle simulations.  _PASA_ 6, 207–210.  10.1017/S1323358000018117
  * Buntemeyer et al. (2016) Buntemeyer, L., Banerjee, R., Peters, T., Klassen, M., and Pudritz, R. E. (2016).  Radiation hydrodynamics using characteristics on adaptive decomposed domains for massively parallel star formation simulations.  _New A_ 43, 49–69.  10.1016/j.newast.2015.07.002
  * Castor (2004) Castor, J. I. (2004).  _Radiation Hydrodynamics_
  * Chan et al. (2021) Chan, T. K., Theuns, T., Bower, R., and Frenk, C. (2021).  Smoothed particle radiation hydrodynamics: two-moment method with local Eddington tensor closure.  _MNRAS_ 505, 5784–5814.  10.1093/mnras/stab1686
  * Chatzopoulos and Weide (2019) Chatzopoulos, E. and Weide, K. (2019).  Gray Radiation Hydrodynamics with the FLASH Code for Astrophysical Applications.  _ApJ_ 876, 148.  10.3847/1538-4357/ab18f9
  * Clark et al. (2012) Clark, P. C., Glover, S. C. O., and Klessen, R. S. (2012).  TreeCol: a novel approach to estimating column densities in astrophysical simulations.  _MNRAS_ 420, 745–756.  10.1111/j.1365-2966.2011.20087.x
  * Clementel et al. (2014) Clementel, N., Madura, T. I., Kruip, C. J. H., Icke, V., and Gull, T. R. (2014).  3D radiative transfer in $\eta$ Carinae: application of the SIMPLEX algorithm to 3D SPH simulations of binary colliding winds.  _MNRAS_ 443, 2475–2491.  10.1093/mnras/stu1287
  * Cohen et al. (1999) Cohen, J., Pratchett, T., and Steward, I. (1999).  _The Science of Discworld_ (Ebury Press)
  * Commerçon et al. (2014) Commerçon, B., Debout, V., and Teyssier, R. (2014).  A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids.  _A &A_ 563, A11.  10.1051/0004-6361/201322858
  * Commerçon et al. (2011) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., and Chabrier, G. (2011).  Radiation hydrodynamics with adaptive mesh refinement and application to prestellar core collapse. I. Methods.  _A &A_ 529, A35.  10.1051/0004-6361/201015880
  * Dale et al. (2007) Dale, J. E., Ercolano, B., and Clarke, C. J. (2007).  A new algorithm for modelling photoionizing radiation in smoothed particle hydrodynamics.  _MNRAS_ 382, 1759–1767.  10.1111/j.1365-2966.2007.12486.x
  * Davis et al. (2012) Davis, S. W., Stone, J. M., and Jiang, Y.-F. (2012).  A Radiation Transfer Solver for Athena Using Short Characteristics.  _ApJS_ 199, 9.  10.1088/0067-0049/199/1/9
  * Delaunay (1934) Delaunay, B. (1934).  Sur la sphère vide. A la mémoire de Georges Voronoï.  _Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et naturelles_ 6, 793–800
  * Dubroca and Feugeas (1999) Dubroca, B. and Feugeas, J. (1999).  Etude théorique et numérique d’une hiérarchie de modèles aux moments pour le transfert radiatif.  _Academie des Sciences Paris Comptes Rendus Serie Sciences Mathematiques_ 329, 915–920.  10.1016/S0764-4442(00)87499-6
  * Dullemond et al. (2012) [Dataset] Dullemond, C. P., Juhasz, A., Pohl, A., Sereshti, F., Shetty, R., Peters, T., et al. (2012).  RADMC-3D: A multi-purpose radiative transfer tool.  Astrophysics Source Code Library, record ascl:1202.015
  * Eddington (1926) Eddington, A. S. (1926).  _The Internal Constitution of the Stars_
  * Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., and Liu, X. W. (2003).  MOCASSIN: a fully three-dimensional Monte Carlo photoionization code.  _MNRAS_ 340, 1136–1152.  10.1046/j.1365-8711.2003.06371.x
  * Flock et al. (2013) Flock, M., Fromang, S., González, M., and Commerçon, B. (2013).  Radiation magnetohydrodynamics in global simulations of protoplanetary discs.  _A &A_ 560, A43.  10.1051/0004-6361/201322451
  * Frostholm et al. (2018) Frostholm, T., Haugbølle, T., and Grassi, T. (2018).  Lampray: Multi-group long characteristics ray tracing for adaptive mesh radiation hydrodynamics.  _arXiv e-prints_ , arXiv:1809.0554110.48550/arXiv.1809.05541
  * Fryer et al. (2006) Fryer, C. L., Rockefeller, G., and Warren, M. S. (2006).  SNSPH: A Parallel Three-dimensional Smoothed Particle Radiation Hydrodynamics Code.  _ApJ_ 643, 292–305.  10.1086/501493
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., et al. (2000).  FLASH: An Adaptive Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear Flashes.  _ApJS_ 131, 273–334.  10.1086/317361
  * Gaches et al. (2023) Gaches, B. A. L., Walch, S., Wünsch, R., and Mackey, J. (2023).  Tree-based solvers for adaptive mesh refinement code FLASH - IV. An X-ray radiation scheme to couple discrete and diffuse X-ray emission sources to the thermochemistry of the interstellar medium.  _MNRAS_ 522, 4674–4690.  10.1093/mnras/stad1206
  * Garcia-Segura and Franco (1996) Garcia-Segura, G. and Franco, J. (1996).  From Ultracompact to Extended H II Regions.  _ApJ_ 469, 171.  10.1086/177769
  * Gehmeyr and Mihalas (1994) Gehmeyr, M. and Mihalas, D. (1994).  Adaptive grid radiation hydrodynamics with TITAN.  _Physica D Nonlinear Phenomena_ 77, 320–341.  10.1016/0167-2789(94)90143-0
  * Gittings et al. (2008) Gittings, M., Weaver, R., Clover, M., Betlach, T., Byrne, N., Coker, R., et al. (2008).  The RAGE radiation-hydrodynamic code.  _Computational Science and Discovery_ 1, 015005.  10.1088/1749-4699/1/1/015005
  * Glover and Mac Low (2007) Glover, S. C. O. and Mac Low, M.-M. (2007).  Simulating the Formation of Molecular Clouds. I. Slow Formation by Gravitational Collapse from Static Initial Conditions.  _ApJS_ 169, 239–268.  10.1086/512238
  * Gnedin and Abel (2001) Gnedin, N. Y. and Abel, T. (2001).  Multi-dimensional cosmological radiative transfer with a Variable Eddington Tensor formalism.  _New A_ 6, 437–455.  10.1016/S1384-1076(01)00068-9
  * González et al. (2007) González, M., Audit, E., and Huynh, P. (2007).  HERACLES: a three-dimensional radiation hydrodynamics code.  _A &A_ 464, 429–435.  10.1051/0004-6361:20065486
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., Wandelt, B. D., Hansen, F. K., Reinecke, M., et al. (2005).  HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere.  _ApJ_ 622, 759–771.  10.1086/427976
  * Grassi et al. (2014) Grassi, T., Bovino, S., Schleicher, D. R. G., Prieto, J., Seifried, D., Simoncini, E., et al. (2014).  KROME - a package to embed chemistry in astrophysical simulations.  _MNRAS_ 439, 2386–2419.  10.1093/mnras/stu114
  * Grond et al. (2019) Grond, J. J., Woods, R. M., Wadsley, J. W., and Couchman, H. M. P. (2019).  TREVR: A general N log2N radiative transfer algorithm.  _MNRAS_ 485, 3681–3695.  10.1093/mnras/stz525
  * Harries et al. (2019) Harries, T. J., Haworth, T. J., Acreman, D., Ali, A., and Douglas, T. (2019).  The TORUS radiation transfer code.  _Astronomy and Computing_ 27, 63.  10.1016/j.ascom.2019.03.002
  * Hasegawa and Umemura (2010) Hasegawa, K. and Umemura, M. (2010).  START: smoothed particle hydrodynamics with tree-based accelerated radiative transfer.  _MNRAS_ 407, 2632–2644.  10.1111/j.1365-2966.2010.17100.x
  * Hayes and Norman (2003) Hayes, J. C. and Norman, M. L. (2003).  Beyond Flux-limited Diffusion: Parallel Algorithms for Multidimensional Radiation Hydrodynamics.  _ApJS_ 147, 197–220.  10.1086/374658
  * Hofmann (1990) Hofmann, G. R. (1990).  Who invented ray tracing?  _The Visual Computer_ 6, 120–124
  * Hopkins and Grudić (2019) Hopkins, P. F. and Grudić, M. Y. (2019).  Numerical problems in coupling photon momentum (radiation pressure) to gas.  _MNRAS_ 483, 4187–4196.  10.1093/mnras/sty3089
  * Iliev et al. (2006) Iliev, I. T., Ciardi, B., Alvarez, M. A., Maselli, A., Ferrara, A., Gnedin, N. Y., et al. (2006).  Cosmological radiative transfer codes comparison project - I. The static density field tests.  _MNRAS_ 371, 1057–1086.  10.1111/j.1365-2966.2006.10775.x
  * Iliev et al. (2009) Iliev, I. T., Whalen, D., Mellema, G., Ahn, K., Baek, S., Gnedin, N. Y., et al. (2009).  Cosmological radiative transfer comparison project - II. The radiation-hydrodynamic tests.  _MNRAS_ 400, 1283–1316.  10.1111/j.1365-2966.2009.15558.x
  * Jaura et al. (2018) Jaura, O., Glover, S. C. O., Klessen, R. S., and Paardekooper, J. P. (2018).  SPRAI: coupling of radiative feedback and primordial chemistry in moving mesh hydrodynamics.  _MNRAS_ 475, 2822–2834.  10.1093/mnras/stx3356
  * Jiang (2021) Jiang, Y.-F. (2021).  An Implicit Finite Volume Scheme to Solve the Time-dependent Radiation Transport Equation Based on Discrete Ordinates.  _ApJS_ 253, 49.  10.3847/1538-4365/abe303
  * Jiang (2022) Jiang, Y.-F. (2022).  Multigroup Radiation Magnetohydrodynamics Based on Discrete Ordinates including Compton Scattering.  _ApJS_ 263, 4.  10.3847/1538-4365/ac9231
  * Jiang et al. (2012) Jiang, Y.-F., Stone, J. M., and Davis, S. W. (2012).  A Godunov Method for Multidimensional Radiation Magnetohydrodynamics Based on a Variable Eddington Tensor.  _ApJS_ 199, 14.  10.1088/0067-0049/199/1/14
  * Kannan et al. (2019) Kannan, R., Vogelsberger, M., Marinacci, F., McKinnon, R., Pakmor, R., and Springel, V. (2019).  AREPO-RT: radiation hydrodynamics on a moving mesh.  _MNRAS_ 485, 117–149.  10.1093/mnras/stz287
  * Kessel-Deynet and Burkert (2000) Kessel-Deynet, O. and Burkert, A. (2000).  Ionizing radiation in smoothed particle hydrodynamics.  _MNRAS_ 315, 713–721.  10.1046/j.1365-8711.2000.03451.x
  * Klahr et al. (1999) Klahr, H. H., Henning, T., and Kley, W. (1999).  On the Azimuthal Structure of Thermal Convection in Circumstellar Disks.  _ApJ_ 514, 325–343.  10.1086/306926
  * Klassen et al. (2014) Klassen, M., Kuiper, R., Pudritz, R. E., Peters, T., Banerjee, R., and Buntemeyer, L. (2014).  A General Hybrid Radiation Transport Scheme for Star Formation Simulations on an Adaptive Grid.  _ApJ_ 797, 4.  10.1088/0004-637X/797/1/4
  * Klepitko et al. (2023) Klepitko, A., Walch, S., Wünsch, R., Seifried, D., Dinnbier, F., and Haid, S. (2023).  Tree-based solvers for adaptive mesh refinement code FLASH - III: a novel scheme for radiation pressure on dust and gas and radiative transfer from diffuse sources.  _MNRAS_ 521, 160–184.  10.1093/mnras/stad385
  * Kley (1989) Kley, W. (1989).  Radiation hydrodynamics of the boundary layer in accretion disks. I \- Numerical methods.  _A &A_ 208, 98–110
  * Kolb et al. (2013) Kolb, S. M., Stute, M., Kley, W., and Mignone, A. (2013).  Radiation hydrodynamics integrated in the PLUTO code.  _A &A_ 559, A80.  10.1051/0004-6361/201321499
  * Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., McKee, C. F., and Bolstad, J. (2007).  Equations and Algorithms for Mixed-frame Flux-limited Diffusion Radiation Hydrodynamics.  _ApJ_ 667, 626–643.  10.1086/520791
  * Kuiper et al. (2010) Kuiper, R., Klahr, H., Dullemond, C., Kley, W., and Henning, T. (2010).  Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation.  _A &A_ 511, A81.  10.1051/0004-6361/200912355
  * Levermore (1984) Levermore, C. D. (1984).  Relating Eddington factors to flux limiters.  _J. Quant. Spec. Radiat. Transf._ 31, 149–160.  10.1016/0022-4073(84)90112-2
  * Levermore and Pomraning (1981) Levermore, C. D. and Pomraning, G. C. (1981).  A flux-limited diffusion theory.  _ApJ_ 248, 321–334.  10.1086/159157
  * Lomax and Whitworth (2016) Lomax, O. and Whitworth, A. P. (2016).  SPAMCART: a code for smoothed particle Monte Carlo radiative transfer.  _MNRAS_ 461, 3542–3551.  10.1093/mnras/stw1568
  * Lucy (1977) Lucy, L. B. (1977).  A numerical approach to the testing of the fission hypothesis.  _AJ_ 82, 1013–1024.  10.1086/112164
  * Lucy (1999) Lucy, L. B. (1999).  Computing radiative equilibria with Monte Carlo techniques.  _A &A_ 344, 282–288
  * Mackey et al. (2021) Mackey, J., Green, S., Moutzouri, M., Haworth, T. J., Kavanagh, R. D., Zargaryan, D., et al. (2021).  PION: simulating bow shocks and circumstellar nebulae.  _MNRAS_ 504, 983–1008.  10.1093/mnras/stab781
  * Mayer et al. (2007) Mayer, L., Lufkin, G., Quinn, T., and Wadsley, J. (2007).  Fragmentation of Gravitationally Unstable Gaseous Protoplanetary Disks with Radiative Transfer.  _ApJ Let_ 661, L77–L80.  10.1086/518433
  * Mellema et al. (2006) Mellema, G., Iliev, I. T., Alvarez, M. A., and Shapiro, P. R. (2006).  C 2-ray: A new method for photon-conserving transport of ionizing radiation.  _New A_ 11, 374–395.  10.1016/j.newast.2005.09.004
  * Melon Fuksman and Mignone (2019) Melon Fuksman, J. D. and Mignone, A. (2019).  A Radiative Transfer Module for Relativistic Magnetohydrodynamics in the PLUTO Code.  _ApJS_ 242, 20.  10.3847/1538-4365/ab18ff
  * Mihalas and Klein (1982) Mihalas, D. and Klein, R. I. (1982).  On the solution of the time-dependent inertial-frame equation of radiative transfer in moving media to O(v/c).  _Journal of Computational Physics_ 46, 97–137.  10.1016/0021-9991(82)90007-9
  * Mihalas and Mihalas (1984) Mihalas, D. and Mihalas, B. W. (1984).  _Foundations of radiation hydrodynamics_
  * Minerbo (1978) Minerbo, G. N. (1978).  Maximum entropy Eddington factors.  _J. Quant. Spec. Radiat. Transf._ 20, 541–545.  10.1016/0022-4073(78)90024-9
  * Moens et al. (2022) Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R. (2022).  Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion.  _A &A_ 657, A81.  10.1051/0004-6361/202141023
  * Murray et al. (1994) Murray, S. D., Castor, J. I., Klein, R. I., and McKee, C. F. (1994).  Accretion Disk Coronae in High-Luminosity Systems.  _ApJ_ 435, 631.  10.1086/174842
  * Nayakshin et al. (2009) Nayakshin, S., Cha, S.-H., and Hobbs, A. (2009).  Dynamic Monte Carlo radiation transfer in SPH: radiation pressure force implementation.  _MNRAS_ 397, 1314–1325.  10.1111/j.1365-2966.2009.15091.x
  * Okamoto et al. (2012) Okamoto, T., Yoshikawa, K., and Umemura, M. (2012).  ARGOT: accelerated radiative transfer on grids using oct-tree.  _MNRAS_ 419, 2855–2866.  10.1111/j.1365-2966.2011.19927.x
  * Paardekooper et al. (2010) Paardekooper, J. P., Kruip, C. J. H., and Icke, V. (2010).  SimpleX2: radiative transfer on an unstructured, dynamic grid.  _A &A_ 515, A79.  10.1051/0004-6361/200913821
  * Pawlik and Schaye (2008) Pawlik, A. H. and Schaye, J. (2008).  TRAPHIC - radiative transfer for smoothed particle hydrodynamics simulations.  _MNRAS_ 389, 651–677.  10.1111/j.1365-2966.2008.13601.x
  * Peters et al. (2010) Peters, T., Banerjee, R., Klessen, R. S., Mac Low, M.-M., Galván-Madrid, R., and Keto, E. R. (2010).  H II Regions: Witnesses to Massive Star Formation.  _ApJ_ 711, 1017–1028.  10.1088/0004-637X/711/2/1017
  * Petkova and Springel (2009) Petkova, M. and Springel, V. (2009).  An implementation of radiative transfer in the cosmological simulation code GADGET.  _MNRAS_ 396, 1383–1403.  10.1111/j.1365-2966.2009.14843.x
  * Price and Bate (2009) Price, D. J. and Bate, M. R. (2009).  Inefficient star formation: the combined effects of magnetic fields and radiative feedback.  _MNRAS_ 398, 33–46.  10.1111/j.1365-2966.2009.14969.x
  * Ramsey and Dullemond (2015) Ramsey, J. P. and Dullemond, C. P. (2015).  Radiation hydrodynamics including irradiation and adaptive mesh refinement with AZEuS. I. Methods.  _A &A_ 574, A81.  10.1051/0004-6361/201424954
  * Razoumov and Scott (1999) Razoumov, A. O. and Scott, D. (1999).  Three-dimensional numerical cosmological radiative transfer in an inhomogeneous medium.  _MNRAS_ 309, 287–298.  10.1046/j.1365-8711.1999.02775.x
  * Rijkhorst et al. (2006) Rijkhorst, E. J., Plewa, T., Dubey, A., and Mellema, G. (2006).  Hybrid characteristics: 3D radiative transfer for parallel adaptive mesh refinement hydrodynamics.  _A &A_ 452, 907–920.  10.1051/0004-6361:20053401
  * Ritzerveld and Icke (2006) Ritzerveld, J. and Icke, V. (2006).  Transport on adaptive random lattices.  _Phys. Rev. E_ 74, 026704.  10.1103/PhysRevE.74.026704
  * Robitaille (2011) Robitaille, T. P. (2011).  HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code.  _A &A_ 536, A79.  10.1051/0004-6361/201117150
  * Rosdahl et al. (2013) Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., and Teyssier, R. (2013).  RAMSES-RT: radiation hydrodynamics in the cosmological context.  _MNRAS_ 436, 2188–2231.  10.1093/mnras/stt1722
  * Rosdahl and Teyssier (2015) Rosdahl, J. and Teyssier, R. (2015).  A scheme for radiation pressure and photon diffusion with the M1 closure in RAMSES-RT.  _MNRAS_ 449, 4380–4403.  10.1093/mnras/stv567
  * Rosen et al. (2017) Rosen, A. L., Krumholz, M. R., Oishi, J. S., Lee, A. T., and Klein, R. I. (2017).  Hybrid Adaptive Ray-Moment Method (HARM2): A highly parallel method for radiation hydrodynamics on adaptive grids.  _Journal of Computational Physics_ 330, 924–942.  10.1016/j.jcp.2016.10.048
  * Roth and Kasen (2015) Roth, N. and Kasen, D. (2015).  Monte Carlo Radiation-Hydrodynamics With Implicit Methods.  _ApJS_ 217, 9.  10.1088/0067-0049/217/1/9
  * Rozyczka (1985) Rozyczka, M. (1985).  Two-dimensional models of stellar wind bubbles. I - Numerical methods and their application to the investigation of outer shell instabilities.  _A &A_ 143, 59–71
  * Sekora and Stone (2010) Sekora, M. D. and Stone, J. M. (2010).  A hybrid Godunov method for radiation hydrodynamics.  _Journal of Computational Physics_ 229, 6819–6852.  10.1016/j.jcp.2010.05.024
  * Shu (1991) Shu, F. H. (1991).  _The physics of astrophysics. Volume 1: Radiation._
  * Sijoy and Chaturvedi (2015) Sijoy, C. D. and Chaturvedi, S. (2015).  TRHD: Three-temperature radiation-hydrodynamics code with an implicit non-equilibrium radiation transport using a cell-centered monotonic finite volume scheme on unstructured-grids.  _Computer Physics Communications_ 190, 98–119.  10.1016/j.cpc.2015.01.019
  * Sądowski et al. (2013) Sądowski, A., Narayan, R., Tchekhovskoy, A., and Zhu, Y. (2013).  Semi-implicit scheme for treating radiation under M1 closure in general relativistic conservative fluid dynamics codes.  _MNRAS_ 429, 3533–3550.  10.1093/mnras/sts632
  * Skinner et al. (2019) Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., and Vartanyan, D. (2019).  FORNAX: A Flexible Code for Multiphysics Astrophysical Simulations.  _ApJS_ 241, 7.  10.3847/1538-4365/ab007f
  * Skinner and Ostriker (2013) Skinner, M. A. and Ostriker, E. C. (2013).  A Two-moment Radiation Hydrodynamics Module in Athena Using a Time-explicit Godunov Method.  _ApJS_ 206, 21.  10.1088/0067-0049/206/2/21
  * Smith et al. (2020) Smith, A., Kannan, R., Tsang, B. T. H., Vogelsberger, M., and Pakmor, R. (2020).  AREPO-MCRT: Monte Carlo Radiation Hydrodynamics on a Moving Mesh.  _ApJ_ 905, 27.  10.3847/1538-4357/abc47e
  * Stamatellos et al. (2007) Stamatellos, D., Whitworth, A. P., Bisbas, T., and Goodwin, S. (2007).  Radiative transfer and the energy equation in SPH simulations of star formation.  _A &A_ 475, 37–49.  10.1051/0004-6361:20077373
  * Steinacker et al. (2013) Steinacker, J., Baes, M., and Gordon, K. D. (2013).  Three-Dimensional Dust Radiative Transfer*.  _ARA & A_ 51, 63–104.  10.1146/annurev-astro-082812-141042
  * Stone et al. (1992) Stone, J. M., Mihalas, D., and Norman, M. L. (1992).  ZEUS-2D: A Radiation Magnetohydrodynamics Code for Astrophysical Flows in Two Space Dimensions. III. The Radiation Hydrodynamic Algorithms and Tests.  _ApJS_ 80, 819.  10.1086/191682
  * Teyssier and Commerçon (2019) Teyssier, R. and Commerçon, B. (2019).  Numerical Methods for Simulating Star Formation.  _Frontiers in Astronomy and Space Sciences_ 6, 51.  10.3389/fspas.2019.00051
  * Turner and Stone (2001) Turner, N. J. and Stone, J. M. (2001).  A Module for Radiation Hydrodynamic Calculations with ZEUS-2D Using Flux-limited Diffusion.  _ApJS_ 135, 95–107.  10.1086/321779
  * Valdivia and Hennebelle (2014) Valdivia, V. and Hennebelle, P. (2014).  A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds.  _A &A_ 571, A46.  10.1051/0004-6361/201423720
  * van der Holst et al. (2011) van der Holst, B., Tóth, G., Sokolov, I. V., Powell, K. G., Holloway, J. P., Myra, E. S., et al. (2011).  CRASH: A Block-adaptive-mesh Code for Radiative Shock Hydrodynamics—Implementation and Verification.  _ApJS_ 194, 23.  10.1088/0067-0049/194/2/23
  * Vandenbroucke and Wood (2018) Vandenbroucke, B. and Wood, K. (2018).  The Monte Carlo photoionization and moving-mesh radiation hydrodynamics code CMACIONIZE.  _Astronomy and Computing_ 23, 40.  10.1016/j.ascom.2018.02.005
  * Vázquez-Semadeni et al. (2017) Vázquez-Semadeni, E., González-Samaniego, A., and Colín, P. (2017).  Hierarchical star cluster assembly in globally collapsing molecular clouds.  _MNRAS_ 467, 1313–1328.  10.1093/mnras/stw3229
  * Viau et al. (2006) Viau, S., Bastien, P., and Cha, S.-H. (2006).  An Implicit Method for Radiative Transfer with the Diffusion Approximation in Smooth Particle Hydrodynamics.  _ApJ_ 639, 559–570.  10.1086/499328
  * Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., Gatto, A., Glover, S. C. O., Wünsch, R., et al. (2015).  The SILCC (SImulating the LifeCycle of molecular Clouds) project - I. Chemical evolution of the supernova-driven ISM.  _MNRAS_ 454, 238–268.  10.1093/mnras/stv1975
  * Whitehouse and Bate (2004) Whitehouse, S. C. and Bate, M. R. (2004).  Smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation.  _MNRAS_ 353, 1078–1094.  10.1111/j.1365-2966.2004.08131.x
  * Whitehouse et al. (2005) Whitehouse, S. C., Bate, M. R., and Monaghan, J. J. (2005).  A faster algorithm for smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation.  _MNRAS_ 364, 1367–1377.  10.1111/j.1365-2966.2005.09683.x
  * Whitworth et al. (1995) Whitworth, A. P., Chapman, S. J., Bhattal, A. S., Disney, M. J., Pongracic, H., and Turner, J. A. (1995).  Binary star formation: accretion-induced rotational fragmentation.  _MNRAS_ 277, 727–746.  10.1093/mnras/277.2.727
  * Wibking and Krumholz (2022) Wibking, B. D. and Krumholz, M. R. (2022).  QUOKKA: a code for two-moment AMR radiation hydrodynamics on GPUs.  _MNRAS_ 512, 1430–1449.  10.1093/mnras/stac439
  * Wise and Abel (2011) Wise, J. H. and Abel, T. (2011).  ENZO+MORAY: radiation hydrodynamics adaptive mesh refinement simulations with adaptive ray tracing.  _MNRAS_ 414, 3458–3491.  10.1111/j.1365-2966.2011.18646.x
  * Wolfire and Cassinelli (1986) Wolfire, M. G. and Cassinelli, J. P. (1986).  The Temperature Structure in Accretion Flows onto Massive Protostars.  _ApJ_ 310, 207.  10.1086/164676
  * Wünsch et al. (2021) Wünsch, R., Walch, S., Dinnbier, F., Seifried, D., Haid, S., Klepitko, A., et al. (2021).  Tree-based solvers for adaptive mesh refinement code FLASH - II: radiation transport module TreeRay.  _MNRAS_ 505, 3730–3754.  10.1093/mnras/stab1482
  * Wünsch et al. (2018) Wünsch, R., Walch, S., Dinnbier, F., and Whitworth, A. (2018).  Tree-based solvers for adaptive mesh refinement code FLASH - I: gravity and optical depths.  _MNRAS_ 475, 3393–3418.  10.1093/mnras/sty015
  * Zhang et al. (2011) Zhang, W., Howell, L., Almgren, A., Burrows, A., and Bell, J. (2011).  CASTRO: A New Compressible Astrophysical Solver. II. Gray Radiation Hydrodynamics.  _ApJS_ 196, 20.  10.1088/0067-0049/196/2/20

## Tables