Numerical Methods for Simulating Star Formation

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-Alfvenic 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 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 (abridged).


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.

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.

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, ∂ρ ∂t + ∇ · (ρv) = 0, where ρ is the mass density and v is the fluid velocity vector. We also have the momentum conservation equation where B is the magnetic field vector and P tot is the total pressure, defined as the sum of the thermal pressure and the magnetic pressure Note that we work here in cgs units, hence the presence of the 4π term in these equations. We have also included the gravitational acceleration vector g as a source term on the right-hand side of the momentum conservation equation. Finally, we have the total energy conservation equation where the total fluid energy is the sum of the kinetic energy, the internal energy and magnetic energy In order to close the system, we need the fluid equation of state, usually given be the ideal gas equation of state P = (γ − 1)e, and the induction equation for the evolution of the magnetic field where we have introduced the electric field E (also known as the electromotive force, emf). We need to add to these equations the most important property of the magnetic field, namely 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.

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 where the vector 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 Similarly, the divergence-free condition, written in integral form as 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 where the operator [A] = A + − A − denotes the jump of a quantity A across the surface element.

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, piecewiseconstant, 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 ∆φ = ∇ · B * .
where φ is a scalar potential. The new magnetic field is then corrected as 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 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 ∇ · B is advected away Frontiers 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 points 1 , 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 ψ that satisfies and is used as a source term in the induction equation.
In the first equation, c h is the divergence wave speed and τ 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 ψ 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). 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.

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;, 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) If λ J = c s √ 4πGρ < 4∆x then refine to level + 1.
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 The ψ i (x) are usually an orthonormal basis of polynomials of degree at most p. The coefficients α i are computed using two different philosophies. The first option, the WENO approach computes the coefficients α 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 α 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 α i , while in the DG case, one needs to store permanently the α 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: 1the 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.

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;. 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); Tomida et al. (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. 3.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 where J = ∇ × B is the current, η Ω , η H , and η AD are the Ohmic, Hall, and ambipolar resistivities, respectively (in units of cm 2 s −1 ), and v is the velocity of the neutrals. The notation ||B|| represents the norm of the magnetic field vector 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 η Ω ∆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 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 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.

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 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 STS sub-timesteps ∆t STS which are based on Chebyshev polynomials to guarantee the stability conditions over the super-timestep ∆t tot = N STS 1 ∆t 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  and anisotropic diffusion (Meyer et al., 2012;Vaidya et al., 2017).

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 where ρ i (resp. v i ) and ρ n (v) are the mass density (resp. velocity) of ions and neutrals, and γ AD is the collisional coupling constant, with γ AD ∼ 3 × 10 13 cm 3 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 The induction equation thus reads From Eq. 22, we have η AD = ||B|| 2 /ρ i ρ n γ AD . In the dense ISM, the heavy charged particles are ions and charged dust grains, and the expression of the resistivity is more complex . 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  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. 3.6).

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 where we introduce the "Hall speed" v H = − η H ||B|| 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) The whistler waves are dispersive, the phase speed c w = ω/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 The time step is then constrained as 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 n+1 x of the magnetic fields is given by B n y and B n z at time n, explicit-implicit for the second dimension, B n+1 y is given by B n+1 x and B n z , and finally implicit for third dimension, B n+1 z is given by B n+1 x and B n+1 y . 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 ∇ · 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 The current J = ∇ × 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.

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 Marchand et al., 2018), the SPH code by Tsukamoto et al. (2017), PHANTOM , ZeusTW , 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).

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(Balsara, , 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.

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 g , the force is where v g is the dust grain velocity, v the gas velocity, and t 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 s ∼ s (Epstein, 1924). The Stokes number St ≡ t s /t dyn characterizes the coupling of the dust grain relative to the gas, where well coupled dust grains follow St < 1. t 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.
Several implementations of dust and gas mixtures have been developed in the literature to deal with the various Stokes regime. For 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 St ≤ 1. This monofluid formalism has been implemented with success in SPH and 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.

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 Requiring that the Jeans length of the gas becomes optically thick to infrared radiation leads to Many MHD simulations use this critical density to define a polytropic equation of state of the form 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;. In conclusion, radiation plays a central role during the formation and the evolution of the first Larson core (see Sec. 5), but later on acts as a self-regulating mechanism for star formation.

Ray tracing and long characteristic methods
The radiative transfer equation written in full generality using the radiation specific intensity I ν (x, n, t) bf where j ν is the emission coefficient and α ν 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 ds = cdt. The equation of radiative transfer can be written as a Lagrangian time derivative following the trajectory of the photons as 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).

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 .
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).

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 where E ν is the radiation energy, F ν is the radiation flux and P = DE ν is the radiation pressure tensor and D is the Eddington tensor. For slowly varying spectral absorption and emission coefficients, for which one has ν ∂α ν ∂ν α ν and ν ∂j ν ∂ν j ν , one can write relativistic correction to first order in (v/c) as 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 Since E ν and F ν are still expressed in the laboratory frame but α 0 ν and j 0 ν are expressed in the comoving frame, these equations are often called mixed frame radiation moments equations. We can also Lorentztransform the radiation moments from the laboratory frame to the comoving frame, which gives, to first order in (v/c) Injecting these relations only in the right-hand side of the moments equations and neglecting all (v/c) terms gives the simpler form 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 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 M1 closure : 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 D n ⊗ n. The unit vector n and the parameter χ are determined using the radiation energy and the radiation flux as 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.

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 D 1 3 I and the radiation flux can be approximated by where λ(R) is the flux limiter, a function that has to connect the diffusion limit with λ 1/3 for R 0 to the free-streaming regime where λ 1/R for R → +∞. A possible form for the flux limiter has been proposed by Minerbo (1978) and reads while the Levermore (1984) flux limiter reads 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.

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 R form within dense cores of sizes 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.

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 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 H 2 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 10 −8 g cm −3 until stellar densities are reached within 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 10 AU, a mass 0.01 M , and a lifetime 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 H 2 dissociation. Masunaga and Inutsuka (2000) integrated their models throughout the Class 0 and Class I phases up to an age of 1.3 × 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(Vaytet et al., , 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).

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. 2, 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. 4). 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 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 γ which varies with the temperature. At low temperature (T 150 K), the H 2 molecule behaves like a monoatomic gas (γ = 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 γ = 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 H 2 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. (2013Vaytet et al. ( , 2014; Vaytet and Haugbølle (2017); Vaytet et al. (2018), and also in the SPH public code PHANTOM .
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,C V ≡ 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 M 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 × 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;Vaytet et al., 2013Vaytet et al., , 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;. 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); Tomida et al. (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. 3. 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 ∆x 6.6 × 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; and AMR (Vaytet et al., 2018). The typical resolution in SPH calculations currently achievable with limited CPU resources is of 3 × 10 6 equal-mass particles per M . 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  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  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.

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.6.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.

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 acc of the sink particle. This ensures that the sink particle is formed from at least N neigh . r 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 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 α and β of the thermal energy and rotational energy to the gravitational energy must satisfy Additionally, the total energy must be negative, and finally, the system must be contracting (∇ · 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 ∼ N 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 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 ρ exceeds the maximum Jeans density ρ J that can be resolved where J < 0.25 is the Jeans length resolution criterion, c s the gas sound speed, and ∆x min the minimum cell size. The initial mass of the sink is m sink = (ρ − ρ J )∆x 3 min . 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 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 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 acc = 4∆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 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 grav ∼ 1/(r 2 + 2 ) where = 2∆x 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 sink ∆t < C∆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.

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 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 acc exceeds the density threshold ρ res for sink creation, the mass increment ∆M = (ρ i − ρ res )∆x 3 i 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 . 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.

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 ∇ · B = 0. (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  who derived a magnetic Truelove criterion. They derive simple steady state accretion rates as a function of the plasma beta β and the Mach number M.

Price and Bate
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 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).

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 disk /M 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  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.

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.

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 (ρ, T ) path tracks reported by Winkler and Newman (1980) from 1D spherical calculations. The protostar mass M and equatorial radius R were then used to compute the accretion luminosity L acc = 0.5GM Ṁ /R , whereṀ is the mass flux crossing the inner boundary. Bodenheimer et al. (1990) further assumed that the radiative energy input L∆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 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) 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,  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;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;Ali et al., 2018;Geen et al., 2018).

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., . 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 ≈ 2 pc box). Each star injects in the ambient medium a momentum that is proportional to the stellar mass M 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 η 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 a continuous momentum injection ∆P = P ∆M where P is a proportionality constant which depends on the stellar mass ∼ M 1/2 ). If the momentum injection occurs in a cell with a too low density, they take 10% of ∆M out of the sink particle and redistribute it evenly between the injection cells in order to avoid too large outflow speeds (and too small integration timestep). Along the same line, Cunningham et al. (2011);Hansen et al. (2012) extended the work of Krumholz et al. ( , 2007 and presented for the first time a combined sub-grid model which includes protostellar outflows (Matzner and McKee, 1999) and protostellar radiative heating , ignoring magnetic fields though. Their sub-grid model accounts for pre-main sequence evolution as in  and the calculated protostellar radius is used in the outflow and radiative heating models. They parameterized their outflow model with dimensionless parameters which, for instance, sets the fraction of mass accreted by the protostar or launched in the outflow, and the ejection velocity as a function of the Kepler speed at the surface. Contrary to the case of radiative heating where energy is put within the accretion volume, they inject outflows at a distance comprised between 4∆x and 8∆x from the sink particle, and used angular dependency from Matzner and McKee (1999). We note that they not only inject momentum, but also mass and thermal energy. They set the wind temperature to 10 4 K, which is appropriate for an ionised wind. Using the same tool, Myers et al. (2014) added magnetic field in their simulations with ORION, but because of the low resolution they used (> 20 AU), they where not able to launch self-consistent outflows and used the sub-grid model we just mentioned. The outflows they produce may form strong shocks at very high temperature, higher than the dust sublimation temperature so that the dust opacity drops. Their radiative transfer scheme, primarily designed for dust thermal emission does not allow the gas to cool efficiently. They thus needed to add a line cooling function for the thermal budget of the gas.
Federrath et al. (2014) implemented a sub-grid-scale outflow model upon the sink particles algorithm they implemented in FLASH (Federrath et al., 2010) for MHD collapse. Their model combines different features of the Li and Nakamura (2006); Nakamura and Li (2007) and Cunningham et al. (2011) implementations. They used a normalised velocity profile for momentum injection, which reproduces the two components outflow/jet with opening angles of 30 • and 5 • respectively. They also improved the algorithm for the outflow orientation by recording the angular momentum transfer from the accreted gas to the sink particle. Their algorithm also conserves mass and momentum exactly in the sink accretion step. Murray et al. (2018) proposed a similar implementation in RAMSES. They follow previous works (Cunningham et al., 2011;Myers et al., 2014) for the outflow injection as well as the protostellar evolution. Outflows are there injected within a conical volume about the spin axis of the sink particle in a radial extent comprised between 4 and 8 cells (at the highest refinement level) away from the sink particle.  performed a parameter study on the number of cells for the radial outflow direction and found that convergence on the maximum outflow velocity requires 16 cells per outflow radius. Convergence on the mass, linear and angular momentum of the outflow is achieved though with 8 cells.
In the context of massive star formation, another class of outflows models are based on the luminosity and ionisation. Sub-grids models of outflows generated by massive star luminosity have been developed (e.g. Peters et al., 2014), but the physical mechanisms driving this type of outflows are not yet well understood. No work in the context of massive star formation has accounted for the combined effects of magnetic fields and radiative force with a resolution sufficient to resolve the magnetic outflow launch. There is certainly a combination of magneto-centrifugal process and from radiative pressure, depending on the mass of the forming protostar.

Second collapse and pre-main sequences
Ideally, the sub-grid models should be designed to reproduce the results of high resolution studies, which resolve the protostar scales. The short time integration after stellar core formation, dominated by the adiabatic evolution, does not allow to design sub-grid models for protostellar feedback from full 3D RMHD calculations with a high level of confidence. For these reason, all recent works mostly rely on the results of 1D spherical symmetry results to set their protostellar core evolution models. For instance, Kuiper and Yorke (2013) use the pre-main sequence evolution code STELLAR (Bodenheimer et al., 2007;Yorke and Bodenheimer, 2008) to determine the protostellar luminosity of the non-resolved protostars. Pelupessy et al. (2013) developed the open source Astrophysical Multi-purpose Software Environment (AMUSE), a component library of simulations involving different physical domains and scales. For instance, it enables one to couple hydrodynamical codes (GADGET, ATHENA, AMRVAC) with pre-main sequence and binary evolution codes such as MESA (Paxton et al., 2011). Recently, Wall et al. (2019) used AMUSE to perform a star cluster formation model that includes individual star formation from self-gravitating, magnetized gas, coupled to collisional stellar dynamics. It couples different tools: the AMR code FLASH (Fryxell et al., 2000), the N-body code ph4 (McMillan et al., 2012), and the stellar evolution code SeBa (Portegies Zwart and Verbunt, 1996). We anticipate to see more works using this strategy of coupling tools and scales in the near future.
Recent work from Vaytet et al. (2018) have managed to resolve the second core formation with nonideal MHD (Ohm diffusion and ambiploar diffusion) from collapse 1 M dense core (see Figure 3). Interestingly, they show that the magnetic field lines topology is very different at the first and second core scales. Ambipolar diffusion is dominating at the first core scale. The magnetic fields direction is essentially vertical and the magnetic fields and the gas evolution are decoupled. At the second core scale, the ionisation rises due to high temperature and the system is back to ideal MHD, but with a high plasma beta β = P therm /P mag . The weak magnetic field is efficiently wrapped up by the gas rotation and the resulting fields direction is essentially toroidal. On the opposite, when ideal MHD is preserved all the way from the dense core to the stellar core, the magnetic fields lines are strongly pinched and the field remains poloidal at the second core scale. The energy balance at the first and second core accretion shocks is also different. While the accretion shock at the first core has been classified as a radiative supercritical shock, with all the incident kinetic energy radiated away (or often referred to as cold accretion), in 1D spherical core collapse studies (Commerçon et al., 2011a;Vaytet et al., 2012)), recent 3D works have shown that the first core accretion experiences both cold and hot accretion at the same time at its surface, depending on the accretion flows morphology. On the surface of the stellar core, the accretion shock has been found to be subcritical (Vaytet et al., 2012;Tomida et al., 2013;Vaytet et al., 2018) with a significant part of the incident kinetic energy transferred to the protostars. Although these results have been investigated only for the very first stages of the protostar evolution (only one month in Vaytet et al. (2018)), they indicate that it is not easy to derive general properties of small scales phenomena to design robust sub-grid models. In particular,  show that the magnetic fields evolution of the protostars they formed in their SPH models is largely affected by numerical diffusion. Last, the radiative efficiency of the stellar accretion shock has been shown to be a key process for the pre-main sequence evolution and sets the radius of the young protostars (Baraffe et al., 2012). We recall that the sub-grids models designed for the radiative and dynamical feedback take the protostellar radius as an input. The launching speed of the protostellar jets as well as the accretion luminosity depends on R −1 .

