# Dust growth, fragmentation and self-induced dust traps in Phantom

Arnaud Vericel1, Jean-François Gonzalez1, Daniel J. Price2, Guillaume Laibe1
and Christophe Pinte2,3
1Univ Lyon, Univ Claude Bernard Lyon 1, Ens de Lyon, CNRS, Centre de Recherche
Astrophysique de Lyon, UMR5574, F-69230, Saint-Genis-Laval, France
2School of Physics & Astronomy, Monash University, Clayton VIC 3800, Australia
3Univ. Grenoble Alpes, CNRS, IPAG, F-38000, Grenoble, France E-mail:
<EMAIL_ADDRESS>
(Accepted 2021 July 30. Received 2021 July 30; in original form 2020 September
17)

###### Abstract

We present the implementation of a dust growth and fragmentation module in the
public Smoothed Particle Hydrodynamics (SPH) code Phantom. This module is made
available for public use with this paper. The coagulation model considers
locally monodisperse dust size distributions around single values that are
carried by the SPH particles. Along with the presentation of the model,
implementation and tests, we showcase growth and fragmentation in a few
typical circumstellar disc simulations and revisit previous results. The
module is also interfaced with the radiative transfer code mcfost, which
facilitates the comparison between simulations and ALMA observations by
generating synthetic maps. Circumstellar disc simulations with growth and
fragmentation reproduce the ‘self-induced dust trap’ mechanism first proposed
by Gonzalez et al., which supports its existence. Synthetic images of discs
featuring this mechanism suggest it would be detectable by ALMA as a bright
axisymmetric ring at several tens of au from the star. With this paper, our
aim is to provide a public tool to be able to study and explore dust growth in
a variety of applications related to planet formation.

###### keywords:

protoplanetary discs — hydrodynamics — methods: numerical — planets and
satellites: formation

††pubyear: 2021††pagerange: Dust growth, fragmentation and self-induced dust
traps in Phantom–B.2

##  1 Introduction

High resolution imaging of protoplanetary discs with the Atacama Large
Millimetre/Submillimetre Array (ALMA) has revealed rich substructure in mm
continuum emission that includes rings, gaps and spiral arms (Andrews et al.,
2018; Huang et al., 2018). The current state of the art in modelling these
discs is to simulate gas plus one, or more, dust populations of constant size
(e.g. Dipierro et al., 2015; Price et al., 2018b; Mentiplay et al., 2019;
Pérez et al., 2019; Calcino et al., 2019; Veronesi et al., 2020). Gas and dust
are coupled via an aerodynamic drag force (Whipple, 1973). Whereas gas orbits
at a sub-Keplerian speed because of its own pressure support, the pressureless
dust phase feels an azimuthal headwind and loses angular momentum. This
transfer leads to the dust drifting radially towards the star and settling
vertically towards the mid-plane (Safronov, 1969; Adachi et al., 1976;
Weidenschilling, 1977; Dubrulle et al., 1995; Haghighipour, 2005). These
processes are most efficient for mm to cm grain sizes, i.e. for dust
populations detected by instruments such as ALMA. While simulating gas, dust
and their mutual effects, calculations indicate that interaction with embedded
planets is the likely explanation for the observed rings and gaps
(Paardekooper & Mellema, 2004; Dipierro & Laibe, 2017; Mentiplay et al., 2019;
Pérez et al., 2019; Toci et al., 2020). Stellar companions are also invoked to
explain spiral arms, cavities or misalignments (Kley & Nelson, 2012; Zhu et
al., 2015; Muñoz & Lai, 2016; Dong et al., 2016; Price et al., 2018b; Dong et
al., 2018; Calcino et al., 2019; Cuello et al., 2019, 2020; Nealon et al.,
2020; Gonzalez et al., 2020). While this general simulation approach has
proven successful, particularly for planet-induced kinematic signatures (Pinte
et al., 2018, 2019, 2020), it fails to explain how planets were formed in
these discs in the first place, or how the dust grows and evolves on longer
time scales.

From a theoretical point of view, dust evolution has been studied extensively
in 3D without dust growth to understand dust radial drift, settling or
concentrations in pressure bumps (e.g. Weidenschilling, 1977; Dubrulle et al.,
1995; Haghighipour, 2005; Fromang & Papaloizou, 2006). Efficient dust
concentrations have been shown to lead to the Streaming Instability — a local
hydrodynamic instability leading to direct planetesimal formation (Youdin &
Goodman, 2005; Johansen & Youdin, 2007; Yang et al., 2017; Schäfer et al.,
2017; Auffinger & Laibe, 2018; Abod et al., 2019; Li et al., 2019b). This
mechanism could be the cornerstone of planet formation. Understanding how
particles can reach such high concentration levels in discs is thus of utmost
importance. One mechanism that could provide the required initial conditions
is the ‘self-induced dust trap’ proposed by Gonzalez et al. (2017a) and
Vericel & Gonzalez (2020), where dust growth/fragmentation, back-reaction and
large scale gradients leads to dust pile-ups and subsequent formation of a
local pressure maximum in the disc.

Dust growth is a complex and computationally challenging process in which
coagulation of small dust particles forms larger bodies by a ‘snowball effect’
(Lissauer & Stewart, 1993; Dominik & Tielens, 1997; Blum & Wurm, 2008). Since
the pioneering work of Weidenschilling (1997), the numerical challenges of
dust growth have led the community to follow different approaches to tackle
this problem. Grid-based codes have considered the resolution of the
Smoluchowski equation (Smoluchowski, 1916) in vertically and azimuthally
averaged discs (Brauer et al., 2008; Birnstiel et al., 2009; Pinilla et al.,
2012; Birnstiel et al., 2012; Drążkowska et al., 2016; Drążkowska & Alibert,
2017) where gas undergoes Keplerian shear and the vertical equilibrium is
reached instantaneously. In most of these models the gas is considered as a
static background and back-reaction, i.e. the drag force exerted by dust onto
the gas, is neglected, although it can be of significant importance during the
disc evolution (Gonzalez et al., 2017a; Dipierro et al., 2018). However, a few
2D grid-based codes compute the Smoluchowski equation at the same time as the
evolution of both gas and dust, more particularly for the study of planet
hosting discs (Li et al., 2005, 2019a; Drążkowska et al., 2019; Laune et al.,
2020). Alternatively, dust growth has also been modelled with Monte Carlo
methods in grid-based codes (e.g. Zsom & Dullemond, 2008; Ormel & Spaans,
2008; Drążkowska et al., 2014), although these need a large number of size
bins to avoid overdiffusion.

While the Eulerian nature of a grid-based code handles dust size distributions
naturally in each cell, particle tracking and complex geometries present
challenges. As a result, alternative methods have also been developed with a
static gas background on top of which dust super-particles are set (Krijt et
al., 2016; Krijt & Ciesla, 2016; Schoonenberg et al., 2018; Misener et al.,
2019). In these studies, dust growth is modelled using the ‘single-size’
approximation, i.e. considering that every dust super-particle carries a
single size that can evolve because of the local conditions.

This method was pioneered by Laibe et al. (2008, hereafter L08) following
Stepinski & Valageas (1997), using a 3D SPH code that solves the equations of
motion self-consistently, including gas drag and back-reaction. This code,
which we will refer to as LyonSPH, initially developed by Barrière-Fouchet et
al. (2005) has been further developed over the years to include fragmentation
(Gonzalez et al., 2015; Gonzalez et al., 2017a), porosity (Garcia, 2018;
Garcia & Gonzalez, 2020) and snow line (Vericel & Gonzalez, 2020) effects. In
this paper, we extend and adapt this growth algorithm to Phantom: a public,
fast, modular and optimised SPH code widely used by the community (Price et
al., 2017, 2018a, hereafter P18). By doing so, our aim is to provide a tool to
study dust growth in a range of planet-forming environments, including discs
around multiple stars, flybys, and in misaligned and warped discs.

The paper is organised as follows. Section 2 describes the growth and
fragmentation model while Section 3 treats its numerical implementation. We
present two tests in Section 4 as well as circumstellar disc simulations in
Section 5. Finally, we discuss our findings and conclude in Sections 6 and 7.

##  2 Model

Dust coagulation is the result of the collision between particles, the outcome
of which is determined by their relative velocity, noted $V_{\mathrm{rel}}$.
The relative velocity can come from multiple sources including disc
turbulence, Brownian motion, dust radial and azimuthal drift and dust
settling. Since we use the ‘single size’ approximation locally, particles in
that vicinity feel the same drag force, which cancels the radial, azimuthal
and vertical component of their relative velocity (see Vericel & Gonzalez,
2020, for a discussion). Moreover, the Brownian motion only plays a
significant role for very small grains (micrometre and smaller, Birnstiel et
al., 2010), which is a regime that is in practice never met in our
simulations. For these reasons, we only consider the turbulence-driven
component of the relative velocity, which is transmitted to the dust from the
gas via drag. We use the form proposed by Stepinski & Valageas (1997), that is

$V_{\mathrm{rel}}=\sqrt{2^{3/2}\mathrm{Ro}\alpha}\dfrac{\sqrt{\mathrm{Sc}-1}}{\mathrm{Sc}}c_{\mathrm{s}}=\sqrt{2}V_{\mathrm{t}}\dfrac{\sqrt{\mathrm{Sc}-1}}{\mathrm{Sc}},$
(1)

where $\mathrm{Ro}$ is the Rossby number that we consider constant and equal
to 3, $\alpha$ is the viscosity parameter defined by Shakura & Sunyaev (1973)
and $c_{\mathrm{s}}$ the gas sound speed. Various turbulent relative velocity
formulations have been proposed, we refer the reader to L08 and Laibe (2014)
for a discussion about their differences. In Equation (1),
$V{{}_{\mathrm{t}}}$ is called the turbulent velocity and $\mathrm{Sc}$ is the
Schmidt number, expressed as

$\mathrm{Sc}=\left(1+\mathrm{St}\right)\sqrt{1+\dfrac{\Delta\bm{v}^{2}}{V_{\mathrm{t}}^{2}}},$
(2)

where $\mathrm{St}$ is the Stokes number and $\Delta\bm{v}$ is the
differential velocity between the gas and dust phases. The dust relative
velocity (Equation 1) is plotted for different values of $\Delta v$ in Fig. 1.
For small differential velocities with respect to $V_{\mathrm{t}}$, the
Schmidt number approximates as $1+\mathrm{St}$, which reduces the relative
velocity to

$V_{\mathrm{rel}}\simeq\frac{\sqrt{2\mathrm{St}}}{1+\mathrm{St}}V_{\mathrm{t}}.$
(3)

The Stokes number is defined as the product of the stopping time
$t_{\mathrm{s}}$ set by the gas drag and the Keplerian frequency
$\Omega_{\mathrm{k}}$ set by the star’s gravitational field

Figure 1: Evolution of the relative velocity (Equation 1) as a function of the
Stokes number for 4 different values of the ratio $\Delta v/V_{\mathrm{t}}$.
In the simulations performed for this paper, the differential velocity is
always small compared to the turbulent velocity (and by extension compared to
the sound speed), thus we only show values in the range $0\leq\Delta v/V_{\rm
t}\leq 1$.

$\mathrm{St}=t_{\mathrm{s}}\Omega_{\mathrm{k}},$ (4)

where $t_{\mathrm{s}}$ depends on the drag regime considered (see Laibe &
Price, 2012b). In the usual Epstein regime (Epstein, 1924), this gives

$\mathrm{St}=\sqrt{\dfrac{\pi\gamma}{8}}\dfrac{\rho_{\mathrm{s}}s}{f\left(\rho_{\mathrm{g}}+\rho_{\mathrm{d}}\right)c_{\mathrm{s}}}\Omega_{\mathrm{k}},$
(5)

