© U. Lebreuilly et al. 2019

1. Introduction

Dust grains are thought to represent on average only 1% of the mass of the diffuse interstellar medium (ISM, Mathis et al. 1977; Weingartner & Draine 2001). Their typical size distribution is usually modeled as a power law (Mathis et al. 1977, hereafter MRN). However, recent works indicate that this is probably not the case in the dense parts of the ISM (Pagani et al. 2010; Compiègne et al. 2011; Harsono et al. 2018). In addition, a dynamical size sorting may operate in molecular clouds (Hopkins & Lee 2016; Tricco et al. 2017) or during the protostellar collapse (Bate & Lorén-Aguilar 2017). This sorting leads to important variations in the dust-to-gas ratio, especially for large dust grains.

The dust budget plays an essential role during the star and disk formation processes. On the one hand, it regulates the thermal budget of the medium through its blackbody emission, its absorption, and the photoelectric effect. On the other hand, the charge on the grain surface controls, through their number, the ambipolar, ohmic and hall resistivities during the protostellar collapse (Marchand et al. 2016). Zhao et al. (2016) show that removing grains smaller than ∼100 nm from the dust distribution enhances the effects of ambipolar diffusion, triggering the formation of a rotationally supported disk. As we need to know the dust concentration, a treatment of its dynamics is therefore essential to improve our understanding of stellar, disk, and planet formation.

Bifluid models of strongly coupled dust and gas mixtures are difficult to integrate numerically. For Lagrangian-based methods they tend to produce spurious dust aggregates when the grains are accumulated below the resolution length of the gas (Ayliffe et al. 2012). They also require extremely high spatial resolution in order to capture correctly the dephasing between gas and dust, i.e., the dissipation of energy through the drag (Laibe & Price 2012). Some efficient bifluid treatments were developed to cope with these issues (e.g., semi implicit methods, Lorén-Aguilar & Bate 2014, 2015; Stoyanovskaya et al. 2017, 2018). As an alternative way to model correctly these mixtures, Laibe & Price (2014a) proposed an algorithm for smoothed particle hydrodynamic codes (SPH, Lucy 1977; Gingold & Monaghan 1977) using a monofluid formalism. It has been proven to be efficient in the so-called diffusion approximation valid for strongly coupled mixtures (Price & Laibe 2015). More recently, Hutchison et al. (2018) improved this method and adapted it to treat simultaneously several dust species. Lin & Youdin (2017) and Chen & Lin (2018) have applied the monofluid formalism on a cylindrical grid adapted to protoplanetary disk simulations with an Eulerian approach.

Adaptive mesh refinement (AMR, Berger & Oliger 1984) codes are a complementary approach to SPH in astrophysical or cosmological simulations. Commerçon et al. (2008) show in the context of protostellar collapse that the two methods lead to converged results, but AMR codes capture details more efficiently thanks to the versatility of the refinement criteria. The adaptive mesh is a high performance tool for star formation simulations where a wide range of physical and dynamical scales are considered (e.g., Hennebelle 2018; Vaytet et al. 2018). In this paper, we present a numerical Eulerian approach of dust dynamics in the diffusion approximation on AMR grids and its implementation in the RAMSES code (Teyssier 2002), as well as an application to the first collapse (Larson 1969).

The paper is organized as follows. In Sect. 2, we recall the theory of gas and dust mixtures and the diffusion approximation. Then the dust diffusion algorithm implemented in RAMSES is presented in Sect. 3 followed by the validation tests in Sect. 4. The application on the dynamics of dust during the first collapse in the spherical case is shown in Sect. 5. In Sect. 6, we present our conclusions and perspectives. A companion paper, paper II (Lebreuilly et al. in prep.), will discuss in detail the dust dynamics during the first collapse with initial solid-body rotation and turbulence, and with multiple grains species.

2. Gas and dust mixtures in the diffusion approximation

2.1. Grain dynamics

First, we examine the fate of a single grain within a gas of density ρ_g. Let us consider a grain of radius s_grain, mass m_grain, and intrinsic density ρ_grain. It exchanges momentum with the surrounding gas molecules through microscopic collisions. Macroscopically, this leads to a drag force on the grain that writes

$$\begin{array}{c}\hfill \mathbf{F}=-\frac{m_{\mathrm{grain}}}{t_{\mathrm{s},\mathrm{grain}}}\mathbf{\Delta }\mathbf{v},\end{array}$$(1)

where the stopping time of the grain t_s, grain, is the typical time required by the dust grain to adjust its dynamics to a change in the gas velocity, and Δv is the differential velocity between the grain and the gas. In astrophysical regimes the mean free path of the gas is much larger than the size of the grain. In this case the grain stopping time is given by (Epstein 1924)

$$\begin{array}{c}\hfill t_{\mathrm{s},\mathrm{grain}}\equiv \sqrt{\frac{\pi \gamma }{8}}\frac{{\rho }_{\mathrm{grain}}}{{\rho }_{\mathrm{g}}}\frac{s_{\mathrm{grain}}}{c_{\mathrm{s}}},\end{array}$$(2)

where c_s and γ denote the sound speed and the adiabatic index of the gas, respectively. This is a particular case of a linear drag regime.

Large and dense grains are less coupled to the gas. If c_s and ρ_g increase, the coupling of the two phases becomes more important as the collision rate between gas and dust particles increases. We note that the opposite drag force applies on the gas particles for momentum conservation.

2.2. The monofluid formalism

2.2.1. Gas and dust bifluids

A usual model for gas and dust mixtures consists of two fluids interacting via the drag term (Saffman 1962). Due to the momentum exchanges between grains and gas particles, collective effects of the dust grains play an important role in the mixture dynamics. A continuous fluid description of the dust is more practical than treating the grains independently although this description might be inappropriate for very large particles. The dust fluid is characterized by a density ρ_d and an advection velocity v_d. The dust density ρ_d must not be confused with the intrinsic grain density ρ_grain. This fluid is considered pressureless and inviscid since the collisions with other grains are much less frequent than with gas particles. Through the drag, dust feels the forces exerted on the gas. The surrounding gas is described, as usual, as a fluid with a density ρ_g and an advection velocity v_g.

The bifluid formulation of gas and dust mass and momentum conservation writes

$$\begin{array}{cc}& \frac{{\mathrm{d}}_{\mathrm{g}}{\rho }_{\mathrm{g}}}{\mathrm{d}t}=-{\rho }_{\mathrm{g}}(\mathbf{\nabla }·{\mathbf{v}}_{\mathrm{g}}),\hfill \\ \hfill & \frac{{\mathrm{d}}_{\mathrm{d}}{\rho }_{\mathrm{d}}}{\mathrm{d}t}=-{\rho }_{\mathrm{d}}(\mathbf{\nabla }·{\mathbf{v}}_{\mathrm{d}}),\hfill \\ \hfill & {\rho }_{\mathrm{g}}\frac{{\mathrm{d}}_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}}{\mathrm{d}t}=-\mathbf{\nabla }{P}_{\mathrm{g}}+{\rho }_{\mathrm{g}}{\mathbf{f}}_{\mathrm{g}}+{\rho }_{\mathrm{g}}\mathbf{f}+K\mathbf{\Delta }\mathbf{v},\hfill \\ \hfill & {\rho }_{\mathrm{d}}\frac{{\mathrm{d}}_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}}{\mathrm{d}t}={\rho }_{\mathrm{d}}{\mathbf{f}}_{\mathrm{d}}+{\rho }_{\mathrm{d}}\mathbf{f}-K\mathbf{\Delta }\mathbf{v},\hfill \\ \hfill & \frac{{\mathrm{d}}_{\mathrm{g}}{e}_{\mathrm{g}}}{\mathrm{d}t}=-\frac{P_{\mathrm{g}}}{{\rho }_{\mathrm{g}}}\mathrm{\nabla }·{\mathbf{v}}_{\mathrm{g}}+\frac{K}{{\rho }_{\mathrm{g}}}\mathbf{\Delta }\mathbf{v}·\mathbf{\Delta }\mathbf{v},\hfill \end{array}$$(3)

where Δv ≡ v_d − v_g is the differential velocity between the two fluids; f is the specific force acting on both gas and dust; and f_g and f_d are the specific forces, drag and pressure force excluded, affecting the gas or the dust, respectively. In the rest of the article, we assume the hydrodynamical case where f_g and f_d are zero, and f is either zero or the gravity force. The gas pressure P_g is given by P_g = (γ − 1)ρ_ge_g, e_g being the gas specific internal energy. The gas derivative $\frac{{\mathrm{d}}_{\mathrm{g}}}{\mathrm{d}t}\equiv \frac{\partial }{\partial t}+({\mathbf{v}}_{\mathrm{g}}·\mathrm{\nabla })$ differs from the dust derivative $\frac{{\mathrm{d}}_{\mathrm{d}}}{\mathrm{d}t}\equiv \frac{\partial }{\partial t}+({\mathbf{v}}_{\mathrm{d}}·\mathrm{\nabla })$. We also define the stopping time as

$$\begin{array}{c}\hfill t_{\mathrm{s}}\equiv \sqrt{\frac{\pi \gamma }{8}}\frac{{\rho }_{\mathrm{grain}}}{\rho }\frac{s_{\mathrm{grain}}}{c_{\mathrm{s}}}=\frac{{\rho }_{\mathrm{g}}}{\rho }{t}_{\mathrm{s},\mathrm{grain}},\end{array}$$(4)

and the drag coefficient K as

$$\begin{array}{c}\hfill K\equiv \frac{{\rho }_{\mathrm{d}}{\rho }_{\mathrm{g}}}{\rho t_{\mathrm{s}}}·\end{array}$$(5)

2.2.2. Two-phases mixture

Alternatively, it is possible to model the gas and dust mixture as a single fluid made of two phases of total density ρ ≡ ρ_d + ρ_g (Laibe & Price 2014a,b). It is moving at the mixture barycentric velocity v:

$$\begin{array}{c}\hfill \mathbf{v}\equiv \frac{{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}+{\rho }_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}}}·\end{array}$$

We also define the dust ratio

$$\begin{array}{c}\hfill ϵ\equiv \frac{{\rho }_{\mathrm{d}}}{\rho },\end{array}$$

which must not be confused with the dust-to-gas ratio $\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}$. Using these new quantities, it is straightforward to write

