# Direct Formation of Planetary Embryos in Self-Gravitating Disks

Hans Baehr Department of Physics and Astronomy, University of Nevada, Las
Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA  Nevada Center
for Astrophysics, University of Nevada, Las Vegas, 4505 South Maryland
Parkway, Las Vegas, NV 89154, USA  Zhaohuan Zhu (朱照寰) Department of Physics
and Astronomy, University of Nevada, Las Vegas, 4505 South Maryland Parkway,
Las Vegas, NV 89154, USA  Nevada Center for Astrophysics, University of
Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA
Chao-Chin Yang (楊朝欽) Department of Physics and Astronomy, University of
Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154, USA
Nevada Center for Astrophysics, University of Nevada, Las Vegas, 4505 South
Maryland Parkway, Las Vegas, NV 89154, USA

###### Abstract

Giant planets have been discovered at large separations from the central star.
Moreover, a striking number of young circumstellar disks have gas and/or dust
gaps at large orbital separations, potentially driven by embedded planetary
objects. To form massive planets at large orbital separations through core
accretion within disk lifetime, however, an early solid body to seed pebble
and gas accretion is desirable. Young protoplanetary disks are likely self-
gravitating, and these gravitoturbulent disks may efficiently concentrate
solid material at the midplane driven by spiral waves. We run 3D local
hydrodynamical simulations of gravitoturbulent disks with Lagrangian dust
particles to determine whether particle and gas self-gravity can lead to the
formation of dense solid bodies, seeding later planet formation. When self-
gravity between dust particles is included, solids of size $\mathrm{St}=0.1$
to $1$ concentrate within the gravitoturbulent spiral features and collapse
under their own self-gravity into dense clumps up to several $M_{\oplus}$ in
mass at wide orbits. Simulations with dust that drift most efficiently,
$\mathrm{St}=1$, form the most massive clouds of particles, while simulations
with smaller dust particles, $\mathrm{St}=0.1$, have clumps with masses an
order of magnitude lower. When the effect of dust backreaction onto the gas is
included, dust clumps become smaller by a factor of a few but more numerous.
The existence of large solid bodies at an early stage of the disk can
accelerate the planet formation process, particularly at wide orbital
separations, and potentially explain planets distant from the central stars
and young protoplanetary disks with substructures.

protoplanetary disks — planets and satellites: formation — planets and
satellites: gaseous planets — hydrodynamics

††journal: ApJ††software: Matplotlib (Hunter, 2007), SciPy & NumPy (Virtanen
et al., 2020; van der Walt et al., 2011), IPython (Pérez & Granger, 2007)

##  1 Introduction

Observing young stars and disks still embedded in their natal cores is
challenging since the progenitor envelope is opaque at most wavelengths. This
makes it difficult to constrain the initial conditions of planet formation,
including the early sedimentation and concentration of dust. During this early
time, the disk is potentially massive enough to be gravitationally unstable.
Depending on the cooling efficiency, gravitationally unstable disks can either
fragment into large, bound, gas-dominated companions (Boss, 1997; Gammie,
2001; Meru & Bate, 2011b) or become marginally stable to form large scale
spiral structures (Mejía et al., 2005; Zhu et al., 2012; Kratter & Lodato,
2016). While observing the dynamics and processes of gravitationally unstable
disks is difficult due to their rapid accretion and thus short lifespans, they
can play an important role on subsequent planet formation (Forgan, 2019).

Table 1: List of Simulations:

Model | $\beta$ |  | Particle
---
backreaction?
$\mathrm{St}$ | $Q$ | $\alpha_{R}$ | $\alpha_{G}$ | $H_{d}[H_{g}]$ | $\delta_{d,x}$ | $\delta_{d,z}$ | $\mathrm{Sc}$ |  | $\sigma_{d,x}$
---
$[c_{s}]$
| $\sigma_{d,z}$
---
$[c_{s}]$
noBR_S_10 | $10$ | No | $0.1$ | $1.2$ | $4.9\times 10^{-3}$ | $0.03$ | $9.8\times 10^{-2}$ | $4.5\times 10^{-3}$ | $9.6\times 10^{-4}$ | $36$ | $0.36$ | $0.09$
noBR_L_10 | $10$ | No | $1$ | $1.3$ | $6.5\times 10^{-3}$ | $0.044$ | $6.7\times 10^{-3}$ | $2.5\times 10^{-2}$ | $4.5\times 10^{-5}$ | $113$ | $0.56$ | $0.08$
BR_S_10 | $10$ | Yes | $0.1$ | $1.3$ | $7.9\times 10^{-3}$ | $0.029$ | $9.1\times 10^{-2}$ | $5.6\times 10^{-3}$ | $8.4\times 10^{-4}$ | $44$ | $0.43$ | $0.11$
BR_L_10 | $10$ | Yes | $1$ | $1.3$ | $6.3\times 10^{-3}$ | $0.039$ | $9.3\times 10^{-3}$ | $2.1\times 10^{-2}$ | $8.6\times 10^{-5}$ | $406$ | $0.72$ | $0.08$
BR_XL_10 | $10$ | Yes | $10$ | $1.3$ | $6.0\times 10^{-3}$ | $0.021$ | $5.8\times 10^{-2}$ | $1.1\times 10^{-2}$ | $1.2\times 10^{-4}$ | $225$ | $0.87$ | $0.07$
BR_L_5 | $5$ | Yes | $1$ | $1.4$ | $9.7\times 10^{-3}$ | $0.11$ | $1.1\times 10^{-2}$ | $1.9\times 10^{-2}$ | $3.4\times 10^{-2}$ | $3.5$ | $0.94$ | $0.09$
noBR_L_10_HR | $10$ | No | $1$ | $1.4$ | $6.8\times 10^{-3}$ | $0.048$ | $1.1\times 10^{-2}$ | $1.4\times 10^{-2}$ | $1.2\times 10^{-4}$ | $458$ | $0.66$ | $0.09$

. Simulations are initially marginally gravitationally stable such that the
gas $Q_{0}=1.02$ and the solid-to-gas mass ratio $Z_{0}=M_{d}/M_{g}=0.01$. For
measured quantities within the particle clouds, see Table 2.

Note. — Simulation parameters and steady-state values of the stability
parameter $Q$, Reynolds’ stress $\alpha_{R}$, gravitational stress
$\alpha_{G}$, dust scale height $H_{d}$, dimensionless particle diffusion
constants $\delta$ and particle velocity dispersions $\sigma$. We define the
Schmidt number $\mathrm{Sc}\equiv(\alpha_{R}+\alpha_{G})/\delta_{d,z}$. All
simulations have spatial resolution $512^{2}\times 256$ with box lengths
$L_{x}=L_{y}=(80/\pi)H_{g}$ and $L_{z}=(40/\pi)H_{g}$, except for the high
resolution simulation which used $1024^{2}\times 512$ cells

Recent ALMA protoplanetary disk observations revealed that substructure forms
early and often, such that rings are nearly ubiquitous and spirals are
occasionally present as well (Andrews et al., 2018; Long et al., 2018; Clarke
et al., 2018; Segura-Cox et al., 2020; Sheehan et al., 2020). Although the
mechanisms to form rings and gaps are heavily debated, the presence of
undetected planets which can carve open these gaps is a notable option (Zhang
et al., 2018; Lodato et al., 2019; Wang et al., 2021), a scenario being tested
by gas kinematic observations (Pinte et al., 2018; Teague et al., 2021). On
the other hand, the existence of systems with multiple giant planets on wide
orbits, such as HR 8799 (Marois et al., 2010; Maire et al., 2016) suggests
that massive planets can indeed form far away from the central star. Thus, it
is natural to ask how to form planets at large orbital separations where core
and pebble accretion models are less efficient. Planetesimals and even
planetary cores may have formed before the Class II stage and could have begun
accreting gaseous envelopes and opening gaps in disks.

The fragmentation of gravitationally unstable disks is sometimes invoked as a
mechanism to directly form giant planets at distant orbital separations
(Cheetham et al., 2018; Bonnefoy et al., 2018; Morales et al., 2019; Janson et
al., 2021), with initial masses expected to be at least a few Jupiter masses
(Boss, 1997; Baehr et al., 2017). However, the planets that potentially
inhabit the gaps of young disks are typically sub-Jovian, likely too low to be
explained by direct fragmentation of the gas disk (Zhang et al., 2018; Lodato
et al., 2019). Tidal downsizing of companions formed by disk fragmentation
could explain the smaller planets at closer radii, but, at orbital separations
greater than $\sim$10 au, the tidal force cannot strip enough material
(Nayakshin, 2010). Thus, the planets may still grow in the traditional “core-
accretion” fashion, starting with dense solid embryos or cores that only
accrete substantial gas envelopes later. Self-gravitating disks which do not
fragment but are instead marginally unstable (a.k.a. gravitoturbulent) could
potentially concentrate solid material enough so that the dust clouds are
gravitationally bound and form planetary cores/embryos directly (Rice et al.,
2004; Boley & Durisen, 2010).

Marginally gravitationally unstable disks occur when the gravitational
stability parameter (Safronov, 1960; Toomre, 1964; Goldreich & Lynden-Bell,
1965)

$Q=\frac{c_{\mathrm{s}}\Omega}{\pi G\Sigma},$ (1)

is slightly above unity. With $G$ being the gravitational constant, stability
from thermal pressure and rotational shear is quantified through the gas sound
speed $c_{\mathrm{s}}$ and Keplerian orbital frequency $\Omega$, respectively.
These counteract gravitational collapse of perturbations in a disk with gas
surface density $\Sigma$. In addition, if the cooling timescale
$t_{\mathrm{c}}$ is long enough, $t_{\mathrm{c}}>\beta\Omega^{-1}$ (Gammie,
2001), fragmentation of the disk into dense gas structures may be prevented.
There remains some uncertainty about the value of $\beta$ and whether this is
a sufficient criterion (Paardekooper et al., 2011; Meru & Bate, 2012; Brucy &
Hennebelle, 2021). In 3D simulations like those in this work, a criterion of
$\beta=3$ has been shown to converge with resolution (Baehr et al., 2017).
Thus, we focus on gravitoturbulent disks with $\beta>3$.

