# Numerical Methods for Simulating Star Formation

Romain Teyssier 1,∗ and Benoît Commerçon 2

###### Abstract

We review the numerical techniques for ideal and non-ideal magneto-
hydrodynamics (MHD) used in the context of star formation simulations. We
outline the specific challenges offered by modeling star forming environments,
which are dominated by supersonic and super-Alfvénic turbulence in a
radiative, self-gravitating fluid. These conditions are rather unique in
physics and engineering and pose particularly severe restrictions on the
robustness and accuracy of numerical codes. One striking aspect is the
formation of collapsing fluid elements leading to the formation of
singularities that represent point-like objects, namely the proto-stars.
Although a few studies have attempted to resolve the formation of the first
and second Larson’s cores, resolution limitations force us to use sink
particle techniques, with sub-grid models to compute the accretion rates of
mass, momentum and energy, as well as their ejection rate due to radiation and
jets from the proto-stars.

We discuss the most popular discretisation techniques used in the community,
namely smoothed particle hydrodynamics, finite difference and finite volume
methods, stressing the importance to maintain a divergence-free magnetic
field. We discuss how to estimate the truncation error of a given numerical
scheme, and its importance in setting the magnitude of the numerical
diffusion. This can have a strong impact on the outcome of these MHD
simulations, where both viscosity and resistivity are implemented at the grid
scale. We then present various numerical techniques to model non-ideal MHD
effects, such as Ohmic and ambipolar diffusion, as well as the Hall effect.
These important physical ingredients are posing strong challenges in term of
resolution and time stepping. For the latter, several strategies are discussed
to overcome the limitations due to prohibitively small time steps.

An important aspect of star formation simulations is the radiation field. We
discuss the current state-of-the-art, with a variety of techniques offering
pros and cons in different conditions. Finally, we present more advanced
strategies to mitigate the adverse effect of finite numerical resolution,
which are very popular in the context of supersonic, self-gravitating fluids,
namely adaptive mesh refinement, moving meshes, Smoothed Particle
Hydrodynamics and high-order methods. Advances in these three directions are
likely to trigger immense progress in the future of our field. We then
illustrate the different aspects of this review by presenting recent results
on supersonic MHD turbulence and magnetised collapse calculations.

##  1 Keywords:

star formation, numerical techniques, MHD: ideal, MHD: non-ideal,
astrophysical fluid dynamics, radiation fields, sink particles

##  2 Introduction

Star formation is one of the main unsolved problems in astrophysics. Although
our view of this fundamental process, at the nexus of galaxy formation,
planetary science and stellar evolution, has considerably changed over the
past decades, thanks to new observations and theoretical progress, many dark
corners remain to be explored. One of the reasons why the true origin of stars
still eludes us is the highly non-linear nature of the governing equations,
describing self-gravitating, compressible, magnetized dust and gas fluids
interacting with radiation. In this review, we present these main governing
equations, focusing on ideal magneto-hydrodynamics, radiation hydrodynamics,
non-ideal effects and sink particle techniques. We describe the most popular
discretisation schemes, focusing on possible sources of errors and their
consequences on clouding our conclusions. We finally give a short description
of the landscape in term of star formation simulations, with different set-up
strategies addressing particular scales in this fundamentally multi-scale
problem.

Star formation is believed to be the consequence of the collapse of mildly
supersonic or transonic to subsonic molecular cores emerging out of supersonic
turbulent flows in the interstellar medium (Mac Low and Klessen, 2004; McKee
and Ostriker, 2007). The source of turbulence is probably to be found on large
scales, as a consequence of large galactic shearing or colliding flows, but
also on small scales, because of various sources of stellar feedback (McKee,
1989; Federrath et al., 2017). In this context, gravitational or thermal
instabilities lead to the formation of dense gas clumps that undergo a more or
less violent gravitational collapse, leading to the formation of a proto-star
surrounded by a proto-stellar disk. Describing these processes using only
simple analytical methods is almost impossible. Moreover, the traditional
engineering methods in Computational Fluid Dynamics (CFD) are usually not
robust enough to sustain highly supersonic and super-Alfvénic flows. Self-
gravity, magnetic fields and radiation fields, taken together, define a very
unique system of equations that has no equivalent in the industry. This has
led astrophysicists to develop their own methods, largely inspired by
traditional methods designed in the applied mathematics community, but
significantly adapted to the specific problem at hand. In this context,
Smoothed Particle Hydrodynamics (SPH) and Adaptive Mesh Refinement (AMR)
techniques turned out to be particularly suitable to the problem of star
formation. Both techniques have their own pros and cons, and comparing the two
allows us to assess the robustness of our numerical results. New techniques
have also been developed, that are specific to star formation, such as sink
particles, a commonly adopted recipe to replace an otherwise collapsing fluid
element by a collision-less particle, saving computational times and
increasing the realism of the simulation.

In this review, we pay attention to the description of the equations, without
necessarily discussing their physical foundations, such as the ideal MHD limit
or the non-ideal diffusion processes. We describe the various numerical
techniques, from low-order schemes to modern high-order methods, as well as
from non-zero divergence schemes to exact divergence-free methods, etc. We
refer to the corresponding literature, including all references that are
relevant to the historical description of the discipline and that give a fair
snapshot of the present state of the field. We apologize in advance for not
having included all possible references on the topic.

##  3 Ideal MHD: numerical techniques

Developing new computational methods for solving the ideal MHD equations has
generated a lot of interest within the applied mathematics and computational
physics communities. Quite naturally, because of their success in solving the
Euler equations, grid-based methods with flux upwinding, also known as
Godunov’s method, were applied to the MHD equations in the framework of
finite-volume discretisation. In parallel, Smoothed Particle Hydrodynamics
(SPH) generated a lot of interest in the astrophysics community because of its
strict Galilean invariance. Both methods, however, quickly ran into
difficulties trying to maintain the divergence-free property of the MHD
equations.

###  3.1 The ideal MHD equations

Before we review the various improvements of MHD numerical solvers over the
past decades, we briefly recall the ideal MHD equations, shown here in
conservative form. We have first the mass conservation equation,

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho{\bf v}\right)=0,$ (1)

where $\rho$ is the mass density and ${\bf v}$ is the fluid velocity vector.
We also have the momentum conservation equation

$\frac{\partial}{\partial t}\left(\rho{\bf v}\right)+\nabla\cdot\left(\rho{\bf
v}\otimes{\bf v}+P_{\rm tot}{\mathbb{I}}-\frac{1}{4\pi}{\bf B}\otimes{\bf
B}\right)=\rho{\bf g},$ (2)

where ${\bf B}$ is the magnetic field vector and $P_{\rm tot}$ is the total
pressure, defined as the sum of the thermal pressure and the magnetic pressure

$P_{\rm tot}=P+\frac{1}{8\pi}B^{2}.$ (3)

Note that we work here in cgs units, hence the presence of the $4\pi$ term in
these equations. We have also included the gravitational acceleration vector
${\bf g}$ as a source term on the right-hand side of the momentum conservation
equation. Finally, we have the total energy conservation equation

$\frac{\partial E_{\rm tot}}{\partial t}+\nabla\cdot\left(\left(E_{\rm
tot}+P_{\rm tot}\right){\bf v}-\frac{1}{4\pi}{\bf B}({\bf B}\cdot{\bf
v})\right)=\rho{\bf g}\cdot{\bf v},$ (4)

where the total fluid energy is the sum of the kinetic energy, the internal
energy and magnetic energy

$E_{\rm tot}=\frac{1}{2}\rho v^{2}+e+\frac{1}{8\pi}B^{2}.$ (5)

In order to close the system, we need the fluid equation of state, usually
given be the ideal gas equation of state

$P=(\gamma-1)e,$ (6)

and the induction equation for the evolution of the magnetic field

$\frac{\partial{\bf B}}{\partial t}=-\nabla\times{\bf
E}{\rm~{}~{}~{}with~{}~{}~{}}{\bf E}=-{\bf v}\times{\bf B},$ (7)

where we have introduced the electric field ${\bf E}$ (also known as the
electromotive force, emf). We need to add to these equations the most
important property of the magnetic field, namely

${\nabla}\cdot{\bf B}=0,$ (8)

also known as the solenoidal constraint or the divergence-free condition.
Since we consider now in this section ideal MHD conditions, we have no
dissipation terms in the induction equation, as well as in the fluid
equations.

###  3.2 Some important mathematical properties

The fluid equations appear to be in conservative form, which naturally calls
for a finite volume representation, which will ensure exact conservation of
mass, momentum and total energy, by construction, owing to the divergence
theorem

$\int_{V}\nabla\cdot{\bf F}{\rm d}V=\int_{S}{\bf F}\cdot{\bf n}{\rm d}S,$ (9)

where the vector ${\bf F}$ can be the flux of mass, momentum or energy.
Designing a numerical scheme boils down to computing the flux through the
surface of the volume elements. On the other hand, the induction equation
naturally calls for a finite surface representation, owing to the Stoke’s
theorem

$\frac{\partial}{\partial t}\int_{S}{\bf B}\cdot{\bf n}{\rm
d}S=\int_{S}\nabla\times{\bf E}\cdot{\bf n}{\rm d}S=\int_{L}{\bf
E}\cdot{\bf{\rm d}l}.$ (10)

Similarly, the divergence-free condition, written in integral form as

$\int_{V}\nabla\cdot{\bf B}{\rm d}V=\int_{S}{\bf B}\cdot{\bf n}{\rm d}S,$ (11)

also calls for defining the magnetic field as a surface averaged quantity.
This naturalness argument, together with the fact that the divergence-free
condition can be maintained to machine precision accuracy, has led to the
design of the constrained transport scheme (Evans and Hawley, 1988). This very
popular method for grid-based techniques comes however with a price: the flow
variables are not all co-located at cell centers, like volume-weighted fluid
variables, but also at face centers, for the magnetic field components, and at
edge centers for the electric field components.

A very important mathematical property of the magnetic field that emerges from
the divergence-free condition is that the normal component of the field should
be continuous across cell faces. The $x$-component of the field, for example,
can be discontinuous in the $y$ and $z$-directions, but has to vary smoothly
in the $x$-direction. This property also naturally holds in the constrained
transport method. It can be written as

$\int_{S}\left[{\bf B}\cdot{\bf n}\right]{\rm d}S=0,$ (12)

where the operator $[A]=A^{+}-A^{-}$ denotes the jump of a quantity $A$ across
the surface element.

###  3.3 Preserving the divergence-free condition

Historically, one of the first ideal MHD, general purpose codes, is ZEUS-2D
developed specifically for astrophysics by Stone and Norman (Stone et al.,
1992a). It is a finite-difference code using constrained transport and
artificial viscosity to handle shocks, The continuity and divergence-free
conditions are therefore satisfied by construction, but since artificial
viscosity is used to handle shocks, instabilities can occur in fast
rarefaction waves (Falle, 2002).

The other popular strategy for grid-based methods is to maintain all MHD
variables as volume-averaged quantities, allowing for discontinuities across
all cell faces. In the late 90s, a series of papers presented finite-volume
MHD codes with proper upwinding of numerical fluxes using Riemann solvers (Dai
and Woodward, 1994; Ryu et al., 1995; Tóth, 1996; Balsara, 1998; Lee and
Deane, 2009). These methods are very powerful, because they are strictly
conservative and because they satisfy the so-called entropy condition, meaning
the entropy can only increase, a key property to maintain the stability of the
numerical solution. These finite-volume codes all considered the magnetic
field as a volume-averaged, piecewise-constant, cell-centered quantity, which
violates the continuity condition of the normal component of the field. The
resulting schemes are therefore not necessarily divergence-free anymore. In
order to maintain the divergence as small as possible, an additional step is
required that modifies the magnetic field components and decreases or
nullifies the divergence: this operation is called divergence cleaning. In a
seminal paper, Tóth (2000) has compared various schemes and showed that they
all passed with mixed success a variety of MHD test problems. We will now
review these grid-based methods that are using divergence cleaning.

The first method we discuss here is the projection scheme, introduced by
Brackbill and Barnes (1980). The idea is to measure the spurious divergence
after the main update of the MHD variables, and solve for a Poisson equation
defined as

$\Delta\phi=\nabla\cdot{\bf B}^{*}.$ (13)

where $\phi$ is a scalar potential. The new magnetic field is then corrected
as

${\bf B}^{n+1}={\bf B}^{*}-\nabla\phi,$ (14)

so that by construction the divergence is now zero. This method works very
well in many applications. It suffers however from two main issues: first, it
is quite expensive as it requires to solve for a Poisson equation. Another
consequence is that the correction is non-local: very localized divergence
errors can be instantaneously transported across the grid to enforce the
solenoidality condition. Second, the correction process modifies the magnetic
field, without modifying the total energy. As a consequence, the resulting
temperature can be modified, and the entropy condition might be violated. An
easy fix is to remove the magnetic energy before the correction step and add
the new magnetic energy after the correction, but the resulting scheme is not
conservative anymore.

The second popular method is the so-called 8-waves formulation or Powell’s
method (Powell et al., 1999). The idea is to write more general ideal MHD
equations allowing for the presence of magnetic monopoles. This results in the
following non-conservative form

$\displaystyle\frac{\partial}{\partial t}\left(\rho{\bf
v}\right)+\nabla\cdot\left(\rho{\bf v}\otimes{\bf v}+P_{\rm
tot}{\mathbb{I}}-\frac{1}{4\pi}{\bf B}\otimes{\bf B}\right)$
$\displaystyle=-\left(\nabla\cdot{\bf B}\right){\bf B}+\rho{\bf g},$ (15)
$\displaystyle\frac{\partial E_{\rm tot}}{\partial
t}+\nabla\cdot\left(\left(E_{\rm tot}+P_{\rm tot}\right){\bf
v}-\frac{1}{4\pi}{\bf B}({\bf B}\cdot{\bf v})\right)$
$\displaystyle=-\left(\nabla\cdot{\bf B}\right){\bf B}\cdot{\bf v}+\rho{\bf
g}\cdot{\bf v},$ (16) $\displaystyle\frac{\partial{\bf B}}{\partial
t}-\nabla\times\left({\bf v}\times{\bf B}\right)$
$\displaystyle=-\left(\nabla\cdot{\bf B}\right){\bf v}.$ (17)

This method proved very useful and robust for many applications. It is still
widely used in astrophysics nowadays. The success of the method originates
from the property that the spurious $\nabla\cdot{\bf B}$ is advected away with
the flow, using the so-called eighth MHD wave, so that divergence errors do
not accumulate. There are however two problems: First, in stagnation
points111Stagnation points are regions where the flow is brought to rest in
the frame of the system under study. More generally, even if the velocity does
not vanish, stagnation points are regions in the flow where the streamlines
are converging and the flow becomes compressive., the divergence errors will
accumulate because the flow trajectories are converging. Second, the scheme is
not strictly conservative. Shock waves lead to solutions that do not converge
to the correct conservation laws anymore.

A very elegant solution to the first problem was proposed in Dedner et al.
(2002) using the so-called hyperbolic-parabolic divergence cleaning technique,
also known as a Generalized Lagrange Multiplier (GLM) formulation of the ideal
MHD equations, in short, Dedner’s scheme. The idea is to add a ninth wave to
the previous Powell’s modified MHD equations, introducing the cleaning field
$\psi$ that satisfies

$\frac{\partial\psi}{\partial t}+{\bf
v}\cdot\nabla\psi+c_{h}^{2}\nabla\cdot{\bf B}=-\frac{\psi}{\tau},$ (18)

and is used as a source term in the induction equation.

$\frac{\partial{\bf B}}{\partial t}-\nabla\times\left({\bf v}\times{\bf
B}\right)=-\nabla\psi.$ (19)

In the first equation, $c_{h}$ is the divergence wave speed and $\tau$ is the
divergence damping time scale. The latter is chosen equal to (or larger than)
the fast magnetosonic wave speed, while the former is usually equal to (or
larger than) the fast magnetosonic cell crossing time. At first sight, this
new nine waves scheme can be seen as a combination of advection and damping of
divergence errors, thus a clever combination of the projection method and the
Powell’s scheme. As shown recently by Derigs et al. (2017), it is in fact much
more than that: this new divergence field $\psi$ allows to restore the entropy
condition, as this field can be interpreted as a divergence cleaning energy
which stores temporarily the magnetic energy lost during the damping step. It
is also conservative in a general sense, but still violates the Rankine-
Hugoniot shock relations locally.

In parallel, accurate and stable MHD solvers have been the topics of many
studies in the SPH framework. Probably because of its strict Galilean
invariance, early SPH MHD solvers were quite oscillatory (see e.g. Dolag et
al., 1999). Truncation errors associated to the non-zero divergence cannot be
damped by numerical diffusion as they are advected away, as in the case of
grid-based codes. Many regularization techniques have been proposed to provide
stable SPH solvers for the ideal MHD equations: using the vector potential
(e.g. Kotarba et al., 2009) or using artificial diffusivity (Price and
Monaghan, 2005; Dolag and Stasyszyn, 2009). Interestingly, the work of Price
and Monaghan (2005) was based on exploring various divergence cleaning methods
used in grid-based codes for SPH. The authors concluded that this works in
most cases, but in difficult cases, such as supersonic turbulent flows, this
can lead to spurious effects, such as the violation of the entropy condition
or wrong shock relations. More worrysome, these authors noticed that
divergence cleaning can lead to an increase of the local magnetic energy, a
price to pay to redistribute the truncation errors in the divergence. This
results in spurious dynamo effects. Recently, however, Tricco and Price (2012)
revisited the Dedner’s scheme for SPH and found a formulation that guarantees
that divergence cleaning leads to a decrease of the magnetic energy and an
increase of the entropy, in the spirit of Derigs et al. (2017) for grid-based
codes.

In light of these rather complex developments of divergence cleaning methods,
the constrained transport (CT) approach we have introduced already seems quite
appealing. Note however that this approach requires to have well defined cell-
interfaces, which is not necessarily the case for SPH or the recently
developed moving-mesh AREPO code (Springel, 2010). The CT scheme, introduced
for astrophysical fluid flows by Evans and Hawley (1988) and used in the
ZEUS-2D code (Stone et al., 1992a), features the nice property that the
divergence of the magnetic field, defined in integral form as the net magnetic
flux across the 6 faces of each cell, will always remain constant, So if it is
initially zero, it will remain equal to zero within numerical round-off
errors. This however requires to define the electric field on the cell edges,
and this is the main difficulty of this approach (see Figure 1). Indeed, in
the Godunov methodology, fluxes, defined at cell faces, are computed as the
solution to the Riemann problem of the two neighboring cells meeting at their
common cell interface. For the electric field, defined at cell edges, 4
neighboring cells meet and in order to properly upwind the electric field, we
need to solve a two-dimensional Riemann problem, a rather daunting task (see
Teyssier et al., 2006, for a discussion). Several solutions have been found to
design 2D MHD Riemann solvers (Londrillo and del Zanna, 2004; Balsara et al.,
2014), and this has been the main characteristic of several simulations codes
with proper upwinding of both fluxes and electric fields, using the CT scheme
within the Godunov methodology (Londrillo and del Zanna, 2004; Fromang et al.,
2006; Stone et al., 2008; Li et al., 2012; Mocz et al., 2014a).

Figure 1: Schematic showing the geometry of a Cartesian cell in the
Constrained Transport approach. The finite volume cell is shown as a grey
cube. It is labelled $i,j,k$. The magnetic field components are defined
perpendicular to the faces of the cube. The are shown in red, only in the
rightermost face in all three directions. The electric field components are
defined on the cell edges. They are shown in blue and only in one face for
sake of simplicity.

###  3.4 Minimizing numerical errors: higher order versus mesh refinement

When performing self-gravitating magnetized fluid simulations, it is of
primary importance to quantify and understand numerical errors. These errors
are often called truncation errors because they arise as leading order terms
in the Taylor expansion of the adopted discrete numerical method. Usually,
this Taylor expansion leads to the so-called modified equation, which encodes
how the original ideal MHD equations have been transformed into a new set of
equations with additional terms coming from the truncation errors. For first-
order schemes, these terms are identical to a diffusion process, and are
therefore called numerical viscosity or numerical diffusivity. In SPH, quite
often, these terms are added explicitly to the equations as artificial
viscosity and artificial diffusivity, while for grid-based Godunov solvers,
these terms are only implicitly present, through this Taylor expansion of the
discrete scheme. In both cases, however, these diffusion processes control how
shock heating and magnetic reconnection proceed in the solution. They play a
fundamental role in preserving the entropy condition, and in regulating the
flow close to the grid scale. Unfortunately, many complex MHD processes, such
as the small scale dynamo (Brandenburg and Subramanian, 2005) or the magneto-
rotational instability (Balbus and Hawley, 1991) depend crucially on the so-
called Prandtl number (Fromang et al., 2007; Federrath et al., 2014), which is
the ratio of the real magnetic diffusivity to the viscosity. In most case,
this ratio is always close to 1 if one uses the numerical Prandtl number
(Fromang et al., 2007; Federrath et al., 2011a), while in nature, it can vary
widely. It is therefore crucial to adopt in some cases explicit viscosity and
diffusivity coefficients and model the flow including these non-ideal
processes (Fromang et al., 2007; Federrath, 2016).

As explained in the next section, in order to model these non-ideal effects,
it is crucial to control and minimize the numerical diffusion as much as
possible. There are two possible strategies to achieve this: refine the grid
or increase the order of accuracy of the method. The first approach leads to
the so-called Adaptive Mesh Refinement method, a very popular and successful
technique in the context of star and galaxy formation. Popular AMR codes are
available to the star formation community, such as RAMSES (Teyssier, 2002;
Fromang et al., 2006), ORION (Krumholz et al., 2007a; Klein, 1999), ENZO
(Bryan et al., 2014) or FLASH (Fryxell et al., 2000). Other codes, that used
to be only unigrid, now propose an adapted grid or an adaptive grid extension,
such as PLUTO (Mignone et al., 2007) and ATHENA (Stone et al., 2008). In all
these codes, cells are adaptively refined according to various refinement
criteria. In the context of star formation, the most popular approach is to
always resolve the local Jeans length with 4 cells or more, the so-called
Truelove criterion (Truelove et al., 1997)

${\rm If~{}}\lambda_{\rm J}=\frac{c_{s}}{\sqrt{4\pi G\rho}}<4\Delta
x_{\ell}~{}~{}~{}{\rm then~{}refine~{}to~{}level~{}}{\ell+1}.$ (20)

This corresponds to a level-dependent density threshold that triggers new
refinements. SPH methods are Lagrangian in nature, so they cannot refine as
much as grid-based codes. Usually, much more SPH particles are needed in the
initial conditions to reach a certain target Jeans mass, corresponding to the
maximum resolution level of the corresponding AMR simulation. Particle
splitting is an interesting alternative to classical SPH but is still under
development (Kitsionas and Whitworth, 2002; Chiaki and Yoshida, 2015), the
difficulty being to handle sharp transitions in particle mass and its
interaction with the smoothing kernel. Note that if similar resolution
requirements are met, AMR and SPH methods largely agree on the quantitative
predictions on how the collapse proceeds (Commerçon et al., 2008). For
magnetized flows, the Jeans length-based refinement strategy has to be more
conservative, of the order of 30 cells per Jeans length, in order to capture
properly the magnetic field amplification in collapsing cores (Sur et al.,
2010; Federrath et al., 2011b; Turk et al., 2012).

A difficulty arises when one uses AMR for Constrained Transport. In this case,
it is mandatory to be able to interpolate the magnetic field from the coarser
level to the finer level and still satisfy the divergence-free condition. For
this, divergence preserving interpolation schemes have been developed
(Balsara, 2001; Tóth and Roe, 2002) and play an important role in the
viability of the CT approach in the context of AMR.

To reduce truncation errors, the other option is to use higher order schemes.
The solution inside each cell is not piecewise constant as in the traditional
first-order Godunov method, but it is reconstructed using high-order
polynomials of degree $p$ as

$\rho(x)=\sum_{i=0,p}\alpha_{i}\psi_{i}(x).$ (21)

The $\psi_{i}(x)$ are usually an orthonormal basis of polynomials of degree at
most $p$. The coefficients $\alpha_{i}$ are computed using two different
philosophies. The first option, the WENO approach computes the coefficients
$\alpha_{i}$ at each time step using neighboring cells. The higher the
polynomial degree, the more neighbors must be included in the stencil of the
method (Jiang and Wu, 1999; Boscheri and Dumbser, 2014). The second option,
the Discontinuous Galerkin (DG) approach, considers that the $\alpha_{i}$ are
new variables, defined for each MHD variables, and updated using a projection
of the ideal MHD equation onto the polynomial basis (Li and Xu, 2012;
Klingenberg et al., 2017; Mocz et al., 2014b; Guillet et al., 2018). In the
WENO case, one needs to access neighbouring but possibly distant cells to
compute the $\alpha_{i}$, while in the DG case, one needs to store permanently
the $\alpha_{i}$ and solve for their evolution at each time step. Note that
other high-order methods are also being developed that do not strictly
correspond to neither WENO nor DG (Felker and Stone, 2018).

In the context of high-order methods, the divergence free condition is
particularly challenging. One can either implement high-order version of one
of the above mentioned divergence cleaning techniques (Derigs et al., 2016),
or one can try to preserve the divergence-free constraint within the method.
The second approach could be seen as a generalisation of the CT scheme to
higher order. The key ingredients are: 1- the use of a divergence-free
polynomial basis for the magnetic field (Li and Shu, 2005; Guillet et al.,
2018), 2- enforcing the continuity of the normal field across cells boundaries
(Li and Xu, 2012). Although some very promising solutions have been found
recently, this is still a very active field of research. Interestingly, the
traditional CT scheme with face-centred magnetic field variables can be re-
interpreted as a DG scheme where each magnetic field component is piecewise-
linear and continuous in the normal direction and piecewise constant in the
transverse direction, so that two cell-centred coefficients are required
(instead of one) for each field component.

##  4 Non-ideal MHD: numerical techniques

###  4.1 Equations and basic concepts

Star formation takes place in molecular clouds, which are made of a mixture of
dust and gas, implying dust and gas collisions, and both constituents are far
from being fully ionised. In the previous section, we presented the work done
in the ideal MHD framework, which do not appear to be well suited to the
ionisation state in collapsing cores, in particular at the onset of disk
formation. Recent works have emphasized the imperfect coupling of the dust and
gas mixture with the magnetic fields at the transition between the envelop and
the disk in collapsing cores (see the review by Wurster and Li (2018) for the
work done in the context of protostellar disk formation) and a lot of effort
has been devoted over the past ten years to include the so called non-ideal
effects: the ambipolar diffusion, the Ohmic diffusion, and the Hall effect.
The ambipolar diffusion is the common name to describe the interaction between
neutrals and charged particles. It can be seen as a friction term, it enables
the neutral field to respond to the magnetic forces, via collisions with
charged particles. The Ohmic diffusion results from the collision of electrons
with the neutrals. Last, the Hall effect is due to the drift between the
positively and negatively charged species. As shown in Marchand et al. (2016)
and Tsukamoto et al. (2015a), all these three terms can be dominant over the
others at different scales within the envelop of collapsing dense cores. For a
classical dust size distribution, Mathis et al. (1977) and Marchand et al.
(2016) find that ampipolar diffusion and the Hall effect dominate at densities
$<10^{12}$ cm-3 and that Ohmic diffusion is the stronger resistive effect at
higher densities.

The exact scale and density at which these resistive effects become dominant
over the other dynamical processes depend on the chemistry, the ionisation
intensity and the dust grain size distribution (Zhao et al., 2016; Dzyurkevich
et al., 2017; Zhao et al., 2018, 2018). Hennebelle et al. (2016) have shown
that ambipolar diffusion regulates the flow dynamics over the other dynamical
processes (induction, rotation and free-fall) at scale of a few 10s AU, which
sets the initial size of protostellar disks. In addition, Masson et al. (2016)
have shown that scales of a few 10s AU exhibit magnetic Reynolds numbers less
than unity (see also Tomida et al. (2013, 2015); Vaytet et al. (2018) for
studies including the Ohmic diffusion). Recently, Koga et al. (2019) performed
a similar analysis considering the Hall effect only and found coherent results
as for ambipolar diffusion.

Before we describe the numerical implementation for the three aforementioned
resistive effects, let us recall the main equations and define the necessary
quantities. Ideally, the system should account for the different behaviour of
the neutrals, as well as negatively and positively charged particles. Among
them the different constituents that participate to the ionisation balance, we
have neutral (molecule, atoms, dust grains), ions (molecular and atomic),
electrons, and charged dust grains. The latter can be positively or negatively
charged, and can hold multiple charges (Draine and Sutin, 1987). Most current
works follow a one-fluid approximation to describe the evolution of this
complex mixture, where the Ohm’s law accounts for the non-ideal effects with
some assumptions. In the following, we focus on the work done in the single
and two-fluid approximations. We review briefly the work done towards a multi-
fluid treatment of the charged and neutral particles in Sec. 4.6.

In the low ionisation limit, each kind of charged particle is scarce, so that
the pressure terms, the gravity, and the inertia term can be neglected in
their momentum evolution equation in comparison to the Lorentz force and the
frictional force exerted by the neutrals. The generalised Ohm’s law reads

