# MULTIGRAIN: a smoothed particle hydrodynamic algorithm for multiple small dust grains and gas

MULTIGRAIN: a smoothed particle hydrodynamic algorithm for multiple small dust grains and gas Abstract We present a new algorithm, multigrain, for modelling the dynamics of an entire population of small dust grains immersed in gas, typical of conditions that are found in molecular clouds and protoplanetary discs. The multigrain method is more accurate than single-phase simulations because the gas experiences a backreaction from each dust phase and communicates this change to the other phases, thereby indirectly coupling the dust phases together. The multigrain method is fast, explicit and low storage, requiring only an array of dust fractions and their derivatives defined for each resolution element. hydrodynamics, methods: numerical, protoplanetary discs, dust, extinction, ISM: kinematics and dynamics 1 INTRODUCTION Modelling the interaction of multiple dust grains simultaneously with the gas is a ‘grand challenge’ in protoplanetary disc modelling (Haworth et al. 2016), since discs involve grains with sizes spanning several orders of magnitude, from sub-micron grains to km-sized planetesimals. Grains of different sizes experience different dynamics since small grains are lighter and more easily influenced by the gas compared to larger, heavier grains. The usual approach to dusty gas dynamics is to model the gas and dust as separate fluids. The gas is modelled either on a grid (Paardekooper & Mellema 2004; Youdin & Johansen 2007; Balsara et al. 2009; Bai & Stone 2010a; Miniati 2010; Yang & Johansen 2016) or on a set of Lagrangian particles (Monaghan & Kocharyan 1995; Monaghan 1997; Barrière-Fouchet et al. 2005; Laibe & Price 2012a,b; Lorén-Aguilar & Bate 2014, 2015); similarly for each dust phase (although the discretization method often defaults to the one used by the gas). During simulation, the gas and dust fluids are evolved independently, but interact via a common drag force (e.g. Saffman 1962; Garaud & Lin 2004). Although grid- and particle-based methods each have their own distinct advantages/disadvantages (e.g. Price & Federrath 2010), they both require prohibitively small time-steps or implicit methods at high drag. Furthermore, Laibe & Price (2012a,b) discovered a drag resolution criterion that becomes increasingly restrictive with smaller grain sizes and applies generally to any method that models dust on a grid or on a set of particles that is not collocated with the gas at all times. While Laibe & Price (2012a,b), and later Lorén-Aguilar & Bate (2014), tested this spatial criterion using smoothed particle hydrodynamics (SPH), Youdin & Johansen (2007) inferred a similarly high-resolution requirement in hybrid grid-particle simulations. Failing to meet this criterion may explain the first-order convergence rate in high drag regimes observed by Miniati (2010), Bai & Stone (2010a), and Yang & Johansen (2016). To address the restrictive temporal and spatial restrictions that exist for high drag regimes, Laibe & Price (2014a,b,c, hereafter LP14a; LP14b; LP14c) and Price & Laibe (2015, hereafter PL15) developed a single-fluid formulation appropriate for small grains – similar to earlier formulations by Johansen & Klahr (2005). The dust–gas mixture is advected at the barycentric velocity and whose density is equal to the total density of the mixture. In the context of SPH, this means the mixture is represented by a single set of SPH particles with an evolution equation for the dust fraction (LP14b; PL15). While the above methods provide a means of modelling discs or molecular clouds with a single embedded dust phase, the challenge is to span the observed range of grain sizes. The typical approach is the one we recently used in Dipierro et al. (2015), where a series of single-phase simulations were stitched together in post-processing to interpret the dark structures observed at millimetre wavelengths by the ALMA interferometer in the disc surrounding the star HL Tau. In that paper, the method from Laibe & Price (2012a) was used to model the dynamics of mm-sized grains and larger, while the smaller grains were modelled using the method from PL15. Besides being tedious, the procedure used by Dipierro et al. (2015) is slow and, more importantly, neglects the indirect coupling between dust phases caused by the ‘backreaction’ of individual phases on the gas, which in turn influences the grain dynamics. Neglecting this backreaction misses important effects such as outward migration of dust particles (Bai & Stone 2010a) and/or modification of the linear growth rate of the streaming instability LP14c. Also, backreaction in individual grain size simulations is both annoying and wrong – annoying because the different response of the gas makes stacking of different dust grain distributions difficult (Tricco, Price & Laibe 2017); wrong because the gas should respond to the entire dust mixture, rather than each grain size individually. In this paper, we develop a new multigrain algorithm for modelling the dynamics of multiple dust phases, based on the analytical work presented in LP14c. Because much of the opacity and accompanying scattering/emission in astrophysical environments stems from the presence of dust grains that can be considered ‘small’ (i.e. where the terminal velocity approximation is valid), we focus on deriving and implementing the SPH versions of the continuum equations for the multiphase, terminal velocity approximation – generalizing the single-dust-phase method developed in PL15. 2 THE DIFFUSION APPROXIMATION FOR MULTIPLE DUST SPECIES 2.1 Continuum equations We consider a system consisting of a mixture of a single gas phase and N strongly coupled dust phases. Throughout this paper, we use the indices a, b, and c to refer to individual simulation particles that move at the barycentric velocity of the mixture. Subscript or superscripts g and d are used for gas and dust properties, respectively. Finally, we identify the fluid quantities for each of the N different dust phases using the index j. 2.1.1 General equations LP14c derived the general continuum fluid equations for a mixture of gas and N coupled dust species moving in a barycentric reference frame. They further showed that in strongly coupled regimes – i.e. first order in tj/T, where tj is a drag time-scale specific to each grain type (see equation 16; note the difference in notation from that of LP14c) and T is the time-scale for a sound wave to propagate over a typical distance L (commonly referred to as the terminal velocity approximation; see e.g. Youdin & Goodman 2005; Chiang 2008; Barranco 2009; Lee et al. 2010; Jacquet, Balbus & Latter 2011) – the fluid equations reduce to   \begin{eqnarray} \frac{{\rm d} \rho }{{\rm d} t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (1)  \begin{eqnarray} \frac{{\rm d} \epsilon _{j}}{{\rm d} t} = - \frac{1}{\rho } \nabla \cdot \left[ \rho \epsilon _{j}\left( \Delta \boldsymbol {v}_{j}- \epsilon \Delta \boldsymbol {v}\right) \right], \end{eqnarray} (2)  \begin{eqnarray} \frac{{\rm d} \boldsymbol {v}}{{\rm d} t} = \left( 1 - \epsilon \right) \boldsymbol {f}_{\mathrm{g}}+ \sum _{j}\epsilon _{j}\boldsymbol {f}_{\mathrm{d}j}+ \boldsymbol {f}, \end{eqnarray} (3)  \begin{eqnarray} \frac{{\rm d} u}{{\rm d} t} = - \frac{P}{\rho _{\mathrm{g}}} \nabla \cdot \boldsymbol {v}+ \epsilon \Delta \boldsymbol {v}\cdot \nabla u, \end{eqnarray} (4)  \begin{eqnarray} \Delta \boldsymbol {v}_{j} = \big[\Delta \boldsymbol {f}_{j}- \sum _{k}\epsilon _{k}\Delta \boldsymbol {f}_{k}\big] \epsilon _{j}t_{j}, \end{eqnarray} (5)where d/dt is the convective derivative using the barycentric velocity $$\boldsymbol {v}$$,   \begin{eqnarray} \boldsymbol {v}\equiv \displaystyle \frac{\rho _{\mathrm{g}}\boldsymbol {v}_{\mathrm{g}}+ \displaystyle \sum _{j}\rho _{\mathrm{d}j}\boldsymbol {v}_{\mathrm{d}j}}{\rho } = \frac{\rho _{\mathrm{g}}\boldsymbol {v}_{\mathrm{g}}+ \rho _{\mathrm{d}}\boldsymbol {v}_{\mathrm{d}}}{\rho }, \end{eqnarray} (6)ρ is the total density of the mixture,   \begin{eqnarray} \rho \equiv \rho _{\mathrm{g}}+ \rho _{\mathrm{d}}= \rho _{\mathrm{g}}+ \sum _{j}\rho _{\mathrm{d}j}, \end{eqnarray} (7)εj and ε are the mass fractions (relative to the mixture) of the individual and combined dust phases, respectively,   \begin{eqnarray} \epsilon _{j} \equiv \frac{\rho _{\mathrm{d}j}}{\rho }, \end{eqnarray} (8)  \begin{eqnarray} \epsilon \equiv \sum _{j}\epsilon _{j}= \frac{\rho _{\mathrm{d}}}{\rho }, \end{eqnarray} (9)$$\Delta \boldsymbol {v}$$ is the weighted sum of the differential velocities $$\Delta \boldsymbol {v}_{j}\equiv \boldsymbol {v}_{\mathrm{d}j}-\boldsymbol {v}_{\mathrm{g}}$$,   \begin{eqnarray} \Delta \boldsymbol {v}\equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}\Delta \boldsymbol {v}_{j}, \end{eqnarray} (10)$$\boldsymbol {f}$$ represents accelerations acting on both components of the fluid, while $$\boldsymbol {f}_{\mathrm{g}}$$ and $$\boldsymbol {f}_{\mathrm{d}j}$$ represent the accelerations acting on the gas and dust components, respectively, $$\Delta \boldsymbol {f}_{j}\equiv \boldsymbol {f}_{\mathrm{d}j}- \boldsymbol {f}_{\mathrm{g}}$$ is the differential force between the gas and each dust phase, u is the specific thermal energy of the gas, and P is the gas pressure. 2.1.2 Drag time-scales When N = 1, the drag time-scale is unambiguously set by the drag stopping time,   \begin{eqnarray} t_{\mathrm{s}}^{N=1}\equiv \displaystyle\frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}}{K \rho }, \end{eqnarray} (11)where K is a drag coefficient that, in general, depends on local properties of the gas and dust. We assume that K is either constant or in the linear Epstein regime, suitable for small dust grains with low Mach numbers (Epstein 1924, also e.g. Laibe & Price 2012b). In the latter case,   \begin{eqnarray} K = \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}}{\rho _{\rm grain}s} \sqrt{\frac{8}{\pi \gamma }} c_\mathrm{s}= \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}c_\mathrm{s}}{\rho _\mathrm{eff}s}, \end{eqnarray} (12)where we assume spherical grains with radius s, with uniform intrinsic dust density ρgrain, or equivalently, an effective density $$\rho _\mathrm{eff}\equiv \rho _{\rm grain}\sqrt{\pi \gamma / 8}$$. As usual, γ is the adiabatic constant. The stopping time for N = 1 in the Epstein regime can therefore be written as   \begin{eqnarray} t_{\mathrm{s}}^{N=1}= \frac{\rho _\mathrm{eff}s}{\rho c_\mathrm{s}}. \end{eqnarray} (13) Generalizing the stopping time to N > 1 is conceptually simple, but difficult in practice. Each dust-type equilibrates with the gas at a different rate depending on both the intrinsic properties of the dust grains and the local properties of the gas. Although we assume dust grains of different species do not interact, they are indirectly coupled by their mutual backreaction on the gas. One approach is to derive time-scales using the eigenvalues of the drag matrix LP14c, but the derivations and the expressions become increasingly unwieldy as N increases (i.e. there is no general algebraic expression as a function of N). The eigenvalues help aid in interpreting results, but they are not needed to evolve the fluid equations numerically. The only potential impact the eigenvalues have is through their influence on the time-step. Even then, LP14c found fixed upper/lower bounds to the eigenvalues of the N × N drag matrix, effectively removing any need for the eigenvalues during computation. For convenience, we define the following time-scales to help simplify our numerical implementation:   \begin{eqnarray}{2} T_{\mathrm{s}j} \equiv \frac{\rho _\mathrm{eff}s_j}{\rho c_\mathrm{s}} = \epsilon _{j}\left( 1 - \epsilon \right) t_{j}, \end{eqnarray} (14)  \begin{eqnarray}{2} \widetilde{T}_{\mathrm{s}j} \equiv \frac{T_{\mathrm{s}j}- \sum _{k}\epsilon _{k}T_{\mathrm{s}k}}{1 - \epsilon } \quad && = \epsilon _{j}t_{j}- \sum _{k}\epsilon _{k}^2 t_{k}, \end{eqnarray} (15)where   \begin{eqnarray} t_{j}\equiv \frac{\rho }{K_j}, \end{eqnarray} (16)and where Kj is the drag coefficient for each dust phase, e.g.   \begin{eqnarray} K_j = \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}j}c_\mathrm{s}}{\rho _\mathrm{eff}s_j}. \end{eqnarray} (17)Note that the weighted sums of equations (14) and (15) happen to be equivalent, i.e.   \begin{eqnarray} \frac{1}{\epsilon } \sum _{j}\epsilon _{j}T_{\mathrm{s}j}\, = \, \frac{1}{\epsilon } \sum _{j}\epsilon _{j}\widetilde{T}_{\mathrm{s}j}\, = \, \frac{1 - \epsilon }{\epsilon } \sum _{j}\epsilon _{j}^2 t_{j}. \end{eqnarray} (18)This new quantity carries physical significance, but its interpretation is clearer if we first define an effective grain size for the mixture,   \begin{eqnarray} s \equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}s_j, \end{eqnarray} (19)such that equation (18) can be written in a more familiar form:   \begin{eqnarray} T_{\mathrm{s}}\equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}T_{\mathrm{s}j}= \frac{\rho _\mathrm{eff}s}{\rho c_\mathrm{s}}. \end{eqnarray} (20)Comparing this to equation (13), one may observe that Ts acts like an effective stopping time for the mixture. The benefit of using Tsj and $$\widetilde{T}_{\mathrm{s}j}$$ in lieu of tj is that they allow us to use our existing codebase with only a few additional lines of code, namely to assemble $$\widetilde{T}_{\mathrm{s}j}$$ (Tsj is calculated identically to $$t_{\mathrm{s}}^{N=1}$$ with s replaced by sj). In return, the form of the evolution equations are unchanged from the N = 1 case, as evidenced in the following sections. 2.1.3 Hydrodynamics For the simple case of hydrodynamics, the only force is the pressure gradient, i.e.   \begin{eqnarray} \boldsymbol {f}_{\mathrm{d}j} = 0, \end{eqnarray} (21)  \begin{eqnarray} \boldsymbol {f}_{\mathrm{g}} = - \frac{\nabla P}{\rho _{\mathrm{g}}}, \end{eqnarray} (22)  \begin{eqnarray} \Delta \boldsymbol {f}_{j} = \displaystyle\frac{\nabla P}{\rho _{\mathrm{g}}}. \end{eqnarray} (23)Using equations (14) and (23) to simplify equation (5), we get   \begin{eqnarray} \Delta \boldsymbol {v}_{j}= \frac{\epsilon _{j}t_{j}\nabla P}{\rho } = \frac{T_{\mathrm{s}j}\nabla P}{\rho _{\mathrm{g}}}, \end{eqnarray} (24)while equations (10), (20), and (24) allow us to write   \begin{eqnarray} \Delta \boldsymbol {v}= \frac{T_{\mathrm{s}}\nabla P}{\rho _{\mathrm{g}}}. \end{eqnarray} (25)As promised, when these last two expressions for $$\Delta \boldsymbol {v}_{j}$$ and $$\Delta \boldsymbol {v}$$ are inserted into equations (1), (2), (3), and (4), we obtain the same form of the fluid equations as reported in PL15 for the N = 1 terminal velocity approximation, namely   \begin{eqnarray} \frac{{\rm d} \rho }{{\rm d} t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (26)  \begin{eqnarray} \frac{{\rm d} \epsilon _{j}}{{\rm d} t} = - \frac{1}{\rho } \nabla \cdot \left( \epsilon _{j}\widetilde{T}_{\mathrm{s}j}\nabla P \right), \end{eqnarray} (27)  \begin{eqnarray} \frac{{\rm d} \boldsymbol {v}}{{\rm d} t} = -\frac{\nabla P}{\rho } + \boldsymbol {f}, \end{eqnarray} (28)  \begin{eqnarray} \frac{{\rm d} \tilde{u}}{{\rm d} t} = - \frac{P}{\rho } \nabla \cdot \boldsymbol {v}, \end{eqnarray} (29)where for convenience we have defined $$\tilde{u} \equiv (1 - \epsilon ) u$$ instead of evolving u directly as in PL15. The corresponding energy equation in terms of u would be   \begin{eqnarray} \frac{{\rm d} u}{{\rm d} t} = - \frac{P}{\rho _{\mathrm{g}}} \nabla \cdot \boldsymbol {v}+ \frac{\epsilon T_{\mathrm{s}}}{\rho _{\mathrm{g}}} \nabla P \cdot \nabla u. \end{eqnarray} (30) In order to recover the special case of a single dust phase, we need only collapse the sums in equations (14) and (15) and set sj → s and εj → ε. It is simple to check that in this limit, $$T_{\mathrm{s}j}= \widetilde{T}_{\mathrm{s}j}= T_{\mathrm{s}}= t_{\mathrm{s}}^{N=1}$$, thereby recovering the N = 1 fluid equations from PL15 exactly. 2.1.4 Equation of state The set of equations above is closed by assuming the usual equation of state, which constrains the gas pressure P in terms of the gas density and temperature. Unless otherwise specified in this paper, we assume an adiabatic equation of state, i.e.   \begin{eqnarray} P = \left( \gamma - 1 \right) \rho _{\mathrm{g}}u = \left( \gamma -1 \right) \left( 1 - \epsilon \right) \rho u, \end{eqnarray} (31)or simply   \begin{eqnarray} P = (\gamma -1) \rho \tilde{u}. \end{eqnarray} (32) 2.2 Time-stepping As pointed out by PL15, the addition of the dust evolution equation adds a further constraint on the time-step that becomes limiting when the diffusion coefficient is large. We can derive this time-step constraint more rigorously than that presented by PL15, albeit with the same result for N = 1, by discretizing the set of equations in time using a forward Euler method   \begin{eqnarray} \frac{\rho ^{n+1} - \rho ^n}{\Delta t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (33)  \begin{eqnarray} \frac{\epsilon _j^{n+1} - \epsilon _j^{n}}{\Delta t} = - \frac{1}{\rho } \nabla \cdot \left( \epsilon _{j}\widetilde{T}_{\mathrm{s}j}\nabla P \right), \end{eqnarray} (34)  \begin{eqnarray} \frac{\boldsymbol {v}^{n+1} - \boldsymbol {v}^n}{\Delta t} = -\frac{\nabla P}{\rho }, \end{eqnarray} (35)and performing a Von Neumann stability analysis on the above semidiscrete equations. That is, we solve the linear system that results from assuming plane wave solutions of the form   \begin{eqnarray} \rho = D {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (36)  \begin{eqnarray} \boldsymbol {v} = \boldsymbol {V} {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (37)  \begin{eqnarray} \epsilon _j = E_j {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (38)where D, $$\boldsymbol {V}$$, and Ej are perturbation amplitudes, k is the wavenumber, $$\boldsymbol {x}$$ is the position vector, and ω is the angular frequency. This analysis generically produces a time-step criterion of the form   \begin{eqnarray} \Delta t < C_0\frac{1}{k c_{\max }}, \end{eqnarray} (39)where C0 is a dimensionless safety factor of order unity and cmax  is the maximum wave speed according to the dispersion relation for linear waves. The wavelength of maximum growth usually occurs on the resolution scale, giving the usual Courant criterion   \begin{eqnarray} \Delta t < C_0\frac{h}{c_{\max }}, \end{eqnarray} (40)where h is the SPH smoothing length. For N = 1, the dispersion relation to first order in ωts is given by LP14a,  \begin{eqnarray} \omega = \pm \tilde{c}_{\rm s} k - \frac{i}{2} t_{\rm s} k^2 c_{\rm s}^2 \epsilon , \end{eqnarray} (41)where $$\tilde{c}^2_{\rm s} \equiv c_{\rm s}^2 (1 - \epsilon )$$ is the modified sound speed (squared). The maximum wave speed is therefore   \begin{eqnarray} c_{\max } =\left|\frac{\omega }{k}\right|= \sqrt{\tilde{c}_{\rm s}^2 + \frac{1}{4} \epsilon ^2 t_{\rm s}^2 k^2 c_{\rm s}^4}, \end{eqnarray} (42)and the time-step constraint appropriate for SPH is   \begin{eqnarray} \Delta t < C_0\frac{h}{\sqrt{\tilde{c}^2_{\rm s} + \epsilon ^2 t^2_{\rm s} c_{\rm s}^4 / h^2}}. \end{eqnarray} (43)This is similar to the time-step criterion proposed by PL15 except that the above combines the usual Courant–Friedrichs–Lewy (CFL) condition ($$\Delta t < h/\tilde{c}_{\rm s}$$) and the additional constraint from the dust evolution ($$\Delta t < h^2 / (\epsilon t_{\rm s} c_{\rm s}^2)$$) into a single criterion. When generalizing to multiple dust phases, we find the same result but with the effective stopping time replacing the N = 1 stopping time, giving   \begin{eqnarray} \Delta t < C_0\frac{h}{\sqrt{\tilde{c}^2_{\rm s} + \epsilon ^2 T_{\rm s}^2 c_{\rm s}^4 / h^2}}. \end{eqnarray} (44) As expected, with Ts in the denominator, restricting ourselves to strong drag regimes weakens the constraint on the time-step. More specifically, the time-step is limited when the grain-size distribution is dominated by large grains (or, alternatively, high dust fraction), such that   \begin{eqnarray} \frac{\epsilon T_{\mathrm{s}}}{1-\epsilon } > \Delta t_{\rm CFL}, \end{eqnarray} (45)where $$\Delta t_{\rm CFL}\equiv h/\tilde{c}_{\rm s}$$ is the CFL time-step. The added advantage of the criterion in equation (44) is that it is less stringent than the explicit time-step for either the full multigrain one-fluid formalism or the multifluid method (equations 79 and 80 of LP14c, respectively):   \begin{eqnarray} \Delta t_{\mathrm{one{\rm -}fluid}} < C \Bigg[ \max _{j}\left(\frac{1}{\epsilon _{j}t_{j}} \right) + \frac{1}{\left(1 - \epsilon \right)} \sum _{j}t_{j}^{-1} \Bigg]^{-1}, \end{eqnarray} (46)  \begin{eqnarray} \Delta t_{\mathrm{multi{\rm -}fluid}} < C \min _{j} \left[ \frac{1}{t_{j}} \left( \frac{1}{\epsilon _{j}} + \frac{1}{1-\epsilon } \right) \right]^{-1}, \end{eqnarray} (47)where C is another safety factor. Thus, as long as the cut-off to our dust distribution is ≲cm (see PL15), our global time-step should be of the order of ΔtCFL. 3 SPH FORMULATION When formulating the discretized SPH fluid equations, we can take advantage of the fact that (i) the only equations that were altered by having multiple dust phases were the dust fraction and energy equations and (ii) we have written the continuum equations in the same form as PL15. The first point allows us to adopt the discretized density and momentum equations from PL15 without any changes (thereby guaranteeing exact conservation of linear and angular momentum),   \begin{eqnarray} \rho _{a} = \sum\limits_{b} m_{b} W_{ab} (h_{a}), \end{eqnarray} (48)  \begin{eqnarray} \frac{{\rm d}\boldsymbol {v}_a}{{\rm d}t} & =& -\sum _{b} m_{b} \Big[ \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \nabla _{a} W_{ab} (h_{a}) \nonumber \\ &&+\,\frac{P_{b} + q^{\rm AV}_{ab, b}}{\Omega _{b} \rho _{b}^{2}} \nabla _{a} W_{ab} (h_{b}) \Big] + \boldsymbol {f}_a, \end{eqnarray} (49)where Wab is the usual SPH kernel, h is the smoothing length, Ω is the usual term to account for smoothing length gradients   \begin{eqnarray} \Omega _{a} = 1 - \frac{\mathrm{\partial} h_{a}}{\mathrm{\partial} \rho _{a}} \sum _{b} m_{b} \frac{\mathrm{\partial} W_{ab} (h_{a})}{\mathrm{\partial} h_{a}}, \end{eqnarray} (50)and h is related to ρ in the usual manner [which requires an iterative procedure to solve equation (48); see Price & Monaghan 2004, 2007; LP14b]. The second point allows us to write down the generalized diffusion equation for εj by inspection. Comparing equation (27) to (12) in PL15 suggests that we can use either of their discretized diffusion equations, provided we make the substitutions $$t_{\mathrm{s}}^{N=1}\rightarrow \widetilde{T}_{\mathrm{s}j}$$ and ε → εj (although in the latter case, care must be taken to leave any instances of the gas fraction, 1 − ε, untouched). Furthermore, because evolving the dust fraction directly can in some instances result in negative values, we prefer to use the positive definite formulation prescribed in appendix B of PL15 by defining $$S_j \equiv \sqrt{\rho \epsilon _{j}}$$ (not to be confused with the grain size sj). The corresponding evolution equation in terms of Sj is   \begin{eqnarray} \frac{\mathrm{d} S_{j,a}}{\mathrm{d}t} & =& - \frac{1}{2} \sum _{b} \frac{m_{b} S_{j,b}}{\rho _{b}} \left( \frac{\widetilde{T}_{\mathrm{s}j,a}}{\rho _{a}} + \frac{\widetilde{T}_{\mathrm{s}j,b}}{\rho _{b}} \right) \left(P_{a} - P_{b} \right) \frac{\overline{F}_{ab}}{\vert r_{ab} \vert } \nonumber \\ &&+\, \frac{S_{j,a}}{2\rho _{a}\Omega _{a}} \sum _{b} m_{b} \boldsymbol {v}_{ab} \cdot \nabla _a W_{ab} (h_{a}), \end{eqnarray} (51)where $$\overline{F}_{ab} \equiv \frac{1}{2}[ F_{ab}(h_{a}) + F_{ab}(h_{b}) ]$$ and Fab is defined such that $$\nabla _{a} W_{ab} \equiv F_{ab} \hat{\boldsymbol {r}}_{ab}$$. In writing the diffusion equation in this form, we have implicitly chosen to use the faster, easier-to-implement ‘direct second derivative’ method; however, the evolution equation for the ‘two first derivatives’ method can be obtained in the same fashion (see PL15 for a comparison of these two methods). 3.1 Conservation of energy This leaves only the energy equation to be determined. It is tempting to simply generalize the equation for the energy in a similar manner to the above, but conservation of energy puts an additional constraint on the form of the equation that is not immediately obvious. Instead, we derive the energy equation using the already discretized fluid equations above and by enforcing exact conservation of energy. The total energy E of the system in the terminal velocity approximation can be expressed as   \begin{eqnarray} E = \sum _a m_a \left( \frac{1}{2} v_a^2 + \tilde{u}_a \right), \end{eqnarray} (52)where $$\tilde{u}_a \equiv (1 - \epsilon _a) u_a$$ as previously. Conservation of energy requires that   \begin{eqnarray} \frac{\mathrm{d} E}{\mathrm{d}t} = \sum _{a} m_{a} \left[ \boldsymbol {v}_{a} \cdot \frac{\mathrm{d} \boldsymbol {v}_{a}}{\mathrm{d} t} + \frac{\mathrm{d} \tilde{u}_{a}}{\mathrm{d} t} \right]= 0, \end{eqnarray} (53)where   \begin{eqnarray} \frac{\mathrm{d} \tilde{u}_{a}}{\mathrm{d} t} = \frac{\rho ^\mathrm{g}_{a}}{\rho _a} \frac{\mathrm{d} u_{a}}{\mathrm{d} t} - u_{a} \sum _{j}\left( \frac{2 S_{j,a}}{\rho _{a}} \frac{\mathrm{d} S_{j,a}}{\mathrm{d} t} - \frac{S^2_{j,a}}{\rho ^2_{a}} \frac{\mathrm{d} \rho _{a}}{\mathrm{d}t} \right). \end{eqnarray} (54) Inserting the different expressions from equations (48), (49), and (51) and solving for the time derivative of the energy dictates that the discretized energy equation should be   \begin{eqnarray} \frac{{\rm d}\tilde{u}_{a}}{{\rm d}t} = \sum _{b} m_{b} \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \left( \boldsymbol {v}_{a} - \boldsymbol {v}_{b}\right) \cdot \nabla _{a} W_{ab} (h_{a}), \end{eqnarray} (55)or, if one evolves u directly as in PL15  \begin{eqnarray} \frac{{\rm d}u_{a}}{{\rm d}t} & =& \frac{1}{1 - \epsilon _{a}} \sum _{b} m_{b} \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \left( \boldsymbol {v}_{a} - \boldsymbol {v}_{b}\right) \cdot \nabla _{a} W_{ab} (h_{a}) \nonumber \\ &&- \, \frac{\rho _{a}}{2 \rho ^\mathrm{g}_{a}} \sum _{j}\sum _{b} m_{b} \frac{S_{j,a} S_{j,b}}{\rho _{a} \rho _{b}} \left( \frac{\widetilde{T}_{\mathrm{s}j,a}}{\rho _{a}} + \frac{\widetilde{T}_{\mathrm{s}j,b}}{\rho _{b}} \right) \nonumber \\ &&\times \, \left( u_{a} - u_{b}\right) \left( P_{a} - P_{b} \right) \frac{\overline{F}_{ab}}{\vert r_{ab} \vert }. \end{eqnarray} (56) 3.2 Shock-capturing terms We include the artificial viscosity and conductive terms below for completeness, but note that they are unchanged by the addition of more dust phases. 3.2.1 Artificial viscosity The artificial viscosity term is computed as follows:   \begin{eqnarray} q^{\rm AV}_{ab, a} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-\frac{1}{2} \left( 1 - \epsilon _{a} \right) v_{\mathrm{sig},a} \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab}, & \qquad \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} < 0\\ 0, & \qquad \mathrm{otherwise}, \end{array}\right. \end{eqnarray} (57)where $$\boldsymbol {v}_{ab} \equiv \boldsymbol {v}_{a} - \boldsymbol {v}_{b}$$ (similarly for $$\hat{\boldsymbol {r}}_{ab}$$) and the signal speed $$v$$sig corresponds to the usual choice for hydrodynamics, i.e.   \begin{eqnarray} v_{\mathrm{sig},a} = \alpha ^{\rm AV}_{a} c_{\mathrm{s},a} + \beta ^{\rm AV} \vert \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} \vert , \end{eqnarray} (58)where $$\alpha ^{\rm AV}_{a} \in [0,1]$$ is the linear dimensionless viscosity parameter (the index implying that αAV can be unique to each particle; see e.g. Morris & Monaghan 1997; Cullen & Dehnen 2010) and βAV (typically βAV = 2) is the von Neumann–Richtmyer viscosity parameter. 3.2.2 Artificial conductivity In order to correctly treat contact discontinuities, an artificial conductivity term must be added to the energy equations (see Price 2008),   \begin{eqnarray} \left( \frac{\mathrm{d} u_{a}}{\mathrm{d} t} \right)_{\rm cond} = \frac{1}{1-\epsilon _{a}} \sum _{b} m_{b} \left[ \frac{Q_{ab, a}}{\Omega _{a} \rho _{a}^{2}} F_{ab} (h_{a}) + \frac{Q_{ab, b}}{\Omega _{b} \rho _{b}^{2}} F_{ab} (h_{b}) \right],\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (59)where   \begin{eqnarray} Q_{ab,a} = \frac{1}{2} \alpha _u \rho _{a} v_{\mathrm{sig},u} \left( u_{a} - u_{b} \right), \end{eqnarray} (60)with αu ∈ [0, 1] the dimensionless conductivity parameter and $$v_{\mathrm{sig},u} = \vert \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} \vert$$ (Price 2008; Wadsley, Veeravalli & Couchman 2008). 4 NUMERICAL TESTS Given the similarity of the SPH equations in Section 3 to those in PL15 and the existing implementation of the latter in our SPH code phantom (e.g. Dipierro et al. 2015; Price et al. 2017), the generalization to N dust phases was straightforward. phantom is well tested (see Price et al. 2017) and we are confident that the implementation of the N = 1 terminal velocity approximation from which we started was correct. Therefore, the tests in this section are less focused on the code as a whole and more focused on specific aspects of our implementation. 4.1 Recovering the N = 1 case By far, the most difficult part of implementing the multigrain method was expanding the relevant arrays in the code to accommodate the N − 1 additional dust phases. To ensure that our new array structures cause no ill effects, we performed 3D versions of the dustyshock, dustywave, and dustydiffuse tests from PL15. Interested readers can find the set-up details for these tests in Price et al. (2017). The results from these tests are shown in the top row of Fig. 1. Importantly, we found that the results calculated with and without our new array structure matched to within machine precision. Note that this agreement shows that equations (48), (49), (51), and (56) reduce numerically to the N = 1 case, analogous to what we observed with the continuum equations. Figure 1. View largeDownload slide Results from the dustyshock (left column), dustywave (middle column), and dustydiffuse (right column) tests as performed by PL15 and Price et al. (2017), but with our new array structure capable of handling multiple dust phases. The top row shows results when N = 1, while the bottom row contains simulations where the same dust phase has been split into N = 10 equal mass bins. As desired, the two cases are identical. Moreover, they match the results run prior to the multigrain implementation. Figure 1. View largeDownload slide Results from the dustyshock (left column), dustywave (middle column), and dustydiffuse (right column) tests as performed by PL15 and Price et al. (2017), but with our new array structure capable of handling multiple dust phases. The top row shows results when N = 1, while the bottom row contains simulations where the same dust phase has been split into N = 10 equal mass bins. As desired, the two cases are identical. Moreover, they match the results run prior to the multigrain implementation. We then added an additional layer of complexity by splitting the single dust phase in each of the above tests into N equal bins and evolving them as if they were N different dust phases. This new set-up can be achieved by setting εj = ε/N and Kj = K. Separating the fluid into mass bins does not alter the physics of the problem, just the numerical method by which it is modelled. Therefore, we should recover the N = 1 solution (modulo numerical errors from calculating and combining quantities differently). Fig. 1 (bottom row) shows the results from the multigrain calculations. Again we found that the dustyshock, dustywave, and dustydiffuse tests agreed with the N = 1 cases to within machine precision. 4.2 Testing the general case It seems like the next logical test would be to extend one or more of the tests above to the general case of N different dust phases. However, there is a fundamental difference in the way the drag is calculated for these tests and the way we have assumed the drag will be calculated when using the equations derived in this paper. Whereas the tests above use a constant drag coefficient K for the entire fluid, the equations in Sections 2 and 3 are optimized for physical dust grains in the Epstein drag regime where the equivalent drag constant (17) changes with grain size. We could reformulate the tests and their solutions to accept a unique value of K for each dust phase, but this would require altering equations (51) and (56) – the very equations we are trying to verify. Therefore, for the general case, we need a test requiring physical grain sizes and drag coefficients. 4.3 Dust settling in a protoplanetary disc The dust settling test from PL15 is an ideal candidate for testing the general case because it mimics one of the environments the multigrain method is designed to simulate, namely the settling of small dust grains in protoplanetary discs. 4.3.1 Initial conditions We simulate a disc-like environment at a radius r = 50 au using a thin, vertical (Cartesian) column of gas in near-hydrostatic equilibrium with an external acceleration in the form of   \begin{eqnarray} \boldsymbol {a}_\mathrm{ext} = - \frac{\mathcal {G} M z}{\left(r^2+z^2\right)^{3/2}}\boldsymbol {\hat{z}}, \end{eqnarray} (61)where $$\mathcal {G}$$ is Newton's gravitational constant, M is the stellar mass, and z is the ‘vertical’ coordinate along the length of the column (x and y represent the two shorter dimensions of the column). The gas density of the column is given by   \begin{eqnarray} \rho _{\mathrm{g}}(z) = \rho _{\mathrm{g},0}\, \mathrm{exp} \left[-\frac{z^2}{2H^2}\right], \end{eqnarray} (62)where we choose H/r = 0.05, giving a disc scale height of H = 2.5 au. We assume an isothermal equation of state with $$P = c_\mathrm{s}^2 \rho _{\mathrm{g}}$$, where cs ≡ HΩ and $$\Omega \equiv \sqrt{\mathcal {G} M / r^3}$$, corresponding to an orbital time $$t_\mathrm{orb} \equiv 2\pi /\Omega \approx 353\, \mathrm{yr}$$. We adopt code units with a distance unit of 10 au, mass in solar masses and time units such that $$\mathcal {G} = 1$$. These choices give an orbital time of ≈70.2 in code units. The particles are initially placed on a close-packed lattice using 100 × 86 × 78 = 670 800 particles in the domain [x, y, z] ∈ [ ± 1, ±0.75, ±0.65]. We then stretch the particles in z using the method described in Price (2004) to give the density profile given in equation (62). We set ρg, 0 to 10−3 in code units (≈6 × 10−13 g cm−3 in physical units), corresponding to a particle mass in code units of 2.42 × 10−9. We use periodic boundary conditions in all directions, but set the boundary in z at ±10H in order to avoid periodicity in the vertical direction. We relaxed the density profile by running the code for 15 orbits with artificial viscosity, at which point we added N = 10 distinct dust phases to the system. We created a cell-edge, logarithmic grid from smin to smax with grid cells of width $$\Delta \log s = \frac{1}{N} \log _{10}\left( s_{\rm max}/s_{\rm min}\right)$$. Then we assigned sj by taking the square root of the product of the cell's endpoints, – thereby skewing the ‘typical’ grain size for each cell towards the smaller, more numerous dust grains. Each dust phase was distributed throughout the disc with an initially uniform dust fraction. We constrained the total dust fraction to be ε = 1/101 (corresponding to a dust-to-gas ratio of 0.01) and set the magnitudes of εj according to the differential power-law distribution   \begin{eqnarray} \mathrm{d} \epsilon = \epsilon _0 s^{3-p} \mathrm{d} s, \qquad \mathrm{for} \quad s_{\rm min}\le s \le s_{\rm max}, \end{eqnarray} (63)where dε is the differential dust fraction with respect to grain size, ε0 is a normalization factor, and p is the usual power-law index for number density as a function of grain size (e.g. Mathis, Rumpl & Nordsieck 1977). In particular, εj is determined by integrating equation (63) across each grain-size cell and then normalizing their combined sum via equation (9). Assuming p = 3.5, we set $$s_{\rm min}\approx 0.0599\, \mu \mathrm{m}$$ and smax ≈ 1.67 mm such that the smallest simulated grain size is $$0.1\, \mu \mathrm{m}$$ and the largest simulated grain size is 1 mm. The initial values for sj and εj in this test are listed in Table 1. Table 1. The initial values for sj and εj used in the settling test assuming a power-law distribution in grain sizes ranging from smin = 0.1 μm to smax = 1 mm with a power-law index of p = 3.5. j  sj [cm]  εj  1  1.00 × 10−5  3.99 × 10−5  2  2.78 × 10−5  6.65 × 10−5  3  7.74 × 10−5  1.11 × 10−4  4  2.15 × 10−4  1.85 × 10−4  5  5.99 × 10−4  3.09 × 10−4  6  1.67 × 10−3  5.15 × 10−4  7  4.64 × 10−3  8.59 × 10−4  8  1.29 × 10−2  1.43 × 10−3  9  3.59 × 10−2  2.39 × 10−3  10  1.00 × 10−1  3.99 × 10−3  j  sj [cm]  εj  1  1.00 × 10−5  3.99 × 10−5  2  2.78 × 10−5  6.65 × 10−5  3  7.74 × 10−5  1.11 × 10−4  4  2.15 × 10−4  1.85 × 10−4  5  5.99 × 10−4  3.09 × 10−4  6  1.67 × 10−3  5.15 × 10−4  7  4.64 × 10−3  8.59 × 10−4  8  1.29 × 10−2  1.43 × 10−3  9  3.59 × 10−2  2.39 × 10−3  10  1.00 × 10−1  3.99 × 10−3  View Large 4.3.2 Results After adding the dust, we ran the simulation for an additional 15 orbits. The resulting dust density for each of the different phases is shown in Fig. 2. As expected, the settling efficiency is proportional to the size of the dust grains, thus enhancing the mid-plane density of the larger grains. However, visually separating this density enhancement is difficult in Fig. 2 because the initial density distribution increases by a factor of 100 from smin to smax. Although the continuum density distribution is a decreasing function with respect to grain size ($$\propto s^{3-p} = 1/\sqrt{s}$$), integrating over each cell to include the mass from non-simulated grains steepens the power law by an additional power of s such that $$\rho _{\mathrm{d}j}\propto \sqrt{s}$$. Figure 2. View largeDownload slide Ten dust densities from a multigrain simulation after having settled for 15 orbits in a 3D vertical column of a protoplanetary disc at r = 50 au (assuming H/r = 0.05; so H = 2.5 au) using 100 × 86 × 78 = 670 800 simulation particles. The grain size and initial dust fraction for each phase is listed in Table 1. Large dust grains efficiently settle towards the disc mid-plane, but still have a much lower density than the smaller dust grains because the global number density of the larger grains is lower. Our multigrain simulation is ∼5 × faster to run than 10 single-phase simulations run serially (see Section 4.3.3). Figure 2. View largeDownload slide Ten dust densities from a multigrain simulation after having settled for 15 orbits in a 3D vertical column of a protoplanetary disc at r = 50 au (assuming H/r = 0.05; so H = 2.5 au) using 100 × 86 × 78 = 670 800 simulation particles. The grain size and initial dust fraction for each phase is listed in Table 1. Large dust grains efficiently settle towards the disc mid-plane, but still have a much lower density than the smaller dust grains because the global number density of the larger grains is lower. Our multigrain simulation is ∼5 × faster to run than 10 single-phase simulations run serially (see Section 4.3.3). In order to better show how the mid-plane density is affected by settling, we set-up and ran a second simulation where the dust fractions were all equal, i.e. εj = ε/N. Fig. 3 shows the resulting time evolution of the dust density for phases j = [1, 9, 10]. This time we clearly see that settling increases the dust density relative to its initial state and at a rate that is commensurate with its settling efficiency. These results are in good agreement with previous settling tests performed in the literature (PL15; Hutchison et al. 2016; Price et al. 2017), albeit with a smaller initial dust fraction. Figure 3. View largeDownload slide Time evolution of the densities of three dust phases (j = [1, 9, 10]). The initial conditions in this simulation were the same as in Fig. 2, except with equal dust fractions (εj = ε/N) to make the relative density enhancement within and between dust phases more visible. We have also adjusted the colourbar in order to allow direct comparison with the settling tests performed by PL15 and Price et al. (2017). Note that the density enhancement due to settling has a shallower dynamic range than the built-in density gamut created by our grains-size distribution (see Fig. 2). Figure 3. View largeDownload slide Time evolution of the densities of three dust phases (j = [1, 9, 10]). The initial conditions in this simulation were the same as in Fig. 2, except with equal dust fractions (εj = ε/N) to make the relative density enhancement within and between dust phases more visible. We have also adjusted the colourbar in order to allow direct comparison with the settling tests performed by PL15 and Price et al. (2017). Note that the density enhancement due to settling has a shallower dynamic range than the built-in density gamut created by our grains-size distribution (see Fig. 2). As a further benchmarking exercise, we ran 10 single-phase simulations using the initial conditions from Table 1 and compared them to the results from the multiphase test above (see Fig. 4). Although the two scenarios are not strictly equivalent – the single-phase simulations do not include the backreaction from the N − 1 other phases – the global solutions still match because (i) the majority of the disc mass resides in the gas and (ii) the gas is essentially motionless throughout the simulation (see e.g. the top row in Fig. 3, which can be used as a proxy for the gas). In the limit of zero backreaction and a stationary gas phase, the system can be modelled analytically and numerically using a simplified set of fluid equations. Section A gives the full analysis. Figure 4. View largeDownload slide Comparison of dust fractions after 15 orbits when calculated by 10 single-phase simulations (black points) versus 1 multigrain simulation (coloured points). Not only does the multigrain method recover the correct solution, but the dispersion in εj is equal to or better than the single-phase simulations. Figure 4. View largeDownload slide Comparison of dust fractions after 15 orbits when calculated by 10 single-phase simulations (black points) versus 1 multigrain simulation (coloured points). Not only does the multigrain method recover the correct solution, but the dispersion in εj is equal to or better than the single-phase simulations. The large-scale agreement, we see in Fig. 4 does not extend to smaller scales. In Fig. 5, we zoom in on the $$s=0.1\, \mu \mathrm{m}$$ grains to illustrate the substructure in εj that develops as a result of the backreaction included from other dust phases. These differences between the single-phase and multigrain simulations continue to grow with time. Therefore, single-phase simulations should be used with caution in situations involving turbulent gas dynamics and/or long time-scales. Figure 5. View largeDownload slide A zoom in of the $$s=0.1\, \mu \mathrm{m}$$ grains in Fig. 4, highlighting the non-linear coupling between dust phases captured in a multigrain simulation (blue points) compared to the single-phase simulation (black points). The location of the peaks (resp. troughs) seen in ε1 correlate with the outer edges (resp. density peaks) of the other phases. We have added semi-transparent lines to help identify the location of the other phases in the figure. With the exception of the largest grain size, the remaining phases exhibit similar discrepancies with their single-phase counterpart. Figure 5. View largeDownload slide A zoom in of the $$s=0.1\, \mu \mathrm{m}$$ grains in Fig. 4, highlighting the non-linear coupling between dust phases captured in a multigrain simulation (blue points) compared to the single-phase simulation (black points). The location of the peaks (resp. troughs) seen in ε1 correlate with the outer edges (resp. density peaks) of the other phases. We have added semi-transparent lines to help identify the location of the other phases in the figure. With the exception of the largest grain size, the remaining phases exhibit similar discrepancies with their single-phase counterpart. The differences we observe in Fig. 5 are small, but prevent the test from being truly rigorous. One way of making the indirect coupling between dust phases vanishingly small is to concentrate all of the dust mass into the smallest grains that remain fixed to the gas, i.e. stationary. We found that by using a power-law index of p = 6.5, we could concentrate $${\sim } 99\, {\rm per \, cent}$$ of the dust mass in the two smallest phases (with j = 1 accounting for $${\sim } 92\, {\rm per \, cent}$$). Under these new conditions, our single-phase and multigrain simulations were a near perfect match at all scales. The only visible difference between the two scenarios was a minor reduction in dispersion in some of the multigrain phases, similar in magnitude to what is seen in Fig. 4. As an interesting aside, the steeper power-law index of p = 6.5 produces a 10 order-of-magnitude gap between the mid-plane densities of the largest and smallest dust grains. Happily, roundoff errors do not appear to corrupt the results in this situation, which we attribute to the fact that each εj is evolved separately. While the dust fractions are combined to calculate the gas properties, the gas–dust interaction depends only on the ratio of their masses. That is, the gas is not sensitive to tiny fluctuations in ε that may be introduced by loss of precision when combining εj of very different magnitudes. So far we have relied on comparing our multigrain results with single-phase simulations. In Fig. 6, we return to using the initial conditions from Table 1 and compare our multigrain solution to a grid-based numerical solution described in Section A2. The settling fronts in our SPH simulations match the simplified solutions to better than a few per cent for all except the largest grain sizes (see Table 2). As pointed out by PL15, the resolution follows the total mass rather than the dust mass, so it tends to over smooth the density peaks in the dust. Despite the smoothing, the locations of the settling fronts and the densities within the disc match very well. The L2 errors scale with the grain size (see column 2 in Table 2) and are a result of the oversmoothed dust peaks and the increased dispersion in the density at larger grain sizes. Figure 6. View largeDownload slide Dust densities for dust settling test in Fig. 4. The coloured points are from our multigrain simulation, while the black solution curves come from solving the equations (A3) and (A4) on a grid. The oversmoothing of the dust fronts is an artefact of tracking the total mass rather than just the dust mass. However, the location of the dust fronts and the enhancement of the mid-plane densities are well captured by the multigrain method. Figure 6. View largeDownload slide Dust densities for dust settling test in Fig. 4. The coloured points are from our multigrain simulation, while the black solution curves come from solving the equations (A3) and (A4) on a grid. The oversmoothing of the dust fronts is an artefact of tracking the total mass rather than just the dust mass. However, the location of the dust fronts and the enhancement of the mid-plane densities are well captured by the multigrain method. Table 2. L2 errors, computed by splash (Price 2007), between the multigrain dusty settling test from Section 4.3 and the analytic/numerical solutions from Section A. Column 2 is obtained using the numeric dust densities from Section A2 and column 3 using the analytic dust velocities from equation (A5). Large L2 errors are caused mainly by inconsistencies between the analytical and numerical models. That our density errors are lowest where our velocity errors are largest (and vice versa) indicates that our multigrain solution is valid. s [cm]  L2 density errors  L2 velocity errors  1.00 × 10−5  4.93 × 10−3  1.47  2.78 × 10−5  4.90 × 10−3  5.27 × 10−1  7.74 × 10−5  4.81 × 10−3  1.89 × 10−1  2.15 × 10−4  5.10 × 10−3  6.80 × 10−2  5.99 × 10−4  5.06 × 10−3  2.45 × 10−2  1.67 × 10−3  6.44 × 10−2  8.90 × 10−3  4.64 × 10−3  1.24 × 10−2  3.50 × 10−3  1.29 × 10−2  3.44 × 10−2  1.92 × 10−3  3.59 × 10−2  7.89 × 10−2  1.61 × 10−3  1.00 × 10−1  9.43 × 10−2  1.57 × 10−3  s [cm]  L2 density errors  L2 velocity errors  1.00 × 10−5  4.93 × 10−3  1.47  2.78 × 10−5  4.90 × 10−3  5.27 × 10−1  7.74 × 10−5  4.81 × 10−3  1.89 × 10−1  2.15 × 10−4  5.10 × 10−3  6.80 × 10−2  5.99 × 10−4  5.06 × 10−3  2.45 × 10−2  1.67 × 10−3  6.44 × 10−2  8.90 × 10−3  4.64 × 10−3  1.24 × 10−2  3.50 × 10−3  1.29 × 10−2  3.44 × 10−2  1.92 × 10−3  3.59 × 10−2  7.89 × 10−2  1.61 × 10−3  1.00 × 10−1  9.43 × 10−2  1.57 × 10−3  View Large In Fig. 7, we compare our multigrain simulation to the analytic solution in equation (A5). While we find L2 errors of order 0.1–1 per cent for grain sizes $$> {}10\, \mu$$m, there is a steady decline in accuracy as grain size decreases (see column 3 in Table 2). This progressive departure from the analytic model is a reflection of the fact that the gas is not completely stationary. Fluctuations in the gas velocity create a size-dependent velocity dispersion in the dust that primarily affects the smaller grain sizes. The larger dust grains, which are less susceptible to these fluctuations in the gas, exhibit less dispersion and better agreement with the analytic solution. Figure 7. View largeDownload slide Settling velocities for the dust settling test in Fig. 4. Again, the coloured points are results from our multigrain calculations, but the black curves are now the single-phase analytic solutions from equation (A5). Slight discrepancies are visible due to the coupling between dust phases and the fact that the gas is not stationary, both of which are ignored in the analytic solution. Figure 7. View largeDownload slide Settling velocities for the dust settling test in Fig. 4. Again, the coloured points are results from our multigrain calculations, but the black curves are now the single-phase analytic solutions from equation (A5). Slight discrepancies are visible due to the coupling between dust phases and the fact that the gas is not stationary, both of which are ignored in the analytic solution. 4.3.3 Performance We ran each of the test simulations above using OpenMP on eight cores from a single node. We found that our multigrain simulations with N = 10 dust phases were a factor of 2 slower than one single-phase simulation, thus making the multigrain simulations five times faster than their single-phase equivalent. This scaling improves as N increases, provided there is enough memory to handle the large array sizes. For example, we found the multigrain method to be ≈13 times faster when N = 100. We expect even better performance ratios relative to multifluid simulations, because multifluid methods require N times more simulation particles and often an added overhead for implicit time-stepping (explicit multifluid methods are impossibly slow for most of the grain sizes considered in this study; see PL15). Finally, because the multigrain method reuses the same simulation particles for all N dust phases, it requires less post-processing and, when N = 10, uses 55 per cent less disc space than an equivalent set of single-phase simulations (65 per cent less when N = 100). Files in which $$\Delta \boldsymbol {v}_{j}$$ is not written to disc1 are reduced by an additional 15 per cent. 4.4 Radial drift in a protoplanetary disc The dusty settling test in the previous section remains well approximated by single-phase methods. To demonstrate that our algorithm also works in regimes of strong backreaction, we computed the radial drift velocities for two dust phases in a protoplanetary disc with conditions such that the inward migration of the larger grains induces a discernable outward migration of the smaller grains. 4.4.1 Analytic solution An analytic solution for multiple dust phases migrating in an inviscid disc was derived by Bai & Stone (2010b). Neglecting vertical gravity, they show that the hydrostatic equilibrium equations can be written in block matrix form as follows:   \begin{eqnarray} {\left(\begin{array}{cc} {\bf I} + \boldsymbol{\Gamma } & -2 \boldsymbol{\Lambda }\\ \boldsymbol {\Lambda }/2 & {\bf I} + \boldsymbol {\Gamma } \end{array}\right)} {\left(\begin{array}{cc}\boldsymbol {\mathcal {V}}_r\\ \boldsymbol {\mathcal {V}}_\phi \end{array}\right)} = - \eta v_K {\left(\begin{array}{cc}\boldsymbol {0}\\ \boldsymbol {1} \end{array}\right)}, \end{eqnarray} (64)where $${\bf I}$$ is the identity matrix, $$\boldsymbol {\mathcal {V}}_r \equiv \left( v_{1r},v_{2r}, \ldots ,v_{nr} \right)^\intercal$$, and $$\boldsymbol {\mathcal {V}}_\phi \equiv \left( v_{1\phi },v_{2\phi }, \ldots ,v_{n\phi } \right)^\intercal$$ are the radial and azimuthal velocities for each dust phase, respectively. The matrix $$\boldsymbol {\Lambda } \equiv {\rm diag}\left\lbrace {\rm St}_1, {\rm St}_2, \ldots , {\rm St}_n \right\rbrace$$ is a diagonal matrix of the Stokes numbers for uncoupled dust phases (i.e. $$\mathrm{St}= t_{\mathrm{s}}^{N=1}\Omega _\mathrm{K}$$), while $$\boldsymbol {\Gamma } \equiv \left(\boldsymbol {\mathcal {E}}, \boldsymbol {\mathcal {E}}, \ldots , \boldsymbol {\mathcal {E}} \right)^\intercal$$ is a matrix made up of the dust-to-gas ratios, where $$\boldsymbol {\mathcal {E}} \equiv \left(\mathcal {E}_1, \mathcal {E}_2, \ldots , \mathcal {E}_n \right)^\intercal$$ and $$\mathcal {E}_j \equiv \rho _{\mathrm{d}j}/\rho _{\mathrm{g}}= \epsilon _{j}/(1-\epsilon )$$. Bai & Stone (2010b) provide a closed-form solution to equation (64); however, we found it more convenient to solve it numerically. 4.4.2 Set-up We set-up a 3D, locally isothermal gas disc using the following power-law parametrizations (see e.g. Laibe, Gonzalez & Maddison 2012)   \begin{eqnarray} c_{{\rm s}}(r) = c_{{\rm s},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{-q/2}, \end{eqnarray} (65)  \begin{eqnarray} H_{\rm g}(r) = H_{{\rm g},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{3/2- q/2}, \end{eqnarray} (66)  \begin{eqnarray} \Sigma _{\rm g}(r) = \Sigma _{{\rm g},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{-p}, \end{eqnarray} (67)  \begin{eqnarray} \rho _{\rm g}(r,z) = \displaystyle\frac{\Sigma _{\rm g}}{\sqrt{2\pi }H} \exp {\left[ -\displaystyle\frac{z^2}{2 H^2}\right]}, \end{eqnarray} (68)where Σg is the local surface density for the gas, quantities with the subscript ‘1 au’ are reference values measured at r = 1 au, and the parameters p = 1 and q = 0.5 are power-law exponents controlling the density and temperature (i.e. flaring) of the disc, respectively. We set the radial velocity to zero and correct the orbital velocities from pure Keplerian rotation to account for the radial pressure gradient in the disc,   \begin{eqnarray} v_\phi = v_K(1 - \eta ), \end{eqnarray} (69)where the pressure gradient parameter η is given by   \begin{eqnarray} \eta \approx \frac{1}{4} \left(\frac{H_{\rm g}}{r} \right)^2 \left[ 3 + 2 p+ q - \left( 3 - q \right) \left( \frac{z}{H_{\rm g}} \right)^2 \right], \end{eqnarray} (70)to order z2/r2 (see e.g. Takeuchi & Lin 2002). We add the dust by assigning a uniform dust fraction to all of the particles. Because the diffusion approximation assumes the stopping time is much shorter than the dynamical time-scale, we do not give the dust a separate azimuthal velocity. The analytic solution from Bai & Stone (2010b) is 2D and assumes that dust resides in the mid-plane of the disc. As both gas density and gravity decrease with increasing z, we expect the dust at high altitudes to migrate slower than dust in the mid-plane. To compensate, we only compare migration rates for |z| < Hg(r), we bin the particles radially into 50 logarithmically spaced bins (using the same binning method described for the grain-size distribution previously), and we average the radial velocities both azimuthally and vertically within each bin. One final caveat remains: the steady-state analytic model assumes the disc is inviscid, whereas SPH disc simulations are inherently viscous. Normally, we would relax our disc into a quasi-steady state and use our instantaneous gas and dust mid-plane densities as initial conditions for the model – thereby allowing us to account for any non-steady-state processes like settling and/or migration. However, the lack of viscidity in the model produces rigid assumptions about the gas velocity that are not met in our viscous SPH simulations. As a result, we find that our simulation relaxes into a steady state that is substantially different to the analytic solution. To our knowledge, there is currently no analytic solution for radial velocities in viscous discs. Deriving such a solution goes beyond the scope of this paper; therefore, we will revisit the problem in a future study. In the meantime, we can circumvent this incompatibility in the present study by using the initial state of the system, where we have full control over the velocities and we can mimic the conditions of an inviscid disc. Testing the initial conditions, albeit unorthodox, still yields valuable information about our method for two reasons. First, the terminal velocity approximation breaks down when the time-step is smaller than a few stopping times. Because the typical time for drift to relax is on the order of a few stopping times, we do not need to wait for the dust velocities to equilibrate. In other words, the full asymptotic radial velocities are obtained after the very first loop over the particles (what we call t = 0) when P, Tsj, Ts, $$\Delta \boldsymbol {v}_{j}$$, and $$\Delta \boldsymbol {v}$$ are all calculated – the quantities we use to construct the velocity profiles of the gas and dust. Secondly, the individual ga s and dust properties are calculated (as opposed to being evolved). Therefore, the test is more sensitive to how we calculate the forces than how we evolve the mixture. Since our force prescription does not vary with time, the test is almost as useful at t = 0, as it would be once the system has reached a steady state. 4.4.3 Results Using the same grain-size distribution from Section 4.3.1, we set-up a dusty protoplanetary disc around a solar mass star with an inner and outer radius of rin = 1 au and rout = 300 au, respectively. The gas disc has the following reference values: $$c_{{\rm s},1{\rm au}} \approx 1.5 \, {\rm km\,s}^{-1}$$, Hg, 1au = 0.05 au, and $$\Sigma _{{\rm g},1{\rm au}} \approx 166 \, {\rm g\,cm}^{-2}$$. The dust disc for each grain size is set equal to the gas in size and shape, but scaled in mass by the dust fractions listed in Table 3, such that, when all of the dust phases are included, the total dust mass comprises one-third of the total mass of the system (i.e. a dust-to-gas ratio of $$\mathcal {E} = 0.5$$). For the 10 single-phase calculations, it is not possible to simultaneously match the dust fractions, dust-to-gas ratios, and the gas/dust densities of the multigrain case. Therefore, we chose to keep the respective dust fraction and the total surface density of the disc the same, while allowing the dust-to-gas ratio for each dust phase to change as needed. This discrepancy with the multigrain calculations results in different surface densities for the gas and dust, but the effects are unimportant in this context since outward migration of dust in a single-phase simulation can only be achieved by drastically changing the structure of the disc (e.g. with a radially increasing pressure profile). Table 3. The initial grain sizes (s), dust fractions (ε), and dust-to-gas ratios ($$\mathcal {E}$$) used in the outward migration test in Fig. 8. Variables with/without the subscript j indicate multigrain/single-phase values, respectively. The final row gives the sum of the individual dust fractions and dust-to-gas ratios, highlighting the inherent discrepancies between setups of single- and multi-phased simulations. s and sj [cm]  ε and εj  $$\mathcal {E}$$  $$\mathcal {E}_j$$  1.00 × 10−5  1.34 × 10−3  1.34 × 10−3  2.01 × 10−3  2.78 × 10−5  2.24 × 10−3  2.24 × 10−3  3.36 × 10−3  7.74 × 10−5  3.74 × 10−3  3.75 × 10−3  5.60 × 10−3  2.15 × 10−4  6.23 × 10−3  6.27 × 10−3  9.35 × 10−3  5.99 × 10−4  1.04 × 10−2  1.05 × 10−2  1.56 × 10−2  1.67 × 10−3  1.73 × 10−2  1.76 × 10−2  2.60 × 10−2  4.64 × 10−3  2.89 × 10−2  2.98 × 10−2  4.34 × 10−2  1.29 × 10−2  4.83 × 10−2  5.07 × 10−2  7.24 × 10−2  3.59 × 10−2  8.05 × 10−2  8.76 × 10−2  1.21 × 10−1  1.00 × 10−1  1.34 × 10−1  1.55 × 10−1  2.01 × 10−1  Sum total:  1/3  3.65 × 10−1  1/2  s and sj [cm]  ε and εj  $$\mathcal {E}$$  $$\mathcal {E}_j$$  1.00 × 10−5  1.34 × 10−3  1.34 × 10−3  2.01 × 10−3  2.78 × 10−5  2.24 × 10−3  2.24 × 10−3  3.36 × 10−3  7.74 × 10−5  3.74 × 10−3  3.75 × 10−3  5.60 × 10−3  2.15 × 10−4  6.23 × 10−3  6.27 × 10−3  9.35 × 10−3  5.99 × 10−4  1.04 × 10−2  1.05 × 10−2  1.56 × 10−2  1.67 × 10−3  1.73 × 10−2  1.76 × 10−2  2.60 × 10−2  4.64 × 10−3  2.89 × 10−2  2.98 × 10−2  4.34 × 10−2  1.29 × 10−2  4.83 × 10−2  5.07 × 10−2  7.24 × 10−2  3.59 × 10−2  8.05 × 10−2  8.76 × 10−2  1.21 × 10−1  1.00 × 10−1  1.34 × 10−1  1.55 × 10−1  2.01 × 10−1  Sum total:  1/3  3.65 × 10−1  1/2  View Large The left-hand and right-hand panels in Fig. 8 show the mean radial dust velocities for the individual and combined cases, respectively. Coloured points are the velocities calculated from the SPH mixture, while the solid black lines show the corresponding analytic solutions. Although the two largest grain sizes exhibit only minor changes to their velocities after the addition of the other dust species, the eight smaller sizes experience a complete reversal in migration direction. This change in sign is caused by the exchange of angular momentum as the larger grains drag the sub-Keplerian gas into faster orbits, thereby pushing the gas radially outwards. The smaller dust grains, who are more sensitive to changes in the gas flow, are then carried outwards along with the gas. Outward migration of dust in a disc with a radially decreasing pressure gradient cannot be replicated with only one dust phase; conservation of angular momentum requires one or more phases to radially contract as the others expand. Importantly, the multigrain formalism correctly resolves the velocities for both outward and inward migrating species. Figure 8. View largeDownload slide Mean radial drift velocities near the mid-plane of a protoplanetary disc for 10 dust phases of different grain size, individually coupled with the gas (left) and simultaneously coupled (right). Coloured points (connected by either dashed or dotted lines for clarity) show the mean radial dust velocity calculated for each dust phase from the barycentric simulation data and solid black lines indicate the corresponding analytic solutions. The complete lack of outward migration in the left-hand panel accentuates how single-phase simulations can miss dynamical effects caused by the presence of other dust phases. The right-hand panel illustrates that the multigrain formalism is capable of resolving outward/inward migration velocities of different dust species within a single set of simulation particles. Figure 8. View largeDownload slide Mean radial drift velocities near the mid-plane of a protoplanetary disc for 10 dust phases of different grain size, individually coupled with the gas (left) and simultaneously coupled (right). Coloured points (connected by either dashed or dotted lines for clarity) show the mean radial dust velocity calculated for each dust phase from the barycentric simulation data and solid black lines indicate the corresponding analytic solutions. The complete lack of outward migration in the left-hand panel accentuates how single-phase simulations can miss dynamical effects caused by the presence of other dust phases. The right-hand panel illustrates that the multigrain formalism is capable of resolving outward/inward migration velocities of different dust species within a single set of simulation particles. The relative angular velocity between the gas and dust varies with height. In fact, η changes sign at z ≈ 1.5Hg, meaning that dust particles rotate slower than the gas above this height. Because phantom is a 3D code, our calculations systematically underestimate the 2D analytic solution, which assumes all of the dust is rotating in the mid-plane of the disc. We can reduce this offset by only considering particles near the mid-plane, but having fewer particles to average can make the data more noisy. In making Fig. 8, we used 107 particles in the disc, and discard all particles with z > Hg (∼1/3 of the particles). Even with so many particles remaining, the inner ∼20 bins are very noisy (note the first 14 are not shown), with values ranging between −1.5 and 0.3 η$$v$$K. Also note that the standard deviation for most bins is larger than the size of the plotting window, with typical magnitudes ranging from tens to hundreds [ηvK]. Thus, we should not take the discrepancy between the numerical and analytic solutions too seriously. 5 DISCUSSION AND CONCLUSIONS We have derived and implemented a numerical scheme using SPH that is capable of simulating multiple dust phases composed of small dust grains coupled to the gas in the terminal velocity approximation (i.e. when the stopping time is short compared to the computational time-step). Our method simulates dust using a dimensionless dust fraction, as opposed to traditional methods that employ additional sets of simulation particles. By expanding the scalar dust fraction into an array of N dust fractions that are independently evolved and coupled to the gas, we obtain a method that scales better in terms of computational time and resources as N becomes large. Another benefit of evolving the mixture is that the multigrain method circumvents having to resolve the prohibitive temporal and spatial resolution criteria for small dust grains that usually choke multifluid simulations with separate gas and dust particles. We have demonstrated that the multigrain continuum and discretized equations correctly reduce to the equations described by PL15 when N = 1 and that there is no loss in accuracy when simulating a single phase using our multigrain framework – even when that dust phase is divided into multiple mass bins. On the other hand, when simulating multiple unique dust phases, the multigrain method is superior to using multiple single-phase simulations, not only in terms of computational speed and efficiency as discussed above, but also in terms of accuracy as a result of capturing the indirect coupling (via the gas) between dust phases. Although the deviations between our multigrain and single-phase simulations were small for the select test cases we performed in Section 4.3 (∼ few per cent), there are a few additional points to consider for general applications: (i) perturbations from other dust phases accumulate over time, (ii) perturbations from concentrated dust grains (or equivalently, higher dust-to-gas ratios) are larger in magnitude than for dispersed grains, and (iii) perturbations between phases can further be accentuated by motion of the gas (as opposed to the stationary gas phase in our settling tests). In light of these concerns, we caution against using single-phase simulations where possible and encourage the adoption of the more accurate and efficient multigrain method we present here. Finally, the present multigrain algorithm can only be used for small dust grains within the terminal velocity approximation, which is accurate only when the stopping time is shorter than the dynamical time-scale (LP14a). To extend to larger grains, we would need to either implement the full multiphase one-fluid equations with implicit time-stepping from LP14c or develop a hybrid between the one- and multi-fluid methods. Both have advantages and disadvantages, but are beyond the scope of this paper. Presently, we do not account for the evolution in grain size through growth and fragmentation. However, incorporating grain size evolution into the multigrain framework would be straightforward because the mass and number of the simulation particles does not have to change with time. Acknowledgements We thank the anonymous referee whose careful review helped improve this paper significantly. We would also like to acknowledge Matthew Bate, Giovanni Dipierro, and Christophe Pinte for useful discussions. D.J.P. is grateful for funding via an Australian Research Council (ARC) Future Fellowship, FT130100034. M.A.H. and D.J.P. were funded by ARC Discovery Project grant DP130102078. G.L. acknowledges financial support from PNP, PNPS, PCMI of CNRS/INSU, CEA, and CNES, France. This work has in part been carried out within the framework of the National Centre for Competence in Research PlanetS, supported by the Swiss National Science Foundation. Finally, computations have been done on both the SwinSTAR supercomputer hosted at Swinburne University of Technology and the Piz Daint supercomputer hosted at the Swiss National Computational Centre. Footnotes 1 In the diffusion approximation, $$\Delta \boldsymbol {v}_{j}$$ is a calculated quantity needed for recovering the gas and dust velocities during post-processing. However, as $$\Delta \boldsymbol {v}_{j}$$ is not required in any of the evolution equations, we often omit writing it to disc in order to save space. REFERENCES Bai X.-N., Stone J. M., 2010a, ApJS , 190, 297 https://doi.org/10.1088/0067-0049/190/2/297 CrossRef Search ADS   Bai X.-N., Stone J. M., 2010b, ApJ , 722, 1437 https://doi.org/10.1088/0004-637X/722/2/1437 CrossRef Search ADS   Balsara D. S., Tilley D. A., Rettig T., Brittain S. D., 2009, MNRAS , 397, 24 https://doi.org/10.1111/j.1365-2966.2009.14606.x CrossRef Search ADS   Barranco J. A., 2009, ApJ , 691, 907 https://doi.org/10.1088/0004-637X/691/2/907 CrossRef Search ADS   Barrière-Fouchet L., Gonzalez J.-F., Murray J. R., Humble R. J., Maddison S. T., 2005, A&A , 443, 185 CrossRef Search ADS   Chiang E., 2008, ApJ , 675, 1549 https://doi.org/10.1086/527354 CrossRef Search ADS   Cullen L., Dehnen W., 2010, MNRAS , 408, 669 https://doi.org/10.1111/j.1365-2966.2010.17158.x CrossRef Search ADS   Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS , 453, L73 https://doi.org/10.1093/mnrasl/slv105 CrossRef Search ADS   Epstein P. S., 1924, Phys. Rev. , 23, 710 https://doi.org/10.1103/PhysRev.23.710 CrossRef Search ADS   Garaud P., Lin D. N. C., 2004, ApJ , 608, 1050 https://doi.org/10.1086/420839 CrossRef Search ADS   Haworth T. J. et al.  , 2016, PASA , 33, e053 https://doi.org/10.1017/pasa.2016.45 Hutchison M. A., Price D. J., Laibe G., Maddison S. T., 2016, MNRAS , 461, 742 https://doi.org/10.1093/mnras/stw1126 CrossRef Search ADS   Jacquet E., Balbus S., Latter H., 2011, MNRAS , 415, 3591 https://doi.org/10.1111/j.1365-2966.2011.18971.x CrossRef Search ADS   Johansen A., Klahr H., 2005, ApJ , 634, 1353 https://doi.org/10.1086/497118 CrossRef Search ADS   Laibe G., Price D. J., 2012a, MNRAS , 420, 2345 https://doi.org/10.1111/j.1365-2966.2011.20202.x CrossRef Search ADS   Laibe G., Price D. J., 2012b, MNRAS , 420, 2365 https://doi.org/10.1111/j.1365-2966.2011.20201.x CrossRef Search ADS   Laibe G., Price D. J., 2014a, MNRAS , 440, 2136 (LP14a) https://doi.org/10.1093/mnras/stu355 (LP2014a) CrossRef Search ADS   Laibe G., Price D. J., 2014b, MNRAS , 440, 2147 (LP14b) https://doi.org/10.1093/mnras/stu359 (LP2014b) CrossRef Search ADS   Laibe G., Price D. J., 2014c, MNRAS , 444, 1940 (LP14c) https://doi.org/10.1093/mnras/stu1367 (LP2014c) CrossRef Search ADS   Laibe G., Gonzalez J.-F., Maddison S. T., 2012, A&A , 537, A61 CrossRef Search ADS   Lee A. T., Chiang E., Asay-Davis X., Barranco J., 2010, ApJ , 718, 1367 https://doi.org/10.1088/0004-637X/718/2/1367 CrossRef Search ADS   Lorén-Aguilar P., Bate M. R., 2014, MNRAS , 443, 927 https://doi.org/10.1093/mnras/stu1173 CrossRef Search ADS   Lorén-Aguilar P., Bate M. R., 2015, MNRAS , 454, 4114 https://doi.org/10.1093/mnras/stv2262 CrossRef Search ADS   Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ , 217, 425 https://doi.org/10.1086/155591 CrossRef Search ADS   Miniati F., 2010, J. Comput. Phys. , 229, 3916 https://doi.org/10.1016/j.jcp.2010.01.034 CrossRef Search ADS   Monaghan J. J., 1997, J. Comput. Phys. , 138, 801 https://doi.org/10.1006/jcph.1997.5846 CrossRef Search ADS   Monaghan J. J., Kocharyan A., 1995, Comput. Phys. Commun. , 87, 225 https://doi.org/10.1016/0010-4655(94)00174-Z CrossRef Search ADS   Morris J. P., Monaghan J. J., 1997, J. Comput. Phys. , 136, 41 https://doi.org/10.1006/jcph.1997.5690 CrossRef Search ADS   Paardekooper S.-J., Mellema G., 2004, A&A , 425, L9 CrossRef Search ADS   Price D. J., 2004, PhD thesis , Univ. Cambridge, Cambridge Price D. J., 2007, PASA , 24, 159 https://doi.org/10.1071/AS07022 CrossRef Search ADS   Price D. J., 2008, J. Comput. Phys. , 227, 10040 https://doi.org/10.1016/j.jcp.2008.08.011 CrossRef Search ADS   Price D. J., Federrath C., 2010, MNRAS , 406, 1659 Price D. J., Laibe G., 2015, MNRAS , 451, 813 (PL15) https://doi.org/10.1093/mnras/stv996 (PL15) CrossRef Search ADS   Price D. J., Monaghan J. J., 2004, MNRAS , 348, 139 https://doi.org/10.1111/j.1365-2966.2004.07346.x CrossRef Search ADS   Price D. J., Monaghan J. J., 2007, MNRAS , 374, 1347 https://doi.org/10.1111/j.1365-2966.2006.11241.x CrossRef Search ADS   Price D. J. et al.  , 2017, PASA, preprint (arXiv:1702.03930) Saffman P. G., 1962, J. Fluid Mech. , 13, 120 https://doi.org/10.1017/S0022112062000555 CrossRef Search ADS   Takeuchi T., Lin D. N. C., 2002, ApJ , 581, 1344 https://doi.org/10.1086/344437 CrossRef Search ADS   Tricco T. S., Price D. J., Laibe G., 2017, MNRAS , 471, L52 https://doi.org/10.1093/mnrasl/slx096 CrossRef Search ADS   Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS , 387, 427 https://doi.org/10.1111/j.1365-2966.2008.13260.x CrossRef Search ADS   Yang C.-C., Johansen A., 2016, ApJS , 224, 39 https://doi.org/10.3847/0067-0049/224/2/39 CrossRef Search ADS   Youdin A. N., Goodman J., 2005, ApJ , 620, 459 https://doi.org/10.1086/426895 CrossRef Search ADS   Youdin A., Johansen A., 2007, ApJ , 662, 613 https://doi.org/10.1086/516729 CrossRef Search ADS   APPENDIX A: SOLUTIONS TO THE SETTLING TEST In the limit of very small dust-to-gas ratios, we can neglect the backreaction of the dust on the gas. The dust can then be treated as N independent phases, moving inside a static gas background, and governed by the following one-dimensional equations   \begin{eqnarray} \frac{\mathrm{\partial} \rho _{\mathrm{d}}}{\mathrm{\partial} t} + v_{\mathrm{d}}\frac{\mathrm{\partial} \rho _{\mathrm{d}}}{\mathrm{\partial} z} = - \rho _{\mathrm{d}}\frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} z}, \end{eqnarray} (A1)  \begin{eqnarray} \frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} t} + v_{\mathrm{d}}\frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} z} = - \frac{v_{\mathrm{d}}}{t_{\mathrm{s}}^{N=1}} - \frac{\mathcal {G} M z}{\left( r^2 + z^2 \right)^{3/2}}, \end{eqnarray} (A2)where we have dropped the subscript j to emphasize that the phases are no longer coupled. To aid our analysis, we define the dimensionless variables $$\bar{\rho } \equiv \rho /\rho _{\mathrm{g},0}$$, $$\bar{v} \equiv v/v_\mathrm{K}$$, $$\bar{z} \equiv z/r$$, and $$\bar{t} \equiv t \, \Omega _\mathrm{K}$$, where $$v_\mathrm{K}= \sqrt{\mathcal {G} M/r}$$ and ΩK = $$v$$K/r are the Keplerian velocity and frequency, respectively. Substituting these quantities into equations (A1) and (A2), we obtain the corresponding non-dimensionalized equations in the form   \begin{eqnarray} \frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{v}_\mathrm{d}\frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{z}} = - \bar{\rho }_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}}, \end{eqnarray} (A3)  \begin{eqnarray} \frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{v}_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}} = - \frac{\bar{v}_\mathrm{d}}{\mathrm{St}} - \frac{\bar{z}}{\left(1 + \bar{z}^2 \right)^{3/2}}, \end{eqnarray} (A4)where $$\mathrm{St}= t_{\mathrm{s}}^{N=1}\Omega _\mathrm{K}$$ is the Stokes number. A1 Analytic solution to settling problem Importantly, equation (A4) is independent of $$\bar{\rho }_\mathrm{d}$$. As a first-order partial differential equation, a solution for $$\bar{v}_\mathrm{d}$$ could potentially be obtained via the method of characteristics. We took a simpler approach by solving the Lagrangian form of the equation with a convective derivative. In this form, the equation is a first-order ordinary differential equation that can be solved using an integrating factor. Assuming the dust starts from rest, the solution for the dust velocity is   \begin{eqnarray} \bar{v}_\mathrm{d}= \frac{\mathrm{St}\, \bar{z}}{\left(1+\bar{z}^2\right)^{3/2}} \left( \mathrm{e}^{-\bar{t}/\mathrm{St}} - 1\right). \end{eqnarray} (A5) The same procedure can be used to obtain an equation for the dust density. However, the resulting solution does not conserve mass, since it assumes that an infinitely extended dust distribution continuously rains down on to the disc. The density solution nevertheless correctly predicts the location of the incoming dust front and also the interior density profile so long as t is less than the settling time-scale. This is not a problem in the velocity solution because all of our dust grains settle at their terminal velocity within the disc, effectively erasing any built-up momentum gained at higher altitudes. A2 1D numerical solution to the settling problem To remove the assumptions imposed by the analytic solution, we also compared our multigrain results with a numerical solution to equations (A3) and (A4). We solved the equations on a one-dimensional grid using an implicit Crank–Nicolson algorithm, using forward differences in time and centred differences in space. Because the temporal derivative is centred half a time-step in the future, we replace all of the other terms with time averages centred about the same time. Grouping terms based on their location in time, we can then write each equation as a linear system in the form   \begin{eqnarray} \bf {A} \boldsymbol {x}^{n+1} = \bf {B} \boldsymbol {x}^n + \boldsymbol {C}, \end{eqnarray} (A6)where superscripts designate the time level, $$\bf {A}$$ and $$\bf {B}$$ are square sparse matrices, $$\boldsymbol {x}$$ is the fluid variable for which we are solving, and the column vector $$\boldsymbol {C}$$ is a placeholder for all terms independent of $$\boldsymbol {x}$$ ($$\boldsymbol {C}=0$$ when $$\boldsymbol {x}$$ represents density). Once the boundary conditions have been accounted for in $$\bf {A}$$ and $$\bf {B}$$, the solution at time n + 1 can be obtained symbolically via   \begin{eqnarray} \boldsymbol {x}^{n+1} = \bf {A}^{-1} \left( \bf {B} \boldsymbol {x}^n + \boldsymbol {C} \right). \end{eqnarray} (A7)The non-linearity in the advection term in equation (A4) keeps us from obtaining the solution using the exact method as outlined above because it would irreversibly mix terms from different time-steps. To overcome this problem, we assume the leading $$\bar{v}_\mathrm{d}$$ in the advection term is known and designate it as $$\tilde{v}_\mathrm{d}$$ to keep it separate from the other velocity terms. We account for $$\tilde{v}_\mathrm{d}$$ and the fact that the fluid equations are coupled by using a predictor-corrector scheme to advance the system forward in time. Designating predicted quantities with asterisks, we advance the system in four steps: $$\bar{\rho }_\mathrm{d}^*$$ is predicted assuming $$\bar{v}_\mathrm{d}$$ is constant (i.e. $$\bar{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^n$$). $$\bar{v}_\mathrm{d}^*$$ is predicted assuming that $$\tilde{v}_\mathrm{d}$$ is constant and equal to $$\bar{v}_\mathrm{d}^n$$. $$\bar{\rho }_\mathrm{d}^{n+1}$$ is corrected assuming $$\bar{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^*$$. $$\bar{v}_\mathrm{d}^{n+1}$$ is corrected assuming $$\tilde{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^*$$ and $$\tilde{v}_\mathrm{d}^n = \bar{v}_\mathrm{d}^n$$. Using the same physical parameters as in Section 4.3, we discretize the region z ∈ [−3H, 3H] with 1002 cell-centred grid points, including ghost points. The boundary condition for the velocity is $$v$$d( ± 3H, t) = 0, which consequently enforces the following boundary condition for the density:   \begin{eqnarray} \frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{\rho }_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}} = 0 \end{eqnarray} (A8)at the same locations. The initial conditions are $$v$$d(z, 0) = 0 and ρd(z, 0) = εjρg. We found that a dimensionless time-step of 1 was sufficient to keep the algorithm stable. As the dust settles, we do get some low-density numerical noise in the wings of the disc, but this noise is always separated from the settling dust layer by a region of zero density. We have verified that our results do not change when we force the density to zero beyond the first encountered zero-density grid point on either side of the mid-plane. The close match between our multigrain results and the analytic and numerical solutions demonstrates that our method works. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# MULTIGRAIN: a smoothed particle hydrodynamic algorithm for multiple small dust grains and gas

, Volume 476 (2) – May 1, 2018
13 pages

/lp/ou_press/multigrain-a-smoothed-particle-hydrodynamic-algorithm-for-multiple-FPePjnXLey
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty367
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present a new algorithm, multigrain, for modelling the dynamics of an entire population of small dust grains immersed in gas, typical of conditions that are found in molecular clouds and protoplanetary discs. The multigrain method is more accurate than single-phase simulations because the gas experiences a backreaction from each dust phase and communicates this change to the other phases, thereby indirectly coupling the dust phases together. The multigrain method is fast, explicit and low storage, requiring only an array of dust fractions and their derivatives defined for each resolution element. hydrodynamics, methods: numerical, protoplanetary discs, dust, extinction, ISM: kinematics and dynamics 1 INTRODUCTION Modelling the interaction of multiple dust grains simultaneously with the gas is a ‘grand challenge’ in protoplanetary disc modelling (Haworth et al. 2016), since discs involve grains with sizes spanning several orders of magnitude, from sub-micron grains to km-sized planetesimals. Grains of different sizes experience different dynamics since small grains are lighter and more easily influenced by the gas compared to larger, heavier grains. The usual approach to dusty gas dynamics is to model the gas and dust as separate fluids. The gas is modelled either on a grid (Paardekooper & Mellema 2004; Youdin & Johansen 2007; Balsara et al. 2009; Bai & Stone 2010a; Miniati 2010; Yang & Johansen 2016) or on a set of Lagrangian particles (Monaghan & Kocharyan 1995; Monaghan 1997; Barrière-Fouchet et al. 2005; Laibe & Price 2012a,b; Lorén-Aguilar & Bate 2014, 2015); similarly for each dust phase (although the discretization method often defaults to the one used by the gas). During simulation, the gas and dust fluids are evolved independently, but interact via a common drag force (e.g. Saffman 1962; Garaud & Lin 2004). Although grid- and particle-based methods each have their own distinct advantages/disadvantages (e.g. Price & Federrath 2010), they both require prohibitively small time-steps or implicit methods at high drag. Furthermore, Laibe & Price (2012a,b) discovered a drag resolution criterion that becomes increasingly restrictive with smaller grain sizes and applies generally to any method that models dust on a grid or on a set of particles that is not collocated with the gas at all times. While Laibe & Price (2012a,b), and later Lorén-Aguilar & Bate (2014), tested this spatial criterion using smoothed particle hydrodynamics (SPH), Youdin & Johansen (2007) inferred a similarly high-resolution requirement in hybrid grid-particle simulations. Failing to meet this criterion may explain the first-order convergence rate in high drag regimes observed by Miniati (2010), Bai & Stone (2010a), and Yang & Johansen (2016). To address the restrictive temporal and spatial restrictions that exist for high drag regimes, Laibe & Price (2014a,b,c, hereafter LP14a; LP14b; LP14c) and Price & Laibe (2015, hereafter PL15) developed a single-fluid formulation appropriate for small grains – similar to earlier formulations by Johansen & Klahr (2005). The dust–gas mixture is advected at the barycentric velocity and whose density is equal to the total density of the mixture. In the context of SPH, this means the mixture is represented by a single set of SPH particles with an evolution equation for the dust fraction (LP14b; PL15). While the above methods provide a means of modelling discs or molecular clouds with a single embedded dust phase, the challenge is to span the observed range of grain sizes. The typical approach is the one we recently used in Dipierro et al. (2015), where a series of single-phase simulations were stitched together in post-processing to interpret the dark structures observed at millimetre wavelengths by the ALMA interferometer in the disc surrounding the star HL Tau. In that paper, the method from Laibe & Price (2012a) was used to model the dynamics of mm-sized grains and larger, while the smaller grains were modelled using the method from PL15. Besides being tedious, the procedure used by Dipierro et al. (2015) is slow and, more importantly, neglects the indirect coupling between dust phases caused by the ‘backreaction’ of individual phases on the gas, which in turn influences the grain dynamics. Neglecting this backreaction misses important effects such as outward migration of dust particles (Bai & Stone 2010a) and/or modification of the linear growth rate of the streaming instability LP14c. Also, backreaction in individual grain size simulations is both annoying and wrong – annoying because the different response of the gas makes stacking of different dust grain distributions difficult (Tricco, Price & Laibe 2017); wrong because the gas should respond to the entire dust mixture, rather than each grain size individually. In this paper, we develop a new multigrain algorithm for modelling the dynamics of multiple dust phases, based on the analytical work presented in LP14c. Because much of the opacity and accompanying scattering/emission in astrophysical environments stems from the presence of dust grains that can be considered ‘small’ (i.e. where the terminal velocity approximation is valid), we focus on deriving and implementing the SPH versions of the continuum equations for the multiphase, terminal velocity approximation – generalizing the single-dust-phase method developed in PL15. 2 THE DIFFUSION APPROXIMATION FOR MULTIPLE DUST SPECIES 2.1 Continuum equations We consider a system consisting of a mixture of a single gas phase and N strongly coupled dust phases. Throughout this paper, we use the indices a, b, and c to refer to individual simulation particles that move at the barycentric velocity of the mixture. Subscript or superscripts g and d are used for gas and dust properties, respectively. Finally, we identify the fluid quantities for each of the N different dust phases using the index j. 2.1.1 General equations LP14c derived the general continuum fluid equations for a mixture of gas and N coupled dust species moving in a barycentric reference frame. They further showed that in strongly coupled regimes – i.e. first order in tj/T, where tj is a drag time-scale specific to each grain type (see equation 16; note the difference in notation from that of LP14c) and T is the time-scale for a sound wave to propagate over a typical distance L (commonly referred to as the terminal velocity approximation; see e.g. Youdin & Goodman 2005; Chiang 2008; Barranco 2009; Lee et al. 2010; Jacquet, Balbus & Latter 2011) – the fluid equations reduce to   \begin{eqnarray} \frac{{\rm d} \rho }{{\rm d} t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (1)  \begin{eqnarray} \frac{{\rm d} \epsilon _{j}}{{\rm d} t} = - \frac{1}{\rho } \nabla \cdot \left[ \rho \epsilon _{j}\left( \Delta \boldsymbol {v}_{j}- \epsilon \Delta \boldsymbol {v}\right) \right], \end{eqnarray} (2)  \begin{eqnarray} \frac{{\rm d} \boldsymbol {v}}{{\rm d} t} = \left( 1 - \epsilon \right) \boldsymbol {f}_{\mathrm{g}}+ \sum _{j}\epsilon _{j}\boldsymbol {f}_{\mathrm{d}j}+ \boldsymbol {f}, \end{eqnarray} (3)  \begin{eqnarray} \frac{{\rm d} u}{{\rm d} t} = - \frac{P}{\rho _{\mathrm{g}}} \nabla \cdot \boldsymbol {v}+ \epsilon \Delta \boldsymbol {v}\cdot \nabla u, \end{eqnarray} (4)  \begin{eqnarray} \Delta \boldsymbol {v}_{j} = \big[\Delta \boldsymbol {f}_{j}- \sum _{k}\epsilon _{k}\Delta \boldsymbol {f}_{k}\big] \epsilon _{j}t_{j}, \end{eqnarray} (5)where d/dt is the convective derivative using the barycentric velocity $$\boldsymbol {v}$$,   \begin{eqnarray} \boldsymbol {v}\equiv \displaystyle \frac{\rho _{\mathrm{g}}\boldsymbol {v}_{\mathrm{g}}+ \displaystyle \sum _{j}\rho _{\mathrm{d}j}\boldsymbol {v}_{\mathrm{d}j}}{\rho } = \frac{\rho _{\mathrm{g}}\boldsymbol {v}_{\mathrm{g}}+ \rho _{\mathrm{d}}\boldsymbol {v}_{\mathrm{d}}}{\rho }, \end{eqnarray} (6)ρ is the total density of the mixture,   \begin{eqnarray} \rho \equiv \rho _{\mathrm{g}}+ \rho _{\mathrm{d}}= \rho _{\mathrm{g}}+ \sum _{j}\rho _{\mathrm{d}j}, \end{eqnarray} (7)εj and ε are the mass fractions (relative to the mixture) of the individual and combined dust phases, respectively,   \begin{eqnarray} \epsilon _{j} \equiv \frac{\rho _{\mathrm{d}j}}{\rho }, \end{eqnarray} (8)  \begin{eqnarray} \epsilon \equiv \sum _{j}\epsilon _{j}= \frac{\rho _{\mathrm{d}}}{\rho }, \end{eqnarray} (9)$$\Delta \boldsymbol {v}$$ is the weighted sum of the differential velocities $$\Delta \boldsymbol {v}_{j}\equiv \boldsymbol {v}_{\mathrm{d}j}-\boldsymbol {v}_{\mathrm{g}}$$,   \begin{eqnarray} \Delta \boldsymbol {v}\equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}\Delta \boldsymbol {v}_{j}, \end{eqnarray} (10)$$\boldsymbol {f}$$ represents accelerations acting on both components of the fluid, while $$\boldsymbol {f}_{\mathrm{g}}$$ and $$\boldsymbol {f}_{\mathrm{d}j}$$ represent the accelerations acting on the gas and dust components, respectively, $$\Delta \boldsymbol {f}_{j}\equiv \boldsymbol {f}_{\mathrm{d}j}- \boldsymbol {f}_{\mathrm{g}}$$ is the differential force between the gas and each dust phase, u is the specific thermal energy of the gas, and P is the gas pressure. 2.1.2 Drag time-scales When N = 1, the drag time-scale is unambiguously set by the drag stopping time,   \begin{eqnarray} t_{\mathrm{s}}^{N=1}\equiv \displaystyle\frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}}{K \rho }, \end{eqnarray} (11)where K is a drag coefficient that, in general, depends on local properties of the gas and dust. We assume that K is either constant or in the linear Epstein regime, suitable for small dust grains with low Mach numbers (Epstein 1924, also e.g. Laibe & Price 2012b). In the latter case,   \begin{eqnarray} K = \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}}{\rho _{\rm grain}s} \sqrt{\frac{8}{\pi \gamma }} c_\mathrm{s}= \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}}c_\mathrm{s}}{\rho _\mathrm{eff}s}, \end{eqnarray} (12)where we assume spherical grains with radius s, with uniform intrinsic dust density ρgrain, or equivalently, an effective density $$\rho _\mathrm{eff}\equiv \rho _{\rm grain}\sqrt{\pi \gamma / 8}$$. As usual, γ is the adiabatic constant. The stopping time for N = 1 in the Epstein regime can therefore be written as   \begin{eqnarray} t_{\mathrm{s}}^{N=1}= \frac{\rho _\mathrm{eff}s}{\rho c_\mathrm{s}}. \end{eqnarray} (13) Generalizing the stopping time to N > 1 is conceptually simple, but difficult in practice. Each dust-type equilibrates with the gas at a different rate depending on both the intrinsic properties of the dust grains and the local properties of the gas. Although we assume dust grains of different species do not interact, they are indirectly coupled by their mutual backreaction on the gas. One approach is to derive time-scales using the eigenvalues of the drag matrix LP14c, but the derivations and the expressions become increasingly unwieldy as N increases (i.e. there is no general algebraic expression as a function of N). The eigenvalues help aid in interpreting results, but they are not needed to evolve the fluid equations numerically. The only potential impact the eigenvalues have is through their influence on the time-step. Even then, LP14c found fixed upper/lower bounds to the eigenvalues of the N × N drag matrix, effectively removing any need for the eigenvalues during computation. For convenience, we define the following time-scales to help simplify our numerical implementation:   \begin{eqnarray}{2} T_{\mathrm{s}j} \equiv \frac{\rho _\mathrm{eff}s_j}{\rho c_\mathrm{s}} = \epsilon _{j}\left( 1 - \epsilon \right) t_{j}, \end{eqnarray} (14)  \begin{eqnarray}{2} \widetilde{T}_{\mathrm{s}j} \equiv \frac{T_{\mathrm{s}j}- \sum _{k}\epsilon _{k}T_{\mathrm{s}k}}{1 - \epsilon } \quad && = \epsilon _{j}t_{j}- \sum _{k}\epsilon _{k}^2 t_{k}, \end{eqnarray} (15)where   \begin{eqnarray} t_{j}\equiv \frac{\rho }{K_j}, \end{eqnarray} (16)and where Kj is the drag coefficient for each dust phase, e.g.   \begin{eqnarray} K_j = \frac{\rho _{\mathrm{g}}\rho _{\mathrm{d}j}c_\mathrm{s}}{\rho _\mathrm{eff}s_j}. \end{eqnarray} (17)Note that the weighted sums of equations (14) and (15) happen to be equivalent, i.e.   \begin{eqnarray} \frac{1}{\epsilon } \sum _{j}\epsilon _{j}T_{\mathrm{s}j}\, = \, \frac{1}{\epsilon } \sum _{j}\epsilon _{j}\widetilde{T}_{\mathrm{s}j}\, = \, \frac{1 - \epsilon }{\epsilon } \sum _{j}\epsilon _{j}^2 t_{j}. \end{eqnarray} (18)This new quantity carries physical significance, but its interpretation is clearer if we first define an effective grain size for the mixture,   \begin{eqnarray} s \equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}s_j, \end{eqnarray} (19)such that equation (18) can be written in a more familiar form:   \begin{eqnarray} T_{\mathrm{s}}\equiv \frac{1}{\epsilon } \sum _{j}\epsilon _{j}T_{\mathrm{s}j}= \frac{\rho _\mathrm{eff}s}{\rho c_\mathrm{s}}. \end{eqnarray} (20)Comparing this to equation (13), one may observe that Ts acts like an effective stopping time for the mixture. The benefit of using Tsj and $$\widetilde{T}_{\mathrm{s}j}$$ in lieu of tj is that they allow us to use our existing codebase with only a few additional lines of code, namely to assemble $$\widetilde{T}_{\mathrm{s}j}$$ (Tsj is calculated identically to $$t_{\mathrm{s}}^{N=1}$$ with s replaced by sj). In return, the form of the evolution equations are unchanged from the N = 1 case, as evidenced in the following sections. 2.1.3 Hydrodynamics For the simple case of hydrodynamics, the only force is the pressure gradient, i.e.   \begin{eqnarray} \boldsymbol {f}_{\mathrm{d}j} = 0, \end{eqnarray} (21)  \begin{eqnarray} \boldsymbol {f}_{\mathrm{g}} = - \frac{\nabla P}{\rho _{\mathrm{g}}}, \end{eqnarray} (22)  \begin{eqnarray} \Delta \boldsymbol {f}_{j} = \displaystyle\frac{\nabla P}{\rho _{\mathrm{g}}}. \end{eqnarray} (23)Using equations (14) and (23) to simplify equation (5), we get   \begin{eqnarray} \Delta \boldsymbol {v}_{j}= \frac{\epsilon _{j}t_{j}\nabla P}{\rho } = \frac{T_{\mathrm{s}j}\nabla P}{\rho _{\mathrm{g}}}, \end{eqnarray} (24)while equations (10), (20), and (24) allow us to write   \begin{eqnarray} \Delta \boldsymbol {v}= \frac{T_{\mathrm{s}}\nabla P}{\rho _{\mathrm{g}}}. \end{eqnarray} (25)As promised, when these last two expressions for $$\Delta \boldsymbol {v}_{j}$$ and $$\Delta \boldsymbol {v}$$ are inserted into equations (1), (2), (3), and (4), we obtain the same form of the fluid equations as reported in PL15 for the N = 1 terminal velocity approximation, namely   \begin{eqnarray} \frac{{\rm d} \rho }{{\rm d} t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (26)  \begin{eqnarray} \frac{{\rm d} \epsilon _{j}}{{\rm d} t} = - \frac{1}{\rho } \nabla \cdot \left( \epsilon _{j}\widetilde{T}_{\mathrm{s}j}\nabla P \right), \end{eqnarray} (27)  \begin{eqnarray} \frac{{\rm d} \boldsymbol {v}}{{\rm d} t} = -\frac{\nabla P}{\rho } + \boldsymbol {f}, \end{eqnarray} (28)  \begin{eqnarray} \frac{{\rm d} \tilde{u}}{{\rm d} t} = - \frac{P}{\rho } \nabla \cdot \boldsymbol {v}, \end{eqnarray} (29)where for convenience we have defined $$\tilde{u} \equiv (1 - \epsilon ) u$$ instead of evolving u directly as in PL15. The corresponding energy equation in terms of u would be   \begin{eqnarray} \frac{{\rm d} u}{{\rm d} t} = - \frac{P}{\rho _{\mathrm{g}}} \nabla \cdot \boldsymbol {v}+ \frac{\epsilon T_{\mathrm{s}}}{\rho _{\mathrm{g}}} \nabla P \cdot \nabla u. \end{eqnarray} (30) In order to recover the special case of a single dust phase, we need only collapse the sums in equations (14) and (15) and set sj → s and εj → ε. It is simple to check that in this limit, $$T_{\mathrm{s}j}= \widetilde{T}_{\mathrm{s}j}= T_{\mathrm{s}}= t_{\mathrm{s}}^{N=1}$$, thereby recovering the N = 1 fluid equations from PL15 exactly. 2.1.4 Equation of state The set of equations above is closed by assuming the usual equation of state, which constrains the gas pressure P in terms of the gas density and temperature. Unless otherwise specified in this paper, we assume an adiabatic equation of state, i.e.   \begin{eqnarray} P = \left( \gamma - 1 \right) \rho _{\mathrm{g}}u = \left( \gamma -1 \right) \left( 1 - \epsilon \right) \rho u, \end{eqnarray} (31)or simply   \begin{eqnarray} P = (\gamma -1) \rho \tilde{u}. \end{eqnarray} (32) 2.2 Time-stepping As pointed out by PL15, the addition of the dust evolution equation adds a further constraint on the time-step that becomes limiting when the diffusion coefficient is large. We can derive this time-step constraint more rigorously than that presented by PL15, albeit with the same result for N = 1, by discretizing the set of equations in time using a forward Euler method   \begin{eqnarray} \frac{\rho ^{n+1} - \rho ^n}{\Delta t} = - \rho \left( \nabla \cdot \boldsymbol {v}\right), \end{eqnarray} (33)  \begin{eqnarray} \frac{\epsilon _j^{n+1} - \epsilon _j^{n}}{\Delta t} = - \frac{1}{\rho } \nabla \cdot \left( \epsilon _{j}\widetilde{T}_{\mathrm{s}j}\nabla P \right), \end{eqnarray} (34)  \begin{eqnarray} \frac{\boldsymbol {v}^{n+1} - \boldsymbol {v}^n}{\Delta t} = -\frac{\nabla P}{\rho }, \end{eqnarray} (35)and performing a Von Neumann stability analysis on the above semidiscrete equations. That is, we solve the linear system that results from assuming plane wave solutions of the form   \begin{eqnarray} \rho = D {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (36)  \begin{eqnarray} \boldsymbol {v} = \boldsymbol {V} {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (37)  \begin{eqnarray} \epsilon _j = E_j {\rm e}^{i (\boldsymbol {k}\cdot \boldsymbol {x} - \omega t)}, \end{eqnarray} (38)where D, $$\boldsymbol {V}$$, and Ej are perturbation amplitudes, k is the wavenumber, $$\boldsymbol {x}$$ is the position vector, and ω is the angular frequency. This analysis generically produces a time-step criterion of the form   \begin{eqnarray} \Delta t < C_0\frac{1}{k c_{\max }}, \end{eqnarray} (39)where C0 is a dimensionless safety factor of order unity and cmax  is the maximum wave speed according to the dispersion relation for linear waves. The wavelength of maximum growth usually occurs on the resolution scale, giving the usual Courant criterion   \begin{eqnarray} \Delta t < C_0\frac{h}{c_{\max }}, \end{eqnarray} (40)where h is the SPH smoothing length. For N = 1, the dispersion relation to first order in ωts is given by LP14a,  \begin{eqnarray} \omega = \pm \tilde{c}_{\rm s} k - \frac{i}{2} t_{\rm s} k^2 c_{\rm s}^2 \epsilon , \end{eqnarray} (41)where $$\tilde{c}^2_{\rm s} \equiv c_{\rm s}^2 (1 - \epsilon )$$ is the modified sound speed (squared). The maximum wave speed is therefore   \begin{eqnarray} c_{\max } =\left|\frac{\omega }{k}\right|= \sqrt{\tilde{c}_{\rm s}^2 + \frac{1}{4} \epsilon ^2 t_{\rm s}^2 k^2 c_{\rm s}^4}, \end{eqnarray} (42)and the time-step constraint appropriate for SPH is   \begin{eqnarray} \Delta t < C_0\frac{h}{\sqrt{\tilde{c}^2_{\rm s} + \epsilon ^2 t^2_{\rm s} c_{\rm s}^4 / h^2}}. \end{eqnarray} (43)This is similar to the time-step criterion proposed by PL15 except that the above combines the usual Courant–Friedrichs–Lewy (CFL) condition ($$\Delta t < h/\tilde{c}_{\rm s}$$) and the additional constraint from the dust evolution ($$\Delta t < h^2 / (\epsilon t_{\rm s} c_{\rm s}^2)$$) into a single criterion. When generalizing to multiple dust phases, we find the same result but with the effective stopping time replacing the N = 1 stopping time, giving   \begin{eqnarray} \Delta t < C_0\frac{h}{\sqrt{\tilde{c}^2_{\rm s} + \epsilon ^2 T_{\rm s}^2 c_{\rm s}^4 / h^2}}. \end{eqnarray} (44) As expected, with Ts in the denominator, restricting ourselves to strong drag regimes weakens the constraint on the time-step. More specifically, the time-step is limited when the grain-size distribution is dominated by large grains (or, alternatively, high dust fraction), such that   \begin{eqnarray} \frac{\epsilon T_{\mathrm{s}}}{1-\epsilon } > \Delta t_{\rm CFL}, \end{eqnarray} (45)where $$\Delta t_{\rm CFL}\equiv h/\tilde{c}_{\rm s}$$ is the CFL time-step. The added advantage of the criterion in equation (44) is that it is less stringent than the explicit time-step for either the full multigrain one-fluid formalism or the multifluid method (equations 79 and 80 of LP14c, respectively):   \begin{eqnarray} \Delta t_{\mathrm{one{\rm -}fluid}} < C \Bigg[ \max _{j}\left(\frac{1}{\epsilon _{j}t_{j}} \right) + \frac{1}{\left(1 - \epsilon \right)} \sum _{j}t_{j}^{-1} \Bigg]^{-1}, \end{eqnarray} (46)  \begin{eqnarray} \Delta t_{\mathrm{multi{\rm -}fluid}} < C \min _{j} \left[ \frac{1}{t_{j}} \left( \frac{1}{\epsilon _{j}} + \frac{1}{1-\epsilon } \right) \right]^{-1}, \end{eqnarray} (47)where C is another safety factor. Thus, as long as the cut-off to our dust distribution is ≲cm (see PL15), our global time-step should be of the order of ΔtCFL. 3 SPH FORMULATION When formulating the discretized SPH fluid equations, we can take advantage of the fact that (i) the only equations that were altered by having multiple dust phases were the dust fraction and energy equations and (ii) we have written the continuum equations in the same form as PL15. The first point allows us to adopt the discretized density and momentum equations from PL15 without any changes (thereby guaranteeing exact conservation of linear and angular momentum),   \begin{eqnarray} \rho _{a} = \sum\limits_{b} m_{b} W_{ab} (h_{a}), \end{eqnarray} (48)  \begin{eqnarray} \frac{{\rm d}\boldsymbol {v}_a}{{\rm d}t} & =& -\sum _{b} m_{b} \Big[ \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \nabla _{a} W_{ab} (h_{a}) \nonumber \\ &&+\,\frac{P_{b} + q^{\rm AV}_{ab, b}}{\Omega _{b} \rho _{b}^{2}} \nabla _{a} W_{ab} (h_{b}) \Big] + \boldsymbol {f}_a, \end{eqnarray} (49)where Wab is the usual SPH kernel, h is the smoothing length, Ω is the usual term to account for smoothing length gradients   \begin{eqnarray} \Omega _{a} = 1 - \frac{\mathrm{\partial} h_{a}}{\mathrm{\partial} \rho _{a}} \sum _{b} m_{b} \frac{\mathrm{\partial} W_{ab} (h_{a})}{\mathrm{\partial} h_{a}}, \end{eqnarray} (50)and h is related to ρ in the usual manner [which requires an iterative procedure to solve equation (48); see Price & Monaghan 2004, 2007; LP14b]. The second point allows us to write down the generalized diffusion equation for εj by inspection. Comparing equation (27) to (12) in PL15 suggests that we can use either of their discretized diffusion equations, provided we make the substitutions $$t_{\mathrm{s}}^{N=1}\rightarrow \widetilde{T}_{\mathrm{s}j}$$ and ε → εj (although in the latter case, care must be taken to leave any instances of the gas fraction, 1 − ε, untouched). Furthermore, because evolving the dust fraction directly can in some instances result in negative values, we prefer to use the positive definite formulation prescribed in appendix B of PL15 by defining $$S_j \equiv \sqrt{\rho \epsilon _{j}}$$ (not to be confused with the grain size sj). The corresponding evolution equation in terms of Sj is   \begin{eqnarray} \frac{\mathrm{d} S_{j,a}}{\mathrm{d}t} & =& - \frac{1}{2} \sum _{b} \frac{m_{b} S_{j,b}}{\rho _{b}} \left( \frac{\widetilde{T}_{\mathrm{s}j,a}}{\rho _{a}} + \frac{\widetilde{T}_{\mathrm{s}j,b}}{\rho _{b}} \right) \left(P_{a} - P_{b} \right) \frac{\overline{F}_{ab}}{\vert r_{ab} \vert } \nonumber \\ &&+\, \frac{S_{j,a}}{2\rho _{a}\Omega _{a}} \sum _{b} m_{b} \boldsymbol {v}_{ab} \cdot \nabla _a W_{ab} (h_{a}), \end{eqnarray} (51)where $$\overline{F}_{ab} \equiv \frac{1}{2}[ F_{ab}(h_{a}) + F_{ab}(h_{b}) ]$$ and Fab is defined such that $$\nabla _{a} W_{ab} \equiv F_{ab} \hat{\boldsymbol {r}}_{ab}$$. In writing the diffusion equation in this form, we have implicitly chosen to use the faster, easier-to-implement ‘direct second derivative’ method; however, the evolution equation for the ‘two first derivatives’ method can be obtained in the same fashion (see PL15 for a comparison of these two methods). 3.1 Conservation of energy This leaves only the energy equation to be determined. It is tempting to simply generalize the equation for the energy in a similar manner to the above, but conservation of energy puts an additional constraint on the form of the equation that is not immediately obvious. Instead, we derive the energy equation using the already discretized fluid equations above and by enforcing exact conservation of energy. The total energy E of the system in the terminal velocity approximation can be expressed as   \begin{eqnarray} E = \sum _a m_a \left( \frac{1}{2} v_a^2 + \tilde{u}_a \right), \end{eqnarray} (52)where $$\tilde{u}_a \equiv (1 - \epsilon _a) u_a$$ as previously. Conservation of energy requires that   \begin{eqnarray} \frac{\mathrm{d} E}{\mathrm{d}t} = \sum _{a} m_{a} \left[ \boldsymbol {v}_{a} \cdot \frac{\mathrm{d} \boldsymbol {v}_{a}}{\mathrm{d} t} + \frac{\mathrm{d} \tilde{u}_{a}}{\mathrm{d} t} \right]= 0, \end{eqnarray} (53)where   \begin{eqnarray} \frac{\mathrm{d} \tilde{u}_{a}}{\mathrm{d} t} = \frac{\rho ^\mathrm{g}_{a}}{\rho _a} \frac{\mathrm{d} u_{a}}{\mathrm{d} t} - u_{a} \sum _{j}\left( \frac{2 S_{j,a}}{\rho _{a}} \frac{\mathrm{d} S_{j,a}}{\mathrm{d} t} - \frac{S^2_{j,a}}{\rho ^2_{a}} \frac{\mathrm{d} \rho _{a}}{\mathrm{d}t} \right). \end{eqnarray} (54) Inserting the different expressions from equations (48), (49), and (51) and solving for the time derivative of the energy dictates that the discretized energy equation should be   \begin{eqnarray} \frac{{\rm d}\tilde{u}_{a}}{{\rm d}t} = \sum _{b} m_{b} \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \left( \boldsymbol {v}_{a} - \boldsymbol {v}_{b}\right) \cdot \nabla _{a} W_{ab} (h_{a}), \end{eqnarray} (55)or, if one evolves u directly as in PL15  \begin{eqnarray} \frac{{\rm d}u_{a}}{{\rm d}t} & =& \frac{1}{1 - \epsilon _{a}} \sum _{b} m_{b} \frac{P_{a} + q^{\rm AV}_{ab, a}}{\Omega _{a} \rho _{a}^{2}} \left( \boldsymbol {v}_{a} - \boldsymbol {v}_{b}\right) \cdot \nabla _{a} W_{ab} (h_{a}) \nonumber \\ &&- \, \frac{\rho _{a}}{2 \rho ^\mathrm{g}_{a}} \sum _{j}\sum _{b} m_{b} \frac{S_{j,a} S_{j,b}}{\rho _{a} \rho _{b}} \left( \frac{\widetilde{T}_{\mathrm{s}j,a}}{\rho _{a}} + \frac{\widetilde{T}_{\mathrm{s}j,b}}{\rho _{b}} \right) \nonumber \\ &&\times \, \left( u_{a} - u_{b}\right) \left( P_{a} - P_{b} \right) \frac{\overline{F}_{ab}}{\vert r_{ab} \vert }. \end{eqnarray} (56) 3.2 Shock-capturing terms We include the artificial viscosity and conductive terms below for completeness, but note that they are unchanged by the addition of more dust phases. 3.2.1 Artificial viscosity The artificial viscosity term is computed as follows:   \begin{eqnarray} q^{\rm AV}_{ab, a} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-\frac{1}{2} \left( 1 - \epsilon _{a} \right) v_{\mathrm{sig},a} \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab}, & \qquad \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} < 0\\ 0, & \qquad \mathrm{otherwise}, \end{array}\right. \end{eqnarray} (57)where $$\boldsymbol {v}_{ab} \equiv \boldsymbol {v}_{a} - \boldsymbol {v}_{b}$$ (similarly for $$\hat{\boldsymbol {r}}_{ab}$$) and the signal speed $$v$$sig corresponds to the usual choice for hydrodynamics, i.e.   \begin{eqnarray} v_{\mathrm{sig},a} = \alpha ^{\rm AV}_{a} c_{\mathrm{s},a} + \beta ^{\rm AV} \vert \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} \vert , \end{eqnarray} (58)where $$\alpha ^{\rm AV}_{a} \in [0,1]$$ is the linear dimensionless viscosity parameter (the index implying that αAV can be unique to each particle; see e.g. Morris & Monaghan 1997; Cullen & Dehnen 2010) and βAV (typically βAV = 2) is the von Neumann–Richtmyer viscosity parameter. 3.2.2 Artificial conductivity In order to correctly treat contact discontinuities, an artificial conductivity term must be added to the energy equations (see Price 2008),   \begin{eqnarray} \left( \frac{\mathrm{d} u_{a}}{\mathrm{d} t} \right)_{\rm cond} = \frac{1}{1-\epsilon _{a}} \sum _{b} m_{b} \left[ \frac{Q_{ab, a}}{\Omega _{a} \rho _{a}^{2}} F_{ab} (h_{a}) + \frac{Q_{ab, b}}{\Omega _{b} \rho _{b}^{2}} F_{ab} (h_{b}) \right],\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (59)where   \begin{eqnarray} Q_{ab,a} = \frac{1}{2} \alpha _u \rho _{a} v_{\mathrm{sig},u} \left( u_{a} - u_{b} \right), \end{eqnarray} (60)with αu ∈ [0, 1] the dimensionless conductivity parameter and $$v_{\mathrm{sig},u} = \vert \boldsymbol {v}_{ab} \cdot \hat{\boldsymbol {r}}_{ab} \vert$$ (Price 2008; Wadsley, Veeravalli & Couchman 2008). 4 NUMERICAL TESTS Given the similarity of the SPH equations in Section 3 to those in PL15 and the existing implementation of the latter in our SPH code phantom (e.g. Dipierro et al. 2015; Price et al. 2017), the generalization to N dust phases was straightforward. phantom is well tested (see Price et al. 2017) and we are confident that the implementation of the N = 1 terminal velocity approximation from which we started was correct. Therefore, the tests in this section are less focused on the code as a whole and more focused on specific aspects of our implementation. 4.1 Recovering the N = 1 case By far, the most difficult part of implementing the multigrain method was expanding the relevant arrays in the code to accommodate the N − 1 additional dust phases. To ensure that our new array structures cause no ill effects, we performed 3D versions of the dustyshock, dustywave, and dustydiffuse tests from PL15. Interested readers can find the set-up details for these tests in Price et al. (2017). The results from these tests are shown in the top row of Fig. 1. Importantly, we found that the results calculated with and without our new array structure matched to within machine precision. Note that this agreement shows that equations (48), (49), (51), and (56) reduce numerically to the N = 1 case, analogous to what we observed with the continuum equations. Figure 1. View largeDownload slide Results from the dustyshock (left column), dustywave (middle column), and dustydiffuse (right column) tests as performed by PL15 and Price et al. (2017), but with our new array structure capable of handling multiple dust phases. The top row shows results when N = 1, while the bottom row contains simulations where the same dust phase has been split into N = 10 equal mass bins. As desired, the two cases are identical. Moreover, they match the results run prior to the multigrain implementation. Figure 1. View largeDownload slide Results from the dustyshock (left column), dustywave (middle column), and dustydiffuse (right column) tests as performed by PL15 and Price et al. (2017), but with our new array structure capable of handling multiple dust phases. The top row shows results when N = 1, while the bottom row contains simulations where the same dust phase has been split into N = 10 equal mass bins. As desired, the two cases are identical. Moreover, they match the results run prior to the multigrain implementation. We then added an additional layer of complexity by splitting the single dust phase in each of the above tests into N equal bins and evolving them as if they were N different dust phases. This new set-up can be achieved by setting εj = ε/N and Kj = K. Separating the fluid into mass bins does not alter the physics of the problem, just the numerical method by which it is modelled. Therefore, we should recover the N = 1 solution (modulo numerical errors from calculating and combining quantities differently). Fig. 1 (bottom row) shows the results from the multigrain calculations. Again we found that the dustyshock, dustywave, and dustydiffuse tests agreed with the N = 1 cases to within machine precision. 4.2 Testing the general case It seems like the next logical test would be to extend one or more of the tests above to the general case of N different dust phases. However, there is a fundamental difference in the way the drag is calculated for these tests and the way we have assumed the drag will be calculated when using the equations derived in this paper. Whereas the tests above use a constant drag coefficient K for the entire fluid, the equations in Sections 2 and 3 are optimized for physical dust grains in the Epstein drag regime where the equivalent drag constant (17) changes with grain size. We could reformulate the tests and their solutions to accept a unique value of K for each dust phase, but this would require altering equations (51) and (56) – the very equations we are trying to verify. Therefore, for the general case, we need a test requiring physical grain sizes and drag coefficients. 4.3 Dust settling in a protoplanetary disc The dust settling test from PL15 is an ideal candidate for testing the general case because it mimics one of the environments the multigrain method is designed to simulate, namely the settling of small dust grains in protoplanetary discs. 4.3.1 Initial conditions We simulate a disc-like environment at a radius r = 50 au using a thin, vertical (Cartesian) column of gas in near-hydrostatic equilibrium with an external acceleration in the form of   \begin{eqnarray} \boldsymbol {a}_\mathrm{ext} = - \frac{\mathcal {G} M z}{\left(r^2+z^2\right)^{3/2}}\boldsymbol {\hat{z}}, \end{eqnarray} (61)where $$\mathcal {G}$$ is Newton's gravitational constant, M is the stellar mass, and z is the ‘vertical’ coordinate along the length of the column (x and y represent the two shorter dimensions of the column). The gas density of the column is given by   \begin{eqnarray} \rho _{\mathrm{g}}(z) = \rho _{\mathrm{g},0}\, \mathrm{exp} \left[-\frac{z^2}{2H^2}\right], \end{eqnarray} (62)where we choose H/r = 0.05, giving a disc scale height of H = 2.5 au. We assume an isothermal equation of state with $$P = c_\mathrm{s}^2 \rho _{\mathrm{g}}$$, where cs ≡ HΩ and $$\Omega \equiv \sqrt{\mathcal {G} M / r^3}$$, corresponding to an orbital time $$t_\mathrm{orb} \equiv 2\pi /\Omega \approx 353\, \mathrm{yr}$$. We adopt code units with a distance unit of 10 au, mass in solar masses and time units such that $$\mathcal {G} = 1$$. These choices give an orbital time of ≈70.2 in code units. The particles are initially placed on a close-packed lattice using 100 × 86 × 78 = 670 800 particles in the domain [x, y, z] ∈ [ ± 1, ±0.75, ±0.65]. We then stretch the particles in z using the method described in Price (2004) to give the density profile given in equation (62). We set ρg, 0 to 10−3 in code units (≈6 × 10−13 g cm−3 in physical units), corresponding to a particle mass in code units of 2.42 × 10−9. We use periodic boundary conditions in all directions, but set the boundary in z at ±10H in order to avoid periodicity in the vertical direction. We relaxed the density profile by running the code for 15 orbits with artificial viscosity, at which point we added N = 10 distinct dust phases to the system. We created a cell-edge, logarithmic grid from smin to smax with grid cells of width $$\Delta \log s = \frac{1}{N} \log _{10}\left( s_{\rm max}/s_{\rm min}\right)$$. Then we assigned sj by taking the square root of the product of the cell's endpoints, – thereby skewing the ‘typical’ grain size for each cell towards the smaller, more numerous dust grains. Each dust phase was distributed throughout the disc with an initially uniform dust fraction. We constrained the total dust fraction to be ε = 1/101 (corresponding to a dust-to-gas ratio of 0.01) and set the magnitudes of εj according to the differential power-law distribution   \begin{eqnarray} \mathrm{d} \epsilon = \epsilon _0 s^{3-p} \mathrm{d} s, \qquad \mathrm{for} \quad s_{\rm min}\le s \le s_{\rm max}, \end{eqnarray} (63)where dε is the differential dust fraction with respect to grain size, ε0 is a normalization factor, and p is the usual power-law index for number density as a function of grain size (e.g. Mathis, Rumpl & Nordsieck 1977). In particular, εj is determined by integrating equation (63) across each grain-size cell and then normalizing their combined sum via equation (9). Assuming p = 3.5, we set $$s_{\rm min}\approx 0.0599\, \mu \mathrm{m}$$ and smax ≈ 1.67 mm such that the smallest simulated grain size is $$0.1\, \mu \mathrm{m}$$ and the largest simulated grain size is 1 mm. The initial values for sj and εj in this test are listed in Table 1. Table 1. The initial values for sj and εj used in the settling test assuming a power-law distribution in grain sizes ranging from smin = 0.1 μm to smax = 1 mm with a power-law index of p = 3.5. j  sj [cm]  εj  1  1.00 × 10−5  3.99 × 10−5  2  2.78 × 10−5  6.65 × 10−5  3  7.74 × 10−5  1.11 × 10−4  4  2.15 × 10−4  1.85 × 10−4  5  5.99 × 10−4  3.09 × 10−4  6  1.67 × 10−3  5.15 × 10−4  7  4.64 × 10−3  8.59 × 10−4  8  1.29 × 10−2  1.43 × 10−3  9  3.59 × 10−2  2.39 × 10−3  10  1.00 × 10−1  3.99 × 10−3  j  sj [cm]  εj  1  1.00 × 10−5  3.99 × 10−5  2  2.78 × 10−5  6.65 × 10−5  3  7.74 × 10−5  1.11 × 10−4  4  2.15 × 10−4  1.85 × 10−4  5  5.99 × 10−4  3.09 × 10−4  6  1.67 × 10−3  5.15 × 10−4  7  4.64 × 10−3  8.59 × 10−4  8  1.29 × 10−2  1.43 × 10−3  9  3.59 × 10−2  2.39 × 10−3  10  1.00 × 10−1  3.99 × 10−3  View Large 4.3.2 Results After adding the dust, we ran the simulation for an additional 15 orbits. The resulting dust density for each of the different phases is shown in Fig. 2. As expected, the settling efficiency is proportional to the size of the dust grains, thus enhancing the mid-plane density of the larger grains. However, visually separating this density enhancement is difficult in Fig. 2 because the initial density distribution increases by a factor of 100 from smin to smax. Although the continuum density distribution is a decreasing function with respect to grain size ($$\propto s^{3-p} = 1/\sqrt{s}$$), integrating over each cell to include the mass from non-simulated grains steepens the power law by an additional power of s such that $$\rho _{\mathrm{d}j}\propto \sqrt{s}$$. Figure 2. View largeDownload slide Ten dust densities from a multigrain simulation after having settled for 15 orbits in a 3D vertical column of a protoplanetary disc at r = 50 au (assuming H/r = 0.05; so H = 2.5 au) using 100 × 86 × 78 = 670 800 simulation particles. The grain size and initial dust fraction for each phase is listed in Table 1. Large dust grains efficiently settle towards the disc mid-plane, but still have a much lower density than the smaller dust grains because the global number density of the larger grains is lower. Our multigrain simulation is ∼5 × faster to run than 10 single-phase simulations run serially (see Section 4.3.3). Figure 2. View largeDownload slide Ten dust densities from a multigrain simulation after having settled for 15 orbits in a 3D vertical column of a protoplanetary disc at r = 50 au (assuming H/r = 0.05; so H = 2.5 au) using 100 × 86 × 78 = 670 800 simulation particles. The grain size and initial dust fraction for each phase is listed in Table 1. Large dust grains efficiently settle towards the disc mid-plane, but still have a much lower density than the smaller dust grains because the global number density of the larger grains is lower. Our multigrain simulation is ∼5 × faster to run than 10 single-phase simulations run serially (see Section 4.3.3). In order to better show how the mid-plane density is affected by settling, we set-up and ran a second simulation where the dust fractions were all equal, i.e. εj = ε/N. Fig. 3 shows the resulting time evolution of the dust density for phases j = [1, 9, 10]. This time we clearly see that settling increases the dust density relative to its initial state and at a rate that is commensurate with its settling efficiency. These results are in good agreement with previous settling tests performed in the literature (PL15; Hutchison et al. 2016; Price et al. 2017), albeit with a smaller initial dust fraction. Figure 3. View largeDownload slide Time evolution of the densities of three dust phases (j = [1, 9, 10]). The initial conditions in this simulation were the same as in Fig. 2, except with equal dust fractions (εj = ε/N) to make the relative density enhancement within and between dust phases more visible. We have also adjusted the colourbar in order to allow direct comparison with the settling tests performed by PL15 and Price et al. (2017). Note that the density enhancement due to settling has a shallower dynamic range than the built-in density gamut created by our grains-size distribution (see Fig. 2). Figure 3. View largeDownload slide Time evolution of the densities of three dust phases (j = [1, 9, 10]). The initial conditions in this simulation were the same as in Fig. 2, except with equal dust fractions (εj = ε/N) to make the relative density enhancement within and between dust phases more visible. We have also adjusted the colourbar in order to allow direct comparison with the settling tests performed by PL15 and Price et al. (2017). Note that the density enhancement due to settling has a shallower dynamic range than the built-in density gamut created by our grains-size distribution (see Fig. 2). As a further benchmarking exercise, we ran 10 single-phase simulations using the initial conditions from Table 1 and compared them to the results from the multiphase test above (see Fig. 4). Although the two scenarios are not strictly equivalent – the single-phase simulations do not include the backreaction from the N − 1 other phases – the global solutions still match because (i) the majority of the disc mass resides in the gas and (ii) the gas is essentially motionless throughout the simulation (see e.g. the top row in Fig. 3, which can be used as a proxy for the gas). In the limit of zero backreaction and a stationary gas phase, the system can be modelled analytically and numerically using a simplified set of fluid equations. Section A gives the full analysis. Figure 4. View largeDownload slide Comparison of dust fractions after 15 orbits when calculated by 10 single-phase simulations (black points) versus 1 multigrain simulation (coloured points). Not only does the multigrain method recover the correct solution, but the dispersion in εj is equal to or better than the single-phase simulations. Figure 4. View largeDownload slide Comparison of dust fractions after 15 orbits when calculated by 10 single-phase simulations (black points) versus 1 multigrain simulation (coloured points). Not only does the multigrain method recover the correct solution, but the dispersion in εj is equal to or better than the single-phase simulations. The large-scale agreement, we see in Fig. 4 does not extend to smaller scales. In Fig. 5, we zoom in on the $$s=0.1\, \mu \mathrm{m}$$ grains to illustrate the substructure in εj that develops as a result of the backreaction included from other dust phases. These differences between the single-phase and multigrain simulations continue to grow with time. Therefore, single-phase simulations should be used with caution in situations involving turbulent gas dynamics and/or long time-scales. Figure 5. View largeDownload slide A zoom in of the $$s=0.1\, \mu \mathrm{m}$$ grains in Fig. 4, highlighting the non-linear coupling between dust phases captured in a multigrain simulation (blue points) compared to the single-phase simulation (black points). The location of the peaks (resp. troughs) seen in ε1 correlate with the outer edges (resp. density peaks) of the other phases. We have added semi-transparent lines to help identify the location of the other phases in the figure. With the exception of the largest grain size, the remaining phases exhibit similar discrepancies with their single-phase counterpart. Figure 5. View largeDownload slide A zoom in of the $$s=0.1\, \mu \mathrm{m}$$ grains in Fig. 4, highlighting the non-linear coupling between dust phases captured in a multigrain simulation (blue points) compared to the single-phase simulation (black points). The location of the peaks (resp. troughs) seen in ε1 correlate with the outer edges (resp. density peaks) of the other phases. We have added semi-transparent lines to help identify the location of the other phases in the figure. With the exception of the largest grain size, the remaining phases exhibit similar discrepancies with their single-phase counterpart. The differences we observe in Fig. 5 are small, but prevent the test from being truly rigorous. One way of making the indirect coupling between dust phases vanishingly small is to concentrate all of the dust mass into the smallest grains that remain fixed to the gas, i.e. stationary. We found that by using a power-law index of p = 6.5, we could concentrate $${\sim } 99\, {\rm per \, cent}$$ of the dust mass in the two smallest phases (with j = 1 accounting for $${\sim } 92\, {\rm per \, cent}$$). Under these new conditions, our single-phase and multigrain simulations were a near perfect match at all scales. The only visible difference between the two scenarios was a minor reduction in dispersion in some of the multigrain phases, similar in magnitude to what is seen in Fig. 4. As an interesting aside, the steeper power-law index of p = 6.5 produces a 10 order-of-magnitude gap between the mid-plane densities of the largest and smallest dust grains. Happily, roundoff errors do not appear to corrupt the results in this situation, which we attribute to the fact that each εj is evolved separately. While the dust fractions are combined to calculate the gas properties, the gas–dust interaction depends only on the ratio of their masses. That is, the gas is not sensitive to tiny fluctuations in ε that may be introduced by loss of precision when combining εj of very different magnitudes. So far we have relied on comparing our multigrain results with single-phase simulations. In Fig. 6, we return to using the initial conditions from Table 1 and compare our multigrain solution to a grid-based numerical solution described in Section A2. The settling fronts in our SPH simulations match the simplified solutions to better than a few per cent for all except the largest grain sizes (see Table 2). As pointed out by PL15, the resolution follows the total mass rather than the dust mass, so it tends to over smooth the density peaks in the dust. Despite the smoothing, the locations of the settling fronts and the densities within the disc match very well. The L2 errors scale with the grain size (see column 2 in Table 2) and are a result of the oversmoothed dust peaks and the increased dispersion in the density at larger grain sizes. Figure 6. View largeDownload slide Dust densities for dust settling test in Fig. 4. The coloured points are from our multigrain simulation, while the black solution curves come from solving the equations (A3) and (A4) on a grid. The oversmoothing of the dust fronts is an artefact of tracking the total mass rather than just the dust mass. However, the location of the dust fronts and the enhancement of the mid-plane densities are well captured by the multigrain method. Figure 6. View largeDownload slide Dust densities for dust settling test in Fig. 4. The coloured points are from our multigrain simulation, while the black solution curves come from solving the equations (A3) and (A4) on a grid. The oversmoothing of the dust fronts is an artefact of tracking the total mass rather than just the dust mass. However, the location of the dust fronts and the enhancement of the mid-plane densities are well captured by the multigrain method. Table 2. L2 errors, computed by splash (Price 2007), between the multigrain dusty settling test from Section 4.3 and the analytic/numerical solutions from Section A. Column 2 is obtained using the numeric dust densities from Section A2 and column 3 using the analytic dust velocities from equation (A5). Large L2 errors are caused mainly by inconsistencies between the analytical and numerical models. That our density errors are lowest where our velocity errors are largest (and vice versa) indicates that our multigrain solution is valid. s [cm]  L2 density errors  L2 velocity errors  1.00 × 10−5  4.93 × 10−3  1.47  2.78 × 10−5  4.90 × 10−3  5.27 × 10−1  7.74 × 10−5  4.81 × 10−3  1.89 × 10−1  2.15 × 10−4  5.10 × 10−3  6.80 × 10−2  5.99 × 10−4  5.06 × 10−3  2.45 × 10−2  1.67 × 10−3  6.44 × 10−2  8.90 × 10−3  4.64 × 10−3  1.24 × 10−2  3.50 × 10−3  1.29 × 10−2  3.44 × 10−2  1.92 × 10−3  3.59 × 10−2  7.89 × 10−2  1.61 × 10−3  1.00 × 10−1  9.43 × 10−2  1.57 × 10−3  s [cm]  L2 density errors  L2 velocity errors  1.00 × 10−5  4.93 × 10−3  1.47  2.78 × 10−5  4.90 × 10−3  5.27 × 10−1  7.74 × 10−5  4.81 × 10−3  1.89 × 10−1  2.15 × 10−4  5.10 × 10−3  6.80 × 10−2  5.99 × 10−4  5.06 × 10−3  2.45 × 10−2  1.67 × 10−3  6.44 × 10−2  8.90 × 10−3  4.64 × 10−3  1.24 × 10−2  3.50 × 10−3  1.29 × 10−2  3.44 × 10−2  1.92 × 10−3  3.59 × 10−2  7.89 × 10−2  1.61 × 10−3  1.00 × 10−1  9.43 × 10−2  1.57 × 10−3  View Large In Fig. 7, we compare our multigrain simulation to the analytic solution in equation (A5). While we find L2 errors of order 0.1–1 per cent for grain sizes $$> {}10\, \mu$$m, there is a steady decline in accuracy as grain size decreases (see column 3 in Table 2). This progressive departure from the analytic model is a reflection of the fact that the gas is not completely stationary. Fluctuations in the gas velocity create a size-dependent velocity dispersion in the dust that primarily affects the smaller grain sizes. The larger dust grains, which are less susceptible to these fluctuations in the gas, exhibit less dispersion and better agreement with the analytic solution. Figure 7. View largeDownload slide Settling velocities for the dust settling test in Fig. 4. Again, the coloured points are results from our multigrain calculations, but the black curves are now the single-phase analytic solutions from equation (A5). Slight discrepancies are visible due to the coupling between dust phases and the fact that the gas is not stationary, both of which are ignored in the analytic solution. Figure 7. View largeDownload slide Settling velocities for the dust settling test in Fig. 4. Again, the coloured points are results from our multigrain calculations, but the black curves are now the single-phase analytic solutions from equation (A5). Slight discrepancies are visible due to the coupling between dust phases and the fact that the gas is not stationary, both of which are ignored in the analytic solution. 4.3.3 Performance We ran each of the test simulations above using OpenMP on eight cores from a single node. We found that our multigrain simulations with N = 10 dust phases were a factor of 2 slower than one single-phase simulation, thus making the multigrain simulations five times faster than their single-phase equivalent. This scaling improves as N increases, provided there is enough memory to handle the large array sizes. For example, we found the multigrain method to be ≈13 times faster when N = 100. We expect even better performance ratios relative to multifluid simulations, because multifluid methods require N times more simulation particles and often an added overhead for implicit time-stepping (explicit multifluid methods are impossibly slow for most of the grain sizes considered in this study; see PL15). Finally, because the multigrain method reuses the same simulation particles for all N dust phases, it requires less post-processing and, when N = 10, uses 55 per cent less disc space than an equivalent set of single-phase simulations (65 per cent less when N = 100). Files in which $$\Delta \boldsymbol {v}_{j}$$ is not written to disc1 are reduced by an additional 15 per cent. 4.4 Radial drift in a protoplanetary disc The dusty settling test in the previous section remains well approximated by single-phase methods. To demonstrate that our algorithm also works in regimes of strong backreaction, we computed the radial drift velocities for two dust phases in a protoplanetary disc with conditions such that the inward migration of the larger grains induces a discernable outward migration of the smaller grains. 4.4.1 Analytic solution An analytic solution for multiple dust phases migrating in an inviscid disc was derived by Bai & Stone (2010b). Neglecting vertical gravity, they show that the hydrostatic equilibrium equations can be written in block matrix form as follows:   \begin{eqnarray} {\left(\begin{array}{cc} {\bf I} + \boldsymbol{\Gamma } & -2 \boldsymbol{\Lambda }\\ \boldsymbol {\Lambda }/2 & {\bf I} + \boldsymbol {\Gamma } \end{array}\right)} {\left(\begin{array}{cc}\boldsymbol {\mathcal {V}}_r\\ \boldsymbol {\mathcal {V}}_\phi \end{array}\right)} = - \eta v_K {\left(\begin{array}{cc}\boldsymbol {0}\\ \boldsymbol {1} \end{array}\right)}, \end{eqnarray} (64)where $${\bf I}$$ is the identity matrix, $$\boldsymbol {\mathcal {V}}_r \equiv \left( v_{1r},v_{2r}, \ldots ,v_{nr} \right)^\intercal$$, and $$\boldsymbol {\mathcal {V}}_\phi \equiv \left( v_{1\phi },v_{2\phi }, \ldots ,v_{n\phi } \right)^\intercal$$ are the radial and azimuthal velocities for each dust phase, respectively. The matrix $$\boldsymbol {\Lambda } \equiv {\rm diag}\left\lbrace {\rm St}_1, {\rm St}_2, \ldots , {\rm St}_n \right\rbrace$$ is a diagonal matrix of the Stokes numbers for uncoupled dust phases (i.e. $$\mathrm{St}= t_{\mathrm{s}}^{N=1}\Omega _\mathrm{K}$$), while $$\boldsymbol {\Gamma } \equiv \left(\boldsymbol {\mathcal {E}}, \boldsymbol {\mathcal {E}}, \ldots , \boldsymbol {\mathcal {E}} \right)^\intercal$$ is a matrix made up of the dust-to-gas ratios, where $$\boldsymbol {\mathcal {E}} \equiv \left(\mathcal {E}_1, \mathcal {E}_2, \ldots , \mathcal {E}_n \right)^\intercal$$ and $$\mathcal {E}_j \equiv \rho _{\mathrm{d}j}/\rho _{\mathrm{g}}= \epsilon _{j}/(1-\epsilon )$$. Bai & Stone (2010b) provide a closed-form solution to equation (64); however, we found it more convenient to solve it numerically. 4.4.2 Set-up We set-up a 3D, locally isothermal gas disc using the following power-law parametrizations (see e.g. Laibe, Gonzalez & Maddison 2012)   \begin{eqnarray} c_{{\rm s}}(r) = c_{{\rm s},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{-q/2}, \end{eqnarray} (65)  \begin{eqnarray} H_{\rm g}(r) = H_{{\rm g},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{3/2- q/2}, \end{eqnarray} (66)  \begin{eqnarray} \Sigma _{\rm g}(r) = \Sigma _{{\rm g},1{\rm au}} \left(\displaystyle\frac{r}{1\, {\rm au}} \right)^{-p}, \end{eqnarray} (67)  \begin{eqnarray} \rho _{\rm g}(r,z) = \displaystyle\frac{\Sigma _{\rm g}}{\sqrt{2\pi }H} \exp {\left[ -\displaystyle\frac{z^2}{2 H^2}\right]}, \end{eqnarray} (68)where Σg is the local surface density for the gas, quantities with the subscript ‘1 au’ are reference values measured at r = 1 au, and the parameters p = 1 and q = 0.5 are power-law exponents controlling the density and temperature (i.e. flaring) of the disc, respectively. We set the radial velocity to zero and correct the orbital velocities from pure Keplerian rotation to account for the radial pressure gradient in the disc,   \begin{eqnarray} v_\phi = v_K(1 - \eta ), \end{eqnarray} (69)where the pressure gradient parameter η is given by   \begin{eqnarray} \eta \approx \frac{1}{4} \left(\frac{H_{\rm g}}{r} \right)^2 \left[ 3 + 2 p+ q - \left( 3 - q \right) \left( \frac{z}{H_{\rm g}} \right)^2 \right], \end{eqnarray} (70)to order z2/r2 (see e.g. Takeuchi & Lin 2002). We add the dust by assigning a uniform dust fraction to all of the particles. Because the diffusion approximation assumes the stopping time is much shorter than the dynamical time-scale, we do not give the dust a separate azimuthal velocity. The analytic solution from Bai & Stone (2010b) is 2D and assumes that dust resides in the mid-plane of the disc. As both gas density and gravity decrease with increasing z, we expect the dust at high altitudes to migrate slower than dust in the mid-plane. To compensate, we only compare migration rates for |z| < Hg(r), we bin the particles radially into 50 logarithmically spaced bins (using the same binning method described for the grain-size distribution previously), and we average the radial velocities both azimuthally and vertically within each bin. One final caveat remains: the steady-state analytic model assumes the disc is inviscid, whereas SPH disc simulations are inherently viscous. Normally, we would relax our disc into a quasi-steady state and use our instantaneous gas and dust mid-plane densities as initial conditions for the model – thereby allowing us to account for any non-steady-state processes like settling and/or migration. However, the lack of viscidity in the model produces rigid assumptions about the gas velocity that are not met in our viscous SPH simulations. As a result, we find that our simulation relaxes into a steady state that is substantially different to the analytic solution. To our knowledge, there is currently no analytic solution for radial velocities in viscous discs. Deriving such a solution goes beyond the scope of this paper; therefore, we will revisit the problem in a future study. In the meantime, we can circumvent this incompatibility in the present study by using the initial state of the system, where we have full control over the velocities and we can mimic the conditions of an inviscid disc. Testing the initial conditions, albeit unorthodox, still yields valuable information about our method for two reasons. First, the terminal velocity approximation breaks down when the time-step is smaller than a few stopping times. Because the typical time for drift to relax is on the order of a few stopping times, we do not need to wait for the dust velocities to equilibrate. In other words, the full asymptotic radial velocities are obtained after the very first loop over the particles (what we call t = 0) when P, Tsj, Ts, $$\Delta \boldsymbol {v}_{j}$$, and $$\Delta \boldsymbol {v}$$ are all calculated – the quantities we use to construct the velocity profiles of the gas and dust. Secondly, the individual ga s and dust properties are calculated (as opposed to being evolved). Therefore, the test is more sensitive to how we calculate the forces than how we evolve the mixture. Since our force prescription does not vary with time, the test is almost as useful at t = 0, as it would be once the system has reached a steady state. 4.4.3 Results Using the same grain-size distribution from Section 4.3.1, we set-up a dusty protoplanetary disc around a solar mass star with an inner and outer radius of rin = 1 au and rout = 300 au, respectively. The gas disc has the following reference values: $$c_{{\rm s},1{\rm au}} \approx 1.5 \, {\rm km\,s}^{-1}$$, Hg, 1au = 0.05 au, and $$\Sigma _{{\rm g},1{\rm au}} \approx 166 \, {\rm g\,cm}^{-2}$$. The dust disc for each grain size is set equal to the gas in size and shape, but scaled in mass by the dust fractions listed in Table 3, such that, when all of the dust phases are included, the total dust mass comprises one-third of the total mass of the system (i.e. a dust-to-gas ratio of $$\mathcal {E} = 0.5$$). For the 10 single-phase calculations, it is not possible to simultaneously match the dust fractions, dust-to-gas ratios, and the gas/dust densities of the multigrain case. Therefore, we chose to keep the respective dust fraction and the total surface density of the disc the same, while allowing the dust-to-gas ratio for each dust phase to change as needed. This discrepancy with the multigrain calculations results in different surface densities for the gas and dust, but the effects are unimportant in this context since outward migration of dust in a single-phase simulation can only be achieved by drastically changing the structure of the disc (e.g. with a radially increasing pressure profile). Table 3. The initial grain sizes (s), dust fractions (ε), and dust-to-gas ratios ($$\mathcal {E}$$) used in the outward migration test in Fig. 8. Variables with/without the subscript j indicate multigrain/single-phase values, respectively. The final row gives the sum of the individual dust fractions and dust-to-gas ratios, highlighting the inherent discrepancies between setups of single- and multi-phased simulations. s and sj [cm]  ε and εj  $$\mathcal {E}$$  $$\mathcal {E}_j$$  1.00 × 10−5  1.34 × 10−3  1.34 × 10−3  2.01 × 10−3  2.78 × 10−5  2.24 × 10−3  2.24 × 10−3  3.36 × 10−3  7.74 × 10−5  3.74 × 10−3  3.75 × 10−3  5.60 × 10−3  2.15 × 10−4  6.23 × 10−3  6.27 × 10−3  9.35 × 10−3  5.99 × 10−4  1.04 × 10−2  1.05 × 10−2  1.56 × 10−2  1.67 × 10−3  1.73 × 10−2  1.76 × 10−2  2.60 × 10−2  4.64 × 10−3  2.89 × 10−2  2.98 × 10−2  4.34 × 10−2  1.29 × 10−2  4.83 × 10−2  5.07 × 10−2  7.24 × 10−2  3.59 × 10−2  8.05 × 10−2  8.76 × 10−2  1.21 × 10−1  1.00 × 10−1  1.34 × 10−1  1.55 × 10−1  2.01 × 10−1  Sum total:  1/3  3.65 × 10−1  1/2  s and sj [cm]  ε and εj  $$\mathcal {E}$$  $$\mathcal {E}_j$$  1.00 × 10−5  1.34 × 10−3  1.34 × 10−3  2.01 × 10−3  2.78 × 10−5  2.24 × 10−3  2.24 × 10−3  3.36 × 10−3  7.74 × 10−5  3.74 × 10−3  3.75 × 10−3  5.60 × 10−3  2.15 × 10−4  6.23 × 10−3  6.27 × 10−3  9.35 × 10−3  5.99 × 10−4  1.04 × 10−2  1.05 × 10−2  1.56 × 10−2  1.67 × 10−3  1.73 × 10−2  1.76 × 10−2  2.60 × 10−2  4.64 × 10−3  2.89 × 10−2  2.98 × 10−2  4.34 × 10−2  1.29 × 10−2  4.83 × 10−2  5.07 × 10−2  7.24 × 10−2  3.59 × 10−2  8.05 × 10−2  8.76 × 10−2  1.21 × 10−1  1.00 × 10−1  1.34 × 10−1  1.55 × 10−1  2.01 × 10−1  Sum total:  1/3  3.65 × 10−1  1/2  View Large The left-hand and right-hand panels in Fig. 8 show the mean radial dust velocities for the individual and combined cases, respectively. Coloured points are the velocities calculated from the SPH mixture, while the solid black lines show the corresponding analytic solutions. Although the two largest grain sizes exhibit only minor changes to their velocities after the addition of the other dust species, the eight smaller sizes experience a complete reversal in migration direction. This change in sign is caused by the exchange of angular momentum as the larger grains drag the sub-Keplerian gas into faster orbits, thereby pushing the gas radially outwards. The smaller dust grains, who are more sensitive to changes in the gas flow, are then carried outwards along with the gas. Outward migration of dust in a disc with a radially decreasing pressure gradient cannot be replicated with only one dust phase; conservation of angular momentum requires one or more phases to radially contract as the others expand. Importantly, the multigrain formalism correctly resolves the velocities for both outward and inward migrating species. Figure 8. View largeDownload slide Mean radial drift velocities near the mid-plane of a protoplanetary disc for 10 dust phases of different grain size, individually coupled with the gas (left) and simultaneously coupled (right). Coloured points (connected by either dashed or dotted lines for clarity) show the mean radial dust velocity calculated for each dust phase from the barycentric simulation data and solid black lines indicate the corresponding analytic solutions. The complete lack of outward migration in the left-hand panel accentuates how single-phase simulations can miss dynamical effects caused by the presence of other dust phases. The right-hand panel illustrates that the multigrain formalism is capable of resolving outward/inward migration velocities of different dust species within a single set of simulation particles. Figure 8. View largeDownload slide Mean radial drift velocities near the mid-plane of a protoplanetary disc for 10 dust phases of different grain size, individually coupled with the gas (left) and simultaneously coupled (right). Coloured points (connected by either dashed or dotted lines for clarity) show the mean radial dust velocity calculated for each dust phase from the barycentric simulation data and solid black lines indicate the corresponding analytic solutions. The complete lack of outward migration in the left-hand panel accentuates how single-phase simulations can miss dynamical effects caused by the presence of other dust phases. The right-hand panel illustrates that the multigrain formalism is capable of resolving outward/inward migration velocities of different dust species within a single set of simulation particles. The relative angular velocity between the gas and dust varies with height. In fact, η changes sign at z ≈ 1.5Hg, meaning that dust particles rotate slower than the gas above this height. Because phantom is a 3D code, our calculations systematically underestimate the 2D analytic solution, which assumes all of the dust is rotating in the mid-plane of the disc. We can reduce this offset by only considering particles near the mid-plane, but having fewer particles to average can make the data more noisy. In making Fig. 8, we used 107 particles in the disc, and discard all particles with z > Hg (∼1/3 of the particles). Even with so many particles remaining, the inner ∼20 bins are very noisy (note the first 14 are not shown), with values ranging between −1.5 and 0.3 η$$v$$K. Also note that the standard deviation for most bins is larger than the size of the plotting window, with typical magnitudes ranging from tens to hundreds [ηvK]. Thus, we should not take the discrepancy between the numerical and analytic solutions too seriously. 5 DISCUSSION AND CONCLUSIONS We have derived and implemented a numerical scheme using SPH that is capable of simulating multiple dust phases composed of small dust grains coupled to the gas in the terminal velocity approximation (i.e. when the stopping time is short compared to the computational time-step). Our method simulates dust using a dimensionless dust fraction, as opposed to traditional methods that employ additional sets of simulation particles. By expanding the scalar dust fraction into an array of N dust fractions that are independently evolved and coupled to the gas, we obtain a method that scales better in terms of computational time and resources as N becomes large. Another benefit of evolving the mixture is that the multigrain method circumvents having to resolve the prohibitive temporal and spatial resolution criteria for small dust grains that usually choke multifluid simulations with separate gas and dust particles. We have demonstrated that the multigrain continuum and discretized equations correctly reduce to the equations described by PL15 when N = 1 and that there is no loss in accuracy when simulating a single phase using our multigrain framework – even when that dust phase is divided into multiple mass bins. On the other hand, when simulating multiple unique dust phases, the multigrain method is superior to using multiple single-phase simulations, not only in terms of computational speed and efficiency as discussed above, but also in terms of accuracy as a result of capturing the indirect coupling (via the gas) between dust phases. Although the deviations between our multigrain and single-phase simulations were small for the select test cases we performed in Section 4.3 (∼ few per cent), there are a few additional points to consider for general applications: (i) perturbations from other dust phases accumulate over time, (ii) perturbations from concentrated dust grains (or equivalently, higher dust-to-gas ratios) are larger in magnitude than for dispersed grains, and (iii) perturbations between phases can further be accentuated by motion of the gas (as opposed to the stationary gas phase in our settling tests). In light of these concerns, we caution against using single-phase simulations where possible and encourage the adoption of the more accurate and efficient multigrain method we present here. Finally, the present multigrain algorithm can only be used for small dust grains within the terminal velocity approximation, which is accurate only when the stopping time is shorter than the dynamical time-scale (LP14a). To extend to larger grains, we would need to either implement the full multiphase one-fluid equations with implicit time-stepping from LP14c or develop a hybrid between the one- and multi-fluid methods. Both have advantages and disadvantages, but are beyond the scope of this paper. Presently, we do not account for the evolution in grain size through growth and fragmentation. However, incorporating grain size evolution into the multigrain framework would be straightforward because the mass and number of the simulation particles does not have to change with time. Acknowledgements We thank the anonymous referee whose careful review helped improve this paper significantly. We would also like to acknowledge Matthew Bate, Giovanni Dipierro, and Christophe Pinte for useful discussions. D.J.P. is grateful for funding via an Australian Research Council (ARC) Future Fellowship, FT130100034. M.A.H. and D.J.P. were funded by ARC Discovery Project grant DP130102078. G.L. acknowledges financial support from PNP, PNPS, PCMI of CNRS/INSU, CEA, and CNES, France. This work has in part been carried out within the framework of the National Centre for Competence in Research PlanetS, supported by the Swiss National Science Foundation. Finally, computations have been done on both the SwinSTAR supercomputer hosted at Swinburne University of Technology and the Piz Daint supercomputer hosted at the Swiss National Computational Centre. Footnotes 1 In the diffusion approximation, $$\Delta \boldsymbol {v}_{j}$$ is a calculated quantity needed for recovering the gas and dust velocities during post-processing. However, as $$\Delta \boldsymbol {v}_{j}$$ is not required in any of the evolution equations, we often omit writing it to disc in order to save space. REFERENCES Bai X.-N., Stone J. M., 2010a, ApJS , 190, 297 https://doi.org/10.1088/0067-0049/190/2/297 CrossRef Search ADS   Bai X.-N., Stone J. M., 2010b, ApJ , 722, 1437 https://doi.org/10.1088/0004-637X/722/2/1437 CrossRef Search ADS   Balsara D. S., Tilley D. A., Rettig T., Brittain S. D., 2009, MNRAS , 397, 24 https://doi.org/10.1111/j.1365-2966.2009.14606.x CrossRef Search ADS   Barranco J. A., 2009, ApJ , 691, 907 https://doi.org/10.1088/0004-637X/691/2/907 CrossRef Search ADS   Barrière-Fouchet L., Gonzalez J.-F., Murray J. R., Humble R. J., Maddison S. T., 2005, A&A , 443, 185 CrossRef Search ADS   Chiang E., 2008, ApJ , 675, 1549 https://doi.org/10.1086/527354 CrossRef Search ADS   Cullen L., Dehnen W., 2010, MNRAS , 408, 669 https://doi.org/10.1111/j.1365-2966.2010.17158.x CrossRef Search ADS   Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS , 453, L73 https://doi.org/10.1093/mnrasl/slv105 CrossRef Search ADS   Epstein P. S., 1924, Phys. Rev. , 23, 710 https://doi.org/10.1103/PhysRev.23.710 CrossRef Search ADS   Garaud P., Lin D. N. C., 2004, ApJ , 608, 1050 https://doi.org/10.1086/420839 CrossRef Search ADS   Haworth T. J. et al.  , 2016, PASA , 33, e053 https://doi.org/10.1017/pasa.2016.45 Hutchison M. A., Price D. J., Laibe G., Maddison S. T., 2016, MNRAS , 461, 742 https://doi.org/10.1093/mnras/stw1126 CrossRef Search ADS   Jacquet E., Balbus S., Latter H., 2011, MNRAS , 415, 3591 https://doi.org/10.1111/j.1365-2966.2011.18971.x CrossRef Search ADS   Johansen A., Klahr H., 2005, ApJ , 634, 1353 https://doi.org/10.1086/497118 CrossRef Search ADS   Laibe G., Price D. J., 2012a, MNRAS , 420, 2345 https://doi.org/10.1111/j.1365-2966.2011.20202.x CrossRef Search ADS   Laibe G., Price D. J., 2012b, MNRAS , 420, 2365 https://doi.org/10.1111/j.1365-2966.2011.20201.x CrossRef Search ADS   Laibe G., Price D. J., 2014a, MNRAS , 440, 2136 (LP14a) https://doi.org/10.1093/mnras/stu355 (LP2014a) CrossRef Search ADS   Laibe G., Price D. J., 2014b, MNRAS , 440, 2147 (LP14b) https://doi.org/10.1093/mnras/stu359 (LP2014b) CrossRef Search ADS   Laibe G., Price D. J., 2014c, MNRAS , 444, 1940 (LP14c) https://doi.org/10.1093/mnras/stu1367 (LP2014c) CrossRef Search ADS   Laibe G., Gonzalez J.-F., Maddison S. T., 2012, A&A , 537, A61 CrossRef Search ADS   Lee A. T., Chiang E., Asay-Davis X., Barranco J., 2010, ApJ , 718, 1367 https://doi.org/10.1088/0004-637X/718/2/1367 CrossRef Search ADS   Lorén-Aguilar P., Bate M. R., 2014, MNRAS , 443, 927 https://doi.org/10.1093/mnras/stu1173 CrossRef Search ADS   Lorén-Aguilar P., Bate M. R., 2015, MNRAS , 454, 4114 https://doi.org/10.1093/mnras/stv2262 CrossRef Search ADS   Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ , 217, 425 https://doi.org/10.1086/155591 CrossRef Search ADS   Miniati F., 2010, J. Comput. Phys. , 229, 3916 https://doi.org/10.1016/j.jcp.2010.01.034 CrossRef Search ADS   Monaghan J. J., 1997, J. Comput. Phys. , 138, 801 https://doi.org/10.1006/jcph.1997.5846 CrossRef Search ADS   Monaghan J. J., Kocharyan A., 1995, Comput. Phys. Commun. , 87, 225 https://doi.org/10.1016/0010-4655(94)00174-Z CrossRef Search ADS   Morris J. P., Monaghan J. J., 1997, J. Comput. Phys. , 136, 41 https://doi.org/10.1006/jcph.1997.5690 CrossRef Search ADS   Paardekooper S.-J., Mellema G., 2004, A&A , 425, L9 CrossRef Search ADS   Price D. J., 2004, PhD thesis , Univ. Cambridge, Cambridge Price D. J., 2007, PASA , 24, 159 https://doi.org/10.1071/AS07022 CrossRef Search ADS   Price D. J., 2008, J. Comput. Phys. , 227, 10040 https://doi.org/10.1016/j.jcp.2008.08.011 CrossRef Search ADS   Price D. J., Federrath C., 2010, MNRAS , 406, 1659 Price D. J., Laibe G., 2015, MNRAS , 451, 813 (PL15) https://doi.org/10.1093/mnras/stv996 (PL15) CrossRef Search ADS   Price D. J., Monaghan J. J., 2004, MNRAS , 348, 139 https://doi.org/10.1111/j.1365-2966.2004.07346.x CrossRef Search ADS   Price D. J., Monaghan J. J., 2007, MNRAS , 374, 1347 https://doi.org/10.1111/j.1365-2966.2006.11241.x CrossRef Search ADS   Price D. J. et al.  , 2017, PASA, preprint (arXiv:1702.03930) Saffman P. G., 1962, J. Fluid Mech. , 13, 120 https://doi.org/10.1017/S0022112062000555 CrossRef Search ADS   Takeuchi T., Lin D. N. C., 2002, ApJ , 581, 1344 https://doi.org/10.1086/344437 CrossRef Search ADS   Tricco T. S., Price D. J., Laibe G., 2017, MNRAS , 471, L52 https://doi.org/10.1093/mnrasl/slx096 CrossRef Search ADS   Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS , 387, 427 https://doi.org/10.1111/j.1365-2966.2008.13260.x CrossRef Search ADS   Yang C.-C., Johansen A., 2016, ApJS , 224, 39 https://doi.org/10.3847/0067-0049/224/2/39 CrossRef Search ADS   Youdin A. N., Goodman J., 2005, ApJ , 620, 459 https://doi.org/10.1086/426895 CrossRef Search ADS   Youdin A., Johansen A., 2007, ApJ , 662, 613 https://doi.org/10.1086/516729 CrossRef Search ADS   APPENDIX A: SOLUTIONS TO THE SETTLING TEST In the limit of very small dust-to-gas ratios, we can neglect the backreaction of the dust on the gas. The dust can then be treated as N independent phases, moving inside a static gas background, and governed by the following one-dimensional equations   \begin{eqnarray} \frac{\mathrm{\partial} \rho _{\mathrm{d}}}{\mathrm{\partial} t} + v_{\mathrm{d}}\frac{\mathrm{\partial} \rho _{\mathrm{d}}}{\mathrm{\partial} z} = - \rho _{\mathrm{d}}\frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} z}, \end{eqnarray} (A1)  \begin{eqnarray} \frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} t} + v_{\mathrm{d}}\frac{\mathrm{\partial} v_{\mathrm{d}}}{\mathrm{\partial} z} = - \frac{v_{\mathrm{d}}}{t_{\mathrm{s}}^{N=1}} - \frac{\mathcal {G} M z}{\left( r^2 + z^2 \right)^{3/2}}, \end{eqnarray} (A2)where we have dropped the subscript j to emphasize that the phases are no longer coupled. To aid our analysis, we define the dimensionless variables $$\bar{\rho } \equiv \rho /\rho _{\mathrm{g},0}$$, $$\bar{v} \equiv v/v_\mathrm{K}$$, $$\bar{z} \equiv z/r$$, and $$\bar{t} \equiv t \, \Omega _\mathrm{K}$$, where $$v_\mathrm{K}= \sqrt{\mathcal {G} M/r}$$ and ΩK = $$v$$K/r are the Keplerian velocity and frequency, respectively. Substituting these quantities into equations (A1) and (A2), we obtain the corresponding non-dimensionalized equations in the form   \begin{eqnarray} \frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{v}_\mathrm{d}\frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{z}} = - \bar{\rho }_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}}, \end{eqnarray} (A3)  \begin{eqnarray} \frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{v}_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}} = - \frac{\bar{v}_\mathrm{d}}{\mathrm{St}} - \frac{\bar{z}}{\left(1 + \bar{z}^2 \right)^{3/2}}, \end{eqnarray} (A4)where $$\mathrm{St}= t_{\mathrm{s}}^{N=1}\Omega _\mathrm{K}$$ is the Stokes number. A1 Analytic solution to settling problem Importantly, equation (A4) is independent of $$\bar{\rho }_\mathrm{d}$$. As a first-order partial differential equation, a solution for $$\bar{v}_\mathrm{d}$$ could potentially be obtained via the method of characteristics. We took a simpler approach by solving the Lagrangian form of the equation with a convective derivative. In this form, the equation is a first-order ordinary differential equation that can be solved using an integrating factor. Assuming the dust starts from rest, the solution for the dust velocity is   \begin{eqnarray} \bar{v}_\mathrm{d}= \frac{\mathrm{St}\, \bar{z}}{\left(1+\bar{z}^2\right)^{3/2}} \left( \mathrm{e}^{-\bar{t}/\mathrm{St}} - 1\right). \end{eqnarray} (A5) The same procedure can be used to obtain an equation for the dust density. However, the resulting solution does not conserve mass, since it assumes that an infinitely extended dust distribution continuously rains down on to the disc. The density solution nevertheless correctly predicts the location of the incoming dust front and also the interior density profile so long as t is less than the settling time-scale. This is not a problem in the velocity solution because all of our dust grains settle at their terminal velocity within the disc, effectively erasing any built-up momentum gained at higher altitudes. A2 1D numerical solution to the settling problem To remove the assumptions imposed by the analytic solution, we also compared our multigrain results with a numerical solution to equations (A3) and (A4). We solved the equations on a one-dimensional grid using an implicit Crank–Nicolson algorithm, using forward differences in time and centred differences in space. Because the temporal derivative is centred half a time-step in the future, we replace all of the other terms with time averages centred about the same time. Grouping terms based on their location in time, we can then write each equation as a linear system in the form   \begin{eqnarray} \bf {A} \boldsymbol {x}^{n+1} = \bf {B} \boldsymbol {x}^n + \boldsymbol {C}, \end{eqnarray} (A6)where superscripts designate the time level, $$\bf {A}$$ and $$\bf {B}$$ are square sparse matrices, $$\boldsymbol {x}$$ is the fluid variable for which we are solving, and the column vector $$\boldsymbol {C}$$ is a placeholder for all terms independent of $$\boldsymbol {x}$$ ($$\boldsymbol {C}=0$$ when $$\boldsymbol {x}$$ represents density). Once the boundary conditions have been accounted for in $$\bf {A}$$ and $$\bf {B}$$, the solution at time n + 1 can be obtained symbolically via   \begin{eqnarray} \boldsymbol {x}^{n+1} = \bf {A}^{-1} \left( \bf {B} \boldsymbol {x}^n + \boldsymbol {C} \right). \end{eqnarray} (A7)The non-linearity in the advection term in equation (A4) keeps us from obtaining the solution using the exact method as outlined above because it would irreversibly mix terms from different time-steps. To overcome this problem, we assume the leading $$\bar{v}_\mathrm{d}$$ in the advection term is known and designate it as $$\tilde{v}_\mathrm{d}$$ to keep it separate from the other velocity terms. We account for $$\tilde{v}_\mathrm{d}$$ and the fact that the fluid equations are coupled by using a predictor-corrector scheme to advance the system forward in time. Designating predicted quantities with asterisks, we advance the system in four steps: $$\bar{\rho }_\mathrm{d}^*$$ is predicted assuming $$\bar{v}_\mathrm{d}$$ is constant (i.e. $$\bar{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^n$$). $$\bar{v}_\mathrm{d}^*$$ is predicted assuming that $$\tilde{v}_\mathrm{d}$$ is constant and equal to $$\bar{v}_\mathrm{d}^n$$. $$\bar{\rho }_\mathrm{d}^{n+1}$$ is corrected assuming $$\bar{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^*$$. $$\bar{v}_\mathrm{d}^{n+1}$$ is corrected assuming $$\tilde{v}_\mathrm{d}^{n+1} = \bar{v}_\mathrm{d}^*$$ and $$\tilde{v}_\mathrm{d}^n = \bar{v}_\mathrm{d}^n$$. Using the same physical parameters as in Section 4.3, we discretize the region z ∈ [−3H, 3H] with 1002 cell-centred grid points, including ghost points. The boundary condition for the velocity is $$v$$d( ± 3H, t) = 0, which consequently enforces the following boundary condition for the density:   \begin{eqnarray} \frac{\mathrm{\partial} \bar{\rho }_\mathrm{d}}{\mathrm{\partial} \bar{t}} + \bar{\rho }_\mathrm{d}\frac{\mathrm{\partial} \bar{v}_\mathrm{d}}{\mathrm{\partial} \bar{z}} = 0 \end{eqnarray} (A8)at the same locations. The initial conditions are $$v$$d(z, 0) = 0 and ρd(z, 0) = εjρg. We found that a dimensionless time-step of 1 was sufficient to keep the algorithm stable. As the dust settles, we do get some low-density numerical noise in the wings of the disc, but this noise is always separated from the settling dust layer by a region of zero density. We have verified that our results do not change when we force the density to zero beyond the first encountered zero-density grid point on either side of the mid-plane. The close match between our multigrain results and the analytic and numerical solutions demonstrates that our method works. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial