# The role of the drag force in the gravitational stability of dusty planet-
forming disc – II. Numerical simulations

Cristiano Longarini,1,2,3,4 Philip J. Armitage,2,3 Giuseppe Lodato,1 Daniel J.
Price4, Simone Ceppi1,4
1Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16,
20133 Milano, Italy
2Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY
11794, USA
3Center for Computational Astrophysics, Flatiron Institute, New York, NY
10010, USA
4School of Physics and Astronomy, Monash University, Clayton, VIC 3800,
Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Young protostellar discs are likely to be both self-gravitating, and to
support grain growth to sizes where the particles decoupled from the gas. This
combination could lead to short-wavelength fragmentation of the solid
component in otherwise non-fragmenting gas discs, forming Earth-mass solid
cores during the Class 0 / I stages of Young Stellar Object evolution. We use
three-dimensional smoothed particle hydrodynamics simulations of two-fluid
discs, in the regime where the Stokes number of the particles ${\rm St}>1$, to
study how the formation of solid clumps depends on the disc-to-star mass
ratio, the strength of gravitational instability, and the Stokes number.
Gravitational instability of the simulated discs is sustained by local
cooling. We find that the ability of the spiral structures to concentrate
solids increases with the cooling time, and decreases with the Stokes number,
while the relative dynamical temperature between gas and dust of the particles
decreases with the cooling time and the disc-to-star mass ratio, and increases
with the Stokes number. Dust collapse occurs in a subset of high disc mass
simulations, yielding clumps whose mass is close to linear theory estimates,
namely $1$–$10~{}\text{M}_{\oplus}$. Our results suggest that if planet
formation occurs via this mechanism, the best conditions correspond to near
the end of the self-gravitating phase, when the cooling time is long and the
Stokes number close to unity.

###### keywords:

accretion, accretion discs — turbulence — planets and satellites: formation

††pubyear: 2023††pagerange: The role of the drag force in the gravitational
stability of dusty planet-forming disc – II. Numerical simulations–18

##  1 Introduction

Direct imaging of protostellar discs, in scattered light and especially at
sub-mm wavelengths, has opened debate as to how the timescales for planet
formation align with the established evolutionary sequence for Young Stellar
Objects (Adams et al., 1987). Many discs show substructure (Andrews, 2020)
that is consistent with — and often interpreted as — the theoretically
expected signature of planet-disc interaction (Dong et al., 2015; Zhang et
al., 2018; Lodato et al., 2019). A subset of these discs, including HL Tau
(ALMA Partnership et al., 2015), IRS 63 (Segura-Cox et al., 2020) and GY91
(Sheehan & Eisner, 2018) are young, with ages estimated to be less than $1\
\text{Myr}$. Assuming a planet origin for the substructure, a substantial part
of the planet formation process must overlap with the time when protostellar
discs at 10–100 au scales are likely self-gravitating. How this occurs is
unclear. Classical planetesimal-only variants of Core Accretion model
(Safronov, 1969; Pollack et al., 1996) are not viable; meeting the timescale
constraint at 10 au is already non-trivial, and at substantially larger radii
sufficiently massive planet formation is not possible at all. Models with a
dominant pebble component (e.g. Rosenthal & Murray-Clay, 2018; Morbidelli,
2020; Guilera et al., 2020; Chambers, 2021; Cummins et al., 2022; Jiang &
Ormel, 2022), or those with large-scale core migration (Levison et al., 2010),
are more promising, but the ability of these processes to form the inferred
large-orbital-radius population of planets has not been fully established.

Alternatively, planets may form via gravitational instability (GI) (Boss,
1997) in a self-gravitating protostellar disk. Self-gravity is expected to be
an important process during the Class 0/I phase of protostellar evolution, and
is most vigorous at large orbital radii where the cooling time is short
(Clarke, 2009; Rafikov, 2009) and infall continues. GI in a gaseous disc,
however, is a poor candidate for forming the ALMA-inferred planet population,
because it substantially overshoots the inferred masses. The mass of a gaseous
fragment created by GI is of the order of the Jeans mass of the spiral
perturbation, and for typical protoplanetary discs its value is, assuming
Toomre $Q=1$ (Toomre, 1964; Kratter & Lodato, 2016)

$M_{\rm J}=\frac{4c_{\rm s}^{4}}{G^{2}\Sigma}\approx 12M_{\rm
Jup}\left(\frac{M_{*}}{1M_{\odot}}\right)\left(\frac{H/R}{0.1}\right)^{3},$
(1)

typically between $1$–$10~{}\text{M}_{\text{J}}$ (where $M_{\rm Jup}$ is the
mass of Jupiter) depending on the disc aspect ratio. However, this is the
initial mass of the object: subsequently the fragment would start accreting
material belonging to the accretion disc, increasing its mass and likely
becoming a brown dwarf (Kratter et al., 2010; Kratter & Lodato, 2016).

A mechanism that may allow lower mass planets to form as a consequence of GI
was proposed by Rice et al. (2005). Concentration of solid particles in spiral
arms, together with vertical settling, can lead to gravitational collapse in
the solid component even in a non-fragmenting gravitationally unstable gaseous
disc (Gibbons et al., 2012, 2014, 2015). It is known that gas spiral arms act
as dust traps, since they are pressure maxima (Rice et al., 2005; Shi et al.,
2016). Solid particles collect inside them (Dipierro et al., 2015a), reaching
dust to gas ratio that can be of the order of unity. Conversely, the
interaction between gas spiral arms and dust particles can excite them,
imparting random motions that reduce the peak density and potential for
collapse (Riols et al., 2020). Walmswell et al. (2013) found that large dust
particles experience gravitational scattering by the spiral arms, while Booth
& Clarke (2016) related the level of excitation of solid particles to the
aerodynamical coupling and the cooling factor. Marginally coupled solid
particles are less excited by spiral arms, while, in rapidly cooled discs, the
level of dust excitation is higher. Baehr & Zhu (2021) confirmed this trend
through 3D shearing box simulations.

In this work, we perform three-dimensional smoothed particle hydrodynamics
(SPH) simulations using the phantom code (Price et al., 2018) of gas and dust
in gravitationally unstable protostellar discs. The aim is to investigate the
interplay between GI and the dynamics of particles that are aerodynamically
coupled to the gas, and to understand under which conditions it is possible to
form solid cores inside GI spirals. In Section 2, we recall the basics of the
aerodynamical coupling between gas and dust, and the role of dust in GI. In
Section 3 we briefly describe the SPH code phantom (Price et al., 2018) and
the set of simulations we have performed for this work. In Section 4 we
analyse the simulations, show our main results, and compare our findings with
previous works. In Section 5, we discuss our findings in an evolutionary
perspective. Finally, in Section 6 we present our conclusions.

##  2 Drag force and gravitational instability

###  2.1 Radial drift

To zeroth order, gas dynamics is dictated by the hydrodynamical equations.
Neglecting the role of the disc self gravity, in centrifugal equilibrium, the
azimuthal component of the gas velocity is

$v_{\phi,g}^{2}=v_{k}^{2}+\frac{1}{\rho_{g}}\frac{\text{d}P}{\text{d}R},$ (2)

where $v_{k}^{2}=GM_{\star}/R$ is the Keplerian speed and the second term is
the negative contribution of the pressure gradient. It is possible to
explicitly compute the pressure gradient knowing the disc temperature
structure: the gas azimuthal speed is slightly sub-Keplerian, and the
correction term is of the order of $(H/R)^{2}\sim 0.01$, where $H$ is the disc
thickness. In addition, the gas has also a negative radial velocity, given by
viscosity.

Conversely, solid particle dynamics is dictated by both gravitational and drag
force. Since dust is not supported by pressure, the basic velocity of gas and
dust is different, but they interact through a drag force, that modifies the
relative velocity. Writing $\Delta\mathbf{u}=\mathbf{v}_{\rm
d}-\mathbf{v}_{\rm g}$, the drag force $\mathbf{F}_{\rm d}$ depends upon the
relative velocity as

$|\textbf{F}_{\rm d}|=\frac{1}{2}C_{\mathrm{D}}\pi s^{2}\rho\Delta u^{2},$ (3)

where $C_{\rm D}$ is a coefficient that depends on the drag regime, $s$ is the
dust particle size and $\rho$ is the total density. The different
aerodynamical coupling regimes depend upon the size of dust particles relative
to the gas mean free path. It is possible to define a dimensionless parameter,
called the Knudsen number (Paardekooper & Mellema, 2006), that measures this
property

$\text{Kn}=\frac{9\lambda_{g}}{4s},$ (4)

where $s$ is the dust particle size and $\lambda_{\rm g}$ is the gas
particles’ mean free path, given by

$\lambda_{\rm g}=\frac{5\pi}{64\sqrt{2}}\frac{\mu
m_{p}}{\rho_{\text{g}}\sigma_{\text{coll}}},$ (5)

where $\sigma_{\text{coll}}=2.367\times 10^{-15}\text{cm}^{2}$ is the cross-
section of gas particles, $\mu=2.1$ is the mean molecular weight, $\rho_{\rm
g}$ is the gas density and $m_{p}$ is the mass of the proton. Typically, in
protoplanetary discs, $\text{Kn}>1$, meaning that the dust particles’ size is
smaller than the gas mean free path: this is called the Epstein regime.
However, when the disc is very massive, there is a transition between Epstein
and Stokes regime ($\text{Kn}<1)$. The main impact of this transition is that
in the inner denser region of the disc, where particles’ size is comparable
with the gas mean free path, the Stokes number is higher compared to Epstein
regime. The $C_{\text{D}}$ coefficient is given by (Fassio & Probstein, 1970)