${\bf E_{NIMHD}}={\bf v}\times{\bf B}-\eta_{\Omega}{\bf
J}-\frac{\eta_{\mathrm{H}}}{||{\bf B}||}{\bf J}\times{\bf
B}+\frac{\eta_{\mathrm{AD}}}{||{\bf B}||^{2}}{\bf J}\times{\bf B}\times{\bf
B},$ (22)

where ${\bf J}=\nabla\times{\bf B}$ is the current, $\eta_{\Omega}$,
$\eta_{\mathrm{H}}$, and $\eta_{\mathrm{AD}}$ are the Ohmic, Hall, and
ambipolar resistivities, respectively (in units of cm2 s-1), and ${\bf v}$ is
the velocity of the neutrals. The notation $||{\bf B}||$ represents the norm
of the magnetic field vector ${\bf B}$. The electric field is then replaced in
Eq. (7). It is worth noticing that all the resistive terms lead to parabolic
partial differential equations to discretise. We note that only the Ohmic term
can lead to reconnection here. It can be rewritten as a Laplacian operator
$\eta_{\Omega}\Delta{\bf J}$, whereas the two others cannot. The ambipolar
diffusion and the Hall effect can be assimilated to drift velocities of the
magnetic fields with respect to the neutral speed ${\bf v}$. The total energy
equation has to be modified accordingly, to account for the heating resulting
from the friction terms, i.e., the ambipolar diffusion and the Ohmic diffusion

$\frac{\partial E_{\rm tot}}{\partial t}+\nabla\cdot\left(\left(E_{\rm
tot}+P_{\rm tot}\right){\bf v}-\frac{1}{4\pi}{\bf B}({\bf B}\cdot{\bf
v})\right)=\rho{\bf g}\cdot{\bf v}+\eta_{\Omega}||{\bf
J}||^{2}+\eta_{\mathrm{AD}}\frac{||{\bf J}\times{\bf B}||^{2}}{||{\bf
B}||^{2}}.$ (23)

We refer readers to Shu et al. (1987), Nakano et al. (2002) and Balbus (2009)
for more details on the derivation of the generalised Ohm’s law in the one-
fluid approximation. In the context of star formation, where a huge dynamical
range has to be considered, the estimate of the resistivities is challenging.
It requires to know the ionisation state of the gas and dust mixture. The
coupled chemistry of gas and dust has to be considered, which involves a huge
variety of physical processes: gas phase and gas-dust chemistry, cosmic ray
ionisation, UV radiation, (in-)elastic collisions, grain growth. This domain
remains the subject of intense research works that we do not detail in this
review. We refer readers to Sec. 4.0.1 in the review by Wurster and Li (2018)
for a scan of the works that tackle the resistivity calculations.

Ambipolar diffusion and Ohmic diffusion were the first to be considered in
star formation applications. Indeed, these two terms do not introduce new MHD
waves in the system (they do not change the spectral properties of the
hyberpoblic system), so that they do not require heavy modifications of the
MHD solver. The induction equation integration is most of the time splitted in
two steps: first the ideal MHD and second the resistive terms. If an energy
equation is considered, the heating terms due to resistive effect in Eq. 23
are integrated as source terms in most numerical implementations.

###  4.2 Ohmic diffusion

The Ohmic diffusion is perhaps the simplest resistive term to introduce
because of its diffusive nature, similar to the introduction of artificial
viscosity and artificial resistivity to prevent any numerical artefact due to
numerical diffusion. Assuming a constant resistivity, the Ohmic term can be
rewritten as a Laplacian (which preserves the solenoidal constraint), leading
to the following stability condition

$\Delta t_{\Omega}=\frac{\Delta x^{2}}{\eta_{\Omega}}.$ (24)

It corresponds to the stability condition required in numerical schemes
integrating parabolic equations such as the heat equation with an explicit
scheme.

Numerous implementations for the Ohmic diffusion have been developed in the
past ten years, in all different kind of codes. Masson et al. (2012) present a
fully explicit implementation in the AMR code RAMSES, where they update the
emf to account for the resistive terms and then use the CT scheme as for ideal
MHD. Because of the restrictive stability condition for the resistive term
which can be smaller that the MHD one, some authors considered schemes that
enable to relax the time step constraint. For SPH, Bonafede et al. (2011)
implemented Ohmic diffusion with a constant resistivity in GADGET following
the method used for artificial dissipation (e.g., Price and Monaghan, 2005;
Dolag and Stasyszyn, 2009). Tomida et al. (2013) implemented an explicit
scheme in a nested-grid code following the implementation of Matsumoto (2011)
in the AMR code SFUMATO. They adopt the Super-Time-Stepping (STS) algorithm
proposed by Alexiades et al. (1996) to accelerate the integration of the
diffusion term and to relax the stringent explicit CFL condition.

The STS consists in a discretisation as a classical explicit scheme, which
does not require to satisfy the stability in every step but after a series of
integration cycles which are made of $N_{\mathrm{STS}}$ sub-timesteps $\Delta
t_{\mathrm{STS}}$ which are based on Chebyshev polynomials to guarantee the
stability conditions over the super-timestep $\Delta
t_{\mathrm{tot}}=\sum_{1}^{N_{\mathrm{STS}}}\Delta t_{\mathrm{STS}}$. We also
note the work of Mignone et al. (2007); O’Sullivan and Downes (2006) where
Ohmic diffusion is also integrated using the STS technique. Similarly,
Tsukamoto et al. (2013) proposes a SPH discretisation of Ohmic dissipation
using STS. Wurster et al. (2016) and Price et al. (2017) (in the PHANTOM code)
use similar implementations on different SPH codes. Last, Marinacci et al.
(2018) propose two implementations in AREPO: one using Powell divergence
cleaning and another using CT. Both implementations can be integrated using an
explicit or implicit scheme. Almost all the implementations mentioned so far
for the Ohm diffusion are combined with ambipolar diffusion (see below).

Apart from resistive MHD, the STS technique is becoming increasingly popular
in the computational astrophysics community to accelerate parabolic
timestepping advancements in the context of radiation hydrodynamics with the
FLD (Commerçon et al., 2011) and anisotropic diffusion (Meyer et al., 2012;
Vaidya et al., 2017).

###  4.3 Ambipolar diffusion

Ambipolar diffusion is not as straightforward to include since the associated
term in the induction equation can not be rewritten as a Laplacian
(Brandenburg and Zweibel, 1994). For the diffuse ISM, the drift between ions
and neutrals is the effect of the ambipolar diffusion, that can be written as
a friction force

$\bf{F}_{in}=\rho_{i}\rho_{n}\gamma_{\mathrm{AD}}(\bf{v}-\bf{v}_{i}),$ (25)

where $\rho_{i}$ (resp. $\bf{v}_{i}$) and $\rho_{n}$ ($\bf{v}$) are the mass
density (resp. velocity) of ions and neutrals, and $\gamma_{\mathrm{AD}}$ is
the collisional coupling constant, with $\gamma_{\mathrm{AD}}\sim 3\times
10^{13}$ cm3 s-1 g-1 in ISM mixture (Draine et al., 1983). In the one-fluid
and strong coupling limit, the inertia of ions is neglected, as well as the
pressure and gravitational forces (the ionisation degree is very low). As a
consequence, the Lorentz force equals the drag force so that the ion drift
velocity is

$\bf{v}_{d}\equiv\bf{v}_{i}-\bf{v}=\frac{1}{\rho_{i}\rho_{n}\gamma_{\mathrm{AD}}}\left(\nabla\times\bf{B}\right)\times\bf{B}.$
(26)

The induction equation thus reads

$\frac{\partial\bf{B}}{\partial
t}=\nabla\times[\bf{v}_{i}\times\bf{B}]=\nabla\times\left[\bf{v}\times\bf{B}+\frac{1}{\rho_{i}\rho_{n}\gamma_{\mathrm{AD}}}\left(\nabla\times\bf{B}\right)\times\bf{B}\times\bf{B}\right].$
(27)

From Eq. 22, we have
$\eta_{\mathrm{AD}}=||\bf{B}||^{2}/\rho_{i}\rho_{n}\gamma_{\mathrm{AD}}$. In
the dense ISM, the heavy charged particles are ions and charged dust grains,
and the expression of the resistivity is more complex (Marchand et al., 2016).
In the two-fluid approximation, the ion inertia is not neglected and the
system of equation accounts for both the ions and the neutrals momenta
equation.

The first implementation of ambipolar diffusion in a grid code was presented
in Black and Scott (1982) using a two-fluid approximation. They integrated the
induction equation using the ion velocity and treated the collision between
neutrals and ions as a friction force in the momenta equations. Mac Low et al.
(1995) present a first 3D implementation of one-fluid ambipolar diffusion in
the ZEUS code using a fully explicit scheme as well as CT. They propose to
subcycle the magnetic fields evolution in the case when ambipolar diffusion
timestep is much shorter than the dynamical one. The first SPH implementation
can be found in Hosking and Whitworth (2004), where they used a two-fluid
approximation, but their scheme did not preserve the solenoidal constraint.
Duffin and Pudritz (2008) present a one-fluid implementation in the FLASH code
using Powell’s method (Powell et al., 1999). Masson et al. (2012) also
implemented a one-fluid scheme for ambipolar diffusion in RAMSES with the same
philosophy as for the Ohmic diffusion. A similar implementation using CT,
operator splitting, and STS in the ATHENA code can be found in Bai and Stone
(2011) and Chen and Ostriker (2014). Last, Wurster et al. (2014) propose a
one-fluid implementation for SPH codes, which they implemented with success in
PHANTOM.

Other one-fluid implementations in the strong coupling approximation can be
found in Padoan et al. (2000) in a 3D explicit, finite difference grid code,
Choi et al. (2009) in the MHD TVD code using a flux-interpolated CT scheme
(Ryu et al., 1998; Balsara and Spicer, 1999), Christie et al. (2017) in the
ENZO’s Dedner MHD solver. Oishi and Mac Low (2006) and Li et al. (2006)
presented two simultaneous but independent implementations in the ZEUS code of
a two-fluid solver in the heavy ion approximation.

For star formation applications, most of these methods are usually second-
order accurate and at best second-order accurate in time (with a few
exceptions though (e.g., Meyer et al., 2014)). Clearly, further work needs to
be done to improve the accuracy of resistive MHD solvers. Alongside, multi-
fluid approaches are also needed to account for more detailed physics (see
Sec. 4.6).

###  4.4 Hall effect

The last resistive effect that has been implemented in star formation models
is becoming increasingly popular. The Hall effect originates from the drift
between ions and electrons, or, more generally the positively and negatively
charged particles. The derivation of the generalised Ohm law (22) is not as
straightforward as in the case of ambipolar diffusion only. We shall consider
the momentum equations of ions and electrons, accounting for the Lorentz force
as well as collisions between ions, electrons, and neutrals. If we restrict
the resistive effects to the Hall effect only, the induction equation reads

$\frac{\partial\bf{B}}{\partial
t}=\nabla\times\left[\bf{v}\times\bf{B}-\frac{\eta_{\mathrm{H}}}{||{\bf
B}^{|}|}{\bf J}\times{\bf B}\right]=\nabla\times\left[(\bf{v}+{\bf
v}_{\mathrm{H}})\times\bf{B}\right],$ (28)

where we introduce the ”Hall speed” ${\bf
v}_{\mathrm{H}}=-\frac{\eta_{\mathrm{H}}}{||{\bf B}||}{\bf J}$. Contrary to
the two previous resistive terms, the Hall effect introduces new waves in the
MHD systems: the whistler waves. The dispersion relation for whistler waves
reads (Balbus and Terquem, 2001)

$\omega=\pm\frac{\eta_{\mathrm{H}}k^{2}}{2}+\sqrt{\left(\frac{\eta_{\mathrm{H}}k^{2}}{2}\right)^{2}+k^{2}c^{2}_{\mathrm{A}}}.$
(29)

The whistler waves are dispersive, the phase speed $c_{\mathrm{w}}=\omega/k$
depends on the wavenumber $k$ and tends to infinity as the wavenumber tends to
infinity. Since any discrete formulation cannot extend to infinity, a
numerical solver cannot follow the whistler waves with very high frequencies.
In practice, the whistler wave speed is chosen to be the one at the grid scale

$c_{\mathrm{w}}=\left|\frac{\eta_{\mathrm{H}}\pi}{2\Delta
x}\right|+\sqrt{\left(\frac{\eta_{\mathrm{H}}\pi}{2\Delta
x}\right)^{2}+c^{2}_{\mathrm{A}}}.$ (30)

The time step is then constrained as

$\Delta t\leq\frac{\Delta x}{c_{\mathrm{w}}}\approx\frac{\Delta
x^{2}}{\eta_{\mathrm{H}}}.$ (31)

The first way to implement the Hall effect is to use an operator split method,
similarly to the Ohm and ambipolar diffusion. A pioneer implementation is
presented in Sano and Stone (2002) in a second-order Godunov finite-difference
code (Sano et al., 1999), in the local shearing box framework. They use the CT
method to update the induction equation with emf computed with the Hall term.
The method has been implemented in 3D in ZEUS in Krasnopolsky et al. (2011).
Falle (2003) presents an implicit scheme for multi-fluid MHD which alleviates
the stringent time step conditions, where the time step goes to zero as the
strength of the Hall effect increases. O’Sullivan and Downes (2006) designed a
Hall Diffusion Scheme (HDS) in an explicit 3D code for multi-fluid MHD. They
split the Hall diffusion in two parts, where the first part uses the STS time
integration up to a critical value of the Hall resistivity (i.e. to satisfy
the stability condition), and the second part, the HDS, diffuses the excess
Hall resistivity with a different discretisation scheme. The balance between
the ambipolar diffusion and the Hall resistivity determines the excess Hall
resistivity. The HDS integrates the induction equation for the excess Hall
using a dimensional splitting which is explicit for the first dimension, i.e.
the $x$ component $B_{x}^{n+1}$ of the magnetic fields is given by $B_{y}^{n}$
and $B_{z}^{n}$ at time $n$, explicit-implicit for the second dimension,
$B_{y}^{n+1}$ is given by $B_{x}^{n+1}$ and $B_{z}^{n}$, and finally implicit
for third dimension, $B_{z}^{n+1}$ is given by $B_{x}^{n+1}$ and
$B_{y}^{n+1}$. Both STS and HDS can be subcycled to reach the time step limit
given by the hyperbolic system. The HDS scheme has been also implemented in
ATHENA by Bai (2014) using CT.

