Small dust grain dynamics on adaptive mesh refinement grids

Context . Small dust grains are essential ingredients of star, disk and planet formation. 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 ﬁnite volume codes. We use it to investigate dust dynamics during the protostellar collapse. Methods . We present a method to solve the monoﬂuid equations of gas and dust mixtures with several dust species in the di ﬀ usion approximation implemented in the adaptive-mesh-reﬁnement code RAMSES . It uses a ﬁnite volume second-order Godunov method with a predictor-corrector MUSCL scheme to estimate the ﬂuxes 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 ﬁrst-order on non-uniform grids. We apply our method on various dustycollapse simulations of 1 M (cid:12) cores composed of gas and dust. Conclusions . We developed an e ﬃ cient approach to treat gas and dust dynamics in the di ﬀ usion 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 (cid:39) 100 µ m fall signiﬁcantly faster than the gas. ISM: kinematics and dynamics – hydrodynamics – stars: formation – protoplanetary disks – methods: numerical


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 2014Bate , 2015;;Stoyanovskaya et al. 2017Stoyanovskaya et al. , 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.

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 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) t s,grain ≡ πγ 8 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.

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 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)ρ g e g , e g being the gas specific internal energy.The gas derivative ).We also define the stopping time as and the drag coefficient K as

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: We also define the dust ratio which must not be confused with the dust-to-gas ratio ρ d ρ g .Using these new quantities, it is straightforward to write After some developments, we write the conservation of mass and momentum for the mixture and of the gas internal energy Here only one Lagrangian derivative d dt ≡ ∂ ∂t + (v • ∇) 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).

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, where St is the Stokes number 1 .System (6) reduces to The errors caused by the diffusion approximation are marginal as long as St 1.We note that the first-order term − P g ρ(1− ) ∇ • ( ∆v) 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) 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 • ∇ e g .
1 ||.|| indicates the L2 norm of either a tensor or a vector.
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 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.