$C_{\mathrm{D}}=\left\\{\begin{array}[]{ccc}\frac{8}{3}\frac{c_{g}}{\Delta
u}&\text{Kn}>1&\\\ 24\text{Re}^{-1}&\text{Kn}<1,&\text{Re}<1\\\
24\text{Re}^{-0.6}&\text{Kn}<1,&1<\text{Re}<800\\\
0.44&\text{Kn}<1,&\text{Re}>800\end{array}\right.$ (6)

where $c_{g}$ is the gas sound speed and the Reynolds number based upon the
difference of velocity is given by

$\text{Re}=\frac{2s\Delta u}{\nu}=4\frac{s\Delta u}{\lambda_{g}c_{g}},$ (7)

and the last equivalence is true for collisional viscosity. We also define the
stopping time, i.e. the time needed to modify the relative velocity between
gas and dust. The longer it is, the less the particles are coupled. For
spherical grains, it is given by

$t_{\text{s}}=\frac{m_{d}\Delta
u}{|\mathbf{F}_{D}|}=\frac{4\pi\rho_{0}s^{3}\Delta
u}{3|\mathbf{F}_{D}|}=\frac{8\rho_{0}s}{3C_{D}\rho\Delta u},$ (8)

where $m_{d}$ is the dust grain mass and $\rho_{0}$ is the dust grain’s
intrinsic density. In order to measure the strength of the aerodynamical
coupling, we define a dimensionless Stokes number as the ratio between the
stopping time and the dynamical one: in general, its expression is

$\text{St}=t_{\text{s}}\Omega=\frac{8}{3C_{\mathrm{D}}}\left(\frac{\rho_{\mathrm{0}}}{\rho}\right)\left(\frac{v_{\mathrm{k}}}{\Delta
u}\right)\left(\frac{s}{r}\right),$ (9)

and, in the Epstein regime, its value is

$\text{St}=1\left(\frac{\Sigma}{0.2\text{g}/\text{cm}^{2}}\right)^{-1}\left(\frac{\rho_{0}}{3\text{g}/\text{cm}^{3}}\right)\left(\frac{s}{1\text{mm}}\right).$
(10)

The smaller the Stokes number is, the more tightly the particles are coupled:
so, for $\text{St}\to 0$, dust dynamics follows the gas dynamics, while for
$\text{St}\to\infty$, the two fluids do not influence each other
aerodynamically. In a smooth axisymmetric disc, the drag force has dramatic
physical consequences. Because of the azimuthal difference of speed between
gas and dust, solid particles experience a headwind that slows them down,
giving them a negative radial velocity. This effect is called radial drift,
and its timescale for ${\rm St}\sim 1$ is of the order of $\sim 100$yr,
several orders of magnitude shorter than the disc lifetime. This is the “metre
sized barrier", so called because ${\rm St}\sim 1$ corresponds to roughly this
physical size in the inner disc (Weidenschilling, 1977). Radial drift can be
stopped in discs where the presence of a gas pressure maximum creates a “dust
trap”. Indeed, in a gas pressure maximum, the pressure gradient is zero, and
hence the velocity difference between gas and dust vanishes and no radial
drift occurs. Several mechanisms have been proposed to form dust traps, such
as gaps made by planets (Paardekooper & Mellema, 2004, 2006; Ayliffe et al.,
2012; Pinilla et al., 2012), zonal flows (Johansen et al., 2009; Simon &
Armitage, 2014), or spirals induced by gravitational instabilities (Dipierro
et al., 2015a). Concentration of solid material by aerodynamic effects must be
considered in assessing models for planet formation in early protoplanetary
stages.

###  2.2 Gravitational instability

In an axisymmetric thin disc composed of a single fluid, the dispersion
relation for linear, tightly-wound density waves is

$D(\omega,k,m)=(\omega-m\Omega)^{2}-c^{2}k^{2}+2\pi G\Sigma{|k|}-\kappa^{2},$
(11)

where $\omega$ is the perturbation frequency, $\Omega(R)$ is the angular
velocity, $m/R$ is the azimuthal wavenumber of the perturbation, $c(R)$ the
sound speed of the fluid, $k$ the perturbation radial wavenumber, $\kappa(R)$
the epicyclic frequency and $\Sigma(R)$ the disc surface density. The well-
known instability criterion for axisymmetric disturbances is (Toomre, 1964)

$Q=\frac{c\kappa}{\pi G\Sigma}<1,$ (12)

that identifies the parameter space for which $\omega^{2}<0$. It is possible
to generalize the instability taking into account a second fluid component,
i.e. particles in the protostellar case. In the context of galactic dynamics,
the gravitational role of a second component has been investigated by Jog &
Solomon (1984) and Bertin & Romeo (1988). Longarini et al. (2023) applied this
method to protoplanetary discs, taking into account also the role of the drag
force between gas and dust. When the second component is considered, the
outcome of GI can significantly change. In the protostellar scenario, the
three fundamental parameters are the relative concentration of the dust
component to that of the gas

$\epsilon=\frac{\Sigma_{d}}{\Sigma_{g}},$ (13)

the dust relative temperature

$\xi=\left(\frac{c_{d}}{c_{g}}\right)^{2},$ (14)

where $c_{d}$ is the dust dispersion velocity, and the Stokes number

$\text{St}=t_{s}\Omega.$ (15)

Two regimes of instability can be identified, depending on whether the
instability is triggered by the gas or the dust. When the instability is
driven by the gas component, the most unstable wavelength is close to the gas-
only one, and the role of dust is negligible. Conversely, when the instability
is controlled by the dust, the perturbation wavelength is much smaller. This
different kind of instability happens when the dust is sufficiently cold,
decoupled and abundant: the threshold between these two regimes is given by
(Bertin & Romeo, 1988)

$\epsilon\xi^{-1/2}>1,$ (16)

or, taking into account the drag interaction (Longarini et al., 2023)

$\epsilon\xi^{-1/2}(1+0.72\text{St}^{-1.36})>1.$ (17)

Dust driven instability could have important consequences for planet formation
because when the most unstable wavelength (i.e. the Jeans wavelength) is
smaller, the Jeans mass will be small as well. In the gas-only model, the
Jeans length is

$\lambda^{\text{1f}}_{J}=\frac{2c_{g}^{2}}{G\Sigma_{g}}=\sqrt{\frac{2}{\pi}}\frac{c_{g}\Omega}{G\rho_{g}},$
(18)

where we used $\rho_{g}=\Sigma_{g}/\sqrt{2\pi}H$, and the Jeans mass is

$M^{\text{1f}}_{J}=\frac{4}{3}\pi\rho_{g}\left(\lambda^{\text{1f}}_{J}\right)^{3}=\frac{4\pi}{3}\left(\frac{2}{\pi}\right)^{3/2}\frac{c_{g}^{3}\Omega^{3}}{G^{3}\rho_{g}^{2}}.$
(19)

When the instability becomes dust driven, the value of the most unstable
wavelength changes dramatically. It is possible to compute it through the
dispersion relation of the two fluid component model with drag force
(Longarini et al., 2023), and its value depends on the Stokes number, the
relative temperature and the dust to gas ratio. In particular, the gas-only
limit is obtained for $\text{St}\to 0$, $\epsilon\to 0$ and $\xi\to 1$. For
convenience, we define the ratio between the Jeans wavelength of one fluid
component model and of the two fluid one with drag force as
$\Lambda(\epsilon,\xi,\text{St})$: the value of $\Lambda$ is between 0 and 1,
where $\Lambda=1$ is the one fluid limit. Indeed, by computing the most
unstable wavelength, it is possible to extract $\Lambda$. Hence, the Jeans
mass of the two-component fluid model is

$M^{\text{2f}}_{J}=M^{\text{1f}}_{J}\Lambda^{3}.$ (20)

In Appendix A, we show plots of $\Lambda^{3}$ as a function of
$(\epsilon,\xi,\text{St})$.

##  3 Numerical simulations

In this work, we perform numerical SPH simulations of gas and dust
protostellar discs using the code phantom (Price et al., 2018). This code is
widely used in the astrophysical community to study gas and dust dynamics in
accretion discs (Ragusa et al., 2017; Ceppi et al., 2022), both in a single
fluid mixture (Veronesi et al., 2020) or dust-as-particles approach (Aly et
al., 2021). In this work, we use the dust-as-particles formulation.

###  3.1 Two-fluid gas and dust mixtures

The dust formulation is based on the continuum fluid equations in the form

$\frac{\partial\rho_{\mathrm{g}}}{\partial
t}+\left(\mathbf{v}_{\mathrm{g}}\cdot\nabla\right)\rho_{\mathrm{g}}=-\rho_{\mathrm{g}}\left(\nabla\cdot\mathbf{v}_{\mathrm{g}}\right),$
(21)

$\frac{\partial\rho_{\mathrm{d}}}{\partial
t}+\left(\mathbf{v}_{\mathrm{d}}\cdot\nabla\right)\rho_{\mathrm{d}}=-\rho_{\mathrm{d}}\left(\nabla\cdot\mathbf{v}_{\mathrm{d}}\right),$
(22)

$\frac{\partial\mathbf{v}_{\mathrm{g}}}{\partial
t}+\left(\mathbf{v}_{\mathrm{g}}\cdot\nabla\right)\mathbf{v}_{\mathrm{g}}=-\frac{\nabla
P}{\rho_{\mathrm{g}}}+\frac{K}{\rho_{\mathrm{g}}}\left(\mathbf{v}_{\mathrm{d}}-\mathbf{v}_{\mathrm{g}}\right),$
(23)

$\frac{\partial\mathbf{v}_{\mathrm{d}}}{\partial
t}+\left(\mathbf{v}_{\mathrm{d}}\cdot\nabla\right)\mathbf{v}_{\mathrm{d}}=-\frac{K}{\rho_{\mathrm{d}}}\left(\mathbf{v}_{\mathrm{d}}-\mathbf{v}_{\mathrm{g}}\right),$
(24)

where the subscripts $g$ and $d$ refer to gas and dust properties. We define
the stopping time, that is given by

$t_{s}\equiv\frac{\rho_{\text{g}}\rho_{\text{d}}}{K(\rho_{\text{g}}+\rho_{\text{d}})},$
(25)

where the drag coefficient $K$ depends on the physical drag regime the system
is in. In general, it is given by

$K=\rho_{\mathrm{g}}\rho_{\mathrm{d}}\frac{1}{2}C_{\mathrm{D}}\frac{\pi
s^{2}}{m_{\text{d }}}|\Delta v|,$ (26)

where $C_{D}$ is defined as before.

In the phantom implementation, the two phases are modelled as two distinct
sets of particles: hereafter, we adopt the convention from Monaghan &
Kocharyan (1995), and refer to gas particles with the subscripts $a,b,c$ and
to dust particles with $i,j,k$. Gas and dust densities are computed by
weighted summation over the particles of the same type according to

$\rho_{a}=\sum_{b}m_{b}W_{ab}\left(h_{a}\right);\quad
h_{a}=h_{\mathrm{fact}}\left(\frac{m_{a}}{\rho_{a}}\right)^{1/3},$ (27)

$\rho_{i}=\sum_{j}m_{j}W_{ij}\left(h_{i}\right);\quad
h_{i}=h_{\mathrm{fact}}\left(\frac{m_{i}}{\rho_{i}}\right)^{1/3},$ (28)

where the kernel $W$ is the same for gas and dust, $h$ is the smoothing length
and $h_{\text{fact}}=1/2$. The drag terms of the equations of motion are
discretized by using a “double hump" kernel (Fulk & Quinn, 1996) that, instead
of the bell-shaped kernel, goes to zero at $r=0$ and has a peak at
$r/h\lesssim 1$. 111Laibe & Price (2012a) showed that using a double-hump
kernel gives a factor of 10 better accuracy at no additional computational
cost. For further information, see Section $\mathbf{2.13.4}$ of Price et al.
(2018). In these simulations, we are using the velocity reconstruction
procedure presented in Price & Laibe (2020).

###  3.2 Heating and cooling

In order to take into account the effect of cooling or heating phenomena, we
write the complete equation for the evolution of gas internal energy $e$

$\frac{\partial e}{\partial
t}+\left(\mathbf{v}_{\mathrm{g}}\cdot\nabla\right)e=-\frac{P}{\rho_{\mathrm{g}}}\left(\nabla\cdot\mathbf{v}_{\mathrm{g}}\right)+\Lambda_{\text{shock}}-\frac{\Lambda_{\text{cool}}}{\rho_{g}}+\frac{\Lambda_{\mathrm{drag}}}{\rho_{g}},$
(29)

where the first term on the RHS is the $P\text{d}V$ work, the second is a
heating term due to the shock viscosity, the third is the cooling of the disc
and the last term is the drag heating term. In this work, we assume an
adiabatic equation of state. For an ideal gas, it is possible to link pressure
and density as follows

$P=(\gamma-1)\rho_{g}e=\frac{c_{g}^{2}\rho_{g}}{\gamma},$ (30)

where $\gamma=5/3$ and $c_{g}$ is the gas sound speed, that is initialized as
a power law $c_{g}\propto R^{-0.25}$.

The shock viscosity term can be written as

$\Lambda_{\text{shock}}=c_{1}\alpha^{\text{AV}}\frac{h_{g}}{H}+c_{2}\beta^{\text{AV}}\left(\frac{h_{g}}{H}\right)^{2},$
(31)

where $c_{1},c_{2}$ are two numerical factors, whose dimension is a specific
energy per unit time, $\alpha^{\text{AV}}$ and $\beta^{\text{AV}}$ are
respectively the linear and the quadratic viscosity coefficients, and
$h_{g}/H_{g}$ is related to the numerical resolution222This quantity tells how
many smoothing lengths are included in the disc thickness.. The viscosity term
is dissipative, so it heats the disc. In the phantom simulations we are not
using the disc-viscosity flag, meaning that the shock capturing viscosity is
not described by an $\alpha_{\text{SS}}$ (Shakura & Sunyaev, 1973)
prescription. We did so since in these systems the main driver of angular
momentum transport is GI.

For the cooling we use the prescription from Gammie (2001) and Rice et al.
(2004), in which the cooling time $t_{\text{cool}}$ is proportional to the
dynamical time, with a factor of proportionality $\beta_{\text{cool}}$

$t_{\text{cool}}=\beta_{\text{cool}}\Omega^{-1}$ (32)

Under the assumption that the transfer of angular momentum driven by gravito-
turbulence occurs locally (Lodato & Rice, 2004; Béthune et al., 2021), we can
relate the cooling parameter to an effective $\alpha-$viscosity parameter

$\alpha_{\text{GI}}=\frac{4}{9}\frac{1}{\gamma(\gamma-1)\beta_{\text{cool}}}.$
(33)

We need to introduce a cooling prescription for the system in order to trigger
GI. Finally, the drag heating term is

$\Lambda_{\text{drag}}=K|\mathbf{v}_{d}-\mathbf{v}_{g}|^{2}.$ (34)

Currently, phantom neglects any thermal coupling between the dust and the gas,
aside from the drag heating.

###  3.3 Numerical setup

We perform simulations with three different disc-to-star mass ratios
$M_{d}/M_{\star}=\\{0.05,0.1,0.2\\}$, three different cooling times
$\beta_{\text{cool}}=\\{8,10,15\\}$ and two different dust particle sizes, for
a total of 18 simulations. See below for further details. We first initialize
a gas-only disc around a solar mass star, with $R_{\text{in}}=0.25$au,
$R_{\text{out}}=25$au, and $\Sigma_{g}\propto R^{-1}$. We set the aspect ratio
so that $Q_{\text{ext}}=2$ initially and because of the cooling, it decreases,
eventually reaching $Q=1$. The shock viscosity coefficients
$\alpha^{\text{AV}}=0.1$, $\beta^{\text{AV}}=0.2$ as in Rice et al. (2005). As
a test, we performed a simulation with $\beta^{\text{AV}}=2$, and we did not
find any differences in terms of the relevant quantities in this work. The
self-gravity of both gas and dust is taken into account, as well as the dust
back-reaction. We performed simulations at two different numerical
resolutions: the standard runs are performed with $N_{g}=10^{6}$ and the high
resolution ones with $N_{g}=2\times 10^{6}$. In both cases, $N_{d}=N_{g}/5$,
where $N_{g}$ and $N_{d}$ are the number of gas and dust particles
respectively. We verify that the results are consistent with the different
resolutions.

We let the system evolve for an outer thermal time
($\beta_{\text{cool}}\Omega^{-1}$ at the outer radius), and then we add dust
particles and evolve for a further 5 outer dynamical times (i.e. $5\times
10^{3}$ inner dynamical times). The dust particles are added proportional to
the gas distribution.333We benchmarked our simulations with the ones of Rice
et al. (2004), as they started with an initial uniform dust distribution, with
a fixed thickness. We obtained the same results, since dust trapping is very
efficient in these systems. Something that should be pointed out is that Rice
et al. (2004) did not account for self-gravity acting on the solid particles.
Distributing dust particles proportional to the gas is a valid assumption only
for St <10\. Shi et al. (2016) and Baehr & Zhu (2021) showed that for
uncoupled particles gravitational interactions and stirring become quite
relevant, and hence dust distribution is closer to a uniform one. This is
particularly clear for $\text{St}\sim 100$, where. In our simulations of large
dust grains, although the initial distribution is proportional to the gas one,
we observe that the spiral is less prominent compared to smaller grains.

Dust back-reaction and dust self-gravity are always taken into account. Since
the aim of the work is to study the effect of the drag force in GI, we use two
different dust sizes: a larger one, to reproduce weakly coupled solid
particles, and a smaller one, to study marginally coupled particles. Since the
disc mass is different across our set of simulation, we chose to adapt the
particles’ size in order to obtain the same Stokes number distribution. To do
so, we computed the radially averaged Stokes number as a function of the
particle size for different disc-to-star mass ratio, starting from the initial
conditions of the gas disc, taking into account the transition between Epstein
and Stokes regime. We decided to choose the small particles’ size so that the
radially averaged Stokes number $\langle\text{St}\rangle=8$, and for the large
ones $\langle\text{St}\rangle=40$: in this way we effectively cover a Stokes
number range from 1 to 10 with smaller grains, and from 10 to 100 for larger
ones. One exception is that the small dust particles for the highest disc-to-
star mass ratio simulations are chosen to have $\langle\text{St}\rangle=16$
for computational reasons. In addition, as we will show in the next section, a
higher disc-to-star mass ratio makes the dust unstable, and collapse can
happen: for that reason, dust in simulations S16, S17 and S18 is evolved only
for an outer dynamical time. In every simulation, dust intrinsic density is
fixed $\rho_{0}=$5\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$$.

Self-gravitating discs may be more radially extended than our models, which
can be rescaled. If we rescale the outer radius of a factor $\lambda$, how
does the dust particles’ size need to be rescaled? Since
$\text{St}\propto{s}/{\Sigma}\propto sR_{\text{out}}^{2}$, if we change the
outer radius according to $R_{\text{out}}^{\prime}=\lambda R_{\text{out}}$, in
order to maintain the same Stokes number, the corresponding rescaling for dust
particles size should be $s^{\prime}=\lambda^{-2}s$. Hence, if we consider a
larger disc with $R_{\text{out}}^{\prime}=10R_{\text{out}}=250$ au, the
corresponding particle sizes should be rescaled as $s^{\prime}=s/100$. Dust
properties are summarized in Table 1, where we also include the rescaled dust
particles’ size. Snapshots of the hydrodynamical simulations are shown in
Figures 1, 2, 3.

Table 1: Parameters of simulations: disc-to-star mass ratio $M_{d}/M_{\star}$, cooling factor $\beta_{\text{cool}}$, size of dust particles $s$, average Stokes number $\langle\text{St}\rangle$ and corresponding dust particles size in a 10 times bigger disc $s_{10}$. Simulation | $M_{d}/M_{\star}$ | $\beta_{\text{cool}}$ | s [cm] | $\langle\text{St}\rangle$ | s10 [cm]
---|---|---|---|---|---
S1 | 0.05 | 8 | 300 | 40 | 3
S2 | 0.05 | 10 | 300 | 40 | 3
S3 | 0.05 | 15 | 300 | 40 | 3
S4 | 0.05 | 8 | 60 | 8 | 0.6
S5 | 0.05 | 10 | 60 | 8 | 0.6
S6 | 0.05 | 15 | 60 | 8 | 0.6
S7 | 0.1 | 8 | 600 | 40 | 6
S8 | 0.1 | 10 | 600 | 40 | 6
S9 | 0.1 | 15 | 600 | 40 | 6
S10 | 0.1 | 8 | 120 | 8 | 1.2
S11 | 0.1 | 10 | 120 | 8 | 1.2
S12 | 0.1 | 15 | 120 | 8 | 1.2
S13 | 0.2 | 8 | 1500 | 40 | 15
S14 | 0.2 | 10 | 1500 | 40 | 15
S15 | 0.2 | 15 | 1500 | 40 | 15
S16 | 0.2 | 8 | 600 | 16 | 6
S17 | 0.2 | 10 | 600 | 16 | 6
S18 | 0.2 | 15 | 600 | 16 | 6
Figure 1: Dust dynamics in gas spiral arms: large and small dust particles
surface density for different cooling factor and $M_{d}/M_{\star}=0.05$.
Figure 2: Dust dynamics in gas spiral arms: large and small dust particles
surface density for different cooling factor and $M_{d}/M_{\star}=0.1$. Figure
3: Dust dynamics in gas spiral arms: large and small dust particles surface
density for different cooling factor and $M_{d}/M_{\star}=0.2$.

##  4 Analysis and Results

In this section, we present the analysis of the numerical simulations. Since
the simulations use a two-fluid approach (Laibe & Price, 2012a, b), gas and
dust particles are treated as two different sets of particles, thus they
occupy different locations and carry their own physical information. In order
to obtain properties that depend on both gas and dust, such as dust to gas
ratio, or Stokes number, we interpolate gas properties to the location of dust
particles. In addition, since dust is modelled as a pressureless fluid, it has
no internal energy and no thermal sound speed. However, stirring phenomena
induce a velocity dispersion onto dust particles (Youdin & Lithwick, 2007): to
obtain this quantity, we compute it with an SPH interpolation over
neighbouring dust particles, via

$c_{d,i}^{2}=\sum_{j=1}^{N_{\text{neigh
}}}m_{j}\frac{\left(v_{d,i}-v_{d,j}\right)^{2}}{\rho_{j}}W_{ij}\left(h_{i}\right).$
(35)

For our analysis, we will mainly focus our attention on the dust to gas ratio
$\epsilon$, the relative temperature between gas and dust
$\xi=(c_{d}/c_{g})^{2}$, the Stokes number St, the cooling factor
$\beta_{\text{cool}}$ and the disc-to-star mass ratio $M_{d}/M_{\star}$. These
parameters are related to different physical phenomena: the dust to gas ratio
and the relative temperature trace dust trapping and dust excitation
respectively, the cooling factor is linked to the gas spiral amplitude (and
hence to the strength of GI), the Stokes number determines the power of the
aerodynamical coupling and the disc-to-star mass ratio is connected to the
spiral morphology.

###  4.1 Dust trapping: the $\epsilon$ parameter

Figure 4: Distribution of the dust to gas ratio for different values of
cooling factor (top panel) and disc to star mass ratio (bottom panel). The
simulations shown in these plots are S1,S2,S3,S7,S13.

GI spiral arms trap dust particles (Dipierro et al., 2015a), since they are
both pressure maxima and gravitational potential minima. We use the
simulations to quantify how this phenomenon depends on the model parameters.
Figure 4 shows the distribution of the dust to gas ratio for different
$\beta_{\text{cool}}$ and $M_{d}/M_{\star}$, for a set of simulations with
large dust particles. The initial value of the dust to gas ratio is $10^{-2}$.
In Figure 4 the higher tail of the distributions reaches values of $\gtrsim
10^{-1}$, implying that dust concentration by up to approximately an order of
magnitude is happening. Figure 5 shows a comparison between the dust to gas
ratio distributions of large and small particles. As expected, dust
concentration in spiral arms is stronger for smaller particles, and it can
approach values of the order of unity, since their aerodynamical coupling with
gas is stronger. The strength of dust trapping is determined by both the
aerodynamic coupling between gas and dust and the gravitational potential of
gas spiral arms. The combined effect of gravitational and drag interaction is
maximised when $\text{St}\simeq Q\simeq 1$ (Baehr & Zhu, 2021), thus, in our
simulations, smaller particles reach higher values of the dust to gas ratio.
No particular correlations are found between $\epsilon$ and the disc to star
mass ratio, while there is a slight dependence on the cooling factor. In order
to understand this relationship, Figure 6 shows the quantity
$\delta\Sigma/\Sigma_{0}$ for gas (orange dots) and dust (large particles -
blue dots, small particles - green dots) as a function of the cooling factor
$\beta_{\text{cool}}$. The quantity $\Sigma_{0}$ is the azimuthally averaged
surface density at a fiducial radius of 10 au, and the quantity $\delta\Sigma$
is its standard deviation. For the gas, it is known that
$\delta\Sigma_{g}/\Sigma_{g0}\propto\beta_{\text{cool}}^{-1/2}$ (Cossins et
al., 2009), and we recover this behaviour in our simulations. For the dust,
the situation is different. We do not find any evident correlation between the
density contrast and the cooling factor. So, if we assume that

$\frac{\delta\Sigma_{g}}{\Sigma_{g0}}\propto\beta_{\text{cool}}^{-1/2},\quad\frac{\delta\Sigma_{d}}{\Sigma_{d0}}=\text{const},$
(36)

the ratio between these two quantities is

$\frac{\delta\Sigma_{g}}{\Sigma_{g0}}\frac{\Sigma_{d0}}{\delta\Sigma_{d}}=0.01\epsilon^{-1}\propto\beta_{\text{cool}}^{-1/2},$
(37)

where $\Sigma_{d0}/\Sigma_{g0}=1/100$, meaning that
$\epsilon\propto\beta_{\text{cool}}^{1/2}$. This is the positive correlation
we found before. Why do we not find any evident correlation between the dust
density contrast and the cooling factor? Physically, the dust experiences the
effect of the gas cooling through gravitational and drag forces. When
$\text{St}\ll 1$, dust and gas particles are indistinguishable, and so
$\delta\Sigma_{g}/\Sigma_{g0}=\delta\Sigma_{d}/\Sigma_{d0}\propto\beta_{\text{cool}}^{-1/2}$.
For higher Stokes number, the drag force is weaker, and the dust is less
influenced by the gas cooling. In this case, we expect the relationship
between $\delta\Sigma_{d}/\Sigma_{d0}$ and $\beta_{\text{cool}}$ to be flatter
than $\beta_{\text{cool}}^{-1/2}$: if this condition is respected, the
correlation between $\epsilon$ and $\beta_{\text{cool}}$ will be positive. In
general, $\delta\Sigma_{d}/\Sigma_{d0}$ is a function of both the cooling
factor and the Stokes number.

Figure 5: Comparison of the distribution of dust to gas ratio of large (orange
line) and small (blue line) dust particles. The simulations shown in this plot
are S1 and S4. Figure 6: Density contrast $\delta\Sigma/\Sigma_{0}$ of gas
(orange), large dust grains (blue) and small dust grains (green) as a function
of the cooling factor, for $M_{d}/M_{\star}=0.05$.

###  4.2 Dust excitation: the $\xi$ parameter

Figure 7: Distribution of the relative temperature for different values of
cooling factor (top panel) and disc to star mass ratio (bottom panel). The
simulations shown in these plots are S1,S2,S3,S7,S13.

To investigate dust excitation by spiral arms, we study how the relative
temperature $\xi=(c_{d}/c_{g})^{2}$ varies as a function of the simulation
parameters. Figure 7 shows the distribution of the dust relative temperature
$\xi$ for different values of $\beta_{\text{cool}}$ and $M_{d}/M_{\star}$ for
a set of simulations with large dust particles. We observe that for the
simulations in which dust collapse is not happening, the dust dispersion
velocity reaches very quickly ($<1$ outer orbital time) a steady value. The
relative temperature shows a negative correlation with both the disc-to-star
mass ratio and the cooling factor. For the disc-to-star mass ratio, this can
be understood by considering the relationship with the spiral morphology:
$M_{d}/M_{\star}$ is inversely proportional to the azimuthal wavenumber $m$
(Cossins et al., 2009), hence massive discs have fewer spiral arms. Since dust
is excited because of “spiral kicks” (Walmswell et al., 2013), the lower the
number of spiral arms, the less the dust is excited. For the cooling factor,
the negative correlation can be understood in two ways. First, the cooling
rate $\beta_{\text{cool}}$ is linked to the amplitude of the spiral
perturbation according to eq. 36. Since gas spiral arms excite dust particles
by kicking them every passage, the higher is the perturbation, the stronger is
the kick, and so the excitation. Second, in gravito-turbulent regime,
transport of angular momentum is driven by the spiral perturbation, and it is
possible to define an $\alpha-$viscosity coefficient related to the cooling
rate (eq. 33). The height of a dust layer is determined by the interaction
with the gas and by the vertical diffusion. In the hypothesis that the
vertical diffusion coefficient is equal to the azimuthal one, we can obtain
the height of the dust layer as

$H_{d}=H_{g}\sqrt{\frac{\alpha}{\alpha+\text{St}}}.$ (38)

If we assume that $\alpha=\alpha_{\text{GI}}$, the dust layer height, and thus
the dust dispersion velocity, is inversely proportional to
$\beta_{\text{cool}}$.

Figure 8 compares the relative temperature distributions of large and small
dust particles. Small particles are colder than the large ones, and their
relative temperature is almost completely enclosed in the interval
$\xi_{\text{small}}\in[10^{-2},1]$, meaning that their random motions are
subsonic. On average, the distribution of $\xi_{\text{small}}$ is shifted by
one order of magnitude compared to $\xi_{\text{large}}$. This behaviour is in
agreement with what Booth & Clarke (2016) found: larger particles tend to be
dynamically hotter, because the kicks of the gas spiral are more effective,
while if the coupling with the gas is stronger, the kicks are damped because
of the drag force.

In principle, the observed trend would imply that increasing
$\beta_{\text{cool}}$ would lead to an arbitrarily thin dust layer, eventually
causing gravitational collapse. This would be a direct analogue of the
classical Goldreich & Ward (1973) mechanism for planetesimal formation in a
(weakly) self-gravitating context, and in a more realistic model it likely be
limited in a similar way by the excitation of shear turbulence (Cuzzi et al.,
1993). However, there is an upper limit for $\beta_{\text{cool}}$ set by the
value above which the transfer of angular momentum would be driven by some
process other than gravitational instability. To compute an estimate of the
maximum cooling time, we require $\alpha_{\text{GI}}$ to be larger than
$10^{-3}$. Observations of protoplanetary discs that are not expected to be
self-gravitating suggest that this is a reasonable upper limit to the strength
of turbulence (Flaherty et al., 2017). Thus, we obtain

$\beta_{\text{cool}}^{\text{max}}=400\left(\frac{\alpha}{10^{-3}}\right)^{-1},$
(39)

that corresponds to a minimum density perturbation

$\left.\frac{\delta\Sigma_{g}}{\Sigma_{g0}}\right|_{\text{min}}=0.05\left(\frac{\alpha}{10^{-3}}\right)^{1/2}.$
(40)

Table 2: Table that summarizes the correlations between $\epsilon,$ $\xi$ and $\beta_{\text{cool}},$ $M_{d}/M_{\star}$ and St. | $\beta_{\text{cool}}$ | $M_{d}/M_{\star}$ | St
---|---|---|---
$\epsilon$ | Positive | None | Negative
$\xi$ | Negative | Negative | Positive

Table 2 summarizes the relationships we have discussed in these paragraphs. A
broader collection of histograms can be found in Appendix C.

Figure 8: Comparison of the distribution of relative temperature of large
(orange line) and small (blue line) dust particles. The simulations shown in
this plot are S1 and S4.

Booth & Clarke (2016) studied the relationship between the dust excitation and
both the cooling and the Stokes number. They found that
$c_{d}\propto\beta^{-1/2}\text{St}^{1/2}v_{k}$, where $v_{k}$ is the Keplerian
speed. To compare with Booth & Clarke (2016), we use our data to reproduce
Figures 7 and 13 of their paper, that show a relationship between the dust
velocity dispersion and $\beta_{\text{cool}}$ and St. To do so, we divided the
particles into equally spaced intervals of Stokes number and, for each
particle, we computed the mean value of $c_{d}/c_{g}=\sqrt{\xi}$. The
comparison is shown in Figure 9, where we show the results of simulations with
$M_{d}/M_{\star}=0.05$, for both standard and high resolution cases. The
previously derived relationships with the Stokes number (the left panel) and
with the cooling factor (the right panel) are well recovered. Using our two
fluid algorithm, it is too computationally expensive to analyse properly the
case of $\text{St}<1$. However, in this case we expect that as the
aerodynamical coupling with the gas is stronger, $c_{d}/c_{g}$ should
increase, eventually reaching $c_{d}=c_{g}$ for $\text{St}\to 0$. This growth
for $\text{St}<1$ has been shown by Booth & Clarke (2016).

Figure 9: Comparison with Booth & Clarke (2016). The left panel shows how dust
dispersion velocity depends on the Stokes number, for different values of
cooling factor, compared to the expected relationship
$\propto\text{St}^{1/2}$. The right panel, shows how dust dispersion velocity
depends on the cooling factor, for different values of Stokes number, compared
to the expected relationship $\propto\beta_{\text{cool}}^{-1/2}$. All the
simulations shown in these plots have $M_{d}/M_{\star}=0.05$. The diamonds are
the values obtained from the standard resolution simulations, while the
squares from the high resolution ones.

###  4.3 Two fluid instability

In this section, we apply the two fluid instability theory presented in
Section 2. The gas-only and gas-and-dust models for gravitational instability
have both been developed within a linear framework; hence, in principle, the
quantities $\epsilon$, $\xi$ and St should be evaluated in the unperturbed
state. However, in this work, we are evaluating them in the perturbed one.
Although not completely self-consistent, it gives us an idea of the most
unstable regions of the disc. Figure 10 shows the distribution of large (blue)
and small (orange) dust particles in the $(\xi,\epsilon)$ diagram: the black
lines corresponds to the dust driven GI threshold for $\text{St}=\infty$ (eq.
16, solid line) and for $\text{St}=0.5$ (eq. 17, dashed line). We choose
$\text{St}=0.5$ as a minimum value since the number of particles with Stokes
number lower than this is negligible. The particles above the region are in a
dust driven GI regime. We find that the number of small particles for which
the instability is dust driven is greater compared to large ones: this is
because small particles have both larger dust to gas ratio and lower
dispersion velocity. In addition, the number of dust driven particles
increases with the cooling factor and with the disc to star mass ratio, as
already discussed in previous sections. To understand the spatial location of
these particles in the disc, Figure 11 shows the particles that satisfy
condition 16 superimposed on the total density map. As expected, the most
unstable regions of the disc are not randomly distributed, but correspond with
the spiral arms.

Figure 12 shows the value of the Jeans mass for $M_{d}/M_{\star}=0.05$ as a
function of the Stokes number. The curve can be divided into two parts: for
$\text{St}\in(0,5)$, the Jeans mass is decreasing with the Stokes number,
reaching its minimum at about $\text{St}\sim 3$. This happens because for
$\text{St}\to 0$, the particles are indistinguishable, and the instability is
gas driven. It is important to notice that the number of particles with small
Stokes number is low, hence the binning in Stokes number presents a
considerable scatter. By increasing the Stokes number, $\xi$ decreases and the
two fluids behave more and more differently. When the Stokes number is
approximately 1, the relative temperature is a minimum and the dust to gas
ratio is high, so the instability becomes dust driven. Otherwise, for
$\text{St}>5$, the Jeans mass increases with the Stokes number. Indeed, the
relative temperature increases, and the dust to gas ratio decreases. Hence,
the system transitions from dust into gas driven instability, again,
eventually reaching the gas-only component model value.

Figure 10: Distribution of large (blue) and small (orange) dust particles in
the $(\epsilon,\xi)$ diagram, for different values of disc-to-star mass ratio
and cooling factor. The black lines corresponds to the dust driven GI
threshold for $\text{St}=\infty$ (eq. 16, solid line) and for $\text{St}=0.5$
(eq. 17, dashed line). We choose $\text{St}=0.5$ as a minimum value since the
number of particles with Stokes number lower than this is negligible. The
particles above the region are in a dust driven GI regime. Figure 11: Total
density maps for small dust particles simulations
(S4,S5,S6,S10,S11,S12,S16,S17,S18 in order). Yellow dots correspond to
particles for which the instability is dust driven. Figure 12: Jeans mass
given by eq. (20) evaluated at each dust particle location as a function of
the Stokes number, for $M_{d}/M_{\star}=0.05$ and $\beta_{\text{cool}}=8$.

##  5 Discussion

###  5.1 Dust collapse

The simulations with $M_{d}/M_{\star}=0.2$ and small dust particles do not
reach 5 outer orbits, since the simulation stops due to the onset of dust
collapse. This happens because the stopping time of collapsing dust particles
becomes smaller than the time step of the code. Indeed, the stopping time is
inversely proportional to the total density
$\rho_{\text{tot}}=\rho_{g}+\rho_{d}$, and since the dust density is
increasing, because of the collapse, the stopping time tends to zero. The top
panel of Figure 13 shows the maximum dust density as a function of time for
small and large dust particles, in the run with $M_{d}/M_{\star}=0.2$ and
$\beta=15$. While large dust particles reach a quasi-steady state, the small
particle density exponentially increases in the first orbit. This is the
signature of dust collapse. This phenomenon is visible in simulations S16, S17
and S18, and happens only in the dust component. The bottom panel of Figure 13
shows a comparison between gas and small dust averaged density as a function
of time, for $M_{d}/M_{\star}=0.2$ and $\beta=15$. At $t=0$,
$\langle\rho_{d}\rangle=10^{-2}\langle\rho_{g}\rangle$. Whereas the gas
average density is constant with time, the dust density increases because of
dust trapping, eventually exceeding that of the gas. This means that any
clumps forming from this mechanism would be substantially made up of solids,
and would likely be identified with the rocky core of a giant planet. However,
in this work we do not want to characterize the outcome of this collapse,
which is a complex topic. Indeed, simulations of planetesimal collapse
(Nesvorný et al., 2021) show that a rotating self-gravitating cloud of dust
does not monolithically collapse, meaning that it is not possible to directly
equate the cloud mass with the planetary core one.

To identify and analyse dust clumps in more detail, we define the numerical
conditions that should be respected for a clump to be physical and not
affected by resolution. For a clump radius $r_{\text{clump}}$ the smoothing
length of the dust particle $h_{i}$ should be less than a fraction of the
clump radius, in order to be resolved. This condition translates into
$h_{i}<\eta r_{\text{clump}}$, where $\eta$ is less than unity. We take
$\eta=1/2$. Physically, a gravitationally bound clump is a collection of
particles whose thermal support does not balance the gravitational one. We
define the thermal and gravitational energy of particles inside
$r_{\text{clump}}$ as follows:

$e_{\text{th}}=\sum_{i}^{i\in r_{\text{clump}}}m_{i}c_{i}^{2},$ (41)
$e_{\text{gr}}=-\frac{1}{2}\sum_{i}^{i\in r_{\text{clump}}}\sum_{j=i}^{j\in
r_{\text{clump}}}m_{i}m_{j}[\varphi(r_{ij},h_{i})+\varphi(r_{ij},h_{j})],$
(42)

where $\varphi$ is the gravitational softening kernel and
$r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$. We used a cubic spline softening
kernel, and the detailed expression can be found in the appendix of Price &
Monaghan (2007). According to the virial theorem, if the force between any two
particles can be described in terms of a potential energy $\Phi\propto r^{n}$,
where $r$ is the distance between two particles, the equilibrium state
respects the following condition

$2\langle T\rangle=n\langle\Phi\rangle,$ (43)

where $T$ is the kinetic energy. In the case of gravitational interaction, the
virial theorem reads $\langle T\rangle/\langle\Phi\rangle=-1/2$. We define a
clump as a region of the space where the dimensionless quantity
$\alpha_{J}=-e_{\text{th}}/e_{\text{gr}}<1/2$. Then, in order to be sure that
the collapse is physical, and not artificial, we verify that $h_{g}<h_{d}$ in
the region where there is the dust clump. The last condition requires that the
resolution of the gas in the region where there is a possible clump should be
higher compared to the dust resolution.

To summarize, the conditions under which we define a dust clump are the
following

  * •

$h_{g,i}\leq h_{d,i}$,

  * •

$h_{d,i}<\eta r_{\text{clump}}$,

  * •

$\alpha_{J}<1/2$.

Table 3: Comparison between the expected and the observed Jeans Mass in the simulations where dust collapse happens. # | Expected $M_{J}$ [M⊕] | Simulation [M⊕]
---|---|---
S16 | 2 | 0.6
S17 | 2.2 | 1.
S18 | 4.5 | 3.2

In simulations S16, S17 and S18 there are particles that respect the previous
conditions, implying that dust collapse has happened. Table 3 shows the mass
of the clumps, obtained by summing the mass of each particle that is
gravitationally bound and the one obtained from analytical theory from Eq. 20.
The masses are broadly in agreement with analytic expectations. For S16, for
example, the derived Jeans mass is of the order of an Earth mass. In general,
the mass of the clump computed from the simulation is smaller compared to the
one expected from the analytical theory: this is not surprising, since with
the simulation we are only able to appreciate the initial phase of the
collapse. Indeed, as soon as the stopping time is smaller than the time step,
the simulation stops. To avoid this problem, one could decrease the time step
of the code, but it is computationally expensive. Otherwise, one could not
consider the dust density when computing the stopping time: this procedure has
been applied in previous works to increase the velocity of the simulations
(Poblete et al., 2019; Longarini et al., 2021) but this approximation is valid
only for small dust to gas ratios. While there this was justified, here this
is not possible, since dust is collapsing and the dust to gas ratio becomes
higher than unity. Finally, to study the early evolution of clumps with an SPH
code, it would be possible to simulate them as sink particles: so far, the
creation of dust sink particles is not possible in phantom.

Figure 13: Maximum dust density as a function of time for simulations S15
($M_{d}/M_{\star}=0.2,\beta=15$, large dust particles, blue line) and S18
($M_{d}/M_{\star}=0.2,\beta=15$, small dust particles, blue line).

####  5.1.1 Gas-dust coupling during the collapse

The bottom panel of Figure 13 shows that the collapse happens only in the dust
component, and the gas is not influenced. This could sound surprising: indeed,
in high density regions, we expect the stopping time (eq. 25) to be small. So,
why is the dust collapse not influencing the gas? The degree of coupling is
measured with the Stokes number, that compares the strength of the drag force
with the ones that are acting on the particle. Usually, in a flat
protoplanetary disc, the Stokes number is computed as the ratio of the
stopping time and the dynamical time, that is the typical timescale of a
particle orbiting around a central object at a distance $R$. However, in this
situation, the dust clump is driving the dynamics of the surrounding
particles, and hence we should compare the stopping time with the free fall
time in order to understand the degree of coupling of particles. The typical
timescale of the infall of a spherically-symmetric distribution of mass is

$t_{\text{ff}}=\sqrt{\frac{3\pi}{32G\rho}}.$ (44)

Comparing the stopping time and the free fall time in the simulations where
dust collapse is happening (Table 4), we obtain that the particles in the
collapsing region are uncoupled, since the ratio between the two timescales is
higher than one.

Table 4: Stopping, dynamical and free fall timescales for the simulations that show dust collapse # |  $t_{s}$ [yr] |  $t_{\text{dyn}}$ [yr] |  $t_{\text{ff}}$ [yr] | $t_{s}/t_{\text{dyn}}$ | $t_{s}/t_{\text{ff}}$
---|---|---|---|---|---
S16 | 2.7 | 7.9 | 0.4 | 0.3 | 6
S17 | 5.6 | 14.1 | 0.7 | 0.4 | 8
S18 | 7.2 | 17.7 | 1.6 | 0.4 | 4

It is possible to quantify the critical density a clump should reach so that
$t_{\text{ff}}<t_{\text{dyn}}$, that is

$\rho_{\text{crit}}=\frac{3\pi}{32}\frac{M_{\star}}{R^{3}}:$ (45)

when a clump reaches this density, the evolution of the surrounding particles
is determined by the clump, and not by the star anymore. From that point, the
aerodynamical coupling should be quantified by taking the ratio between the
stopping time and the free fall time. Hence, for $\rho<\rho_{\text{crit}}$,
$\text{St}\propto\rho^{-1}$, since the dynamical time does not depend on the
density. For $\rho>\rho_{\text{crit}}$, the scaling changes since the free
fall time depends on the density $\text{St}\propto\rho^{-1/2}$, and so does
the degree of coupling.

###  5.2 Gravitational instability in the context of protoplanetary disc
evolution

In this work, we focused on cooling driven GI, but it is also possible to
trigger it through infall. Kratter et al. (2010) found that in the infall-
driven case the strength of the spiral perturbation is controlled by two
dimensionless parameters, a thermal one, that relates the infall mass
accretion rate $\dot{M}_{\text{inf}}$ to the characteristic sound speed of the
disc, and a rotational one, that compares the relative strength of rotation
and gravity in the core. Obviously, the higher the accretion rate, the
stronger is the spiral perturbation in the disc. We can naively associate the
dust evolution we have modelled in the limit of fast (slow) cooling with high
(low) infall rate. This association is expected to be qualitatively correct,
however a detailed comparison between these two regimes is given by Kratter &
Lodato (2016).

By studying the non-linear evolution of gas and dust in protoplanetary discs,
we found that the instability conditions for the two components are different.
It is well established that spiral fragmentation occurs in fast cooling gas
discs, and it is possible to define a critical value of $\beta_{\text{cool}}$
below which fragmentation occurs. Simulations of cooling-driven fragmentation
(Gammie, 2001; Rice et al., 2005; Lodato & Clarke, 2011; Meru & Bate, 2012)
currently suggest that $\beta_{\text{min}}\simeq 3$ (Deng et al., 2017). For
the dust, Booth & Clarke (2016) found that dust becomes more unstable for
higher $\beta_{\text{cool}}$, and we confirm this trend with 3D simulations.
The differing behaviour of gas and dust suggests an interesting evolution in
the outcome of gravitational instability within protoplanetary discs. During a
first stage, at the beginning of the disc lifetime, we expect a very massive
disc system characterized by strong GI, caused by the high infall rate from
the molecular cloud. If conditions allow gas fragmentation, because of the
high Jeans mass, low mass stellar companions can be formed (Kratter et al.,
2010). A second stage is characterized by a less massive disc, with lower
infall rate, or equivalently longer cooling time. If conditions during this
epoch trigger gravitational instability, it will lead to dust-driven
fragmentation that could be responsible for the formation of rocky cores of
giant planets. Then, in the third stage, there is a protostar surrounded by a
planet hosting disc, characterized by substructures such as gaps, rings or
planetary spirals. In this stage, GI is not effective anymore because the disc
mass is small and the transport of angular momentum is controlled by disc
winds (Tabone et al., 2022) or other, non-self-gravitating, sources of
turbulence (Lesur et al., 2022). A schematic view of these stages is given by
figure 14.

###  5.3 Application to an actual case: HL Tau

HL Tau is a young ($<1\text{Myr}$) protostellar system that shows axisymmetric
structures (gaps and rings) in dust continuum emission (ALMA Partnership et
al., 2015). The origin of rings and gaps is usually attributed to planet disc
interaction (Lin & Papaloizou, 1986), with Dipierro et al. (2015b) finding
that three protoplanets with masses $M_{p1}=$61\text{\,}\mathrm{M}_{\oplus}$$,
$M_{p2}=$83\text{\,}\mathrm{M}_{\oplus}$$,
$M_{p3}=$170\text{\,}\mathrm{M}_{\oplus}$$ at $R_{p1}=13.2$au, $R_{p2}=32.3$au
and $R_{p3}=68.8$au best match the observations. The formation of super-Earth
mass planets at large radii, in such a young system, is a challenge for core
accretion theory. If planets can form via the dust-induced collapse mechanism
we have discussed in this work, HL Tau is a plausible example of what the
resulting planetary system might look like. We note that the inferred trend of
increasing mass with radius can be interpreted within our model. In a
simplistic way, we assume that $\text{St}\propto\Sigma^{-1}\propto R$, and
that $\xi\propto\beta_{\text{cool}}^{-1}\text{St}\propto R^{4.5}R=R^{5.5}$,
where we have supposed a realistic cooling law, according to which
$\beta_{\text{cool}}\propto R^{-4.5}$ (Rafikov, 2005; Clarke & Lodato, 2009).
Supposing that the dust to gas ratio is constant, the Jeans mass of the gas-
and-dust fluid model will be an increasing function of the radius, since the
relative temperature is smaller for particles closer to the central star. So,
within this hypothesis, we expect that the mass of the dust driven GI
protoplanets will be an increasing function of the radius. We want to point
out that this is just a qualitative argument. Indeed, to thoroughly
investigate whether HL Tau planets can be formed through dust driven GI, one
should properly model the system. In addition, planet migration and planet
accretion should also be considered.

Figure 14: Scheme of the three stages of disc lifetime

##  6 Conclusions

Self-gravitating gas discs may be ubiquitous during the Class 0/I phases of
YSO evolution (e.g. Lin & Pringle, 1990; Xu & Kunz, 2021). Observations
suggest that the self-gravity can in some systems be strong enough as to
trigger fragmentation (Tobin et al., 2016), while in other cases, such as the
massive disc of the more evolved IM Lup system (Lodato et al., 2023) the
instability is expected to be more gentle. If that fragmentation occurs in the
gas, as in the L1448 IRS3B system, the outcome is typically star or brown
dwarf formation. Lower mass objects can be formed if the fragmenting fluid is
instead the solid component of a two-fluid self-gravitating disc, given disc
conditions that allows the collisional growth of dust to small macroscopic
dimensions while the disc remains self-gravitating. Analytic estimates suggest
that planetary cores of $\sim 1-10\ M_{\oplus}$ could form from this mechanism
at large orbital radius, with properties that could be identified with the
population of ALMA-inferred disc-embedded planets (Andrews, 2020).

In this paper, we presented the results of SPH simulations of gas and dust in
protoplanetary discs, studying the role of aerodynamic coupling in the context
of gravitational instability. We analysed our results in the framework of two
fluid gravitational instability, and compared our findings with previous
numerical works, finding generally good agreement.

Our main results can be summarized as follows:

  1. 1.

We studied the relationship between the dust to gas ratio $\epsilon$, the
relative temperature between gas and dust $\xi$ and the cooling factor
$\beta_{\text{cool}}$, the disc to star mass ratio and the Stokes number. We
found that the dust to gas ratio increases with the cooling factor and
decreases with the Stokes number, and that the relative temperature increases
with the Stokes number and decreases with the cooling factor and the disc to
star mass ratio. It is possible to explain these relationships by considering
the interaction between dust particles and gas spiral arms. We compared our
findings with Booth & Clarke (2016) and found good agreement.

  2. 2.

We investigated the role of dust in gravitational instability, and found that
the most unstable regions of the disc are the spiral arms, where the
instability tends to be dust driven. In addition, we studied the relationship
between the theoretical Jeans mass and the Stokes number. We found that the
Jeans mass — when the instability is dust driven — can reach values of the
order of the Earth mass.

  3. 3.

We observe three cases of dust collapse in our set of simulation, which occur
(as expected) in the high disc mass models. The values of the clump masses
obtained numerically are close to those predicted by linear theory.

In applying our results to the possibility of early planet formation in Class
0/I discs, the main prerequisite is the requirement that dust is able to grow
via coagulation to a large enough Stokes number, with a top-heavy particle
mass function, in a short enough time. Our simulations that explicitly
exhibited dust collapse had solid particles with an average Stokes number
$\langle{\rm St}\rangle=16$, which would correspond (rescaling our simulations
to a disc size of $R_{\text{out}}=250$au) to particles with sizes between a
few cm and a few metres. Fragmentation and radial drift pose barriers to
growth to the required sizes (Birnstiel et al., 2010), and further work is
needed to assess whether there are circumstances where the required Stokes
numbers can be reached. Simulations with a constant Stokes number, along with
runs with ${\rm St}\sim 1$ (a regime which is numerically difficult to access
using our code), would also help to better define the regime where dust can
fragment in a gas disc that is itself stable against fragmentation.

In an evolutionary context, our results imply that planet formation — if it is
able to occur via the mechanism of dust-dominated gravitational disc
instability — is likely to occur toward the end of the self-gravitating phase.
This is when the competing effects of particle trapping and particle
excitation are jointly most favourable for collapse. The masses and orbital
radii of the planets formed via dust collapse are qualitatively in agreement
with those inferred for the HL Tau system, and we speculate that they may form
from the collapse of solids in the spiral arms of a formerly self-gravitating
protostellar disc.

## Acknowledgements

The authors thank the referee for useful comments and suggestions that
significantly improved the quality of the work. In this work we used splash to
create figures of the hydrodynamical simulations (Price, 2007). This project
and the authors have received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement N. 823823 (DUSTBUSTERS RISE project). CL acknowledges support from
Fulbright Commission through VRS scholarship. PJA acknowledges support from
NASA TCAN award 80NSSC19K0639, and from award 644616 from the Simons
Foundation. DJP acknowledges Australian Research Council funding via
DP220103767 and DP180104235. The authors thank Sahl Rowther, Benedetta
Veronesi, Cathie Clarke and Richard Booth for useful discussions. We thank
Stony Brook Research Computing and Cyberinfrastructure, and the Institute for
Advanced Computational Science at Stony Brook University for access to the
SeaWulf computing system, supported by National Science Foundation grant #
1531492.

## Data Availability

The data presented in this article will be shared on reasonable request to the
corresponding author.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, ApJ, 808, L3
  * Adams et al. (1987) Adams F. C., Lada C. J., Shu F. H., 1987, ApJ, 312, 788
  * Aly et al. (2021) Aly H., Gonzalez J.-F., Nealon R., Longarini C., Lodato G., Price D. J., 2021, MNRAS, 508, 2743
  * Andrews (2020) Andrews S. M., 2020, ARA&A, 58, 483
  * Ayliffe et al. (2012) Ayliffe B. A., Laibe G., Price D. J., Bate M. R., 2012, MNRAS, 423, 1450
  * Baehr & Zhu (2021) Baehr H., Zhu Z., 2021, ApJ, 909, 135
  * Bertin & Romeo (1988) Bertin G., Romeo A. B., 1988, A&A, 195, 105
  * Béthune et al. (2021) Béthune W., Latter H., Kley W., 2021, A&A, 650, A49
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Booth & Clarke (2016) Booth R. A., Clarke C. J., 2016, MNRAS, 458, 2676
  * Boss (1997) Boss A. P., 1997, Science, 276, 1836
  * Ceppi et al. (2022) Ceppi S., Cuello N., Lodato G., Clarke C., Toci C., Price D. J., 2022, MNRAS, 514, 906
  * Chambers (2021) Chambers J., 2021, ApJ, 914, 102
  * Clarke (2009) Clarke C. J., 2009, MNRAS, 396, 1066
  * Clarke & Lodato (2009) Clarke C. J., Lodato G., 2009, MNRAS, 398, L6
  * Cossins et al. (2009) Cossins P., Lodato G., Clarke C. J., 2009, MNRAS, 393, 1157
  * Cummins et al. (2022) Cummins D. P., Owen J. E., Booth R. A., 2022, MNRAS, 515, 1276
  * Cuzzi et al. (1993) Cuzzi J. N., Dobrovolskis A. R., Champney J. M., 1993, Icarus, 106, 102
  * Deng et al. (2017) Deng H., Mayer L., Meru F., 2017, ApJ, 847, 43
  * Dipierro et al. (2015a) Dipierro G., Pinilla P., Lodato G., Testi L., 2015a, MNRAS, 451, 974
  * Dipierro et al. (2015b) Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015b, MNRAS, 453, L73
  * Dong et al. (2015) Dong R., Zhu Z., Whitney B., 2015, ApJ, 809, 93
  * Fassio & Probstein (1970) Fassio F., Probstein R. F., 1970, AIAA Journal, 8, 772
  * Flaherty et al. (2017) Flaherty K. M., et al., 2017, ApJ, 843, 150
  * Fulk & Quinn (1996) Fulk D. A., Quinn D. W., 1996, Journal of Computational Physics, 126, 165
  * Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
  * Gibbons et al. (2012) Gibbons P. G., Rice W. K. M., Mamatsashvili G. R., 2012, MNRAS, 426, 1444
  * Gibbons et al. (2014) Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2014, MNRAS, 442, 361
  * Gibbons et al. (2015) Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2015, MNRAS, 453, 4232
  * Goldreich & Ward (1973) Goldreich P., Ward W. R., 1973, ApJ, 183, 1051
  * Guilera et al. (2020) Guilera O. M., Sándor Z., Ronco M. P., Venturini J., Miller Bertolami M. M., 2020, A&A, 642, A140
  * Jiang & Ormel (2022) Jiang H., Ormel C. W., 2022, arXiv e-prints, p. arXiv:2207.13002
  * Jog & Solomon (1984) Jog C. J., Solomon P. M., 1984, ApJ, 276, 114
  * Johansen et al. (2009) Johansen A., Youdin A., Klahr H., 2009, ApJ, 697, 1269
  * Kratter & Lodato (2016) Kratter K., Lodato G., 2016, ARA&A, 54, 271
  * Kratter et al. (2010) Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 708, 1585
  * Laibe & Price (2012a) Laibe G., Price D. J., 2012a, MNRAS, 420, 2345
  * Laibe & Price (2012b) Laibe G., Price D. J., 2012b, MNRAS, 420, 2365
  * Lesur et al. (2022) Lesur G., et al., 2022, arXiv e-prints, p. arXiv:2203.09821
  * Levison et al. (2010) Levison H. F., Thommes E., Duncan M. J., 2010, AJ, 139, 1297
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Lin & Pringle (1990) Lin D. N. C., Pringle J. E., 1990, ApJ, 358, 515
  * Lodato & Clarke (2011) Lodato G., Clarke C. J., 2011, MNRAS, 413, 2735
  * Lodato & Rice (2004) Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630
  * Lodato et al. (2019) Lodato G., et al., 2019, MNRAS, 486, 453
  * Lodato et al. (2023) Lodato G., et al., 2023, MNRAS, 518, 4481
  * Longarini et al. (2021) Longarini C., Lodato G., Toci C., Aly H., 2021, MNRAS, 503, 4930
  * Longarini et al. (2023) Longarini C., Lodato G., Bertin G., Armitage P. J., 2023, MNRAS, 519, 2017
  * Meru & Bate (2012) Meru F., Bate M. R., 2012, MNRAS, 427, 2022
  * Monaghan & Kocharyan (1995) Monaghan J. J., Kocharyan A., 1995, Computer Physics Communications, 87, 225
  * Morbidelli (2020) Morbidelli A., 2020, A&A, 638, A1
  * Nesvorný et al. (2021) Nesvorný D., Li R., Simon J. B., Youdin A. N., Richardson D. C., Marschall R., Grundy W. M., 2021, PSJ, 2, 27
  * Paardekooper & Mellema (2004) Paardekooper S. J., Mellema G., 2004, A&A, 425, L9
  * Paardekooper & Mellema (2006) Paardekooper S. J., Mellema G., 2006, A&A, 453, 1129
  * Pinilla et al. (2012) Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L., Testi L., Natta A., 2012, A&A, 538, A114
  * Poblete et al. (2019) Poblete P. P., Cuello N., Cuadra J., 2019, MNRAS, 489, 2204
  * Pollack et al. (1996) Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62
  * Price (2007) Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159
  * Price & Laibe (2020) Price D. J., Laibe G., 2020, MNRAS, 495, 3929
  * Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  * Price et al. (2018) Price D. J., et al., 2018, Publ. Astron. Soc. Australia, 35, e031
  * Rafikov (2005) Rafikov R. R., 2005, ApJ, 621, L69
  * Rafikov (2009) Rafikov R. R., 2009, ApJ, 704, 281
  * Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
  * Rice et al. (2004) Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2004, MNRAS, 355, 543
  * Rice et al. (2005) Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56
  * Riols et al. (2020) Riols A., Roux B., Latter H., Lesur G., 2020, MNRAS, 493, 4631
  * Rosenthal & Murray-Clay (2018) Rosenthal M. M., Murray-Clay R. A., 2018, ApJ, 864, 66
  * Safronov (1969) Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka.
  * Segura-Cox et al. (2020) Segura-Cox D. M., et al., 2020, Nature, 586, 228
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Sheehan & Eisner (2018) Sheehan P. D., Eisner J. A., 2018, ApJ, 857, 18
  * Shi et al. (2016) Shi J.-M., Zhu Z., Stone J. M., Chiang E., 2016, MNRAS, 459, 982
  * Simon & Armitage (2014) Simon J. B., Armitage P. J., 2014, ApJ, 784, 15
  * Tabone et al. (2022) Tabone B., Rosotti G. P., Cridland A. J., Armitage P. J., Lodato G., 2022, MNRAS, 512, 2290
  * Tobin et al. (2016) Tobin J. J., et al., 2016, Nature, 538, 483
  * Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
  * Veronesi et al. (2020) Veronesi B., et al., 2020, MNRAS, 495, 1913
  * Walmswell et al. (2013) Walmswell J., Clarke C., Cossins P., 2013, MNRAS, 431, 1903
  * Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
  * Xu & Kunz (2021) Xu W., Kunz M. W., 2021, MNRAS, 508, 2142
  * Youdin & Lithwick (2007) Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588
  * Zhang et al. (2018) Zhang S., et al., 2018, ApJ, 869, L47

##  Appendix A Value of $\Lambda$

In this appendix we discuss the value of $\Lambda$ used to compute the Jeans
mass for the two fluid components model. The value of $\Lambda$ can be
obtained through the dispersion relation of the two fluid component model with
drag force (Longarini et al., 2023), by computing the most unstable wavelength
having fixed $(\epsilon,\xi,St)$. Figure 15 shows the value of $\Lambda^{3}$
as a function of $(\epsilon,\xi,St)$. The curves show a sharp fall: this is
the point where instability becomes dust driven. When the instability is gas
driven, the value of $\Lambda^{3}$ decreases with the dust to gas ratio. Then,
it starts increasing since the mass of the dust becomes higher and higher,
significantly contributing to the total density. Conversely, the value of
$\Lambda^{3}$ increases with the relative temperature both in gas and dust
driven regimes. Finally, $\Lambda^{3}$ decreases with the Stokes number.

Figure 15: $\Lambda^{3}$ as a function of $(\epsilon,\xi,St)$, from top to
bottom panel.

##  Appendix B Resolution requirements

In this appendix, we compute the minimum number of dust particles to resolve
the expected Jeans mass: to study this problem, we refer to Lodato & Clarke
(2011). In order to resolve the Jeans length, we require the smoothing length
to be smaller than the disc height: so, the condition for dust particles is

$\frac{h_{d}}{H_{d}}<1.$ (46)

Dust smoothing length and height can be written as a function of gas
properties as follows

$H_{d}=H_{g}\xi^{1/2},$ (47)
$h_{d}=h_{g}\left(\frac{1}{100}\frac{N_{g}\rho_{g}}{N_{d}\rho_{d}}\right)^{1/3},$
(48)

where we assumed $M_{\text{dust}}/M_{\text{gas}}=1/100$. By considering the
gas disc to be marginally unstable ($Q_{g}=1$), condition 46 becomes

$\frac{N_{d}}{2.5\times
10^{5}}>\xi^{-3/2}\left(\frac{\epsilon}{0.01}\right)^{-1}\left(\frac{m(r)}{5\times
10^{-3}}\right)^{-3}\left(\frac{M_{d}/M\star}{0.1}\right),$ (49)

where $\epsilon$ is the dust to gas ratio and
$m(r)=\Sigma_{g}r^{2}/M_{\star}$. If we write the disc mass as $M_{d}=\pi
r_{\text{out}}^{2}\Sigma_{g}$, $m(r)\simeq 1/\pi(M_{d}/M_{\star})$, thus the
resolution requirement is

$N_{d}>4\times
10^{4}\left(\frac{\xi}{0.1}\right)^{-3/2}\left(\frac{\epsilon}{0.01}\right)^{-1}\left(\frac{M_{d}/M_{\star}}{0.1}\right)^{-2}.$
(50)

Figure 16 shows the minimum number of dust particles required to resolve the
Jeans length as a function of the relative temperature $\xi$ for
$M_{d}/M_{\star}=\\{0.05,0.1,0.2\\}$, for a fixed dust to gas ratio
$\epsilon=0.1$: even for extreme values of $\xi$, the standard resolution we
have chosen ($N_{d}=2\times 10^{5}$) is sufficient to resolve the Jeans
length. One should note that we never reach $\xi\sim 10^{-2}$, and in addition
we have underestimated the dust to gas ratio: indeed, inside spiral arms it
could reach values $>0.1$, making our $N_{d}$ choice even safer.

Figure 16: Minimum number of dust particles as a function of the relative
temperature $\xi$ for different values of the disc-to-star mass ratio. The
black line represents the value of this work

##  Appendix C Large and small dust particles comparison

Figure 18 show the complete set of histograms of the dust to gas ratio and the
relative temperature, for small and large dust particles.

Figure 17: Comparison between large and small dust particles. Figure 18:
Comparison between large and small dust particles.