An alternative for finite volume methods is to include the Hall effect into
the conservative integration scheme. Tóth et al. (2008) proposed a first
implementation in the AMR code BATSRUS (Tóth et al., 2006). It uses block-
adaptive grids with both explicit and implicit time discretisation, and
various methods to control numerical errors in $\nabla\cdot{\bf B}$ (8-waves
scheme with diffusion, projection scheme). They achieve spatial second-order
convergence by using symmetric slope limiters (like monotonized central)
instead of asymmetric limiters (like minmod). In addition, they show that the
first-order Lax-Friedrich Riemann solver is inconsistent for the Hall MHD
equations. This comes from the fact that the whistler wave speed is
proportional to the grid spacing, which makes a one order loss in space in the
truncation errors of the scheme. The same applies for MUSCL schemes and
asymmetric slope limiters which are only first-order accurate. Lesur et al.
(2014) built upon Tóth et al. (2008) a Hall MHD solver in PLUTO for
protoplanetary disc evolution studies. They implement a HLL solver to estimate
the flux and then update the magnetic fields using CT. When written in a
conservative form, the Hall induction equation reads

$\frac{\partial{\bf B}}{\partial t}=\nabla\cdot\left[{\bf vB}-{\bf
Bv}-\frac{\eta_{\mathrm{H}}}{||{\bf B}||}({\bf JB}-{\bf BJ})\right].$ (32)

The current ${\bf J}=\nabla\times{\bf B}$ involves spatial derivatives of the
magnetic fields, so that the flux depends on the conserved quantities, and on
its derivative. This leads to a ill-defined Riemann problem. Lesur et al.
(2014) solved this issue by assuming that J is an external parameter, which is
computed at the cell interfaces using the predicted states. The left and right
fluxes take this unique value as an input, with the predicted B. The Godunov
fluxes are then computed using a whistler-modified HLL solver, which accounts
for the whistler waves speed (30) in the characteristics speeds.

Last, Marchand et al. (2018) recently extended the Lesur et al. (2014)
implementation to the CT scheme of RAMSES, which uses 2D Riemann solver to
integrate the induction equation. They designed a 2D whistler-modified HLL
solver to estimate the emf on cell border, assuming a uniform current in the
2D Riemann solver, again computed from the predicted states.

In SPH codes, the Hall effect has been implemented and used with success for
star formation applications in Wurster et al. (2016), Tsukamoto et al. (2017)
and Price et al. (2017). The standard method is based on an operator splitting
and either a full explicit or a STS scheme based on the implementation for
ambipolar diffusion (Wurster et al., 2014).

To date, only Krasnopolsky et al. (2011), Wurster et al. (2016), Tsukamoto et
al. (2017) and Marchand et al. (2018) have applied their methods for
protostellar collapse. While the SPH methods seem to lead to a relatively
accurate conservation of the angular momentum, the grid-based methods using CT
are suffering from severe non-conservation of the angular momentum. The origin
of this non-conservation is unclear. Krasnopolsky et al. (2011) invoked Alfvén
waves going out of the simulation volume, while Marchand et al. (2018) did not
find evidence of this. Instead, they propose that the non-conservation takes
place after the accretion shock formation. Shocks indeed generate strong
gradients and thus large Hall velocities. Further investigation is clearly
needed to cope with this fundamental issue.

###  4.5 Full resistive implementation

Currently, a handful of 3D codes benefit of a full implementation of the three
resistive terms: PLUTO (Lesur et al., 2014), ATHENA (Bai, 2014), RAMSES
(Masson et al., 2012; Marchand et al., 2018), the SPH code by Tsukamoto et al.
(2017), PHANTOM (Price et al., 2017), ZeusTW (Li et al., 2011), and GIZMO
(Hopkins, 2017). A very recent full implementation for solar physics can be
found in González-Morales et al. (2018) in the MANCHA3D grid code. It combines
STS for ambipolar and ohmic diffusion, and HDS for the Hall effect.

Currently, the ATHENA, PLUTO, RAMSES, and ZeusTW implementations rely on CT
algorithms for the resistive MHD equations integration. All other works use
divergence cleaning algorithms. Given the variety of the methods developed in
the literature, a quantitative comparison (in standard tests as well as in
star formation applications) of all these implementations would be welcome in
the coming years.

Last but not least, it is not clear which non-ideal process dominates over the
others. The physical conditions in star-forming regions are so wide that there
is room for regions where all three effects dominate. In addition, every
effect needs to be tested and quantified before any conclusions on the
relative importance of each of these effects (see the review by Wurster and
Li, 2018).

###  4.6 Multifluid approach

