Enforcing dust mass conservation in 3D simulations of tightlycoupled grains with the Phantom SPH code
Abstract
We describe a new implementation of the onefluid method in the SPH code Phantom to simulate the dynamics of dust grains in gas protoplanetary discs. We revise and extend previously developed algorithms by computing the evolution of a new fluid quantity that produces a more accurate and numerically controlled evolution of the dust dynamics. Moreover, by limiting the stopping time of uncoupled grains that violate the assumptions of the terminal velocity approximation, we avoid fatal numerical errors in mass conservation. We test and validate our new algorithm by running 3D SPH simulations of a large range of disc models with tightly and marginallycoupled grains.
keywords:
hydrodynamics, dust dynamics – methods: numerical – one fluid – accretion, accretion discs.1 Introduction
Protoplanetary discs are composed of a mixture of gas and dust. While gas usually dominates the mass, and hence the hydrodynamics of the system, dust is the dominant source of opacity in the bulk of the disc. As a result, the optical appearance of discs is strongly influenced by the dust distribution (Testi et al., 2014; Birnstiel et al., 2016). Recent highresolution observations of protoplanetary discs have revealed a wealth of asymmetric structures in both gas and dust phases (e.g. Casassus, 2016; Boehler et al., 2018). The physical mechanisms driving the formation of these structures are best understood using 3D hydrodynamical simulations that accurately model the coupling between gas and dust for a wide range of grain sizes (Haworth et al., 2016).
Solid particles embedded in a gas fluid are often treated using a continuous fluid description (Garaud et al., 2004). The macroscopic properties of the dust (e.g. density and velocity) are evolved on a set of grid points or particles that represent a volume large enough to be statistically meaningful, but sufficiently small as to ignore variations of the fluid quantities within that volume. In Smoothed Particle Hydrodynamic (SPH) simulations, the dust dynamics can be computed using two different approaches: the twofluid algorithm described in Laibe & Price (2012), typically used for large dust grains in weak drag regimes, and the onefluid algorithm (Price & Laibe, 2015) based on the socalled terminal velocity approximation (Youdin & Goodman, 2005), which is better suited for simulating dust phases that are tightly coupled with the gas. In terms of gasdust modelling, the twofluid implementation treats the gas and the dust as two separate sets of simulation particles, coupled by a drag force. In contrast, the SPH particles in the onefluid approach represent the mixture, whose composition is determined by the dust fraction, that is evolved as a local property of the mixture.
Since it is numerically difficult to simulate all of the physical drag regimes that occur in nature with a single algorithm, methods/studies are often distinguished by the degree of coupling between phases, usually quantified by the socalled Stokes number: particles with the same Stokes number are aerodynamically identical – regardless of their shape, size, and/or density. The Stokes number is found by comparing the typical dynamical timescale of the system, , to the typical stopping timescale, , i.e. the time it takes for drag to significantly modify the relative velocity between a single grain and the gas. When the grains size is smaller than the mean free path of the gas (Epstein, 1924) – which is generally the case for mmcm size grains in protoplanetary discs (Garaud et al., 2004) –, the Stokes number is given by (Price & Laibe, 2015)
(1) 
where is the grain size, is the intrinsic grain density, is the sound speed, is the adiabatic index, is the Keplerian angular velocity (), is a correction factor for supersonic drag and is the total density.
1.1 The onefluid method
The onefluid equations can be derived by rewriting the fluid equations for the gas and the dust in the barycentric reference frame of the mixture (Laibe & Price, 2014a). In doing so, we substitute out the individual velocities of the gas and dust phases in favour of the new barycentric velocity of the mixture,
(2) 
and the differential velocity between the two phases,
(3) 
where is the density, is the velocity, and the subscripts g and d identify gas and dust quantities, respectively. Similar to the velocities, we replace the gas and dust densities by the total density, , and the dust fraction, , such that
(4)  
(5) 
The equations describing the evolution of a dustgas mixture can be therefore written in the form (Laibe & Price, 2014a)
(6)  
(7)  
(8)  
(9)  
(10) 
where is the gas pressure and represents the external forces acting on both components (e.g. gravity). Moreover, for convenience we have parametrised the thermal energy as . The stopping time, , is given by
(11) 
where is the drag coefficient, which regulates the aerodynamical coupling between the two phases (Weidenschilling, 1977). The equations of the mixture are closed by the equation of state, such as the adiabatic one, i.e.
(12) 
There are several advantages to using the onefluid formulation over the twofluid approach (see Price & Laibe 2015), particularly for small dust grains. For example, since the gas and dust are colocated in the onefluid approach, it does not require (or can easily circumvent) the prohibitive temporal and spatial resolution requirements at high drag (needed in twofluid simulations by the interpolation of fluid quantities between different phases, Laibe & Price 2012). Furthermore, the onefluid method prevents artificial trapping of dust beneath the resolution length of the gas. Finally, the onefluid formalism naturally generalises to account for multiple dust species coupled to the same gas phase (Laibe & Price, 2014c; Hutchison et al., 2018).
Terminal Velocity Approximation
The fluid equations in the onefluid formalism can be simplified when the stopping time is small compared to the typical hydrodynamic timescale, i.e. the time required for a sound wave to propagate over a characteristic distance. In the context of SPH, we can write this condition as, , where is the local smoothing length of the particles. In this regime, usually referred to as terminal velocity regime (Youdin & Goodman, 2005), the relative velocity between the two phases rapidly reach a terminal velocity due to the balancing of the drag and pressure forces. As a consequence, the time dependence of the differential velocity between the gas and the dust can be ignored,
(13) 
Neglecting terms of second order in , Eqs. (6)(10) reduce to
(14)  
(15)  
(16)  
(17) 
Apart from the additional evolution equation for the dust fraction, the equations in the terminal velocity approximation bear striking resemblance to the usual hydrodynamic equations for the gas without the dust. Therefore, the SPH discretisation of the continuity and momentum equations are identical to that of a regular gasonly simulation while the dust fraction and the thermal energy are discretised directly as shown in Eq. 43 in Price & Laibe (2015) and Eq. 55 in Hutchison et al. (2018).
1.2 Timestepping
The addition of the diffusion equation for the dust fraction (Eq. 15) leads to an additional constraint on the timestep. Assuming a constant density and an isothermal equation of state, , Eq. (15) can be rewritten as
(18) 
where is the diffusion coefficient. A new constraint on the timestep is needed when the diffusion coefficient is larger. Indeed, Price & Laibe (2015) provide a stability criterion of the form
(19) 
which implies that the timestep needs to be constrained when the stopping time is long – the opposite of the twofluid case where the timestep is constrained for short stopping times. It is worth remarking that the terminal velocity approximation is only strictly valid when the stopping time is less than the computational timestep. Actually, a more general timestep condition can be derived, taking into account possible gradients in (see Appendix A). This is given as:
(20) 
where and . It can be easily seen that Eq. (20) reduces to Eq. (19) for constant . This condition is safer than Eq. (19) in regions of strong gradients of , but it is more difficult to implement (since it requires an additional loop over the particles to obtain the gradient of ) and can lead to severe timestep restrictions in certain practical applications (see Sect. 4). As a result, we default back to Eq. (19) for our timestep control in this work.
2 Enforcing positivity of the dust fraction
The onefluid approach does not put any constraint on the positivity of the dust fraction. This problem can arise in regions where, for example, particles containing a finite amount of dust are adjacent to pure gas particles (i.e. ). As the particles evolve in time, the infinite gradient in created at this interface leads the pure gas particles to develop a negative dust fraction. We can avert this problem by parameterising and evolving the dust fraction using a new variable, . The positivity of the physical variable is now guaranteed since
(21) 
The corresponding diffusion equation for the new variable is
(22) 
We note that the first term on the right hand side of Eq. 2 is written so as to prevent an infinite gradient in when (i.e. ). The usual method for discretising Eq. 14,
(23) 
trivially conserves the total mass of the mixture, but does nothing to conserve the mass of each of the components. Formally, mass conservation of the dust and gas also holds as long as the energy equation is modified appropriately Price et al. (2017), i.e. such that
(24) 
which, in terms of the new variable , requires that
(25) 
The SPH discretisation for the evolution of is shown in Eq. 280 of Price et al. (2017). Although the formulation prevents from going negative, it does not guarantee that the dust fraction will remain smaller than unity. Numerical artefacts can appear in regions where the gradient of the dust fraction is steep, resulting in a spontaneous increase in dust mass. These artefacts are most severe when or and, at least in some instances, quickly drive the dust fraction to values larger than unity.
3 A new implementation
In this section we propose a new parametrization of the dust fraction similar to that used by Price & Laibe (2015), but that enforces the constraint by mapping the dust fraction to a function whose codomain is only defined from , thereby preventing from becoming unphysical. A promising parametrization that meets the above criterion is given by
(26) 
In this new formulation the variable is then related simply to the ratio of dust to gas densities, . We calculate the time derivative as
(27) 
Substituting Eq. (7) and manipulating the term on the right hand side of Eq. (27), we obtain
(28) 
The SPH discretisation is implemented in the form
(29) 
where . Like the previous implementation (Sect. 2), our new expressions conserve linear and angular momentum, energy, and mass — at least up to the accuracy of the timestepping algorithm. Although it is true that the total mass is trivially conserved by virtue of Eq. (23), this attribute is not bequeathed to the individual phases due to their dependence on , an evolved quantity. This contingency on the timeevolution accuracy of plays an important role in the discussion that follows.
We implemented the above formalism into the SPH code Phantom (Lodato & Price, 2010; Price et al., 2017) and tested it using Phantom’s standard nightly test suite (described in Sect. 5.1 of Price et al. 2017), which includes (among others) the dustywave, dustyshock, and dustydiffuse tests described in Price & Laibe (2015). The new implementation not only passed within the ‘acceptable’ tolerances set for each test, it outperformed the existing algorithm. As a specific example, when compared with the previous method, the ‘derivatives test’ (Sect. 5.1.1 of Price et al. 2017) showed an improvement in the accuracy of the time derivative of by a factor of five while the total energy conservation improved by a factor of .
Next we looked at some typical configurations involving the interaction of an embedded protoplanet with its parent disc. Comparing the two parametrizations discussed in this paper, Fig. 1 follows the evolution of the dust surface density (initial power law profile with index ) in a 3D simulation of a dusty protostellar disc with a radial extent of au and an embedded planet of mass 0.5 located at a distance of 60 au from the central star. The temperature profile drops as a power law with and the disc aspect ratio is , at au. We embed the planet in order to further investigate diffusivity gradients that arise due to planetdisc interactions. The planet also alters the relative dust fractions in the inner and outer parts of the disc with time. The simulation describes the evolution of a millimeter grain population. Particles with a nonnegligible dust fraction exhibit Stokes numbers in the range , which safely correspond to stopping times below .
As time progresses, the viscous and pressure forces in the disc cause the gas to expand radially outward, creating a strong gradient in the dust fraction (and hence diffusivity) at the edge of the dusty disc. Fig. 1 shows that the numerical artefacts that occurred with the old implementation are removed with the new parametrization. This improved accuracy is thanks to the more accurate timeevolution of in regions with steep gradients in the dust diffusivity (i.e., at the outer edge of the dusty disc).
Again comparing the two implementations, Fig. 2 shows the time evolution of the total dust mass, i.e (dashed line) and (solid line). While with the old implementation the dust mass increases in time (starting from a value of and reaching , after years), the new implementation better computes the evolution of the dust density, avoiding most of the numerical artefacts occurring at the edge of the dusty disc due to the strong gradients in the dust fraction. Moreover, our tests show that the computation of the dust fraction and the thermal energy in our new implementation is faster than the parametrization described in Sect. 2.
4 Limiting the stopping time
As mentioned earlier, despite the conservation ensured by the spatial discretisation of the fluid equations, nonconservation may still arise due to timestepping errors. Nonconservation of gas/dust mass are particularly vulnerable in regions of small where the dust fraction tends to relax the timestep (see Eq. 19). However, since these regions are usually occupied by dust grains with large stopping times, they are the very regions that need a small timestep in order to be accurate. This breakdown of our timestep criterion is most likely due to the violation of the assumptions used to derive Eq. (19), and in particular to the fact that it was derived neglecting gradents in the dust fraction, as discussed already in Sect. 2 above. In theory, we should be able to reduce our timestep (by adopting the full timestep condition, Eq. 20, or by reducing ) to maintain mass conservation. We have verified that maintaining a ‘sufficiently small’ timestep for these problematic particles preserves mass conservation for the system, but at the cost of impossibly slow simulations when, e.g., very small amounts of dust get flung out and stranded in the lowdensity outer disc. Therefore, in practice we seek a more viable option that can circumvent these problem particles while still conserving gas/dust mass for the system. It is rather vexing that such violations most likely occur in ‘peripherial’ particles that often have little influence on the simulation at large. From experience, numerical artefacts are mostly likely to occur in the upper/outer regions of discs with high aspect ratio, , and low (in absolute value) radial powerlaw index for the temperature, . The dust diffusion, i.e. , in these regions is strong due to the steep gradients in the pressure and for particles with large stopping time.
To prevent the numerical inaccuracies we see when such strong gradients are present in the disc for particles with large stopping time, we propose moderating the rapid dust diffusion for problematic particles by enforcing the following limit on the stopping time
(30) 
that results in limiting the flux of the mass embodied in large particles. Limiting the flux of dust mass through the stopping time (as opposed to the pressure gradient or the diffusion coefficient as a whole) has the advantage that it is localised strictly to particles that violate the terminal velocity approximation and requires no prior knowledge about the dynamical state of the system.
Fig. 3 compares the evolution of the dust surface density using our new dust implementation presented in Sect. 3 (lower panels) and the same implementation, but limiting the stopping time (upper panels). The protoplanetary disc used in these two simulations has a radial extent of au and it is thicker than the one used in Fig. 1, with an aspect ratio, , at au and a temperature profile index . The gaseous disc mass is , with a dusttogas ratio of . The initial gas and dust surface densities are given by a power law (index ) with an exponential cutoff at au. The dust grain size is . In this case, we include two planets at au and au, of and , respectively. The outer planet is deliberately placed so as to fling dust into regions where we know the terminal velocity approximation has difficulty. Importantly, the flux limited simulations conserve the dust mass to machine precision while our other simulations do not. The spurious increase in dust mass in our unmodified simulation takes place in the outer disc where the gradients in the dust diffusivity are large.
It is important to note that by limiting the stopping time we are artificially modifying the Stokes number. Rewriting Eq. 30 in terms of the Stokes number (for discs in vertical hydrostatic equilibrium, i.e. ) yields
(31) 
Since in typical SPH simulations , this new implementation affects the dust density evolution of large dust grains, even with moderately low . Since the radial dust velocity increases with for (Nakagawa et al., 1986), limiting the stopping time leads to an underestimate of the radial flux of large grains towards disc regions of high pressure. Consequently, the newfound mass conservation afforded by limiting the flux is not an excuse to apply our method in every situation. In particular, care should be taken when simulating protoplanetary discs with high aspect ratio and low , where it is more likely to find dust grains with both large and small Stokes number. For these discs, a correct physical description of the system may only be attainable with the full onefluid approach (Laibe & Price, 2014a), hybrid method combining the one and twofluid approaches or semianalytical twofluid methods (e.g. LorénAguilar & Bate, 2014).
In summary, limiting the stopping time conserves dust mass and prevents numerical artefacts from developing in particles in the outer disc where Stokes numbers are large and dust mass is negligible. The evolution of the dust density in these situations can be considered reliable. However, when the stopping time of particles are being limited in the bulk of the disc, where mass fractions are still high, we recommend using a different approach.
5 Conclusions
We introduce a new algorithm to compute the dynamics of tightlycoupled dust grains in the context of the one fluid approach described in Laibe & Price (2014a). Our algorithm avoids certain numerical artefacts that arise in the previous formalism (Price & Laibe, 2015), rendering our method both faster and more accurate. We do this by

