Small dust grain dynamics on adaptive mesh-refinement grids. I. Methods

TL;DR

This paper introduces a novel Eulerian numerical method for simulating small dust grain dynamics in star and disk formation, implemented in the RAMSES code, and demonstrates its accuracy and application to protostellar collapse scenarios.

Contribution

It presents a new finite volume, second-order Godunov method for dust-gas mixtures with multiple dust species in adaptive mesh refinement codes.

Findings

The scheme is second-order accurate on uniform grids.

Dust grains larger than ~100 micrometers decouple from gas during collapse.

The method successfully passes various benchmark tests.

Abstract

Small dust grains are essential ingredients of star, disk and planet formation. We present an Eulerian numerical approach to study small dust grains dynamics in the context of star and protoplanetary disk formation. It is designed for finite volume codes. We use it to investigate dust dynamics during the protostellar collapse. We present a method to solve the monofluid equations of gas and dust mixtures with several dust species in the diffusion approximation implemented in the adaptive-mesh-refinement code RAMSES. It uses a finite volume second-order Godunov method with a predictor-corrector MUSCL scheme to estimate the fluxes between the grid cells. We benchmark our method against six distinct tests, dustyadvect, dustydiffuse, dustyshock, dustywave, settling and dustycollapse. We show that the scheme is second-order accurate in space on uniform grids and intermediate between second and first-order on non-uniform grids. We apply our method on various dustycollapse simulations of 1 solar mass cores composed of gas and dust. We developed an efficient approach to treat gas and dust dynamics in the diffusion regime on grid based codes. The canonical tests have been successfully passed. In the context of protostellar collapse, we show that dust is less coupled to the gas in the outer regions of the collapse where grains larger than approximately 100 micrometers fall significantly faster than the gas.

Tables

Table 1. Table 1: Dust distribution and stopping time at z = 0 𝑧 0 z=0 for the multigrain settling test

$j$	$s_{\mathrm{grain},j}$ (cm)	${ϵ}_j$	$t_s$ (s) at $z=0$
1	$1.0\times {10}^{-5}$	$3.99\times {10}^{-5}$	$40$
2	$2.78\times {10}^{-5}$	$6.65\times {10}^{-5}$	$111$
3	$7.74\times {10}^{-5}$	$1.11\times {10}^{-4}$	$310$
4	$2.15\times {10}^{-4}$	$1.85\times {10}^{-4}$	$860$
5	$5.99\times {10}^{-4}$	$3.09\times {10}^{-4}$	$2396$
6	$1.67\times {10}^{-3}$	$5.15\times {10}^{-4}$	$6680$
7	$4.64\times {10}^{-3}$	$8.59\times {10}^{-4}$	$18560$
8	$1.29\times {10}^{-2}$	$1.43\times {10}^{-3}$	$51600$
9	$3.59\times {10}^{-2}$	$2.39\times {10}^{-3}$	$143600$
10	$1.0\times {10}^{-1}$	$3.99\times {10}^{-3}$	$400000$

Equations

$$\mathbf{F}=-\frac{m_{\mathrm{grain}}}{t_{\mathrm{s},\mathrm{grain}}}\mathbf{\Delta v},$$

$$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}}},$$

$$\frac{{\mathrm{d}}_{\mathrm{g}}{\rho }_{\mathrm{g}}}{\mathrm{d}t}$$

$$\frac{{\mathrm{d}}_{\mathrm{d}}{\rho }_{\mathrm{d}}}{\mathrm{d}t}$$

$${\rho }_{\mathrm{g}}\frac{{\mathrm{d}}_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}}{\mathrm{d}t}$$

$${\rho }_{\mathrm{d}}\frac{{\mathrm{d}}_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}}{\mathrm{d}t}$$

$$\frac{{\mathrm{d}}_{\mathrm{g}}{e}_{\mathrm{g}}}{\mathrm{d}t}$$

$$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}},$$

$$K\equiv \frac{{\rho }_{\mathrm{d}}{\rho }_{\mathrm{g}}}{\rho t_{\mathrm{s}}}.$$

$$\mathbf{v}\equiv \frac{{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}+{\rho }_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}}}.$$

$$ϵ\equiv \frac{{\rho }_{\mathrm{d}}}{\rho },$$

$${\mathbf{v}}_{\mathrm{g}}=\mathbf{v}-ϵ\mathbf{\Delta v},$$

$${\mathbf{v}}_{\mathrm{d}}=\mathbf{v}+(1-ϵ)\mathbf{\Delta v},$$