Multifluid approaches are designed to describe multiphase or multifluid flows.
In the context of star formation, different fluids or phases are at stack:
neutrals, ions, electrons, and dust grains. A few attempts have been made to
describe all or a part of these components using multi-fluid approach. In this
section, we consider only the approaches which do account for $N$ momentum
equations, where $N$ is the number of phases or fluids. For instance, a
popular two-fluid framework considers ions and neutrals and describes the
friction between these two fluids, i.e. the ambipolar diffusion. Toth (1994)
presented a first 2D implementation of two-fluid ions and neutrals dynamics,
which is based on a flux-corrected transport finite difference method. Inoue
and Inutsuka (2008) propose an unconditionally stable numerical method to
solve the coupling between two fluids ions-neutrals (frictional
forces/heatings, ionization, and recombination, (Draine, 1986) for an ISM gas
subject to the thermal instability. They split time integration as follows: 1/
the ideal hydrodynamical part for neutrals, 2/ the ideal MHD part for ions.
The first part is integrated using a second-order Godunov scheme. The ideal
MHD part is solved in two steps, similarly to what is done for one-fluid ideal
MHD solvers (e.g. Stone et al., 1992a): the magnetic pressure terms are solved
using a second-order Godunov method, and the magnetic tension term and the
induction equation are solved using the method of characteristics for Alfvèn
waves with the CT algorithm. Tilley and Balsara (2008) present a semi-implicit
method for two-fluid ion-neutral ambipolar drift in the RIEMANN code (Balsara,
2001, 2004). Their scheme is second-order accurate in space and time and uses
a mixed implicit-explicit formulation to deal with strong friction regimes.
Pinto et al. (2008) consider a system composed of three fluids: positively-
charged particles, negatively-charged particles, and neutrals. The charged
particles include ions, electrons, as well as charged dust grains, depending
on the main charge carrier. The self-consistent set of equations they designed
remains to be implemented in 3D codes.

Last but not least, a complete multi-fluid approach should include dust grains
dynamics on top of the neutrals, ions, and electrons. Falle (2003) presented a
scheme for multi-fluid hydrodynamics in the limit of small mass densities of
the charged particles, i.e., where the inertia of the charged particles are
neglected. This is a similar to the heavy ion approximation but for neutrals,
ions, electrons, and dust grains of various sizes. The HD equations are first
integrated for the neutrals using a second-order Godunov scheme, and then the
induction equation and the charged species velocities are updated using an
implicit scheme. This scheme captures the three non-ideal effects and can deal
with regime where the Hall effect dominates. We note that similar approaches
were introduced in Ciolek and Roberge (2002) and O’Sullivan and Downes (2006).

The most complete work on multi-fluid designed for star formation was
presented in Kunz and Mouschovias (2009). They considered a six-fluid set of
equations (neutrals, electrons, molecular and atomic ions, positively charged,
negatively charged, and neutral grains) and implemented it in a 2D version of
the radiation-MHD ZEUS-MP code (Hayes et al., 2006). They modified the MHD
solver to account for non-ideal effects (Ohmic dissipation and ambipolar
diffusion). The abundances of the 6 species are calculated using a reduced
chemical-equilibrium network. They applied their methods to protostellar
collapse in Kunz and Mouschovias (2010). Unfortunately, the CPU cost of such a
complete set of physics is prohibitive for 3D models, and no work has taken
over since this first application.

###  4.7 What is next? Dust and gas mixture dynamics

Currently, numerous studies report on the effect of non-ideal MHD in the star
formation process. The wide majority of these works use the single-fluid
approximation, assuming furthermore that the dust grains are perfectly coupled
to the gas via collisions. Nevertheless, dust grains are observed to cover a
wide range of sizes, from nanometer to micrometer in the ISM (Mathis et al.,
1977) and even millimeter in protostellar disks environments (Chiang et al.,
2012). Generally, it assumes that dust grains of different sizes coexist and
follow a power-law dust distribution. Grains react to changes in the gas
velocity via a drag force. For a given dust grain of mass $m_{\mathrm{g}}$,
the force is

${\bf F}_{\mathrm{drag}}=-\frac{m_{\mathrm{g}}}{t_{\mathrm{s}}}({\bf
v}_{\mathrm{g}}-{{\bf v}}),$ (33)

where ${\bf v}_{\mathrm{g}}$ is the dust grain velocity, ${\bf v}$ the gas
velocity, and $t_{\mathrm{s}}$ the stopping time, i.e., the characteristic
decay timescale for the dust grain velocity relative to the gas. In star
forming regions, the mean free path of the gas molecules is larger than the
dust particles radius and the stopping time depends linearly on the dust grain
size $t_{\mathrm{s}}\sim s$ (Epstein, 1924). The Stokes number ${\rm St}\equiv
t_{\mathrm{s}}/t_{\mathrm{dyn}}$ characterizes the coupling of the dust grain
relative to the gas, where well coupled dust grains follow ${\rm St}<1$.
$t_{\mathrm{dyn}}$ represents a characteristic dynamical time, e.g., the eddy
turbulent crossing time for molecular clouds or the orbital time for
prostostellar discs. The Stokes number thus depends on the dust grain size so
that part of the dust grains of a given dust size distribution will be coupled
to the gas and the other part will be decoupled.

Figure 2: SPH calculations of dust dynamics during the collapse of a 1 ${\rm
M}_{\odot}$ dense core. The left column shows the gas column density (edge-on
view) and from left to right, the three other columns show the dust column
density for dust grain sizes of 10 $\mu$m, 100 $\mu$m, and 1 mm. Large dust
grains decouple from the gas and settle in the mid-plane faster into a large
dusty disc with larger dust-to-gas ratio than in the envelop. Figure adapted
from Bate and Lorén-Aguilar (2017) with permission from the authors.

Several implementations of dust and gas mixtures have been developed in the
literature to deal with the various Stokes regime. For ${\rm St}>1$, the bi-
fluid formalism seems to be the most adapted (e.g. Bai and Stone, 2010; Meheut
et al., 2012; Laibe and Price, 2012; Lorén-Aguilar and Bate, 2014; Booth et
al., 2015). However, bi-fluid algorithms encounter sever limitations in the
strong drag regime (small stopping time), where prohibitive spatial and time
discretisations may be required (Laibe and Price, 2012). In recent years,
Laibe and Price (2014a) proposed a monofluid formalism for the dust and gas
mixture which is well adapted to regimes with ${\rm St}\leq 1$. This monofluid
formalism has been implemented with success in SPH and grid-based codes in the
past five years (Laibe and Price, 2014b; Lebreuilly et al., 2019). Last, Price
and Laibe (2015) and Lin and Youdin (2017) present methods for simulating the
dynamics of small grains with only one additional equation on the dust
concentration on top of the Euler equations. These last methods have been
recently implemented with success in PHANTOM (Price et al., 2017), as well as
in PLUTO (Chen and Lin, 2018) and in RAMSES (Lebreuilly et al., 2019).

Recent works have emphasized the possible decoupling of micrometer dust grains
in molecular clouds (Hopkins and Lee, 2016; Tricco et al., 2017) as well as of
millimeter grains in collapsing dense cores (Bate and Lorén-Aguilar, 2017, see
Fig. 2). This decoupling is at the origin of a dynamical sorting which changes
the dust size distribution. Other mechanisms may also affect the size
distribution: coagulation, fragmentation, sublimation, etc. Further work
should investigate the relative importance of each of these processes on the
dust size distribution. In addition, the shape of the dust size distribution
plays a critical role in the calculations of the resistivity as well as of the
opacity. We can anticipate that coupled dust and gas dynamics will consider
the evolution of the dust size distribution, coupled to non-ideal MHD and RHD
algorithms. Last but not least, all current works on this field only account
for the hydrodynamical (the pressure gradient) and gravitational forces.
Adding the dynamics of charged dust grains is the natural next step.

##  5 Radiative transfer: numerical techniques

###  5.1 A poor’s man approach: the polytropic equation of state

Radiative transfer plays a fundamental role in star formation and cannot be
ignored, even for a review on magnetic fields. The importance of radiation
fields is two-fold. First, when the molecular cores collapse and reach a high
enough density, the dust opacity becomes large enough to absorb the cooling
radiation and the gas becomes optically thick. The traditional value for the
dust opacity is

$\kappa_{\rm dust}\simeq 0.1~{}{\rm cm}^{2}/{\rm g}~{}~{}~{}{\rm
at}~{}~{}~{}T\simeq 10{\rm K}.$ (34)

Requiring that the Jeans length of the gas becomes optically thick to infrared
radiation leads to

$\rho\kappa_{\rm dust}\lambda_{\rm J}>1~{}~{}~{}{\rm
or}~{}~{}~{}\rho>\rho_{\rm crit}=\frac{G}{\pi\kappa_{\rm dust}^{2}c_{\rm
s}^{2}}\simeq 10^{-14}~{}{\rm g/cm}^{3}.$ (35)

Many MHD simulations use this critical density to define a polytropic equation
of state of the form

$T(\rho)=10{\rm K}\left(1+\left(\frac{\rho}{\rho_{\rm
crit}}\right)^{2/3}\right),$ (36)

to capture the transition from an optically thin, isothermal gas at low
density to an optically thick, adiabatic gas at high density. Although this
simplified approach is very useful for many already expensive simulations of
collapse of turbulent clouds, this does not model properly the physics of
radiative transfer (e.g. Commerçon et al., 2010). More importantly, once
protostars have formed, higher energy radiation will emerge from the accretion
shocks on the surface of the protostars, or from the main sequence stars after
they ignited nuclear reactions. This high energy radiation in the optical or
UV bands will then interact with the parent cloud, providing radiation
feedback and influence the formation of the next generation of stars (Krumholz
et al., 2007b; Offner et al., 2009; Myers et al., 2013). In conclusion,
radiation plays a central role during the formation and the evolution of the
first Larson core (see Sec. 6), but later on acts as a self-regulating
mechanism for star formation.

###  5.2 Ray tracing and long characteristic methods

The radiative transfer equation written in full generality using the radiation
specific intensity $I_{\nu}({\bf x},{\bf n},t)$

$\frac{1}{c}\frac{\partial}{\partial t}I_{\nu}+{\bf n}\cdot\nabla
I_{\nu}=j_{\nu}-\alpha_{\nu}I_{\nu}.$ (37)

bf where $j_{\nu}$ is the emission coefficient and $\alpha_{\nu}$ is the
absorption coefficient. Note that the radiative transfer equation is written
here in the laboratory frame. The absorption and emission coefficients are
however defined in the comoving frame. Properly introducing the relativistic
corrections to the radiative transfer equation and its various moments is a
very important aspect of the problem we will only briefly touch upon in this
review.

The most natural numerical technique to solve the radiative transfer equation
is ray-tracing. The idea is to shoot a discrete set of light rays from a point
source, solving the previous equation along the path of the light ray using
the curvilinear coordinate $s$ defined as ${\rm d}s=c{\rm d}t$. The equation
of radiative transfer can be written as a Lagrangian time derivative following
the trajectory of the photons as

$\frac{{\rm d}}{{\rm d}s}I_{\nu}=j_{\nu}-\alpha_{\nu}I_{\nu}.$ (38)

The light ray is discretised along its path, matching the underlying gas
distribution, whether it is a grid or a set of SPH particles. The difficulty
is then to choose the appropriate number of rays, so that the angular
distribution of the radiation field is properly sampled, but also so that
enough light rays intersect the gas elements. In Abel and Wandelt (2002), an
adaptive method was proposed, that splits or merges light rays adaptively,
using the Healpix tessellation to sample the sphere. This method was first
developed for the ENZO code (Wise and Abel, 2011) and used with success to
describe photo-ionization regions at the Epoch of Reionization (Wise et al.,
2012), as well as in turbulent molecular clouds (Shima et al., 2017). A
similar scheme based on the Monte-Carlo approach has been developed for SPH in
Altay et al. (2008) and Pawlik and Schaye (2008).

A similar approach is provided by the long characteristics method. This method
can describe light rays emanating from a single point source like the ray-
tracing method, but also diffuse radiation by adopting for each direction a
set of parallel rays whose spacing matches the local grid resolution. This
method has been implemented in the FLASH AMR code by Rijkhorst et al. (2006),
with however an important addition we will discuss later. Diffuse radiation
using long characteristics was included later by Buntemeyer et al. (2016) in
the FLASH code, with an emphasis of describing both the optically thin and
optically thick regimes. These methods exploit the particular property of the
AMR implementations of FLASH and ENZO, namely large rectangular patches of
various sizes and resolution. The long characteristic method has been
developed for the graded octree AMR code RAMSES only recently by Frostholm et
al. (2018).

###  5.3 Monte Carlo and short characteristic methods

Although long characteristics offer the advantage of solving the exact
attenuation along light rays, an especially accurate approach for point
sources, they are tricky to parallelise because light rays propagate over
multiple MPI domains. An alternative technique is the short characteristics
method, for which light rays as short as the hydrodynamics cell size are
considered. The radiation intensity from neighboring light rays is
interpolated between neighboring cells, avoiding the use of long
characteristics. Rijkhorst et al. (2006) use an hybrid long- and short-
characteristics method in the FLASH code to overcome the difficulty to
parallelise the long characteristics method. The short characteristic method
also offers better scaling properties. The short characteristics method has
also been used in the ATHENA code to solve the stationary radiative transfer
equation by Davis et al. (2012) and in the C2-ray code for ionization problems
(Mellema et al., 2006).

A third method which shares with the two previous ones the property of being
accurate and able to capture complex radiation geometry is the Monte Carlo
method. The Monte Carlo approach samples the radiation field with photon
packets carrying information such as angle of propagation and frequency. This
is as close as one can be from true photons propagating in the fluid and
allows for complex scattering processes to be included in the models (Ercolano
et al., 2003; Altay et al., 2008; Harries, 2011; Roth and Kasen, 2015). The
method can be however quite expensive, as large particle numbers are required
to avoid that the Poisson noise dominates the computed signal. In the
diffusion limit, the mean free path becomes much smaller than the length scale
of the system and the number of interactions becomes prohibitively large: It
scales as the square of the optical depth (Nayakshin et al., 2009). It has
also been mostly developed for stationary problems but many codes able to deal
with non-stationary problems are now emerging (Nayakshin et al., 2009;
Harries, 2011; Lomax and Whitworth, 2016).

###  5.4 Two moments methods

The common aspect in these three methods is that they can be quite expensive.
This is probably not true for the short characteristics method, although the
number of angular domains adopted plays a role in the final accuracy of the
solution, and affects the cost. An alternative approach, called moment-based
radiative transfer, can potentially solve this cost problem by integrating the
radiative transfer equation over the angles. We obtain the radiation energy
and the radiation momentum equation as

$\frac{\partial}{\partial t}E_{\nu}+\nabla\cdot{\bf F}_{\nu}=4\pi
j_{\nu}-\alpha_{\nu}cE_{\nu},$ (39) $\frac{1}{c}\frac{\partial}{\partial
t}{\bf F}_{\nu}+c\nabla\cdot{\mathbb{P}}_{\nu}=-\alpha_{\nu}{\bf F}_{\nu},$
(40)

where $E_{\nu}$ is the radiation energy, ${\bf F}_{\nu}$ is the radiation flux
and ${\mathbb{P}}={\mathbb{D}}E_{\nu}$ is the radiation pressure tensor and
${\mathbb{D}}$ is the Eddington tensor. For slowly varying spectral absorption
and emission coefficients, for which one has

$\nu\frac{\partial\alpha_{\nu}}{\partial\nu}\ll\alpha_{\nu}~{}~{}~{}{\rm
and}~{}~{}~{}\nu\frac{\partial j_{\nu}}{\partial\nu}\ll j_{\nu},$ (41)

one can write relativistic correction to first order in $(v/c)$ as

$\alpha_{\nu}({\bf n})=\alpha^{0}_{\nu}-({\bf n}\cdot{\bf
v}/c)\alpha^{0}_{\nu},$ (42) $j_{\nu}({\bf n})=j^{0}_{\nu}+2({\bf n}\cdot{\bf
v}/c)j^{0}_{\nu},$ (43)

where superscript $0$ refers to quantities measured in the comoving frame (see
Mihalas and Klein, 1982; Mihalas and Mihalas, 1984; Krumholz et al., 2007a,
for a detailed discussion). Injecting this in the radiative transfer equation
and taking the moments over the angles leads to new terms accounting for the
relativistic corrections

$\frac{\partial}{\partial t}E_{\nu}+\nabla\cdot{\bf F}_{\nu}=4\pi
j^{0}_{\nu}-\alpha^{0}_{\nu}cE_{\nu}+\alpha^{0}_{\nu}{\bf
F}_{\nu}\cdot\frac{{\bf v}}{c},$ (44) $\frac{1}{c}\frac{\partial}{\partial
t}{\bf F}_{\nu}+c\nabla\cdot{\mathbb{P}}_{\nu}=-\alpha^{0}_{\nu}{\bf
F}_{\nu}+\alpha^{0}_{\nu}{\mathbb{P}}~{}{\bf v}+4\pi j^{0}_{\nu}\frac{{\bf
v}}{c}.$ (45)

Since $E_{\nu}$ and ${\bf F}_{\nu}$ are still expressed in the laboratory
frame but $\alpha^{0}_{\nu}$ and $j^{0}_{\nu}$ are expressed in the comoving
frame, these equations are often called mixed frame radiation moments
equations. We can also Lorentz-transform the radiation moments from the
laboratory frame to the comoving frame, which gives, to first order in $(v/c)$

$E^{0}_{\nu}=E_{\nu}-\frac{2}{c}{\bf F}_{\nu}\cdot\frac{{\bf v}}{c},$ (46)
${\bf F}^{0}_{\nu}={\bf F}_{\nu}-E_{\nu}{\bf v}-{\mathbb{P}}~{}{\bf v}.$ (47)

Injecting these relations only in the right-hand side of the moments equations
and neglecting all $(v/c)$ terms gives the simpler form

$\frac{\partial}{\partial t}E_{\nu}+\nabla\cdot{\bf F}_{\nu}=\left(4\pi
j^{0}_{\nu}-\alpha^{0}_{\nu}cE^{0}_{\nu}\right)-\alpha^{0}_{\nu}{\bf
F}^{0}_{\nu}\cdot\frac{{\bf v}}{c},$ (48) $\frac{1}{c}\frac{\partial}{\partial
t}{\bf F}_{\nu}+c\nabla\cdot{\mathbb{P}}_{\nu}=-\alpha^{0}_{\nu}{\bf
F}^{0}_{\nu}+\left(4\pi
j^{0}_{\nu}-\alpha^{0}_{\nu}cE^{0}_{\nu}\right)\frac{{\bf v}}{c},$ (49)

where all right-hand side terms are now expressed in the comoving frame. In
practice, after the radiation transport step, encoded in the left-hand side,
one then converts the laboratory frame radiation variables into the comoving
frame and solve the thermo-chemistry encoded in the right-hand side in the
comoving frame. Local thermodynamical equilibrium conditions can then be
reached in the comoving frame, eventually reaching the diffusion limit (see
below). At the end of the thermo-chemistry step, one finally converts back the
radiation variables into the laboratory frame before entering the next
transport step (see Rosdahl and Teyssier, 2015, for implementation details).

These equations are quite handy because of their simplicity: We have no
angular dependence anymore. There is a catch however, since everything depends
on the Eddington tensor, for which we have no equation at this order of the
moment’s hierarchy. At this point, two different strategies have been explored
in the star formation literature: 1- compute an accurate Eddington tensor by
solving exactly the stationary radiative transfer equation using one of the
methods presented above (short or long characteristics, ray-tracing or Monte
Carlo), 2- compute an approximate form of the Eddington tensor based on a
particular closure of the moment’s hierarchy.

The first approach is often referred to a the Variable Eddington Tensor (VET)
method. In Gnedin and Abel (2001), the authors developed a moment-based solver
coupled to a cosmological hydrodynamics code, called Optically Thin Variable
Eddington Tensor or OTVET. A similar technique has also been implemented for
the SPH code GADGET (Petkova and Springel, 2011). The idea is to solve the
stationary radiative transfer equation in the full domain assuming an
optically thin medium, so that the problem boils down to a collection of
$1/r^{2}$ sources combined together. This can be achieved efficiently using
any gravity solver. Once the corresponding Eddington tensor is computed, the
moments equations are solved using a finite difference scheme. In Jiang et al.
(2012), on the other hand, the authors developed a moment-based solver for the
ATHENA MHD code, using VET together with the short characteristic method of
Davis et al. (2012) to compute the Eddington tensor. Obviously, this leads to
significantly more accurate results in optically thicker environments. An
intermediate solution has been implemented in the TreeCol code (Clark et al.,
2012), for which a gravity tree code is used to collect column densities
between particles.

The second approach relies on simple but approximate models for the Eddington
tensor. The simplest one, called M0, assumes the Eddington tensor is
isotropic, with

${\rm M0~{}closure:}~{}~{}~{}{\mathbb{D}}=\frac{1}{3}{\mathbb{I}},$ (50)

as it is the case for example in the diffusion limit (see below). A slightly
more elaborate model is called the M1 closure (Dubroca and Feugeas, 1999;
Ripoll et al., 2001). It assumes that the radiation intensity can be fitted by
a Lorentz-boosted Planckian distribution, in other words an ellipse in angular
space. The parameters of this distribution are determined by matching the
radiation energy and the radiation flux. This leads to

${\rm
M1~{}closure:}~{}~{}~{}{\mathbb{D}}=\frac{1-\chi}{2}{\mathbb{I}}+\frac{3\chi-1}{2}{\bf
n}\otimes{\bf n}.$ (51)

This closure captures the optically thick regime, but also optically thin
conditions, in case there is one dominant source of radiation. Indeed, in this
case, one gets the free-streaming regime for radiation with
${\mathbb{D}}\simeq{\bf n}\otimes{\bf n}$. The unit vector ${\bf n}$ and the
parameter $\chi$ are determined using the radiation energy and the radiation
flux as

${\bf n}=\frac{{\bf F_{\nu}}}{\left|{\bf
F_{\nu}}\right|},~{}~{}~{}\chi(f)=\frac{3+4f^{2}}{5+2\sqrt{4-3f^{2}}},~{}~{}~{}f=\frac{\left|{\bf
F_{\nu}}\right|}{cE_{\nu}}$ (52)

A particularly interesting property of the M1 model is that it leads to a
system of hyperbolic conservation laws, that can be numerically integrated
using Godunov schemes. The M1 model has also some major caveats, such as
radiation shocks, in case multiple strong sources are present, or incorrect
radiation geometry, in case of sharp shadows. This model was introduced for
the first time in astrophysics by González et al. (2007) in the HERACLES code.
It was then adapted to cosmic reionization by Aubert and Teyssier (2008) in
the ATON code, and ported to GPU into the RAMSES-CUDATON code (Aubert and
Teyssier, 2010), used now routinely for simulations of the Epoch of
Reionization (EoR) (e.g. Ocvirk et al., 2016). Later on, the M1 method has
been ported to the AMR code RAMSES, leading to the design of the RAMSES-RT
solver (Rosdahl et al., 2013; Rosdahl and Teyssier, 2015). The M1 method has
also been ported to ATHENA by Skinner and Ostriker (2013) and recently to
AREPO by Kannan et al. (2018), the latter based on a new higher-order
implementation of the M1 closure.

###  5.5 Flux-limited diffusion methods

A definitive advantage of moment-based radiative transfer methods is the
possibility to model both optically thin and optically thick regimes. For
optically thick conditions, the radiation field follows the diffusion limit of
the radiative transfer equation. In that limit, the radiation field in the
comoving frame is quasi-isotropic, so that
${\mathbb{D}}\simeq\frac{1}{3}{\mathbb{I}}$ and the radiation flux can be
approximated by

${\bf F}^{0}_{\nu}\simeq-\frac{c}{3\alpha^{0}_{\nu}}\nabla E^{0}_{\nu}.$ (53)

Since now the flux in the comoving frame is a direct function of the energy
density in the comoving frame, we do not need the radiation flux equation
anymore and we are left with only the radiation energy equation. The Lorentz
transform can be written to leading order as

$E_{\nu}\simeq E^{0}_{\nu}$ (54) ${\bf
F}_{\nu}\simeq-\frac{c}{3\alpha^{0}_{\nu}}\nabla
E^{0}_{\nu}+\frac{4}{3}E^{0}_{\nu}{\bf v}.$ (55)

Injecting these into the radiation energy equation leads to a fully covariant
formulation in the diffusion limit

$\frac{\partial}{\partial
t}E^{0}_{\nu}+\nabla\cdot\left(\frac{4}{3}E^{0}_{\nu}{\bf
v}\right)-\nabla\cdot\left(\frac{c}{3\alpha^{0}_{\nu}}\nabla
E^{0}_{\nu}\right)=\left(4\pi
j^{0}_{\nu}-\alpha^{0}_{\nu}cE^{0}_{\nu}\right)+\frac{1}{3}\nabla
E^{0}_{\nu}\cdot{\bf v},$ (56)

which can be written in the familiar form

$\frac{\partial}{\partial t}E^{0}_{\nu}+\nabla\cdot\left(E^{0}_{\nu}{\bf
v}\right)+\frac{1}{3}E^{0}_{\nu}\nabla\cdot{\bf
v}=\nabla\cdot\left(\frac{c}{3\alpha^{0}_{\nu}}\nabla
E^{0}_{\nu}\right)+\left(4\pi
j^{0}_{\nu}-\alpha^{0}_{\nu}cE^{0}_{\nu}\right).$ (57)

This is the radiation energy equation in the diffusion limit using only
comoving radiation variables. If one wants to extend this equation outside its
validity range (high optical depth), then one can use Flux Limited Diffusion
(FLD), for which the flux function is modified as

${\bf F}^{0}_{\nu}\simeq-\frac{c\lambda(R)}{\alpha^{0}_{\nu}}\nabla
E^{0}_{\nu}~{}~{}~{}{\rm and}~{}~{}~{}R=\frac{\left|\nabla
E^{0}_{\nu}\right|}{\alpha^{0}_{\nu}E^{0}_{\nu}},$ (58)

where $\lambda(R)$ is the flux limiter, a function that has to connect the
diffusion limit with $\lambda\simeq 1/3$ for $R\simeq 0$ to the free-streaming
regime where $\lambda\simeq 1/R$ for $R\rightarrow+\infty$. A possible form
for the flux limiter has been proposed by Minerbo (1978) and reads

$\displaystyle\lambda=\frac{2}{3+\sqrt{9+12R^{2}}}~{}~{}~{}$
$\displaystyle{\rm for}$ $\displaystyle~{}~{}~{}0\leq R\leq\frac{3}{2},$ (59)
$\displaystyle\lambda=\frac{1}{1+R+\sqrt{1+2R^{2}}}~{}~{}~{}$
$\displaystyle{\rm for}$ $\displaystyle~{}~{}~{}\frac{3}{2}\leq R\leq+\infty.$
(60)

while the Levermore (1984) flux limiter reads

$\displaystyle\lambda=\frac{1}{R}\left(\coth R-\frac{1}{R}\right)$ (61)

Historically, FLD has been the first method implemented in radiation
hydrodynamics codes. Here again, the ZEUS 2D code was a precursor (Stone et
al., 1992b). For AMR, the ORION code was developed specifically for FLD in
collapsing star forming core by Krumholz et al. (2007a), using the mixed frame
formulation outlined above or in Mihalas and Klein (1982). A version of the
RAMSES code with FLD was developed later by Commercon et al. (2011), with a
particular emphasis on adaptive time stepping together with implicit time
integration described in Commercon et al. (2014). Because of its simplicity,
FLD has been used for many detailed studies of gravitational collapse with the
effect of radiation included. The current frontier for FLD is probably the
multigroup treatment of radiation, allowing for more accurate models of the
full spectral energy distribution (Shestakov and Offner, 2008; Zhang et al.,
2013; González et al., 2015). In the context of multifrequency radiative
transfer, hybrid solutions have also been proposed, with ray-tracing
techniques dedicated to the high-energy radiation field, and FLD for the
infrared, dust-absorbed radiation field (Kuiper et al., 2010a; Klassen et al.,
2014).

FLD has also been developed for SPH by Whitehouse and Bate (2004), Whitehouse
et al. (2005) and Mayer et al. (2007), where most of the difference with the
many grid-based versions lies in the radiation energy gradient that uses the
SPH kernel. A simplified version based on a local estimate of the column
density has also been proposed by Stamatellos et al. (2007) which avoids
entirely the need for an implicit radiation solver. This method can be
considered as an intermediate solution between the polytropic equation of
state that we introduced at the beginning of this section and FLD. The main
caveat here is the complete lack of any proper transport mechanism for
radiation. In a similar spirit, Dale et al. (2007) also developed an
approximate method for ionizing radiation, with Stromgren sphere iteratively
grown around each star to locate photo-ionized, photo-heated gas.

##  6 Second collapse

From first principles, the ultimate goal of star formation is the formation of
the protostar itself. The major difficulty sits in the huge dynamical range in
physical and temporal scales that have to be described: protostars with radii
of about 1 ${\rm R}_{\odot}$ form within dense cores of sizes $\simeq 0.1$ pc.
So there are two ways to deal with protostars in star formation. The first one
is to try to follow the collapse down to the stellar scales with the best
numerical model accounting for the complex physical processes at stack
combined with a very high numerical resolution at a cost of very short horizon
of predictability. The second is the opposite: the physics and numerical
resolutions are degraded, but the models are integrated over longer dynamical
timescales to study the impact of protostellar feedback on the ISM dynamics.
In the next three sections, we review the work that has been done using these
two approaches, as well as attempts to account for protostellar evolution in
large-scale models.

###  6.1 Historical work and 1D studies

Contemporary numerical star formation studies started in the late 60s with the
outstanding pioneering work of Larson (1969) who first computed numerically
the collapse of a dense core down to the formation of the protostar. He used a
1D spherically symmetric model, accounting for coupled gas dynamics and
radiative transfer in a modified Eulerian scheme. Larson identified two
distinct stages during the protostellar collapse, called the first and the
second collapse. Each of these two stages are followed by the formation of a
hydrostatic object, commonly referred to as the first and second Larson cores
in post Larson (1969) studies. Larson’s work established an empirical
evolutionary sequence as follows. Dense cores first collapse isothermally
because dust thermal emission is very efficient at radiating away the gas
compression energy (dust and gas are thermally coupled within dense cores, see
for instance Galli et al., 2002). The first collapse stops at densities of the
order of $\simeq 10^{-13}$ g cm-3 at which dust grains become opaque to their
own radiation so that the gas begins to heat up quasi-adiabatically. This is
the formation of the first Larson core. The first collapse lasts roughly a
free-fall time. The first core accretes matter and its temperature increases
up about 2000K where H2 dissociation starts. This endothermic reaction enables
the gas to behave as a fluid with an effective polytropic index smaller than
the critical value $4/3$ for collapse. The second collapse thus starts at
typical densities of $\simeq 10^{-8}$ g cm-3 until stellar densities are
reached within $\simeq 100$ yr to form the second core, i.e., the protostar.
The typical properties of first Larson cores at the start of the second
collapse are a size $\simeq 10$ AU, a mass $\simeq 0.01$ ${\rm M}_{\odot}$,
and a lifetime $\simeq 1000$ yr. We note that first Larson cores are predicted
by theoretical and numerical studies, but there is no observational
confirmation of such objects even though a few sources are considered as good
candidates (Tsitali et al., 2013; Gerin et al., 2015; Maureira et al., 2017).

Since Larson’s work numerous 1D calculations using spherical symmetry have
been conducted. We note the work in the 80s by Stahler et al. (1980) and
Winkler and Newman (1980) in which Larson predictions were confirmed to
establish the current empirical evolution sequence of collapsing low-mass
dense cores. This evolutionary picture is still currently the commonly
admitted scenario for low-mass star formation. For massive star formation,
recent work indicates that first cores may not have time to form because of
the large accretion rates (Bhandare et al., 2018). Note that in the reminder
of the section dedicated to second collapse, we mention only the work which is
consistent with the Larson evolutionary picture, i.e., where the gas thermal
and chemical budgets are taken into account, and we do not mention the work
that has been done using an isothermal approximation.

Modern studies began with the calculations performed in Masunaga and Inutsuka
(2000) who incorporated a more accurate radiation transport scheme as well as
a realistic gas equation of state which accounts for H2 dissociation. Masunaga
and Inutsuka (2000) integrated their models throughout the Class 0 and Class I
phases up to an age of $1.3\times 10^{5}$ yr. Nowadays, 1D models are still
used either as a first step towards describing more accurately the physics of
the protostellar collapse, e.g., with multigroup radiative transfer (Vaytet et
al., 2013, 2014) or with a view to provide quantitative predictions for the
first and second Larson core properties (Vaytet and Haugbølle, 2017; Bhandare
et al., 2018). These models are not pushed in time as far as the ones by
Masunaga and Inutsuka (2000).

###  6.2 From 1D to 3D

Pioneered studies acknowledged up front that 1D spherical models were limited
because the effect of rotation, turbulence, and magnetic fields cannot be
taken into account. Three decades after Larson work, the first 3D numerical
simulation able to describe the formation of the protostar was performed by
Bate (1998) using SPH. The model accounted for initial rotation, but magnetic
fields were neglected and radiative transfer was crudely mimicked by a
piecewise polytropic equation of state. It is worth to notice at this point
that grid-based codes took some time to reproduce Bate (1998) results for
almost one decade. SPH is indeed naturally well-suited to collapse problems
thanks to its Lagrangian nature while innovative techniques, such as AMR and
nested grids, were not mature enough to capture the entire dynamical range and
to cover eight orders of magnitude in physical length given the computer
capabilities at that time. As mentioned in the Sec. 3, ideal MHD was
introduced early in protostellar collapse models using grid-based codes (e.g.,
Dorfi, 1982; Phillips and Monaghan, 1985).

The first calculations of second collapse incorporating magnetic fields were
thus done on a grid-based code by Tomisaka (2002) using a nested grid with
cylindrical coordinates assuming axisymmetry (Tomisaka and Bregman, 1993).
Tomisaka (2002) code solved the ideal MHD equations using a second-order
“monotonic scheme” (van Leer, 1977; Norman and Winkler, 1986), a constrained
transport scheme for the induction equation (Evans and Hawley, 1988), as well
as an angular momentum preserving scheme. To capture the huge range in spatial
scales, Tomisaka (2002) used up to 16 levels of refinement. He found that two
different types of outflows were launched at the first and second cores
scales. He suggested that the second core outflow corresponds to the optical
jets observed in young stellar objects (YSOs), while the first core outflow
corresponds to the molecular bipolar outflow. The first full 3D MHD models
were presented in Machida et al. (2006) using nested grid with Cartesian
coordinates and 21 levels of refinement. They compared the results obtained
with resistive (Ohmic diffusion) and ideal MHD. They confirmed the results of
Tomisaka (2002) even with non-ideal MHD included. After Machida et al. (2006)
work, the numerical challenge of forming a protostar in 3D with MHD was faced,
but the horizon of predictability remained very short since they had to stop
their calculations 20 days after the protostar formation because of extremely
small timesteps resulting from the high velocity of the jet propagating from
the finest level. Banerjee and Pudritz (2006) also performed 3D collapse
calculations down to stellar scales, using ideal MHD as well as molecular
cooling, and they confirmed the results of Tomisaka (2002) and Machida et al.
(2006).

All the multidimensional studies we mentioned so far were done assuming a
piecewise polytropic equation of state to mimic the thermal behaviour of the
gas throughout the evolutionary sequence. Such barotropic equation of states
are parameterized using the results of 1D spherical symmetry simulations
previously mentioned, so that they cannot account for multidimensional effects
such as the optical depth drop in the vertical direction that allows efficient
cooling of the nascent protostellar disks and thus leads to incorrect results
for the fragmentation (e.g., Boss et al., 2000; Commerçon et al., 2010; Tomida
et al., 2010).

The next level of complexity is to add a more accurate model for the thermal
balance of the gas and dust mixture, i.e., to model radiative transfer (see
Sec. 5). The addition of a more accurate model for radiation transport can be
crudely resumed to an additional heating/cooling term in the gas internal
energy equation due to matter/radiation interactions. For the models which are
restricted to the first collapse and first core formation, i.e., for
temperatures less than $\simeq 2000$ K, the addition is conceptually
straightforward for the gas evolution equations. An usual ideal gas equation
of state can be used, with the only modification on the adiabatic index
$\gamma$ which varies with the temperature. At low temperature ($T\lesssim
150$ K), the H2 molecule behaves like a monoatomic gas ($\gamma=5/3$) because
only the vibrational degrees of freedom are excited. At higher temperature,
the rotational degrees start to be excited and the adiabatic index decreases
to $\gamma=7/5$. Accounting for this non-constant adiabatic index is important
regarding the fragmentation properties (Commerçon et al., 2010) and the first
core properties (Tomida et al., 2010; Lee and Hennebelle, 2018a). The driver
of the second collapse is the endothermic dissociation of the H2 molecule,
which modifies considerably the gas equation of state.

Currently, there are two ways to integrate non-ideal EOS in second collapse
calculations. The first one is to compute on-the-fly the thermodynamic
quantities from equilibrium chemical abundances and to use an ideal gas
equation of state, including the effects of dissociation and ionisation, to
compute the gas pressure and the internal energy (Black and Bodenheimer,
1975). Most of the second-collapse including RHD are done with this
approximation (Masunaga and Inutsuka, 2000; Whitehouse and Bate, 2006;
Stamatellos et al., 2007; Forgan et al., 2009; Tomida et al., 2013; Bhandare
et al., 2018). The second one consists in using tabulated EOS table coming
from detailed studies to account for non-ideal effects (e.g., ionisation by
pressure, interaction between particles) which are important at high
density/pressure (Saumon et al., 1995). Saumon et al. (1995) EOS is used in
the series of papers by Vaytet et al. (2013, 2014); Vaytet and Haugbølle
(2017); Vaytet et al. (2018), and also in the SPH public code PHANTOM (Price
et al., 2017).

We note that these developments imply that the dependency of the specific heat
capacity as a function of the physical conditions has to be taken into
account. Caution should thus be taken within the radiation solvers which are
often integrated using implicit schemes. The most common and simple way to
deal with it is to compute a heat capacity-like factor,
$\tilde{C_{\mathrm{V}}}\equiv e/T$ where $e$ is the gas internal energy
density, which is kept constant in the radiation solve. Readers can refer to
Tomida et al. (2013) for a detailed implementation. Unfortunately, details on
the numerical implementations are usually not reported in the literature.

The first full 3D RHD calculations of the second collapse of an initial 1
${\rm M}_{\odot}$ dense core were performed by Whitehouse and Bate (2006)
using the FLD approximation in a SPH code, followed by Stamatellos et al.
(2007) who used a local radiative cooling approximation and SPH. Bate (2010)
extended the work of Whitehouse and Bate (2006) beyond the formation of the
stellar core for about 50 yr, using $3\times 10^{6}$ equal-mass SPH particles
to satisfy the Bate and Burkert (1997) mass resolution criterion throughout
the entire collapse. Interestingly, there is to date no 3D RHD calculations of
the second collapse, i.e., without magnetic field, performed with a grid-based
code.

Last but not least, it has been demonstrated in Commerçon et al. (2010) that
the interplay between magnetic fields and radiative transfer is of primary
importance for the thermal budget of the collapsing gas and of the first and
second core accretion shocks (Commerçon et al., 2011a; Vaytet et al., 2012,
2013, 2018). The infall velocity is indeed greatly modified in the presence of
magnetic fields which tends to focus the collapsing gas on a smaller area
compared to the case without magnetic fields (Commerçon et al., 2011b; Myers
et al., 2013). The magnetic braking mechanism transports angular momentum from
the inner parts of the collapsing cloud to the envelop. As a consequence, the
infall velocity increases along with the incident kinetic energy, which then
modifies the thermal budget through the accretion shocks. The next step in the
increasing complexity is to perform full 3D RMHD calculations. The first SPH
RMHD calculations of protostellar core formation were performed in Bate et al.
(2014), extending their previous work to ideal MHD. Tomida et al. (2013, 2015)
performed the first 3D full RMHD calculations using a grid-based code,
extending the work of Machida et al. (2006) to account for radiative transfer
with the FLD approximation, as well as the Ohmic dissipation and the ambipolar
diffusion. They used a finite-volume nested-grid code with standard MUSCL
predictor-corrector scheme to achieve second-order accuracy.

More details on Tomida et al. (2015) implementation for non-ideal MHD are
given in Sec. 4. For the implicit radiation update, they used a linear system
solver with a combination of the BiCGStab solver and the incomplete LU
decomposition preconditionner. Tomida et al. (2013) used a nested-grid
constructed to ensure that the Jeans length is always resolved by at least 16
cells, with typical spatial resolution at the maximum level of refinement of
$\Delta x\simeq 6.6\times 10^{-5}$ AU (23 levels). When this resolution
condition cannot be satisfied, they stopped the calculations. This is
particularly the case when discs form and grow with time to spread over
different level of refinements. As a consequence, they cannot follow the
second collapse in their model with ambipolar diffusion since a large disc is
formed. More flexible refinement techniques, such as AMR or SPH are thus
required to continue the calculations with a reasonable numerical resolution.
Currently, the state-of-the-art includes 3D RMHD calculations with resistive
effects, done with both SPH (Tsukamoto et al., 2015b; Wurster et al., 2018)
and AMR (Vaytet et al., 2018). The typical resolution in SPH calculations
currently achievable with limited CPU resources is of $3\times 10^{6}$ equal-
mass particles per ${\rm M}_{\odot}$. In AMR, Vaytet et al. (2018) used a
coarse grid of $64^{3}$ with 21 additional levels of refinement, and a
refinement criterion of 32 points per Jeans length. While the SPH calculations
are able to integrate the three non-ideal effects (Wurster et al., 2018) down
to stellar scales, there is no AMR work describing the Hall effect in second
collapse calculations (Vaytet et al. (2018) accounts for the ambipolar
diffusion and the Ohmic diffusion). The calculations by Wurster et al. (2018)
were integrated 17 yr after the stellar core formation, while Vaytet et al.
(2018) could only integrate for 24 days. Beside the major achievement which
represents the formation of a protostar on a computer, we must admit that
these kind of studies are not yet capable to give long term predictability. In
addition, the CPU time required to integrate full 3D RMHD models is relatively
big.

Clearly, second collapse study cannot be done on a frame covering a large
parameter space. In addition, even at the first core scale, the integration
timestep becomes so short in the adiabatic fragments that long time evolution
calculations becomes prohibitive. In that view, the sink particles are a good
alternative to reach a long horizon of predictability and to cover a wide
parameter space.

##  7 Sink particles

Sink particles are Lagrangian particles, and can accrete matter as well as
angular momentum in order to ensure mass and angular momentum conservation.
Sink particles should not be confounded with sink cells which are usually
introduced in grid-based codes which use a spherical grid to deal with the
singularity in the center (e.g., Boss and Black, 1982; Kuiper et al., 2010b;
Li et al., 2011). The differences between sink particles and sink cells are
discussed in Sec.7.4. Sink particles are often erroneously referred as
protostars or stars in the literature whereas the resolution at which they are
introduced is generally much larger than the stellar scales, typically $1-10$
AU for collapse calculations, and a few $100$ AU for cluster formation and ISM
evolution calculations. Multiple systems (small clusters, binaries), as well
as systems made of a disk around a star can thus be encompassed within a
single sink particle depending on the resolution.

###  7.1 Standard implementations

Sink particles were introduced in star formation calculations by Bate et al.
(1995) to enable long time integration to study the fragmentation of
collapsing clouds. Sink particles were introduced to mimic the second collapse
and the protostars at scales that could not be reached in numerical
simulations. They are treated as accreting non-collisional point masses and
enable one to follow the evolution through the main accretion phase, with a
view to compare with observations of clusters in which most of the gas has
been accreted. Their original implementation was done in a SPH code, with
three main parts in the algorithm: i) sink particle creation, ii) sink
particle accretion, iii) boundary conditions for sink particles. A SPH gaseous
particle is tagged for sink particle creation if its density exceeds a
specific density threshold (defined by the user) and satisfied a number of
tests. First, the smoothing length of the SPH particle that is a candidate for
sink particle creation has to be less than half the accretion radius
$r_{\mathrm{acc}}$ of the sink particle. This ensures that the sink particle
is formed from at least $N_{\mathrm{neigh}}$. $r_{\mathrm{acc}}$ is defined
before the calculations start and remains fixed. Its value is chosen by the
user depending on the required level of resolution. It sets the smallest scale
that can be resolved in the calculations. The knowledge of the flow within
$r_{\mathrm{acc}}$ is lost and the gas is assumed to collapse beyond this
scale to form protostars. Then, a series of tests is performed on the system
composed of the particle and its neighbours to decide if it should create a
sink particle. The tests ensure that the gas particles would continue to
collapse if the sink particle was not created. First, the ratios $\alpha$ and
$\beta$ of the thermal energy and rotational energy to the gravitational
energy must satisfy

$\alpha\leq 1/2~{};~{}\alpha+\beta\leq 1.$ (62)

Additionally, the total energy must be negative, and finally, the system must
be contracting ($\nabla\cdot\vec{v}<0$). If all the tests are passed, a sink
particle is formed from the system at the centre of mass. Initially, the sink
particle has thus a mass of $\sim N_{\mathrm{neigh}}$ times the standard SPH
particle mass. Afterwards, sink particles interact only with the gas particles
via gravity. The sink particle is allowed to accrete mass at each integration
timestep if a gas particle enters $r_{\mathrm{acc}}$ and passes several
criteria: the particle must be the most bound to the sink particle (and not to
another more distant but more massive sink particle), and the particle must
have a moderate specific angular momentum. When a gas particle is accreted,
its mass and linear momentum are added to the sink particle, as well as the
angular momentum which modifies the spin of the sink particle. The sink
particle is moved to the center of mass of the sink and the accreted gas
particle. Sink particles also contribute to the computation of the total
gravitational potential in the same way as standard gaseous particles. Bate et
al. (1995) pointed out that the introduction of sink particles in SPH
calculations may affect the gas outside the accretion radius, in particular
because there is a discontinuity in the number of particles across the
accretion radius. They proposed different types of boundary conditions to
account for missing neighbours. Note that the use of boundary conditions in
modern sink particles implementation is not used anymore thanks to the
improvements that have been made, in particular in the case of optically thick
flows (e.g., Hubber et al., 2013).

