Sec. Stellar and Solar Physics
Numerical Methods for Simulating Star Formation
- 1Centre for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Zurich, Switzerland
- 2CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, Univ Lyon, Ens de Lyon, Univ Lyon1, Lyon, France
We review the numerical techniques for ideal and non-ideal magneto-hydrodynamics (MHD) used in the context of star formation simulations. We outline the specific challenges offered by modeling star forming environments, which are dominated by supersonic and super-Alfvénic turbulence in a radiative, self-gravitating fluid. These conditions are rather unique in physics and engineering and pose particularly severe restrictions on the robustness and accuracy of numerical codes. One striking aspect is the formation of collapsing fluid elements leading to the formation of singularities that represent point-like objects, namely the proto-stars. Although a few studies have attempted to resolve the formation of the first and second Larson's cores, resolution limitations force us to use sink particle techniques, with sub-grid models to compute the accretion rates of mass, momentum and energy, as well as their ejection rate due to radiation and jets from the proto-stars. We discuss the most popular discretisation techniques used in the community, namely smoothed particle hydrodynamics, finite difference and finite volume methods, stressing the importance to maintain a divergence-free magnetic field. We discuss how to estimate the truncation error of a given numerical scheme, and its importance in setting the magnitude of the numerical diffusion. This can have a strong impact on the outcome of these MHD simulations, where both viscosity and resistivity are implemented at the grid scale. We then present various numerical techniques to model non-ideal MHD effects, such as Ohmic and ambipolar diffusion, as well as the Hall effect. These important physical ingredients are posing strong challenges in term of resolution and time stepping. For the latter, several strategies are discussed to overcome the limitations due to prohibitively small time steps. An important aspect of star formation simulations is the radiation field. We discuss the current state-of-the-art, with a variety of techniques offering pros and cons in different conditions. Finally, we present more advanced strategies to mitigate the adverse effect of finite numerical resolution, which are very popular in the context of supersonic, self-gravitating fluids, namely adaptive mesh refinement, moving meshes, Smoothed Particle Hydrodynamics and high-order methods. Advances in these three directions are likely to trigger immense progress in the future of our field. We then illustrate the different aspects of this review by presenting recent results on supersonic MHD turbulence and magnetized collapse calculations.
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.
2. 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.
2.1. The Ideal MHD Equations
Before we review the various improvements of MHD numerical solvers over the past decades, we briefly recall the ideal MHD equations, shown here in conservative form. We have first the mass conservation equation,
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 Ptot 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
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.
2.2. Some Important Mathematical Properties
The fluid equations appear to be in conservative form, which naturally calls for a finite volume representation, which will ensure exact conservation of mass, momentum and total energy, by construction, owing to the divergence theorem
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.
2.3. Preserving the Divergence-Free Condition
Historically, one of the first ideal MHD, general purpose codes, is ZEUS-2D developed specifically for astrophysics by Stone et al. (1992b). It is a finite-difference code using constrained transport and artificial viscosity to handle shocks, The continuity and divergence-free conditions are therefore satisfied by construction, but since artificial viscosity is used to handle shocks, instabilities can occur in fast rarefaction waves (Falle, 2002).
The other popular strategy for grid-based methods is to maintain all MHD variables as volume-averaged quantities, allowing for discontinuities across all cell faces. In the late 90s, a series of papers presented finite-volume MHD codes with proper upwinding of numerical fluxes using Riemann solvers (Dai and Woodward, 1994; Ryu et al., 1995; Tóth, 1996; Balsara, 1998; Lee and Deane, 2009). These methods are very powerful, because they are strictly conservative and because they satisfy the so-called entropy condition, meaning the entropy can only increase, a key property to maintain the stability of the numerical solution. These finite-volume codes all considered the magnetic field as a volume-averaged, piecewise-constant, cell-centered quantity, which violates the continuity condition of the normal component of the field. The resulting schemes are therefore not necessarily divergence-free anymore. In order to maintain the divergence as small as possible, an additional step is required that modifies the magnetic field components and decreases or nullifies the divergence: this operation is called divergence cleaning. In a seminal paper, Tóth (2000) has compared various schemes and showed that they all passed with mixed success a variety of MHD test problems. We will now review these grid-based methods that are using divergence cleaning.
The first method we discuss here is the projection scheme, introduced by Brackbill and Barnes (1980). The idea is to measure the spurious divergence after the main update of the MHD variables, and solve for a Poisson equation defined as
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 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 points1, 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, ch 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., 1992b), features the nice property that the divergence of the magnetic field, defined in integral form as the net magnetic flux across the 6 faces of each cell, will always remain constant, So if it is initially zero, it will remain equal to zero within numerical round-off errors. This however requires to define the electric field on the cell edges, and this is the main difficulty of this approach (see Figure 1). Indeed, in the Godunov methodology, fluxes, defined at cell faces, are computed as the solution to the Riemann problem of the two neighboring cells meeting at their common cell interface. For the electric field, defined at cell edges, 4 neighboring cells meet and in order to properly upwind the electric field, we need to solve a two-dimensional Riemann problem, a rather daunting task (see Teyssier et al., 2006, for a discussion). Several solutions have been found to design 2D MHD Riemann solvers (Londrillo and del Zanna, 2004; Balsara et al., 2014), and this has been the main characteristic of several simulations codes with proper upwinding of both fluxes and electric fields, using the CT scheme within the Godunov methodology (Londrillo and del Zanna, 2004; Fromang et al., 2006; Stone et al., 2008; Li et al., 2012; Mocz et al., 2014a).
Figure 1. Schematic showing the geometry of a Cartesian cell in the Constrained Transport approach. The finite volume cell is shown as a gray cube. It is labeled 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.
2.4. Minimizing Numerical Errors: Higher Order vs. Mesh Refinement
When performing self-gravitating magnetized fluid simulations, it is of primary importance to quantify and understand numerical errors. These errors are often called truncation errors because they arise as leading order terms in the Taylor expansion of the adopted discrete numerical method. Usually, this Taylor expansion leads to the so-called modified equation, which encodes how the original ideal MHD equations have been transformed into a new set of equations with additional terms coming from the truncation errors. For first-order schemes, these terms are identical to a diffusion process, and are therefore called numerical viscosity or numerical diffusivity. In SPH, quite often, these terms are added explicitly to the equations as artificial viscosity and artificial diffusivity, while for grid-based Godunov solvers, these terms are only implicitly present, through this Taylor expansion of the discrete scheme. In both cases, however, these diffusion processes control how shock heating and magnetic reconnection proceed in the solution. They play a fundamental role in preserving the entropy condition, and in regulating the flow close to the grid scale. Unfortunately, many complex MHD processes, such as the small scale dynamo (Brandenburg and Subramanian, 2005) or the magneto-rotational instability (Balbus and Hawley, 1991) depend crucially on the so-called Prandtl number (Fromang et al., 2007; Federrath et al., 2014a), 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., 2011b), 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 (Klein, 1999; Krumholz et al., 2007b), 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)
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., 2011c; 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; Mocz et al., 2014b; Klingenberg et al., 2017; Guillet et al., 2018). In the WENO case, one needs to access neighboring 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 generalization of the CT scheme to higher order. The key ingredients are: 1- the use of a divergence-free polynomial basis for the magnetic field (Li and Shu, 2005; Guillet et al., 2018), 2- enforcing the continuity of the normal field across cells boundaries (Li and Xu, 2012). Although some very promising solutions have been found recently, this is still a very active field of research. Interestingly, the traditional CT scheme with face-centered 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-centered coefficients are required (instead of one) for each field component.
3. Non-ideal MHD: Numerical Techniques
3.1. Equations and Basic Concepts
Star formation takes place in molecular clouds, which are made of a mixture of dust and gas, implying dust and gas collisions, and both constituents are far from being fully ionized. In the previous section, we presented the work done in the ideal MHD framework, which do not appear to be well suited to the ionization 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 10 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), Marchand et al. (2016) find that amipolar diffusion and the Hall effect dominate at densities < 1012 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 ionization intensity and the dust grain size distribution (Zhao et al., 2016, 2018a; Zhao et al., 2018b; 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 behavior of the neutrals, as well as negatively and positively charged particles. Among them the different constituents that participate to the ionization 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 toward a multi-fluid treatment of the charged and neutral particles in section 3.6.
In the low ionization 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 generalized Ohm's law reads
where J = ∇ × B is the current, ηΩ, ηH, and ηAD are the Ohmic, Hall, and ambipolar resistivities, respectively (in units of cm2 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 Equation (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 generalized 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 ionization 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 ionization, 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 section 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 Equation (23) are integrated as source terms in most numerical implementations.
3.2. Ohmic Diffusion
The Ohmic diffusion is perhaps the simplest resistive term to introduce because of its diffusive nature, similar to the introduction of artificial viscosity and artificial resistivity to prevent any numerical artifact 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 10 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 NSTS sub-timesteps ΔtSTS which are based on Chebyshev polynomials to guarantee the stability conditions over the super-timestep . We also note the work of Mignone et al. (2007) and O'Sullivan and Downes (2006) where Ohmic diffusion is also integrated using the STS technique. Similarly, Tsukamoto et al. (2013) proposes a SPH discretisation of Ohmic dissipation using STS. Wurster et al. (2016) and Price et al. (2017) (in the PHANTOM code) use similar implementations on different SPH codes. Last, Marinacci et al. (2018) propose two implementations in AREPO: one using Powell divergence cleaning and another using CT. Both implementations can be integrated using an explicit or implicit scheme. Almost all the implementations mentioned so far for the Ohm diffusion are combined with ambipolar diffusion (see below).
Apart from resistive MHD, the STS technique is becoming increasingly popular in the computational astrophysics community to accelerate parabolic timestepping advancements in the context of radiation hydrodynamics with the FLD (Commerçon et al., 2011) and anisotropic diffusion (Meyer et al., 2012; Vaidya et al., 2017).
3.3. Ambipolar Diffusion
Ambipolar diffusion is not as straightforward to include since the associated term in the induction equation can not be rewritten as a Laplacian (Brandenburg and Zweibel, 1994). For the diffuse ISM, the drift between ions and neutrals is the effect of the ambipolar diffusion, that can be written as a friction force
where ρi (resp. vi) and ρn (v) are the mass density (resp. velocity) of ions and neutrals, and γAD is the collisional coupling constant, with cm3 s−1 g−1 in ISM mixture (Draine et al., 1983). In the one-fluid and strong coupling limit, the inertia of ions is neglected, as well as the pressure and gravitational forces (the ionization 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 Equation (22), we have . In the dense ISM, the heavy charged particles are ions and charged dust grains, and the expression of the resistivity is more complex (Marchand et al., 2016). In the two-fluid approximation, the ion inertia is not neglected and the system of equation accounts for both the ions and the neutrals momenta equation.
The first implementation of ambipolar diffusion in a grid code was presented in Black and Scott (1982) using a two-fluid approximation. They integrated the induction equation using the ion velocity and treated the collision between neutrals and ions as a friction force in the momenta equations. Mac Low et al. (1995) present a first 3D implementation of one-fluid ambipolar diffusion in the ZEUS code using a fully explicit scheme as well as CT. They propose to subcycle the magnetic fields evolution in the case when ambipolar diffusion timestep is much shorter than the dynamical one. The first SPH implementation can be found in Hosking and Whitworth (2004), where they used a two-fluid approximation, but their scheme did not preserve the solenoidal constraint. Duffin and Pudritz (2008) present a one-fluid implementation in the FLASH code using Powell's method (Powell et al., 1999). Masson et al. (2012) also implemented a one-fluid scheme for ambipolar diffusion in RAMSES with the same philosophy as for the Ohmic diffusion. A similar implementation using CT, operator splitting, and STS in the ATHENA code can be found in Bai and Stone (2011) and Chen and Ostriker (2014). Last, Wurster et al. (2014) propose a one-fluid implementation for SPH codes, which they implemented with success in PHANTOM.
Other one-fluid implementations in the strong coupling approximation can be found in Padoan et al. (2000) in a 3D explicit, finite difference grid code, Choi et al. (2009) in the MHD TVD code using a flux-interpolated CT scheme (Ryu et al., 1998; Balsara and Spicer, 1999), Christie et al. (2017) in the ENZO's Dedner MHD solver. Oishi and Mac Low (2006) and Li et al. (2006) present 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 section 3.6).
3.4. Hall Effect
The last resistive effect that has been implemented in star formation models is becoming increasingly popular. The Hall effect originates from the drift between ions and electrons, or, more generally the positively and negatively charged particles. The derivation of the generalized 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” . 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 cw = ω/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 (Miyama 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 of the magnetic fields is given by and at time n, explicit-implicit for the second dimension, is given by and , and finally implicit for third dimension, is given by and . 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; Marchand et al., 2018) 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.
3.5. Full Resistive Implementation
Currently, a handful of 3D codes benefit of a full implementation of the three resistive terms: PLUTO (Lesur et al., 2014), ATHENA (Bai, 2014), RAMSES (Masson et al., 2012; Marchand et al., 2018), the SPH code by Tsukamoto et al. (2017), PHANTOM (Price et al., 2017), ZeusTW (Li et al., 2011), and GIZMO (Hopkins, 2017). A very recent full implementation for solar physics can be found in González-Morales et al. (2018) in the MANCHA3D grid code. It combines STS for ambipolar and ohmic diffusion, and HDS for the Hall effect.
Currently, the ATHENA, PLUTO, RAMSES, and ZeusTW implementations rely on CT algorithms for the resistive MHD equations integration. All other works use divergence cleaning algorithms. Given the variety of the methods developed in the literature, a quantitative comparison (in standard tests as well as in star formation applications) of all these implementations would be welcome in the coming years.
Last but not least, it is not clear which non-ideal process dominates over the others. The physical conditions in star-forming regions are so wide that there is room for regions where all three effects dominate. In addition, every effect needs to be tested and quantified before any conclusions on the relative importance of each of these effects (see the review by Wurster and Li, 2018).
3.6. Multifluid Approach
Multifluid approaches are designed to describe multiphase or multifluid flows. In the context of star formation, different fluids or phases are at stack: neutrals, ions, electrons, and dust grains. A few attempts have been made to describe all or a part of these components using multi-fluid approach. In this section, we consider only the approaches which do account for N momentum equations, where N is the number of phases or fluids. For instance, a popular two-fluid framework considers ions and neutrals and describes the friction between these two fluids, i.e., the ambipolar diffusion. Toth (1994) presented a first 2D implementation of two-fluid ions and neutrals dynamics, which is based on a flux-corrected transport finite difference method. Inoue and Inutsuka (2008) propose an unconditionally stable numerical method to solve the coupling between two fluids ions-neutrals (frictional forces/heatings, ionization, and recombination, Draine, 1986) for an ISM gas subject to the thermal instability. They split time integration as follows: 1/ the ideal hydrodynamical part for neutrals, 2/ the ideal MHD part for ions. The first part is integrated using a second-order Godunov scheme. The ideal MHD part is solved in two steps, similarly to what is done for one-fluid ideal MHD solvers, i.e., Stone et al. (1992b): the magnetic pressure terms are solved using a second-order Godunov method, and the magnetic tension term and the induction equation are solved using the method of characteristics for Alfvèn waves with the CT algorithm. Tilley and Balsara (2008) present a semi-implicit method for two-fluid ion-neutral ambipolar drift in the RIEMANN code (Balsara, 2001, 2004). Their scheme is second-order accurate in space and time and uses a mixed implicit-explicit formulation to deal with strong friction regimes. Pinto et al. (2008) consider a system composed of three fluids: positively-charged particles, negatively-charged particles, and neutrals. The charged particles include ions, electrons, as well as charged dust grains, depending on the main charge carrier. The self-consistent set of equations they designed remains to be implemented in 3D codes.
Last but not least, a complete multi-fluid approach should include dust grains dynamics on top of the neutrals, ions, and electrons. Falle (2003) presented a scheme for multi-fluid hydrodynamics in the limit of small mass densities of the charged particles, i.e., where the inertia of the charged particles are neglected. This is a similar to the heavy ion approximation but for neutrals, ions, electrons, and dust grains of various sizes. The HD equations are first integrated for the neutrals using a second-order Godunov scheme, and then the induction equation and the charged species velocities are updated using an implicit scheme. This scheme captures the three non-ideal effects and can deal with regime where the Hall effect dominates. We note that similar approaches were introduced in Ciolek and Roberge (2002) and O'Sullivan and Downes (2006).
The most complete work on multi-fluid designed for star formation was presented in Kunz and Mouschovias (2009). They considered a six-fluid set of equations (neutrals, electrons, molecular and atomic ions, positively charged, negatively charged, and neutral grains) and implemented it in a 2D version of the radiation-MHD ZEUS-MP code (Hayes et al., 2006). They modified the MHD solver to account for non-ideal effects (Ohmic dissipation and ambipolar diffusion). The abundances of the 6 species are calculated using a reduced chemical-equilibrium network. They applied their methods to protostellar collapse in Kunz and Mouschovias (2010). Unfortunately, the CPU cost of such a complete set of physics is prohibitive for 3D models, and no work has taken over since this first application.
3.7. What Is Next? Dust and Gas Mixture Dynamics
Currently, numerous studies report on the effect of non-ideal MHD in the star formation process. The wide majority of these works use the single-fluid approximation, assuming furthermore that the dust grains are perfectly coupled to the gas via collisions. Nevertheless, dust grains are observed to cover a wide range of sizes, from nanometer to micrometer in the ISM (Mathis et al., 1977) and even millimeter in protostellar disks environments (Chiang et al., 2012). Generally, it assumes that dust grains of different sizes coexist and follow a power-law dust distribution. Grains react to changes in the gas velocity via a drag force. For a given dust grain of mass mg, the force is
where vg is the dust grain velocity, v the gas velocity, and ts 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 ts ~ s (Epstein, 1924). The Stokes number St ≡ ts/tdyn characterizes the coupling of the dust grain relative to the gas, where well coupled dust grains follow St < 1. tdyn 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; Laibe and Price, 2012; Meheut et al., 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) (Laibe and Price, 2012), where prohibitive spatial and time discretisations may be required. 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 grid-based codes in the past 5 years (Laibe and Price, 2014b; Lebreuilly et al., 2019). Last, Price and Laibe (2015); Lin and Youdin (2017) present methods for simulating the dynamics of small grains with only one additional equation on the dust concentration on top of the Euler equations. These last methods have been recently implemented with success in PHANTOM (Price et al., 2017), as well as in PLUTO (Chen and Lin, 2018) and in RAMSES (Lebreuilly et al., 2019).
Recent works have emphasized the possible decoupling of micrometer dust grains in molecular clouds (Hopkins and Lee, 2016; Tricco et al., 2017) as well as of millimeter grains in collapsing dense cores (Bate and Lorén-Aguilar, 2017, see Figure 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.
Figure 2. SPH calculations of dust dynamics during the collapse of a 1 M⊙ dense core. The left column shows the gas column density (edge-on view) and from left to right, the three other columns show the dust column density for dust grain sizes of 10 μm, 100 μm, and 1 mm. Large dust grains decouple from the gas and settle in the mid-plane faster into a large dusty disc with larger dust-to-gas ratio than in the envelop. Figure adapted from Bate and Lorén-Aguilar (2017) with permission from the authors.
4. Radiative Transfer: Numerical Techniques
4.1. A Poor's Man Approach: The Polytropic Equation of State
Radiative transfer plays a fundamental role in star formation and cannot be ignored, even for a review on magnetic fields. The importance of radiation fields is 2-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., 2007a; Myers et al., 2013). In conclusion, radiation plays a central role during the formation and the evolution of the first Larson core (see section 5), but later on acts as a self-regulating mechanism for star formation.
4.2. Ray Tracing and Long Characteristic Methods
The radiative transfer equation written in full generality using the radiation specific intensity Iν(x, n, t)
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).
4.3. Monte Carlo and Short Characteristic Methods
Although long characteristics offer the advantage of solving the exact attenuation along light rays, an especially accurate approach for point sources, they are tricky to parallelise because light rays propagate over multiple MPI domains. An alternative technique is the short characteristics method, for which light rays as short as the hydrodynamics cell size are considered. The radiation intensity from neighboring light rays is interpolated between neighboring cells, avoiding the use of long characteristics. Rijkhorst et al. (2006) use an hybrid long- and short-characteristics method in the FLASH code to overcome the difficulty to parallelise the long characteristics method. The short characteristic method also offers better scaling properties. The short characteristics method has also been used in the ATHENA code to solve the stationary radiative transfer equation by Davis et al. (2012) and in the C2-ray code for ionization problems (Mellema et al., 2006).
A third method which shares with the two previous ones the property of being accurate and able to capture complex radiation geometry is the Monte Carlo method. The Monte Carlo approach samples the radiation field with photon packets carrying information such as angle of propagation and frequency. This is as close as one can be from true photons propagating in the fluid and allows for complex scattering processes to be included in the models (Ercolano et al., 2003; Altay et al., 2008; Harries, 2011; Roth and Kasen, 2015). The method can be however quite expensive, as large particle numbers are required to avoid that the Poisson noise dominates the computed signal. In the diffusion limit, the mean free path becomes much smaller than the length scale of the system and the number of interactions becomes prohibitively large: It scales as the square of the optical depth (Nayakshin et al., 2009). It has also been mostly developed for stationary problems but many codes able to deal with non-stationary problems are now emerging (Nayakshin et al., 2009; Harries, 2011; Lomax and Whitworth, 2016).
4.4. Two Moments Methods
The common aspect in these three methods is that they can be quite expensive. This is probably not true for the short characteristics method, although the number of angular domains adopted plays a role in the final accuracy of the solution, and affects the cost. An alternative approach, called moment-based radiative transfer, can potentially solve this cost problem by integrating the radiative transfer equation over the angles. We obtain the radiation energy and the radiation momentum equation as
where Eν is the radiation energy, Fν is the radiation flux and ℙ = 𝔻Eν is the radiation pressure tensor and 𝔻 is the Eddington tensor. For slowly varying spectral absorption and emission coefficients, for which one has
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., 2007b, 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 and are expressed in the comoving frame, these equations are often called mixed frame radiation moments equations. We can also Lorentz-transform the radiation moments from the laboratory frame to the comoving frame, which gives, to first order in (v/c)
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 a detailed implementation).
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/r2 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
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 𝔻 ≃ 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.
4.5. Flux-Limited Diffusion Methods
A definitive advantage of moment-based radiative transfer methods is the possibility to model both optically thin and optically thick regimes. For optically thick conditions, the radiation field follows the diffusion limit of the radiative transfer equation. In that limit, the radiation field in the comoving frame is quasi-isotropic, so that and the radiation flux can be approximated by
Since now the flux in the comoving frame is a direct function of the energy density in the comoving frame, we do not need the radiation flux equation anymore and we are left with only the radiation energy equation. The Lorentz transform can be written to leading order as
Injecting these into the radiation energy equation leads to a fully covariant formulation in the diffusion limit
which can be written in the familiar form
This is the radiation energy equation in the diffusion limit using only comoving radiation variables. If one wants to extend this equation outside its validity range (high optical depth), then one can use Flux Limited Diffusion (FLD), for which the flux function is modified as
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., 1992a). For AMR, the ORION code was developed specifically for FLD in collapsing star forming core by Krumholz et al. (2007b), 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 Commerçon 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., 2010b; 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.
5. 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.
5.1. Historical Work and 1D Studies
Contemporary numerical star formation studies started in the late 60s with the outstanding pioneering work of Larson (1969) who first computed numerically the collapse of a dense core down to the formation of the protostar. He used a 1D spherically symmetric model, accounting for coupled gas dynamics and radiative transfer in a modified Eulerian scheme. Larson identified two distinct stages during the protostellar collapse, called the first and the second collapse. Each of these two stages are followed by the formation of a hydrostatic object, commonly referred to as the first and second Larson cores in post (Larson, 1969) studies. Larson's work established an empirical evolutionary sequence as follows. Dense cores first collapse isothermally because dust thermal emission is very efficient at radiating away the gas compression energy (dust and gas are thermally coupled within dense cores, see for instance Galli et al., 2002). The first collapse stops at densities of the order of ≃ 10−13 g cm−3 at which dust grains become opaque to their own radiation so that the gas begins to heat up quasi-adiabatically. This is the formation of the first Larson core. The first collapse lasts roughly a free-fall time. The first core accretes matter and its temperature increases up about 2000K where H2 dissociation starts. This endothermic reaction enables the gas to behave as a fluid with an effective polytropic index smaller than the critical value 4/3 for collapse. The second collapse thus starts at typical densities of ≃ 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 H2 dissociation. Masunaga and Inutsuka (2000) integrated their models throughout the Class 0 and Class I phases up to an age of 1.3 × 105 yr. Nowadays, 1D models are still used either as a first step toward describing more accurately the physics of the protostellar collapse, e.g., with multigroup radiative transfer (Vaytet et al., 2013, 2014) or with a view to provide quantitative predictions for the first and second Larson core properties (Vaytet and Haugbølle, 2017; Bhandare et al., 2018). These models are not pushed in time as far as the ones by Masunaga and Inutsuka (2000).
5.2. From 1D to 3D
Pioneered studies acknowledged up front that 1D spherical models were limited because the effect of rotation, turbulence, and magnetic fields cannot be taken into account. Three decades after Larson work, the first 3D numerical simulation able to describe the formation of the protostar was performed by Bate (1998) using SPH. The model accounted for initial rotation, but magnetic fields were neglected and radiative transfer was crudely mimicked by a piecewise polytropic equation of state. It is worth to notice at this point that grid-based codes took some time to reproduce (Bate, 1998) results for almost one decade. SPH is indeed naturally well-suited to collapse problems thanks to its Lagrangian nature while innovative techniques, such as AMR and nested grids, were not mature enough to capture the entire dynamical range and to cover eight orders of magnitude in physical length given the computer capabilities at that time. As mentioned in the section 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; Winkler and Norman, 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 behavior 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 section 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 H2 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, 2018b). The driver of the second collapse is the endothermic dissociation of the H2 molecule, which modifies considerably the gas equation of state. Currently, there are two ways to integrate non-ideal EOS in second collapse calculations. The first one is to compute on-the-fly the thermodynamic quantities from equilibrium chemical abundances and to use an ideal gas equation of state, including the effects of dissociation and ionization, 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., ionization 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 work by Vaytet et al. (2013), Vaytet et al. (2014), Vaytet and Haugbølle (2017), and Vaytet et al. (2018), and also in the SPH public code PHANTOM (Price et al., 2017). We note that these developments imply that the dependency of the specific heat capacity as a function of the physical conditions has to be taken into account. Caution should thus be taken within the radiation solvers which are often integrated using implicit schemes. The most common and simple way to deal with it is to compute a heat capacity-like factor, 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 × 106 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., 2013, 2018). The infall velocity is indeed greatly modified in the presence of magnetic fields which tends to focus the collapsing gas on a smaller area compared to the case without magnetic fields (Commerçon et al., 2011b; Myers et al., 2013). The magnetic braking mechanism transports angular momentum from the inner parts of the collapsing cloud to the envelop. As a consequence, the infall velocity increases along with the incident kinetic energy, which then modifies the thermal budget through the accretion shocks. The next step in the increasing complexity is to perform full 3D RMHD calculations. The first SPH RMHD calculations of protostellar core formation were performed in Bate et al. (2014), extending their previous work to ideal MHD. Tomida et al. (2013) and 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 section 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; Wurster et al., 2018) and AMR (Vaytet et al., 2018). The typical resolution in SPH calculations currently achievable with limited CPU resources is of 3 × 106 equal-mass particles per M⊙. In AMR, Vaytet et al. (2018) used a coarse grid of 643 with 21 additional levels of refinement, and a refinement criterion of 32 points per Jeans length. While the SPH calculations are able to integrate the three non-ideal effects (Wurster et al., 2018) down to stellar scales, there is no AMR work describing the Hall effect in second collapse calculations (Vaytet et al., 2018 accounts for the ambipolar diffusion and the Ohmic diffusion). The calculations by Wurster et al. (2018) were integrated 17 yr after the stellar core formation, while Vaytet et al. (2018) could only integrate for 24 days. Beside the major achievement which represents the formation of a protostar on a computer, we must admit that these kind of studies are not yet capable to give long term predictability. In addition, the CPU time required to integrate full 3D RMHD models is relatively big. Clearly, second collapse study cannot be done on a frame covering a large parameter space. In addition, even at the first core scale, the integration timestep becomes so short in the adiabatic fragments that long time evolution calculations becomes prohibitive. In that view, the sink particles are a good alternative to reach a long horizon of predictability and to cover a wide parameter space.
6. 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., 2010a; Li et al., 2011). The differences between sink particles and sink cells are discussed in section 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.
6.1. Standard Implementations
Sink particles were introduced in star formation calculations by Bate et al. (1995) to enable long time integration to study the fragmentation of collapsing clouds. Sink particles were introduced to mimic the second collapse and the protostars at scales that could not be reached in numerical simulations. They are treated as accreting non-collisional point masses and enable one to follow the evolution through the main accretion phase, with a view to compare with observations of clusters in which most of the gas has been accreted. Their original implementation was done in a SPH code, with three main parts in the algorithm: (i) sink particle creation, (ii) sink particle accretion, (iii) boundary conditions for sink particles. A SPH gaseous particle is tagged for sink particle creation if its density exceeds a specific density threshold (defined by the user) and satisfied a number of tests. First, the smoothing length of the SPH particle that is a candidate for sink particle creation has to be less than half the accretion radius racc of the sink particle. This ensures that the sink particle is formed from at least Nneigh. racc 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 racc 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 neighbors 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 (). If all the tests are passed, a sink particle is formed from the system at the center of mass. Initially, the sink particle has thus a mass of ~Nneigh 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 racc 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 neighbors. 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, cs the gas sound speed, and Δxmin the minimum cell size. The initial mass of the sink is . 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 racc. 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 racc. 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 racc = 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 racc 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 83 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., where ϵ = 2Δxmin 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 (vsinkΔt < CΔx, with C < 1).
Sink particles algorithms have been complemented in the past 10 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.
6.2. Contemporary Implementations
In all current implementations, the sink creation and accretion are the most sensitive algorithms. Of course, the sink particles can be introduced in the initial conditions. Otherwise, they are created after the creation checks.
Each of the commonly used codes in star formation has its own implementation of sink particles. For the grid-based codes, we note the work of Wang et al. (2010) in the ENZO code, which is very similar to Krumholz et al. (2004) but with a different merging scheme which involves a maximum mass sink mass for merging. Federrath et al. (2010) presented an implementation in the AMR code FLASH, where the creation checks were similar to the ones by Bate et al. (1995) for SPH with two additional checks: the cell tagged for sink creation should be at the maximum level of refinement, and it should sit in a local gravitational potential minimum. They also account for magnetic energy in the energy budget done in the check algorithm. The checks are performed on a spherical region of radius racc 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 racc exceeds the density threshold ρres for sink creation, the mass increment 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. (2011a) 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 neighbors), (iv) the gaseous particle is outside the Hill sphere of fragments harboring another already formed sink.
Last, Bleuler and Teyssier (2014) presented an innovative sink creation scheme based on a clump finder algorithm used on-the-fly in the RAMSES code. The clump finder is performed in a first step to identify peaks and their associated regions. Second, the peaks are considered for sink creation. A virial theorem type analysis is performed on the clumps, accounting for surface pressure as well as tidal forces. The virial check is passed if the second derivative of the moment of inertia in the center of mass frame is negative. Then a collapse check is performed, which corresponds to the contraction check of Krumholz et al. (2004) and Federrath et al. (2010), but adapted according to the virial analysis. Last, the final check is a proximity check. We also note that Bleuler and Teyssier (2014) offer also the possibility to use the alternative sink creation and accretion schemes reported in the literature. They provided an excellent review of the different creation checks, accretion, sink merging and trajectories schemes. In particular, they tested different schemes to perform sink accretion: Bondi-Hoyle accretion, flux accretion, and density threshold accretion. They show that all methods give good results for spherical Bondi accretion (i.e., when the flow is supersonic). In the subsonic regime, only the Bondi formula gives a correct accretion rate. They recommend to automatically adapt the accretion scheme depending on the surrounding flow properties.
A number of other sink particle implementations have been reported in the literature that we cannot detail in this review. We give in the following a non-exhaustive list, restricted to the field of star formation. In SPH, sink particles have been implemented in the codes GADGET (Jappsen et al., 2005), DRAGON (Goodwin et al., 2004), and SEREN (Hubber et al., 2011). For grid-based codes, sink particle are used in the uniform grid code ATHENA by Gong and Ostriker (2013), in RAMSES by Padoan and Nordlund (2011), and in the ORION2 code by Lee et al. (2014). Last, Greif et al. (2011) also proposed a sink particle algorithm within the AREPO moving-mesh code. Their implementation is very similar to the one used in SPH codes. We also note that sink particle merging is now widely used in most of the implementations.
Twenty years after their introduction, sink particles remain a sensible 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 behaviors 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.
6.3. Sink Particles and Magnetic Fields
In the previous section, we did not discuss the sink particle algorithms in the presence of magnetic fields. We try to describe here some attempts done to provide a good description of MHD flows for sink creation and accretion. Sink particles can of course be introduced in MHD calculations, but the implementation is not as straightforward as in the hydrodynamical case. While it is relatively easy to accrete gas and momentum, magnetic flux accretion cannot bypass the solenoidal constraint .
Price and Bate (2008) and Wang et al. (2010) performed star formation calculations in turbulent and magnetized clouds, but do not account for any contribution of the magnetic fields in their sink algorithms. Federrath et al. (2010) have introduced the magnetic energy in their negative total energy check for sink creation, but not in the accretion scheme. Lee et al. (2014) have extended the work of Krumholz et al. (2004) to deal with magnetized flows. For sink creation, they follow the work of Myers et al. (2013) who derived a magnetic Truelove criterion. They derive simple steady state accretion rates as a function of the plasma beta β and the Mach number .
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; Masson 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 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).
6.4. Sink Particles vs. 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 outflow (Yorke et al., 2002; Kuiper et al., 2010a; 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, Mdisk/Mstar > 1/3 (Krumholz et al., 2007a; Kuiper et al., 2011; Sigalotti et al., 2018).
Last, the use of “outflow” boundary conditions for the mass accretion on the sink cell is very different from the algorithms we just presented for the sink particles. Machida et al. (2014) have shown that applying similar criteria for the sink accretion than the ones employed in studies with a sink cell (Li et al., 2011) may lead to very different results, i.e., no disk formation, compared to the same calculations done with a classical sink particle algorithm, i.e., formation of a disk.
7. 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, ionization), and dynamical (jets/outflow, winds). In any case, the source terms introduced by sub-grid physics are treated using operator splitting schemes.
7.1. Radiative Feedback
The first class of sub-grid models were designed to account for the radiative feedback of the nascent protostars. The radiative feedback term includes the accretion, the internal, and the ionization 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 gray 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 Lacc = 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 Lint. They realized 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 gray model. In addition, Krumholz et al. (2005) showed how protostellar outflows allow radiation to escape in the outflow cavities to allow a continued accretion onto the central massive stars. These results have been further confirmed by 3D dynamical simulations of Krumholz et al. (2007a), Kuiper et al. (2011, 2016), and Rosen 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 gray FLD model for the radiation transport where the protostellar luminosity is simply a source term in the radiative energy equation (Krumholz et al., 2007a; Offner et al., 2009; Stamatellos et al., 2011; Fontani et al., 2018; Jones and Bate, 2018). Frequency dependent protostellar irradiation modules have been developed using ray-tracing to compute the radiation fields, combined with moment models for radiation hydrodynamics (Kuiper et al., 2010a; 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 (Klassen et al., 2016; Rosen et al., 2016).
Last, we note the developments made to capture ionizing 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 (Dale and Bonnell, 2011; Peters et al., 2011; Dale et al., 2012; Geen et al., 2015, 2018; Ali et al., 2018).
7.2. Dynamical Feedback
The most advanced MHD collapse calculations have shown that outflows and jets naturally develop at different scales during the collapse (Tomisaka, 2002; Banerjee and Pudritz, 2006; Machida et al., 2006; Hennebelle and Fromang, 2008; Ciardi and Hennebelle, 2010; Commerçon et al., 2010; Tomida et al., 2010; Price et al., 2012). In addition, outflows and jets are commonly observed in YSO's (it is for instance a selection criterion for Class 0 sources, Andre et al., 1993) and they are considered as possible sources of turbulence driving in star forming clouds and participate in the regulation of the star formation rate (Matzner, 2007; Nakamura and Li, 2007; Krumholz, 2014; Federrath, 2015). Outflows and jets launching scales are nevertheless not captured in most studies because it requires a very high numerical resolution to reach sub-AU scales. Besides, studies neglecting magnetic fields cannot generate MHD winds driven centrifugally or by toroidal magnetic pressure. Consequently, very few numerical calculations have been able to launch self-consistently MHD outflows (e.g., Hennebelle et al., 2011). Sub-grid models have thus been developed to account for outflows/jets in numerical works with unresolved launching scales and missing physics, similarly to what has just been presented for the radiative feedback. The numerical implementation are based on analytical works dedicated to MHD wind launching (Blandford and Payne, 1982; Pelletier and Pudritz, 1992; Shu et al., 1994; Ferreira, 1997; Matzner and McKee, 1999).
The first works on sub-grid models have been reported in Li and Nakamura (2006) and Nakamura and Li (2007) who employed ideal MHD and self-gravity. They attached a protostellar outflow model on sink particles following (Matzner and McKee, 2000) because they were unable to describe the self-consistent launch of outflows with the crude resolution they use (1283 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 ). 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. (2005) and Krumholz et al. (2007a) and presented for the first time a combined sub-grid model which includes protostellar outflows (Matzner and McKee, 1999) and protostellar radiative heating (Offner et al., 2009), ignoring magnetic fields though. Their sub-grid model accounts for pre-main sequence evolution as in Offner et al. (2009) 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 104 K, which is appropriate for an ionized 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. (2014b) 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 normalized 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 (Offner et al., 2009; Cunningham et al., 2011; Federrath et al., 2014b; 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. Federrath et al. (2014b) 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.
All these protostellar outflow sub-grid models are currently widely used in the community for astrophysical applications on the star formation rate/efficiency (e.g., Federrath, 2015; Kuiper et al., 2016; Offner and Chaban, 2017; Murray et al., 2018), stellar initial mass function (Krumholz et al., 2012), origin of turbulence driving (Nakamura and Li, 2011; Offner and Liu, 2018), as well as cluster and massive star formation (Wang et al., 2010; Cunningham et al., 2011; Kuiper et al., 2015; Li et al., 2018).
In the context of massive star formation, another class of outflows models are based on the luminosity and ionization. Sub-grids models of outflows generated by massive star luminosity have been developed (e.g., Krumholz et al., 2005; 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.
7.3. 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; Beuther et al., 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 (Capuzzo-Dolcetta 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 non-ideal 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 ionization rises due to high temperature and the system is back to ideal MHD, but with a high plasma beta β = Ptherm/Pmag. 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; Vaytet et al., 2018; Tomida et al., 2013) 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 1 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, Wurster et al. (2018) 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 .
Figure 3. Volume rendering from collapsing 1 M⊙ dense cores at the protostellar core scale. The left panel shows the results of an ideal MHD 3D simulations, and the right panel the comparison with non-ideal MHD including ambipolar diffusion and Ohmic diffusion. The color lines indicate the 3D magnetic field lines, and the color coding the magnetic field amplitude. Figure reproduced from Vaytet et al. (2018) with permission of A&A.
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.
Modeling 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 modeled individually, at least for massive galaxies like the Milky Way (see Figure 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 × 103 M⊙ with RAMSES. Renaud et al. (2013) and 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) analogs 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, Hopkins et al. (2018) 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 modeled individually: Hu et al. (2016); Hu (2019) using SPH and a resolution of 1 − 5 M⊙ by gas particle, Emerick et al. (2018) and Emerick et al. (2019) using ENZO with a maximum resolution of 1.8 pc.
Figure 4. Simulation of the Milky Way: column density of the gas disc in a sub-parsec resolution simulation. The color table only applies to the main panel: the table has been changed in each zoom-in view to enhance contrast. Figure adapted from Renaud et al. (2013) with permission from the authors.
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), and Colling et al. (2018) as well as the work of Martizzi et al. (2016) and Butler et al. (2015, 2017) using RAMSES, by Kim et al. (2011), Kim et al. (2013), Kim and Ostriker (2015), Kim and Ostriker (2017), and Kim and Ostriker (2018) using ATHENA, the SILCC project papers using FLASH (Walch et al., 2015; Girichidis et al., 2016, 2018; Gatto et al., 2017; Peters et al., 2017).
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 103 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 modeled. 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, 2012) and magnetic fields (Price and Bate, 2008; Price and Bate, 2009). The work of Dale et al. (2007), Dale and Bonnell (2011), and 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. (2011), Girichidis et al. (2012b), and Girichidis et al. (2012a) using FLASH who studied the importance of isolated initial conditions in isothermal cloud core collapse (without stellar feedback). Federrath et al. (2014b) also used FLASH to perform simulations of isolated cores with protostellar feedback (jets). Krumholz et al. (2007a), Krumholz et al. (2012), Hansen et al. (2012), Myers et al. (2013), Myers et al. (2014), Cunningham et al. (2011), Li et al. (2018), and 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: Hennebelle et al. (2011) and Commerçon et al. (2011b) for 100 M⊙ isolated core collapse models with (R)MHD, Lee and Hennebelle (2018a), Lee and Hennebelle (2018b), and Lee and Hennebelle (2018c) focused on the peak of the IMF using (M)HD 1000 M⊙ isolated core collapse models, Geen et al. (2015), Geen et al. (2016), and Gavagnin et al. (2017) 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) and 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 (Zamora-Avilés et al., 2018). They used FLASH with a maximum resolution of 25 AU and were able to cover scales from 256 pc to 25 AU (see Figure 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.
Figure 5. Simulation of protostellar disk formation within molecular clouds. The column density maps show the entire molecular (upper left panel) and successive zoom within a star-forming dense core. Figure adapted from Zhang et al. (2018) with permission from the authors.
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 40963 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 vs. box size, statistics vs. internal structure, physical realism vs. 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 (Krumholz and McKee, 2005; Hennebelle and Chabrier, 2011; Padoan and Nordlund, 2011; Federrath and Klessen, 2012) but the methodology, in the context of the full spectrum of required physical processes, is still at its infancy.
RT has written the introduction, the ideal MHD part and the radiative transfer part, as well as the discussion. BC has written the non-ideal MHD part, the second collapse part, the sink particle part and the subgrid models part.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to warmly thank our two referees, Richard Klein and Christoph Federrath, for their help in reviewing the manuscript. We thank Matthew Bate, Florent Renaud, and Shangjia Zhang for providing adapted figures from their work.
1. ^Stagnation points are regions where the flow is brought to rest in the frame of the system under study. More generally, even if the velocity does not vanish, stagnation points are regions in the flow where the streamlines are converging and the flow becomes compressive.
Ali, A., Harries, T. J., and Douglas, T. A. (2018). Modelling massive star feedback with Monte Carlo radiation hydrodynamics: photoionization and radiation pressure in a turbulent cloud. Month. Notices R. Astron. Soc. 477, 5422–5436. doi: 10.1093/mnras/sty1001
Altay, G., Croft, R. A. C., and Pelupessy, I. (2008). SPHRAY: a smoothed particle hydrodynamics ray tracer for radiative transfer. Month. Notices R. Astron. Soc. 386, 1931–1946. doi: 10.1111/j.1365-2966.2008.13212.x
Andre, P., Ward-Thompson, D., and Barsony, M. (1993). Submillimeter continuum observations of rho ophiuchi A: the candidate protostar VLA 1623 and prestellar clumps. Astrophys. J. 406:122. doi: 10.1086/172425
Aubert, D., and Teyssier, R. (2008). A radiative transfer scheme for cosmological reionization based on a local Eddington tensor. Month. Notices R. Astron. Soc. 387, 295–307. doi: 10.1111/j.1365-2966.2008.13223.x
Aubert, D., and Teyssier, R. (2010). Reionization simulations powered by graphics processing units. I. On the structure of the ultraviolet radiation field. Astrophys. J. 724, 244–266. doi: 10.1088/0004-637X/724/1/244
Bai, X.-N., and Stone, J. M. (2011). Effect of ambipolar diffusion on the nonlinear evolution of magnetorotational instability in weakly ionized disks. Astrophys. J. 736:144. doi: 10.1088/0004-637X/736/2/144
Balsara, D. S., Dumbser, M., and Abgrall, R. (2014). Multidimensional HLLC Riemann solver for unstructured meshes - With application to Euler and MHD flows. J. Comput. Phys. 261, 172–208. doi: 10.1016/j.jcp.2013.12.029
Balsara, D. S., and Spicer, D. S. (1999). A staggered mesh algorithm using high order godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations. J. Comput. Phys. 149, 270–292. doi: 10.1006/jcph.1998.6153
Bate, M. R. (2010). Collapse of a molecular cloud core to stellar densities: the radiative impact of stellar core formation on the circumstellar disc. Month. Notices R. Astron. Soc. 404, 79–83. doi: 10.1111/j.1745-3933.2010.00839.x
Bate, M. R. (2012). Stellar, brown dwarf and multiple star properties from a radiation hydrodynamical simulation of star cluster formation. Month. Notices R. Astron. Soc. 419, 3115–3146. doi: 10.1111/j.1365-2966.2011.19955.x
Bate, M. R., Bonnell, I. A., and Bromm, V. (2003). The formation of a star cluster: predicting the properties of stars and brown dwarfs. Month. Notices R. Astron. Soc. 339, 577–599. doi: 10.1046/j.1365-8711.2003.06210.x
Bate, M. R., and Burkert, A. (1997). Resolution requirements for smoothed particle hydrodynamics calculations with self-gravity. Month. Notices R. Astron. Soc. 288, 1060–1072. doi: 10.1093/mnras/288.4.1060
Bate, M. R., Tricco, T. S., and Price, D. J. (2014). Collapse of a molecular cloud core to stellar densities: stellar-core and outflow formation in radiation magnetohydrodynamic simulations. Month. Notices R. Astron. Soc. 437, 77–95. doi: 10.1093/mnras/stt1865
Beuther, H., Linz, H., and Henning, T. (2008). “Massive star formation: observations confront theory,” in ASP Conference Series, Vol. 387 (San Francisco, CA: Astronomical Society of the Pacific), 189.
Bhandare, A., Kuiper, R., Henning, T., Fendt, C., Marleau, G.-D., and Kölligan, A. (2018). First core properties: from low- to high-mass star formation. arXiv:1807.06597. doi: 10.1051/0004-6361/201832635
Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., and Borgani, S. (2011). A non-ideal magnetohydrodynamic GADGET: simulating massive galaxy clusters. Month. Notices R. Astron. Soc. 418, 2234–2250. doi: 10.1111/j.1365-2966.2011.19523.x
Boscheri, W., and Dumbser, M. (2014). A direct Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume scheme on unstructured tetrahedral meshes for conservative and non-conservative hyperbolic systems in 3D. J. Comput. Phys. 275, 484–523. doi: 10.1016/j.jcp.2014.06.059
Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., and Puerari, I. (2010). ISM properties in hydrodynamic galaxy simulations: turbulence cascades, cloud formation, role of gravity and feedback. Month. Notices R. Astron. Soc. 409, 1088–1099. doi: 10.1111/j.1365-2966.2010.17370.x
Brackbill, J. U., and Barnes, D. C. (1980). The effect of nonzero product of magnetic gradient and B on the numerical solution of the magnetohydrodynamic equations. J. Comput. Phys. 35, 426–430. doi: 10.1016/0021-9991(80)90079-0
Bryan, G. L., Norman, M. L., O'Shea, B. W., Abel, T., Wise, J. H., Turk, M. J., et al. (2014). ENZO: an adaptive mesh refinement code for astrophysics. Astrophys. J. Suppl. 211:19. doi: 10.1088/0067-0049/211/2/19
Buntemeyer, L., Banerjee, R., Peters, T., Klassen, M., and Pudritz, R. E. (2016). Radiation hydrodynamics using characteristics on adaptive decomposed domains for massively parallel star formation simulations. New Astron. 43, 49–69. doi: 10.1016/j.newast.2015.07.002
Butler, M. J., Tan, J. C., Teyssier, R., Rosdahl, J., Van Loo, S., and Nickerson, S. (2017). Kiloparsec-scale simulations of star formation in disk galaxies. IV. regulation of galactic star formation rates by stellar feedback. Astrophys. J. 841:82. doi: 10.3847/1538-4357/aa7054
Butler, M. J., Tan, J. C., and Van Loo, S. (2015). Kiloparsec-scale simulations of star formation in disk galaxies III. Structure and dynamics of filaments and clumps in giant molecular clouds. Astrophys. J. 805:1. doi: 10.1088/0004-637X/805/1/1
Capuzzo-Dolcetta, R., Limongi, M., and Tornambè, A. (2012). “Advances in computational astrophysics: methods, tools, and outcome,” in ASP Conference Proceedings, Vol. 453 (San Francisco, CA: Astronomical Society of the Pacific), 129.
Chiang, H.-F., Looney, L. W., and Tobin, J. J. (2012). The envelope and embedded disk around the class 0 protostar L1157-mm: dual-wavelength interferometric observations and modeling. Astrophys. J. 756:168. doi: 10.1088/0004-637X/756/2/168
Ciardi, A., and Hennebelle, P. (2010). Outflows and mass accretion in collapsing dense cores with misaligned rotation axis and magnetic field. Month. Notices R. Astron. Soc. 409, L39–L43. doi: 10.1111/j.1745-3933.2010.00942.x
Ciolek, G. E. L. E. C., and Roberge, W. A. G. R. (2002). Time-dependent, multifluid, magnetohydrodynamic shock waves with grain dynamics I. Formulation and numerical tests. Astrophys. J. 567, 947–961. doi: 10.1086/338591
Clark, P. C., Glover, S. C. O., and Klessen, R. S. (2012). TreeCol: a novel approach to estimating column densities in astrophysical simulations. Month. Notices R. Astron. Soc. 420, 745–756. doi: 10.1111/j.1365-2966.2011.20087.x
Colling, C., Hennebelle, P., Geen, S., Iffrig, O., and Bournaud, F. (2018). Impact of galactic shear and stellar feedback on star formation. Astron. Astrophys. 620:A21. doi: 10.1051/0004-6361/201833161
Commercon, B., Debout, V., and Teyssier, R. (2014). A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids. Astron. Astrophys. 563:A11. doi: 10.1051/0004-6361/201322858
Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., and Teyssier, R. (2008). Protostellar collapse: a comparison between smoothed particle hydrodynamics and adaptative mesh refinement calculations. Astron. Astrophys. 482:371. doi: 10.1051/0004-6361:20078591
Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., and Teyssier, R. (2010). Protostellar collapse: radiative and magnetic feedbacks on small-scale fragmentation. Astron. Astrophys. 510:L3. doi: 10.1051/0004-6361/200913597
Commerçon, B., Hennebelle, P., and Henning, T. (2011b). Collapse of massive magnetized dense cores using radiation magnetohydrodynamics: early fragmentation inhibition. Astrophys. J. Lett. 742:L9. doi: 10.1088/2041-8205/742/1/L9
Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., and Chabrier, G. (2011). Radiation hydrodynamics with adaptive mesh refinement and application to prestellar core collapse. I. Methods. Astron. Astrophys. 529:A35. doi: 10.1051/0004-6361/201015880
Cunningham, A. J., Klein, R. I., Krumholz, M. R., and McKee, C. F. (2011). Radiation-hydrodynamic simulations of massive star formation with protostellar outflows. Astrophys. J. 740:107. doi: 10.1088/0004-637X/740/2/107
Dale, J. E., and Bonnell, I. (2011). Ionizing feedback from massive stars in massive clusters: fake bubbles and untriggered star formation. Month. Notices R. Astron. Soc. 414, 321–328. doi: 10.1111/j.1365-2966.2011.18392.x
Dale, J. E., Ercolano, B., and Bonnell, I. A. (2012). Ionizing feedback from massive stars in massive clusters - II. Disruption of bound clusters by photoionization. Month. Notices R. Astron. Soc. 424, 377–392. doi: 10.1111/j.1365-2966.2012.21205.x
Dale, J. E., Ercolano, B., and Clarke, C. J. (2007). A new algorithm for modelling photoionizing radiation in smoothed particle hydrodynamics. Month. Notices R. Astron. Soc. 382, 1759–1767. doi: 10.1111/j.1365-2966.2007.12486.x
de Avillez, M. A., and Breitschwerdt, D. (2005). Global dynamical evolution of the ISM in star forming galaxies. I. High resolution 3D simulations: effect of the magnetic field. Astron. Astrophys. 436, 585–600. doi: 10.1051/0004-6361:20042146
Dedner, A., Kemm, F., Kröner, D., Munz, C. D., Schnitzer, T., and Wesenberg, M. (2002). Hyperbolic divergence cleaning for the MHD equations. J. Comput. Phys. 175, 645–673. doi: 10.1006/jcph.2001.6961
Derigs, D., Winters, A. R., Gassner, G. J., and Walch, S. (2016). A novel high-order, entropy stable, 3D AMR MHD solver with guaranteed positive pressure. J. Comput. Phys. 317, 223–256. doi: 10.1016/j.jcp.2016.04.048
Derigs, D., Winters, A. R., Gassner, G. J., Walch, S., and Bohm, M. (2017). Ideal GLM-MHD: about the entropy consistent nine-wave magnetic field divergence diminishing ideal magnetohydrodynamics equations. arXiv, 420–467. doi: 10.1016/j.jcp.2018.03.002
Dobbs, C. L., Pettitt, A. R., Corbelli, E., and Pringle, J. E. (2018). Simulations of the flocculent spiral M33: what drives the spiral structure? Month. Notices R. Astron. Soc. 478, 3793–3808. doi: 10.1093/mnras/sty1231
Dubroca, B., and Feugeas, J. (1999). Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif. Academie des Sciences Paris Comptes Rendus Serie Sciences Mathematiques 329, 915–920. doi: 10.1016/S0764-4442(00)87499-6
Duffin, D. F., and Pudritz, R. E. (2008). Simulating hydromagnetic processes in star formation: introducing ambipolar diffusion into an adaptive mesh refinement code. Month. Notices R. Astron. Soc. 391, 1659–1673. doi: 10.1111/j.1365-2966.2008.14026.x
Dzyurkevich, N., Commerçon, B., Lesaffre, P., and Semenov, D. (2017). Magnetic diffusivities in 3D radiative chemo-hydrodynamic simulations of protostellar collapse. Astron. Astrophys. 603:A105. doi: 10.1051/0004-6361/201628995
Emerick, A., Bryan, G. L., and Mac Low, M.-M. (2018). Stellar radiation is critical for regulating star formation and driving outflows in low-mass dwarf galaxies. Astrophys. J. 865:L22. doi: 10.3847/2041-8213/aae315
Emerick, A., Bryan, G. L., and Mac Low, M.-M. (2019). Simulating an isolated dwarf galaxy with multichannel feedback and chemical yields from individual stars. Month. Notices R. Astron. Soc. 482, 1304–1329. doi: 10.1093/mnras/sty2689
Ercolano, B., Barlow, M. J., Storey, P. J., and Liu, X. W. (2003). MOCASSIN: a fully three-dimensional Monte Carlo photoionization code. Month. Notice R. Astron. Soc. 340, 1136–1152. doi: 10.1046/j.1365-8711.2003.06371.x
Federrath, C., Banerjee, R., Clark, P., and Klessen, R. (2010). Modeling collapse and accretion in turbulent gas clouds: implementation and comparison of sink particles in AMR and SPH. Astrophys. J. 713, 269–290. doi: 10.1088/0004-637X/713/1/269
Federrath, C., Banerjee, R., Seifried, D., Clark, P. C., and Klessen, R. S. (2011a). “Implementing and comparing sink particles in AMR and SPH,” in Computational Star Formation, Proceedings of the International Astronomical Union, IAU Symposium, vol. 270, 425–428. doi: 10.1017/S1743921311000755
Federrath, C., Chabrier, G., Schober, J., Banerjee, R., Klessen, R. S., and Schleicher, D. R. (2011b). Mach number dependence of turbulent magnetic field amplification: solenoidal versus compressive flows. Phys. Rev. Lett. 107:114504. doi: 10.1103/PhysRevLett.107.114504
Federrath, C., and Klessen, R. S. (2012). The star formation rate of turbulent magnetized clouds: comparing theory, simulations, and observations. Astrophys. J. 761:156. doi: 10.1088/0004-637X/761/2/156
Federrath, C., Rathborne, J. M., Longmore, S. N., Kruijssen, J. M. D., Bally, J., Contreras, Y., et al. (2017). The link between solenoidal turbulence and slow star formation in G0.253+0.016. Multi-Messenger Astrophys. Galactic Centre 322, 123–128. doi: 10.1017/S1743921316012357
Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., and Klessen, R. S. (2011c). A new jeans resolution criterion for (M)HD simulations of self-gravitating gas: application to magnetic field amplification by gravity-driven turbulence. Astrophys. J. 731:62. doi: 10.1088/0004-637X/731/1/62
Fontani, F., Commerçon, B., Giannetti, A., Beltrán, M. T., Sánchez-Monge, Á., Testi, L., et al. (2018). Fragmentation properties of massive protocluster gas clumps: an ALMA study. Astron. Astrophys. 615:A94. doi: 10.1051/0004-6361/201832672
Forgan, D., Rice, K., Stamatellos, D., and Whitworth, A. (2009). Introducing a hybrid radiative transfer method for smoothed particle hydrodynamics. Month. Notices R. Astron. Soc. 394, 882–891. doi: 10.1111/j.1365-2966.2008.14373.x
Fromang, S., Hennebelle, P., and Teyssier, R. (2006). A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics. Astron. Astrophys. 457, 371–384. doi: 10.1051/0004-6361:20065371
Fromang, S., Papaloizou, J., Lesur, G., and Heinemann, T. (2007). MHD simulations of the magnetorotational instability in a shearing box with zero net flux. Astron. Astrophys. 476, 1123–1132. doi: 10.1051/0004-6361:20077943
Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., et al. (2000). FLASH: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. Astrophys. J. Suppl. Ser. 131, 273–334. doi: 10.1086/317361
Gatto, A., Walch, S., Naab, T., Girichidis, P., Wünsch, R., Glover, S. C. O., et al. (2017). The SILCC project - III. Regulation of star formation and outflows by stellar winds and supernovae. Month. Notices R. Astron. Soc. 466, 1903–1924. doi: 10.1093/mnras/stw3209
Gavagnin, E., Bleuler, A., Rosdahl, J., and Teyssier, R. (2017). Star cluster formation in a turbulent molecular cloud self-regulated by photoionization feedback. Month. Notices R. Astron. Soc. 472, 4155–4172. doi: 10.1093/mnras/stx2222
Geen, S., Hennebelle, P., Tremblin, P., and Rosdahl, J. (2015). Photoionization feedback in a self-gravitating, magnetized, turbulent cloud. Month. Notices R. Astron. Soc. 454, 4484–4502. doi: 10.1093/mnras/stv2272
Geen, S., Hennebelle, P., Tremblin, P., and Rosdahl, J. (2016). Feedback in Clouds II: UV photoionization and the first supernova in a massive cloud. Month. Notices R. Astron. Soc. 463, 3129–3142. doi: 10.1093/mnras/stw2235
Geen, S., Watson, S. K., Rosdahl, J., Bieri, R., Klessen, R. S., and Hennebelle, P. (2018). On the indeterministic nature of star formation on the cloud scale. Month. Notices R. Astron. Soc. 481, 2548–2569. doi: 10.1093/mnras/sty2439
Gerin, M., Pety, J., Fuente, A., Cernicharo, J., Commerçon, B., and Marcelino, N. (2015). Nascent bipolar outflows associated with the first hydrostatic core candidates Barnard 1b-N and 1b-S. Astron. Astrophys. 577:L2. doi: 10.1051/0004-6361/201525777
Girichidis, P., Federrath, C., Allison, R., Banerjee, R., and Klessen, R. S. (2012a). Importance of the initial conditions for star formation - III. Statistical properties of embedded protostellar clusters. Month. Notices R. Astron. Soc. 420, 3264–3280. doi: 10.1111/j.1365-2966.2011.20250.x
Girichidis, P., Federrath, C., Banerjee, R., and Klessen, R. S. (2011). Importance of the initial conditions for star formation - I. Cloud evolution and morphology. Month. Notices R. Astron. Soc. 413, 2741–2759. doi: 10.1111/j.1365-2966.2011.18348.x
Girichidis, P., Federrath, C., Banerjee, R., and Klessen, R. S. (2012b). Importance of the initial conditions for star formation - II. Fragmentation-induced starvation and accretion shielding. Month. Notices R. Astron. Soc. 420, 613–626. doi: 10.1111/j.1365-2966.2011.20073.x
Girichidis, P., Seifried, D., Naab, T., Peters, T., Walch, S., Wünsch, R., et al. (2018). The SILCC project - V. The impact of magnetic fields on the chemistry and the formation of molecular clouds. Month. Notices R. Astron. Soc. 480, 3511–3540. doi: 10.1093/mnras/sty2016
Girichidis, P., Walch, S., Naab, T., Gatto, A., Wünsch, R., Glover, S. C. O., et al. (2016). The SILCC (SImulating the lifeCycle of molecular Clouds) project - II. Dynamical evolution of the supernova-driven ISM and the launching of outflows. Month. Notices R. Astron. Soc. 456, 3432–3455. doi: 10.1093/mnras/stv2742
González, M., Vaytet, N., Commercon, B., and Masson, J. (2015). Multigroup radiation hydrodynamics with flux-limited diffusion and adaptive mesh refinement. Astron. Astrophys. 578:A12. doi: 10.1051/0004-6361/201525971
González-Morales, P., Khomenko, E., Downes, T., and De Vicente, A. (2018). MHDSTS: a new explicit numerical scheme for simulations of partially ionised solar plasma. Astron. Astrophys. 615:67. doi: 10.1051/0004-6361/201731916
Goodwin, S. P., Whitworth, A. P., and Ward-Thompson, D. (2004). Astronomy & Astrophysics Simulating star formation in molecular cloud cores I. The influence of low levels of turbulence on fragmentation and multiplicity. Astron. Astrophys. 414, 633–650. doi: 10.1051/0004-6361:20031594
Greif, T. H., Springel, V., White, S. D., Glover, S. C., Clark, P. C., Smith, R. J., et al. (2011). Simulations on a moving mesh: the clustered formation of population III protostars. Astrophys. J. 737:75. doi: 10.1088/0004-637X/737/2/75
Guillet, T., Pakmor, R., Springel, V., Chandrashekar, P., and Klingenberg, C. (2018). High-order magnetohydrodynamics for astrophysics with an adaptive mesh refinement discontinuous galerkin scheme. arXiv. doi: 10.1093/mnras/stz314
Harries, T. J. (2015). Radiation-hydrodynamical simulations of massive star formation using Monte Carlo radiative transfer - I. Algorithms and numerical methods. Month. Notices R. Astron. Soc. 448, 3156–3166. doi: 10.1093/mnras/stv158
Hayes, J. C., Norman, M. L., Fiedler, R. A., Bordner, J. O., Li, P. S., Clark, S. E., et al. (2006). Simulating radiating and magnetized flows in multiple dimensions with ZEUS-MP. Astrophys. J. Suppl. Ser. 165, 188–228. doi: 10.1086/504594
Hennebelle, P., Commerçon, B., Joos, M., Klessen, R. S., Krumholz, M., Tan, J. C., et al. (2011). Collapse, outflows and fragmentation of massive, turbulent and magnetized prestellar barotropic cores. Astron. Astrophys. 528:A72. doi: 10.1051/0004-6361/201016052
Hopkins, P. F., Kereš, D., Oñorbe, J., Faucher-Giguère, C.-A., Quataert, E., Murray, N., et al. (2014). Galaxies on FIRE (feedback in realistic environments): stellar feedback explains cosmologically inefficient star formation. Month. Notices R. Astron. Soc. 445, 581–603. doi: 10.1093/mnras/stu1738
Hopkins, P. F., Quataert, E., and Murray, N. (2012). The structure of the interstellar medium of star-forming galaxies. Month. Notices R. Astron. Soc. 421, 3488–3521. doi: 10.1111/j.1365-2966.2012.20578.x
Hopkins, P. F., Wetzel, A., Kereš, D., Faucher-Giguère, C.-A., Quataert, E., Boylan-Kolchin, M., et al. (2018). FIRE-2 simulations: physics versus numerics in galaxy formation. Month. Notices R. Astron. Soc. 480, 800–863. doi: 10.1093/mnras/sty1690
Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., and Clark, P. C. (2016). Star formation and molecular hydrogen in dwarf galaxies: a non-equilibrium view. Month. Notices R. Astron. Soc. 458, 3528–3553. doi: 10.1093/mnras/stw544
Hubber, D. A., Batty, C. P., McLeod, A., and Whitworth, A. P. (2011). SEREN - a new SPH code for star and planet formation simulations. Algorithms and tests. Astron. Astrophys. 529:A27. doi: 10.1051/0004-6361/201014949
Iffrig, O., and Hennebelle, P. (2017). Structure distribution and turbulence in self-consistently supernova-driven ISM of multiphase magnetized galactic discs. arXiv 70, 1–22. doi: 10.1051/0004-6361/201630290
Inoue, T., and Inutsuka, S. (2008). Two-fluid magnetohydrodynamic simulations of converging H I flows in the interstellar medium. I. Methodology and basic results. Astrophys. J. 687, 303–310. doi: 10.1086/590528
Jappsen, A. K., Klessen, R. S., Larson, R. B., Li, Y., and Mac Low, M. M. (2005). The stellar mass spectrum from non-isothermal gravoturbulent fragmentation. Astron. Astrophys. 435, 611–623. doi: 10.1051/0004-6361:20042178
Jiang, Y.-F., Stone, J. M., and Davis, S. W. (2012). A godunov method for multidimensional radiation magnetohydrodynamics based on a variable eddington tensor. Astrophys. J. Suppl. 199:14. doi: 10.1088/0067-0049/199/1/14
Jones, M. O., and Bate, M. R. (2018). Sink particle radiative feedback in smoothed particle hydrodynamics models of star formation. Month. Notices R. Astron. Soc. 480, 2562–2577. doi: 10.1093/mnras/sty1969
Kim, C.-G., Kim, W.-T., and Ostriker, E. C. (2011). Regulation of star formation rates in multiphase galactic disks: numerical tests of the thermal/dynamical equilibrium model. Astrophys. J. 743:25. doi: 10.1088/0004-637X/743/1/25
Kim, C.-G., and Ostriker, E. C. (2015). Vertical equilibrium, energetics, and star formation rates in magnetized galactic disks regulated by momentum feedback from supernovae. Astrophys. J. 815:67. doi: 10.1088/0004-637X/815/1/67
Kim, C.-G., and Ostriker, E. C. (2017). Three-phase interstellar medium in galaxies resolving evolution with star formation and supernova feedback (TIGRESS): algorithms, fiducial model, and convergence. Astrophys. J. 846:133. doi: 10.3847/1538-4357/aa8599
Kim, C.-G., and Ostriker, E. C. (2018). Numerical simulations of multiphase winds and fountains from star-forming galactic disks. I. Solar neighborhood TIGRESS model. Astrophys. J. 853:173. doi: 10.3847/1538-4357/aaa5ff
Kim, C.-G., Ostriker, E. C., and Kim, W.-T. (2013). Three-dimensional hydrodynamic simulations of multiphase galactic disks with star formation Feedback. I. Regulation of star formation rates. Astrophys. J. 776:1. doi: 10.1088/0004-637X/776/1/1
Kim, W.-T., Ostriker, E. C., and Stone, J. M. (2002). Three-dimensional simulations of parker, magneto-jeans, and swing instabilities in shearing galactic gas disks. Astrophys. J. 581, 1080–1100. doi: 10.1086/344367
Kitsionas, S., and Whitworth, A. P. (2002). Smoothed particle hydrodynamics with particle splitting, applied to self-gravitating collapse. Month. Notices R. Astron. Soc. 330, 129–136. doi: 10.1046/j.1365-8711.2002.05115.x
Klassen, M., Kuiper, R., Pudritz, R. E., Peters, T., Banerjee, R., and Buntemeyer, L. (2014). A general hybrid radiation transport scheme for star formation simulations on an adaptive grid. Astrophys. J. 797:4. doi: 10.1088/0004-637X/797/1/4
Klassen, M., Pudritz, R., Kuiper, R., Peters, T., and Banerjee, R. (2016). Simulating the formation of massive protostars. I. Radiative feedback and accretion disks. Astrophys. J. 823:28. doi: 10.3847/0004-637X/823/1/28
Klingenberg, C., Pörner, F., and Xia, Y. (2017). An efficient implementation of the divergence free constraint in a discontinuous galerkin method for magnetohydrodynamics on unstructured meshes. Commun. Comput. Phys. 21, 423–442. doi: 10.4208/cicp.180515.230616a
Koga, S., Tsukamoto, Y., Okuzumi, S., and Machida, M. N. (2019). Dependence of Hall coefficient on grain size and cosmic ray rate and implication for circumstellar disc formation. Month. Notices R. Astron. Soc. 484, 2119–2136. doi: 10.1093/mnras/sty3524
Korpi, M. J., Brandenburg, A., Shukurov, A., Tuominen, I., and Nordlund, Å. (1999). A supernova-regulated interstellar medium: simulations of the turbulent multiphase medium. Astrophys. J. 514, L99–L102. doi: 10.1086/311954
Kotarba, H., Lesch, H., Dolag, K., Naab, T., Johansson, P. H., and Stasyszyn, F. A. (2009). Magnetic field structure due to the global velocity field in spiral galaxies. Month. Notices R. Astron. Soc. 397, 733–747. doi: 10.1111/j.1365-2966.2009.15030.x
Kraljic, K., Renaud, F., Bournaud, F., Combes, F., Elmegreen, B., Emsellem, E., et al. (2014). The role of turbulence in star formation laws and thresholds. Astrophys. J. 784:112. doi: 10.1088/0004-637X/784/2/112
Krumholz, M. R., Klein, R. I., and McKee, C. F. (2007a). Radiation-hydrodynamic simulations of collapse and fragmentation in massive protostellar cores. Astrophys. J. 656, 959–979. doi: 10.1086/510664
Krumholz, M. R., Klein, R. I., and McKee, C. F. (2012). Radiation-hydrodynamic simulations of the formation of orion-like star clusters. II. The initial mass function from winds, turbulence, and radiation. Astrophys. J. 754:71. doi: 10.1088/0004-637X/754/1/71
Krumholz, M. R., Klein, R. I., McKee, C. F., and Bolstad, J. (2007b). Equations and algorithms for mixed-frame flux-limited diffusion radiation hydrodynamics. Astrophys. J. 667, 626–643. doi: 10.1086/520791
Kuffmeier, M., Frimann, S., Jensen, S. S., and Haugbølle, T. (2018). Episodic accretion: the interplay of infall and disc instabilities. Month. Notices R. Astron. Soc. 475, 2642–2658. doi: 10.1093/mnras/sty024
Kuiper, R., and Hosokawa, T. (2018). First hydrodynamics simulations of radiation forces and photoionization feedback in massive star formation. Astron. Astrophys. 616:A101. doi: 10.1051/0004-6361/201832638
Kuiper, R., Klahr, H., Beuther, H., and Henning, T. (2010a). Circumventing the radiation pressure barrier in the formation of massive stars via disk accretion. Astrophys. J. 722, 1556–1576. doi: 10.1088/0004-637X/722/2/1556
Kuiper, R., Klahr, H., Beuther, H., and Henning, T. (2011). Three-dimensional simulation of massive star formation in the disk accretion scenario. Astrophys. J. 732:20. doi: 10.1088/0004-637X/732/1/20
Kuiper, R., Kuiper, R., Klahr, H., Klahr, H., Dullemond, C., Dullemond, C., et al. (2010b). Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation. Astron. Astrophys. 511:81. doi: 10.1051/0004-6361/200912355
Kuiper, R., Turner, N. J., and Yorke, H. W. (2016). Protostellar outflows and radiative feedback from massive stars. II. Feedback, star-formation efficiency, and outflow broadening. Astrophys. J. 832:40. doi: 10.3847/0004-637X/832/1/40
Kunz, M. W., and Mouschovias, T. C. (2009). The nonisothermal stage of magnetic star formation. I. Formulation of the problem and method of solution. Astrophys. J. 693, 1895–1911. doi: 10.1088/0004-637X/693/2/1895
Lee, Y.-N., and Hennebelle, P. (2018b). Stellar mass spectrum within massive collapsing clumps. II. Thermodynamics and tidal forces of the first Larson core. A robust mechanism for the peak of the IMF. Astron. Astrophys. 611:A89. doi: 10.1051/0004-6361/201731523
Lesur, G., Kunz, M. W., and Fromang, S. (2014). Thanatology in protoplanetary discs. The combined influence of Ohmic, Hall, and ambipolar diffusion on dead zones. Astron. Astrophys. 566:A56. doi: 10.1051/0004-6361/201423660
Li, P. S., Klein, R. I., and McKee, C. F. (2018). Formation of stellar clusters in magnetized, filamentary infrared dark clouds. Month. Notices R. Astron. Soc. 473, 4220–4241. doi: 10.1093/mnras/stx2611
Li, P. S., Martin, D. F., Klein, R. I., and McKee, C. F. (2012). A stable, accurate methodology for high mach number, strong magnetic field MHD turbulence with adaptive mesh refinement: resolution and refinement studies. Astrophys. J. 745:139. doi: 10.1088/0004-637X/745/2/139
Li, Z.-Y., Krasnopolsky, R., Shang, H., and Zhao, B. (2014). On the role of pseudodisk warping and reconnection in protostellar disk formation in turbulent magnetized cores. Astrophys. J. 793:130. doi: 10.1088/0004-637X/793/2/130
Londrillo, P., and del Zanna, L. (2004). On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method. J. Comput. Phys. 195, 17–48. doi: 10.1016/j.jcp.2003.09.016
Lorén-Aguilar, P., and Bate, M. R. (2014). Two-fluid dust and gas mixtures in smoothed particle hydrodynamics: a semi-implicit approach. Month. Notices R. Astron. Soc. 443, 927–945. doi: 10.1093/mnras/stu1173
Machida, M. N., Inutsuka, S.-i., and Matsumoto, T. (2006). Second core formation and high-speed jets: resistive magnetohydrodynamic nested grid simulations. Astrophys. J. Lett. 2, 151–154. doi: 10.1086/507179
Machida, M. N., Inutsuka, S.-i., and Matsumoto, T. (2014). Conditions for circumstellar disc formation: effects of initial cloud configuration and sink treatment. Month. Notices R. Astron. Soc. 438, 2278–2306. doi: 10.1093/mnras/stt2343
Marchand, P., Commerçon, B., and Chabrier, G. (2018). Impact of the Hall effect in star formation and the issue of angular momentum conservation. arXiv:1808.08731, 37. doi: 10.1051/0004-6361/201832907
Marchand, P., Masson, J., Chabrier, G., Hennebelle, P., Commerçon, B., and Vaytet, N. (2016). Chemical solver to compute molecule and grain abundances and non-ideal MHD resistivities in prestellar core-collapse calculations. Astron. Astrophys. 592:A18. doi: 10.1051/0004-6361/201526780
Marinacci, F., Vogelsberger, M., Kannan, R., Mocz, P., Pakmor, R., and Springel, V. (2018). Non-ideal magnetohydrodynamics on a moving mesh. Month. Notices R. Astron. Soc. 476, 2476–2492. doi: 10.1093/mnras/sty397
Martizzi, D., Fielding, D., Faucher-Giguère, C.-A., and Quataert, E. (2016). Supernova feedback in a local vertically stratified medium: interstellar turbulence and galactic winds. Month. Notices R. Astron. Soc. 459, 2311–2326. doi: 10.1093/mnras/stw745
Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., and Commerçon, B. (2016). Ambipolar diffusion in low-mass star formation. I. General comparison with the ideal magnetohydrodynamic case. Astron. Astrophys. 587:A32. doi: 10.1051/0004-6361/201526371
Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., and Chabrier, G. (2012). Incorporating ambipolar and ohmic diffusion in the AMR MHD code RAMSES. Astrophys. J. Suppl. 201:24. doi: 10.1088/0067-0049/201/2/24
Masunaga, H., and Inutsuka, S.-i. (2000). A radiation hydrodynamic model for protostellar collapse. II. The second collapse and the birth of a protostar. Astrophys. J. 531, 350–365. doi: 10.1086/308439
Matsumoto, T., Machida, M. N., and ichiro Inutsuka, S. (2017). Circumstellar disks and outflows in turbulent molecular cloud cores: possible formation mechanism for misaligned systems. Astrophys. J. 839:69. doi: 10.3847/1538-4357/aa6a1c
Maureira, M. J., Arce, H. G., Dunham, M. M., Pineda, J. E., Fernández-López, M., Chen, X., et al. (2017). Kinematics of a young low-mass star-forming core: understanding the evolutionary state of the first-core candidate L1451-mm. Astrophys. J. 838:60. doi: 10.3847/1538-4357/838/1/60
Mayer, L., Lufkin, G., Quinn, T., and Wadsley, J. (2007). Fragmentation of gravitationally unstable gaseous protoplanetary disks with radiative transfer. Astrophys. J. 661, L77–L80. doi: 10.1086/518433
Mellema, G., Iliev, I. T., Alvarez, M. A., and Shapiro, P. R. (2006). C 2-ray: a new method for photon-conserving transport of ionizing radiation. New Astron. 11, 374–395. doi: 10.1016/j.newast.2005.09.004
Meyer, C. D., Balsara, D. S., and Aslam, T. D. (2012). A second-order accurate Super TimeStepping formulation for anisotropic thermal conduction. Month. Notices R. Astron. Soc. 422, 2102–2115. doi: 10.1111/j.1365-2966.2012.20744.x
Meyer, C. D., Balsara, D. S., and Aslam, T. D. (2014). A stabilized Runge-Kutta-Legendre method for explicit super-time-stepping of parabolic and mixed equations. J. Comput. Phys. 257, 594–626. doi: 10.1016/j.jcp.2013.08.021
Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., and Ferrari, A. (2007). PLUTO: a numerical code for computational astrophysics. Astrophys. J. Suppl. Ser. 170, 228–242. doi: 10.1086/513316
Mihalas, D., and Klein, R. I. (1982). On the solution of the time-dependent inertial-frame equation of radiative transfer in moving media to O(vc). J. Comput. Phys. 46, 97–137. doi: 10.1016/0021-9991(82)90007-9
Mocz, P., Burkhart, B., Hernquist, L., McKee, C. F., and Springel, V. (2017). Moving-mesh Simulations of Star-forming Cores in Magneto-gravo-turbulence. Astrophys. J. 838:40. doi: 10.3847/1538-4357/aa6475
Mocz, P., Vogelsberger, M., and Hernquist, L. (2014a). A constrained transport scheme for MHD on unstructured static and moving meshes. Month. Notices R. Astron. Soc. 442, 43–55. doi: 10.1093/mnras/stu865
Mocz, P., Vogelsberger, M., Sijacki, D., Pakmor, R., and Hernquist, L. (2014b). A discontinuous Galerkin method for solving the fluid and magnetohydrodynamic equations in astrophysical simulations. Month. Notices R. Astron. Soc. 437, 397–414. doi: 10.1093/mnras/stt1890
Myers, A. T., Klein, R. I., Krumholz, M. R., and McKee, C. F. (2014). Star cluster formation in turbulent, magnetized dense clumps with radiative and outflow feedback. Month. Notices R. Astron. Soc. 439, 3420–3438. doi: 10.1093/mnras/stu190
Myers, A. T., McKee, C. F., Cunningham, A. J., Klein, R. I., and Krumholz, M. R. (2013). The fragmentation of magnetized, massive star-forming cores with radiative feedback. Astrophys. J. 766:97. doi: 10.1088/0004-637X/766/2/97
Nayakshin, S., Cha, S.-H., and Hobbs, A. (2009). Dynamic Monte Carlo radiation transfer in SPH: radiation pressure force implementation. Month. Notices R. Astron. Soc. 397, 1314–1325. doi: 10.1111/j.1365-2966.2009.15091.x
Ocvirk, P., Gillet, N., Shapiro, P. R., Aubert, D., Iliev, I. T., Teyssier, R., et al. (2016). Cosmic Dawn (CoDa): the first radiation-hydrodynamics simulation of reionization and galaxy formation in the local universe. Month. Notices R. Astron. Soc. 463, 1462–1485. doi: 10.1093/mnras/stw2036
Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., and Timmes, F. (2011). Modules for experiments in stellar astrophysics (MESA). Astrophys. J. Suppl. Ser. 192:3. doi: 10.1088/0067-0049/192/1/3
Pelupessy, F. I., van Elteren, A., de Vries, N., McMillan, S. L. W., Drost, N., and Portegies Zwart, S. F. (2013). The astrophysical multipurpose software environment. Astron. Astrophys. 557:A84. doi: 10.1051/0004-6361/201321252
Peters, T., Banerjee, R., Klessen, R. S., and Mac Low, M.-M. (2011). The interplay of magnetic fields, fragmentation, and ionization feedback in high-mass star formation. Astrophys. J. 729:72. doi: 10.1088/0004-637X/729/1/72
Peters, T., Banerjee, R., Klessen, R. S., Mac Low, M.-M., Galván-Madrid, R., and Keto, E. R. (2010). H II regions: witnesses to massive star formation. Astrophys. J. 711, 1017–1028. doi: 10.1088/0004-637X/711/2/1017
Peters, T., Klaassen, P. D., Mac Low, M.-M., Schrön, M., Federrath, C., Smith, M. D., et al. (2014). Collective outflow from a small multiple stellar system. Astrophys. J. 788:14. doi: 10.1088/0004-637X/788/1/14
Peters, T., Naab, T., Walch, S., Glover, S. C. O., Girichidis, P., Pellegrini, E., et al. (2017). The SILCC project - IV. Impact of dissociating and ionizing radiation on the interstellar medium and Hα emission as a tracer of the star formation rate. Month. Notices R. Astron. Soc. 466, 3293–3308. doi: 10.1093/mnras/stw3216
Petkova, M., and Springel, V. (2011). A novel approach for accurate radiative transfer in cosmological hydrodynamic simulations. Month. Notices R. Astron. Soc. 415, 3731–3749. doi: 10.1111/j.1365-2966.2011.18986.x
Phillips, G. J., and Monaghan, J. J. (1985). A numerical method for three-dimensional simulations of collapsing, isothermal, magnetic gas clouds. Month. Notices R. Astron. Soc. 216, 883–895. doi: 10.1093/mnras/216.4.883
Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., and De Zeeuw, D. L. (1999). A solution-adaptive upwind scheme for ideal magnetohydrodynamics. J. Comput. Phys. 154, 284–309. doi: 10.1006/jcph.1999.6299
Price, D. J., and Bate, M. R. (2009). Inefficient star formation: the combined effects of magnetic fields and radiative feedback. Month. Notices R. Astron. Soc. 398, 33–46. doi: 10.1111/j.1365-2966.2009.14969.x
Price, D. J., and Laibe, G. (2015). A fast and explicit algorithm for simulating the dynamics of small dust grains with smoothed particle hydrodynamics. Month. Notices R. Astron. Soc. 451, 813–826. doi: 10.1093/mnras/stv996
Price, D. J., and Monaghan, J. J. (2005). Smoothed particle magnetohydrodynamics — III. multidimensional tests and the div B= 0 constraint. Month. Notices R. Astron. Soc. 364, 384–406. doi: 10.1111/j.1365-2966.2005.09576.x
Price, D. J., Wurster, J., Tricco, T. S., Nixon, C., Toupin, S., Pettitt, A., et al. (2017). Phantom: a smoothed particle hydrodynamics and magnetohydrodynamics code for astrophysics. arXiv:1702.03930. doi: 10.1017/pasa.2018.25
Ramsey, J. P., and Dullemond, C. P. (2015). Radiation hydrodynamics including irradiation and adaptive mesh refinement with AZEuS. I. Methods. Astron. Astrophys. 574:A81. doi: 10.1051/0004-6361/201424954
Renaud, F., Bournaud, F., Emsellem, E., Elmegreen, B., Teyssier, R., Alves, J., et al. (2013). A sub-parsec resolution simulation of the Milky Way: global structure of the interstellar medium and properties of molecular clouds. Month. Notices R. Astron. Soc. 436, 1836–1851. doi: 10.1093/mnras/stt1698
Rijkhorst, E. J., Plewa, T., Dubey, A., and Mellema, G. (2006). Hybrid characteristics: 3D radiative transfer for parallel adaptive mesh refinement hydrodynamics. Astron. Astrophys. 452, 907–920. doi: 10.1051/0004-6361:20053401
Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., and Teyssier, R. (2013). RAMSES-RT: radiation hydrodynamics in the cosmological context. Month. Notices R. Astron. Soc. 436, 2188–2231. doi: 10.1093/mnras/stt1722
Rosen, A. L., Krumholz, M. R., Oishi, J. S., Lee, A. T., and Klein, R. I. (2017). Hybrid adaptive ray-moment method (HARM2): a highly parallel method for radiation hydrodynamics on adaptive grids. J. Comput. Phys. 330, 924–942. doi: 10.1016/j.jcp.2016.10.048
Ryan, B. R., Gammie, C. F., Fromang, S., and Kestener, P. (2017). Resolution dependence of magnetorotational turbulence in the isothermal stratified shearing box. Astrophys. J. 840:6. doi: 10.3847/1538-4357/aa6a52
Sano, T., and Stone, J. M. (2002). The effect of the hall term on the nonlinear evolution of the magnetorotational instability. I. Local axisymmetric simulations. Astrophys. J. 570, 314–328. doi: 10.1086/339504
Seifried, D., Walch, S., Girichidis, P., Naab, T., Wünsch, R., Klessen, R. S., et al. (2017). SILCC-Zoom: the dynamic and chemical evolution of molecular clouds. Month. Notices R. Astron. Soc. 472, 4797–4818. doi: 10.1093/mnras/stx2343
Shestakov, A. I., and Offner, S. S. R. (2008). A multigroup diffusion solver using pseudo transient continuation for a radiation-hydrodynamic code with patch-based AMR. J. Comput. Phys. 227, 2154–2186. doi: 10.1016/j.jcp.2007.09.019
Shima, K., Tasker, E. J., and Habe, A. (2017). Does feedback help or hinder star formation? The effect of photoionization on star formation in giant molecular clouds. Month. Notices R. Astron. Soc. 467, 512–523. doi: 10.1093/mnras/stw3279
Shu, F., Najita, J., Ostriker, E., Wilkin, F., Ruden, S., and Lizano, S. (1994). Magnetocentrifugally driven flows from young stars and disks. 1: a generalized model. Astrophys. J. 429, 781–796. doi: 10.1086/174363
Sigalotti, L. D. G., Cruz, F., Gabbasov, R., Klapp, J., and Ramírez-Velasquez, J. (2018). From large-scale to protostellar disk fragmentation into close binary stars. Astrophys. J. 857:40. doi: 10.3847/1538-4357/aab619
Smith, M. C., Sijacki, D., and Shen, S. (2018). Supernova feedback in numerical simulations of galaxy formation: separating physics from numerics. Month. Notices R. Astron. Soc. 478, 302–331. doi: 10.1093/mnras/sty994
Spruit, H. C., Stehle, R., and Papaloizou, J. C. B. (1995). Interchange instability in and accretion disc with a poloidal magnetic field. Month. Notices R. Astron. Soc. 275, 1223–1231. doi: 10.1093/mnras/275.4.1223
Stamatellos, D., Whitworth, A., Bisbas, T., and Goodwin, S. (2007). Radiative transfer and the energy equation in SPH simulations of star formation. Astron. Astrophys. 475, 37–49. doi: 10.1051/0004-6361:20077373
Stone, J. M., Mihalas, D., and Norman, M. L. (1992a). ZEUS-2D: a radiation magnetohydrodynamics code for astrophysical flows in two space dimensions. III - The radiation hydrodynamic algorithms and tests. Astrophys. J. Suppl. Ser. 80, 819–845. doi: 10.1086/191682
Stone, J. M., Stone, J. M., and Norman, M. L. (1992b). ZEUS-2D: a radiation magnetohydrodynamics code for astrophysical flows in two space dimensions. II. The magnetohydrodynamic algorithms and tests. Astrophys. J. Suppl. 80:791. doi: 10.1086/191681
Sur, S., Schleicher, D. R. G., Banerjee, R., Federrath, C., and Klessen, R. S. (2010). The generation of strong magnetic fields during the formation of the first stars. Astrophys. J. Lett. 721, L134–L138. doi: 10.1088/2041-8205/721/2/L134
Takasao, S., Tomida, K., Iwasaki, K., and Suzuki, T. K. (2018). A three-dimensional simulation of a magnetized accretion disk: fast funnel accretion onto a weakly magnetized star. Astrophys. J. 857:4. doi: 10.3847/1538-4357/aab5b3
Tasker, E. J. (2011). Star formation in disk galaxies. II. The effect Of star formation and photoelectric heating on the formation and evolution of giant molecular clouds. Astrophys. J. 730:11. doi: 10.1088/0004-637X/730/1/11
Teyssier, R., Fromang, S., and Dormy, E. (2006). Kinematic dynamos using constrained transport with high order Godunov schemes and adaptive mesh refinement. J. Comput. Phys. 218, 44–67. doi: 10.1016/j.jcp.2006.01.042
Tomida, K., Okuzumi, S., and Machida, M. N. (2015). Radiation magnetohydrodynamic simulations of protostellar collapse: nonideal magnetohydrodynamic effects and early formation of circumstellar disks. Astrophys. J. 801:117. doi: 10.1088/0004-637X/801/2/117
Tomida, K., Tomisaka, K., Matsumoto, T., Hori, Y., Okuzumi, S., Machida, M. N., et al. (2013). Radiation magnetohydrodynamic simulations of protostellar collapse: protostellar core formation. Astrophys. J. 763:6. doi: 10.1088/0004-637X/763/1/6
Tomida, K., Tomisaka, K., Matsumoto, T., Ohsuga, K., MacHida, M. N., and Saigo, K. (2010). Radiation magnetohydrodynamics simulation of proto-stellar collapse: two-component molecular outflow. Astrophys. J. Lett. 714(1 Pt 2):58–63. doi: 10.1088/2041-8205/714/1/L58
Tóth, G., De Zeeuw, D. L., Gombosi, T. I., and Powell, K. G. (2006). A parallel explicit/implicit time stepping scheme on block-adaptive grids. J. Comput. Phys. 217, 722–758. doi: 10.1016/j.jcp.2006.01.029
Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H. I., Howell, L. H., and Greenough, J. A. (1997). The jeans condition: a new constraint on spatial resolution in simulations of isothermal self-gravitational hydrodynamics. Astrophys. J. Lett. 489:L179. doi: 10.1086/310975
Tsitali, A. E., Belloche, A., Commerçon, B., and Menten, K. M. (2013). The dynamical state of the first hydrostatic core candidate Chamaeleon-MMS1. Astron. Astrophys. 557:A98. doi: 10.1051/0004-6361/201321204
Tsukamoto, Y., Iwasaki, K., and Inutsuka, S.-i. (2013). An explicit scheme for ohmic dissipation with smoothed particle magnetohydrodynamics. Month. Notices R. Astron. Soc. 434, 2593–2599. doi: 10.1093/mnras/stt1205
Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., and Inutsuka, S. (2015a). Bimodality OF circumstellar disk evolution induced by the hall current. Astrophys. J. Lett. 810:L26. doi: 10.1088/2041-8205/810/2/L26
Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., and Inutsuka, S. (2015b). Effects of Ohmic and ambipolar diffusion on formation and evolution of first cores, protostars, and circumstellar discs. Month. Notices R. Astron. Soc. 452, 278–288. doi: 10.1093/mnras/stv1290
Tsukamoto, Y., Okuzumi, S., Iwasaki, K., Machida, M. N., and Inutsuka, S.-I. (2017). The impact of the Hall effect during cloud core collapse: implications for circumstellar disk evolution. Publ. Astron. Soc. Jpn. 69:95. doi: 10.1093/pasj/psx113
Vaidya, B., Prasad, D., Mignone, A., Sharma, P., and Rickler, L. (2017). Scalable explicit implementation of anisotropic diffusion with Runge- Kutta-Legendre super-time stepping. Month. Notices R. Astron. Soc. 472, 3147–3160. doi: 10.1093/mnras/stx2176
Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., and Masson, J. (2012). Simulations of protostellar collapse using multigroup radiation hydrodynamics: I. the first collapse. Astron. Astrophys. 543:60. doi: 10.1051/0004-6361/201219427
Vaytet, N., Chabrier, G., Audit, E., Commerçon, B., Masson, J., Ferguson, J., et al. (2013). Simulations of protostellar collapse using multigroup radiation hydrodynamics. II. The second collapse. Astron. Astrophys. 557:A90. doi: 10.1051/0004-6361/201321423
Walch, S., Girichidis, P., Naab, T., Gatto, A., Glover, S. C. O., Wünsch, R., et al. (2015). The SILCC (SImulating the LifeCycle of molecular Clouds) project - I. Chemical evolution of the supernova-driven ISM. Month. Notices R. Astron. Soc. 454, 238–268. doi: 10.1093/mnras/stv1975
Wall, J. E., McMillan, S. L. W., Mac Low, M.-M., Klessen, R. S., and Portegies Zwart, S. (2019). Collisional N-body dynamics coupled to self-gravitating magnetohydrodynamics reveals dynamical binary formation. arXiv:1901.01132.
Wang, P., Li, Z.-Y., Abel, T., and Nakamura, F. (2010). Outflow feedback regulated massive star formation in parsec-scale cluster-forming clumps. Astrophys. J. 709, 27–41. doi: 10.1088/0004-637X/709/1/27
Whitehouse, S. C., and Bate, M. R. (2004). Smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation. Month. Notices R. Astron. Soc. 353, 1078–1094. doi: 10.1111/j.1365-2966.2004.08131.x
Whitehouse, S. C., and Bate, M. R. (2006). The thermodynamics of collapsing molecular cloud cores using smoothed particle hydrodynamics with radiative transfer. Month. Notices R. Astron. Soc. 367, 32–38. doi: 10.1111/j.1365-2966.2005.09950.x
Whitehouse, S. C., Bate, M. R., and Monaghan, J. J. (2005). A faster algorithm for smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation. Month. Notices R. Astron. Soc. 364, 1367–1377. doi: 10.1111/j.1365-2966.2005.09683.x
Winkler, K. A., and Norman, M. L. (1986). “Astrophysical radiation hydrodynamics,” in NATO Advanced Research Workshop on Astrophysical Radiation Hydrodynamics (Dordrecht; Holland; Boston, MA: D. Reidel Publishing Co.), 187.
Wise, J. H., and Abel, T. (2011). ENZO+MORAY: radiation hydrodynamics adaptive mesh refinement simulations with adaptive ray tracing. Month. Notices R. Astron. Soc. 414, 3458–3491. doi: 10.1111/j.1365-2966.2011.18646.x
Wurster, J., Bate, M. R., and Price, D. J. (2018). The collapse of a molecular cloud core to stellar densities using radiation non-ideal magnetohydrodynamics. Month. Notices R. Astron. Soc. 475, 1859–1880. doi: 10.1093/mnras/stx3339
Yorke, H. W., and Bodenheimer, P. (1999). The formation of protostellar disks. III. The influence of gravitationally induced angular momentum transport on disk structure and appearance. Astrophys. J. 525, 330–342. doi: 10.1086/307867
Zamora-Avilés, M., Vázquez-Semadeni, E., Körtgen, B., Banerjee, R., and Hartmann, L. (2018). Magnetic suppression of turbulence and the star formation activity of molecular clouds. Month. Notices R. Astron. Soc. 474, 4824–4836. doi: 10.1093/mnras/stx3080
Zhang, S., Hartmann, L., Zamora-Avilés, M., and Kuznetsova, A. (2018). On estimating angular momenta of infalling protostellar cores from observations. Month. Notices R. Astron. Soc. 480, 5495–5503. doi: 10.1093/mnras/sty2244
Zhang, W., Howell, L., Almgren, A., Burrows, A., Dolence, J., and Bell, J. (2013). Castro: a new compressible astrophysical solver. III. Multigroup radiation hydrodynamics. Astrophys. J. Suppl. Ser. 204:7. doi: 10.1088/0067-0049/204/1/7
Zhao, B., Caselli, P., and Li, Z.-Y. (2018a). Effect of grain size on differential desorption of volatile species and on non-ideal MHD diffusivity. Month. Notices R. Astron. Soc. 478, 2723–2736. doi: 10.1093/mnras/sty1165
Zhao, B., Caselli, P., Li, Z.-Y., Krasnopolsky, R., Shang, H., and Nakamura, F. (2016). Protostellar disc formation enabled by removal of small dust grains. Month. Notices R. Astron. Soc. 460, 2050–2076. doi: 10.1093/mnras/stw1124
Zhao, B., Caselli, P., Li, Z. Y., and Krasnopolsky, R. (2018b). Decoupling of magnetic fields in collapsing protostellar envelopes and disc formation and fragmentation. Month. Notices R. Astron. Soc. 473, 4868–4889. doi: 10.1093/mnras/stx2617
Keywords: star formation, numerical techniques, MHD: ideal, MHD: non-ideal, astrophysical fluid dynamics, radiation fields, sink particles
Citation: Teyssier R and Commerçon B (2019) Numerical Methods for Simulating Star Formation. Front. Astron. Space Sci. 6:51. doi: 10.3389/fspas.2019.00051
Received: 22 January 2019; Accepted: 02 July 2019;
Published: 24 July 2019.
Edited by:Christopher F. McKee, University of California, Berkeley, United States
Reviewed by:Christoph Federrath, Australian National University, Australia
Richard Klein, University of California, Berkeley, United States
Copyright © 2019 Teyssier and Commerçon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Romain Teyssier, email@example.com