Spiral arms are generated in these gravitoturbulent disks (Cossins et al.,
2009). These spirals can concentrate solids and potentially produce
planetesimals, as shown in both local (Gibbons et al., 2012, 2014; Shi et al.,
2016) and global simulations (Rice et al., 2004; Boley & Durisen, 2010; Booth
& Clarke, 2016; Cadman et al., 2020). The studies using local simulations were
all two-dimensional which did not consider how dust settling affects the
formation of a thin dust layer. The global studies are all in 3D, but did not
include self-gravity of the dust, which is necessary for the formation of
self-bound dust clouds. Even so, these studies showed that gravitoturbulent
disks could potentially provide young disks with the early planetesimals or
embryos. Among the key questions that need to be addressed with
gravitoturbulent planetesimal formation is: 1) what are the size of bound
clumps of dust and 2) will these bound clumps collapse into solid objects or
be disrupted by high gas or particle velocities?

In this paper, we use 3D local hydrodynamical simulations to model the
interaction of particles in a gravitationally unstable but non-fragmenting
disk. We include particle backreaction and self-gravity for a self-consistent
treatment of particle dynamics. In Sections 2 and 3 we detail the necessary
theory of dust dynamics in marginally gravitationally unstable disks and the
numerics of the Pencil code111http://pencil-code.nordita.org/, respectively.
Section 4 details the identification of clumps and the scaling relations to
determine clump masses. In Section 5 we discuss the results, focusing on
particle and gas velocities and the masses of the resulting particle clumps.
We continue in Section 6 with a discussion of the implications on
protoplanetary disk evolution and planet formation and summarize with our main
conclusions in Section 7.

##  2 Gravitational Collapse

When the density of a region is high enough in a diffuse medium, the high
density region can undergo gravitational collapse. Whereas the collapse of an
adiabatic gas is resisted by the thermal pressure of the gas (Jeans, 1902;
Gammie, 2001; Kratter & Murray-Clay, 2011), this thermal pressure does not
provide dust particles any extra stability against collapse. Instead, the
random motion of the particles determines the resistance to gravitational
collapse. Collapse may be prevented on scales closer to the final solid object
due to the terminal velocity of dust collapsing in a gas (Wahlberg Jansson &
Johansen, 2017; Visser et al., 2021), but that is beyond the resolution
capabilities of our simulations. Here, we only consider that the collapse of
dust in a self-gravitating disk relies on concentrating enough dust locally to
overcome internal turbulent diffusion processes which may arise through
coupling to turbulent motions of the gas (Klahr & Schreiber, 2020, 2021).
Without considering any particle movement, a cloud of particles at orbital
radius $R$ around a star of mass $M_{*}$ should be unstable to collapse when
the cloud has a density higher than the Hill density

$\rho_{\mathrm{Hill}}=\frac{9}{4\pi}\frac{M_{*}}{R^{3}}.$ (2)

This should not be confused with the Roche density of a body (Chandrasekhar,
1963; Shi & Chiang, 2013),

$\rho_{R}=3.5\frac{M_{*}}{R^{3}},$ (3)

which is the minimum density a body needs to remain bound. In this paper, the
former is used to establish the collapse criterion of diffusive particles in a
disk and the latter is used to identify persistent dense objects.

Turbulent gas flows can simultaneously facilitate and disrupt the local
concentration of dust. Turbulent motions can distribute dust between eddies or
to local gas density maxima (Squires & Eaton, 1991; Cuzzi et al., 2001;
Johansen et al., 2006), however at the same time impart diffusive particle
motions, keeping particles from remaining at high densities for long. The
strength of this dust diffusion depends on how well particles couple to the
gas turbulent motions (Youdin & Lithwick, 2007) and is defined in one
direction by the diffusion constant $D$ as the average distance a particle
moves from its original position per unit time (i.e. Johansen & Youdin, 2007;
Yang et al., 2009):

$D_{d,x}\equiv\frac{1}{2}\frac{d\langle|x(t)-x(0)|^{2}\rangle}{dt}\,,$ (4)

where $x(t)$ is the position of a particle at time $t$.

Turbulence in gravitoturbulent disks is largely subsonic, but regions of
supersonic turbulence can be found at various heights above the disk midplane
(Shi & Chiang, 2014), particularly near the midplane of spiral density
features (Cossins et al., 2009; Riols et al., 2020). The motion of larger dust
particles is largely determined by the gravitational interaction with the
gaseous spirals, leading to particle concentration within the spirals. For a
dust particle, the gravitational force from the spiral is stronger than the
aerodynamic drag force when (Shi et al., 2016; Baehr & Zhu, 2021a)

$\frac{Q}{\mathrm{St}}\lesssim 1\,,$ (5)

where $\mathrm{St}$ is the particle’s Stokes number, normally proportional to
the particle’s size. For a marginally stable disk $Q\approx 1$, this condition
is met for $\mathrm{St}\gtrsim 1$. Dust particles settle more efficiently in
gravitoturbulent disks than in disks which do not include self-gravity, due to
the gas self-gravity and the anisotropic turbulence that has weaker vertical
particle diffusion compared to radial transport (Riols et al., 2020; Baehr &
Zhu, 2021b). Including the drift of particles to density maxima, this results
in a dense particle layer predominantly at the spiral midplane.

Smaller particles which are well-coupled to the gas are more affected by the
aerodynamic drag force than the gravitational force from the gas, which leads
to lower velocity dispersions (Booth & Clarke, 2016). The velocity dispersion
$\sigma$ for a collection of particles is defined as

$\sigma=\sqrt{\frac{1}{N_{\mathrm{par}}}\sum_{j}^{N_{\mathrm{par}}}|\bm{w}_{j}-\langle\bm{w}\rangle|^{2}},$
(6)

where $\bm{w}_{j}$ is the velocity of particle $j$ and $\langle\bm{w}\rangle$
is the average velocity of all $N_{\mathrm{par}}$ particles in a clump or in
the whole simulation domain. As with particle diffusion, greater dispersion of
the particle velocities makes it more difficult for particles to
gravitationally collapse.

Figure 1: The vertically-integrated dust particle density $\Sigma_{d}$ at
$t=60\,\Omega^{-1}$ for simulations without dust backreaction (left) and with
dust backreaction (right), normalized by the initial dust particle surface
density $\Sigma_{d,0}=0.01\Sigma_{g,0}$. All particles are size
$\mathrm{St}=1$, which migrate efficiently towards regions of high gas
density, triggering collapse of high density particle clouds. Red circles
indicate locations where the particle density at one or more cells is above
Roche density (3) and gravitationally bound. The size of the circle indicates
the Hill radius (26) of each dust clump.

Including particle backreaction changes the picture slightly. In this case, a
collection of particles can push the gas around as they move through the gas,
but this effect is small when local dust-to-gas mass ratios are less than one.
Only when the dust strongly accumulates does it begin to affect gas
velocities.

When one includes diffusive movement of particles, the condition for particle
cloud collapse can be derived by balancing the gravitational force and the
particle diffusion. A cloud of particles with uniform density
$\rho_{\mathrm{d}}$ and size $r$ can be prevented from collapsing at large
scales by tidal shear and at small scales by the turbulent diffusion generated
by a particle-gas instability, such as the streaming instability (SI; Johansen
& Youdin, 2007; Yang & Zhu, 2021) or the Kelvin-Helmholtz instability (KHI;
Weidenschilling, 1980). A stability parameter $Q_{\mathrm{d}}$ is derived in
Klahr & Schreiber (2021); Gerbig et al. (2020):

$Q_{\mathrm{d}}=\frac{3}{2}\frac{Q}{\epsilon
Z}\sqrt{\frac{\delta_{x}}{\mathrm{St}}}<1.$ (7)

In other words, the gas gravitational stability parameter $Q$ is modified by
dimensionless radial diffusion coefficient $\delta_{x}=D_{d,x}/(c_{s}H_{g})$,
the particle’s Stokes number $\mathrm{St}$, the overall dust-to-gas mass ratio
of the clump $Z$ and gas scale height $H_{g}$. The factor $\epsilon$ is a
measure of the local enhancement of the dust surface density as a result of
radial-azimuthal turbulent concentration of the dust
$\epsilon\equiv\Sigma_{d,max}/\langle\Sigma_{d}\rangle$.

Figure 2: Simulations similar to those in Figure 1 but with $\mathrm{St}=0.1$.
The panel on the left is the simulation without dust backreaction and the
panel on the right includes backreaction. Since dust of this size does not
drift towards the gas structures as efficiently, clumps are fewer and less
massive, but still of considerable mass.

In the simulations of this study, we test the above criterion for scales where
GI turbulence dominates dust diffusion and SI and KHI are not resolved.
Furthermore, no radial pressure gradient is included in the simulations of
this paper, so even with sufficient resolution SI does not develop.

##  3 Model

Figure 3: Top row: On the left, the histogram of particle velocity magnitudes
in simulations with $\mathrm{St}=1$ particles at $t=60\,\Omega^{-1}$. Sharp
peaks in the velocity space correspond to concentrations of particles that
have low relative velocities (within the clump), but may have a high
collective motion. These peaks thus show the velocity of the center of mass of
each particle clump. On the right, the histogram of gas velocities at the same
point in time. Bottom row: particle velocity (left) and gas velocity (right)
histograms of simulations with other particle sizes $\mathrm{St}=0.1$ and
$\mathrm{St}=10$.

For this study we conduct 3D hydrodynamic shearing box simulations of a self-
gravitating disk with Lagrangian super-particles embedded in the Eulerian mesh
using the Pencil code (Brandenburg, 2003; Pencil Code Collaboration et al.,
2021). Both gas and particles are treated as self-gravitating and particle
backreaction is calculated on the gas by mapping the change of particle
momentum due to the gas drag back to the grid with triangular-shaped clouds
(Youdin & Johansen, 2007). Local simulations that allow the Toomre wavelength
$\sim 2\pi H_{g}$ to be well-resolved in the radial and azimuthal coordinates
($x$ and $y$ in the linearized coordinates, respectively) should avoid
spurious fragmentation (Truelove et al., 1997; Nelson, 2006). All simulations
use $512^{2}\times 256$ grid cells with box lengths
$L_{x}=L_{y}=(80/\pi)H_{g}$ and $L_{z}=(40/\pi)H_{g}$ such that $\Delta
x=\Delta y=\Delta z\simeq 0.05\,H_{g}$, where $H_{g}$ is the vertical scale
height of the gas at the initial uniform temperature. At this grid resolution,
we run the simulations with 512 processors up to a simulation time of
$t=80\,\Omega^{-1}$. This takes around 100-200 hours ($\sim$50,000 to 100,000
cpu-hours per simulation) depending on the particle size, where smaller
particle sizes require smaller timesteps to resolve high dust-to-gas ratios.
We run a high resolution simulation with $\Delta x=\Delta y=\Delta z\simeq
0.025\,H_{g}$ up to $t=30\,\Omega^{-1}$ when an initial distribution of clumps
has formed to check for convergence.

Shearing box simulations use hydrodynamic equations which are linearized and
transformed into Cartesian coordinates co-rotating in a Keplerian disk, where
$q=-d\mathrm{ln}\Omega/d\mathrm{ln}R=3/2$ is the shear parameter:

Figure 4: Same as Figure 3, but at $t=72\,\Omega^{-1}$. The top panels show
the simulations with particle size $\mathrm{St}=1$, and have the same general
features as the earlier snapshot. Both simulations with particle backreaction
have more dust and gas at trans- to supersonic velocities, compared to the
simulations without backreaction.

$\displaystyle\frac{\partial{\rho_{\mathrm{g}}}}{\partial t}$
$\displaystyle-q\Omega x\frac{\partial{\rho_{\mathrm{g}}}}{\partial
y}+\nabla\cdot(\rho_{\mathrm{g}}\bm{u})=f_{D}(\rho_{\mathrm{g}})$ (8)
$\displaystyle\frac{\partial\bm{u}}{\partial t}$ $\displaystyle-q\Omega
x\frac{\partial\bm{u}}{\partial y}+\bm{u}\cdot\nabla\bm{u}=-\frac{\nabla
p}{\rho_{\mathrm{g}}}+q\Omega u_{x}\bm{\hat{y}}$
$\displaystyle-2\bm{\Omega}\times\bm{u}-\nabla\Phi-\mathbf{g}-\frac{\varepsilon}{\tau_{s}}(\bm{u}-\bm{w})+f_{\nu}(\bm{u})$
(9) $\displaystyle\frac{\partial s}{\partial t}$ $\displaystyle-q\Omega
x\frac{\partial s}{\partial
y}+(\bm{u}\cdot\nabla)s=-\frac{\Lambda}{\rho_{\mathrm{g}}T}+f_{\chi}(s).$ (10)

In equations (8) - (10), $\mathbf{u}$ is the gas velocity deviation from the
background shear velocity in the local box, $\bm{w}$ is the particle velocity
which imparts a backreaction onto the gas proportional to the local dust-to-
gas ratio $\varepsilon\equiv\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$,
$\rho_{\mathrm{g}}$ is the gas density, $\rho_{\mathrm{d}}$ is the dust
density. The thermodynamic variable is the specific entropy $s$, while $p$ is
the gas pressure and $T$ is the gas temperature. The vertical gravitational
acceleration $\bm{g}=g\hat{\bm{z}}$ due to the central potential is a linear
profile modified with zero acceleration near the $z$-boundary to avoid an
abrupt discontinuity at the periodic vertical boundary.

We use an ideal equation of state such that

$p=(\gamma-1)\rho_{g}e,$ (11)

where $\gamma\equiv c_{p}/c_{v}\equiv 5/3$ is the adiabatic index, $c_{p}$ and
$c_{v}$ are the specific heat capacities at constant pressure and volume,
respectively, and $e=c_{v}T$ is the specific internal energy. Heat can be
generated through the dissipation of shocks via a dimensionless shock
viscosity of $\nu_{sh}=5.0$ (see Lyra et al., 2008b) and through compression
of the gas. The relationship between the gas temperature $T$, the sound speed
$c_{s}$, and entropy is given by

$c_{s}^{2}=(\gamma-1)c_{p}T=c_{s,0}^{2}\exp[\gamma
s/c_{p}+(\gamma-1)\ln(\rho/\rho_{0})]$ (12)

where $c_{s,0}$ is the initial uniform speed of sound. Heat is lost via the
simple $\beta$-cooling prescription

$\Lambda=\frac{\rho_{g}(c_{\textnormal{s}}^{2}-c_{\textnormal{s,irr}}^{2})}{(\gamma-1)t_{\textnormal{c}}}$
(13)

with the cooling timescale $t_{\mathrm{c}}$ parametrized as
$t_{\mathrm{c}}=\beta\Omega^{-1}$. We include a background irradiation term
which is different from most previous simulations. This background irradiation
is important at the outer regions of protoplanetary disks where GI prevails.
The background temperature mimics the effect of stellar irradiation (D’Alessio
et al., 1998), keeping the gas from dropping below the initial $Q_{0}$ due to
a low local gas temperature. We choose $c_{\textnormal{s,irr}}=c_{s,0}$ such
that the disk cools towards the initial uniform temperature. This cooling
prescription has no dependence on local variations in optical depth and thus
all regions cool with the same timescale. With a more realistic treatment of
thermodynamics, the opacity is dominated by small dust grains and an increase
of the particle density can increase the local cooling timescale.

The term $\epsilon(\bm{u}-\bm{w})/\tau_{s}$ in Equation (9) is the
backreaction of particles with stopping time $\tau_{s}$ (see Equation (19)
below), included or removed depending on whether this effect is active or
inactive. Hyperdissipation is applied with the terms
$f_{D}(\rho_{\mathrm{g}})$, $f_{\nu}(\bm{u})$, $f_{\chi}(s)$ which for each
has the form

$f(\xi)=\nu_{3}(\nabla^{6}\xi),$ (14)

with hyperdissipation constant $\nu_{3}=2.5\,H_{g}^{6}\Omega$, which leads to
a mesh Reynolds number of about 0.15 (Yang & Krumholz, 2012; Lyra et al.,
2017) during fully developed gravitoturbulent state.

The gravitational potential of the gas and dust is solved in Fourier space by
transforming the density to find the potential at wavenumber $k$ and
transforming the solution back into real space. The solution is shifted in the
$y$-direction by applying a phase in the Fourier space such that the shear
periodic boundary conditions are accounted for (i.e. Johansen et al., 2007).
The solution to the Poisson equation in Fourier space at sheared wavenumber
$k$ is

$\Phi(\bm{k},t)=-\frac{4\pi G\rho(\bm{k},t)}{\bm{k}^{2}},$ (15)

where the density of the particles and gas is combined
$\rho=\rho_{\mathrm{g}}+\rho_{\mathrm{d}}$ to produce a potential
$\Phi=\Phi_{\mathrm{g}}+\Phi_{\mathrm{d}}$. This means that self-gravity
includes both gas and particles and affects both components.

Figure 5: A pair of $x$-$z$ slices at $t=60\Omega^{-1}$ through simulations
with $\mathrm{St}=1$ dust particles (left) and $\mathrm{St}=0.1$ dust
particles (right), both without dust backreaction. Vectors indicate the gas
flow direction and magnitude in the $x-z$ plane, while contours in yellow and
purple indicate regions where the 3D gas velocity magnitude is above Mach
number equals 1 and 2, respectively. Blue colors show the gas density
structure and gray particles plotted as triangular-shaped clouds, which are
settled close to $z=0$.

Every particle in our simulations, also known as a superparticle, represents a
collection of solids such that the $i$-th superparticle has position
$\bm{x}^{(i)}$ and velocity $\bm{w}^{(i)}$ as in Youdin & Johansen (2007);
Yang & Johansen (2016). It follows that

$\displaystyle\frac{d\bm{w}^{(i)}}{dt}$
$\displaystyle=-2\Omega\times\bm{w}^{(i)}+q\Omega
w^{(i)}_{x}\bm{\hat{y}}-\nabla\Phi$ (16)
$\displaystyle-\frac{1}{\tau_{s}}\left(\bm{w}^{(i)}-\bm{u}(\bm{x}^{(i)})\right),$
(17) $\displaystyle\frac{d\bm{x}^{(i)}}{dt}$
$\displaystyle=\bm{w}^{(i)}-q\Omega x^{(i)}\bm{\hat{y}}.$ (18)

In the Epstein regime, the stopping time is proportional to the particle size
as (Weidenschilling, 1977)

$\tau_{\mathrm{s}}=\frac{a\rho_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}}{c_{\mathrm{s}}\rho_{\mathrm{g}}}$
(19)

where $a$ is the particle diameter and
$\rho_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}}$
is the material density of an individual dust particle. The dimensionless
friction time $\tau_{f}$ (also referred to as the Stokes number $\mathrm{St}$)
is the particle stopping time normalized by the orbital time $\Omega^{-1}$

$\mathrm{St}=\tau_{f}=\tau_{\mathrm{s}}\Omega.$ (20)

Larger particles have higher Stokes numbers, are less coupled to small scale
gas motions, and retain their initial motion for longer. As a corollary,
smaller particles with low Stokes numbers are well-coupled and closely move
with the gas. Particles are added such that the initial distribution maintains
a physically motivated metallicity of $Z=0.01$, roughly that of the
interstellar medium (ISM). The vertical dust distribution follows a vertical
Gaussian profile with the same width as the gas. Particle mass is calculated
from the total gas mass and the specified metallicity (see also Section 4).

Since the radial pressure gradient is set as zero, particles do not drift in
the disk. The contribution of particles to the gravitational potential is
initially minuscule and the potential is dominated by the gas distribution,
but with rapid dust settling to the midplane, the entire dust layer is covered
by 5 to 10 grid cells ($\sim 0.25-0.5\,H_{g}$) (Baehr & Zhu, 2021b). Since our
simulations include $1.5\times 10^{6}$ particles we resolve the midplane layer
($\sim 5$ cells thick) with roughly $1.1$ particles per cell.