$${\rho }_{\mathrm{g}}=\left(1-ϵ\right)\rho ,$$

$${\rho }_{\mathrm{d}}=ϵ\rho .$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}ϵ}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{\Delta v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}ϵ}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{\Delta v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}$$

$$\mathbf{\Delta v}=t_{\mathrm{s}}\frac{\nabla P_{\mathrm{g}}}{(1-ϵ)\rho }.$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}ϵ}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}$$

$$\frac{\partial \rho ϵ}{\partial t}+\nabla \cdot \left[\rho ϵ\left(\mathbf{v}+\frac{t_{\mathrm{s}}\nabla P_{\mathrm{g}}}{\rho }\right)\right]=0.$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{ϵ}_k}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{e}_{\mathrm{g}}}{\mathrm{d}t}$$

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

$$\mathcal{E}\equiv \sum _{l=1}^{\mathcal{N}}{ϵ}_l,$$

$${\mathcal{T}}_{\mathrm{s}}\equiv \frac{1}{\mathcal{E}}\sum _{l=1}^{\mathcal{N}}{ϵ}_l{T}_{\mathrm{s},l}.$$

$${\rho }_{\mathrm{g}}$$

Full text

*Aims. *We present an Eulerian numerical approach to study small dust grain dynamics in the context of star and protoplanetary disk formation. It is designed for finite volume codes. We use it to investigate dust dynamics during the protostellar collapse.

*Methods. * We present a method to solve the monofluid equations of gas and dust mixtures with several dust species in the diffusion approximation implemented in the adaptive-mesh-refinement code RAMSES. It uses a finite volume second-order Godunov method with a predictor-corrector MUSCL scheme to estimate the fluxes between the grid cells.

*Results. *We benchmark our method against six distinct tests, dustyadvect, dustydiffuse, dustyshock, dustywave, settling and dustycollapse. We show that the scheme is second-order accurate in space on uniform grids and intermediate between second- and first-order on non-uniform grids. We apply our method on various dustycollapse simulations of $1M_{\odot}$ cores composed of gas and dust.

*Conclusions. *We developed an efficient approach to treat gas and dust dynamics in the diffusion regime on grid-based codes. The canonical tests were successfully passed. In the context of protostellar collapse, we show that dust is less coupled to the gas in the outer regions of the collapse where grains larger than $\simeq 100\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{m}$ fall significantly faster than the gas.

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 $\sim 100\leavevmode\nobreak\ \mathrm{n}\mathrm{m}$ 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 $\rho_{\mathrm{g}}$ . Let us consider a grain of radius $s_{\mathrm{grain}}$, mass $m_{\mathrm{grain}}$, and intrinsic density $\rho_{\mathrm{grain}}$. It exchanges momentum with the surrounding gas molecules through microscopic collisions. Macroscopically, this leads to a drag force on the grain that writes

[TABLE]

where the stopping time of the grain $t_{\mathrm{s,grain}}$, is the typical time required by the dust grain to adjust its dynamics to a change in the gas velocity, and $\mathbf{\Delta 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)

[TABLE]