$$\begin{array}{cc}\hfill {\mathbf{v}}_{\mathrm{g}}=& \mathbf{v}-ϵ\mathbf{\Delta }\mathbf{v},\hfill \\ \hfill {\mathbf{v}}_{\mathrm{d}}=& \mathbf{v}+(1-ϵ)\mathbf{\Delta }\mathbf{v},\hfill \\ \hfill {\rho }_{\mathrm{g}}=& (1-ϵ)\rho ,\hfill \\ \hfill {\rho }_{\mathrm{d}}=& ϵ\rho .\hfill \end{array}$$

After some developments, we write the conservation of mass and momentum for the mixture and of the gas internal energy

$$\begin{array}{cc}& \frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho (\mathbf{\nabla }·\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }+\mathbf{f}-\frac{1}{\rho }\mathbf{\nabla }·(ϵ(1-ϵ)\rho \mathbf{\Delta }\mathbf{v}\otimes \mathbf{\Delta }\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}ϵ}{\mathrm{d}t}=-\frac{1}{\rho }\mathbf{\nabla }·(ϵ(1-ϵ)\rho \mathbf{\Delta }\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{\Delta }\mathbf{v}}{\mathrm{d}t}=\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{(1-ϵ)\rho }-\frac{\mathbf{\Delta }\mathbf{v}}{t_{\mathrm{s}}}-(\mathbf{\Delta }\mathbf{v}·\mathbf{\nabla })\mathbf{v}+\frac{1}{2}\mathrm{\nabla }((2ϵ-1)\mathbf{\Delta }\mathbf{v}·\mathbf{\Delta }\mathbf{v})\hfill \\ \hfill & +(1-ϵ)\mathbf{\Delta }\mathbf{v}\times (\mathrm{\nabla }\times (1-ϵ)\mathbf{\Delta }\mathbf{v})-ϵ\mathbf{\Delta }\mathbf{v}\times (\mathrm{\nabla }\times ϵ\mathbf{\Delta }\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}=-\frac{P_{\mathrm{g}}}{\rho (1-ϵ)}\mathrm{\nabla }·(\mathbf{v}-ϵ\mathrm{\Delta }\mathbf{v})+(ϵ\mathrm{\Delta }\mathbf{v}·\mathrm{\nabla })e_{\mathrm{g}}+ϵ\frac{\mathbf{\Delta }\mathbf{v}·\mathbf{\Delta }\mathbf{v}}{t_{\mathrm{s}}}·\hfill \end{array}$$(6)

Here only one Lagrangian derivative $\frac{\mathrm{d}}{\mathrm{d}t}\equiv \frac{\partial }{\partial t}+(\mathbf{v}·\mathrm{\nabla })$ has to be defined. This system must be closed by the gas equation of state. We note the presence of the last terms in the differential velocity equation that were omitted in previous studies.

At this point, the monofluid approach is a dual reformulation of the dust and gas fluid Eqs. (3). However, it may be much more convenient for numerical simulations since it requires only one resolution scale (the mixture resolution) and it avoids the artificial trapping of dust particles (Laibe & Price 2014a).

2.2.3. Diffusion approximation

Laibe & Price (2014a) show that, for a strong drag regime, i.e., when the stopping time is short compared with the dynamical timescale t_dyn of the system, two simplifications can be made.

In the so-called diffusion approximation, ||Δv²/v²||, ||Δv⊗Δv/v²||, ||Δv×(∇ × (1−ϵ)Δv)/v²|| and ||Δv×(∇ × ϵΔv)/v²|| are second-order in St ≡ (t_s/t_dyn), where St is the Stokes number¹. System (6) reduces to

$$\begin{array}{cc}& \frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho (\mathbf{\nabla }·\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }+\mathbf{f},\hfill \\ \hfill & \frac{\mathrm{d}ϵ}{\mathrm{d}t}=-\frac{1}{\rho }\mathbf{\nabla }·(\rho ϵ(\mathbf{\Delta }\mathbf{v}-ϵ\mathbf{\Delta }\mathbf{v})),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{\Delta }\mathbf{v}}{\mathrm{d}t}=\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{(1-ϵ)\rho }-\frac{\mathbf{\Delta }\mathbf{v}}{t_{\mathrm{s}}}-(\mathbf{\Delta }\mathbf{v}·\mathbf{\nabla })\mathbf{v},\hfill \\ \hfill & \frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}=-\frac{P_{\mathrm{g}}}{\rho (1-ϵ)}\mathrm{\nabla }·\mathbf{v}+(ϵ\mathrm{\Delta }\mathbf{v}·\mathrm{\nabla })e_{\mathrm{g}}.\hfill \end{array}$$(7)

The errors caused by the diffusion approximation are marginal as long as St ≪ 1. We note that the first-order term $-\frac{P_{\mathrm{g}}}{\rho (1-ϵ)}\mathrm{\nabla }·\left(ϵ\mathrm{\Delta }\mathbf{v}\right)$ in the energy equation also simplifies owing to energy conservation, as explained by Price & Laibe (2015) in their second footnote. Another approximation known as the terminal velocity approximation (Youdin & Goodman 2005; Chiang 2008) can be made when St ≪ 1. In this case $\left||\frac{\mathrm{d}\mathbf{\Delta }\mathbf{v}}{\mathrm{d}t}/\left(\frac{\mathbf{\Delta }\mathbf{v}}{t_{\mathrm{s}}}\right)|\right|$ and $\left||(\mathbf{\Delta }\mathbf{v}·\mathbf{\nabla })\mathbf{v}/\left(\frac{\mathbf{\Delta }\mathbf{v}}{t_{\mathrm{s}}}\right)|\right|$ are transitory terms, and are thus neglected. A consequence of the terminal velocity approximation for linear drag regimes is that the differential velocity depends directly on the force balance

$$\begin{array}{c}\hfill \mathbf{\Delta }\mathbf{v}=t_{\mathrm{s}}\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{(1-ϵ)\rho }·\end{array}$$

It should be noted that this expression can change if other forces apply on the dust or the gas, e.g., the radiation or the Lorentz forces. In the terminal velocity approximation, we finally obtain

$$\begin{array}{cc}& \frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho (\mathbf{\nabla }·\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }+\mathbf{f},\hfill \\ \hfill & \frac{\mathrm{d}ϵ}{\mathrm{d}t}=-\frac{1}{\rho }\mathbf{\nabla }·\left(ϵt_{\mathrm{s}}\mathbf{\nabla }{P}_{\mathrm{g}}\right),\hfill \\ \hfill & \frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}=-\frac{P_{\mathrm{g}}}{\rho (1-ϵ)}\mathrm{\nabla }·\mathbf{v}+(ϵt_{\mathrm{s}}\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{(1-ϵ)\rho }·\mathrm{\nabla })e_{\mathrm{g}}.\hfill \end{array}$$

The two first equations are identical to pure hydrodynamics when only gas is involved. The third equation describes the evolution of the dust ratio. Using a conservative formulation, it can be expressed as an advection equation

$$\begin{array}{c}\hfill \frac{\partial \rho ϵ}{\partial t}+\mathbf{\nabla }·\left[\rho ϵ(\mathbf{v}+\frac{t_{\mathrm{s}}\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho })\right]=0.\end{array}$$(8)

In this formulation, the equation is almost identical to mass conservation but with a different advection velocity due to the dephasing between the dust and the barycenter. The specific internal energy equation is similar to pure hydrodynamics with an additional term that accounts for the back-reaction of dust on the gas.

2.2.4. (𝒩 + 1) phase mixtures

In the diffuse interstellar medium a grain size distribution commonly modeled as a power law (Mathis et al. 1977) is observed. It is therefore more realistic to consider several dust species in order to reproduce the observed dust size distribution.

In this perspective, the previous monofluid formalism has been extended to (𝒩+1) phase mixtures with 𝒩 distinct dust species and a gas phase (Laibe & Price 2014c; Hutchison et al. 2018). They show that, in the diffusion approximation the monofluid set of equations writes

$$\begin{array}{cc}& \frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho (\mathbf{\nabla }·\mathbf{v}),\hfill \\ \hfill & \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }+\mathbf{f},\hfill \\ \hfill & \frac{\mathrm{d}{ϵ}_k}{\mathrm{d}t}=-\frac{1}{\rho }\mathbf{\nabla }·\left({ϵ}_k{T}_{\mathrm{s},k}\mathbf{\nabla }{P}_{\mathrm{g}}\right),\forall k\in [1,\mathcal{N}],\hfill \\ \hfill & \frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}=-\frac{P_{\mathrm{g}}}{\rho (1-\mathcal{E})}\mathrm{\nabla }·\mathbf{v}+(\mathcal{E}{\mathcal{T}}_{\mathrm{s}}\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{(1-\mathcal{E})\rho }·\mathrm{\nabla })e_{\mathrm{g}},\hfill \end{array}$$(9)

where ϵ_k is the dust ratio of the phase k and T_s, k is the effective stopping time of the dust phase k defined as

$$\begin{array}{c}\hfill T_{\mathrm{s},k}\equiv \frac{t_{\mathrm{s},k}}{1-{ϵ}_k}-{\displaystyle \sum _{l=1}^{\mathcal{N}}}\frac{{ϵ}_l}{1-{ϵ}_l}{t}_{\mathrm{s},l},\end{array}$$

where t_s, k is the individual stopping time of the phase k. The $\sum _{l=1}^{\mathcal{N}}\frac{{ϵ}_l}{1-{ϵ}_l}{t}_{\mathrm{s},l}$ term accounts for the interaction between dust species that is due to their cumulative back-reaction on the gas. We also introduces the total dust ratio

$$\begin{array}{c}\hfill \mathcal{E}\equiv {\displaystyle \sum _{l=1}^{\mathcal{N}}}{ϵ}_l,\end{array}$$

and the mean stopping time

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathrm{s}}\equiv \frac{1}{\mathcal{E}}{\displaystyle \sum _{l=1}^{\mathcal{N}}}{ϵ}_l{T}_{\mathrm{s},l}.\end{array}$$

The gas and dust densities are simply given by

$$\begin{array}{cc}& {\rho }_{\mathrm{g}}\equiv (1-\mathcal{E})\rho ,\hfill \\ \hfill & {\rho }_{\mathrm{d},k}\equiv {ϵ}_k\rho .\hfill \end{array}$$

