Skip to content

Dust modelled as Lagrangian particles

Solver and general options

When dust is modelled as Lagrangian particles, we treat them as test particles, which means that they feel the gravity of the star (direct and indirect terms), of the planets (if any), of the gas disc (if gas self-gravity is included) and the drag acceleration from the gas (both the Epstein and the Stokes drag laws are implemented with a simple interpolation between both regimes as in Paardekooper 2007 (A&A, 462 ; ADS). However, dust particles do not feel each other (no dust self-gravity, no collisions, no growth nor fragmentation). Dust drag onto the gas (dust feedback) is actually implemented in the present version of the code, but only in a beta version which is valid only for same-sized particles (it will not work for a size distribution). The effects of gas turbulence are modelled by applying stochastic kicks to the orbital radius and azimuth of the particles at each hydrodynamical timestep. Kicks follow a Gaussian distribution with mean and standard deviation as in Charnoz et al. (2011, ApJ, 737 ; ADS – see their equation 21, where we discard the terms with spatial derivatives of the dust’s turbulent diffusion, which are generally small). The coefficient of dust’s turbulent diffusion, \(D\), is taken to be \(D = \nu(1 + 4St^2 )/(1 + St^2 )^2\) , with \(\nu\) the local kinematic viscosity of the gas, and \(St\) the local Stokes number of the dust (Youdin & Lithwick 2007, Icarus, 192 ; ADS). Note that it is possible to discard any of the above forces in the code’s parameter file (parameters DustFeelPlanets, DustFeelSG, DustFeelDisk, DustFeedback and DustFeelTurb). The particles number is set by NbPart in the input parameter file.

The equations of motion for the particles are solved with a semi-implicit first-order Euler integrator which is similar to, but different from, that implemented in the code Athena (Zhu et al. 2014, ApJ, 785 ; ADS). Differences between both integrators are found to be marginal, except when the particles friction time is much shorter than the hydrodynamical timestep. When the friction time is much shorter than the hydro timestep, the so-called short-friction time approximation may be used instead of the particles integrator (e.g., Johansen & Klahr 2005, ApJ, 634 ; ADS – see SFTApprox in the input parameter file, which can be deactivated, for instance for testing purposes). We have checked that the particles integrator behaves well for dust particles that have a friction time much shorter than the hydro timestep, and gives very similar results to the short-friction time approximation. Note also that the use of a leapfrog integrator (second-order in time) has a minor impact on the results for typical applications of the code (gas/dust simulations over a few thousand orbits at most).

In the calculation of the gas drag acceleration, the gas quantities need to be interpolated at the particles location. By default, a Triangular-Shaped Cloud interpolation scheme is used, but other interpolation schemes are also available : Cloud-In-Cell and Nearest-Grid Point – usually, Cloud-In-Cell provides good enough results. The same interpolation scheme is used when calculating the feedback acceleration of the dust particles on the gas. The interpolation scheme can be changed via the Interpolation parameter in the input parameter file. The functions that take care of the particles update are SemiUpdateDustPositions and UpdateDustVelocities in src/DustUpdate.c. The interpolation scheme is implemented in function interpolation in src/Dsys.c. Note that by default, particles that enter the planets Hill radius are removed from there and put back in the disc at a different location. This avoids a large accumulation of particles immediately near the planets. If you are interested in the particles’ behavior near the planets, you can deactivate this feature by setting RemoveDustFromPlanetsHillRadius to No in the input parameter file.

Initial conditions

The initial orbital radius of the dust particles is sorted out according to a power-law probability distribution, whose power-law exponent is set by DustSlope in the input parameter file (DustSlope should be set to SigmaSlope-1 if you would like a uniform dust-to-gas surface density ratio initially). The minimum and maximum orbital radii where particles are inserted in the grid are also defined in the input parameter file (RMinDust and RMaxDust). The initial azimuth of the particles is sorted out randomly between \(0\) and \(2\pi\). Their initial radial velocity is zero1, and their initial azimuthal velocity accounts for the star’s gravity and the disc’s self-gravity (if included). This ensures that the only difference in initial azimuthal velocities between gas and dust is due to the gas pressure gradient. The function that takes care of the particles initialization is InitDustSystem in the src/Dsys.c source file.

Size distribution

The particles have a size distribution such that the number of particles in the size interval \([s, s + ds]\), which we denote by \(n(s)ds\), is a power-law function of \(s\). The negative value of the power-law exponent is set by SizePartSlope in the input parameter file (it therefore takes positive values). The minimum and maximum sizes of the particles are set in physical units (meters) in the input parameter file (SizeMinPart and SizeMaxPart). Given that NbPart is usually up to about \(10^6\) to maintain a tractable computational cost, it may not be possible to use a realistic, ISM-like value for SizePartSlope of \(3.5\). This is not a problem for observational predictions of the dust’s thermal emission (see Baruteau et al. 2019, MNRAS, 486 ; ADS). Taking SizePartSlope = 1 ensures that there is roughly the same particles number per decade of size.

Internal density and code units

The internal density of the particles is set in physical units by the Rhopart parameter in the input parameter file (in \(\rm{g.cm^{−3}}\)). This implies that the code’s units of length and mass need to be specified in the input parameter file. By default, the code’s unit of length is \(1\) Astronomical Unit (AU) and the code’s unit of mass is \(1\) Solar mass (the units of time and of temperature follow by the convention that the gravitational constant and the reduced ideal gas constant are equal to unity in code units ; see function ComputeCodeUnits in src/LowTasks.c). The code’s units of length and mass can be changed via the parameters FactorUnitLength and FactorUnitMass in the input parameter file. For example, in in/template.par, FactorUnitLength is set to \(10\) so that the code’s unit of length is \(10\) AU (and the grid thus extends from \(4\) to \(22\) AU).

Outputs

There are two kinds of outputs the code writes for the dust particles. On the one hand, the code writes dustsystatX.dat ascii files, which contain \(6\) columns : the particles orbital radius, azimuth, radial velocity, azimuthal velocity, Stokes number and physical size in meters (this requires the keyword WriteDustSystem set to Yes in the input parameter file). The first four columns are in code units. On the other hand, the code may write dustX.dat binary files, which contain 2D arrays of the total dust surface density (requires the keyword WriteDustDensity set to Yes in the input parameter file).

Parallelization

The particles evolution is parallelized with MPI using FARGO’s domain decomposition into rings, which means that the particles’ position and velocity are communicated from one CPU to another when particles cross the radial boundary between two CPUs. You should be aware that, on some very rare occasions, this may cause the code’s execution to suddenly ’freeze’ if at some point there are too many particles that cross the boundary between two CPUs : the MPI messages to be sent from one CPU to another are too big and cause MPI Wait instructions to fail, with the consequence that the code’s execution is put on hold indefinitely. Should this happen, please restart the code by changing the number of CPUs (and check that the RestartWithNewDust parameter is set to No in the input parameter file). I am well aware of this problem but still have not found a smart way to solve it.

  1. This could be improved in the future by using the initial gas radial velocity interpolated at the particles position. ↩