##  4 Analysis

Our simulation parameters are summarized in Table 1. Simulations without dust
backreaction are indicated by the ‘noBR’ segment in the name and ‘BR’ if dust
backreaction is included. Three sizes of particles are used, denoted by ‘S’
($\mathrm{St=0.1}$), ‘L’ ($\mathrm{St=1}$) and ‘XL’ ($\mathrm{St=10}$). A
final simulation at a higher resolution is indicated by ‘HR’. Included in
Table 1 are time and space averaged diagnostics of the gas gravitational
stability $Q$, Reynolds and gravitational stresses $\alpha_{R}$ and
$\alpha_{G}$, defined as

$\alpha_{R}=\frac{2}{3}\frac{\langle\rho u_{x}u_{y}\rangle}{\langle\rho
c_{\mathrm{s}}^{2}\rangle},$ (21)

and

$\alpha_{G}=\frac{2}{3}\frac{\langle g_{x}g_{y}\rangle}{4\pi G\langle\rho
c_{\mathrm{s}}^{2}\rangle},$ (22)

where $g_{x}$ and $g_{y}$ are the gravitational accelerations in the radial
and azimuthal directions, respectively. As in Baehr & Zhu (2021b), these
stresses are calculated as simple volume averages over the time range
$t=50\,\Omega^{-1}$ to $t=80\,\Omega^{-1}$ when the gas gravitoturbulence has
been established. The sum of these stresses yields the total $\alpha$ stress
in the disk (Shakura & Sunyaev, 1973; Gammie, 2001).

Also measured in Table 1 are the dimensionless particle diffusion constants
$\delta\equiv D_{d}/c_{s}H_{g}$ and velocity dispersions $\sigma$ in the
vertical and radial directions. Particle diffusions are calculated in Pencil
using the scheme from Yang et al. (2009) where the particle displacement
$x(t)-x(0)$ of all particles is averaged as in Equation (4). The particle
velocity dispersion is calculated using all particles according to Equation
(6) and averaged in time between $t=50\,\Omega^{-1}$ and $t=80\,\Omega^{-1}$.
While diffusion may sometimes be assumed equal in each direction, turbulent
processes in protoplanetary disks can be non-isotropic (Zhu et al., 2015) and
may result in particle concentration in certain structures (Yang et al.,
2018).

We identify dense clumps based on the cells in which the dust density is at or
above Roche surface density

$\Sigma_{R}\approx\sqrt{2\pi}H_{g}\rho_{R}=8.8\frac{M_{*}H_{g}}{R^{3}}=8.8\frac{\Omega^{2}H_{g}}{G}.$
(23)

Since the dust scale height $H_{d}$ is much smaller than the gas scale height
$H_{g}$, this serves as a criterion to identify the largest dust clouds,
although smaller bound clouds may exist. The Roche surface density may not be
a definitive criterion since a dense particle layer may be susceptible to
disruption from shockwaves that can transverse the particle concentration
faster than it can collapse (Shi & Chiang, 2013). Nevertheless, we find the
Roche density to be a suitable threshold for the identification of persistent
clumps.

We search for clumps with a $3\times 3$ grid cell footprint over the
vertically integrated particle density. The particles are so settled to the
midplane that the volumetric density in the midplane and surface density are
closely correlated. Within each $3\times 3$ patch, cells which are at or above
the Roche surface density are used to calculate the total dust mass of each
clump $M_{c}$. From the clump mass, we compute its Hill radius $R_{H}$, the
region around the clump where its gravity dominates the tidal force. Most
clumps are small enough that the Hill radius is within the $3\times 3$
footprint, but for larger clumps, the mass estimate might neglect some
particles that are within the Hill radius, but outside the $3\times 3$ grid
cells, where particle density is usually low compared to the central density
(Figures 1 and 2).

$R_{H}=\left(\frac{M_{c}}{3M_{*}}R^{3}\right)^{1/3}$ (24)

We can rewrite $M_{*}/R^{3}$ from our adopted Toomre stability parameter
$Q_{0}=1.02$ and initial midplane density $\rho_{m,0}$ as

$\frac{M_{*}}{R^{3}}=\frac{\Omega^{2}}{G}=\sqrt{2\pi}\pi\rho_{m,0}Q_{0}$ (25)

which yields an expression for the Hill radius of a clump

$R_{H}=\left(\frac{M_{c}}{3\sqrt{2\pi}\pi\rho_{m,0}Q_{0}}\right)^{1/3}\approx\left(\frac{M_{c}}{24\rho_{m,0}}.\right)^{1/3}$
(26)

We plot the Hill radius of each identified bound clump in Figures 1 and 2
which range from two to eight grid cells in radius.

Figure 6: A 3D render of the snapshot at the left in Figure 1. Dense particle
clumps show as bright red regions while dense particle streams that have not
collapsed are orange to yellow. The most massive clouds $>1M_{\oplus}$ are
circled. Values for gas and particle densities are normalized by the initial
mid-plane gas density.

Similar to Schäfer et al. (2017), we provide the following scaling relations,
using a mean molecular weight of $\mu=2.33$ at a reference distance of
$50\,au$ around a $1\,M_{\odot}$ star, to convert our results in physical
units for the sound speed $c_{s}$, orbital frequency $\Omega$, orbital period
$P$, and gas scale height $H_{g}$.

$c_{s}=199\left(\frac{T}{11.25\,K}\right)^{1/2}\,m\,s^{-1}$ (27)
$\Omega=1.78\times
10^{-2}\left(\frac{R}{50\,au}\right)^{-3/2}\left(\frac{M_{*}}{1M_{\odot}}\right)^{1/2}\,yr^{-1}$
(28)
$P=353\left(\frac{R}{50\,au}\right)^{3/2}\left(\frac{M_{*}}{1M_{\odot}}\right)^{-1/2}\,yr$
(29)
$H_{g}=2.36\left(\frac{T}{11.25\,K}\right)^{1/2}\left(\frac{R}{50\,au}\right)^{3/2}\left(\frac{M_{*}}{1M_{\odot}}\right)^{-1/2}\,au.$
(30)

The scale height changes as the local temperature fluctuates, so the $H_{g}$
we define here and use throughout the paper is based on the initial uniform
temperature, i.e., with the speed of sound at
$c_{s}=c_{s,\mathrm{irr}}=c_{s,0}$. The above definitions mean the initial
surface density is

$\begin{split}\Sigma_{0}&=\frac{c_{s,0}\Omega}{\pi GQ_{0}}\\\
&=53\left(\frac{T}{11.25\,K}\right)^{1/2}\left(\frac{R}{50\,au}\right)^{-3/2}\left(\frac{M_{*}}{1M_{\odot}}\right)^{1/2}\,g\,cm^{-2}.\end{split}$
(31)

Figure 7: Normalized cumulative mass distribution of particle concentrations
in the five simulations where particles reach Roche density and are
gravitationally bound. Plotted distributions are at a time when there is
little variation in subsequent snapshots and is considered stable. No particle
clumps were identified in the simulation using particles of size
$\mathrm{St}=10$. Mass is displayed in code units (bottom axis), physical
units assuming the shearing box is at 50 au (top axis; see Equation (33)) and
fit with an exponentially-tapered power law (dotted curves).

In our models, we adopt a length of $H_{g}=\pi$, a unit time of
$\Omega^{-1}=P/2\pi$, and $G=1$, so our unit volume density is

$\begin{split}\hat{\rho}_{0}&=\frac{\Omega^{2}}{G}\\\ &=4.7\times
10^{-12}\left(\frac{R}{50\,au}\right)^{-3}\left(\frac{M_{*}}{1M_{\odot}}\right)\,g\,cm^{-3},\end{split}$
(32)

and our unit mass is

$\begin{split}\hat{M}_{0}&=\frac{4}{\pi}\frac{H_{g}^{3}}{GP^{2}}\\\
&=6.8\times
10^{27}\left(\frac{T}{11.25\,K}\right)^{3/2}\left(\frac{R}{50\,au}\right)^{3/2}\\\
&\times\left(\frac{M_{*}}{1M_{\odot}}\right)^{-1/2}\,g=1.1\,M_{\oplus}\end{split}$
(33)

With this value for code unit of mass, we derive the total gas mass in the
entire simulation box (assuming at values from above at 50 au)

$\begin{split}M_{\mathrm{total,gas}}&=\Sigma_{0}L_{x}L_{y}=80^{2}\hat{M}_{0}/Q_{0}\\\
&=4.3\times 10^{31}\,g=0.021\,M_{\odot}.\end{split}$ (34)

With this total gas mass, and the dust mass assumed is to be at a 1:100 ratio,
the total dust mass is

$\begin{split}M_{\mathrm{total,dust}}&=64\hat{M}_{0}/Q_{0}\\\ &=4.3\times
10^{29}\,g=72\,M_{\oplus}.\end{split}$ (35)

Dividing it by the number of super particles ($1.5\times 10^{6}$) yields the
mass per super particle

$\begin{split}M_{\mathrm{sp}}=2.8\times 10^{23}\,g=4.8\times
10^{-5}\,M_{\oplus}.\end{split}$ (36)

Using the above scaling relations, we can translate the code units to physical
units at any disk radius of interest. The surface density we adopt results in
a disk to star mass ratio around 0.1, which is appropriate for self-
gravitating disks and within the constraints of observations, e.g. GM Aurigae
(Schwarz et al., 2021). Additionally, many protoplanetary disks could be
optically thick in submillimeter observations, resulting in a larger disk mass
than estimates based on the optically thin assumption (Zhu et al., 2019). In
the next section, we proceed with the simulation results, starting with the
dust and gas velocities, and how they depend on our choice of particle size
and dust backreaction, followed by the mass of the dust concentrations that
form.

Figure 8: Evolution of the clump mass distribution in select simulations. The
progression from lighter to darker colors indicates the progression of time in
increments of $4\Omega^{-1}$ until the snapshot which represents the converged
distribution which is indicated with a black outline. Colors in each panel
correspond the simulations of the same color in Figure 7.

##  5 Results