When 𝒩 = 1, Eqs. (9) reduce to the formulation for a single dust phase.

Combining the different equations with mass conservation in the Lagrangian frame, co-moving with the barycenter leads to the formulation of the system of equations in a conservative form

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\left[\rho \mathbf{v}\right]=0,\hfill \\ \hfill & \frac{\partial \rho \mathbf{v}}{\partial t}+\mathbf{\nabla }·[P_{\mathrm{g}}\mathbb{I}+\rho (\mathbf{v}\otimes \mathbf{v})]=\rho \mathbf{f},\hfill \\ \hfill & \frac{\partial {\rho }_{\mathrm{d},k}}{\partial t}+\mathbf{\nabla }·\left[{\rho }_{\mathrm{d},k}(\mathbf{v}+\frac{T_{\mathrm{s},k}\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho })\right]=0,\hfill \\ \hfill & \frac{\partial E}{\partial t}+\mathbf{\nabla }·[(E+P_{\mathrm{g}})\mathbf{v}]=\mathrm{\nabla }·\left[\frac{\mathcal{E}{\mathcal{T}}_{\mathrm{s}}}{1-\mathcal{E}}\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }\frac{P_{\mathrm{g}}}{\gamma -1}\right],\hfill \end{array}$$(10)

where ρ_d, k is the density of the dust phase k, $E\equiv \frac{1}{2}\rho {\mathbf{v}}^2+\rho (1-\mathcal{E})e_{\mathrm{g}}$ is the total energy of the mixture and 𝕀 is the identity matrix.

3. Numerical methods

3.1. RAMSES

3.1.1. Basic features

Before we present our method, we recall the main features of the RAMSES code (Teyssier 2002) in order to facilitate the understanding of our implementation of dust dynamics in the diffusion approximation.

RAMSES is a finite volume Eulerian code that integrates the equation of hydrodynamics in their conservative form on an AMR grid (Berger & Oliger 1984). Among others applications, it has been extended to magnetohydrodynamics (Teyssier et al. 2006; Fromang et al. 2006; Masson et al. 2012; Marchand et al. 2018), radiation hydrodynamics (Commerçon et al. 2011, 2014; Rosdahl et al. 2013; Rosdahl & Teyssier 2015; González et al. 2015), and cosmic ray and anisotropic heat conduction (Dubois & Commerçon 2016).

For simplicity, the hydrodynamical version of RAMSES is presented here. It solves the Euler equations of a pure gas in their hyperbolic form

$$\begin{array}{c}\hfill \frac{\partial \mathbb{U}}{\partial t}+\mathbf{\nabla }·\mathbb{F}(\mathbb{U})=0,\end{array}$$(11)

the state 𝕌 vector and flux 𝔽 vector being

$$\begin{array}{cc}& \mathbb{U}\equiv ({\rho }_{\mathrm{g}},{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}},E_{\mathrm{g}}),\hfill \\ \hfill & \mathbb{F}(\mathbb{U})\equiv ({\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}},{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}\otimes {\mathbf{v}}_{\mathrm{g}}+P_{\mathrm{g}}\mathbb{I},{\mathbf{v}}_{\mathrm{g}}(E_{\mathrm{g}}+P_{\mathrm{g}})),\hfill \end{array}$$

where $E_{\mathrm{g}}\equiv \frac{1}{2}{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}^2+{\rho }_{\mathrm{g}}{e}_{\mathrm{g}}$ is the total energy of the gas.

3.1.2. Godunov scheme

Finite volume methods are based on the estimation of the average of 𝕌 over the cells. For a 1D problem integrated over time (for t ∈ [t₁, t₂]) and space (for x ∈ [x₁, x₂]), the previous system writes

$$\begin{array}{c}\hfill {\displaystyle {\int }_{x_1}^{x_2}}(\mathbb{U}(t_2)-\mathbb{U}(t_1))\mathrm{d}x+{\displaystyle {\int }_{t_1}^{t_2}}(\mathbb{F}(x_2)-\mathbb{F}(x_1))\mathrm{d}t=0.\end{array}$$(12)

Hence, the evolution of the state vector is constrained by the fluxes at the interfaces of the cells. RAMSES uses a second-order predictor-corrector Godunov method (Godunov 1959) to update its value.

At the timestep n for a cell i, any physical quantity A is discretized as

$$\begin{array}{c}\hfill A(x,t)\to A_i^n.\end{array}$$

Cell interfaces are denoted with half integer subscripts, e.g, i + 1/2 is to the interface between the cell i and i + 1. Similarly, half timesteps are also denoted with half integer superscripts, e.g, n + 1/2.

For a cell of length Δx and a timestep Δt, the scheme writes

$$\begin{array}{c}\hfill {\mathbb{U}}_i^{n+1}={\mathbb{U}}_i^n-({\mathbb{F}}_{i+1/2}^{n+1/2}-{\mathbb{F}}_{i-1/2}^{n+1/2})\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x},\end{array}$$(13)

where the discretized fluxes ${\mathbb{F}}_{i\pm 1/2}^{n+1/2}$, which are the solutions of the Riemann problem at the interfaces, are approximated with a MUSCL predictor-corrector scheme (van Leer 1979) that ensures the second-order accuracy in space and in time for the Godunov method. We note that, in RAMSESΔx = Δy = Δz.

3.1.3. AMR grid and adaptive timestep

The AMR grid enables an accurate description of the regions of interest in the simulation box. Cells are tagged with a refinement level ℓ with a value between ℓ_min for the coarsest cells and ℓ_max for the finest. The length of a cell of level ℓ is

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\ell }=\frac{L_{\mathrm{box}}}{2^{\ell }},\end{array}$$

where L_box is the size of the simulation box. For a uniform grid with a unique level ℓ, the effective resolution is simply 2^ℓ.

RAMSES uses an adaptive timestep to reduce the computation time while maintaining a good accuracy for the refined cells. When a level ℓ is updated with a timestep Δt_ℓ, the level ℓ + 1 is updated twice with two timesteps, Δt_{1, ℓ + 1} and Δt_{2, ℓ + 1}, each verifying the CourantFriedrichsLewy condition (Courant et al. 1928, hereafter CFL). The two levels are synchronized with the following condition

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\ell }=\mathrm{\Delta }{t}_{1,\ell +1}+\mathrm{\Delta }{t}_{2,\ell +1}.\end{array}$$

The state vector is updated from fine to coarse levels. For a cell of level ℓ, three types of interfaces are possible

• The fine-to-coarse interfaces, when the neighbor cell is at level ℓ − 1;

• the fine-fine interfaces, when the neighbor cell is also at level ℓ;

• the coarse-to-fine interfaces, when the neighbor cell is at level ℓ + 1.

In RAMSES, the levels are updated from fine to coarse. During the update of level ℓ + 1, the fluxes at the interface with level ℓ are used to update the level ℓ + 1 and stored for the later update of the level ℓ. During the update of the level, the state vector of the coarser levels is held constant. The fine-to-coarse and fine-fine fluxes are computed during the update of the level ℓ. The coarse neighbor cells at level ℓ − 1 are virtually refined to compute the flux at the boundaries with level ℓ.

3.2. Dust diffusion scheme

3.2.1. Operator splitting

In the previous section, we explain how RAMSES solves the equations of hydrodynamics in their conservative form. The monofluid formalism in the diffusion approximation has the same structure for the mass, momentum and energy conservation equation. A total of 𝒩 additional equations are required to follow the evolution of the dust ratios. It is still possible to write the system in a hyperbolic form, but we choose to decompose the flux in two distinct terms. The system now writes

$$\begin{array}{c}\hfill \frac{\partial \mathbb{U}}{\partial t}+\mathbf{\nabla }·{\mathbb{F}}_{\mathrm{H}}(\mathbb{U})+\mathbf{\nabla }·{\mathbb{F}}_{\mathrm{\Delta }}(\mathbb{U})=0.\end{array}$$(14)

The new state vector 𝕌 writes

$$\begin{array}{c}\hfill \mathbb{U}\equiv (\rho ,\rho \mathbf{v},E,{\rho }_{\mathrm{d},k});\end{array}$$

the flux 𝔽_H is similar to the flux of pure hydrodynamics and is given by

$$\begin{array}{c}\hfill {\mathbb{F}}_{\mathrm{H}}(\mathbb{U})\equiv (\rho \mathbf{v},\rho \mathbf{v}\otimes \mathbf{v}+P_{\mathrm{g}}\mathbb{I},\mathbf{v}(E+P_{\mathrm{g}}),{\rho }_{\mathrm{d},k}\mathbf{v});\end{array}$$

the flux 𝔽_Δ accounts for dust diffusion and is given by

$$\begin{array}{c}\hfill {\mathbb{F}}_{\mathrm{\Delta }}(\mathbb{U})\equiv (0,0,\frac{P_{\mathrm{g}}}{\gamma -1}{\mathbf{w}}_{\mathrm{g},\mathrm{b}},{\rho }_{\mathrm{d},k}{\mathbf{w}}_{\mathbf{k}}),\end{array}$$

where w_k, the differential advection velocity between the dust species k and the barycenter, is

$$\begin{array}{c}\hfill {\mathbf{w}}_{\mathbf{k}}\equiv \frac{T_{\mathrm{s},k}\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho };\end{array}$$(15)

and w_g, b is the differential velocity between the gas and the barycenter that writes

$$\begin{array}{c}\hfill {\mathbf{w}}_{\mathrm{g},\mathrm{b}}=-\frac{\mathcal{E}{\mathcal{T}}_{\mathrm{s}}}{1-\mathcal{E}}\frac{\mathbf{\nabla }{P}_{\mathrm{g}}}{\rho }\end{array}$$(16)

or

$$\begin{array}{c}\hfill {\mathbf{w}}_{\mathrm{g},\mathrm{b}}=-{\displaystyle \sum _k}\frac{{\rho }_{\mathrm{d},k}}{\rho -{\sum }_j{\rho }_{\mathrm{d},j}}{\mathbf{w}}_{\mathbf{k}}.\end{array}$$(17)

The classical second-order Godunov method presented above is used to update the state vector. An operator splitting method is performed to solve the system in two steps. In the so-called hydrodynamical step, we only compute 𝔽_H and update the state vector accordingly. An important difference between this step and the classical hydrodynamical version of RAMSES is that the fluid density ρ is different from the gas density. As a consequence, the pressure and the wave speeds must be computed using ρ_g = ρ(1 − ℰ) instead of ρ.