where $\gamma$ is the adiabatic index, $\rho_{\rm s}$ is the dust intrinsic
density, and $\rho_{\rm g}$ and $\rho_{\rm d}$ are the volume densities of the
gas and dust phases. The definition of the Stokes number corresponds to the
mixture (Dipierro & Laibe, 2017; Kanagawa et al., 2017). $f$ is the Kwok
parameter (Kwok, 1975) correcting drag for supersonic differential velocities
and expressed as

$f=\sqrt{1+\dfrac{9\pi}{128}\dfrac{\Delta\bm{v}^{2}}{c_{\mathrm{s}}^{2}}}.$
(6)

The radial drift process is fastest when $\mathrm{St}\sim 1$. Considering
this, we rewrite the Stokes number as the ratio between the dust size and its
optimal drift size for which $\mathrm{St}=1$, called $s_{\rm opt}$, which can
be expressed as

$\displaystyle s_{\rm opt}$
$\displaystyle=\sqrt{\dfrac{8}{\pi\gamma}}\dfrac{f(\rho_{\mathrm{g}}+\rho_{\mathrm{d}})c_{\mathrm{s}}}{\rho_{\mathrm{s}}\Omega_{\mathrm{k}}},$
(7)
$\displaystyle=\sqrt{\dfrac{8}{\pi\gamma}}f(1+\varepsilon)\dfrac{\rho_{\mathrm{g}}}{\rho_{\mathrm{s}}}H,$
(8)

where $\varepsilon=\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ is the dust-to-gas
ratio and $H=c_{\mathrm{s}}/\Omega_{\mathrm{k}}$ is the scale height of the
gas disc.

###  2.1 Growth

In the model, we consider a single dust SPH particle to represent a swarm of
equal-sized physical grains. This means that at the location of the SPH
particle, the size distribution is considered highly peaked around a certain
value: the one that the particle carries. We represent grain growth inside
each SPH particle by assuming perfect coagulation of two grains during a
collision time. Since the distribution is locally monodisperse (i.e. the size
is locally uniform), only same-size coagulation events are modelled and the
grain mass $m_{\mathrm{d}}$ doubles during a collision

$\dfrac{\mathrm{d}m_{\rm d}}{\mathrm{d}t}=m_{\rm d}f_{\mathrm{col}},$ (9)

where $f_{\mathrm{col}}=\sigma_{\mathrm{d}}n_{\mathrm{d}}V_{\mathrm{rel}}$ is
the collision frequency, $\sigma_{\mathrm{d}}$ is the dust cross section and
$n_{\mathrm{d}}$ is the dust numerical density. We consider compact spherical
grains, for which the growth rate is more conveniently expressed using the
size

$\dfrac{\mathrm{d}s}{\mathrm{d}t}=\dfrac{s}{3}f_{\mathrm{col}}.$ (10)

The growth rate is associated with the physical grains within the SPH particle
and not to the SPH particle itself. This is fundamental, because we wish to
conserve the mass of each SPH particle to preserve numerical robustness. This
constraint is satisfied considering a larger size corresponds to fewer
physical grains represented by a given SPH particle, without affecting its
total mass. This representation holds as long as the mass of the SPH particle
itself is much larger than the mass of a single physical grain it ‘contains’.
This translates as the condition

$s\ll\left(\dfrac{M_{\mathrm{d}}}{n^{\mathrm{SPH}}_{\mathrm{d}}}\dfrac{3}{4\pi\rho_{\mathrm{s}}}\right)^{\frac{1}{3}},$
(11)

where $M_{\mathrm{d}}$ is the total dust mass, and
$n^{\mathrm{SPH}}_{\mathrm{d}}$ is the number of dust SPH particles.
Quantitatively speaking, this gives

$s\ll
362\left(\dfrac{M_{\mathrm{d}}}{10^{-4}~{}M_{\odot}}\right)^{\frac{1}{3}}\left(\frac{n^{\mathrm{SPH}}_{\mathrm{d}}}{10^{6}}\right)^{-\frac{1}{3}}\left(\dfrac{\rho_{\mathrm{s}}}{1000~{}\mathrm{kg.m}^{-3}}\right)^{-\frac{1}{3}}\mathrm{km},$
(12)

which is entirely satisfied for simulations with sizes between a few tens of
micrometres to a few tens of kilometres at most. Equation (10) finally gives
the growth rate

$\dfrac{\mathrm{d}s}{\mathrm{d}t}=\dfrac{\rho_{\mathrm{d}}}{\rho_{\mathrm{s}}}V_{\mathrm{rel}}.$
(13)

###  2.2 Fragmentation

While Equation (13) describes the growth process, fragmentation can also be
considered in the model. More specifically, as the grains relative velocity
increase, the kinetic energy of the collision also increases. At sufficiently
high relative velocities, the kinetic energy can destroy the grains chemical
bonds and make them fragment to smaller sizes. A threshold relative velocity
above which grains are fragmenting upon collisions can then be considered.
This threshold is often called the ‘fragmentation velocity’, noted
$V_{\mathrm{frag}}$. In case of fragmentation, the resulting growth rate is
therefore negative and expressed as

$\dfrac{\mathrm{d}s}{\mathrm{d}t}=-\dfrac{\rho_{\mathrm{d}}}{\rho_{\mathrm{s}}}V_{\mathrm{rel}}\psi,$
(14)

where $\psi$ refers to two models of fragmentation:

  1. 1.

‘Hard’, that is symmetrical to the growth case (Gonzalez et al., 2015);

  2. 2.

‘Smooth’, which considers a more progressive loss of mass with increasing
relative velocities (Garcia, 2018).

Defining the ratio of the relative velocity to the fragmentation velocity as
$v=V_{\mathrm{rel}}/V_{\mathrm{frag}}$, the parameter $\psi$ is given by