(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 (N + 1) phase mixtures with 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 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 where t s,k is the individual stopping time of the phase k.
The N l=1 l 1− l t 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 and the mean stopping time The gas and dust densities are simply given by When N = 1, Eqs. (9) reduce to the formulation for a single dust phase.
A96, page 3 of 17 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 where ρ d,k is the density of the dust phase k, E ≡ 1 2 ρv 2 +ρ(1−E)e g is the total energy of the mixture and I is the identity matrix.

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.
For simplicity, the hydrodynamical version of RAMSES is presented here.It solves the Euler equations of a pure gas in their hyperbolic form the state U vector and flux F vector being where E g ≡ 1 2 ρ g v 2 g + ρ g e g is the total energy of the gas.

Godunov scheme
Finite volume methods are based on the estimation of the average of U over the cells.For a 1D problem integrated over time (for t ∈ [t 1 , t 2 ]) and space (for x ∈ [x 1 , x 2 ]), the previous system writes 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 A(x, t) → A n i .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 where the discretized fluxes F n+1/2 i±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.

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 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 ∆t = ∆t 1, +1 + ∆t 2, +1 .
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 .

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 A96, page 4 of 17 same structure for the mass, momentum and energy conservation equation.A total of 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 The new state vector U writes the flux F H is similar to the flux of pure hydrodynamics and is given by the flux F ∆ accounts for dust diffusion and is given by where w k , the differential advection velocity between the dust species k and the barycenter, is and w g,b is the differential velocity between the gas and the barycenter that writes 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 F 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

MUSCL scheme for dust diffusion/advection
The second step of the operator splitting consists in taking into account the second flux vector F ∆ .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 where ρ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 F (ρ d ) ≡ ρ d 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 ρ is constant.

Differential advection velocity.
During the diffusion step, we aim to compute the differential advection velocity w n i, j,k .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 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 secondorder 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 xcomponent of w n i, j,k is given by 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 where ∆ σ ρ d n i, j,k and ∆ σ w σ n i, j,k are the TVD variations of ρ d and w σ in the direction σ.After the prediction of ρ d n+1/2 i, j,k , 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 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 A96, page 5 of 17 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 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 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 F n+1/2 i, j−1/2,k and F n+1/2 i, j,k−1/2 are estimated in a similar way.The dust density is finally updated according to Energy.Similarly and simultaneously with the dust diffusion step, the energy is updated after the hydrodynamical step to account for the energy component of F ∆ .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.

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 where C CFL < 1 is a safety factor.In addition, the hydrodynamical step imposes another stability condition that writes where ν σ is the mixture velocity in the direction σ.Dust diffusion is stable without intervention on the timestep as long as 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.

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.

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 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 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 ∀x ∈ [0, L] 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 −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, ∀x ∈ [0, L], We use a very small timestep ∆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 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.Fig. 2. dustyadvect tests with an initial condition given by the function f 2 (Eq.( 31)).L 2 (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).

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: An isothermal equation of state P g = c 2 s (1 − )ρ is considered.The advection velocity is then given by ρ being constant, the equation can be written as the a non-linear diffusion equation Equation ( 35) has a self similar solution known as the Barenblatt-Pattle solution (Barenblatt 1952) that writes where C is a constant depending on the initial conditions.We consider the following initial profile which is consistent with Eq. (36) if • We additionally set ρ = 1, t s = 0.1, 0 = 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 2 norm.Even though our scheme is fundamentally designed for advection, it is also efficient at handling diffusion problems.

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 cs cs = c s 1 − 0 , 0 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 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 ρ 0 , v 0 , P g 0 , 0 left = (1, 0, 1, 0.5) ρ 0 , v 0 , P g 0 , 0 right = (0.125, 0, 0.1, 0.5) .
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.

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 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 cs 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 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 The adiabatic index of the gas is γ = 1.000001 2 to simulate an isothermal soundwave propagation.The initial perturbation has a relative amplitude δ = 10 −4 .The simulation box has a size L = 1.0 and the grid is taken as uniform with 64 cells.The timestep ∆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 ν x and the perturbation density δρ = ρ − ρ 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 = 5 to = 9, with a small timestep ∆t = 10 −5 and for K = 50.We also perform runs on AMR grids 2 Taking exactly γ = 1 is not possible with RAMSES as the internal energy is computed as 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 −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.hydrostatic equilibrium.As in Price & Laibe (2015), we set an analytic gravitational force where G 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 where H is the local scale height of the disk, which is determined by the ratio H/r = 0.05.We then assume an isothermal equation of state where the imposed soundspeed c s is Ω k being the Keplerian angular velocity at the radius r.Finally, a uniform initial dust ratio 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 which is the limit at low Stokes number (St = Ω k t 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 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 L 2 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).
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 8H.As in Price & Laibe (2015) and Hutchison et al. (2018), the initial mid-plane gas density ρ g ≈ 6 × 10 −13 g cm −3 .We also set an initial dust ratio of 4.99 × 10 −3 of millimeter grains with an intrinsic density of ρ grain = 3 g cm −3 .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).

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 N ∈ [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 N = 5 (top) and 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 N. We see that the multigrain simulations are not very expensive.The CPU time agrees more with a square root scaling with 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.

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

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 A96, page 11 of 17 a resolution of = 8.As before, the intrinsic grain density is 3 g cm −3 , 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 3 , 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 3 s grain ∝ s 4−m grain with m = 3.5.3.99 × 10 −3 400 000 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 ≈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 grain > 10 −2 cm 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 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 onespecies 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 1 at the dust front, which emphasizes the fundamental character equivalent stopping time in multigrain simulations.

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.

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 ρ 0 and a dust ratio 0 .Its radius R 0 is set by the ratio between the thermal and gravitational energies α α = 5 2 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 the evolution becoming adiabatic above the critical density ρ 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 ρ grain = 3 g cm −3 which is a good approximation for a combination of carbonaceous and silicate grains.

Setup
All the models are set with M 0 = 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 −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 (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).

Free-fall timescale for strongly coupled mixtures
A core only composed with gas would collapse at the following free-fall timescale t ff,g ≡ 3π 32Gρ g • For a perfectly coupled gas and dust mixture (t s t ff ) it writes In the context of the Boss and Bodenheimer test, this timescale indirectly depends on the initial dust ratio Physically, this is due to the fact that two cores with the same initial α but different 0 have different initial radius, hence free-fall A96, page 13 of 17 timescale.Indeed, R 0 increases with 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 α = 0.5 are performed.The condition t s t 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 ρ 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.
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 Col-GAS 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; Dust ratio 1 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 onespecies solution.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.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 .
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.
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.

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

Fig. 1 .
Fig.1.dustyadvect tests using the function f 1 (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.

Fig. 4 .Fig. 5 .
Fig. 4. Dustyshock with = 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.

Fig. 7 .
Fig.7.dustywave numerical convergence tests.L 2 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).

Fig. 10 .
Fig. 10.CPU time t CPU of the ten equivalent dustydiffuse tests as a function of the number of species N.

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

Fig. 12 .
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.

Fig. 16 .
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 St 0 .

Fig. 17 .
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%.

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