Our simulations show that particles clump into dense, bound clouds up to a few
Earth masses if placed at $\sim$50 au. This is exemplified in Figure 1, which
shows the locations where particle densities within a single grid cell are
above Roche surface density222A movie of the evolution of these clumps can be
found at https://youtu.be/hqtwMIe1dwk. The red circles drawn around them
indicate the Hill radius of each clump of particles, representing the
approximate region around the clump where the gravity of the clump is stronger
than the stellar gravity. Figure 2 shows similar clumping behavior for dust
with a smaller Stokes number (St=0.1), although clumps are smaller and fewer.

###  5.1 Gas and Particle Velocities

In the top two panels of Figure 3, we compare the dust velocities in
simulations with $\mathrm{St}=1$ particles and the gas velocities of all cells
in the simulation domain. Particles of this size drift towards pressure maxima
with the greatest efficiency. The numerous spikes in the velocity histogram
are where a significant number of particles have the same velocity. Thus, they
roughly represent the center-of-mass velocity of each clump. For the case
without particle backreaction, numerous large clumps form and most remain
together at subsonic speeds.

Dust backreaction, combined with the high local enhancements in dust-to-gas
ratio brought on by the clumping pushes the gas velocity in a number of cells
near the mid-plane further into the supersonic regime. Increasing the strength
of the GI turbulence by decreasing the cooling time $\beta$ has a similar
effect and shifts even more gas towards supersonic speeds. Both effects
increase the number of particles in the supersonic regime, but also accelerate
some clumps. Overall, most particle clumps are still within two times the
sound speed. Figure 4 shows the same simulations at a later time $\Delta
t=12\,\Omega^{-1}$. The simulations in the top panel of both Figures 3 and 4
with $\mathrm{St}=1$ have similar velocity patterns in gas and particle
velocities.

Figure 9: Scatter plot of the total gas mass versus total dust mass within the
Hill radius of each identified clump of particles. Dashed lines indicate
constant dust-to-gas mass ratios. Plotted on top and to the right are the
logarithmically-binned histograms of the data.

At the smaller particle size, $\mathrm{St}=0.1$, particles couple more tightly
to the gas and the flows at smaller scales. As shown in Figure 5, less
efficient settling to the midplane by the $\mathrm{St}=0.1$ particles means
gas flows near the surface can occasionally sweep up particles (Riols et al.,
2020; Baehr & Zhu, 2021b). While clumps form with smaller dust, the particle
layer is not as well settled and the local dust loading is lower. We present
the full 3D view of the gas and dust in Figure 6, where we highlight the most
massive clumps with gray circles.

The velocity distribution of the largest dust species ($\mathrm{St}=10$) is
noticeably shifted towards higher velocities than the smaller particle sizes
studied here. Even though there are only small differences between the gas
velocity distributions in all simulations, particles of the largest dust
species is more likely to be at higher velocities and even a significant
amount in the supersonic regime. Larger dust species are less coupled to gas
motions, and gravitational stirring from the dense gas filaments can excite
the larger dust grains to much higher velocities (Shi et al., 2016).

We measure radial particle diffusion constants using the method described in
Section 2. Our models have the same resolution in the radial and azimuthal
directions as in the 2D simulations of Shi et al. (2016) and we derive a
similar increasing trend in the diffusion constants up to $\mathrm{St}=1$.
However, whereas 2D simulations showed that radial diffusion is highest when
$\mathrm{St}=10$ and decreases with increasing particle size, we notice that
radial diffusion is already beginning to decrease from $\mathrm{St}=1$ to
$\mathrm{St}=10$ in our 3D models.

Similar to Baehr & Zhu (2021b), vertical dust diffusion is inferred from the
dust scale height once the simulation has reached an equilibrium and the dust
scale height is quasi-steady in time. With the ratio of dust scale height and
the equilibrium gas scale height due to the gas self-gravity $H=0.75H_{g}$,
the vertical dust diffusion constant can be determined via (Dubrulle et al.
1995; see also Yang et al. 2018)

$\frac{H_{d}}{H}=\sqrt{\frac{\delta_{d,z}}{\delta_{d,z}+\mathrm{St}}},$ (37)

where $\delta_{d,z}\equiv D_{d,z}/c_{s}H_{g}$ (Section 4). The measured
vertical particle diffusion constants are sometimes around two orders of
magnitude lower that the radial particle diffusion, indicating a high degree
of anisotropy in dust transport driven by gravitoturbulence.

Particle diffusion is measured with a cadence of $\Delta t=2\Omega^{-1}$ and
reported in Table 1, along with the corresponding Schmidt number
$\mathrm{Sc}\equiv(\alpha_{R}+\alpha_{G})/\delta_{d,z}$. The Schmidt number is
a dimensionless number that measures how well coupled the dust diffusion is to
the gas turbulence (Dullemond & Dominik, 2004). The Schmidt numbers derived
here are particularly high, as is often the case for self-gravitating
turbulence (Riols et al., 2020; Baehr & Zhu, 2021b). Schmidt numbers are
notably higher in simulations where more gravitationally bound clumps are
identified. This could in part be due to the large number of particles in
clumps will skew the vertical particle distribution towards a narrower profile
and lower measured diffusion constant. The simulation with the largest dust
species has the lowest Schmidt number. However, low vertical diffusion could
instead help cause the collapse of dust rather than be a symptom of particle
collapse.

###  5.2 Cloud Masses

All simulations carried out here show a similar size distribution that is
dominated by a larger number of smaller clumps with a small number of clumps
that are larger by over an order of magnitude in mass. How each population
arrives at this distribution differs both by size of the dust and whether
backreaction is included.

We characterize the mass distribution of solids using an exponentially-tapered
power law, modified to use minimum mass $M_{\mathrm{min}}$ as a fitting
parameter rather than the characteristic mass $M_{\mathrm{pow}}$ as in Schäfer
et al. (2017). We therefore define the number of clumps above a particular
mass $N_{>}(M)$ as

$\frac{N_{>}(M)}{N_{\mathrm{tot}}}=\left(\frac{M}{M_{\mathrm{min}}}\right)^{-\alpha}\mathrm{exp}\left[\left(\frac{M_{\mathrm{min}}}{M_{\mathrm{exp}}}\right)^{\beta}-\left(\frac{M}{M_{\mathrm{exp}}}\right)^{\beta}\right],$
(38)

and use it to fit the distributions at a time $t=60\,\Omega^{-1}$ when gas
stability parameter $Q$ is steady in time. This can be compared with other
studies of planetesimal formation via gravitational collapse of particles
(e.g. Johansen et al., 2015; Schäfer et al., 2017; Abod et al., 2019; Li et
al., 2019). We find that an exponentially tapered power law fits the data
best, but not as well at higher masses, similar to the results of Johansen et
al. (2015), Schäfer et al. (2017) and Abod et al. (2019). Li et al. (2019)
finds the best fit distribution varies with simulation parameters, but a
truncated broken power law, broken cumulative power law or three segment power
law can all match the data.

We find a noticeable difference between clump formation with particle
backreaction versus those without when particle size $\mathrm{St}=1$.
Backreaction from particles causes more rapid concentration of the dust at the
earliest onset of dust sedimentation and collapse. This creates a large number
of smaller dust clumps that will only grow larger through hierarchical mergers
of clumps. Without backreaction however, particles do not segregate into
individual clumps as rapidly, instead forming long filament-like structures
which will later collapse into fewer, but larger objects. This is shown in
Figure 8, where the simulations with dust backreaction have considerably more
identified clumps at the lower mass end to the left, but fewer at the higher
end of the distribution to the right.

Table 2: Particle clump properties:

model | $\overline{N}$ |  | $\overline{\langle M_{\mathrm{dust}}\rangle}$
---
$[\hat{M}_{0}]$
| $M_{\mathrm{dust,max}}$
---
$[\hat{M}_{0}]$
| $\overline{\langle M_{\mathrm{gas}}\rangle}$
---
$[\hat{M}_{0}]$
$\overline{\langle Z\rangle}$ |  | $\overline{\langle\sigma\rangle}$
---
$[c_{s}]$
$\alpha$ | $\beta$ |  | $M_{\mathrm{min}}$
---
$[\hat{M}_{0}]$
| $M_{\mathrm{exp}}$
---
$[\hat{M}_{0}]$
noBR_S_10 | $37$ | $0.11$ | $0.69$ | $0.12$ | $0.70$ | $3.2\times 10^{-2}$ | $-5.9\times 10^{-1}$ | $0.35$ | $4\times 10^{-3}$ | $9.4\times 10^{-4}$
noBR_L_10 | $37$ | $1.01$ | $7.8$ | $0.15$ | $4.8$ | $5.0\times 10^{-2}$ | $3.0\times 10^{-2}$ | $0.58$ | $5\times 10^{-3}$ | $0.74$
BR_S_10 | $16$ | $0.11$ | $0.20$ | $0.14$ | $0.52$ | $1.5\times 10^{-1}$ | $6.0\times 10^{-1}$ | $2.5$ | $1\times 10^{-2}$ | $0.2$
BR_L_10 | $115$ | $0.29$ | $1.9$ | $1.8$ | $10$ | $5.5\times 10^{-2}$ | $-1.1\times 10^{-1}$ | $0.51$ | $4\times 10^{-3}$ | $8.8\times 10^{-4}$
BR_XL_10 | $0$ | – | – | – | – | – | – | – | – | –
BR_L_5 | $164$ | $0.18$ | $3.1$ | $0.16$ | $1.3$ | $8.0\times 10^{-2}$ | $4.6\times 10^{-2}$ | $0.97$ | $5\times 10^{-3}$ | $0.21$
noBR_L_10_HR | $82$ | $0.31$ | $1.7$ | $0.11$ | $3.2$ | $1.0\times 10^{1}$ | $-4.5\times 10^{-1}$ | $0.36$ | $2\times 10^{-3}$ | $0.55$