Sink particles were introduced in Eulerian calculations for star formation
purposes by Krumholz et al. (2004) in the AMR code ORION. They introduce a
Jeans length criterion for sink creation, based on the resolution study of
Truelove et al. (1997). A sink is created within a cell at the maximum level
of refinement if the density $\rho$ exceeds the maximum Jeans density
$\rho_{\mathrm{J}}$ that can be resolved

$\rho_{\mathrm{J}}<J^{2}\frac{\pi c_{\mathrm{s}}^{2}}{G\Delta
x^{2}_{\mathrm{min}}},$ (63)

where $J<0.25$ is the Jeans length resolution criterion, $c_{\mathrm{s}}$ the
gas sound speed, and $\Delta x_{\mathrm{min}}$ the minimum cell size. The
initial mass of the sink is $m_{\mathrm{sink}}=(\rho-\rho_{\mathrm{J}})\Delta
x_{\mathrm{min}}^{3}$. In Krumholz et al. (2004), no additional checks are
performed to validate the sink creation. Their algorithm is coupled to a sink
merging scheme based on a friends-of-friends (FOF) algorithm to handle
situations where a block of contiguous cells create sink particles in a single
time step. The FOF is performed after each time step to group all the sink
particles (old and new), with a linking length equal to $r_{\mathrm{acc}}$.
All the groups of sink particles found by the FOF are then merged and replaced
by a single particle at the center of mass of the group. All the merged sink
particles quantities (mass, momentum, angular momentum) are also added
conservatively to the new sink particle. Krumholz et al. (2004) sink accretion
scheme was designed to handle situations where the flow onto the sink particle
is subsonic (which is analogous to the case where boundary conditions have to
be introduced in standard SPH implementation). They set the accretion rate
using an approximate formula from Bondi-Hoyle accretion which handles also the
regime where the flow is supersonic. The accretion rate is estimated using
average properties in the region within $r_{\mathrm{acc}}$. An important point
in Krumholz et al. (2004) implementation is the description of the accretion
zone. The mass accreted by the sink particle comes from all the cells in this
accretion zone. Based on their experiments, Krumholz et al. (2004) use a value
of $r_{\mathrm{acc}}=4\Delta x$. In the case where the Bondi-Hoyle radius is
smaller than the accretion zone, it is incorrect to estimate the accretion
rate from all the cells in the entire accretion zone. Krumholz et al. (2004)
thus use an accretion kernel where each cell within $r_{\mathrm{acc}}$ is
assigned a weight from a Gaussian-like function. The kernel is used to compute
the average quantities to determine the mass accretion rate. The mass to be
accreted is redistributed within the accretion zone onto virtual cloud
particles (each cell in the accretion zone is divided into $8^{3}$ cloud
particles). Additional checks can be then performed to avoid to accrete gas
that has for instance too much angular momentum to be gravitationally bound.
When a parcel of gas is accreted, it also transfers its linear and angular
momentum to the sink. An absolute mass cap is also used to limit to 25% the
amount of mass accreted from a unique cell in a single time step.

After accretion, the position of the sink particle is changed in two steps
according to its momentum after gas accretion and through gravity. First the
momentum is updated to account for the gravitational interaction with the gas.
The total force on a sink particle is computed by a direct summation over all
the cells in the computational domain, except the ones in the accretion zone.
This method remains practicable as long as the number of sink particles
created remains limited. For the accretion zone, the interaction is computed
between the sink and the cloud particles created during the accretion step. A
Plummer law is used to soften the gravitational interaction at short distance,
i.e., $F_{\mathrm{grav}}\sim 1/(r^{2}+\epsilon^{2})$ where $\epsilon=2\Delta
x_{\mathrm{min}}$ is the softening length. Second, the acceleration due to
particle-particle interaction is added to the sink particle momentum using a
Bulirsch-Stoer algorithm with adaptive time step and error control. Note that
the stability of the temporal discretisation scheme required that the sink
particle velocity satisfies a CFL-like condition ($v_{\mathrm{sink}}\Delta
t<C\Delta x$, with $C<1$).

Sink particles algorithms have been complemented in the past ten years to
improve their versatility, in particular in the transition between the sub-
and super-sonic regimes. All newer implementations have been built upon these
two pioneering works.

###  7.2 Contemporary implementations

In all current implementations, the sink creation and accretion are the most
sensitive algorithms. Of course, the sink particles can be introduced in the
initial conditions. Otherwise, they are created after the creation checks.

Each of the commonly used codes in star formation has its own implementation
of sink particles. For the grid-based codes, we note the work of Wang et al.
(2010) in the ENZO code, which is very similar to Krumholz et al. (2004) but
with a different merging scheme which involves a maximum mass sink mass for
merging. Federrath et al. (2010) presented an implementation in the AMR code
FLASH, where the creation checks were similar to the ones by Bate et al.
(1995) for SPH with two additional checks: the cell tagged for sink creation
should be at the maximum level of refinement, and it should sit in a local
gravitational potential minimum. They also account for magnetic energy in the
energy budget done in the check algorithm. The checks are performed on a
spherical region of radius $r_{\mathrm{acc}}$ around the cell candidate for
sink creation. Their scheme for accretion is also different from the previous
grid-based works. If a cell $i$ within $r_{\mathrm{acc}}$ exceeds the density
threshold $\rho_{\mathrm{res}}$ for sink creation, the mass increment $\Delta
M=(\rho_{i}-\rho_{\mathrm{res}})\Delta x_{i}^{3}$ is considered for accretion
onto the sink particle. The mass increment has to be gravitationally bound to
the sink, and in the case of overlapping sink particles, it is accreted to the
most gravitationally bound sink particle, the latter being moved to the center
of mass of the particle-gas configuration before the accretion step. They also
improve the treatment of the gas-sink gravitational interaction. Similarly to
Krumholz et al. (2004), the sink-sink interaction is done by direct summation
over all the sink particles. The gravitational acceleration of the sink due to
the gas is estimated from the gravitational potential of the gas which is
handled by the Poisson solver of FLASH to handle the gas-gas interaction. The
acceleration is interpolated using a first-order cloud-in-cell method for each
sink particle. Last, the acceleration of the gas due to the sink is done by a
direct summation. Each sink particle contributes to the acceleration of each
cell in the computational domain. The direct summations for the sink-sink and
sink-gas interactions requires a gravitational softening. Federrath et al.
(2010) use a cubic spline softening similar to what is used in SPH, which is
less aggressive within the softening length than a simple Plummer profile.
Federrath et al. (2011) updated their first implementation to ensure exact
momentum conservation using a direct summation between all cells and sink
particles for the gas-sink interaction as well.

In addition, Federrath et al. (2010) performed a quantitative comparison
between their implementation and the one by Bate et al. (1995) by comparing
the results of two star formation calculations performed with FLASH and with a
SPH code. They find a good quantitative agreement which was quite encouraging
given the fundamental differences between SPH and AMR sink implementations.

Hubber et al. (2013) presented an improved algorithm for the sink particles in
SPH, which can handle situations where sinks are introduced after the gas has
become adiabatic. They propose a different combination of creation criteria
which are based on the studies we just mentioned: i) density threshold, ii) no
overlapping with another sink accretion volume, iii) the gaseous particle
flagged for sink creation sits in a minimum of gravitational potential (on a
volume defined by its neighbours), iv) the gaseous particle is outside the
Hill sphere of fragments harbouring another already formed sink.

Last, Bleuler and Teyssier (2014) presented an innovative sink creation scheme
based on a clump finder algorithm used on-the-fly in the RAMSES code. The
clump finder is performed in a first step to identify peaks and their
associated regions. Second, the peaks are considered for sink creation. A
virial theorem type analysis is performed on the clumps, accounting for
surface pressure as well as tidal forces. The virial check is passed if the
second derivative of the moment of inertia in the center of mass frame is
negative. Then a collapse check is performed, which corresponds to the
contraction check of Krumholz et al. (2004) and Federrath et al. (2010), but
adapted according to the virial analysis. Last, the final check is a proximity
check. We also note that Bleuler and Teyssier (2014) offer also the
possibility to use the alternative sink creation and accretion schemes
reported in the literature. They provided an excellent review of the different
creation checks, accretion, sink merging and trajectories schemes. In
particular, they tested different schemes to perform sink accretion: Bondi-
Hoyle accretion, flux accretion, and density threshold accretion. They show
that all methods give good results for spherical Bondi accretion (i.e., when
the flow is supersonic). In the subsonic regime, only the Bondi formula gives
a correct accretion rate. They recommend to automatically adapt the accretion
scheme depending on the surrounding flow properties.