3.2.2. MUSCL scheme for dust diffusion/advection

The second step of the operator splitting consists in taking into account the second flux vector 𝔽_Δ. For simplicity, we now focus on the update of the dust density keeping in mind that the dust density and energy updates are done simultaneously. In this perspective, the following equation is solved

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{\rho }}_{\mathrm{d},k}}{\partial t}=-\mathbf{\nabla }·\left[{\stackrel{\sim }{\rho }}_{\mathrm{d},k}{\mathbf{w}}_{\mathbf{k}}\right],\end{array}$$(18)

where ${\stackrel{\sim }{\rho }}_{\mathrm{d},k}$ refers to ρ_d, k after its advection as a passive scalar at the velocity v. In the remainder of the section, we omit this symbol and the k index. A MUSCL predictor-corrector scheme is used to compute the dust diffusion fluxes ℱ(ρ_d)≡ρ_dw. The flux storage and the Godunov update are performed as they are the hydrodynamical step.

Apart from the dust densities and the energy of the mixture, the variables are constant during this step. In particular, the total density ρ is constant.

Differential advection velocity. During the diffusion step, we aim to compute the differential advection velocity ${\mathit{w}}_{i,j,k}^n$. To get its value, the pressure gradient is evaluated at the center of the cell i, j, k of length Δx. For example, its component in the x-direction, is

$$\begin{array}{c}\hfill {\mathbf{\nabla }}_x{P_{\mathrm{g}}}_{i,j,k}^n=\frac{{P_{\mathrm{g}}}_{i+1,j,k}^n-{P_{\mathrm{g}}}_{i-1,j,k}^n}{\mathrm{\Delta }{x}_{i+1,j,k}+\mathrm{\Delta }{x}_{i-1,j,k}},\end{array}$$(19)

where Δx_{i − 1, j, k} and Δx_{i + 1, j, k} denote the distance in the x-direction from the center of the cell to the center of the left and right neighbor cells, respectively. They take into account the level of the neighbor cell.

If the neighbor cell is coarser, we interpolate the pressure at the fine level using a slope limiter to compute the gradient. This method is similar to what is used in the hydrodynamical solver of RAMSES, with a loss of one order of accuracy at level interface. In certain conditions this method maintains a second-order accuracy of the scheme in presence of AMR (see Sect. 4.4). The case where the neighbor cell is at a finer level is never considered since fine-to-coarse fluxes are already computed during the update of the finer level. Finally, the expression of the x-component of ${\mathbf{w}}_{i,j,k}^n$ is given by

$$\begin{array}{c}\hfill {w_x}_{i,j,k}^n=\frac{{T_{\mathrm{s}}}_{i,j,k}^n{\mathbf{\nabla }}_x{P_{\mathrm{g}}}_{i,j,k}^n}{{\rho }_{i,j,k}^n}·\end{array}$$(20)

Predictor step. During the predictor step, the dust density is estimated at the left and right interfaces with a simple finite difference method. It uses slope limiters to preserve the total variation diminishing property of the scheme (TVD, Harten 1983). The centered value of the dust density is estimated at half timesteps as

$$\begin{array}{c}\hfill {{\rho }_{\mathrm{d}}}_{i,j,k}^{n+1/2}={{\rho }_{\mathrm{d}}}_{i,j,k}^n-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}{\displaystyle \sum _{\sigma =x,y,z}}({w_{\sigma }}_{i,j,k}^n{\mathrm{\Delta }}_{\sigma }{{\rho }_{\mathrm{d}}}_{i,j,k}^n+{{\rho }_{\mathrm{d}}}_{i,j,k}^n{\mathrm{\Delta }}_{\sigma }{w_{\sigma }}_{i,j,k}^n),\end{array}$$(21)

where ${\Delta }_{\sigma }{\rho }_{\text{d}}{}_{i,j,k}^n$ and ${\Delta }_{\sigma }{w}_{\sigma }{}_{i,j,k}^n$ are the TVD variations of ρ_d and w_σ in the direction σ. After the prediction of ${{\rho }_{\mathrm{d}}}_{i,j,k}^{n+1/2}$, we interpolate ρ_d and w_σ at the interfaces.

Let us consider the interface in the x-direction. The left and right interface values, denoted by the L and R subscripts, are given by

$$\begin{array}{cc}& {{\rho }_L}_{i,j,k}={{\rho }_{\mathrm{d}}}_{i,j,k}^{n+1/2}-\frac{{\mathrm{\Delta }}_x{{\rho }_{\mathrm{d}}}_{i,j,k}^n}{2}·\hfill \\ \hfill & {{\rho }_R}_{i,j,k}={{\rho }_{\mathrm{d}}}_{i,j,k}^{n+1/2}+\frac{{\mathrm{\Delta }}_x{{\rho }_{\mathrm{d}}}_{i,j,k}^n}{2},\hfill \\ \hfill & {{w_x}_L}_{i,j,k}={w_x}_{i,j,k}^n-\frac{{\mathrm{\Delta }}_x{w_x}_{i,j,k}^n}{2},\hfill \\ \hfill & {{w_x}_R}_{i,j,k}={w_x}_{i,j,k}^n+\frac{{\mathrm{\Delta }}_x{w_x}_{i,j,k}^n}{2}·\hfill \end{array}$$(22)

The bottom, top, back, and front states are estimated in a similar way. We do not perform a predictor step in time for the differential advection velocity because the current scheme is sufficient to get the second-order convergence (see Sect. 4.4).

Corrective step. The correction operation consists in computing the fluxes at the interface using the left and right predicted values, respectively. We consider the left interface of the cell i, j, k in the x-direction. To avoid issues with the velocity discontinuities at the interfaces (e.g., due to spurious pressure jumps, Sharma et al. 2009), we impose a unique advection velocity

$$\begin{array}{c}\hfill w_{i-1/2,j,k}^{n+1/2}=\frac{{w_L}_{i,j,k}^{n+1/2}+{w_R}_{i-1,j,k}^{n+1/2}}{2}·\end{array}$$(23)

This average does not appear to be critical for the second-order accuracy of the scheme but leads to a better convergence and is consistent with the upwind implementation in RAMSES. Our scheme uses the upwind method (Courant et al. 1952) to estimate the flux, sufficient to get the second-order accuracy in space. It writes as

$$\begin{array}{c}\hfill {\mathcal{F}}_{i-1/2,j,k}^{n+1/2}=\max (w_{i-1/2,j,k}^{n+1/2}{{\rho }_{\mathrm{d},L}}_{i,j,k}^{n+1/2},0)+\min (w_{i-1/2,j,k}^{n+1/2}{{\rho }_{\mathrm{d},R}}_{i-1,j,k}^{n+1/2},0)·\end{array}$$(24)

We note that several other approximate Riemann solvers can be used in our implementation as well, such as Lax–Wendroff (Lax & Wendroff 1960) or Harten-Lax-van Leer (Harten et al. 1983, hereafter HLL). The fluxes in the other directions ${\mathcal{F}}_{i,j-1/2,k}^{n+1/2}$ and ${\mathcal{F}}_{i,j,k-1/2}^{n+1/2}$ are estimated in a similar way. The dust density is finally updated according to

$$\begin{array}{cc}\hfill {{\rho }_{\mathrm{d}}}_{i,j,k}^{n+1}={{\rho }_{\mathrm{d}}}_{i,j,k}^n& -\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\mathcal{F}}_{i+1/2,j,k}^{n+1/2}-{\mathcal{F}}_{i-1/2,j,k}^{n+1/2})\hfill \\ \hfill & -\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\mathcal{F}}_{i,j+1/2,k}^{n+1/2}-{\mathcal{F}}_{i,j-1/2,k}^{n+1/2})\hfill \\ \hfill & -\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\mathcal{F}}_{i,j,k+1/2}^{n+1/2}-{\mathcal{F}}_{i,j,k-1/2}^{n+1/2}).\hfill \end{array}$$(25)

Energy. Similarly and simultaneously with the dust diffusion step, the energy is updated after the hydrodynamical step to account for the energy component of 𝔽_Δ. It is computed using the same scheme as the dust density with the velocity computed with Eq. (17) and the internal energy instead of the dust density. The fluxes are then added to the state vector similarly to Eq. (25) with the total energy instead of the dust density.

3.3. Time-stepping

Since the dust diffusion step consists in solving an advection equation explicitly, the scheme stability is achieved if the timestep Δt verifies the split CFL condition

$$\begin{array}{c}\hfill \mathrm{\Delta }t<\mathrm{\Delta }{t}_{\mathrm{dust}}\equiv C_{\mathrm{CFL}}\frac{\mathrm{\Delta }x}{{\sum }_{\sigma =x,y,z}|\max (w_{\sigma },w_{\mathrm{g},\mathrm{b}\sigma })|},\end{array}$$(26)

where C_CFL <  1 is a safety factor. In addition, the hydrodynamical step imposes another stability condition that writes

$$\begin{array}{c}\hfill \mathrm{\Delta }t<C_{\mathrm{CFL}}\frac{\mathrm{\Delta }x}{|c_{\mathrm{s}}|+{\sum }_{\sigma =x,y,z}\left|{\nu }_{\sigma }\right|},\end{array}$$(27)

where ν_σ is the mixture velocity in the direction σ. Dust diffusion is stable without intervention on the timestep as long as

$$\begin{array}{c}\hfill |c_{\mathrm{s}}|+{\displaystyle \sum _{\sigma =x,y,z}}\left|{\nu }_{\sigma }\right|>{\displaystyle \sum _{\sigma =x,y,z}}|\max (w_{\sigma },w_{\mathrm{g},\mathrm{b}\sigma })|.\end{array}$$(28)

If the former condition is not verified, which is the case when the pressure gradient is steep, dust diffusion constrains the timestep. In this case, we impose the dust timestep instead of the hydrodynamical one.

4. Validation tests

To benchmark the implementation of dust dynamics in RAMSES, we run the canonical tests for gas and dust mixtures, DUSTYDIFFUSE, DUSTYSHOCK, DUSTYWAVE and the SETTLING test. We also test this advection solver with the DUSTYADVECT test.

4.1. Dustyadvect