Note. — Summary of clump properties and particle properties within bound
clumps in each simulation where $N$ is the number of identified bound objects,
$M_{\mathrm{dust}}$ and $M_{\mathrm{gas}}$ are the dust and gas mass within
the Hill sphere of each clump, $Z=M_{\mathrm{dust}}/M_{\mathrm{gas}}$, $n$ is
the number of particles within a clump, $\sigma$ is the particle velocity
dispersion, and $\alpha$, $\beta$, $M_{\mathrm{min}}$ and $M_{\mathrm{exp}}$
are best fit parameters in Equation (38). Quantities in $\langle\cdot\rangle$
brackets are averaged over all clumps in the simulation. Quantities with a bar
above are averaged over multiple snapshots in time from $t=60\,\Omega^{-1}$ to
$t=80\,\Omega^{-1}$, except for the high resolution run which was calculated
from $t=20\,\Omega^{-1}$ to $t=30\,\Omega^{-1}$. The simulation with
$\mathrm{St}=10$ particles ’xlarge’ did not form any bound clumps, but is
included for completeness. Masses are in units of
$\hat{M}_{0}=(4/\pi)G^{-1}H^{3}_{g}p^{-2}$ (Equation (33)).

The more efficiently drifting particles with $\mathrm{St}=1$ are especially
efficient at clumping and form more clumps with higher masses. These clumps
gradually merge over time, and the total number of clumps decreases while the
number in the high mass end increases. This could suggest that planetesimal-
planetesimal accretion is the more common way to grow than pebble accretion.
However, the opposite appears to be the case when dust is dominated by smaller
species. The size distribution of clumps with smaller particles
($\mathrm{St}=0.1$) evolves in a different way as illustrated in Figure 8.
Bound clouds of dust formed when $\mathrm{St}=0.1$ are initially scarcer and
smaller than those at $\mathrm{St=1}$. The number of bound clumps does
steadily rise across the entire mass range however. Since the number of
objects is lower than in the case of $\mathrm{St}=1$ and more space between
them, planetesimal-planetesimal mergers is less likely to be the cause of
clump growth. Instead, the plentiful dust in and between the filamentary dust
structures is likely accreted by the clumps.

Since the simulations presented here are very nearly 3D analogs of those in
Gibbons et al. (2014), we compare our findings with theirs. Their most massive
dust clouds have notably lower masses, nearly two orders of magnitude of lower
than identified here. Additionally, while the number of clumps is not
identified in their work, by visual inspection there appear to be fewer in the
2D simulations.

In Figure 9, we compare clumps the dust and gas masses at a single snapshot in
all simulations. Dust masses vary by a few orders of magnitude, and tend to be
broadly distributed when numbers are low. When large numbers of dust clumps
form, as in the two simulations where dust particles have size $\mathrm{St}=1$
and include backreaction, clumps tend to have masses around $0.1\hat{M}_{0}$.
Gas masses consistently fall within a more narrow range, which makes sense
when one considers that all were formed within the non-axisymmetric gas
structures which do not vary strongly in peak gas density. The small handful
of low gas masses are the result of clumps that have drifted out of the dense
gas structure in which they formed. Overall, most dense particle regions have
high dust-to-gas ratios between $Z=1$ and $Z=100$.

It is important to consider that resolution and initial conditions may have an
effect on the results. Convergence is not a given, particularly in simulations
that concern the gravitational collapse (Meru & Bate, 2011a; Paardekooper,
2012), potentially due to effects of dimensionality (Baehr et al., 2017) or
other numerical factors (Deng et al., 2017; Klee et al., 2017, 2019). Thus we
performed a simulation with double the grid resolution of the noBR_L_10
simulation and compared the mass distribution. The number of particles was
also increased in the high resolution run to maintain the same number of
particles in a midplane layer assumed to be one cell thick, leading to a four
fold increase in particles. The high concentrations of particles that lead to
clumping could also be a result of the linear collapse phase of the gas.
Therefore, we also conducted a simulation where particles were only added to
the simulation after the gas had reached a turbulent steady state.

Figure 10: Comparison of the mass distribution of clumps in simulations with
the same parameters but with different resolution and initial conditions

. All mass distributions are from the same point in time after particle are
introduced for $t=26\,\Omega^{-1}$.

Figure 10 shows the distributions of all three simulations at the same time
$t=26\,\Omega^{-1}$. In each case, the upper clump mass limit is the same and
the number of identified bound clumps is similar, with more smaller clumps
forming in the simulation where dust particles were added later. Clump and
turbulent diagnostics of the high resolution simulation are reported in Tables
1 and 2 and are within reasonable values compared to the standard resolution.
Simulations of planetesimal formation via the streaming instability also
appear similarly converged with resolution, such that the differential size
distribution is unaffected by grid resolution (Johansen et al., 2015).

It should be noted that this suggests convergence for the simulations where
the dust scale height is resolved ($\mathrm{St}=0.1$ and $\mathrm{St}=1$). For
$\mathrm{St}=10$, increasing the resolution should resolve the critical
wavelength of the dust and change the results at this size (see Section 6.1).

##  6 Discussion

###  6.1 Particle Collapse Criterion

Figure 11: Particle stability parameter $Q_{d}$ for three particle sizes which
include the effect of particle backreaction. The diffusion constant and range
of maximum concentration $\epsilon Z$ in Equation (7) is determined from the
$\Sigma_{d,\mathrm{max}}/\langle\Sigma_{d}\rangle$ across a few snapshots
after concentration in the spiral structures has occurred but before particle
self-gravity dominates the high end of particle densities
($t=12$–$16\,\Omega^{-1}$).

The particle collapse criterion of equation (7) from Section 2 can be tested
for the three particle species $\mathrm{St}=0.1,1,10$. Because particle self-
gravity significantly affects the particle density enhancement ($\epsilon Z$)
once gravitational collapse is triggered, the criterion has to be measured
after concentration into filaments has started, but before particles have
collapsed. For this purpose, we consider the time span
$t\in[12\Omega^{-1},16\Omega^{-1}]$ to measure $\epsilon Z$ as well as the
radial diffusion constants $\delta_{x}$ to compute $Q_{\mathrm{d}}$ for the
simulations that include particle backreaction. The resulting range of
$Q_{\mathrm{d}}$ in this duration are plotted in Figure 11.

The two smaller sizes of dust species reach values of $Q_{\mathrm{d}}$ well
below the stability threshold. Considering the simulations have an equilibrium
gas stability parameter of $Q\sim 1.2-1.4$, suggesting that dust particles can
become gravitationally unstable and collapse even when the gas disk is
marginally stable. If disk becomes more unstable with time through the
accretion of envelope material (Küffmeier et al., 2018; Kuznetsova et al.,
2022), structure in massive/gravitoturbulent disks may be determined not only
by the instability of the gas, but the dense solid concentrations that form
first. On the contrary, due to higher radial dust diffusion, the simulation
with the largest dust particles straddles the stability threshold and the
particles in the simulation do not collapse into dense clumps.

When the simulation reaches a quasi-steady state after gravitational collapse
of dust, the diffusion decreases. The dust diffusion can decrease to the
extent that the critical length scale of the dust is no longer resolved and
gravitational collapse at this size cannot occur. The critically unstable
length (Gerbig et al., 2020; Klahr & Schreiber, 2020)

$r_{\mathrm{crit}}=\frac{1}{3}\sqrt{\frac{\delta_{x}}{\mathrm{St}}}H_{g}$ (39)

for $\mathrm{St}=10$ is $r_{\mathrm{crit},10}=0.01H_{g}$ in the turbulent
steady state, which is smaller than the grid spacing, indicating that there is
not enough resolution to capture collapse at this size. Even at
$\mathrm{St}=1$, $r_{\mathrm{crit},1}=0.05H_{g}$ and hence this dust species
is not well resolved, but the dense clumps have already formed by this point
from the initial collapse. Thus to study the collapse of dust sizes
$\mathrm{St}=1$ and larger, grid spacing will have to be nearly an order of
magnitude smaller. At smaller dust sizes, on the other hand, the critical
length is resolved by a few grid cells, albeit not all unstable lengths are
resolved.

Additional diffusive processes, such as the ones driven by the streaming
instability and Kelvin-Helmholtz instability, can prevent collapse at scales
smaller than a scale height, and the particle collapse criterion can be
applied in the same way (Gerbig et al., 2020). While scales smaller than $\sim
0.1H_{g}$ are not resolved in our simulations to capture SI or KHI, the high
clump metallicities would favor low values of $Q_{\mathrm{d}}$ and the
formation of planetesimals or planetary embryos. Although SI and KHI are
unlikely to compete with GI to determine the stability of dust clumps, the
interaction of three together should be explored in future works.

The stability of dust sizes outside the range of this study, in particular the
smaller sizes, should be explored in future work. Self-gravitating disks are
likely to be young, with less grain growth and coagulation to larger sizes, so
understanding the lower particle size range susceptible to gravitational
collapse will be important. Due to timestep limitations at dust sizes smaller
than $\mathrm{St}=0.1$, we do not include any simulations at these small
sizes.

###  6.2 Clump Survival

As shown, particles not only settle to the midplane, but they drift to the
dense gas structures generated in a $Q\approx 1$ disk and collapse of dust
into gravitationally bound clouds follows shortly thereafter. The power
spectrum of the turbulent eddies of gravitoturbulent disks peaks at scales
around the gas scale height (Cossins et al., 2009; Booth & Clarke, 2019) and
decreases down to the sizes of the clumps. Interactions with the high velocity
gas motions can potentially push some dust clumps into supersonic regime. This
is supported by the high collective motions measured for a few clumps in our
simulations, seen in Figures 3 and 4. While gas velocities can often reach the
supersonic regime close to the midplane but rarely in the midplane, the clumps
are pushed collectively but particle dispersions rarely exceed the escape
velocity.

In Figure 12, we measure the particle velocity dispersion in each clump. The
velocity dispersion with any given clump is generally less than 10% the sound
speed, but a small number of clumps at all masses have dispersions approaching
$0.5c_{s}$. The escape velocity $v_{\mathrm{esc}}$ of a particle from a
massive cloud can be found in terms of the gas sound speed through the
definition of the Bondi radius

$\left(\frac{v_{\mathrm{esc}}}{c_{s}}\right)^{2}=\frac{2GM_{\mathrm{c}}}{c_{s}^{2}r}=\frac{2R_{B}}{r},$
(40)