A number of other sink particle implementations have been reported in the
literature that we cannot detail in this review. We give in the following a
non-exhaustive list, restricted to the field of star formation. In SPH, sink
particles have been implemented in the codes GADGET (Jappsen et al., 2005),
DRAGON (Goodwin et al., 2004), and SEREN (Hubber et al., 2011). For grid-based
codes, sink particle are used in the uniform grid code ATHENA by Gong and
Ostriker (2013), in RAMSES by Padoan and Nordlund (2011), and in the ORION2
code by Lee et al. (2014). Last, Greif et al. (2011) also proposed a sink
particle algorithm within the AREPO moving-mesh code. Their implementation is
very similar to the one used in SPH codes. We also note that sink particle
merging is now widely used in most of the implementations.

Twenty years after their introduction, sink particles remain a sensitive
subject in the community. It is commonly accepted that an accreting sink
particle should have a minimal impact on the collapsing gas around it. The
easiest way to deal with it is to create the sink particle during the (quasi-)
isothermal phases of the collapse (either the first or the second collapse),
when the gas velocity is supersonic. The perturbation created by the sink
particles on the gas dynamics cannot propagate since the CFL timestep is
limited by the supersonic motions. An example of bad behaviours observed in
protostellar collapse SPH calculations using sink particles introduced in the
vicinity of the first cores can be found in Appendix D of Commerçon et al.
(2008). When a barotropic law is used to mimic the gas thermal budget, a
simple method used to introduce sink particles in the vicinity of hydrostatic
object can be found in Price and Bate (2008). They artificially change the gas
equation of state to isothermal at an arbitrary density (e.g. $10^{-11}$ g
cm-3) and set the density threshold for sink creation to a value of more than
two orders of magnitudes larger.

###  7.3 Sink particles and magnetic fields

In the previous section, we did not discuss the sink particle algorithms in
the presence of magnetic fields. We try to describe here some attempts done to
provide a good description of MHD flows for sink creation and accretion. Sink
particles can of course be introduced in MHD calculations, but the
implementation is not as straightforward as in the hydrodynamical case. While
it is relatively easy to accrete gas and momentum, magnetic flux accretion
cannot bypass the solenoidal constraint $\nabla\cdot\vec{B}=0$.

Price and Bate (2008) and Wang et al. (2010) performed star formation
calculations in turbulent and magnetized clouds, but do not account for any
contribution of the magnetic fields in their sink algorithms. Federrath et al.
(2010) have introduced the magnetic energy in their negative total energy
check for sink creation, but not in the accretion scheme. Lee et al. (2014)
have extended the work of Krumholz et al. (2004) to deal with magnetized
flows. For sink creation, they follow the work of Myers et al. (2013) who
derived a magnetic Truelove criterion. They derive simple steady state
accretion rates as a function of the plasma beta $\beta$ and the Mach number
$\mathcal{M}$.

In all implementations, gas is accreted by the sink particles but not magnetic
flux. It stays the gas in the surrounding of the sink particles, and magnetic
flux accumulates as collapse goes on. While at scales of a few 10-100s AU this
flux accumulation is not problematic (Federrath and Klessen, 2012; Li et al.,
2018), it eventually leads to a magnetic explosion and magnetic flux
redistribution in the surrounding of the sink (Zhao et al., 2011) due to the
development of an interchange instability (Spruit and Taam, 1990; Spruit et
al., 1995; Li and McKee, 1996), which is also observed in high resolution
models without sink particles and ideal MHD (Li et al., 2014; Masson et al.,
2016). In addition, since density stays roughly constant while magnetic
intensity increases in the sink accretion zone, the Alfvén speed increases
significantly, leading to non-physical accretion rates, as well as short
integration time steps. Lee et al. (2014) proposes to cap the accretion rate
at a maximum value corresponding to a constant Alfvén speed after mass
accretion. This trick tempers the accumulation of magnetic flux as well.

In the recent years, resistive MHD calculations have shown that all non-ideal
effects leads to a decrease of the magnetic flux at the first core scale.
Magnetic fields decouple from the gas and dust mixture dynamics. The magnetic
flux does not accumulate at the first core border, but rather at a larger
distance which corresponds to the radius at which non-ideal effects start to
dominate (Masson et al., 2016; Hennebelle et al., 2016). Matter is accreted
but leaves its associated magnetic flux outside the decoupling region. Adding
sink particles in restive MHD calculations with resolution capable of
resolving the decoupling region has not yet been investigated in detail
(Machida et al., 2009; Matsumoto et al., 2017; Tomida et al., 2017), but one
can foresee that introducing a sink within the decoupling region should
prevent magnetic flux accumulation. It is expected that non-ideal MHD can
alleviate the conditions of interchange instability, but is is then necessary
to resolve the decoupling scale, i.e., the disk. To date, the most detailed
study of sink particles introduction in MHD calculations can be found in
Machida et al. (2014).

###  7.4 Sink particles versus sink cells

A fundamental difference between sink cells and sink particles is that sink
cells are fixed in position and have physical boundaries. There are pros and
cons for the two methods, which we briefly summarize in the following.

First, the sink cell can be considered as a boundary in codes using a
spherical grid. In collapse calculations, a sink cell is most often introduced
as a passive boundary with zero gradients, or also called outflowing boundary
conditions (Yorke et al., 2002; Kuiper et al., 2010b; Li et al., 2011). The
mass flux across the boundary is added to the central point mass of the sink
cell and used to compute the gravitational acceleration. The sink cell enables
one to make more detailed models for the magnetic field configuration at the
boundary, as for instance in star-disk interaction studies, (e.g. Mellon and
Li, 2008; Takasao et al., 2018). In principle, this is a quite powerful
numerical set-up to overstep the issue of time integration after the second
core formation. The current state-of-the-art is not yet mature to design such
numerical studies, but this should certainly be the purpose of future works.

On the other hand, a major inconvenience of sink cells resides in their fixed
position in time, which forces the center of mass to be equal to the
geometrical center of the disk. As a consequence, the development of non-
axisymmetric azimuthal low wave-number modes may be artificially damped, which
alters the development of eccentric gravitational instabilities (Adams et al.,
1989; Shu et al., 1990). This effect can be particularly problematic in the
case of young stellar objects with disks of mass comparable to the stellar
mass, $M_{\mathrm{disk}}/M_{\mathrm{star}}>1/3$ (Krumholz et al., 2007; Kuiper
et al., 2011; Sigalotti et al., 2018).

Last, the use of ”outflow” boundary conditions for the mass accretion on the
sink cell is very different from the algorithms we just presented for the sink
particles. Machida et al. (2014) have shown that applying similar criteria for
the sink accretion than the ones employed in studies with a sink cell (Li et
al., 2011) may lead to very different results, i.e., no disk formation,
compared to the same calculations done with a classical sink particle
algorithm, i.e., formation of a disk.

##  8 Sub-grid models for sink particles

Beyond the practical advantage of saving computational time, the sink particle
has been shown to be very useful in works where (proto-)stellar feedback is
considered. Similar to cosmological and galactic evolution studies, the sub-
grid models can be associated to the sink particle evolution. We distinguish
two types of sub-grid physics, depending on the feedback origin: radiative
(luminosity, accretion shock, ionisation), and dynamical (jets/outflow,
winds). In any case, the source terms introduced by sub-grid physics are
treated using operator splitting schemes.

###  8.1 Radiative feedback

The first class of sub-grid models were designed to account for the radiative
feedback of the nascent protostars. The radiative feedback term includes the
accretion, the internal, and the ionisation luminosities. When integrating
over discrete timesteps, the luminosity is an energy or radiative flux input.
Radiative feedback was introduced early in 2D models (e.g., Bodenheimer et
al., 1990). They used a 2D axisymmetric hydrodynamical code (Rozyczka, 1985),
which was primarily designed for stellar winds. Bodenheimer et al. (1990) used
the grey FLD approximation for the radiation transport. They combined the
evolution of their central zone, or sink cell, with a pre-main sequence
evolution model based on Maclaurin spheroids in hydrostatic equilibrium. They
computed the structure of the protostar according to the central $(\rho,T)$
path tracks reported by Winkler and Newman (1980) from 1D spherical
calculations. The protostar mass $M_{\star}$ and equatorial radius $R_{\star}$
were then used to compute the accretion luminosity
$L_{\mathrm{acc}}=0.5GM_{\star}\dot{M}/R_{\star}$, where $\dot{M}$ is the mass
flux crossing the inner boundary. Bodenheimer et al. (1990) further assumed
that the radiative energy input $L\Delta t$ was only due to the accretion
luminosity and that all the infall kinetic energy was converted into radiation
at the stellar core accretion shock. They found that most of the energy that
went for heating in the envelope comes from this protostellar luminosity.

Yorke and Bodenheimer (1999) and Yorke et al. (2002) used a more sophisticated
protostellar evolution model which accounts for the intrinsic luminosity of
the protostar $L_{\mathrm{int}}$. They realised that the central luminosity
evolution sets the thermal budget within the envelop, but could not resolve
the innermost regions of the star-disk system. They designed a protostellar
accretion luminosity model which accounts for the luminosity coming from
Keplerian motion of the disk being converted into heat and radiated away at
the disk-core boundary layer following (Adams and Shu, 1986)

$L_{\mathrm{tot}}=L_{\mathrm{int}}+\frac{3GM_{\star}\dot{M}}{4R_{\star}},$
(64)

where they assume that a fourth of the total potential energy of the accreted
material is released into heat within the disk, the remaining being radiated
away within the unresolved inner disk region and at the accretion shock. In
addition, Yorke and Bodenheimer (1999) and Yorke et al. (2002) followed
radiation transport using a frequency dependent ray-tracing algorithm (with 64
frequency bins). In particular, Yorke et al. (2002) demonstrated the
importance to handle the protostellar irradiation using a frequency dependent
irradiation scheme in the context of massive star formation. The flashlight
effect, i.e., radiation escapes in the vertical direction and the radiative
acceleration is reduced within the disk to allow accretion, is enhanced
compared to a simple grey model. In addition, Krumholz et al. (2005) showed
how protostellar outflows allow radiation to escape in the outflow cavities to
allow a continued accretion onto the central massive stars. These results have
been further confirmed by 3D dynamical simulations of several authors
(Krumholz et al., 2007; Kuiper et al., 2011; Rosen et al., 2016; Klassen et
al., 2016; Kuiper et al., 2016).

Currently, most 3D numerical RHD calculations include stellar radiative
feedback sub-grid models attached to sink particles/cells evolution. Most of
the implementations account for both the internal and the accretion
luminosities, with some modulation factors on the amount of potential energy
radiated away and on the efficiency of the mass transfer for the sink
accretion radius to the protostar. The simpler sub-grid models use a grey FLD
model for the radiation transport where the protostellar luminosity is simply
a source term in the radiative energy equation (Krumholz et al., 2007; Offner
et al., 2009; Stamatellos et al., 2011; Fontani et al., 2018; Jones and Bate,
2018). Frequency dependent protostellar irradiation modules have been
developed using ray-tracing to compute the radiation fields, combined with
moment models for radiation hydrodynamics (Kuiper et al., 2010b; Wise and
Abel, 2011; Klassen et al., 2014; Ramsey and Dullemond, 2015; Buntemeyer et
al., 2016; Rosen et al., 2017). These frequency dependent modules have been
used primarily in the context of massive star formation (Rosen et al., 2016;
Klassen et al., 2016).

Last, we note the developments made to capture ionising radiation feedback
from massive stars: Dale et al. (2007) used a Strömgren volume method in SPH,
Peters et al. (2010) ray tracing in the FLASH code, Kuiper and Hosokawa (2018)
an hybrid ray-tracing and FLD irradiation module in the PLUTO code, Geen et
al. (2015) the M1 moment method in the RAMSES code, and Harries (2015) Monte
Carlo radiative transfer in the AMR code TORUS. These schemes are applied
first in isolated collapse calculations (Kuiper and Hosokawa, 2018) and,
mostly, in cluster formation studies (Peters et al., 2011; Dale and Bonnell,
2011; Dale et al., 2012; Geen et al., 2015; Gavagnin et al., 2017; Ali et al.,
2018; Geen et al., 2018).

###  8.2 Dynamical feedback

The most advanced MHD collapse calculations have shown that outflows and jets
naturally develop at different scales during the collapse (Tomisaka, 2002;
Banerjee and Pudritz, 2006; Machida et al., 2006; Hennebelle and Fromang,
2008; Ciardi and Hennebelle, 2010; Commerçon et al., 2010; Tomida et al.,
2010; Price et al., 2012). In addition, outflows and jets are commonly
observed in YSO’s (it is for instance a selection criterion for Class 0
sources, Andre et al. (1993)) and they are considered as possible sources of
turbulence driving in star forming clouds and participate in the regulation of
the star formation rate (Matzner, 2007; Nakamura and Li, 2007; Krumholz, 2014;
Federrath, 2015). Outflows and jets launching scales are nevertheless not
captured in most studies because it requires a very high numerical resolution
to reach sub-AU scales. Besides, studies neglecting magnetic fields cannot
generate MHD winds driven centrifugally or by toroidal magnetic pressure.
Consequently, very few numerical calculations have been able to launch self-
consistently MHD outflows (e.g., Hennebelle et al., 2011). Sub-grid models
have thus been developed to account for outflows/jets in numerical works with
unresolved launching scales and missing physics, similarly to what has just
been presented for the radiative feedback. The numerical implementation are
based on analytical works dedicated to MHD wind launching (Blandford and
Payne, 1982; Pelletier and Pudritz, 1992; Shu et al., 1994; Ferreira, 1997;
Matzner and McKee, 1999).

The first works on sub-grid models have been reported in Li and Nakamura
(2006) and Nakamura and Li (2007) who employed ideal MHD and self-gravity.
They attached a protostellar outflow model on sink particles following Matzner
and McKee (2000) because they were unable to describe the self-consistent
launch of outflows with the crude resolution they use ($128^{3}$ for a
$\approx 2$ pc box). Each star injects in the ambient medium a momentum that
is proportional to the stellar mass $M_{\star}$ with a two-component outflow
model. The volume of injection is given by the sink accretion volume (which
they refer to as ”supercell”). A fraction $\eta$ of the total outflow momentum
is put in a conical collimated jet component to facilitate the transport of
energy and momentum at large distances. The rest of the momentum is put into a
spherical component around the protostar. The direction of the jet is given by
the orientation of the magnetic field in the central cell, and they assume an
opening angle of 30∘ about this direction.

Most of the implementation of outflows and jets have a similar construction,
with flavors depending on the sink particle algorithm. Wang et al. (2010)
implemented in the MHD ENZO code a simplified version for protostellar
outflows where they only account for the collimated jet component. They adopt