DISCUSSION
We have described almost all the numerical methods that are required to model the star formation process within the turbulent interstellar medium. These discretised equations must now be solved for a given set of initial and boundary conditions, over some prescribed time. This defines the numerical set up and the overall simulation strategy. Obviously, given the wide range of spatial and temporal scales involved, one cannot realistically model an entire galaxy all the way down to the formation of brown dwarves and proto-planetary disks. Unfortunately, given the complexity of the star formation process, and the fact that large and small scales are strongly coupled, this is required by the physics of the problem. Indeed, supersonic turbulence is seeded on large scales, where kinetic energy is injected through shearing or colliding flows induced by galactic rotation, spiral waves and extra-galactic accretion flows. Moreover, turbulence and gas ejection can be triggered by stellar feedback processes occurring on very small scales, such as radiation driven flows, jets and ultimately supernovae explosions.
Modelling star formation is therefore fundamentally a multi-scale, multi-physics problem, and a very difficult one. Past or present day simulations always rely on some compromises. We can decompose them into several categories, corresponding roughly to the adopted box size and initial conditions. First, full galaxy simulations, where spiral waves and spiral shocks are used as the seed for turbulence. The box size must be of several tens of kiloparsecs to contain the entire rotating disk and possibly the galactic corona. The gas is compressed through the spiral shocks or various colliding flows. Thermal and gravitational instabilities fragment the gas into small molecular clouds. The spatial resolution never really drops below 0.1 or 1 pc in these simulations, so that stars cannot be modelled individually, at least for massive galaxies like the Milky Way (see Fig. 4 for illustration). Simulators rely on simplified, sub-grid star formation recipes based on stochastically spawning star cluster particles of the same mass, or on the sink particle technique where in this case sink particles represent a star cluster, rather than a single star. Bournaud et al. (2010) studied the properties the ISM substructure and turbulence in galaxy simulations with resolutions up to 0.8 pc and 5 × 10 3 M with RAMSES. (Renaud et al., 2013;Kraljic et al., 2014) were able to capture the transition from turbulence-supported to self-gravitating gas with resolution up to 0.05 pc in simulations of Milky Way-like galaxies using RAMSES. Hopkins et al. (2012) presented SPH simulations dwarfs galaxies and Milky Way (MW) analogues to massive star-forming galactic discs with pc-scale resolution. Later, Hopkins et al. (2014) wrapped their implementation of feedback in the FIRE (Feedback In Realistic Environments) simulations suite, aiming, among others objectives, at resolving the formation of giant molecular clouds and the multi-phase interstellar medium (ISM). More recently,  updated the FIRE implementation in the GIZMO code and performed sub-parsec resolution simulations of galactic discs. We also refer readers to the work of Dobbs et al. (2018) using SPH, Tasker (2011) using ENZO and Smith et al. (2018 using AREPO. In addition, we note that for dwarf galaxies, it becomes possible to model much smaller scales so that at least massive stars can be modelled individually: Hu et al. (2016); Hu (2019) using SPH and a resolution of 1 − 5 M by gas particle, Emerick et al. (2018Emerick et al. ( , 2019 using ENZO with a maximum resolution of 1.8 pc. The second approach relies on simulating only a small portion of vertically stratified galactic disks. The box size is usually around 1 kilo-parsec, with periodic boundary conditions in the direction of the disk plane and outflow conditions in the direction perpendicular to the disk plane. The geometry is clearly heavily constrained by this elongated, vertical and stratified layer of gas but it captures most of the phenomenon at work locally in the disk. Using the shearing box technique, one can also add more realistic shearing conditions, so that turbulence can be maintained both by stellar feedback and a large scale galactic shearing flow (Kim et al., 2002). In this case, the resolution can drop significantly below 0.1 parsec, may be 0.01 parsec or even less for the highest resolution simulations. Unfortunately, it is still impossible to model individual stars at this resolution so that here again most papers are based on the sink particle techniques for representing star clusters. It is however almost possible to resolve massive molecular cores (and their associated sink particles) individually, so that individual supernovae and HII regions can be resolved. These simulations are probably the most realistic models of the interstellar medium of the Galaxy, although still far from resolving the entire stellar population. In standard setups, the models account for the thermal instability and feedback (mechanical by SN and radiative by massive stars). First grid-based models were presented in Korpi et al. (1999) and included magnetic fields, SN heating and radiative cooling. 3D MHD AMR models of a kpc-scale ISM driven by SNe were presented in de Avillez and Breitschwerdt (2005) and Joung and Mac Low (2006). Currently, numerous studies and projects are based on a similar setup, with increasing physics put on over several years (SNe, radiative feedback, cosmic-rays, shear, etc...): the series of papers by Hennebelle and Iffrig (2014); Iffrig and Hennebelle (2017); Colling et al. (2018) as well as the work of Martizzi et al. (2016) and Butler et al. (2015Butler et al. ( , 2017 using RAMSES, by Kim et al. (2011Kim et al. ( , 2013; Kim and Ostriker (2015, 2018 using ATHENA, the SILCC project papers using FLASH (Walch et al., 2015;Girichidis et al., 2016;Gatto et al., 2017;Peters et al., 2017;Girichidis et al., 2018).
The third approach aims at simulating individual molecular cores. The box size ranges from 100 parsec down to 0.1 parsec for the highest resolution simulations. The mass typically ranges from 50 to a few 10 3 M . Boundary conditions are either periodic or isolated. In the latter case, the cloud as a whole can collapse, while the former set of simulations represents a more uniform background with collapsing regions only at very small scales. The smallest box sizes represents internal regions of larger clouds and are the only simulations for which the entire star population, down to the brown dwarf limit, can be modelled. These simulations are directly tackling the origin of the stellar IMF, and the formation of proto-planetary disks. We mention in the following a non-exhaustive list of the work done in this intense field of research. The series of papers by Bate & collaborators have investigated the collapse of molecular clouds of mass ranging from 50 to 500 M with increasing physics and numerical resolution since their pioneer work (Bate et al., 2003) using SPH. In particular, they have investigated the effect of radiative transfer (Bate, 2009(Bate, , 2012 and magnetic fields (Price and Bate, 2008;Price and Bate, 2009). The work of Dale et al. (2007); Dale and Bonnell (2011); Dale et al. (2012) have also used SPH simulations to study the effect of ionizing radiation from the forming massive stars. AMR codes have also been widely used in this context, accounting for a lot of physics and initial and boundary conditions. We note the work of Girichidis et al. (2011Girichidis et al. ( , 2012a using FLASH who studied the importance of isolated initial conditions in isothermal cloud core collapse (without stellar feedback).  also used FLASH to perform simulations of isolated cores with protostellar feedback (jets). Krumholz et al. (2007Krumholz et al. ( , 2012; Hansen et al. (2012); Myers et al. ( , 2014; Cunningham et al. (2011);Li et al. (2018); Offner and Chaban (2017) presented 3D collapse models using ORION with radiative transfer and/or ideal MHD, as well as protostellar feedback (luminosity, outflows). RAMSES has also been extensively used in this approach: Commerçon et al. (2011b) for 100 M isolated core collapse models with (R)MHD, Lee and Hennebelle (2018b,a,c) focused on the peak of the IMF using (M)HD 1000 M isolated core collapse models, Geen et al. (2015Geen et al. ( , 2016;  studied the effect of ionizing radiation using the M1 moment method for radiative transfer. Last, we note the work of Wang et al. (2010) using ENZO with isolated boundary conditions, ideal MHD and protostellar jets. This non-exhaustive list of works done using this approach demonstrates the importance and the utility of such models. These simulations, however, rarely follow the formation of the second Larson core. They still rely on the sink particle technique, and on sub-grid models as well, when it comes to modeling proto-stellar jet launching, energy budget at accretion shocks and their associated radiation.
The last category of simulations deals with isolated collapsing low-mass dense cores, with the goal of following the formation of protostellar disks and ultimately the entire second collapse without relying on any sub-grid model. This ambitious strategy has already been discussed at length in the previous sections, so we will not describe it here again. We also refer the reader to the review by Wurster and Li (2018) on the formation of protostellar disks. An interesting intermediate strategy has been developed recently with the zoom-in technique. This method is used now routinely in galaxy formation simulations, for which only a single dark matter halo can be re-simulated at much higher resolution, keeping the entire cosmological environment alive but at a coarser resolution. Here, the idea is to model the turbulent molecular cloud on large scales, say with a box of size 100 parsecs, to identify a single collapsing core, and to zoom-in on the core with much higher resolution and follow the first and possibly the second collapse and the formation of the proto-planetary disk. For instance, Kuffmeier et al. (2017); Kuffmeier et al. (2018) presented RAMSES AMR zoom-in calculations from an outer scale of 40 pc down to cell sizes of 2 au to study the effect of the environment on the formation and the evolution of protoplanetary disks. Zhang et al. (2018) performed high resolution calculations of the formation and evolution of a star-forming core, obtained by running larger scale calculations of molecular cloud formation . They used FLASH with a maximum resolution of 25 AU and were able to cover scales from 256 pc to 25 AU (see Fig. 5). Mocz et al. (2017) performed moving mesh AREPO calculations from a box size of about 5 pc down to a maximum resolution of 4 AU ( 5 × 10 −5 pc) to study the magnetic fields morphology from large to dense core scales. At larger scale, Padoan et al. (2017) performed a simulation of supernova-driven turbulence using RAMSES, with a box size of 250 pc and a maximum resolution of 0.0076 pc. At kpc scales, Hennebelle (2018) (FRIGG project with RAMSES, simulation of the formation of self-gravitating cores) and Seifried et al. (2017) (SLICC-Zoom project with FLASH, simulations of the formation of molecular clouds) presented zoom calculations of stratified Galactic disks down to resolution of 10 −2 − 10 −3 pc. Thanks to the steady increase in CPU power, we expect to see more and more work using the zoom-in technique in the coming years.
These different strategies are all in different ways very ambitious and address different problems in the theory of star formation. In addition, a compromise needs to be found between physical realism and computational efficiency. On one hand, we need to include magnetic fields, radiation fields and complex chemistry, but on the other hand, we need many small resolution elements and many time steps. Moreover, in the context of high performance computing, modern supercomputers are very hard to use at full efficiency when complex grid geometries is used in conjunction with expensive and demanding algorithms. There is a clear tendency for simulation projects deployed on the largest supercomputers in the world to use simple grid geometry, like Cartesian meshes and periodic boxes, with a highly simplified physical model. These large scale simulations are however very interesting to explore statistical aspects and large inertial range for turbulent flows. For example, Federrath (2013) performed 4096 3 hydrodynamic isothermal turbulence simulations using FLASH with periodic boundary condition for more than 40 000 time-steps on 32 768 CPU cores. In the framework of the magnetorotational turbulence, Fromang (2010); Ryan et al. (2017) performed high resolution simulations of isothermal shearing boxes using, respectively, ZEUS and a GPU version of RAMSES.
In conclusion, all simulations of star formation published in the literature, many of which have been discussed here, explore various corners of the numerical parameter space: resolution versus box size, statistics versus internal structure, physical realism versus computational speed. It is interesting that at the smallest scales of interest here, a "first principle" approach is in principle possible and pursued by several groups. At larger scales, however, star formation simulations share the same kind of limitations as galaxy formation simulations, or climate models and weather forecasting simulations, namely a strong dependence of the results on small, unresolved scales. Although bigger computers with more efficient, higher order codes and more realistic models will certainly help shed light of the mysteries of star formation, a robust methodology to implement sub-grid models and couple them properly to the fluid equations still needs to be invented in our field. Various attempts have been proposed in the context of unresolved turbulence and star formation from analytical works Hennebelle and Chabrier, 2011; Draine BT, Roberge WG, Dalgarno A. Magnetohydrodynamic shock waves in molecular clouds. ApJ 264