where $M_{\mathrm{c}}$ is the mass of the clump and $r$ is the distance away
from the center of the clump. For particles within the Hill radius $r\sim
R_{H}$, we arrive at an expression for the minimum escape velocity in terms of
the clump mass as a fraction of the thermal mass (see Fung et al., 2019) and
our unit mass $M_{\mathrm{th}}=M_{*}(H_{g}/R)^{3}=\pi^{3}\hat{M}_{0}$

$\frac{v_{\mathrm{esc}}}{c_{s}}=2^{1/2}3^{1/6}\left(\frac{M_{\mathrm{c}}}{M_{\mathrm{th}}}\right)^{1/3}=0.54\left(\frac{M_{\mathrm{c}}}{\hat{M}_{0}}\right)^{1/3}.$
(41)

As shown in Figure 12, some low mass clumps are above this threshold and
contain particles that are moving fast enough to fall apart, while the
majority of clumps either have sufficiently low particle velocity dispersions
or high enough masses to remain bound. This indicates that these clumps are
not just transient consequences of turbulent concentration, but persistent and
stable clouds of particles.

Figure 12: Velocity dispersion of particles in a clump versus the dust mass of
the clump. The gray dashed line indicates the threshold where the particle
velocity dispersion within the clump is the local escape velocity. A particle
clump above this line will likely be disrupted if it does not accrete more
particles or the particles in the clump become less excited.

###  6.3 Planetesimal and Embryo Formation

The formation of early planetary embryos and potentially planetary cores could
have important implications for the formation of gas giant planets. Planets
potentially carve out the gas and/or dust gaps on wide orbits observed in ALMA
images of young systems. Systems younger than one million years (i.e. Alves et
al., 2020; Segura-Cox et al., 2020) already show evidence of structures that
could be the result of planet formation. This would work together with pebble
accretion models, which require a solid seed mass of on the order of a
fraction of an Earth mass (Bitsch et al., 2015, 2019; Tanaka et al., 2020).
The accretion of additional pebbles and gas could take a few hundred to a
million years (Andama et al., 2022).

The streaming instability is an efficient mechanism to form planetesimals in
the inner few tens of au in protoplanetary disks, but less efficient at more
distant orbital separations. Thus, gravitoturbulent disks provide a means to
concentrate dust early such that the gravitational collapse could produce
solid bodies up to the size of a few Earth masses at several tens of au early
in the disk lifetime.

Gravitotrubulent disks are not the only way to form early planetesimals.
Besides the aforementioned streaming instability, secular gravitational
instabilities (Tominaga et al., 2020) could fill the same role in a young disk
while dust trapping in vortices via hydrodynamic instabilities could form the
seeds of early planet formation (Lyra et al., 2008a; Raettig et al., 2015;
Lyra & Umurhan, 2019).

###  6.4 Limitations

Without a radial pressure gradient, there is no radial drift included in these
simulations, which would add significant radial velocities to all particles,
but particularly the $\mathrm{St}=1$ particles. A radial pressure gradient
could also induce the streaming instability which would add further dust
diffusion. What effects radial drift ultimately has on the ability of dust
particles to concentrate should be further unexplored in 3D stratified
gravitoturbulent simulations.

While fully able to resolve the required length scales for gravitational
collapse of the gas, ideally these simulations would be able to resolve the
midplane dust layer with more than $\sim 10$ grid cells. Even more helpful
would be to model GI and particle-gas instabilities in the same simulation
such that one can continue the gravitational collapse of the dust down to
smaller scales. Shearing boxes with adaptively refined grids will be the most
useful tool to explore this scenario.

We considered grains all to be of the same size, which are also sizes most
likely to drift into the dense gas structures formed by GI. This is in part
for simplicity, but also because the interaction of multiple concurrent
species through backreaction effects is not yet fully understood (Zhu & Yang,
2021; Yang & Zhu, 2021). In reality, dust should initially have a size
distribution weighted heavily towards smaller grains similar to the
interstellar medium, where not much growth has occurred. The dust size
distribution in the ISM scales with size as a power law to the $-3.5$ (Mathis
et al., 1977). This would likely disfavor the mass distributions that rely on
the more massive dust particles, but dust growth from $\mu$m to mm sizes can
be efficient once the disk has formed (Birnstiel et al., 2016). Future work
will have to determine the mass of particle clouds that can be formed when
using a realistic dust size distribution. On the other hand, gravitationally
unstable disks are most likely to occur directly after star formation
(Küffmeier et al., 2018; Xu & Kunz, 2021), which may limit the amount of time
for the collisional growth of dust up to the sizes included in this work.
Whether collisional growth in these disks is efficient enough to create enough
large dust grains to reach the point where the particle collapse criterion
applies is a topic of continued research (Sengupta et al., 2019; Elbakyan et
al., 2020).

##  7 Conclusion

Gravitationally unstable disks are dominated by their gas content but may
concentrate significant quantities of solid material via gravitoturbulence
which triggers the self-gravitational collapse of concentrated dust. We use 3D
hydrodynamical simulations to compare how dust particles of different sizes
concentrate, collapse, and form bound objects. In these simulations, we
consider both simple dust drag without backreaction from the dust and a self-
consistent drag backreaction from the dust onto the gas. This helps to better
understand how dust and gas velocities are affected by one another and how the
mass distribution of bound solid dust clouds changes in each case. We
summarize the results of the paper below:

  1. 1.

Including self-gravitating dust in 3D simulations of disk gravitoturbulence
produces bound clumps of dust up to several $M_{\oplus}$ in mass at $\sim$50
au. Larger clumps are most efficiently formed when dust is of size
$\mathrm{St}=1$, while for smaller sizes $\mathrm{St}=0.1$, clumps form but
with a maximum mass roughly an order of magnitude lower.

  2. 2.

The formation of clumps at the small and intermediate sizes is consistent with
a particle collapse criterion based on the disruption of internal dust
diffusion. No bound clumps form when particles are size $\mathrm{St}=10$,
which is likely due to insufficient resolution of the critically unstable
length at this size.

  3. 3.

The effect of dust backreaction onto the gas can increase the perturbed
velocities of both dust and gas into the supersonic regime. The clumps as a
whole move at these high velocities, but the particle velocity dispersion
within the clumps remains low and the clumps are stable over extended periods
of time.

  4. 4.

The formation of solid objects as early as the period when a circumstellar
disk is self-gravitating could provide the seeds for planet formation. Growth
of these embryos via pebble accretion models could help to explain the
existence of early ring/gap structure in protoplanetary disks.