where $c_{\mathrm{s}}$ and $\gamma$ 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_{\mathrm{s}}$ and $\rho_{\mathrm{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 $\rho_{\mathrm{d}}$ and an advection velocity $\mathbf{v_{\mathrm{d}}}$. The dust density $\rho_{\mathrm{d}}$ must not be confused with the intrinsic grain density $\rho_{\mathrm{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 $\rho_{\mathrm{g}}$ and an advection velocity $\mathbf{v_{\mathrm{g}}}$.

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

[TABLE]

where $\mathbf{\Delta v}\equiv\mathbf{v_{\mathrm{d}}}-\mathbf{v_{\mathrm{g}}}$ is the differential velocity between the two fluids; $\mathbf{f}$ is the specific force acting on both gas and dust; and $\mathbf{f_{\mathrm{g}}}$ and $\mathbf{f_{\mathrm{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 $\mathbf{f_{\mathrm{g}}}$ and $\mathbf{f_{\mathrm{d}}}$ are zero, and $\mathbf{f}$ is either zero or the gravity force. The gas pressure $P_{\mathrm{g}}$ is given by $P_{\mathrm{g}}=(\gamma-1)\rho_{\mathrm{g}}e_{\mathrm{g}}$, $e_{\mathrm{g}}$ being the gas specific internal energy. The gas derivative $\frac{\mathrm{d}_{\mathrm{g}}}{\mathrm{d}t}\equiv\frac{\partial}{\partial t}+\left(\mathbf{v}_{\mathrm{g}}\cdot\nabla\right)$ differs from the dust derivative $\frac{\mathrm{d}_{\mathrm{d}}}{\mathrm{d}t}\equiv\frac{\partial}{\partial t}+\left(\mathbf{v}_{\mathrm{d}}\cdot\nabla\right)$. We also define the stopping time as

[TABLE]

and the drag coefficient $K$ as

[TABLE]

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 $\rho\equiv\rho_{\mathrm{d}}+\rho_{\mathrm{g}}$ (Laibe & Price, 2014a, b). It is moving at the mixture barycentric velocity $\mathbf{v}$:

[TABLE]

We also define the dust ratio

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

Here only one Lagrangian derivative $\frac{\mathrm{d}}{\mathrm{d}t}\equiv\frac{\partial}{\partial t}+\left(\mathbf{v}\cdot\nabla\right)$ 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. (2.2.1). 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_{\mathrm{dyn}}$ of the system, two simplifications can be made.

In the so-called diffusion approximation, $\left|\left|\mathbf{\Delta v}^{2}/\mathbf{v}^{2}\right|\right|$, $\left|\left|\mathbf{\Delta v}\otimes\mathbf{\Delta v}/\mathbf{v}^{2}\right|\right|$, $\left|\left|\mathbf{\Delta v}\times(\nabla\times(1-\epsilon)\mathbf{\Delta v})/\mathbf{v}^{2}\right|\right|$ and $\left|\left|\mathbf{\Delta v}\times(\nabla\times\epsilon\mathbf{\Delta v})/\mathbf{v}^{2}\right|\right|$ are second-order in $\mathrm{St}\equiv\left(t_{\mathrm{s}}/t_{\mathrm{dyn}}\right)$, where $\mathrm{St}$ is the Stokes number. 111 $\left|\left|.\right|\right|$ indicates the L2 norm of either a tensor or a vector. System (2.2.2) reduces to

[TABLE]

The errors caused by the diffusion approximation are marginal as long as $\mathrm{St}\ll 1$. We note that the first-order term $-\frac{P_{\mathrm{g}}}{\rho(1-\epsilon)}\nabla\cdot\left(\epsilon\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 $\mathrm{St}\ll 1$. In this case $\left|\left|\frac{\mathrm{d}\mathbf{\Delta v}}{\mathrm{d}t}/\left(\frac{\mathbf{\Delta v}}{t_{\mathrm{s}}}\right)\right|\right|$ and $\left|\left|(\mathbf{\Delta v}\cdot\mathbf{\nabla})\mathbf{v}/\left(\frac{\mathbf{\Delta v}}{t_{\mathrm{s}}}\right)\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

[TABLE]

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

[TABLE]

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

[TABLE]

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 ($\mathcal{N}+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 $\left(\mathcal{N}+1\right)$ phase mixtures with $\mathcal{N}$ 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

[TABLE]

where $\epsilon_{k}$ is the dust ratio of the phase $k$ and $T_{\mathrm{s},k}$ is the effective stopping time of the dust phase $k$ defined as

[TABLE]

where $t_{\mathrm{s},k}$ is the individual stopping time of the phase $k$. The $\sum_{l=1}^{\mathcal{N}}\frac{\epsilon_{l}}{1-\epsilon_{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

[TABLE]

and the mean stopping time

[TABLE]

The gas and dust densities are simply given by

[TABLE]

When $\mathcal{N}=1$, Eqs. (2.2.4) 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

[TABLE]

where $\rho_{\mathrm{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 $\mathbb{I}$ 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; Rosdahl et al., 2013; Commerçon et al., 2014; 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

[TABLE]

the state $\mathbb{U}$ vector and flux $\mathbb{F}$ vector being

[TABLE]

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 $\mathbb{U}$ over the cells. For a 1D problem integrated over time (for $t\in[t_{1},t_{2}]$) and space (for $x\in[x_{1},x_{2}]$), the previous system writes

[TABLE]

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

[TABLE]

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 $\Delta x$ and a timestep $\Delta t$, the scheme writes

[TABLE]

where the discretized fluxes $\mathbb{F}^{n+1/2}_{i\pm 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 $\Delta x=\Delta y=\Delta 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 $\ell$ with a value between $\ell_{\mathrm{min}}$ for the coarsest cells and $\ell_{\mathrm{max}}$ for the finest. The length of a cell of level $\ell$ is

[TABLE]

where $L_{\mathrm{box}}$ is the size of the simulation box. For a uniform grid with a unique level $\ell$, the effective resolution is simply $2^{\ell}$.

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

[TABLE]

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

• •

The fine-to-coarse interfaces, when the neighbor cell is at level $\ell-1$;

• •

the fine-fine interfaces, when the neighbor cell is also at level $\ell$;

• •

the coarse-to-fine interfaces, when the neighbor cell is at level $\ell+1$.

In RAMSES, the levels are updated from fine to coarse. During the update of level $\ell+1$, the fluxes at the interface with level $\ell$ are used to update the level $\ell+1$ and stored for the later update of the level $\ell$. 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 $\ell$. The coarse neighbor cells at level $\ell-1$ are virtually refined to compute the flux at the boundaries with level $\ell$.

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 $\mathcal{N}$ 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

[TABLE]

The new state vector $\mathbb{U}$ writes

[TABLE]

the flux $\mathbb{F}_{\mathrm{H}}$ is similar to the flux of pure hydrodynamics and is given by

[TABLE]

the flux $\mathbb{F}_{\mathrm{\Delta}}$ accounts for dust diffusion and is given by

[TABLE]

where $\mathbf{w_{k}}$, the differential advection velocity between the dust species k and the barycenter, is

[TABLE]

and $\mathbf{w_{\mathrm{g,b}}}$ is the differential velocity between the gas and the barycenter that writes

[TABLE]

or

[TABLE]

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 $\mathbb{F}_{\mathrm{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 $\rho$ is different from the gas density. As a consequence, the pressure and the wave speeds must be computed using $\rho_{\mathrm{g}}=\rho(1-\mathcal{E})$ instead of $\rho$.

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 $\mathbb{F}_{\mathrm{\Delta}}$. 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

[TABLE]

where $\tilde{\rho}_{\mathrm{d},k}$ refers to ${\rho}_{\mathrm{d},k}$ after its advection as a passive scalar at the velocity $\mathbf{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 $\mathcal{F}(\rho_{\mathrm{d}})\equiv\rho_{\mathrm{d}}\mathbf{w}$. 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 $\rho$ is constant.

Differential advection velocity:

During the diffusion step, we aim to compute the differential advection velocity ${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 $\Delta x$. For example, its component in the $x$-direction, is

[TABLE]

where $\Delta x_{i-1,j,k}$ and $\Delta 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}^{n}_{i,j,k}$ is given by

[TABLE]

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

[TABLE]

where $\Delta_{\sigma}{\rho_{\mathrm{d}}}_{i,j,k}^{n}$ and $\Delta_{\sigma}{w_{\sigma}}^{n}_{i,j,k}$ are the TVD variations of $\rho_{\mathrm{d}}$ and $w_{\sigma}$ in the direction $\sigma$. After the prediction of ${\rho_{\mathrm{d}}}_{i,j,k}^{n+1/2}$, we interpolate $\rho_{\mathrm{d}}$ and $w_{\sigma}$ 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

[TABLE]

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

[TABLE]

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

[TABLE]

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}^{n+1/2}_{i,j-1/2,k}$ and $\mathcal{F}^{n+1/2}_{i,j,k-1/2}$ are estimated in a similar way. The dust density is finally updated according to

[TABLE]

Energy :

Similarly and simultaneously with the dust diffusion step, the energy is updated after the hydrodynamical step to account for the energy component of $\mathbb{F}_{\mathrm{\Delta}}$. 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 $\Delta t$ verifies the split CFL condition

[TABLE]

where $C_{\mathrm{CFL}}<1$ is a safety factor. In addition, the hydrodynamical step imposes another stability condition that writes

[TABLE]

where $v_{\sigma}$ is the mixture velocity in the direction $\sigma$. Dust diffusion is stable without intervention on the timestep as long as

[TABLE]

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

[TABLE]

Considering an initial condition $\rho_{\mathrm{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 $\rho_{\mathrm{d}}(x,t)=f(x-w_{x}t)$. In the remainder of the section, $w_{x}=1$.

At first, four dustyadvect tests are performed. We impose the initial function $f_{1}$, which writes $\forall x\in\left[0,L\right]$

[TABLE]

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 $\ell=7$ (128 cells) to $\ell=13$ (8192 cells) are considered. Extremely small timesteps compared with the CFL condition $\Delta t=8\times 10^{-6}$ 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.

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_{2}$ is imposed to quantify the truncation error in space of the scheme which writes, $\forall x\in\left[0,L\right]$,

[TABLE]

We use a very small timestep $\Delta t=10^{-8}$ 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_{2}$ norm

[TABLE]

Figure 2 shows the evolution of the $L_{2}$ 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.

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_{\mathrm{s}}$ and $c_{\mathrm{s}}$ are set constant as well. Therefore, we only solve the following equation:

[TABLE]

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

[TABLE]

$\rho$ being constant, the equation can be written as the a non-linear diffusion equation

[TABLE]

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

[TABLE]

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

[TABLE]

which is consistent with Eq. (38) if

[TABLE]

We additionally set $\rho=1$, $t_{\mathrm{s}}=0.1$, $\epsilon_{0}=0.1$ and $c_{\mathrm{s}}=1$ and an AMR grid of $\ell_{\mathrm{min}}=4$ and with $\ell_{\mathrm{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_{2}$ norm. Even though our scheme is fundamentally designed for advection, it is also efficient at handling diffusion problems.

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 $\tilde{c}_{\mathrm{s}}$

[TABLE]

$\epsilon_{0}$ being the initial dust ratio. The dust ratio $\epsilon$ 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

[TABLE]

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 $\ell_{\mathrm{min}}=4$ and $\ell_{\mathrm{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

[TABLE]

Finally, a drag coefficient $K=1000$ and an adiabatic index $\gamma=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.

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. (2.2.2) on the density, the dust ratio, the pressure, and the velocity are imposed

[TABLE]

where the index [math] and the symbol $\delta$ 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 $\tilde{c}_{\mathrm{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

[TABLE]

where $x$ is the position in the box, $L$ is the box length, $c_{\mathrm{s},0}$ is the initial sound speed, and $\delta$ is a parameter that sets the amplitude of the perturbation. The initial uniform state is set such as that

[TABLE]

The adiabatic index of the gas is $\gamma=1.000001$222 Taking exactly $\gamma=1$ is not possible with RAMSES as the internal energy is computed as $\frac{P_{\mathrm{g}}}{\gamma-1}$ to simulate an isothermal soundwave propagation. The initial perturbation has a relative amplitude $\delta=10^{-4}$. The simulation box has a size $L=1.0$ and the grid is taken as uniform with $64$ cells. The timestep $\Delta t=10^{-4}$ is the same for the three tests and respects the stability condition for the considered drag coefficients.

Figures 5 and 6 show the velocity $v_{x}$ and the perturbation density $\delta\rho=\rho-\rho_{0}$ 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.

We perform dustywave tests on uniform grids with resolutions ranging from $\ell=5$ to $\ell=9$, with a small timestep $\Delta t=10^{-5}$ and for $K=50$. We also perform runs on AMR grids with coarse resolutions ranging from $\ell_{\mathrm{min}}=4$ to $\ell_{\mathrm{min}}=8$. For these tests, the cells for $x\in\left[0.25,0.75\right]$ have a level of refinement $\ell_{\mathrm{max}}=\ell_{\mathrm{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 ($\ell=11$) in both time and space ($\Delta t=10^{-6}$). 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_{2}$ 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.

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

[TABLE]

where $\mathcal{G}$ is the gravitational constant, $M_{\star}$ 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|\ll r$ is

[TABLE]

where $\mathcal{H}$ is the local scale height of the disk, which is determined by the ratio $\mathcal{H}/r=0.05$. We then assume an isothermal equation of state where the imposed soundspeed $c_{\mathrm{s}}$ is

[TABLE]

$\Omega_{\mathrm{k}}$ being the Keplerian angular velocity at the radius $r$. Finally, a uniform initial dust ratio $\epsilon_{0}$ 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

[TABLE]

which is the limit at low Stokes number ($\mathrm{St}=\Omega_{\mathrm{k}}t_{\mathrm{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

[TABLE]

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\mathcal{H}$. As in Price & Laibe (2015) and Hutchison et al. (2018), the initial mid-plane gas density $\rho_{\mathrm{g}}\approx 6\times 10^{-13}$g cm*-3*. We also set an initial dust ratio of $4.99\times 10^{-3}$ of millimeter grains with an intrinsic density of $\rho_{\mathrm{grain}}=3$ g cm*-3*. The adiabatic index of the gas is $\gamma=5/3$.

A first test with a uniform grid at level $\ell=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).

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 $\ell=7$. They are operated with one dust phase and with this same phase split with $\mathcal{N}\in\left[2,10\right]$ (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 $\mathcal{N}=5$ (top) and $\mathcal{N}=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.

Figure 10 shows the CPU time as a function of the number of species $\mathcal{N}$. We see that the multigrain simulations are not very expensive. The CPU time agrees more with a square root scaling with $\mathcal{N}$ 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.

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^{-4}$, 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 $\ell=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).

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 $\ell=8$. As before, the intrinsic grain density is $3$ g cm*-3*, but every species is characterized by a proper dust ratio $\epsilon_{j}$ and grain size $s_{\mathrm{grain},j}$. These quantities are determined using a MRN-like distribution 333$\epsilon\left(s_{\mathrm{grain}}\right)\propto s_{\mathrm{grain}}^{4-m}$ with $m=3.5$, 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.

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 $\approx 10^{-19}$ g cm*-3*(Andre et al., 2000) which is approximately $10^{5}-10^{6}$ years. An efficient settling is thus expected to happen during the collapse of this cloud, especially for large grains, e.g., $s_{\mathrm{grain}}>10^{-2}\mathrm{c}\mathrm{m}$ here, for which the typical settling timescale is a few 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 $\epsilon_{1}$ 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 $\epsilon_{1}$ at the dust front, which emphasizes the fundamental character equivalent stopping time in multigrain simulations.

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_{0}$, a uniform density $\rho_{0}$ and a dust ratio $\epsilon_{0}$. Its radius $R_{0}$ is set by the ratio between the thermal and gravitational energies $\alpha$

[TABLE]

where $k_{\mathrm{B}}$ is the Boltzmann constant, $\mu$ is the mean molecular weight, $m_{\mathrm{H}}$ the hydrogen atomic mass, and $T_{\mathrm{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

[TABLE]

the evolution becoming adiabatic above the critical density $\rho_{\mathrm{ad}}=10^{-13}$ g cm*-3* (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 $\rho_{\mathrm{grain}}=3$ g cm*-3* which is a good approximation for a combination of carbonaceous and silicate grains.

5.2 Setup

All the models are set with $M_{0}=1M_{\odot}$, $\mu=2.31$, $\gamma=5/3$, and $T_{\mathrm{g}}=10$ K. AMR grids of $\ell_{\mathrm{min}}=5$ and $\ell_{\mathrm{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 $\mathbf{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*-3*.

5.3 Validity of the diffusion approximation

The diffusion approximation is valid as long as the stopping time $t_{\mathrm{s}}$ is shorter than the dynamical timescale which is, for gravitational collapses, roughly the free-fall time

[TABLE]

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

[TABLE]

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 $\simeq 4\leavevmode\nobreak\ \mathrm{m}\mathrm{m}$. 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

[TABLE]

For a perfectly coupled gas and dust mixture ($t_{\mathrm{s}}\ll t_{\mathrm{ff}}$) it writes

[TABLE]

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

[TABLE]

Physically, this is due to the fact that two cores with the same initial $\alpha$ but different $\epsilon_{0}$ have different initial radius, hence free-fall timescale. Indeed, $R_{0}$ increases with $\epsilon_{0}$ 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 $\alpha=0.5$ are performed. The condition $t_{\mathrm{s}}\ll t_{\mathrm{ff}}$ is ensured by considering very small grains with $s=10^{-7}$ 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 $\rho_{\mathrm{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.

5.5 Core properties

Non-rotating dustycollapse of $\alpha=0.5$ are performed with an initial dust-to-gas ratio of $1\%$ (corresponding to $\epsilon\sim 0.00990099$). Single grain species of size $1\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (Col1), $10\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (Col10) and $100\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (Col100) are considered. The results are compared, when the gas maximum density reaches $10^{-11}$ g cm*-3*, with a fiducial collapse without any treatment of dust (ColGAS).

Figure 16 shows the total mass of the first hydrostatic core $M_{\mathrm{core}}$, the dust mass ratio $M_{\mathrm{d}}/M_{\mathrm{core}}$ and its formation time $t_{\mathrm{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\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{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\geq 100\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{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_{\mathrm{d}}}/{M_{\mathrm{core}}}$.

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 $\mathbf{v_{\mathrm{g}}}\approx\mathbf{v_{\mathrm{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 $\sim 20\%$ for the Col10 test, it goes up to $\sim 300\%$ for the Col100 test. We see in the Col100 case that the 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.

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\simeq 100\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{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\leavevmode\nobreak\ \mathrm{\SIUnitSymbolMicro}\mathrm{m}$). However, they depend on the presence of dust itself as the thermal-to-gravitational energy ratio $\alpha$ 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.

Acknowledgements.

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 n°ANR-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.