The scheme convergence and behavior at discontinuities is examined with 1D advection tests. In these DUSTYADVECT tests, the advection velocity w_x is constant and the hydrodynamical update deactivated. It consists in solving the following equation

$$\begin{array}{c}\hfill \frac{\partial {\rho }_{\mathrm{d}}}{\partial t}=-w_x\frac{\partial {\rho }_{\mathrm{d}}}{\partial x}·\end{array}$$(29)

Considering an initial condition ρ_d(x, 0)=f(x), f(x) having a period of the size of the box L = 1, the analytic solution at time t is simply given by ρ_d(x, t) = f(x − w_xt). In the remainder of the section, w_x = 1.

At first, four DUSTYADVECT tests are performed. We impose the initial function f₁, which writes ∀x ∈ [0,L]

$$\begin{array}{c}\hfill f_1\left(x\right)=\{\begin{array}{cc}0.1\hfill & \text{if}\frac{L}{4}<x<\frac{3L}{4},\hfill \\ 0.01\hfill & \text{otherwise.}\hfill \end{array}\end{array}$$(30)

In the first test, no predictor step is operated. Different slope limiters are then used for the three other tests in the predictor step; Minmod (Roe 1986), Van-Leer (van Leer 1974), and Superbee (Roe 1986). Uniform grids of resolutions ranging from ℓ = 7 (128 cells) to ℓ = 13 (8192 cells) are considered. Extremely small timesteps compared with the CFL condition Δt = 8 × 10⁻⁶ are imposed to ensure that the spatial error dominates.

Figure 1 shows the outcome of these tests after one period, at t = 1. As expected, the quality of the results strongly depends on the slope limiter. Without the predictor step (no slope limiter) the solver is simply a first-order centered upwind scheme. In this case a larger resolution is required to achieve the same accuracy as in the three other tests. As expected, the Van-Leer and Superbee slope limiters give more accurate results than the Minmod test but at the cost of a lack of symmetry. The Minmod slope limiter was therefore chosen as a good compromise to achieve a satisfying precision and to preserve symmetry.

Fig. 1. DUSTYADVECT tests using the function f1 (Eq. (30)). Dust density after one on period various grids (ℓ = 7 in blue, ℓ = 9 in orange, ℓ = 11 in green, ℓ = 13 in red) as a function of the position compared with the analytic solution (black solid line). We present this test using four different slope limiters. Top left: No predictor step. Top right: Minmod slope limiter. Bottom left: Van-Leer slope limiter. Bottom right: superbee slope limiter.

Another series of DUSTYADVECT tests are performed for different resolutions, using the Minmod slope limiter and without a predictor step. A smooth initial Gaussian-like function f₂ is imposed to quantify the truncation error in space of the scheme which writes, ∀x ∈ [0,L],

$$\begin{array}{c}\hfill f_2\left(x\right)=0.01+0.1exp\left[-{\left(\frac{x-L/2}{L/4}\right)}^2\right]·\end{array}$$(31)

We use a very small timestep Δt = 10⁻⁸ to minimize the truncation errors in time compared with spatial errors. The results are compared at the same time t = 0.01 with the analytic solution using the L₂ norm

$$\begin{array}{c}\hfill L_2=\sqrt{\frac{\sum _{i=1}^{N_{\mathrm{cell}}}{|{{\rho }_{\mathrm{d}}}_{\mathrm{RESULTS},i}^n-{{\rho }_{\mathrm{d}}}_{\mathrm{ANALYTIC},i}^n|}^2}{N_{\mathrm{cell}}}}·\end{array}$$(32)

Figure 2 shows the evolution of the L₂ norm with (blue squares) and without (orange diamonds) a prediction step as a function of the size of the cells for the Gaussian-like test. As expected, without prediction, the scheme has only a first-order accuracy in space, while it is a second-order scheme when the full predictor-corrector scheme with slope limiters are used.

Fig. 2. DUSTYADVECT tests with an initial condition given by the function f2 (Eq. (31)). L2 (Eq. (32)) norm as a function of the cell size for the scheme using the Minmod slope limiter (blue squares) and without predictor step (orange diamonds). The results are compared with a first-order slope (dashed line) and a second-order slope (dotted line).

4.2. Dustydiffuse

When the hydrodynamical variables (pressure and dust ratio excluded) remain constant, Eq. (8) behaves as a non-linear diffusion equation. Pure dust diffusion tests can thus be performed. In the DUSTYDIFFUSE tests (Price & Laibe 2015), the hydrodynamical step is deactivated and t_s and c_s are set constant as well. Therefore, we only solve the following equation:

$$\begin{array}{c}\hfill \frac{\partial {\rho }_{\mathrm{d}}}{\partial t}=\frac{\partial {\rho }_{\mathrm{d}}{w}_x}{\partial x}·\end{array}$$(33)

An isothermal equation of state $P_{\text{g}}=c_{\text{s}}^2(1-ϵ)\rho$ is considered. The advection velocity is then given by

$$\begin{array}{c}\hfill w_x=-t_{\mathrm{s}}{c}_{\mathrm{s}}^2\frac{\partial ϵ}{\partial x};\end{array}$$(34)

ρ being constant, the equation can be written as the a non-linear diffusion equation

$$\begin{array}{c}\hfill \frac{\partial ϵ}{\partial t}=t_{\mathrm{s}}{c}_{\mathrm{s}}^2\frac{\partial }{\partial x}\left(ϵ\frac{\partial ϵ}{\partial x}\right)·\end{array}$$(35)

Equation (35) has a self similar solution known as the Barenblatt-Pattle solution (Barenblatt 1952) that writes

$$\begin{array}{c}\hfill ϵ(t,x,C)={(t_{\mathrm{s}}{c}_{\mathrm{s}}^{2}t)}^{-1/3}(C-\frac{1}{6}\frac{x^2}{{(t_{\mathrm{s}}{c}_{\mathrm{s}}^{2}t)}^{2/3}}),\end{array}$$(36)

where C is a constant depending on the initial conditions. We consider the following initial profile

$$\begin{array}{c}\hfill ϵ(t_0,x)={ϵ}_0(1-{\left(\frac{x}{x_{\mathrm{c}}}\right)}^2),\end{array}$$(37)

which is consistent with Eq. (36) if

$$\begin{array}{cc}\hfill t_0=& \frac{C^3}{t_{\mathrm{s}}{c}_{\mathrm{s}}^2{ϵ}_0^3},\hfill \\ \hfill C=& {\left(\frac{{ϵ}_0{x}_{\mathrm{c}}}{\sqrt{6}}\right)}^{2/3}·\hfill \end{array}$$

We additionally set ρ = 1, t_s = 0.1, ϵ₀ = 0.1 and c_s = 1 and an AMR grid of ℓ_min = 4 and with ℓ_max = 10. The refinement criterion, based on the dust density gradient, forbids a variation of more than 5% between two cells.

Figure 3 shows a comparison between the outcome of the tests and the analytic solutions at t = 1, t = 5, t = 10 and t = 20. At each time the numerical results agree with the exact solution to a precision of less then 1% in L₂ norm. Even though our scheme is fundamentally designed for advection, it is also efficient at handling diffusion problems.

Fig. 3. DUSTYDIFFUSE tests. Dust ratio as a function of the position at t = 1 (red), t = 5 (green), t = 10 (blue), t = 20 (purple) compared with the exact solution (black solid lines).

4.3. Dustyshock

Another canonical test for gas and dust mixtures consists of 1D hydrodynamical shocks. For strongly coupled mixtures, these so-called DUSTYSHOCK tests are closely approximated by the same analytic solution as the usual Sod test (Sod 1978), but with the modified sound speed c̃_s

$$\begin{array}{c}\hfill {\stackrel{\sim }{c}}_{\mathrm{s}}=c_{\mathrm{s}}\sqrt{1-{ϵ}_0},\end{array}$$

ϵ₀ being the initial dust ratio. The dust ratio ϵ remains almost constant through the shock; however, pressure bumps must occur where there is either a pressure or density gradient.

The same prescription as Laibe & Price (2014a) for the stopping time is used. It writes

$$\begin{array}{c}\hfill t_{\mathrm{s}}=\frac{ϵ(1-ϵ)\rho }{K},\end{array}$$

where K is the drag coefficient defined previously. This prescription is consistent with the expression of the stopping time presented before is convenient for tests since it allows us simply parametrize the coupling between the gas and dust phase.

A DUSTYSHOCK test is performed on an AMR grid with ℓ_min = 4 and ℓ_max = 13 with a high initial dust ratio. The grid is refined with a criterion that forbids dust density variations of more than 5% between two neighbor cells. Two distinct regions, representing the left and right half of the box, are set with different initial conditions given by

$$\begin{array}{cc}& {({\rho }_0,{\mathbf{v}}_0,{P_{\mathrm{g}}}_0,{ϵ}_0)}_{\mathrm{left}}=(1,0,1,0.5)\hfill \\ \hfill & {({\rho }_0,{\mathbf{v}}_0,{P_{\mathrm{g}}}_0,{ϵ}_0)}_{\mathrm{right}}=(0.125,0,0.1,0.5).\hfill \end{array}$$

Finally, a drag coefficient K = 1000 and an adiabatic index γ = 1.4 are imposed.

Figure 4 shows the gas and dust densities, the velocity, and the gas pressure as a function of the position at t = 0.2. The Sod solution with the modified sound speed is very well recovered.

Fig. 4. DUSTYSHOCK with ϵ = 0.5. Top left: gas density as a function of position. Top right: same but for dust. Bottom left: gas pressure. Bottom right: gas and dust velocities. The AMR level (right axis) is represented with dotted blue lines. The analytic solution is given by the black solid lines. Red circles and blue crosses indicate gas and dust numerical quantities, respectively.

4.4. Dustywave

The 1D DUSTYWAVE test (Laibe & Price 2011) consists in following the evolution of an isothermal sound wave in a gas and dust mixture. Small periodic perturbations on the Eqs. (6) on the density, the dust ratio, the pressure, and the velocity are imposed

$$\begin{array}{cc}& \rho ={\rho }_0+\delta {\rho }_0,\hfill \\ \hfill & ϵ={ϵ}_0+\delta {ϵ}_0,\hfill \\ \hfill & P_{\mathrm{g}}={P_{\mathrm{g}}}_0+\delta {P_{\mathrm{g}}}_0,\hfill \\ \hfill & {\nu }_x=0+\delta {{\nu }_x}_0,\hfill \end{array}$$

where the index 0 and the symbol δ indicate respectively the initial uniform quantities and the perturbations. Laibe & Price (2011) provide the analytic solution for this test. They find that the sound waves propagate with the modified sound speed c̃_s and are damped because of the dust inertia. Large grains damp these sound waves faster than small grains as they are more massive.

Three DUSTYWAVE tests are performed using different drag coefficients K = 50 (K50, strong back-reaction regime), K = 100 (K100) and K = 1000 (K1000,weak back-reaction regime). The initial perturbations are in the form

$$\begin{array}{cc}& \delta {\rho }_0={\rho }_0\delta sin\left(2\pi \frac{x}{L}\right),\hfill \\ \hfill & \delta {{\nu }_x}_0={\nu }_0\delta sin\left(2\pi \frac{x}{L}\right),\hfill \\ \hfill & \delta {P_{\mathrm{g}}}_0=(1-{ϵ}_0)c_{\mathrm{s},0}^2\delta {\rho }_0,\hfill \\ \hfill & \delta {ϵ}_0=0,\hfill \end{array}$$

where x is the position in the box, L is the box length, c_s, 0 is the initial sound speed, and δ is a parameter that sets the amplitude of the perturbation. The initial uniform state is set such as that

$$\begin{array}{cc}& {\rho }_0=2,\hfill \\ \hfill & {ϵ}_0=0.5,\hfill \\ \hfill & {\nu }_0=1,\hfill \\ \hfill & c_{\mathrm{s},0}=1.\hfill \end{array}$$

The adiabatic index of the gas is γ = 1.000001² to simulate an isothermal soundwave propagation. The initial perturbation has a relative amplitude δ = 10⁻⁴. The simulation box has a size L = 1.0 and the grid is taken as uniform with 64 cells. The timestep Δt = 10⁻⁴ is the same for the three tests and respects the stability condition for the considered drag coefficients.

Figures 5 and 6 show the velocity ν_x and the perturbation density δρ = ρ − ρ₀ for both gas and dust at t = 4.5 for these three tests. The amplitude of the damping increases with a decreasing K and the results are increasingly less accurate as the errors due to the diffusion approximation increase, consistently with the theory (Laibe & Price 2011; Price & Laibe 2015). Larger grains, i.e, with smaller K, have more inertia and allow an efficient damping of the gas.

Fig. 5. DUSTYWAVE test. Velocity of the gas (top) and dust (bottom) phase as a function of the position at t = 4.5 for the K50, K100 and K1000 tests (from left to right) compared with the analytical solution (black lines) given by Laibe & Price (2011).

Fig. 6. DUSTYWAVE test. Same as in Fig. 5 but for the density perturbations.

We perform DUSTYWAVE tests on uniform grids with resolutions ranging from ℓ = 5 to ℓ = 9, with a small timestep Δt = 10⁻⁵ and for K = 50. We also perform runs on AMR grids with coarse resolutions ranging from ℓ_min = 4 to ℓ_min = 8. For these tests, the cells for x ∈ [0.25,0.75] have a level of refinement ℓ_max = ℓ_min + 1. We test the order of our scheme and measure the accuracy of numerical solution against a solution of reference obtained at very high resolution (ℓ = 11) in both time and space (Δt = 10⁻⁶). We do not compare the results with the analytic solution presented above as it not exact in the diffusion approximation. Figure 7 shows the L₂ errors obtained when increasing the number of cells. The scheme is first-order in space without correction and second-order when using the Minmod limiter. In the presence of AMR, the scheme keeps a second-order accuracy in space for low resolution but deteriorates at high resolutions. This is due to the estimate of the pressure gradient that is first-order at a refined interface. At high resolution, the error is smaller (almost two orders of magnitude) than the tests without a predictor step.

Fig. 7. DUSTYWAVE numerical convergence tests. L2 norm as a function of the minimum cell size for the scheme using the Minmod slope limiter (blue squares), without predictor step (orange diamonds) and with AMR (green circles). The results are compared with a first-order slope (dashed line) and a second-order slope (dotted line).

4.5. Disk settling

The SETTLING test, introduced by Price & Laibe (2015), is designed to test the algorithm in realistic astrophysical conditions. It simulates the local settling of dust grains in a disk in hydrostatic equilibrium. As in Price & Laibe (2015), we set an analytic gravitational force

$$\begin{array}{c}\hfill f_{\mathrm{grav}}=-\frac{\mathcal{G}{M}_{\star }y}{{(y^2+r^2)}^{3/2}},\end{array}$$

where 𝒢 is the gravitational constant, M_⋆ is the mass of the central star, y is the altitude and r the cylindrical radius at which the disk is simulated. The gas density at hydrostatic equilibrium for |y|≪r is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{g}}={\rho }_{\mathrm{g},0}{\mathrm{e}}^{-\frac{y^2}{2{\mathcal{H}}^2}},\end{array}$$