HB thanks Hubert Klahr and Konstantin Gerbig for valuable discussions and
comments. Simulations were conducted on the Pleiades supercomputer hosted by
the NASA Advanced Supercomputing (NAS) division as a part of the High-End
Computing (HEC) program. This research was supported by NASA Theoretical and
Computational Astrophysics Networks (TCAN) award 80NSSC19K0639 and discussions
with associated collaborators. CCY and ZZ acknowledge the support by NASA via
the Emerging Worlds program (grant number 80NSSC20K0347) and the Astrophysics
Theory Program (grant number 80NSSC21K0141). CCY is also grateful for the
support by NASA TCAN program (grant number 80NSSC21K0497). ZZ acknowledges
support from the National Science Foundation under CAREER Grant Number
AST-1753168.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, The Astrophysical Journal, 883, 192
  * Alves et al. (2020) Alves, F. O., Cleeves, L. I., Girart, J. M., et al. 2020, The Astrophysical Journal, 904, L6
  * Andama et al. (2022) Andama, G., Ndugu, N., Anguma, S. K., & Jurua, E. 2022, Monthly Notices of the Royal Astronomical Society, 510, 1298
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, The Astrophysical Journal Letters, 869, L41
  * Baehr et al. (2017) Baehr, H., Klahr, H., & Kratter, K. M. 2017, The Astrophysical Journal, 848, 40
  * Baehr & Zhu (2021a) Baehr, H., & Zhu, Z. 2021a, The Astrophysical Journal, 909, 135
  * Baehr & Zhu (2021b) —. 2021b, The Astrophysical Journal, 909, 136
  * Birnstiel et al. (2016) Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Science Reviews, 205, 41
  * Bitsch et al. (2019) Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, Astronomy & Astrophysics, 623, A88
  * Bitsch et al. (2015) Bitsch, B., Lambrechts, M., & Johansen, A. 2015, Astronomy & Astrophysics, 582, A112
  * Boley & Durisen (2010) Boley, A. C., & Durisen, R. H. 2010, The Astrophysical Journal, 724, 618
  * Bonnefoy et al. (2018) Bonnefoy, M., Perraut, K., Lagrange, A.-M., et al. 2018, Astronomy & Astrophysics, 618, A63
  * Booth & Clarke (2016) Booth, R. A., & Clarke, C. J. 2016, Monthly Notices of the Royal Astronomical Society, 458, 2676
  * Booth & Clarke (2019) —. 2019, Monthly Notices of the Royal Astronomical Society, 483, 3718
  * Boss (1997) Boss, A. P. 1997, Science, 276, 1836
  * Brandenburg (2003) Brandenburg, A. 2003, in Advances in Non-linear Dynamos, ed. A. Ferriz-Mas & M. Nunez (London: Taylor and Francis), 269
  * Brucy & Hennebelle (2021) Brucy, N., & Hennebelle, P. 2021, Monthly Notices of the Royal Astronomical Society, 503, 4192
  * Cadman et al. (2020) Cadman, J., Hall, C., Rice, K., Harries, T. J., & Klaassen, P. D. 2020, Monthly Notices of the Royal Astronomical Society, 498, 4256
  * Chandrasekhar (1963) Chandrasekhar, S. 1963, The Astrophysical Journal, 138, 1182
  * Cheetham et al. (2018) Cheetham, A. C., Bonnefoy, M., Desidera, S., et al. 2018, Astronomy & Astrophysics, 615, A160
  * Clarke et al. (2018) Clarke, C. J., Tazzari, M., Juhász, A., et al. 2018, The Astrophysical Journal Letters, 866, L6
  * Cossins et al. (2009) Cossins, P., Lodato, G., & Clarke, C. J. 2009, Monthly Notices of the Royal Astronomical Society, 393, 1157
  * Cuzzi et al. (2001) Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, The Astrophysical Journal, 546, 496
  * D’Alessio et al. (1998) D’Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, The Astrophysical Journal, 500, 411
  * Deng et al. (2017) Deng, H., Mayer, L., & Meru, F. 2017, The Astrophysical Journal, 847, 43
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. F. 1995, Icarus, 114, 237
  * Dullemond & Dominik (2004) Dullemond, C. P., & Dominik, C. 2004, Astronomy & Astrophysics, 421, 1075
  * Elbakyan et al. (2020) Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V. V., & Vorobyov, E. I. 2020, Astronomy & Astrophysics, 637, A5
  * Forgan (2019) Forgan, D. H. 2019, Monthly Notices of the Royal Astronomical Society, 485, 4465
  * Fung et al. (2019) Fung, J., Zhu, Z., & Chiang, E. I. 2019, The Astrophysical Journal, 887, 152
  * Gammie (2001) Gammie, C. F. 2001, The Astrophysical Journal, 553, 174
  * Gerbig et al. (2020) Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, The Astrophysical Journal, 895, 91
  * Gibbons et al. (2014) Gibbons, P. G., Mamatsashvili, G. R., & Rice, W. K. M. 2014, Monthly Notices of the Royal Astronomical Society, 442, 361
  * Gibbons et al. (2012) Gibbons, P. G., Rice, W. K. M., & Mamatsashvili, G. R. 2012, Monthly Notices of the Royal Astronomical Society, 426, 1444
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, Monthly Notices of the Royal Astronomical Society, 130, 125
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Janson et al. (2021) Janson, M., Gratton, R., Rodet, L., et al. 2021, Nature, 600, 231
  * Jeans (1902) Jeans, J. H. 1902, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 199, 1
  * Johansen et al. (2006) Johansen, A., Klahr, H., & Henning, T. 2006, The Astrophysical Journal, 636, 1121
  * Johansen et al. (2015) Johansen, A., Low, M.-M. M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. N. 2007, The Astrophysical Journal, 662, 627
  * Klahr & Schreiber (2020) Klahr, H., & Schreiber, A. 2020, The Astrophysical Journal, 901, 54
  * Klahr & Schreiber (2021) —. 2021, The Astrophysical Journal, 911, 9
  * Klee et al. (2017) Klee, J., Illenseer, T. F., Jung, M., & Duschl, W. J. 2017, Astronomy & Astrophysics, 606, A70
  * Klee et al. (2019) —. 2019, Astronomy & Astrophysics, 632, A35
  * Kratter & Lodato (2016) Kratter, K. M., & Lodato, G. 2016, Annual Review of Astronomy and Astrophysics, 54, 271
  * Kratter & Murray-Clay (2011) Kratter, K. M., & Murray-Clay, R. A. 2011, The Astrophysical Journal, 740, 1
  * Küffmeier et al. (2018) Küffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, Monthly Notices of the Royal Astronomical Society, 475, 2642
  * Kuznetsova et al. (2022) Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, arXiv preprint arXiv: 2202.05301v1, arXiv:2202.05301
  * Li et al. (2019) Li, R., Youdin, A. N., & Simon, J. B. 2019, The Astrophysical Journal, 885, 69
  * Lodato et al. (2019) Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, Monthly Notices of the Royal Astronomical Society, 486, 453
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, The Astrophysical Journal, 869, 17
  * Lyra et al. (2008a) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008a, Astronomy & Astrophysics, 491, L41
  * Lyra et al. (2008b) —. 2008b, Astronomy & Astrophysics, 479, 883
  * Lyra et al. (2017) Lyra, W., McNally, C. P., Heinemann, T., & Masset, F. 2017, The Astronomical Journal, 154, 146
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, Publications of the Astronomical Society of the Pacific, 131, 072001
  * Maire et al. (2016) Maire, A.-L., Bonnefoy, M., Ginski, C., et al. 2016, Astronomy & Astrophysics, 587, A56
  * Marois et al. (2010) Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. S. 2010, Nature, 468, 1080
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, The Astrophysical Journal, 217, 425
  * Mejía et al. (2005) Mejía, A. C., Durisen, R. H., & Pickett, M. K. 2005, The Astrophysical Journal, 619, 1098
  * Meru & Bate (2011a) Meru, F., & Bate, M. R. 2011a, Monthly Notices of the Royal Astronomical Society, 411, L1
  * Meru & Bate (2011b) —. 2011b, Monthly Notices of the Royal Astronomical Society, 410, 559
  * Meru & Bate (2012) —. 2012, Monthly Notices of the Royal Astronomical Society, 427, 2022
  * Morales et al. (2019) Morales, J. C., Mustill, A. J., Ribas, I., et al. 2019, Science, 365, 1441
  * Nayakshin (2010) Nayakshin, S. 2010, Monthly Notices of the Royal Astronomical Society Letters, 408, L36
  * Nelson (2006) Nelson, A. F. 2006, Monthly Notices of the Royal Astronomical Society, 373, 1039
  * Paardekooper (2012) Paardekooper, S.-J. 2012, Monthly Notices of the Royal Astronomical Society, 421, 3286
  * Paardekooper et al. (2011) Paardekooper, S.-J., Baruteau, C., & Meru, F. 2011, Monthly Notices of the Royal Astronomical Society, 416, L65
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration, T., Brandenburg, A., Johansen, A., et al. 2021, Journal of Open Source Software, 6, 2807
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * Pinte et al. (2018) Pinte, C., Price, D. J., Menard, F., et al. 2018, The Astrophysical Journal Letters, 860, L13
  * Raettig et al. (2015) Raettig, N., Klahr, H., & Lyra, W. 2015, The Astrophysical Journal, 804, 35
  * Rice et al. (2004) Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, Monthly Notices of the Royal Astronomical Society, 355, 543
  * Riols et al. (2020) Riols, A., Roux, B., Latter, H. N., & Lesur, G. 2020, Monthly Notices of the Royal Astronomical Society, 493, 4631
  * Safronov (1960) Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, Astronomy & Astrophysics, 597, A69
  * Schwarz et al. (2021) Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, The Astrophysical Journal Supplement Series, 257, 20
  * Segura-Cox et al. (2020) Segura-Cox, D. M., Schmiedeke, A., Pineda, J. E., et al. 2020, Nature, 586, 228
  * Sengupta et al. (2019) Sengupta, D., Dodson-Robinson, S. E., Hasegawa, Y., & Turner, N. J. 2019, The Astrophysical Journal, 874, 26
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astronomy & Astrophysics, 24, 337
  * Sheehan et al. (2020) Sheehan, P. D., Tobin, J. J., Federman, S., Megeath, S. T., & Looney, L. W. 2020, The Astrophysical Journal, 902, 141
  * Shi & Chiang (2013) Shi, J.-M., & Chiang, E. I. 2013, The Astrophysical Journal, 764, 20
  * Shi & Chiang (2014) —. 2014, The Astrophysical Journal, 789, 34
  * Shi et al. (2016) Shi, J.-M., Zhu, Z., Stone, J. M., & Chiang, E. I. 2016, Monthly Notices of the Royal Astronomical Society, 459, 982
  * Squires & Eaton (1991) Squires, K. D., & Eaton, J. K. 1991, Physics of Fluids A: Fluid Dynamics, 3, 1169
  * Tanaka et al. (2020) Tanaka, K. E. I., Zhang, Y., Hirota, T., et al. 2020, The Astrophysical Journal, 900, L2
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, The Astrophysical Journal Supplement, 257, 18
  * Tominaga et al. (2020) Tominaga, R. T., Takahashi, S. Z., & Inutsuka, S.-i. 2020, The Astrophysical Journal, 900, 182
  * Toomre (1964) Toomre, A. 1964, The Astrophysical Journal, 139, 1217
  * Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, The Astrophysical Journal Letters, 489, L179
  * van der Walt et al. (2011) van der Walt, S. J., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Visser et al. (2021) Visser, R. G., Dra̧żkowska, J., & Dominik, C. 2021, Astronomy & Astrophysics, 647, A126
  * Wahlberg Jansson & Johansen (2017) Wahlberg Jansson, K., & Johansen, A. 2017, Monthly Notices of the Royal Astronomical Society, 469, S149
  * Wang et al. (2021) Wang, S., Kanagawa, K. D., & Suto, Y. 2021, arXiv preprint arXiv: 2110.04827v1, 1
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, Monthly Notices of the Royal Astronomical Society, 180, 57
  * Weidenschilling (1980) —. 1980, Icarus, 44, 172
  * Xu & Kunz (2021) Xu, W., & Kunz, M. W. 2021, Monthly Notices of the Royal Astronomical Society, 502, 4911
  * Yang & Johansen (2016) Yang, C.-C., & Johansen, A. 2016, The Astrophysical Journal Supplement Series, 224, 39
  * Yang & Krumholz (2012) Yang, C.-C., & Krumholz, M. R. 2012, The Astrophysical Journal, 758, 48
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, The Astrophysical Journal, 868, 27
  * Yang et al. (2009) Yang, C.-C., Mac Low, M.-M., & Menou, K. 2009, The Astrophysical Journal, 707, 1233
  * Yang & Zhu (2021) Yang, C.-C., & Zhu, Z. 2021, Monthly Notices of the Royal Astronomical Society, 508, 5538
  * Youdin & Johansen (2007) Youdin, A. N., & Johansen, A. 2007, The Astrophysical Journal, 662, 613
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, The Astrophysical Journal, 869, L47
  * Zhu et al. (2012) Zhu, Z., Hartmann, L. W., Nelson, R. P., & Gammie, C. F. 2012, The Astrophysical Journal, 746, 110
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, The Astrophysical Journal, 801, 81
  * Zhu & Yang (2021) Zhu, Z., & Yang, C.-C. 2021, Monthly Notices of the Royal Astronomical Society, 501, 467
  * Zhu et al. (2019) Zhu, Z., Zhang, S., Jiang, Y.-f., et al. 2019, The Astrophysical Journal Letters, 877, L18