$\psi=\left\\{\begin{array}[]{cl}1&\mbox{ : Hard fragmentation},\\\
v^{2}/(1+v^{2})&\mbox{ : Smooth fragmentation},\\\ \end{array}\right.$ (15)

The evolution of $\psi$ with the relative to fragmentation velocity ratio is
shown in Fig. 2. The Hard model considers that most of the mass is lost during
fragmentation, whatever the relative velocity. Conversely, the Smooth model
only considers large loss of mass for very high relative velocities. Both
models are implemented in Phantom and can be selected by the user at runtime.

Figure 2: Evolution of the fragmentation parameter $\psi$ as a function of the
ratio $v=V_{\mathrm{rel}}/V_{\mathrm{frag}}$ for both Hard and Smooth
fragmentation models. The grey zone corresponds to growth
($V_{\mathrm{rel}}<V_{\mathrm{frag}}$), where $\psi$ is not defined.

####  2.2.1 Fragmentation sizes and distance to the star

When considering fragmentation, we can estimate the equilibrium size that the
dust particles can reach when their relative velocity is in equilibrium with
the fragmentation velocity. Following Gonzalez et al. (2017a), by writing the
equality $V_{\rm rel}=V_{\rm frag}$, we find more precisely

$\sqrt{2\tilde{\alpha}}c_{\rm
s}\dfrac{\sqrt{\mathrm{St}}}{1+\mathrm{St}}=V_{\rm frag},$ (16)

where $\tilde{\alpha}=2^{1/2}\mathrm{Ro}\alpha$ and Sc is considered equal to
$1+\mathrm{St}$. Using this leads to the following quadratic equation

$\mathrm{St}^{2}+2\left(1-\dfrac{\tilde{\alpha}c_{\rm s}^{2}}{V_{\rm
frag}^{2}}\right)\mathrm{St}+1=0,$ (17)

which has two solutions, that we call $s_{\rm frag}^{\pm}$, expressed as (see
Appendix A of Gonzalez et al., 2017a)

$s_{\rm frag}^{\pm}=s_{\rm opt}\left(\dfrac{\tilde{\alpha}c_{\rm
s}^{2}}{V_{\rm frag}^{2}}-1\pm\dfrac{\sqrt{\tilde{\alpha}}c_{\rm s}}{V_{\rm
frag}}\sqrt{\dfrac{\tilde{\alpha}c_{\rm s}^{2}}{V_{\rm frag}^{2}}-2}\right),$
(18)

where ‘+’ and ‘-’ refer to the cases where $\mathrm{St}>1$ and $\mathrm{St}<1$
respectively. One should note that Equation (18) is not defined when
$\tilde{\alpha}c_{\rm s}^{2}<2V_{\rm frag}^{2}$. This corresponds to the case
where the maximum relative velocity is smaller than the fragmentation
velocity, which means that fragmentation never occurs. Using this inequality
with a power-law prescription for the temperature, we find a corresponding
fragmentation radius $r_{\rm frag}$

$r_{\rm frag}=r_{0}\left(\dfrac{\tilde{\alpha}}{2}\dfrac{c_{\rm
s,0}^{2}}{V_{\rm frag}^{2}}\right)^{\frac{1}{q}},$ (19)

where $q$ is the temperature power-law index and $r_{0}$ is a reference
radius. Beyond $r_{\rm frag}$, fragmentation is not possible.

##  3 Implementation

The dust growth model is introduced in both Phantom’s two-fluid and one-fluid
algorithms (Laibe & Price, 2012a, b, 2014a, 2014b, 2014c; Price & Laibe, 2015;
Ballabio et al., 2018; Price & Laibe, 2020). Firstly, we specify the
implementation with the two-fluid algorithm. In that configuration, the
accuracy of the implementation relies on the evaluation of a few gas-related
quantities on each dust particle. This is needed in particular for $\Delta v$,
$\mathrm{St}$ and $c_{\mathrm{s}}$, which are not defined for a given dust
particle but rather for a gas-dust pair or for gas only.

###  3.1 Two fluid means two smoothing lengths

The SPH method allows one to evaluate any given quantity by means of
interpolations over the particle neighbours. Modelling dust and gas as two
separate fluids comes at the price of increasing the total number of particles
and having to handle two sets of resolution lengths. For a gas or dust
quantity (e.g the density), the typical length over which the neighbours are
used is the gas or dust smoothing length (more precisely a few smoothing
lengths, depending on the kernel), which are expressed as

$h_{\mathrm{gas}}=h_{\mathrm{0}}\left(\dfrac{m_{a}}{\rho_{a}}\right)^{1/3},$
(20)

where $h_{0}$ is a constant setting the average number of neighbours, $m_{a}$
is the mass of a gas particle and $\rho_{a}$ its volume density, and

$h_{\mathrm{dust}}=h_{\mathrm{0}}\left(\dfrac{m_{i}}{\rho_{i}}\right)^{1/3},$
(21)

where $m_{i}$ is the mass of a dust particle and $\rho_{i}$ its volume
density. We used the index $i$ for dust particles and $a$ for gas particles,
as in e.g. Monaghan (1997); Laibe & Price (2012a, b).

As gas and dust do not necessarily have the same behaviour, one needs to be
careful when dealing with two sets of smoothing lengths, especially if an
interpolated quantity for a particle type requires a loop over the other one.
This is typically the case for the drag between gas and dust, and in our case
additionally for $\Delta v$, $\mathrm{St}$ and $c_{\mathrm{s}}$. The risk with
such a method is that if the dust-to-gas ratio increases significantly, the
dust smoothing length can become much smaller than that of the gas. This can
jeopardise further interpolations due to a lack of neighbours found inside the
corresponding smoothing length. Secondly, as gas is pressure-supported but
dust is not, reducing the dust smoothing length under that of the gas can also
locally clump close dust particles together as they would not feel any
repulsive force from the gas. To limit both of these risks, we follow Laibe &
Price (2012a) and use in our cross-species interpolations the maximum of the
gas and dust smoothing lengths, which is slightly more computationally
expensive but safer. We denote this maximum smoothing length
$h_{ia}=\mathrm{max}(h_{i},h_{a})$ for future use.

###  3.2 Sound speed and gas density

The gas density is interpolated using the kernel-weighted sum over the
neighbours

$\rho_{\mathrm{g},i}=\sum_{a}m_{a}W(\bm{r}_{ia},h_{ia}),$ (22)

where $W$ is the kernel function (quintic M6 in our case). The other
quantities are then interpolated using $\rho_{\mathrm{g},i}$, that is for a
function $g$

$\displaystyle g_{i}$
$\displaystyle=\dfrac{\sum_{a}m_{a}g_{a}W(\bm{r}_{ia},h_{ia})}{\sum_{a}m_{a}W(\bm{r}_{ia},h_{ia})},$
(23)
$\displaystyle=\dfrac{1}{\rho_{\mathrm{g},i}}\sum_{a}m_{a}g_{a}W(\bm{r}_{ia},h_{ia}).$
(24)

We found that this form, slightly different than the usual interpolation
(weighted by $m$ instead of $m/\rho$), was generally more accurate than its
counterpart, especially in regions where the dust concentrates. This choice of
implementation does not affect the speed of execution of the code.

The sound speed on each dust particle is interpolated using

$c_{\mathrm{s},i}=\dfrac{1}{\rho_{\mathrm{g},i}}\sum_{a}m_{a}c_{\mathrm{s,a}}W(\bm{r}_{ia},h_{ia}).$
(25)

###  3.3 Differential velocity

We evaluate each component of the differential velocity vector before
computing the norm

$\Delta
v_{x,i}=\dfrac{1}{\rho_{\mathrm{g},i}}\sum_{a}m_{a}\Delta{v}_{x,ia}W(\bm{r}_{ia},h_{ia}),$
(26) $\Delta
v_{y,i}=\dfrac{1}{\rho_{\mathrm{g},i}}\sum_{a}m_{a}\Delta{v}_{y,ia}W(\bm{r}_{ia},h_{ia}),$
(27) $\Delta
v_{z,i}=\dfrac{1}{\rho_{\mathrm{g},i}}\sum_{a}m_{a}\Delta{v}_{z,ia}W(\bm{r}_{ia},h_{ia}),$
(28) $||\Delta\bm{v}_{i}||=\Delta v_{i}=\sqrt{\Delta v_{x,i}^{2}+\Delta
v_{y,i}^{2}+\Delta v_{z,i}^{2}}.$ (29)

One should notice that $\Delta v$ is not evaluated by the double-humped kernel
$D$ used in Laibe & Price (2012a) but rather by the usual density kernel $W$.
Indeed, here we do not need to evaluate $\Delta\bm{v}$ along the line of sight
of a gas-dust pair, thus the scalar product
$\Delta\bm{v}_{ia}\cdot\bm{r}_{ia}$ is avoided and the evaluation is more
accurate with the usual bell-shaped kernel $W$.

###  3.4 Stokes number

The evaluation of the Stokes number depends on the drag regime that the grains
experience. In Phantom, the Epstein and Stokes regimes are automatically
selected depending on the Knudsen number
$K_{\mathrm{n}}=9\lambda_{\mathrm{g}}/4s$ (Weidenschilling, 1977; Stepinski &
Valageas, 1996). To follow this and remain as general as possible in our
implementation, we evaluate the local Stokes number (which corresponds to the
mixture) by computing the stopping time using $\rho_{\mathrm{g},i}$,
$c_{\mathrm{s},i}$ and $\Delta v_{i}$ rather than interpolating the stopping
time between each gas-dust pair

$\mathrm{St}_{i}=t_{\mathrm{s}}\left(\rho_{\mathrm{g},i},c_{\mathrm{s},i},\Delta
v_{i}\right)\Omega_{\mathrm{k},i}.$ (30)

The stopping time $t_{\mathrm{s}}$ is expressed as (Laibe & Price, 2012a)

$t_{\mathrm{s}}=\dfrac{\rho_{\mathrm{g}}\rho_{\mathrm{d}}}{K(\rho_{\mathrm{g}}+\rho_{\mathrm{d}})},$
(31)

where $K$ is the drag coefficient. Considering the Epstein drag regime again,
this gives for example (P18)

$\mathrm{St}_{i}=\sqrt{\dfrac{\pi\gamma}{8}}\dfrac{\rho_{\mathrm{s}}s_{i}}{f_{i}\left(\rho_{\mathrm{g},i}+\rho_{i}\right)c_{\mathrm{s},i}}\Omega_{\mathrm{k},i}.$
(32)

###  3.5 One-fluid implementation

To simulate the most strongly coupled populations of grains, Phantom also has
a ‘one-fluid’ algorithm that represents the gas-dust mixture with only one set
of particles and thus one resolution (Laibe & Price, 2014a, b, c; Price &
Laibe, 2015; Ballabio et al., 2018) We also implemented the growth model with
the one-fluid formalism, although its validity is restricted to relatively
small Stokes numbers in order to respect the terminal velocity approximation
(Youdin & Goodman, 2005). With that in mind, the one-fluid model with dust
growth can be used for a short amount of time and mainly to accelerate the
evolution of tightly coupled grains. We also added to the code a conversion
tool that allows the user to convert to the two-fluid method when the one-
fluid regime reaches its limit.

The implementation itself is straightforward in one-fluid since the cross
species interpolations vanish. To compute the values of interest, we simply
use the values stored onto the mixture particles and transform them out of the
barycentric frame, that is applying the changes of variables

$\displaystyle\rho_{\rm g}$ $\displaystyle=(1-\epsilon)\rho,$ (33)
$\displaystyle\rho_{\rm d}$ $\displaystyle=\epsilon\rho,$ (34)
$\displaystyle{\bm{v}}_{\rm g}$ $\displaystyle=\bm{v}-\epsilon\Delta\bm{v},$
(35) $\displaystyle{\bm{v}}_{\rm d}$
$\displaystyle=\bm{v}+(1-\epsilon)\Delta\bm{v},$ (36)

where $\rho=\rho_{\rm g}+\rho_{\rm d}$ and $\epsilon=\rho_{\rm d}/\rho$ is the
dust fraction.

###  3.6 Timestepping

####  3.6.1 Numerical scheme

The size evaluation by the integrator follows the same predictor-corrector
scheme that is described in P18.

$\displaystyle s^{n+\frac{1}{2}}$ $\displaystyle=s^{n}+\dfrac{1}{2}\Delta
t\left(\dfrac{\mathrm{d}s}{\mathrm{d}t}\right)^{n},$ (37) $\displaystyle
s^{*}$ $\displaystyle=s^{n+\frac{1}{2}}+\dfrac{1}{2}\Delta
t\left(\dfrac{\mathrm{d}s}{\mathrm{d}t}\right)^{n},$ (38)
$\displaystyle\left(\dfrac{\mathrm{d}s}{\mathrm{d}t}\right)^{n+1}$
$\displaystyle=\dfrac{\mathrm{d}s}{\mathrm{d}t}\left(s^{*}\right),$ (39)
$\displaystyle s^{n+1}$ $\displaystyle=s^{n+\frac{1}{2}}+\dfrac{1}{2}\Delta
t\left(\dfrac{\mathrm{d}s}{\mathrm{d}t}\right)^{n+1}.$ (40)

The size integration is also implemented with individual timestepping, which
can be significantly faster if dust grain sizes span a wide range.

####  3.6.2 Constraints

The timestepping for growth and fragmentation is constrained by the Courant-
Friedrichs-Lewy (Courant et al., 1928, CFL) condition, which requires the
timestep to be smaller than the typical growth timescale

$\Delta t\leq C_{\mathrm{cour}}\tau_{\mathrm{g}},$ (41)

where $C_{\mathrm{cour}}$ is a constant typically of order unity and
$\tau_{\mathrm{g}}$ is the growth time scale

$\tau_{\mathrm{g}}=\dfrac{s}{|{\mathrm{d}s}/{\mathrm{d}t}|}.$ (42)

Considering that $\Delta v\ll V_{\mathrm{t}}$ and $c_{\rm s}$ (such that
$\mathrm{Sc}\simeq 1+\mathrm{St}$ and $f\simeq 1$), this gives

$\dfrac{\tau_{\mathrm{g}}}{t_{\mathrm{s}}}=\dfrac{1+\varepsilon}{\varepsilon}\sqrt{\dfrac{8}{\pi\gamma}}\dfrac{\left(1+\mathrm{St}\right)}{\sqrt{2^{3/2}\mathrm{Ro}\alpha\mathrm{St}}}.$
(43)

Figure 3: Ratio $\tau_{\rm g}/t_{\rm s}$ (Equation 43) as a function of the
Stokes number for different values of the dust-to-gas ratio.

As seen in Fig. 3, the minimum timestep required to respect the CFL condition
is much larger than the stopping time. This condition is already satisfied in
Phantom, since the timestep is always at most equal to the stopping time. This
also results in very small timesteps for tightly coupled grains, i.e. the
smallest ones, which can drastically slow down the simulations. To limit this
effect when fragmentation is involved, we incorporate a minimum allowed grain
size, under which dust can not fragment. We check and adjust the grain size
both during the predictor and corrector steps.

##  4 Benchmark tests

Phantom is equipped with a set of tests that are performed frequently to hunt
down potential bugs. Implementing a dust growth algorithm comes with a few
more tests that can be added to these checkups. In this section, we
(re-)present two dust growth-related tests: the existing Dustybox and the new
Farmingbox.

###  4.1 The return of the Dustybox

Dustybox is a simple two-fluid test in which gas and dust have an initial
differential velocity in a uniform box. This test ensures that the gas drag
suppresses the differential velocity between the two phases according to an
analytical solution that is relatively simple (the drag coefficient is
constant). This test has been extensively described in Laibe & Price (2012a)
and we refer the reader to this paper for more details about the setup. We
extended this test by adding the verification of the interpolated differential
velocity (Equation 29). With a constant drag coefficient $K$, the differential
velocity between gas and dust decreases exponentially in time, with a
characteristic timescale that is the stopping time:

$\Delta v=\Delta v_{0}e^{-t/t_{\mathrm{s}}}=\Delta v_{0}e^{-2Kt},$ (44)

where $\Delta v_{0}$ is the initial differential velocity and the factor 2
comes from the fact that $\rho_{\mathrm{g}}=\rho_{\mathrm{d}}$ in the test
(see Equation 31).
We present a few examples of this test with different values of the drag
coefficient $K$ in Fig. 4.

Figure 4: Differential velocity as a function of time for 3 different values
of the drag coefficient $K$ in Dustybox. The solid lines represent the
analytical solutions while the dots represent interpolated values during the
test. All the values are in code units throughout the test.

As the drag coefficient increases (from blue to green), the grains are more
coupled to the gas. As a result the initial differential velocity is more
efficiently damped and can reach a steady state (i.e. $\Delta v=0$) faster.
The differential velocities of all the particles are tested at every timestep
and the results show maximum relative errors of about $10^{-4}$ with the
analytical solution.

###  4.2 Farmingbox: motionless dust growth in the Epstein regime

The main goal of Farmingbox is to simulate the motionless growth or
fragmentation of dust particles and compare the integrated size with its
analytical solution. To attain an analytically solvable growth rate, we fix
the particles at their initial position (which cancels $\Delta v$). Similarly
to Dustybox, we set $\sim 7000$ gas and dust particles on a 3D lattice inside
of a 1 $\times$ 0.5 $\times$ 0.3 au3 box. Gas and dust are set so that the
dust experiences the Epstein regime throughout its growth or fragmentation. We
use this test in both pure growth and pure Hard fragmentation models. The
physical parameters used for each mode are shown in Table 1. We only present
the results of the two-fluid algorithm in this Section, although the one-fluid
version was also tested and gave sensible results.

Table 1: List of parameters used in Farmingbox for both pure growth and pure fragmentation modes. All parameters are uniform throughout the box. Parameter | Unit | Growth | Fragmentation
---|---|---|---
$c_{\mathrm{s}}$ | m.s-1 | $942$ | $942$
$\rho_{\mathrm{g}}$ | kg.m-3 | $10^{-8}$ | $10^{-8}$
$\rho_{\mathrm{d}}$ | kg.m-3 | $10^{-8}$ | $5\times 10^{-9}$
$\rho_{\mathrm{s}}$ | kg.m-3 | $1000$ | $1000$
$s_{0}$ | m | $10^{-4}$ | $10^{-2}$
$\alpha$ | $\varnothing$ | $10^{-2}$ | $2.5\times 10^{-2}$

####  4.2.1 Growth

Firstly, we test the pure growth case. Gas and dust particles being
motionless, the Schmidt number is simply $\mathrm{Sc}=1+\mathrm{St}$.
Following L08, we can rewrite Equation (13) with a change of variables such
that

$T=\dfrac{t}{\tau}+2\sqrt{\mathrm{St}_{0}}\left(1+\dfrac{\mathrm{St}_{0}}{3}\right),$
(45)

where $\mathrm{St}_{0}$ is the initial Stokes number of the considered
particle, of initial size $s_{0}$, and $\tau$ is a timescale defined as

$\tau=\sqrt{\dfrac{8}{\pi\gamma}}\dfrac{1}{\sqrt{2^{3/2}\mathrm{Ro}\alpha}\Omega_{\mathrm{k}}}\dfrac{\rho}{\rho_{\mathrm{d}}}.$
(46)

The Keplerian frequency $\Omega_{k}$ is set as $r^{-3/2}$ in code units, which
emulates a $1$ $M_{\odot}$ central star. Even though there is no such star in
the test, we keep this functional form to maintain the property that growth is
different for particles at different distances to the center of the box. This
is also appropriate in order to have values of sizes and Stokes numbers that
would be found in real simulations.

With these new variables we can rewrite Equation (13) in terms of Stokes
number and dimensionless time

$\dfrac{\mathrm{dSt}}{\mathrm{d}T}=\dfrac{\sqrt{\mathrm{St}}}{1+\mathrm{St}}.$
(47)

The solution to this equation is (L08)

$\mathrm{St}=\dfrac{\sigma}{2}+\dfrac{2}{\sigma}-2,$ (48)

where $\sigma=\left(8+9T^{2}+3T\sqrt{16+9T^{2}}\right)^{1/3}$. Note that this
solution has been adapted from L08 with our definition of the Stokes number
(Equation 46).

Finally, this gives the time evolution of the grain size

$s=\sqrt{\dfrac{8}{\pi\gamma}}\dfrac{\rho
c_{\mathrm{s}}}{\rho_{\mathrm{s}}\Omega_{k}}\left(\dfrac{\sigma^{2}-4\sigma+4}{2\sigma}\right)=s_{0}\left(\dfrac{\sigma^{2}-4\sigma+4}{2\mathrm{St}_{0}\sigma}\right).$
(49)

Figure 5: Dust size (top) and Stokes number (bottom) analytical distributions
as a function of the distance to the center of the Farmingbox for 6 different
times in the pure growth case. Three specific particles are plotted on top
(squares, triangles and stars). The error in L1 norm is approximately
$6.13\times 10^{-3}$.

Figure 6: Size (top) and Stokes number (bottom) evolution of the 3 particles
plotted in Fig. 5 as a function of time. The lines represent the exact
solutions while the dots are the values outputted by the code.

Figs. 5 and 6 show the results of Farmingbox in pure growth mode. Initially,
all the particles have the same size (top panel of Fig. 5, black line), which
results in the Stokes number profile to be the same as the Keplerian frequency
(bottom panel, black line) and spanning between 2 and 3 orders of magnitude
from $10^{-3}$ and up. As growth starts, particles at different radii grow
differently due to their different Stokes number. More specifically, the
growth rate is maximum for Stokes unity and decreases when grains are both
very coupled or decoupled from the gas. As a result, the final size
distribution (top panel, red curve) shows a maximum at $r\sim 0.1$ au,
corresponding to the location where the grains have spent an optimal amount of
time around $\mathrm{St}=1$. Fig. 6 demonstrates that in more details, as
particle 1 (bottom panel, blue dots) crosses $\mathrm{St}=1$ rapidly and
decreases its growth rate (top panel). Conversely, particle 2 (red dots)
starts with a lower Stokes number but spends most of its time close to Stokes
unity, which overall is more efficient at making the particle grow. Finally,
particle 3 (green dots) spends the entire test at low Stokes numbers due to
its larger distance from the center of the box, which is growth inefficient.

####  4.2.2 Fragmentation

We now consider pure fragmentation. Similarly to the pure growth case, we
introduce the same change of variable but with the time going in the opposite
direction

$T=-\dfrac{t}{\tau}+2\sqrt{\mathrm{St}_{0}}\left(1+\dfrac{\mathrm{St}_{0}}{3}\right),$
(50)

which also gives Equation (47) and the subsequent solution (equations 48 and
49). In this particular case of pure fragmentation, $T$ can reach negative
values, which would correspond to a negative size, which is not possible. As
this is a physical limit, we ensure that $T$ stays positive throughout the
test, i.e. we set the physical parameters in the box so that (see Table 1)

$t_{\mathrm{max}}<2\sqrt{\mathrm{St}_{0}}\left(1+\dfrac{\mathrm{St}_{0}}{3}\right)\tau,$
(51)

where $t_{\mathrm{max}}$ is the maximum duration of the test, i.e. the time
after which $T<0$.

Figure 7: Similar to Fig. 5 but for pure fragmentation.

Figure 8: Similar to Fig. 6 but for pure fragmentation. The error in L1 norm
is approximately $1.16\times 10^{-3}$.

Figs. 7 and 8 show the results of Farmingbox in pure fragmentation mode. The
behaviour is very similar to the growth case, in the sense that fragmentation
is also the most efficient near Stokes unity, which is again observed for
grains at $r\sim 0.1$ au.

During Farmingbox, the Stokes number, size, and sound speed of all dust
particles are tested against the analytical values at each timesteps. Fig. 5
to 8 only display a subset of these data points to keep the plots readable.
The relative errors between the quantities and their analytical values
throughout Farmingbox stay in the range of $10^{-3}-10^{-4}$.

##  5 Circumstellar disc simulations

In this Section, we perform several 3D circumstellar disc simulations to
compare the results given by Phantom with previous studies. In that regard, we
will simulate a single disc model, that we will refer to as ‘Standard’ (Std)
and that is similar to the one used in previous papers (e.g. Barrière-Fouchet
et al., 2005; Laibe et al., 2008). We only consider the two-fluid algorithm
from now on.

###  5.1 Numerical setup

We model a $1$ $M_{\odot}$ star around which orbits a $0.01$ $M_{\odot}$ disc
represented by $1.2\times 10^{6}$ particles: 1M gas and 0.2M dust. This ratio
has been chosen to resolve the gas scale height and to limit dust over-
concentrations (see Section 6.2). The gas particles are positioned between
$r_{\rm in}$ and $r_{\rm out}$ (see Table 2), whereas the dust particles only
extend up to $5/6\,r_{\rm out}$. We found having a dust disc smaller than the
gas disc to be quite important in order to avoid numerical artefacts linked to
the gas relaxation in the outer disc. This numerical choice is also justified
by the radial drift, which inevitably tends to reduce the size of the dust
disc. To assist the disc relaxation in the inner parts of the disc, the gas
initially follows a smoothed surface density profile such that

$\Sigma_{\mathrm{g}}(r)=\Sigma_{0}\left(1-\sqrt{\dfrac{r_{\mathrm{in}}}{r}}\right)\left(\dfrac{r}{r_{0}}\right)^{-p},$
(52)

where $r$ refers to the cylindrical radius and $\Sigma_{0}$ is set by the
total disc mass and is then allowed to evolve. The dust disc initially follows
the same law, with an initial and uniform dust-to-gas ratio of 1%. Both gas
and dust particles are removed from the simulation if they cross an accretion
radius $r_{\mathrm{acc}}$ that is set to be equal to $r_{\mathrm{in}}$ given
our interest for the regions where $r>r_{\rm in}$. The gas temperature is
supposed locally isothermal and follows the usual power-law prescription

$T(r)=T_{0}\left(\dfrac{R}{r_{0}}\right)^{-q},$ (53)

where $R$ refers to the spherical radius and $T_{0}$ is set by the disc aspect
ratio at $r_{0}$ (see Table. 2). In order to represent a Shakura & Sunyaev
(1973) viscosity parameter $\alpha_{\rm SS}=5\times 10^{-3}$, we follow the
formalism of Lodato & Price (2010), in which the SPH artificial viscosity
parameter can be related to $\alpha_{\rm SS}$ as

$\alpha^{\rm AV}\sim\dfrac{\alpha_{\rm SS}}{10}\dfrac{\langle h\rangle}{H},$
(54)

where $<h>$ is the average smoothing length on a given portion of the disc,
with $\alpha^{\rm AV}\sim 0.2$. We take the coefficient preventing particle
interpenetration $\beta^{\rm AV}=2$.

The dust initially has a uniform size $s_{0}$ of $30$ $\mu$m throughout the
disc and is allowed to grow or fragment in agreement with the model presented
in Section 2. The influence of the initial grain size is explained in Appendix
A. Considering the low disc mass, we neglect the disc self-gravity but do take
into account the dust back-reaction onto the gas.

Table 2: The disc model used in our simulations, with $r_{0}=1$ $\mathrm{au}$. Name | $p$ | $q$ | $M_{\mathrm{disc}}$ | $H/r_{0}$ | $r_{\mathrm{in}}$ | $r_{\mathrm{out}}$
---|---|---|---|---|---|---
Std | $3/2$ | $3/4$ |  $0.01$ $M_{\odot}$ | $0.0281$ |  $20$ au |  $300$ au

The short list of circumstellar disc simulations performed for this paper is
showed in Table 3.

Table 3: Simulation suite. Label | Fragmentation |  $V_{\rm frag}$ [m.s-1] | Back-reaction
---|---|---|---
noF | No | $\varnothing$ | Yes
F-Hard | Hard | 15 | Yes
F-Hard-noBR | Hard | 15 | No
F-Smooth | Smooth | 15 | Yes

###  5.2 Pure growth

We first consider and analyse a simulation without fragmentation (noF). Pure
dust growth in a similar disc and with a similar growth model has been studied
by L08, therefore our results will be compared to theirs.

####  5.2.1 Radial drift and dust concentration

Figure 9: Left: Pressure radial profile at 5 different times between 1 and 45
kyr in simulation noF. The position of the pressure maximum is showed with the
vertical black dotted line. Right: Dust (solid lines) and gas (dashed lines)
surface density radial profiles for the same times as the left plot. The
position of the pressure maximum is also shown with the vertical black dotted
line at $r\sim 28$ au.

Fig. 9 shows the pressure and surface density radial profiles at different
times. The gas reaches a steady state in a few thousand years, resulting in a
pressure maximum located near 28 au (left panel), which corresponds to the
inner edge of the disc. The negative pressure gradient between $\sim 28$ au
and the outermost parts of the disc induces a headwind on the dust that leads
to its radial drift towards the star. As a result, the dust mass is
transferred towards the pressure maximum (right panel), where the headwind
cancels itself out. At this location the dust concentrates, while the outer
parts of the disc are slowly being emptied. The dust radial drift also reduces
the extent of the dust disc as one can see from the dust surface density
profile (right panel). This whole process is accelerated by dust growth, which
increases radial drift velocities up to the point where $\mathrm{St=1}$.

Figure 10: Trajectories in the ($r,s$) plane of 6 particles initially in the
mid-plane of the disc (solid lines) between 0 and 50 kyr in simulation noF.
Empty circles of the same colour mark each particle’s position every 10 kyr,
with a filled circle for the final state. They grey points in the background
represent the dust distribution at 50 kyr. The black dashed line represents
$s_{\rm opt}$ (Equation 8) at $t=0$ ($\varepsilon=0.01$), i.e. where
$\mathrm{St}=1$.

To study the coupled effects of dust growth and radial drift, we follow in
Fig. 10 the evolution of 6 particles initially in the mid-plane of the disc
with a size $s_{0}=30$ $\mu$m. Their initial distance to the star are 40, 60,
80, 100, 150 and 200 au respectively. We also over-plotted an estimation of
$s_{\rm opt}$ at $t=0$, for which we considered the dust-to-gas ratio to be
uniform and the differential velocity to be subsonic (such that $f\sim 1$, see
Equation 6). While this estimation is limited near local dust concentrations,
it still provides a useful comparative basis for the rest of the simulation.
Dust particles shown in Fig. 10 have similar behaviours that can be decomposed
into 3 steps:

  1. 1.

A growth phase with very little radial drift: $\mathrm{St}\ll 1$.

  2. 2.

A phase of radial drift and growth where radial drift is more efficient:
$\mathrm{St\sim 1}$.

  3. 3.

A slower growth without radial drift: $\mathrm{St}\gg 1$.

The further away the particles are, the longer it takes for them to experience
all three phases. For example, after 50 kyr, particles P1 to P3 experienced
these different steps, whereas P4 is on the verge of entering phase (iii).
Particles P6 is only entering phase (ii).

These different steps shape the dust distribution as we can see on the
background (grey points, shown at 50 kyr). Grains in phase (i) form a tail at
$r\gtrsim 200$ au, while grains in phase (ii) form a diagonal between $\sim
40$ and $\sim 200$ au. Particles that are decoupled (phase (iii)) form a
reservoir of material at the pressure maximum, that is at $\sim 28$ au. The
three specified evolutionary stages are the same as those identified in L08.

The diagonal formed by drifting and growing grains (phase (ii)) is steeper
than the curve corresponding to $s_{\rm opt}$. This result, also found by L08,
shows the effect of dust growth during the radial drift process, which
increases their Stokes number faster than radial drift decreases it. This
behaviour was also found analytically by Laibe (2014).

####  5.2.2 Vertical settling

Another important aspect of dust growth that we want to explore is its effect
on vertical settling. To this end, Fig. 11 shows the evolution of the disc
scale height. Pure growth is very efficient at transforming the initially
flared dust disc into a thin stratified layer in the mid-plane. While the gas
disc has a roughly constant aspect ratio of about 0.05 throughout the
simulation, the dust disc aspect ratio tends to be less than $10^{-3}$ in a
few tens of thousands of years. This process is efficient because dust
decouples from the gas as it grows. More precisely, the aerodynamic drag
efficiently damps the vertical oscillations of the dust particles around the
mid-plane, which tends to form layers with a thickness limited by the disc
turbulence (Dubrulle et al., 1995; Fromang & Papaloizou, 2006). These studies
predict that the dust layer thickness is roughly proportional to
$\mathrm{St^{-1/2}}$. This is consistent with what we find, since the average
Stokes number decreases when the distance to the star increases (see Fig. 10,
particles above $s_{\rm opt}$ are in the inner disc in contrast to the
particles in the outer disc).

Figure 11: Scale height radial profiles of the gas (black dashed line) and
dust (coloured solid lines) discs at 6 different times between 0 and 45 kyr in
simulation noF.

We take advantage of the Lagrangian formalism and follow the vertical settling
process of a few particles in more detail. We select 5 particles, noted P7 to
P11, with an initial distance to the star of 50, 70, 100, 150 and 200 au
respectively. Fig. 12 plots their evolution in altitude and Stokes number.
Particle trajectories in the ($r,z$) plane are first vertical, and then
radial, allowing us to consider both movements as decoupled from each other.
The settling timescale is smaller than the drift timescale by a factor of the
order of $(H/r)^{2}\partial\ln P/\partial\ln r\ll 1$ (Garaud et al., 2004;
Barrière-Fouchet et al., 2005). The evolution of their Stokes number is
another interesting point, as it reveals two opposing processes during the
settling phase:

  1. 1.

The settling of grains to denser layers of gas and dust, which strengthens the
coupling between the two ($\mathrm{St}\propto(\rho_{\rm d}+\rho_{\rm
g})^{-1}\searrow$).

  2. 2.

Dust growth that tends to decouple the dust from the gas ($\mathrm{St}\propto
s\nearrow$).

For P8, P9 and P11, this competition is clear. The settling and growth
dominate the Stokes number evolution one after the other. For P7, its
relatively close proximity to the star and to the mid-plane makes the growth
process very efficient and therefore it always dominates.

Figure 12: Left: Trajectories in the ($r,z$) plane of 5 particles between 0
and 50 kyr in simulation noF. The filled squares, empty squares and filled
circles represent the initial time, the time of fist contact with the mid-
plane and the final time respectively. Right: Same as the left panel but for
the ($r,\mathrm{St}$) plane.

The pure growth simulation has shown the coupled effect of dust growth,
vertical settling and radial drift. Our results highlight the exact same dust
regimes of evolution as those described in L08.

Although pure growth is an interesting academic case in order to understand
the interplay between these different processes, dust growth in discs is more
realistically limited by fragmentation when the collision velocities become
larger than a given threshold. In the next section, we perform more realistic
simulations by taking fragmentation into account.

###  5.3 Adding fragmentation

Considering the same disc, we add fragmentation as presented in Section 2. We
firstly analyse the simulation with the Hard fragmentation model (F-Hard). We
consider a fragmentation threshold of 15 m.s-1, value often used for wet icy
aggregates (e.g. Wada et al., 2009; Gonzalez et al., 2015; Vericel & Gonzalez,
2020). Simulations with fragmentation were evolved twice as long as those with
pure growth — that is $\sim 100$ kyr — because they reach a quasi steady state
more slowly than their counterpart.

####  5.3.1 General evolution

The evolution of the dust size radial distribution is firstly displayed in
Fig. 13, colour-coded with $V_{\mathrm{rel}}/V_{\mathrm{frag}}$. The
fragmentation sizes $s_{\rm frag}^{\pm}$ (Equation 18) are also over-plotted:
they are estimated using the form of $s_{\rm opt}$ computed in Section 5.2
with the power laws describing the initial disc structure and thus do not
reflect the evolution of the gas profile. The curves corresponding to $s_{\rm
frag}^{\pm}$ define a ‘fragmentation zone’ that is hatched in Fig. 13. The
fragmentation zone is shown using increasing dust-to-gas ratios over time,
which are estimated using the surface density profiles. Notice that the
fragmentation zone extends out to $\sim 100$ au, which fits with the
calculation of $r_{\rm frag}$ (Equation 19) that gives $\sim 105$ au.

Figure 13: Dust size radial distribution between 0 and 90 kyr in simulation
F-Hard. Particles are coloured with their ratio $V_{\rm rel}/V_{\rm frag}$.
The hatched zone corresponds to the fragmentation zone, which is delimited by
the fragmentation sizes $s_{\rm frag}^{\pm}$. These sizes are evaluated with
increasing dust-to-gas ratios over time (see legend on each panel).

Dust at $50\lesssim r\lesssim 100$ au reaches fragmentation sizes $s_{\rm
frag}^{-}$ relatively rapidly ($\sim 20$ kyr, third panel). Meanwhile, dust in
the innermost parts of the disc is incapable of growing, which creates a steep
size gradient at $r\sim 50$ au ($40-50$ kyr). This size contrast induces a
gradient of radial drift velocities, i.e. a traffic jam. In other words, dust
concentrates radially near the gradient of sizes, which in the figure can be
seen from 40 kyr onward. A second – but slight – size gradient can also be
observed at the entry of the fragmentation zone ($\sim 100$ au). Nevertheless,
this gradient is less important and does not create a significant traffic jam
in the disc. Note that once dust grains are trapped at $r\sim 50$ au, they
continue to grow even when they appear to be in the hatched area. This is
because the latter represents the initial shape of the fragmentation zone, not
reflecting the changes in the gas profile or the reduction of the relative
velocities of grains due to their confinement in the trap. Finally, the
runaway growth seen after $\sim 80$ kyr results from artificial clumping,
which locally overestimates the growth rate (see Section 6.2.)

To understand in more detail the process taking place in this simulation, we
show the evolution of the pressure and surface density profiles in Fig. 14.
The dust concentration is clearly observable in the surface density profiles
at $\sim 50$ au. Additionally, a pressure maximum develops near the dust
concentration after a few thousand years, which hints at the back-reaction
onto the gas as the source of this feature.

The combination of dust growth/fragmentation, its pile-up and the formation of
a local pressure maximum in the disc might be the manifestation of the ‘self-
induced dust trap’ mechanism proposed by Gonzalez et al. (2017a). The next
section focuses on this hypothesis.

Figure 14: Left: Pressure radial profiles between 0 and 90 kyr in simulation
F-Hard. Right: Same, but for the dust (solid lines) and gas (dashed lines)
surface densities.

####  5.3.2 Self-induced dust trap: the importance of the dust back-reaction

The ‘self-induced dust trap’ mechanism has been proposed as a means for dust
to overcome both the radial drift and fragmentation barriers in discs: by
taking into account dust growth and fragmentation as well as large-scale
gradients, dust particles pile-up and concentrate in the disc, such that the
back-reaction reverses the gas flow locally to create a pressure maximum that
further traps the dust.

To estimate the importance of the back-reaction, Gonzalez et al. (2017a)
considered the stationary solution for the gas radial velocity, which gives
(Dipierro & Laibe, 2017; Kanagawa et al., 2017; Gonzalez et al., 2017a)

$v_{\mathrm{g},r}=-\dfrac{1}{1+\varepsilon}\left[\dfrac{\varepsilon\mathrm{St}}{1+\mathrm{St}^{2}}v_{\rm
drift}-\left(1+\dfrac{\varepsilon\mathrm{St}^{2}}{1+\mathrm{St}^{2}}\right)v_{\rm
visc}\right],$ (55)

where

$v_{\rm drift}=\left(\dfrac{H}{r}\right)^{2}\dfrac{\partial\ln P}{\partial\ln
r}v_{\rm k}$ (56)

corresponds to the drift velocity caused by the pressure gradient (Nakagawa et
al., 1986) and

$v_{\mathrm{visc}}=\dfrac{\dfrac{\partial}{\partial
r}\left(\rho_{\mathrm{g}}\nu r^{3}\dfrac{\partial\Omega_{\mathrm{k}}}{\partial
r}\right)}{r\rho_{\mathrm{g}}\dfrac{\partial}{\partial
r}(r^{2}\Omega_{\mathrm{k}})}$ (57)

is the usual viscosity-induced velocity (Lynden-Bell & Pringle, 1974). Notice
that the form used by Gonzalez et al. (2017a) is slightly different, because
their Stokes number definition concerns the dust only and not the mixture,
which differs by a factor $1+\varepsilon$. The viscosity term tends to push
the gas inward, whereas the drag term does the opposite. If the drag term is
greater than its counterpart, the net motion is towards the outer parts of the
disc, which changes the gas structure locally.

To measure the impact of back-reaction on the gas motion, we use a parameter,
called $x_{\rm br}$, which measures the balance between the viscous and drag
terms in equation (55). Following Gonzalez et al. (2017a, b), i.e. assuming
$v_{\rm visc}/v_{\rm drift}\sim 1/\alpha$, we subsequently find

$x_{\rm
br}\simeq\dfrac{\alpha^{-1}\varepsilon\mathrm{St}}{1+\mathrm{St}^{2}(1+\varepsilon)}.$
(58)

Where $x_{\rm br}>1$, $v_{\mathrm{g},r}>0$ and the gas flows outwards, which
corresponds to the formation of a local pressure maximum. The gas mass flux as
well as the estimation of $x_{\rm br}$ are plotted in Fig. 15. When the dust
concentration develops ($\sim 30$ kyr), the gas starts to be pulled outwards
($r\Sigma_{\mathrm{g}}v_{\mathrm{g},r}>0$, cyan). Moreover, we observe a
tendency where the gas seems to be pulled further away with increasing time.
This is consistent with the nature of the self-induced dust trap mechanism,
because the dust concentration extends at increasing distances to the star
over time. Indeed, the regions where $x_{\mathrm{br}}\geq 1$ becomes more
extended as time goes on. The estimation of $x_{\rm br}$ also shows that the
back-reaction should be able to reverse the gas flow. This is also consistent
with the findings of Gárate et al. (2020).

Figure 15: Left: Gas mass flux radial profiles at different times between 10
and 90 kyr. Right: Same, but for the estimation of $x_{\rm br}$, which is
evaluated with the vertically integrated dust-to-gas ratio, that is
$\varepsilon=\Sigma_{\rm d}/\Sigma_{\rm g}$.

To confirm this hypothesis, we finally performed a simulation without back-
reaction as verification. Fig. 16 compares the corresponding radial profiles
of pressure and surface density. Without back-reaction, the pressure maximum
is the same as the pure growth simulation, i.e. at the inner edge of the disc.
In this simulation, no other pressure maximum develops in the disc, even
though the dust experiences growth and fragmentation. As a result, grains
drift towards the inner edge of the disc and are trapped there, as opposed to
the simulation with back-reaction that shows self-trapping at $\sim 50$ au.
With this comparison, we are confident that the mechanism observed is a self-
induced dust trap. Although this mechanism has been observed several times
(Gonzalez et al., 2017a; Garcia, 2018; Pignatale et al., 2019; Vericel &
Gonzalez, 2020), this is the first time it is observed in another, independent
code, which validates its existence. A face-on rendered view of the gas and
dust density evolution during the self-induced dust trap formation is shown in
Fig. 17.

Figure 16: Same as Fig. 14 but with a comparison between simulations with
(F-Hard, solid lines) and without (F-Hard-noBR, dashed lines) back-reaction
and only showing $\Sigma_{\mathrm{d}}$.

Figure 17: Face-on rendered view of the gas (top) and dust (bottom) column
densities at 4 different times during the formation of the self-induced dust
trap in simulation F-Hard.

####  5.3.3 Comparison with the pure growth case

To further point out the major differences that fragmentation produces, we
compare simulations with (F-Hard) and without (noF) fragmentation in this
section. A first example can be seen in Fig. 18, where we show the difference
in dust settling. With fragmentation, the internal parts of the disc are, as
one would expect, less settled than the pure growth case, mainly because the
dust cannot reach large sizes. More generally, the disc can be separated into
3 zones (see black vertical dashed lines):

  1. 1.

An external one ($r\gtrsim r_{\rm frag}$) where fragmentation is impossible.
In this region, fragmentation has no impact and the dust discs reach similar
aspect ratios.

  2. 2.

An intermediate one ($50\lesssim r\lesssim r_{\rm frag}$ au) where the dust
relative velocity is in equilibrium with the fragmentation velocity (see right
panel). In this region, fragmentation limits dust growth, thus the grains are
more coupled to the gas and the disc is less settled than in the pure growth
case.

  3. 3.

An inner one ($r\lesssim 50$ au) where the grains are growing and fragmenting
on short timescales and do not reach an equilibrium with the fragmentation
threshold. In this region, grains have even smaller sizes and the dust disc is
thinner than the gas disc, but much thicker than in the pure growth case.

Figure 18: Left: Disc aspect ratio radial profiles with (F-Hard, blue) and
without (noF, red) fragmentation at 50 kyr. Right: Radial profile of the
relative-to-fragmentation velocity ratio for simulation with fragmentation at
the same time as the left panel. The vertical lines delimit the zone where the
dust relative velocity is in equilibrium with the fragmentation velocity.

A comparison of the trajectories of similar particles is also shown in Fig.
19. Trajectories outside of the fragmentation zone are extremely similar since
the particles experience the same growth path. On the other hand, when
particles enter this zone (salmon colour), their size first evolve along the
$s_{\rm frag}^{-}$ line, that is they tend to fragment slowly (see P2 for
example). Note that P5 or P6 enter the fragmentation zone after typically 50
kyr or more and follow a fragmentation size corresponding to a larger dust-to-
gas ratio. This makes perfect sense, since dust settling and pile-up tend to
concentrate dust in the mid-plane of the disc and increase the dust-to-gas
ratio over time. The formation of the self-induced dust trap can be seen from
the point of views of P2, P3 and P4, as after some time ($\sim 50$ kyr) their
radial drift stops and their size increases near $\sim 50$ au. They are then
able to grow inside the estimated fragmentation zone because the gas structure
is modified locally and the grains relative velocity decreases as they
concentrate.

Figure 19: Trajectories in the ($r,s$) plane of 6 particles initially at 40,
60, 80, 100, 150 and 200 au from the star for simulations with (F-Hard, solid
lines) and without (noF, dashed lines) fragmentation. The fragmentation zone
(salmon colour, computed using the power laws initially describing the disc
structure) is delimited by the fragmentation sizes $s_{\rm frag}^{\pm}$
(dotted lines) and estimated with $\varepsilon=0.01$. $s_{\rm frag}^{-}$ is
also plotted for larger dust-to-gas ratios ($\varepsilon=0.1,0.25$ and 0.5,
from smaller to larger sizes). The optimal drift size $s_{\rm opt}$ is also
added (dash dotted line). The simulation with fragmentation extends up to 100
kyr, that is twice as long as the pure growth case, which reaches a quasi
steady state at this time. To compare them efficiently, the stars represent
the particles state at 50 kyr, while the filled circles mark their state at
100 kyr.

Finally, the shape of the dust size distribution carries these differences, as
one can see in Fig. 20. Following the previous explanations, it is very clear
that the distribution only differs in the fragmentation zone, which we
estimated rather well.

Figure 20: Dust size radial distribution in simulations with (F-Hard, coloured
particles) and without (noF, grey particles) fragmentation at 50 kyr. The
colour represents the ratio $V_{\rm rel}/V_{\rm frag}$. The fragmentation zone
is hatched and estimated using $s_{\rm frag}^{\pm}$ and $\varepsilon=0.2$ (see
also Fig. 13).

####  5.3.4 Effects of the fragmentation model

We implemented two fragmentation models to measure their impact on dust
evolution. Here, we compare the Hard model already presented to the Smooth one
(simulation F-Smooth). The comparison can be seen in Fig. 21. What we gather
from the comparison is that the Smooth model also forms a self-induced dust
trap in the disc. It also tends to form this dust trap slightly closer to the
star. Indeed, this model is less harsh on the dust fragmentation, which shifts
the size gradient towards the star (it allows slightly larger grains at
similar radii). As a consequence, the pressure maximum formed, which is
intimately linked to the position of the dust size gradient, is also located
closer to the star, finally leading to the observed position of the trap. Even
though this tendency is observed, we stress that the differences are extremely
minor, with a radial shift of a few au only. We thus argue that the two
fragmentation models produce largely similar results.

Figure 21: Left: Dust size radial profiles between 0 and 90 kyr for
fragmentation models Hard (solid lines) and Smooth (dashed lines). Right: Same
but for the dust surface density.

The effects of different disc parameters are considered in Appendix B.

####  5.3.5 Synthetic images

We transformed our simulations with fragmentation at 50 kyr into synthetic
observations using the radiative transfer code mcfost (Pinte et al., 2006;
Pinte et al., 2009). At this time, the traps have developed, but the
artificial clumping has not started yet (see Fig. 19 and Section 6.2.) To
achieve this, we recreated a local dust size distribution on each gas particle
from our single set of dust particles. More precisely, we separated a given
simulation into 20 grain size bins and recomputed using the SPH kernel the
density structure at the location of the gas particles for each bin
independently. For each of these locations and from the 20 dust density
fields, mcfost then interpolated the grain size distributions over 100 bins
uniformly distributed in log and sampling grain sizes from $0.03$ to $1000$
$\mu$m. We assume dust grains smaller than 1 $\mu$m follow the gas. The total
grain size distribution integrated over the whole disc is finally normalised,
assuming a power-law distribution $\mathrm{d}n(s)\propto s^{-3.5}\mathrm{d}s$
and the total disc mass. Grains optical properties were calculated using the
Mie theory, assuming astrosilicates composition (Weingartner & Draine, 2001).

Figure 22: Synthetic images at $\lambda=350,850,1300$ and 3000 $\mu$m (left to
right) of simulations with Hard (top) and Smooth (bottom) fragmentation models
at 50 kyr. Images are convolved using a $0.1\times 0.1$ arcsec2 Gaussian beam
(displayed on the bottom left of each panel). The system is assumed to be at a
distance of $140$ pc and seen pole-on.

The central star was represented by a sphere of radius $R_{\star}=2$
$R_{\odot}$ and effective temperature $T_{\rm eff}=4000$ K radiating as a
black body. The disc temperature structure was first computed using $\sim
10^{6}$ photon packets. The images were then computed at
$\lambda=350,850,1300$ and 3000 $\mu$m using $10^{7}$ photon packets at each
wavelength to sample the specific intensity, followed by a ray-tracing
integration to generate the synthetic maps. The system was assumed to be at a
distance of $140$ pc and seen pole-on. The synthetic images are shown in Fig.
22. The dust concentration at $\sim 50$ au translates into a bright and thin
ring at millimetre wavelengths. The observed ring also shows a thinner
structure as the wavelength increases, which is a manifestation of the dust
size distribution being more peaked for higher sizes (the thermal emission of
grains of size $s$ peaks at $\lambda\sim 2\pi s$). Both fragmentation models
give similar images. The major difference is seen at lower wavelengths, where
the Hard model shows a thicker ring than its counterpart.

Considering these synthetic images, the self-induced dust trap mechanism
observed in our simulations should be detectable by instruments such as ALMA,
which is sensitive to the dust thermal emission in the mid-plane of the disc.
This is particularly interesting given the commonness of axisymmetric rings in
recent disc observations (Andrews et al., 2018; Huang et al., 2018). However,
we stress that it might be hard to disentangle this mechanism from others that
produce bright axisymetric rings (e.g. at the edge of a massive planet’s gap
or because of dust sintering, Rice et al., 2006; Okuzumi et al., 2016; Toci et
al., 2020).

##  6 Discussion

###  6.1 The monodisperse approximation

The dust growth module presented in this paper does not solve the Smoluchowski
coagulation equation (Smoluchowski, 1916) but rather considers locally
monodisperse size distributions. This facilitates the implementation as it
becomes a single differential equation to solve numerically for every particle
of the simulation. The monodisperse approximation considers the local
distribution to be highly peaked around a single value that is carried by the
dust SPH particle. By construction, this makes it impossible to track all
sizes at the same time and therefore we only keep the largest ones in memory,
which creates a top-heavy dust size distribution. This is particularly the
case in regions of efficient growth, as illustrated in Fig. 23 between 50 and
100 au in simulation F-Hard. In that regard, small sizes are less represented,
which physically is an acceptable approximation since most of the dust mass is
carried by the larger populations. However, this becomes more concerning when
simulating synthetic images, as most of the stellar light is absorbed and
scattered by tightly coupled grains (i.e. small grains $\lesssim$ 1 $\mu$m).
This directly impacts the scattered light images, but also the disk thermal
structure (and sub-mm images) as the disk midplane is heated by photons
scattered from the disk upper layers. We alleviate part of this issue at the
stage where mcfost reconstructs the size distribution (see Section 5.3.5) by
adding grains smaller than 1 $\mu$m, considering that they follow the gas,
which limits the impact of this caveat. One might also argue that the SPH
method is already poorly adapted to images in scattered light, since it comes
from the upper layers of the disc, which, by effect of settling and overall
vertical density distributions, have a lower resolution. On the contrary, the
mid-plane of the disc is better resolved and contains mm to cm grain sizes,
which makes the synthetic images at millimeter wavelengths more adapted to
this kind of dust growth model. Adding the sub-micron sized grains ensures
that our disk temperature structure is correct and ALMA synthetic images of
self-induced dust traps, such as shown in Fig. 22, should be mostly unaffected
by our dust treatment, with the possible exception of its shortest
wavelengths. In the mid- to far-infrared, under-representation of the 1–100
$\mu$m grain population would mostly affect synthetic images of regions of
efficient growth, and in particular the contrast between bright and dark
rings. We do not compute images in this wavelength range.

The resolution of the Smoluchoswki equation in a code that directly integrates
the equations of motion self-consistently for both gas and dust is extremely
challenging, although methods suitable for use in phantom are under
development (Lombart & Laibe, 2021). The results of such an implementation
providing a continuous size distribution, will precise the extent of the
validity of the monodisperse approximation.

Figure 23: Grain size distribution in simulation F-Hard at 50 kyr in the whole
disc (black) and in the 50–100 au region only (blue).

Lastly, the growth model presented in this paper and based on Stepinski &
Valageas (1997) requires an explicit formulation of the Stokes number —
presently expressed with the Keplerian frequency. This can be a source of
limitation, as it assumes the disc is Keplerian by nature. Significant non-
Keplerian rotations are a possibility, e.g. when the disc is perturbed by
massive companions, which would limit the use of such growth model. Although
alternatives to this issue are beyond the scope of this paper, the Stokes
number would probably be best computed in that scenario with a dynamical
timescale estimated within the simulation itself.

###  6.2 Numerical limitations

As mentioned in section 3.1, the dust growth model is also limited by the
numerical formalism that we employ. Here, we will discuss the one we used in
our simulations: the two-fluid algorithm. More precisely, the handling of two
sets of resolution can become delicate to manage when the two resolutions stop
being comparable. In practice, the gas structure is relatively steady because
of the pressure support. On the other hand, dust particles are pressureless,
which in case of local dust concentrations can become a problem and lead to
numerical clumps. Dust concentrations arise naturally at pressure bumps since
the radial drift vanishes, which increases locally the dust-to-gas ratio. At a
certain point, the dust smoothing length becomes so small that the gas is
effectively invisible to the dust particles, which underestimates the extent
of the dust structures (see for example Fig. 10 and 13). The code used by L08,
LyonSPH, seems to be less prone to this limitation. We suspect this difference
to come from two main reasons:

  1. 1.

The way the smoothing length is handled. In Phantom, the smoothing length is
carefully computed and adjusted using iterations on $h-\rho$, while in LyonSPH
$h$ is independently estimated and thus changes less over time, which might
smooth out the density in the disc. While the method used by Phantom produces
better resolved density structures, it could also be more affected by
numerical clumping.

  2. 2.

The viscosity treatment. While we used the formalism proposed by Lodato &
Price (2010), LyonSPH uses the older Monaghan (1992) formulation. We have
noticed that in LyonSPH, particles tend to be more accreted by the star, which
increases the effective viscosity of the disc over time. As a result, an
effectively more dissipating disc could smooth out the disc structures and
thus limit local dust over-resolutions.

Numerical clumping is hard to detach from in the two-fluid formalism,
especially when considering dust growth that further decouples dust from gas.
However, to limit its effects as much as possible, the best solution is to
adjust the ratio between dust and gas particles, as we already did in our
setup. For example, we want a high enough number of gas particles to resolve
the disc scale height ($\gtrsim 1$M) while having a large gas-to-dust
particles ratio such that the gas smoothing length is as small as possible
compared to the dust’s ($\sim 5-10$). By simple analytical arguments, if we
consider the condition $h_{\rm g}\lesssim h_{\rm d}$ to be fulfilled to avoid
any numerical clumping, this gives a maximum safely attainable dust-to-gas
ratio of

$\varepsilon_{\rm max}\sim\varepsilon_{0}\dfrac{n_{\rm g}^{\rm SPH}}{n_{\rm
d}^{\rm SPH}},$ (59)

where $\varepsilon_{0}$ is the initial dust-to-gas ratio. With great safely
attainable dust-to-gas ratio comes great computational time.

###  6.3 Perspectives

Beyond the presentation of the model implementation in Phantom and the
associated tests, this paper mainly revisited and found similar results as
previous studies that considered dust growth in protoplanetary discs. The main
objective of this paper is to present the algorithm and stress that it is
public and usable by the community. With Phantom’s modularity, this means that
dust growth can be studied in 3D in any given system that was previously out
of reach, all with better overall performances and a lower memory footprint
than LyonSPH.

For example, multiple systems are of much interest, since they can constitute
between a few and about 60% of systems (Rubenstein & Bailyn, 1997; Bellazzini
et al., 2002; Guerrero et al., 2014). Grasping planet formation in general
surely means that dust growth in these poorly studied systems is needed. This
point is especially of interest considering the recent substructures than have
been theoretically and observationally discovered in circumbinary discs
(neglecting dust growth, e.g. Ragusa et al., 2017; Calcino et al., 2019;
Poblete et al., 2019).

Dust properties are also another important element that can be improved in
this model. While we consider perfectly spherical grains composed of a given
species, we are aware that reality is vastly more complicated. For example,
snow lines can change the dust chemical composition, affect their surface
properties and thus their growth and dynamic (e.g. Wada et al., 2009;
Drążkowska & Alibert, 2017; Vericel & Gonzalez, 2020). Another important
element could be the dust porosity, which can be significant considering for
instance the very small filling factors of about $10^{-3}$ that have been
measured on comet 67P/Churyumov–Gerasimenko (Fulle et al., 2015). Dust
porosity can accelerate dust growth and help overcome planet formation
barriers (Garcia, 2018; Garcia & Gonzalez, 2020), which is promising for
planetesimal formation. The implementation of a porosity model is underway in
Phantom and will be the subject of further investigations in the future.

##  7 Conclusion

We presented a dust growth model that considers local monodisperse size
distributions. After specifying the numerical implementation in Phantom, we
presented two tests that ensure the algorithm works as intended. We then
proceeded to circumstellar disc simulations, while first focusing on pure
growth and adding fragmentation effects later on. We finally performed
radiative transfer calculations to produce synthetic images. Our main findings
can be summarised as follows:

  1. 1.

The tests show that the parameters of interest are evaluated with a maximum
relative error of about $10^{-4}$, a satisfying degree of precision.

  2. 2.

Simulations without fragmentation develop a dust concentration at the inner
edge of the disc, which corresponds to the usual gas pressure maximum. The
dust concentrated at this location drifts from the outer parts of the disc and
ends its course well decoupled from the gas. Results are consistent with those
found by L08.

  3. 3.

Simulations with growth and fragmentation show the formation of a self-induced
dust trap as first proposed by Gonzalez et al. (2017a). This trap forms at
several tens of au from the star and is the result of the dust piling-up in
the disc and modifying the gas structure through back-reaction.

  4. 4.

Simulations in which a self-induced dust trap forms produce synthetic images
showing an axisymmetric bright ring at millimetre wavelengths. These should be
detectable by instruments such as ALMA.

## Acknowledgements

The authors would like to thank the referee, Matías Gárate, for constructive
comments and suggestions that helped improved this work. AV would also like to
thank Hossam Aly and Kieran Hirsh for their creative suggestions. This
research was supported by the École Doctorale PHAST (ED 52) of the Université
de Lyon. The authors acknowledge funding from ANR (Agence Nationale de la
Recherche) of France under contract number ANR-16-CE31-0013 (Planet-Forming-
Disks) and thank the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the
Université de Lyon for its financial support within the programme
‘Investissements d’Avenir’ (ANR-11-IDEX-0007) of the French government
operated by the ANR. This project has received funding from the European
Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 823823. DJP and CP acknowledge funding
from the Australian Research Council via grants DP180104235, FT130100034 and
FT170100040. SPH simulations were run at the Common Computing Facility (CCF)
of LABEX LIO and at the PSMN (Pôle Scientifique de Modélisation Numérique) of
the ENS de Lyon, and analysed with Fantomanalysis, which has been developed by
AV. We acknowledge use of the Ozstar supercomputer, funded by Swinburne
University and the Australian Government. Figures were made with Matplotlib
(Hunter, 2007), except for Fig. 17 which was made with Splash (Price, 2007,
2011).

## Data Availability

The Phantom SPH code is available from
https://github.com/danieljprice/phantom. Mcfost is available for use on a
collaborative basis from
https://ipag.osug.fr/~pintec/mcfost/docs/html/overview.html. Fantomanalysis is
available from https://github.com/arnaudvericel/fantomanalysis. The parameter
files for generating our SPH simulations and radiative transfer models are
available upon request.

## References

  * Abod et al. (2019) Abod C. P., Simon J. B., Li R., Armitage P. J., Youdin A. N., Kretke K. A., 2019, ApJ, 883, 192
  * Adachi et al. (1976) Adachi I., Hayashi C., Nakazawa K., 1976, Progress of Theoretical Physics, 56, 1756
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Auffinger & Laibe (2018) Auffinger J., Laibe G., 2018, MNRAS, 473, 796
  * Ballabio et al. (2018) Ballabio G., Dipierro G., Veronesi B., Lodato G., Hutchison M., Laibe G., Price D. J., 2018, MNRAS, 477, 2766
  * Barrière-Fouchet et al. (2005) Barrière-Fouchet L., Gonzalez J.-F., Murray J. R., Humble R. J., Maddison S. T., 2005, A&A, 443, 185
  * Bellazzini et al. (2002) Bellazzini M., Fusi Pecci F., Messineo M., Monaco L., Rood R. T., 2002, AJ, 123, 1509
  * Birnstiel et al. (2009) Birnstiel T., Dullemond C. P., Brauer F., 2009, A&A, 503, L5
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Blum & Wurm (2008) Blum J., Wurm G., 2008, ARA&A, 46, 21
  * Brauer et al. (2008) Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859
  * Calcino et al. (2019) Calcino J., Price D. J., Pinte C., van der Marel N., Ragusa E., Dipierro G., Cuello N., Christiaens V., 2019, MNRAS, 490, 2579
  * Courant et al. (1928) Courant R., Friedrichs K., Lewy H., 1928, Mathematische Annalen, 100, 32
  * Cuello et al. (2019) Cuello N., et al., 2019, MNRAS, 483, 4114
  * Cuello et al. (2020) Cuello N., et al., 2020, MNRAS, 491, 504
  * Dipierro & Laibe (2017) Dipierro G., Laibe G., 2017, MNRAS, 469, 1932
  * Dipierro et al. (2015) Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73
  * Dipierro et al. (2018) Dipierro G., Laibe G., Alexander R., Hutchison M., 2018, MNRAS, 479, 4187
  * Dominik & Tielens (1997) Dominik C., Tielens A. G. G. M., 1997, ApJ, 480, 647
  * Dong et al. (2016) Dong R., Fung J., Chiang E., 2016, ApJ, 826, 75
  * Dong et al. (2018) Dong R., et al., 2018, ApJ, 860, 124
  * Drążkowska & Alibert (2017) Drążkowska J., Alibert Y., 2017, A&A, 608, A92
  * Drążkowska et al. (2014) Drążkowska J., Windmark F., Dullemond C. P., 2014, A&A, 567, A38
  * Drążkowska et al. (2016) Drążkowska J., Alibert Y., Moore B., 2016, A&A, 594, A105
  * Drążkowska et al. (2019) Drążkowska J., Li S., Birnstiel T., Stammler S. M., Li H., 2019, ApJ, 885, 91
  * Dubrulle et al. (1995) Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237
  * Epstein (1924) Epstein P. S., 1924, Phys. Rev., 23, 710
  * Fromang & Papaloizou (2006) Fromang S., Papaloizou J., 2006, A&A, 452, 751
  * Fulle et al. (2015) Fulle M., et al., 2015, ApJ, 802, L12
  * Gárate et al. (2020) Gárate M., Birnstiel T., Drążkowska J., Stammler S. M., 2020, A&A, 635, A149
  * Garaud et al. (2004) Garaud P., Barrière-Fouchet L., Lin D. N. C., 2004, ApJ, 603, 292
  * Garcia (2018) Garcia A., 2018, Phd thesis, Université de Lyon, https://tel.archives-ouvertes.fr/tel-01977317
  * Garcia & Gonzalez (2020) Garcia A. J. L., Gonzalez J.-F., 2020, MNRAS, 493, 1788
  * Gonzalez et al. (2015) Gonzalez J.-F., Laibe G., Maddison S. T., Pinte C., Ménard F., 2015, Planet. Space Sci., 116, 48
  * Gonzalez et al. (2017a) Gonzalez J.-F., Laibe G., Maddison S. T., 2017a, MNRAS, 467, 1984
  * Gonzalez et al. (2017b) Gonzalez J. F., Laibe G., Maddison S. T., 2017b, MNRAS, 472, 1162
  * Gonzalez et al. (2020) Gonzalez J.-F., et al., 2020, MNRAS, submitted
  * Guerrero et al. (2014) Guerrero C. A., Orlov V. G., Monroy-Rodríguez M. A., Voitsekhovich V. V., 2014, AJ, 147, 28
  * Haghighipour (2005) Haghighipour N., 2005, MNRAS, 362, 1015
  * Huang et al. (2018) Huang J., et al., 2018, ApJ, 869, L42
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Johansen & Youdin (2007) Johansen A., Youdin A., 2007, ApJ, 662, 627
  * Kanagawa et al. (2017) Kanagawa K. D., Ueda T., Muto T., Okuzumi S., 2017, ApJ, 844, 142
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Krijt & Ciesla (2016) Krijt S., Ciesla F. J., 2016, ApJ, 822, 111
  * Krijt et al. (2016) Krijt S., Ormel C. W., Dominik C., Tielens A. G. G. M., 2016, A&A, 586, A20
  * Kwok (1975) Kwok S., 1975, ApJ, 198, 583
  * Laibe (2014) Laibe G., 2014, MNRAS, 437, 3037
  * Laibe & Price (2012a) Laibe G., Price D. J., 2012a, MNRAS, 420, 2345
  * Laibe & Price (2012b) Laibe G., Price D. J., 2012b, MNRAS, 420, 2365
  * 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
  * Laibe et al. (2008) Laibe G., Gonzalez J.-F., Fouchet L., Maddison S. T., 2008, A&A, 487, 265
  * Laune et al. (2020) Laune J., Li H., Li S., Li Y.-P., Walls L. G., Birnstiel T., Drążkowska J., Stammler S., 2020, ApJ, 889, L8
  * Li et al. (2005) Li H., Li S., Koller J., Wendroff B. B., Liska R., Orban C. M., Liang E. P. T., Lin D. N. C., 2005, ApJ, 624, 1003
  * Li et al. (2019a) Li Y.-P., et al., 2019a, ApJ, 878, 39
  * Li et al. (2019b) Li R., Youdin A. N., Simon J. B., 2019b, ApJ, 885, 69
  * Lissauer & Stewart (1993) Lissauer J. J., Stewart G. R., 1993, in Levy E. H., Lunine J. I., eds, Protostars and Planets III. pp 1061–1088
  * Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
  * Lombart & Laibe (2021) Lombart M., Laibe G., 2021, MNRAS, 501, 4298
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Mentiplay et al. (2019) Mentiplay D., Price D. J., Pinte C., 2019, MNRAS, 484, L130
  * Misener et al. (2019) Misener W., Krijt S., Ciesla F. J., 2019, ApJ, 885, 118
  * Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
  * Monaghan (1997) Monaghan J. J., 1997, Journal of Computational Physics, 138, 801
  * Muñoz & Lai (2016) Muñoz D. J., Lai D., 2016, ApJ, 827, 43
  * Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
  * Nealon et al. (2020) Nealon R., Cuello N., Alexander R., 2020, MNRAS, 491, 4108
  * Okuzumi et al. (2016) Okuzumi S., Momose M., Sirono S.-i., Kobayashi H., Tanaka H., 2016, ApJ, 821, 82
  * Ormel & Spaans (2008) Ormel C. W., Spaans M., 2008, ApJ, 684, 1291
  * Paardekooper & Mellema (2004) Paardekooper S.-J., Mellema G., 2004, A&A, 425, L9
  * Pérez et al. (2019) Pérez S., Casassus S., Baruteau C., Dong R., Hales A., Cieza L., 2019, AJ, 158, 15
  * Pignatale et al. (2019) Pignatale F. C., Gonzalez J. F., Bourdon B., Fitoussi C., 2019, MNRAS, 490, 4428
  * Pinilla et al. (2012) Pinilla P., Benisty M., Birnstiel T., 2012, A&A, 545, A81
  * Pinte et al. (2006) Pinte C., Ménard F., Duchêne G., Bastien P., 2006, A&A, 459, 797
  * Pinte et al. (2009) Pinte C., Harries T. J., Min M., Watson A. M., Dullemond C. P., Woitke P., Ménard F., Durán-Rojas M. C., 2009, A&A, 498, 967
  * Pinte et al. (2018) Pinte C., et al., 2018, ApJ, 860, L13
  * Pinte et al. (2019) Pinte C., et al., 2019, Nature Astronomy, 3, 1109
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Poblete et al. (2019) Poblete P. P., Cuello N., Cuadra J., 2019, MNRAS, 489, 2204
  * Price (2007) Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159
  * Price (2011) Price D. J., 2011, SPLASH: An Interactive Visualization Tool for Smoothed Particle Hydrodynamics Simulations, Astrophysics Source Code Library (ascl:1103.004)
  * Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
  * Price & Laibe (2020) Price D. J., Laibe G., 2020, MNRAS, 495, 3929
  * Price et al. (2017) Price D. J., et al., 2017, PHANTOM: Smoothed particle hydrodynamics and magnetohydrodynamics code (ascl:1709.002)
  * Price et al. (2018a) Price D. J., et al., 2018a, Publ. Astron. Soc. Australia, 35, e031
  * Price et al. (2018b) Price D. J., et al., 2018b, MNRAS, 477, 1270
  * Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
  * Rice et al. (2006) Rice W. K. M., Armitage P. J., Wood K., Lodato G., 2006, MNRAS, 373, 1619
  * Rubenstein & Bailyn (1997) Rubenstein E. P., Bailyn C. D., 1997, ApJ, 474, 701
  * Safronov (1969) Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka.
  * Schäfer et al. (2017) Schäfer U., Yang C.-C., Johansen A., 2017, A&A, 597, A69
  * Schoonenberg et al. (2018) Schoonenberg D., Ormel C. W., Krijt S., 2018, A&A, 620, A134
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Smoluchowski (1916) Smoluchowski M. V., 1916, Zeitschrift fur Physik, 17, 557
  * Stepinski & Valageas (1996) Stepinski T. F., Valageas P., 1996, A&A, 309, 301
  * Stepinski & Valageas (1997) Stepinski T. F., Valageas P., 1997, A&A, 319, 1007
  * Toci et al. (2020) Toci C., Lodato G., Fedele D., Testi L., Pinte C., 2020, ApJ, 888, L4
  * Vericel & Gonzalez (2020) Vericel A., Gonzalez J.-F., 2020, MNRAS, 492, 210
  * Veronesi et al. (2020) Veronesi B., et al., 2020, MNRAS, 495, 1913
  * Wada et al. (2009) Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2009, in Henning T., Grün E., Steinacker J., eds, Astronomical Society of the Pacific Conference Series Vol. 414, Cosmic Dust - Near and Far. p. 347
  * Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
  * Weidenschilling (1997) Weidenschilling S. J., 1997, Icarus, 127, 290
  * Weingartner & Draine (2001) Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
  * Whipple (1973) Whipple F. L., 1973, NASA Special Publication, 319, 355
  * Yang et al. (2017) Yang C. C., Johansen A., Carrera D., 2017, A&A, 606, A80
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Zhu et al. (2015) Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88
  * Zsom & Dullemond (2008) Zsom A., Dullemond C. P., 2008, A&A, 489, 931

##  Appendix A Pure growth: influence of the initial grain size

The effects of the initial size has been discussed by L08, where they found
that dust growth was fast between 1 and $10$ $\mu$m and therefore the choice
of the initial size between those values was inconsequential on dust
evolution. The initial size is numerically very important to control, since
small grains dictate the timestepping in the two-fluid dust method. As a
matter of fact, slightly increasing the initial size can lead to a simulation
being a few times faster.

Here, we tested 3 initial sizes with 10, 30 (reference value used throughout
the paper) and 50 $\mu$m. We can see their effects on grain growth and
settling in Figs. 24 and 25.

Figure 24: Dust size radial profiles in pure growth simulations for $s_{0}=$
10 (blue), 30 (green) and 50 (red) $\mu$m at 0 (top left), 10 (top right), 25
(bottom left) and 45 (bottom right) kyr.

Figure 25: Same as Fig. 24 but for the gas (black dashed line) and dust (solid
coloured lines) disc scale heights.

Grains in the innermost parts of the disc reach similar sizes very rapidly (a
few thousand years), since the growth timescale there is the shortest. This
allows the grains to quickly forget their initial size. In terms of vertical
settling, the inner parts of the discs also match exactly as quickly. At later
times, the distance to the star for which simulations have converged to the
same state increases. At 45 kyr for instance, we can notice that the
simulations are similar up to $\sim 100$ au.

It seems clear that two simulations with different initial sizes can reach a
similar state if enough time has passed so that most of the disc can forget
its initial size, that is respecting the condition $\Delta s\gg s_{0}$, where
$\Delta s$ is the size gained through collisions. In practice, the initial
size is therefore constrained upward by this limitation, and downward by the
timestepping constraint we already mentioned. Our choice of $30$ $\mu$m seems
adequate for most cases but can be adapted to the user’s needs.

##  Appendix B Influence of the disc model on self-induced dust trap formation

###  B.1 Inner boundary

Figure 26: Simulation for the disc with $r_{\mathrm{in}}=10$ au. Dust size
radial distribution at the end of the simulation where particles are coloured
with their ratio $V_{\mathrm{rel}}/V_{\mathrm{frag}}$ (left) and radial
profiles of the dust surface density (center) and of $x_{\mathrm{br}}$ (right)
at different times.

One may wonder whether the location of the disc inner boundary can affect the
formation of self-induced dust traps. We ran a new simulation of our standard
disc with a smaller inner radius of $r_{\mathrm{in}}=10$ au, up to the trap
formation and before the onset of runaway growth. Figure 26 shows that the
self-induced dust trap forms at the same location of $\sim 50$ au, where both
the grain size (left) and dust surface density (centre) peak, and where
$x_{\mathrm{br}}$ becomes larger than 1. Due to the smaller inner radius, the
disc regions interior to the trap can supply more material outwards, helping
the trap to from sooner than in the simulation with $r_{\mathrm{in}}=20$ au
presented in Section 5.3.

###  B.2 Surface density and temperature profiles

Figure 27: Same as Fig. 26 for the disc with $p=1$ and $q=1/2$.

We present the case of a disc with $p=1$ and $q=1/2$, all other parameters
being equal to those of our standard disc. This configuration also results in
self-induced dust trap formation, as can be seen in Fig. 27. The different
profiles, and in particular the smaller value of $q$, lead to the trap forming
at larger radii, as reported in Gonzalez et al. (2017a), here at $\sim 120$
au. Since evolution timescales increase with radius, the trap also forms later
than in our standard disc.