where ℋ is the local scale height of the disk, which is determined by the ratio ℋ/r = 0.05. We then assume an isothermal equation of state where the imposed soundspeed c_s is

$$\begin{array}{c}\hfill c_{\mathrm{s}}=\mathcal{H}{\mathrm{\Omega }}_{\mathrm{k}},\end{array}$$(38)

Ω_k being the Keplerian angular velocity at the radius r. Finally, a uniform initial dust ratio ϵ₀ is imposed.

In the terminal velocity approximation, the analytic solution for the dust velocity is directly set by the pressure gradient that compensates the gravitational force, hence

$$\begin{array}{c}\hfill {\nu }_{\mathrm{d},y}=w_Y=-\mathrm{St}{\mathrm{\Omega }}_{\mathrm{k}}\frac{y}{{(1+{(y/r)}^2)}^{3/2}},\end{array}$$(39)

which is the limit at low Stokes number (St = Ω_kt_s in this case) of the expression given by Hutchison et al. (2018). As the gas approximately remains in equilibrium, solving the settling problem consists in solving

$$\begin{array}{c}\hfill \frac{\partial {\rho }_{\mathrm{d}}}{\partial t}+\frac{\partial {\rho }_{\mathrm{d}}{w}_y}{\partial y}=0.\end{array}$$(40)

Even though the density can, in principle, be determined using the same method as the velocity, it relies on the hypothesis of an infinite dust reservoir that is inconsistent with our choice of boundaries. We choose to compare our results with a numerical solution as in Hutchison et al. (2018). We use a similar Crank-Nicholson scheme to get this solution, except that it only solves the dust density equation using the analytic dust velocity.

We perform a 2D SETTLING test at a disk orbiting around a solar mass star at radius of 50 au. The simulation box has periodic boundary conditions and L = 20 au, which is approximately 8ℋ. As in Price & Laibe (2015) and Hutchison et al. (2018), the initial mid-plane gas density ρ_g ≈ 6 × 10⁻¹³ g cm⁻³. We also set an initial dust ratio of 4.99 × 10⁻³ of millimeter grains with an intrinsic density of ρ_grain = 3 g cm⁻³. The adiabatic index of the gas is γ = 5/3.

A first test with a uniform grid at level ℓ = 8 is performed. Figure 8 shows the dust ratio as a function of y. As can be seen, the results are essentially similar to those obtained in Price & Laibe (2015) and Hutchison et al. (2018).

Fig. 8. SETTLING test on a uniform grid at level ℓ = 8 with millimeter-size dust grains. Dust ratio at t = 1 (red circles), 5 (blue circles), and 10 orbits (green circles) as a function of the altitude.

4.6. Multigrain

4.6.1. Example 1: Dustydiffuse

To benchmark the implementation of the multiple dust species in the code, ten DUSTYDIFFUSE tests are performed on uniform grids of ℓ = 7. They are operated with one dust phase and with this same phase split with 𝒩 ∈ [2,10] (MULTIGRAIN). All the dust bins have the same intrinsic properties. As all the tests are equivalent, it is expected that they give the same results. The parameters used are the same as in Sect. 4.2.

Figure 9 shows the total dust density as a function of position with 𝒩 = 5 (top) and 𝒩 = 1 (bottom) at t = 1, t = 5, t = 10, and t = 20. The results agree for the MULTIGRAIN simulation and for the single phase simulation to machine precision.

Fig. 9. Comparison between multiple phase with 𝒩 = 5 (top) and an equivalent single phase (bottom) DUSTYDIFFUSE test. Dust ratio as a function of the position at t = 1 (red), t = 5 (green), t = 10 (blue) and t = 20 s (purple) compared with the exact solution (black solid lines).

Figure 10 shows the CPU time as a function of the number of species 𝒩. We see that the MULTIGRAIN simulations are not very expensive. The CPU time agrees more with a square root scaling with 𝒩 than a linear one. As discussed by Hutchison et al. (2018) with the MULTIGRAIN algorithm in the PHANTOM code (Price et al. 2017), the monofluid formalism is a highly computationally effective tool to treat multiple phases.

Fig. 10. CPU time tCPU of the ten equivalent DUSTYDIFFUSE tests as a function of the number of species 𝒩.

4.6.2. Example 2: Dustywave

In this section, we test the cumulative back-reaction of dust on the gas and the interaction between dust species. To benchmark our MULTIGRAIN implementation in a strong back-reaction regime we run a DUSTYWAVE simulation with two dust species. The initial barycentric velocity is a sinusoidal perturbation of amplitude 10⁻⁴, the other variables are not initially perturbed. We numerically integrate the linearly perturbed equations of gas and dust mixtures assuming solutions of the type A(t)exp(ikx) as in Laibe & Price (2014c) to obtain our reference solution. The test is performed with a uniform grid of ℓ = 9 and timesteps given by the CFL condition. In this test, we consider two dust phases and a total initial dust ratio of 0.5. The first bin has a drag coefficient K = 50 and an initial dust ratio of 0.4, the second has a drag coefficient K = 1000 and an initial dust ratio of 0.1. In this case, we expect the first dust species to damp the gas efficiently.

Figure 11 shows the amplitude of the density perturbations (gas and dust) and velocity (gas) compared with the reference solution. Our results are in excellent agreement with the reference solution in terms of amplitude, period and damping rate. As expected the damping of the gas is significant. This emphasizes the fundamental role of the cumulative back-reaction on the dynamics of dust grains in presence of multiple species as the second phase K = 1000 could not damp the gas as efficiently (see Sect. 4.4).

Fig. 11. MULTIGRAIN DUSTYWAVE test. The three panels show the evolution of the maximum amplitude of the perturbation for the dust densities (top), the gas density (middle), and the gas velocity (bottom). The K = 50 and K = 1000 dust phases are shown with blue and green circles respectively. The red circles represent the gas. The semi-analytic solution is given by the dashed black line.

4.6.3. Example 3: Disk settling

The disk SETTLING test presented in Sect. 4.5 is nicely adapted to test our MULTIGRAIN solver in realistic conditions. Dust grains of various sizes are present in protoplanetary disks and they experience different dynamical evolutions. To test the solver in these conditions, a settling test with ten dust species is performed with a resolution of ℓ = 8. As before, the intrinsic grain density is 3 g cm⁻³, but every species is characterized by a proper dust ratio ϵ_j and grain size s_grain, j. These quantities are determined using a MRN-like distribution³, we use the same values for the minimum and maximum grain size as in Hutchison et al. (2018). The details of the size distribution are summarized in Table 1, which also shows the stopping time of each species in the mid-plane of the disk.

Figure 12 shows the dust ratio and velocity of every species after ten orbits compared with the one-species solution, which is a good approximation as long as the cumulative back-reaction on the gas is small. Again, the solution is very well captured by our solver.

Fig. 12. MULTIGRAIN SETTLING test. Dust ratios and dust velocities (circles) of the ten phases of the SETTLING test after ten orbits compared with the numerical (dotted black lines) and analytic (black lines) one-species solution.

Figure 13 shows the dust densities in the (xy)-plane. There is no dispersion of the values in the x-direction, and the initial symmetry of the problem is conserved to machine precision which originate from the Eulerian nature of the numerical scheme. As in previous studies, we see that ten orbits are enough to efficiently separate the dust phases. We note that an orbit at 50 au is approximately 353 years which is a few hundred times smaller than the free-fall timescale of a typical protostellar cloud of density ≈10⁻¹⁹ g cm⁻³ (Andre et al. 2000) which is approximately 10⁵−10⁶ years. An efficient settling is thus expected to happen during the collapse of this cloud, especially for large grains, e.g., s_grain >  10⁻² cm here, for which the typical settling timescale is a few orbits.

Fig. 13. MULTIGRAIN SETTLING test. Dust density for the species j = 2, 4, 6, 8, and 10 in the (xy)-plane at t = 0, 1, 5, and 10 orbits.

As in Hutchison et al. (2018), we are interested in the effects of the interaction between different dust species. Figure 14 shows a zoom of the vertical profile of ϵ₁ as a function of y. The behavior of this phase is similar to what was observed in Hutchison et al. (2018). The dashed black line shows the one-species solution. As we see, for this species, the discrepancies between the one-species solution and the MULTIGRAIN test are on the order of magnitude of the variations of ϵ₁ at the dust front, which emphasizes the fundamental character equivalent stopping time in MULTIGRAIN simulations.

Fig. 14. MULTIGRAIN SETTLING test. Magnification of the dust front as a function of y for the j = 1 dust species (brown circles) compared with the one-species solution (dotted line). These dust grains are dragged by the gas, which is itself submitted to the cumulative back-reaction of all the other dust species through the gas, hence the difference is the one-species solution.

5. Dusty protostellar collapses

The final test of the algorithm verifies again that both gas and dust are sensitive to gravity and the 3D implementation of the solver.

5.1. Dustycollapse

To study the collapse of a dense core, 3D canonical test cases, introduced by Boss & Bodenheimer (1979), but in the presence of dust (DUSTYCOLLAPSE) are performed to follow dust dynamics during the formation of the first Larson core (Larson 1969). The initial dense core, composed with gas and dust, has a mass M₀, a uniform density ρ₀ and a dust ratio ϵ₀. Its radius R₀ is set by the ratio between the thermal and gravitational energies α

$$\begin{array}{c}\hfill \alpha =\frac{5}{2}\frac{(1-{ϵ}_0)R_0}{\mathcal{G}{M}_0}\frac{k_{\mathrm{B}}{T}_{\mathrm{g}}}{\mu m_{\mathrm{H}}},\end{array}$$(41)

where k_B is the Boltzmann constant, μ is the mean molecular weight, m_H the hydrogen atomic mass, and T_g the initial gas temperature.

A barotropic equation of state is used to describe the optically thick regime at which the evolution is adiabatic. It writes

$$\begin{array}{c}\hfill \frac{P_{\mathrm{g}}}{(1-ϵ)\rho }=c_{\mathrm{s}}^2[1+{\left((1-ϵ)\frac{\rho }{{\rho }_{\mathrm{ad}}}\right)}^{\gamma -1}],\end{array}$$(42)

the evolution becoming adiabatic above the critical density ρ_ad = 10⁻¹³ g cm⁻³ (Larson 1969). For the sake of simplicity, this critical density is taken constant in this paper. A dependency with the dust ratio will be discussed in future works.

As in Tricco et al. (2017), we consider an unique grain density ρ_grain = 3 g cm⁻³ which is a good approximation for a combination of carbonaceous and silicate grains.

5.2. Setup

All the models are set with M₀ = 1 M_⊙, μ = 2.31, γ = 5/3, and T_g = 10 K. AMR grids of ℓ_min = 5 and ℓ_max = 17 are considered. With respect to the criterion given by Truelove et al. (1997) to avoid artificial clumps, the grid is refined to impose at least 15 points per Jeans length. The region surrounding the core has a density that is one hundred time lower than the initial core but the same temperature.

The initial density and pressure gradients between the core and the ambient medium are numerical artifacts. To avoid unphysical variations of the dust ratio in these regions and unnecessary constraints on the timestep, the differential advection velocity w_k is set to zero outside the initial core.

We consider the first hydrostatic core as the region where the gas density is higher than 10^−12.5 g cm⁻³.

5.3. Validity of the diffusion approximation

The diffusion approximation is valid as long as the stopping time t_s is shorter than the dynamical timescale which is, for gravitational collapses, roughly the free-fall time

$$\begin{array}{c}\hfill t_{\mathrm{ff}}\equiv \sqrt{\frac{3\pi }{32\mathcal{G}\rho }}·\end{array}$$

With an initial dust ratio of 0.01 and a temperature of 10 K, the Stokes number is

$$\begin{array}{c}\hfill \mathrm{St}\sim 0.114\left(\frac{M_0}{1M_{\odot }}\right)\left(\frac{{\rho }_{\mathrm{grain}}}{3\mathrm{g}{\mathrm{cm}}^{-3}}\right)\left(\frac{s_{\mathrm{grain}}}{0.05\mathrm{cm}}\right){\left(\frac{\alpha }{0.5}\right)}^{3/2}<1.\end{array}$$(43)

For a wide range of density and grain size the stopping time is short compared to the dynamical timescale. Initially, the stopping time is shorter than the free-fall timescale for grain sizes up to ≃4 mm. However, as the Stokes number may vary, the validity of the approximation needs to be discussed during the collapse (see Sect. 5.5).

5.4. Free-fall timescale for strongly coupled mixtures

A core only composed with gas would collapse at the following free-fall timescale

$$\begin{array}{c}\hfill t_{\mathrm{ff},\mathrm{g}}\equiv \sqrt{\frac{3\pi }{32\mathcal{G}{\rho }_{\mathrm{g}}}}·\end{array}$$

For a perfectly coupled gas and dust mixture (t_s ≪ t_ff) it writes

$$\begin{array}{c}\hfill t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32\mathcal{G}\rho }}·\end{array}$$

In the context of the Boss and Bodenheimer test, this timescale indirectly depends on the initial dust ratio

$$\begin{array}{c}\hfill t_{\mathrm{ff}}=\frac{\pi }{5\sqrt{5}}\mathcal{G}{M}_0{\left(\frac{\alpha \mu m_{\mathrm{H}}}{k_{\mathrm{b}}{T}_{\mathrm{g}}(1-{ϵ}_0)}\right)}^{3/2}·\end{array}$$(44)

Physically, this is due to the fact that two cores with the same initial α but different ϵ₀ have different initial radius, hence free-fall timescale. Indeed, R₀ increases with ϵ₀ so that the gas provides the same thermal support against gravity.

To test this relation, non-rotating DUSTYCOLLAPSE with various dust ratios and the same α = 0.5 are performed. The condition t_s ≪ t_ff is ensured by considering very small grains with s = 10⁻⁷ cm. Figure 15 shows the ratio of the free-fall timescale of a dusty cloud to that of a pure-gas cloud compared with the analytical solution. We define the free-fall timescale as the time at which the peak density reaches ρ_ad. The scaling relation is very well verified. As expected, the gas, and the dust, are sensitive to gravity and fall at the mixture free-fall timescale.

Fig. 15. Ratio of the free-fall timescale of the mixture to the free-fall timescale of the gas. The red circle corresponds to the fiducial simulation with gas only, and the blue circles correspond to gas and dust mixtures for various dust ratios. The black solid line is the analytical solution.

5.5. Core properties

Non-rotating DUSTYCOLLAPSE of α  =  0.5 are performed with an initial dust-to-gas ratio of 1% (corresponding to ϵ  ∼  0.00990099). Single grain species of size 1 μm (Col1), 10 μm (Col10) and 100 μm (Col100) are considered. The results are compared, when the gas maximum density reaches 10⁻¹¹ g cm⁻³, with a fiducial collapse without any treatment of dust (ColGAS).

Figure 16 shows the total mass of the first hydrostatic core M_core, the dust mass ratio M_d/M_core and its formation time t_core as a function of the grain size (bottom x-axis) and the initial Stokes number (top x-axis). The first core properties depend on the amount of dust in the initial dense core. In the fiducial ColGAS test, the gas density is assimilated to the total density, as is often the case in the literature. This results in a lower first-core mass and a shorter formation time. We note that the barotropic equation of state implicitly assumes a constant dust-to-gas ratio of 1%, which is slightly inconsistent. Dust-to-gas ratios might be even higher in core forming regions (Hopkins & Lee 2016; Tricco et al. 2017) which could lead to even more important discrepancies. In terms of dynamics, small grains (s <  100 μm) do not have a significant impact as they fall alongside the gas. However, the differential velocity between the gas and the dust is proportional to the stopping time, hence large grains (s ≥ 100 μm) tend to fall substantially faster than the gas. Their infall provokes a slight acceleration of the first core formation as it changes the balance between the gravity and the thermal support of the gas. The mass of the core decreases as the grain size increases because the gas has less time to be accreted during the collapse. The settling of dust in the central regions of the core, however, leads to an increase in M_d/M_core.