parameterising the dust fraction using the square root of the dusttogas ratio, which enforces ;

limiting the stopping time below a value that ensures the validity of the equations of motion in the terminal velocity approximation, i.e. .
The latter leaves the numericallystable, stronglycoupled dust grains untouched, while limiting the amount of dust that can be transferred between weaklycoupled particles that would otherwise violate the assumptions of the onefluid diffusion approximation. When the flux in these weaklycoupled grains is not constrained, the dust mass can unphysically grow over long times in some regions of the disc, violating mass conservation. We find no adverse effects of limiting the flux of particles with low dust fraction, which are typically found in the upper/outer regions of the disc. However, we caution that the stopping time limiter needs to be used with care, since it can lead to an incorrect computation of the dust dynamics of large decoupled dust grains when the dust fraction is nonnegligible. In these situations, we recommend switching to a twofluid formalism (Laibe & Price, 2012; LorénAguilar & Bate, 2014).
Finally, there are realistic scenarios in which a single grain size can be stronglycoupled in one region of the disc and weaklycoupled in another – with a significant dust mass in each region. In this scenario, neither the onefluid diffusion approximation or the twofluid method would be adequate, but would require a hybrid scheme that marries the two approaches or, alternatively, the full onefluid formalism that allows for a wider range in drag regimes (Laibe & Price, 2014a, b). Alternatively, implicit or semianalytic methods have been proposed to simulate tightlycoupled particles in multifluid simulations with strong drag regimes (LorénAguilar & Bate, 2014; Booth et al., 2015; LorénAguilar & Bate, 2015).
Acknowledgments
We are grateful to an anonymous referee for constructive suggestions that improved our paper. The authors would like to thank Richard Alexander and Chris Nixon for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 681601). This research used the ALICE High Performance Computing Facility and the DiRAC Complexity system, funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. We used SPLASH (Price, 2007).
References
 Birnstiel et al. (2016) Birnstiel T., Fang M., Johansen A., 2016, Space Sci. Rev., 205, 41
 Boehler et al. (2018) Boehler Y., et al., 2018, ApJ, 853, 162
 Booth et al. (2015) Booth R. A., Sijacki D., Clarke C. J., 2015, MNRAS, 452, 3932
 Casassus (2016) Casassus S., 2016, Publ. Astron. Soc. Australia, 33, e013
 Epstein (1924) Epstein P. S., 1924, Physical Review, 23, 710
 Garaud et al. (2004) Garaud P., BarriéreFouchet L., Lin D. N. C., 2004, ApJ, 603, 292
 Haworth et al. (2016) Haworth T. J., et al., 2016, preprint, (arXiv:1608.01315)
 Hutchison et al. (2018) Hutchison M. A., Price D. J., Laibe G., 2018, preprint, (arXiv:1802.03213)
 Laibe & Price (2012) Laibe G., Price D. J., 2012, MNRAS, 420, 2345
 Laibe & Price (2014a) Laibe G., Price D. J., 2014a, MNRAS, 440, 2136
 Laibe & Price (2014b) Laibe G., Price D. J., 2014b, MNRAS, 440, 2147
 Laibe & Price (2014c) Laibe G., Price D. J., 2014c, MNRAS, 444, 1940
 Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
 LorénAguilar & Bate (2014) LorénAguilar P., Bate M. R., 2014, MNRAS, 443, 927
 LorénAguilar & Bate (2015) LorénAguilar P., Bate M. R., 2015, MNRAS, 454, 4114
 Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
 Price (2007) Price D. J., 2007, PASA, 24, 159
 Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
 Price et al. (2017) Price D. J., et al., 2017, preprint, (arXiv:1702.03930)
 Testi et al. (2014) Testi L., et al., 2014, Protostars and Planets VI, pp 339–361
 Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
 Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
Appendix A Time stepping with gradients of
For simplicity, we derive the time step condition with nonzero derivatives of the dust fraction in 1D first, from
(32) 
For the first order backward Euler scheme, the linear expansion of Eq. 32 for modes of the form provides
(33) 
The numerical scheme requires for stability. With the usual substitution , this condition gives
(34) 
where , and . Hence,
(35) 
Expanding the lefthand side of Eq. 35 and dividing by provides finally
(36) 
Putting a safety constant in front of the righthand side of Eq. 36 gives the generic form for the time step condition with gradients of , which can be generalised in 3D accordingly.