Fig. 16. Properties of the first Larson core when ρg ∼ 10−11 g cm−3 for Col1, Col10, and Col100. The core mass, the dust mass ratio, and the formation time are shown as a function of the grain size (blue circles). The results are compared with the fiducial ColGAS simulation (dashed red lines). The top x-axis shows the initial Stokes number in the core St0.

Figures 17 and 18 show the radial profiles for the Col10 and Col100 tests, respectively, of the gas density, the dust ratio, the gas and dust velocities, and the Stokes number. As can be seen, the strong dust depletion in the outer regions of the collapse provides the enrichment of the core. However, as the damping is more efficient in the high density regions (see the velocities v_g ≈ v_d), the dust ratio increases at a slower rate in the inner regions of the collapse, i.e, the first core. If the maximum increase of dust ratio in the outer regions of the collapse is about ∼20% for the Col10 test, it goes up to ∼300% for the Col100 test. We see in the Col100 case that the terminal velocity approximation probably leads to non-negligible errors as the Stokes number is close to one in the outer regions of the collapse . More accurate simulations using a future implementation of the full monofluid formalism (Laibe & Price 2014a) should investigate these errors.

Fig. 17. Radial profiles of the Col10 test. Top left: gas density. Top right: dust ratio. Bottom left: gas and dust velocities. Bottom right: stokes number. The horizontal green line corresponds to a dust-to-gas ratio of 1%.

Fig. 18. Same as Fig.17, but for the Col100 test.

If initial solid-body rotation is considered, it is expected that the large dust grains will settle in the mid-plane of the protoplanetary disk (Bate & Lorén-Aguilar 2017). Future works will take this into account and the presence of multiple grain species as well.

6. Conclusions and perspectives

We have presented an Eulerian approach to treat the dynamics of small dust grains. It is efficient in the diffusion regime and can be employed to treat multiple dust species simultaneously and efficiently. After a summary of the theory, we presented the numerical scheme that we implemented in the AMR code RAMSES. It successfully passed the canonical validation tests for advection schemes and dust dynamics solvers, i.e, DUSTYADVECT, DUSTYDIFFUSE, DUSTYSHOCK, DUSTYWAVE, SETTLING, and DUSTYCOLLAPSE. We show that our scheme has a second-order accuracy in space on uniform grids and intermediate between first- and second-order on AMR grids. The method also appears to efficiently treat a non-linear diffusion problem. The DUSTYSHOCK, DUSTYWAVE, SETTLING and DUSTYCOLLAPSE, show that the waves and shock propagate at the correct velocity, and that the dust phase feels the common forces between gas and dust, e.g., gravity. The DUSTYDIFFUSE test operated with a split dust phase shows that our method is able to efficiently treat multiple dust species simultaneously as the computation time scales in $\sqrt{N}$. We emphasize the importance of the cumulative back-reaction on both the gas and the dust with the MULTIGRAIN DUSTYWAVE tests. We show with a MULTIGRAIN SETTLING test that our solver efficiently treats multiple phases in a realistic astrophyiscal environment. We also observe the same interactions between the dust phases that were uncovered in Hutchison et al. (2018).

The method is applied to study dust dynamics during the formation of the first hydrostatic core. As found in Bate & Lorén-Aguilar (2017), grains larger than s ≃ 100 μm can fall substantially faster than the gas. Above this size, the grains settle in the core, enriching its dust content and speeding up the collapse. We showed that the core properties are not affected very much by the dynamics of dust for small grains (s <  100 μm). However, they depend on the presence of dust itself as the thermal-to-gravitational energy ratio α depends on both the gas and the dust density.

The scheme was presented, for simplicity, in the hydrodynamical case; however, it is implemented in the radiation non-ideal MHD version of RAMSES. Realistic simulations of star formation would require us to consider the impact of dust on magnetic fields and radiative transfer coupling with the mixture. In that perspective, later work will focus on taking into account the impact of the variation of the dust ratio on the opacity and the magnetic resistivities.

In the paper II we will discuss in more detail the grain dynamics during the collapse of rotating clouds with MULTIGRAIN simulations.

Acknowledgments

First, we sincerely thank the referee for providing extremely useful remarks and ideas which helped to significantly improve the quality of our work. We acknowledge financial support from the Programme National de Physique Stellaire (PNPS) of CNRS/INSU, CEA, and CNES, France. This work was granted access to the HPC resources of CINES (Occigen) under the allocation DARI A0020407247 made by GENCI. Computations were also performed at the Common Computing Facility (CCF) of the LABEX Lyon Institute of Origins (ANR-10-LABX-0066). This work took part under the programs ISM3D and Core2disk of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. This project was partly supported by the IDEXLyon project (contract nANR-16-IDEX-0005) under the auspices University of Lyon. Some of the DUSTYCOLLAPSE plots were generated using the very efficient Osiris library developed by Neil Vaytet, Tommaso Grassi and Matthias González whom we thank. In particular, we thank Matthias González again for his useful advice on the paper. We also thank Daniel Price for his useful advice and Mark Hutchison for providing his solver to get the solution of the settling test. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 823823.

References

All Tables

All Figures

Fig. 1. DUSTYADVECT tests using the function f1 (Eq. (30)). Dust density after one on period various grids (ℓ = 7 in blue, ℓ = 9 in orange, ℓ = 11 in green, ℓ = 13 in red) as a function of the position compared with the analytic solution (black solid line). We present this test using four different slope limiters. Top left: No predictor step. Top right: Minmod slope limiter. Bottom left: Van-Leer slope limiter. Bottom right: superbee slope limiter.

In the text

	Fig. 2. DUSTYADVECT tests with an initial condition given by the function f2 (Eq. (31)). L2 (Eq. (32)) norm as a function of the cell size for the scheme using the Minmod slope limiter (blue squares) and without predictor step (orange diamonds). The results are compared with a first-order slope (dashed line) and a second-order slope (dotted line).
In the text

	Fig. 3. DUSTYDIFFUSE tests. Dust ratio as a function of the position at t = 1 (red), t = 5 (green), t = 10 (blue), t = 20 (purple) compared with the exact solution (black solid lines).
In the text

	Fig. 4. DUSTYSHOCK with ϵ = 0.5. Top left: gas density as a function of position. Top right: same but for dust. Bottom left: gas pressure. Bottom right: gas and dust velocities. The AMR level (right axis) is represented with dotted blue lines. The analytic solution is given by the black solid lines. Red circles and blue crosses indicate gas and dust numerical quantities, respectively.
In the text

	Fig. 5. DUSTYWAVE test. Velocity of the gas (top) and dust (bottom) phase as a function of the position at t = 4.5 for the K50, K100 and K1000 tests (from left to right) compared with the analytical solution (black lines) given by Laibe & Price (2011).
In the text

	Fig. 6. DUSTYWAVE test. Same as in Fig. 5 but for the density perturbations.
In the text

	Fig. 7. DUSTYWAVE numerical convergence tests. L2 norm as a function of the minimum cell size for the scheme using the Minmod slope limiter (blue squares), without predictor step (orange diamonds) and with AMR (green circles). The results are compared with a first-order slope (dashed line) and a second-order slope (dotted line).
In the text

	Fig. 8. SETTLING test on a uniform grid at level ℓ = 8 with millimeter-size dust grains. Dust ratio at t = 1 (red circles), 5 (blue circles), and 10 orbits (green circles) as a function of the altitude.
In the text

	Fig. 9. Comparison between multiple phase with 𝒩 = 5 (top) and an equivalent single phase (bottom) DUSTYDIFFUSE test. Dust ratio as a function of the position at t = 1 (red), t = 5 (green), t = 10 (blue) and t = 20 s (purple) compared with the exact solution (black solid lines).
In the text

	Fig. 10. CPU time tCPU of the ten equivalent DUSTYDIFFUSE tests as a function of the number of species 𝒩.
In the text

	Fig. 11. MULTIGRAIN DUSTYWAVE test. The three panels show the evolution of the maximum amplitude of the perturbation for the dust densities (top), the gas density (middle), and the gas velocity (bottom). The K = 50 and K = 1000 dust phases are shown with blue and green circles respectively. The red circles represent the gas. The semi-analytic solution is given by the dashed black line.
In the text

	Fig. 12. MULTIGRAIN SETTLING test. Dust ratios and dust velocities (circles) of the ten phases of the SETTLING test after ten orbits compared with the numerical (dotted black lines) and analytic (black lines) one-species solution.
In the text

	Fig. 13. MULTIGRAIN SETTLING test. Dust density for the species j = 2, 4, 6, 8, and 10 in the (xy)-plane at t = 0, 1, 5, and 10 orbits.
In the text

	Fig. 14. MULTIGRAIN SETTLING test. Magnification of the dust front as a function of y for the j = 1 dust species (brown circles) compared with the one-species solution (dotted line). These dust grains are dragged by the gas, which is itself submitted to the cumulative back-reaction of all the other dust species through the gas, hence the difference is the one-species solution.
In the text

	Fig. 15. Ratio of the free-fall timescale of the mixture to the free-fall timescale of the gas. The red circle corresponds to the fiducial simulation with gas only, and the blue circles correspond to gas and dust mixtures for various dust ratios. The black solid line is the analytical solution.
In the text

	Fig. 16. Properties of the first Larson core when ρg ∼ 10−11 g cm−3 for Col1, Col10, and Col100. The core mass, the dust mass ratio, and the formation time are shown as a function of the grain size (blue circles). The results are compared with the fiducial ColGAS simulation (dashed red lines). The top x-axis shows the initial Stokes number in the core St0.
In the text

	Fig. 17. Radial profiles of the Col10 test. Top left: gas density. Top right: dust ratio. Bottom left: gas and dust velocities. Bottom right: stokes number. The horizontal green line corresponds to a dust-to-gas ratio of 1%.
In the text

	Fig. 18. Same as Fig.17, but for the Col100 test.
In the text