# A solution for the density dichotomy problem of Kuiper Belt objects with
multi-species streaming instability and pebble accretion

Manuel H. Cañas New Mexico State University, Department of Astronomy, PO Box
30001 MSC 4500, Las Cruces, NM 88001, USA  Wladimir Lyra New Mexico State
University, Department of Astronomy, PO Box 30001 MSC 4500, Las Cruces, NM
88001, USA  Daniel Carrera Department of Physics and Astronomy, Iowa State
University, Ames, IA, 50010, USA  Leonardo Krapp 51 Pegasus Fellow  Department
of Astronomy and Steward Observatory, University of Arizona, Tucson, Arizona
85721, USA  Debanjan Sengupta New Mexico State University, Department of
Astronomy, PO Box 30001 MSC 4500, Las Cruces, NM 88001, USA  Jacob B. Simon
Department of Physics and Astronomy, Iowa State University, Ames, IA, 50010,
USA  Orkan M. Umurhan NASA Ames Research Center, Space Sciences Division,
Planetary Sciences Branch, Moffatt Field, CA 94035, USA  Chao-Chin Yang (楊朝欽)
Department of Physics and Astronomy, University of Alabama, Box 870324,
Tuscaloosa, AL 35487-0324  Andrew N. Youdin Department of Astronomy and
Steward Observatory, University of Arizona, Tucson, Arizona 85721, USA  The
Lunar and Planetary Laboratory, University of Arizona

###### Abstract

Kuiper belt objects show an unexpected trend, whereby large bodies have
increasingly higher densities, up to five times greater than their smaller
counterparts. Current explanations for this trend assume formation at constant
composition, with the increasing density resulting from gravitational
compaction. However, this scenario poses a timing problem to avoid early
melting by decay of 26Al. We aim to explain the density trend in the context
of streaming instability and pebble accretion. Small pebbles experience
lofting into the atmosphere of the disk, being exposed to UV and partially
losing their ice via desorption. Conversely, larger pebbles are shielded and
remain more icy. We use a shearing box model including gas and solids, the
latter split into ices and silicate pebbles. Self-gravity is included,
allowing dense clumps to collapse into planetesimals. We find that the
streaming instability leads to the formation of mostly icy planetesimals,
albeit with an unexpected trend that the lighter ones are more silicate-rich
than the heavier ones. We feed the resulting planetesimals into a pebble
accretion integrator with a continuous size distribution, finding that they
undergo drastic changes in composition as they preferentially accrete silicate
pebbles. The density and masses of large KBOs are best reproduced if they form
between 15 and 22 AU. Our solution avoids the timing problem because the first
planetesimals are primarily icy, and 26Al is mostly incorporated in the slow
phase of silicate pebble accretion. Our results lend further credibility to
the streaming instability and pebble accretion as formation and growth
mechanisms.

Pebble accretion, planet formation, streaming instability, Kuiper belt
objects, Pluto

††journal: PSJ

##  1 Introduction

The formation of kilometer-sized planetesimals from $\mu$m-sized dust
particles is a process still relatively poorly understood, despite significant
advances (Chiang & Youdin, 2010; Johansen et al., 2014; Lesur et al., 2022).
Sub-$\mu$m interstellar grains grow by collisional coagulation (Safronov,
1972; Nakagawa et al., 1981; Tominaga et al., 2021), but accumulated evidence
from laboratory experiments (Blum & Wurm, 2008; Güttler et al., 2010),
numerical simulations (Güttler et al., 2009; Geretshauser et al., 2010; Zsom
et al., 2010) and observations (Pérez et al., 2015; Tazzari et al., 2016;
Carrasco-González et al., 2019) suggest that growth is inefficient beyond
millimeter (mm) and centimeter (cm) size (hereafter called “pebbles”) due to
bouncing, fragmentation, and drift (Dullemond & Dominik, 2005; Brauer et al.,
2008; Krijt et al., 2015). Therefore, an alternative route to planetesimal
formation is needed, perhaps via direct gravitational collapse of pebble
clouds Youdin & Shu (2002).

A promising mechanism to overcome collisional barriers, and trigger
gravitational collapse, is the streaming instability (Youdin & Goodman, 2005;
Youdin & Johansen, 2007; Johansen & Youdin, 2007). In this feedback mechanism,
particle overdensities trigger a gas flow that enhances the particle
overdensities (Youdin & Goodman, 2005; Squire & Hopkins, 2020). In the typical
conditions of gas-rich protoplanetary disks, the overdensities caused by
streaming instability have been shown to be massive enough to result in the
formation of gravitationally bound objects (Johansen et al., 2007; Yang &
Johansen, 2014; Carrera et al., 2015; Simon et al., 2016, 2017; Schäfer et
al., 2017; Yang et al., 2017; Li et al., 2019; Abod et al., 2019; Nesvorný et
al., 2019; Li & Youdin, 2021). Growth to planetary masses then continues via
accretion of pebbles, as they lose energy from aerodynamical drag while
passing a protoplanet, significantly increasing its effective accretional
cross-section (Lyra et al., 2008a; Johansen & Lacerda, 2010; Ormel & Klahr,
2010; Lambrechts & Johansen, 2012; Lyra et al., 2023).

Figure 1: Density distribution of Kuiper belt objects. Small KBOs have density
around 0.5 g cm-3 (with some dispersion), implying porosity, whereas larger
bodies have increasingly higher density. Values from Bierson & Nimmo (2019,
and references therein), Morgado et al. (2023), and Pereira et al. (2023).

Given that the current models of planet formation involve the concentration
and accumulation of neighboring pebbles, it is surprising that Kuiper belt
objects (KBOs) demonstrate a large range of densities (from 0.5 g cm-3 to 2.6
g cm-3, Brown 2013), with a trend that the small KBOs ($\lesssim 2\times
10^{-2}$ Pluto masses) have a near-constant density of $0.5$ g cm-3 and larger
KBOs have an increasingly larger density that reaches $\approx 2.5\text{ g
cm}^{-3}$ (Stansberry et al., 2006; Brown & Butler, 2017; Grundy et al., 2007;
Brown, 2012, 2013; Stansberry et al., 2012; Barr & Schwamb, 2016; Ortiz et
al., 2017; Grundy et al., 2019). We show in Fig. 1 the mass vs density
relation; the data for mass and density are from Bierson & Nimmo (2019, and
references therein), except for Quaoar, for which we use the more recently
evaluated mass of $\left(1.20\pm 0.05\right)\times{10^{24}}$ g from Morgado et
al. (2023) and Pereira et al. (2023).

One plausible explanation is that smaller KBOs are more porous, and as a KBO
grows, porosity is removed through gravitational compaction. Assuming water
ice (density $\approx$1 g cm-3) is the lowest density material of substantial
abundance in the composition of the small objects, the bulk density of 0.5 g
cm-3 implies a porosity of at least 50% for them. The larger objects (e.g.
Eris, Pluto, and Triton) should be less porous even in the unlikely case of
100% rocky composition. Interestingly, Baer et al. (2011) find a similar trend
in the asteroid belt, with the largest asteroids showing porosities less than
20%, with an abrupt change for diameters under 300 km where a range of
porosities from 0-70% is seen.

Gravitational compaction at constant composition was explored by Bierson &
Nimmo (2019). Assuming a rock mass fraction of 70%, these authors successfully
match 15 of 18 KBOs within two standard deviations allowed by their model (or
11 within one standard deviation). While promising, this scenario poses a
timing problem requiring KBOs to be formed at least four million years after
the formation of the Calcium-Aluminum-rich Inclusions (CAIs). The reason for
this time constraint is that heat from the decay of the short-lived
radioactive isotope 26Al, with a half-life of 0.7 Myr (Norris et al., 1983),
would act to melt the objects and remove porosity. Thus, if KBOs formed early
in the solar nebula with a high rock fraction, we would expect to see high-
density, low-mass KBOs, which are not seen.

The required time delay, however, is unlikely because four million years is
potentially longer than the disk lifetime (e-folding time 2.5 Myr, e.g. Ribas
et al., 2014), which would contradict evidence for KBO formation by
gravitational collapse in a gas disk (Nesvorný et al., 2010, 2019; Lisse et
al., 2021), and Arrokoth needing nebular drag for the two lobes to come into
contact (McKinnon et al., 2020; Lyra et al., 2021). The special timing also
contradicts indication that planets might form early in protoplanetary disks
(e.g. ALMA Partnership et al., 2015; Manara et al., 2018; Tobin et al., 2020;
Yamato et al., 2023; Sai et al., 2023). Thus we search for another explanation
for the density trend, where the large change in density stems from
compositional differences between large and small KBOs, with large KBOs
containing a higher rock fraction than their smaller counterparts.

We consider formation beyond the water ice line, the region of the disk where
it is cold enough for water to condense into solid ice grains. While
coagulation of ices and silicates beyond the snowline should produce pebbles
of homogeneous composition, a compositional difference in pebbles should be
expected from preferential depletion of ices in small grains. Due to
turbulence and bulk vertical motions in the disk – caused by, e.g., the
vertical shear instability (VSI, Nelson et al., 2013; Lin & Youdin, 2015;
Lesur et al., 2022) or the magnetorotational instability (MRI, Balbus &
Hawley, 1998; Lyra et al., 2008b; Johansen et al., 2009; Bai & Stone, 2011;
Simon et al., 2018) if the disk is sufficiently ionized –, smaller grains are
lofted up in the atmosphere of the disk. There, they are exposed to stellar UV
irradiation, which causes removal of ice by photodesorption (photosputtering,
Westley et al., 1995a, b; Bergin et al., 1995; Andersson et al., 2006;
Andersson & van Dishoeck, 2008; Öberg et al., 2009; Krijt et al., 2016, 2018;
Potapov et al., 2019; Powell et al., 2022). In an optically thin environment,
the removal rate produced by solar UV irradiation is estimated at 4 mm per Myr
(Harrison & Schoen, 1967) at 10 AU (similar to the coagulation rate at steady
state), obliterating small icy grains and perhaps explaining their apparent
absence in debris disks (Grigorieva et al., 2007; Stuber & Wolf, 2022).
Evidence for this process in a gas-rich disk is seen in the presence of water
vapor at very low temperatures in circumstellar disks (Hogerheijde et al.,
2011), which is interpreted as water molecules released by non-thermal
desorption into the gas phase, before quickly recondensing as amorphous ice
(Ciesla, 2014).

The height at which the disk becomes optically thick to UV is estimated at 3
scale heights for a primordial disk where all the dust is in sub-$\mu$m size
(Krijt et al., 2018). Yet, the optical depth at a given height decreases with
time as the sub-$\mu$m grains responsible for the opacity coagulate into
larger grains that settle toward the midplane. As most of the solid masses is
converted into larger grains, the mass remaining in sub-$\mu$m grain species
decreases significantly. When the dust-to-gas ratio of these sub-$\mu$m grain
species decreases to $10^{-4}$, typical of Class I-II disks, the gas becomes
optically thin to UV as close to the midplane as 1.5 scale heights, a layer
where even low levels of turbulence would suffice to loft sub-mm sized grains
into. Bulk motions due to the VSI would make it even easier to get these
particles to the UV-irradiated layer, and cosmic ray desorption (Silsbee et
al., 2021; Sipilä et al., 2021; Rawlings, 2022) would further enhance the loss
of ice coating from grains. As a result of these processes, smaller grains
should have less ice than larger grains that reside near the disk midplane,
shielded from UV and cosmic rays. A bimodal distribution of grains is then
expected, with smaller ice-poor pebbles (hereafter “silicate” pebbles), and
larger pebbles composed of dirty ice (hereafter “icy” pebbles).

Considering the results of Drażkowska & Dullemond (2014); Carrera et al.
(2015); Yang et al. (2017), and Li & Youdin (2021), where the effectiveness of
the streaming instability is tested amongst various grain sizes, we expect the
(larger) icy pebbles to be more conducive to the streaming instability, while
the (smaller) silicate pebbles are more tightly coupled to the gas and hence
participate less (see also Yang & Zhu, 2021). The outcome is that
planetesimals formed beyond the ice line should be mostly icy. If this
hypothesis is correct, there remains the question of how objects in the Kuiper
belt achieve densities beyond $2.0\text{ g cm}^{-3}$. The answer might be
through the subsequent processes of pebble accretion and compaction. Pebble
accretion has a strong dependency on both grain size and planetesimal mass
(Lyra et al., 2023) such that, as a planetesimal grows, the favorable grain
size for pebble accretion also changes. This allows protoplanets to accrete a
wealth of silicates through pebble accretion. The accretion of silicates, in
combination with limited removal of porosity from gravitational compaction,
would result in a natural trend where low-mass objects (i.e., those that
failed to accrete pebbles) have low densities, and the bodies increase in
density as a higher fraction of their mass comes from pebble accretion. This
is the scenario we explore in this work.

This paper is organized as follows. In Sect. 2 we discuss our model, including
the hydrodynamical simulation of the streaming instability and the numerical
integration of polydisperse pebble accretion. In Sect. 3 we describe our
results of both the streaming instability and pebble accretion. In Sect. 4 we
discuss the interpretation of our results and compare it with other works that
sought to explain the density trend of KBOs. In Sect. 5 we discuss the
limitations of the model and future works. We conclude the work in Sect. 6.

Table 1: Symbols used in this work. Symbol | Definition | Description |  | Symbol | Definition | Description
---|---|---|---|---|---|---
$t$ |  | time |  | $r_{0}$ |  | arbitrary origin of shearing box
$x,y,z$ |  | Cartesian coordinates |  | $p$ | Eq. (8) | gas pressure
$\varepsilon$ | $\rho_{d}/\rho_{g}$ | dust-to-gas density ratio |  | $c_{s}$ | Eq. (11) | the sound speed
$Z$ | $\frac{\int\rho_{d}dV}{\int\rho_{g}dV}$ | dust-to-gas mass ratio |  | $dV$ | $dx\,dy\,dz$ | volume infinitesimal
$v_{K}$ | $\varOmega_{0}r_{0}$ | Keplerian velocity |  | $\Delta v$ | $\eta v_{K}$ | gas velocity reduction
$\gamma$ |  | adiabatic index |  | $\tau$ | Eq. (A2) | particle friction time
$v_{\rm th}$ | $\sqrt{\frac{8}{\pi}}c_{s}$ | thermal velocity |  | $\rho_{\bullet}$ |  | pebble material density
$k_{B}$ |  | Boltzmann constant |  | ${\boldsymbol{v}}$ | $(v_{x},v_{y},v_{z})$ | pebble velocity
$T$ |  | temperature |  | $\lambda_{\rm SI}$ | $\eta r$ | streaming instability length scale
$m_{H}$ |  | atomic mass unit |  | $H_{d}$ | Eq. (13) | pebble scale height
$\varOmega_{K}$ | Eq. (12) | Keplerian frequency |  | $N$ | $N_{x},N_{y},N_{z}$ | box resolution
$M$ |  | stellar mass |  | $\widetilde{G}$ | Eq. (18) | dimensionless gravitational constant
$\delta$ |  | dimensionless pebble diffusion |  | $\rho_{d}$ |  | volume density of pebbles
$\Delta$ |  | mesh spacing |  | $a$ |  | pebble radius
$\phi$ |  | porosity |  | $\rm{St}$ | $\varOmega_{0}\tau$ | Stokes number of pebbles
$R$ |  | protoplanet radius |  | $\varPhi$ |  | self-gravitational potential
$\rho_{s}$ |  | silicate density |  | $\mu$ |  | mean molecular weight
$\rho_{i}$ |  | ice density |  | $L$ |  $L_{x}$=$L_{y}$=$L_{z}$ | length of shearing box sides
$H$ | $c_{s}/\varOmega_{0}$ | gas scale height |  | $\alpha$ | Eq. (14) | dimensionless gas viscosity
${\boldsymbol{u}}$ | $(u_{x},u_{y},u_{z})$ | gas velocity in Cartesian coordinates |  | $N_{p}$ |  | number of numerical particles
$\varOmega_{0}$ |  | Keplerian frequency at $r_{0}$ |  | $Q$ | Eq. (19) | Toomre $Q$ parameter
$\rho_{g}$ |  | gas density |  | $\rho_{R}$ | Eq. (20) | Roche density
$h$ | $H/r_{0}$ | disk aspect ratio |  | $\dot{M}$ |  | mass accretion rate
$\eta$ | Eq. (6) | dimensionless pressure support |  | $r$ |  | astrocentric distance
$u_{0}$ | $-\frac{3}{2}\varOmega_{0}x$ | Keplerian shear flow |  | $M_{p}$ |  | protoplanet mass
$\theta$ | Eq. (30) | Heaviside step function |  | $\xi$ | Eq. (31) | boxcar function
$G$ |  | gravitational constant |  | $\Pi$ | $\Delta v/c_{s}$ | dimensionless gas velocity reduction
$V$ |  | protoplanet volume |  | $m_{j}$ |  | mass of protoplanet in constituent $j$
$\rho$ |  | protoplanet density |  | $V_{j}$ | $m_{j}/\rho_{j}$ | volume occupied by constituent $j$
$F_{j}$ | $m_{j}/M_{p}$ | mass fraction of constituent $j$ |  | $f_{j}$ | $V_{j}/V$ | volume fraction of constituent $j$
$t_{\rm settle}$ | Eq. (33) | pebble settling time |  | [26Al] |  | isotopic abundance of 26Al
$\mathcal{H}$ |  | heat production rate per mass |  | $c_{p}$ |  | heat capacity of water ice
$t_{1/2}$ |  | half-life of 26Al |  | $\lambda$ | $\ln(2)/t_{1/2}$ | decay constant of 26Al

##  2 The Model

We use the Pencil Code111http://pencil-code.nordita.org (see Pencil Code
Collaboration et al., 2021, and references therein for details) to model the
formation of the first planetesimals by streaming instability of ices and
silicates. We employ the shearing-box approximation (Brandenburg et al., 1995;
Hawley et al., 1995), centered at an arbitrary position $r_{0}$ and orbiting
at the corresponding Keplerian frequency $\varOmega_{0}$. Our model includes
gas, solids, vertical stellar gravity, and dust self-gravity; the latter
allows a concentration of solids to collapse into planetesimals once the Roche
density is exceeded.

###  2.1 Equations of motion

We consider an isothermal gas disk such that any gain in internal energy is
considered to be radiated away effectively. While the gas is solved on a
uniform mesh, the pebbles are represented by Lagrangian particles. We ignore
the selfgravity of the gas (explained a posteriori in Sect. 2.3), but consider
the gravity of the dust grains. Under this assumption, the equations of motion
for the system are

$\displaystyle\frac{{\mathcal{D}}\rho_{g}}{{\mathcal{D}}t}$ $\displaystyle=$
$\displaystyle-\rho_{g}{\boldsymbol{\nabla}}\cdot{\boldsymbol{u}},$ (1)
$\displaystyle\frac{{\mathcal{D}}{\boldsymbol{u}}}{{\mathcal{D}}t}$
$\displaystyle=$
$\displaystyle-\frac{1}{\rho_{g}}{\boldsymbol{\nabla}}{p}-2\varOmega_{0}\hat{{\boldsymbol{z}}}\times{\boldsymbol{u}}+\frac{3}{2}\varOmega_{0}u_{x}\hat{{\boldsymbol{y}}}-\varOmega_{0}^{2}{\boldsymbol{z}}+{\boldsymbol{f}}_{g}$
(2) $\displaystyle+2\varOmega_{0}\Delta v\hat{{\boldsymbol{x}}},$
$\displaystyle\frac{d{\boldsymbol{x}}}{dt}$ $\displaystyle=$
$\displaystyle{\boldsymbol{v}}-\frac{3}{2}\varOmega x\hat{{\boldsymbol{y}}},$
(3) $\displaystyle\frac{d{\boldsymbol{v}}}{dt}$ $\displaystyle=$
$\displaystyle-2\varOmega_{0}\hat{{\boldsymbol{z}}}\times{\boldsymbol{v}}+\frac{3}{2}\varOmega_{0}v_{x}\hat{{\boldsymbol{y}}}-\varOmega_{0}^{2}{\boldsymbol{z}}+{\boldsymbol{f}}_{d}-{\boldsymbol{\nabla}}{\varPhi},$
(4) $\displaystyle\nabla^{2}{\varPhi}$ $\displaystyle=$ $\displaystyle 4\pi
G\rho_{d}.$ (5)

Here, $\rho_{g}$ is the gas density, $t$ is time, ${\boldsymbol{u}}$ is the
gas velocity, $p$ is the gas pressure, ${\boldsymbol{x}}$ and
${\boldsymbol{v}}$ are the position and velocity vectors of the pebbles,
respectively, $\rho_{d}$ is the volume density of pebbles, $\varPhi$ is the
selfgravity potential, $G$ is the gravitational constant, and ($x$,$y$,$z$)
are the local Cartesian coordinates of the shearing box. The terms $f_{d}$ and
$f_{g}$ are the drag force and its backreaction, respectively, explained in
Appendix A.

The quantity $\Delta v=\eta v_{K}$ is the orbital velocity reduction of the
gas due to the pressure support, and $\eta$ is related to the global-scale
radial pressure gradient (Nakagawa et al., 1986),

$\eta\equiv-\frac{1}{2\rho_{g}\varOmega_{K}^{2}r}\frac{\partial p}{\partial
r}.$ (6)

A table of symbols is provided in

table:symbols. We have defined the operator

$\frac{{\mathcal{D}}\ }{{\mathcal{D}}t}=\frac{\partial}{\partial
t}+u_{0}\frac{\partial}{\partial
y}+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}},$ (7)

where $u_{0}=-\left(3/2\right)\varOmega_{0}x$ is the linearized Keplerian
shear flow. The pressure is related to the sound speed $c_{s}$ by the equation
of state

$p=\rho_{g}c_{s}^{2}/\gamma,$ (8)

which closes the system. Here $\gamma=1$ is the adiabatic index for an
isothermal gas.

The gas is initially in hydrostatic equilibrium between its own pressure and
the vertical gravity from the central star. This results in the gas density
having a vertical profile

$\rho_{g}\left(z\right)=\rho_{0}\exp\left(-\frac{z^{2}}{2H^{2}}\right),$ (9)

where $\rho_{0}$ is the midplane gas density and

$H\equiv\frac{c_{s}}{\varOmega_{0}}$ (10)

is the gas scale height.

###  2.2 Box Domain and Resolution

We choose the length of our box such that the scales of the streaming
instability, $\lambda_{\rm SI}\approx r\eta$, are well contained within the
box. Given Eq. (10), we can rewrite this length in terms of the scale height
as $\lambda_{\rm SI}\approx hH$, where the disk aspect ratio $h\equiv
H/r_{0}$. The sound speed is, for constant adiabatic index $\gamma$ and mean
molecular weight $\mu$, solely dependent on temperature

$c_{s}^{2}=\frac{\gamma k_{B}T}{\mu m_{H}}.$ (11)

We substitute the values $\gamma=1$ for the adiabatic index of an isothermal
flow, $\mu=2.3$ for the mean molecular weight (i.e., 5 $\rm{H}_{2}$ for every
2 He), $k_{B}$ and $m_{H}$ are the Boltzmann constant and atomic mass unit,
respectively. The disk temperature $T$ is found using the temperature relation
of Chiang & Goldreich (1997). In this model at $r=45\text{\,AU}$, the disk
temperature is $T\approx 30\text{ K}$, which results in a sound speed of
$c_{s}\approx 0.3\text{ km s}^{-1}$. The Keplerian frequency is

$\varOmega^{2}_{K}\equiv\frac{GM}{r^{3}},$ (12)

where $M$ is the mass of the host star, and $r$ is the astrocentric distance.
At 45 AU for a solar mass star, the Keplerian frequency is
$\varOmega_{K}=6.596\times 10^{-10}\text{ s}^{-1}$ resulting in a scale height
of $H\approx 3.3\text{\,AU}$. Thus, we have the characteristic length of the
streaming instability being $\eta r\approx 0.07H$, placing an upper limit for
the distance between mesh cells, $\Delta$. Even for low resolutions of 16 mesh
cells along each axis, a cubic box of sides $L=0.2H$ ($0.66\text{\,AU}$ for
$H=3.3\text{\,AU}$) will ensure that the streaming instability is resolved.
Consequently, this is the chosen size of our box.

While the size of our box is largely determined by $\lambda_{\rm SI}$, the
required resolution of our mesh is primarily constrained by the dust scale
height, $H_{d}$, which is determined as a balance between vertical stellar
gravity and turbulent diffusion

$H_{d}=H\left(1+\frac{\mathrm{St}}{\delta}\right)^{-1/2},$ (13)

where $\delta$ is a dimensionless diffusion parameter (Dubrulle et al., 1995;
Johansen & Klahr, 2005; Youdin & Lithwick, 2007; Lyra & Lin, 2013) related to
(in the case of isotropic turbulence, Yang et al. 2018) the dimensionless
Shakura-Sunyaev $\alpha$ viscosity (Shakura & Sunyaev, 1973)

$\alpha\equiv\frac{\langle\rho_{g}\,\delta v_{x}\,\delta v_{y}\rangle}{p}.$
(14)

Here $\delta v_{i}$ is the turbulent deviation from the mean flow in the
direction $i$. The $\alpha$-viscosity parameter is related to $\delta$ by
(Youdin & Lithwick, 2007)

$\delta=\alpha\left(1+\mathrm{St}^{2}\right)^{-1}.$ (15)

Eq. (13) is strictly valid only when particles are completely passive. While
this is not the case for simulations of the streaming instability (as feedback
is essential to this instability), we adopt this formula for simplicity. Also,
Eq. (14) ignores magnetic stresses.

We do not consider external turbulence in this simulation; the only turbulence
present is due to the streaming instability itself, which in turn produces
$\alpha$-viscosity around $\alpha\approx 10^{-5}$. 222While that is
inconsistent with our assumption that external turbulence is necessary to loft
the silicate grains up in the disk atmosphere to remove the ice (see Sect. 5
on limitations of the model), the low turbulence only increases the number of
small grains near the midplane, and is thus a conservative approach to what
upper limit for the silicate fraction we should expect in the resulting
planetesimals. For Stokes number of $\mathrm{St}=0.5$, this results in a dust
scale height of $H_{d}\approx H/250$. We want to resolve this layer with at
least 6 mesh points (the size of a stencil), and to this end, we have set the
resolution of our simulation to be $N_{x}=N_{y}=N_{z}=256$ resulting in a
points-per-scale-height value of $H/\Delta=1280$.

###  2.3 Simulation Parameters

We set the dust-to-gas mass ratio (metallicity) to be $Z=0.03$, slightly above
Solar, which is known to trigger strong clumping by the streaming instability
(Carrera et al., 2015; Yang et al., 2017; Li & Youdin, 2021). The bulk mass of
solids is evenly split into $N_{p}=1.536\times 10^{7}$ numerical super-
particles; each particle’s mass is the cumulative mass of the pebbles it
represents, but the aerodynamical behavior of a particle is identical to that
of a single pebble. We chose Stokes numbers such that ices and silicates are
on cm and sub-mm size scales, respectively. Given

$\displaystyle\mathrm{St}$ $\displaystyle\approx$ $\displaystyle
0.1\left(\frac{a}{1\rm\ cm}\right)\left(\frac{\rho_{\bullet}}{1\rm\ g\
cm^{-3}}\right)\left(\frac{T}{35\rm\
K}\right)^{-1/2}\left(\frac{\rho_{g}}{10^{-13}\rm\ g\ cm^{-3}}\right)^{-1}$

we assign half the particles a Stokes number of $\mathrm{St}=0.5$,
representing ices, and the other half a Stokes number of $\mathrm{St}=5\times
10^{-3}$, representing silicates.

For the pebble drift, we use the dimensionless parameter defined by Bai &
Stone (2010)

$\Pi\equiv\frac{\Delta\,v}{c_{s}},$ (17)

and set it to $\Pi=0.05$ (note also that $\Pi c_{s}=\eta v_{K}$). The relative
strength of self-gravity to tidal shear is determined by another dimensionless
parameter

$\widetilde{G}\equiv\frac{4\pi G\rho_{0}}{\varOmega_{0}^{2}},$ (18)

and is related to the usual Toomre $Q$ parameter (Safronov, 1960; Toomre,
1964) by

$Q=\sqrt{\frac{8}{\pi}}\widetilde{G}^{-1}.$ (19)

The constants on the right-hand side of Eq. (5) are non-dimensionalized and in
code units are set to be $4\pi G=0.1$, consequently setting the gravitational
constant in code units to be $G=0.1/4\pi$. Furthermore, with
$\rho_{0}=\varOmega_{0}=1$, then $\widetilde{G}=0.1$ and $Q\approx 16$,
justifying the exclusion of self-gravity from the equations of motion for gas.
Gravitational collapse of the pebbles ensues once the pebble density within a
mesh cell exceeds the Roche density,

$\rho_{R}\equiv\frac{9\varOmega_{0}^{2}}{4\pi G}=9\rho_{0}\widetilde{G}^{-1}.$
(20)

Once this happens, all the particles within a sphere of radius equal to one
mesh cell are removed and replaced by a sink particle (Johansen et al., 2015;
Schäfer et al., 2017) that has a Stokes number akin to that of a planetesimal
and thus no longer feels gas drag. These sink particles can continue to
accrete pebbles.

###  2.4 The composition of icy pebbles

We expect from sticking velocity arguments that the large pebbles will be
mostly icy. This is because, starting from a bare silicate nucleus, and
letting the grain aggregate icy and silicate monomers, a higher ice sticking
velocity means that the grains would grow progressively icier as the silicate
nucleus is diluted in the ice, forming ice-mantled silicate grains. If a 1 mm
bare silicate grain accretes only ice until cm-size, the volume ratio implies
that the grain will be, per volume (mass), only 0.1% (0.3%) silicate. Allowing
for a smaller bare silicate nucleus will increase the ice ratio, so we
consider the bare silicate to only be a trace in the final composition. The
main constraint will come from the concurrent coagulation of silicates and
ices, which will depend on the coagulation efficiency (Lambrechts & Johansen,
2014).

The sticking velocity for water ice is expected to be about 10 times higher
than silicates (Wada et al., 2009; Gundlach & Blum, 2015). This conclusion was
challenged by Musiolik & Wurm (2019), who find that the high surface energy of
ice grains applies only in a narrow temperature range near the ice line; below
175 K, the surface energy of water ice is akin to that of silicates. However,
the translation from surface energy to sticking velocity is not
straightforward, as collision outcomes also depend strongly on porosity, mass
ratio, and impact parameter, with sticking collisions possible with velocities
up to 100 m/s (Planes et al., 2020, 2021). Also, Musiolik (2021) reports
higher surface energies for irradiated ice due to the development of a liquid
microfilm, depending on the width of the ice crust. We therefore consider the
value suggested by Gundlach & Blum (2015), 10 m/s, for the sticking velocity
for ice, and 1 m/s for silicates. With 10 times higher efficiency of ice
coagulation compared to silicate coagulation, we expect the matrix of the
larger pebbles to be about 90% ice by mass (silicates are denser than ice, by
a factor 3, but ices are more abundant than silicates, by a factor 4, so the
two effects almost neatly cancel).

###  2.5 Pebble Accretion

Figure 2: Internal density (upper panels) and ice volume fraction (lower
panels) of accreted pebbles for the three models. (Left) The power law
distribution in which the smallest grains are pure silicates and the largest
pebbles are pure ice, with the trend that larger solids are icier. (Middle)
The bimodal distribution in which pebbles between $\mu$m and mm are pure
silicates, and pebbles above this size are increasingly icier with pebbles at
one centimeter being pure ice. (Right) The constant distribution, in which
pebbles of all grain sizes have a 0% ice volume fraction and a density of
$\rho_{\bullet}=3$ g cm-3.

Figure 3: The accretion rates for the three dust distributions at $t=0$. The
orange lines are the accretion rates in the Hill regime, the magenta lines are
the accretion rates in the Bondi regime, and the green lines are the accretion
rates in the geometric/focusing regime. The scatter points are color-coded
corresponding to the ice volume fraction of the pebble of most efficient
accretion at the given protoplanet mass. (Left) Accretion rates of the power
law distribution, which accrete no silicates except for a small window in the
transition from geometric to Bondi, where roughly 20% of the material being
accreted is silicate. (Middle) The bimodal distribution begins accreting
roughly 50% ice and 50% silicates in the geometric regime, but then accretes
nearly 100% silicates upon entering the Bondi regime. As a planetesimal grows
through Bondi accretion, it begins to favor larger pebbles, thereby increasing
the ice volume fraction that is being accreted, finally accreting mostly ices
in the Hill regime. (Right) Accretion rates for the constant distribution,
where all pebbles being accreted have a 0% ice volume fraction.

Once gravitationally bound objects are produced, they continue to grow by
accreting pebbles (Lyra et al., 2008a, 2009a, 2009b; Ormel & Klahr, 2010;
Lambrechts & Johansen, 2012). Yet, for the objects produced, of the order of
100 km in radius, the pebble accretion rates are much longer than we can model
in a hydrodynamical simulation. Therefore, we take the planetesimal population
produced in the streaming instability simulation, and feed them into a
separate, relatively inexpensive, pebble accretion integrator, that solves the
pebble accretion equations on evolutionary timescales.

The pebble accretion model we adopt considers a polydisperse distribution of
pebbles, as recently developed by Lyra et al. (2023), who found efficient
pebble accretion on top of the direct products of streaming instability.
Particularly, in the polydisperse model, the pebbles can have different
internal density, which we will use as different composition.

We consider three different models for pebble accretion, each differing in
their respective internal density and ice volume fractions. The first model is
a control, modeled as a power law

$\rho_{\bullet}(a)=\rho_{i}\left(\frac{a}{a_{\rm max}}\right)^{-q_{1}},$ (21)

with

$q_{1}\equiv\frac{\log(\rho_{s}/\rho_{i})}{\log(a_{\rm max}/a_{\rm min})}>0$
(22)

so that the smallest grain has a density $\rho_{s}$, and $\rho_{\bullet}$
decreases in internal density with increasing grain size such that the largest
particle has a density $\rho_{i}$. We choose for this model $a_{\rm
min}=1\,\mu$m and $a_{\rm max}=1\,$cm.

The second model is our physically motivated bimodal distribution

$\rho_{\bullet}(a)=\begin{cases}\rho_{s},&\text{if }a<a_{0}\\\
\rho_{i}\left(\frac{a}{a_{\rm max}}\right)^{-q_{2}},&\text{if }a\geq
a_{0}\end{cases}$ (23)

with

$q_{2}\equiv\frac{\log(\rho_{s}/\rho_{i})}{\log(a_{\rm max}/a_{0})}>0.$ (24)

We choose for this model $a_{0}=1$ mm; that is, particles between the sizes
$1\,\mu\text{m}\leq a\leq 1\,\text{ mm}$ are assigned internal density
$\rho_{s}$. Pebbles beyond 1 mm in size then follow a power law such that the
largest pebble has internal density $\rho_{i}$.

The last model assumes that all pebbles are silicates, with internal density
$\rho_{s}$. It could represent a streaming instability filament of bare
silicate pebbles and, as we will show, is necessary to reproduce Eris.

The models are illustrated in Fig. 2, where each column represents the
different models. The upper row shows the internal density of the pebbles, and
the bottom row is the ice volume fraction of the pebbles as a function of
grain size. The ice volume fraction is calculated according to

$\displaystyle\rho=\rho_{i}f_{i}+\rho_{s}f_{s}$ (25) $\displaystyle
f_{i}+f_{s}+f_{v}=1$ (26)

where $f_{j}\equiv V_{j}/V$ is the volume fraction of component $j$; $V_{j}$
is the volume occupied by the component, and $V$ is the total volume. Here
$f_{v}$ is the fractional volume of empty space, or “porosity”. For compact
pebbles, $f_{v}=0$ (and hence $\rho=\rho_{\bullet}$). Solving for the ice
volume fraction of a pebble

$f_{i}=1-\left(\frac{\rho_{\bullet}-\rho_{i}}{\rho_{s}-\rho_{i}}\right)$ (27)

which is the quantity plotted in the bottom row of Fig. 2. We use $\rho_{i}=1$
g cm-3 and $\rho_{s}=3$ g cm-3.

The accretion rates of each model at 20 AU are illustrated in Fig. 3, where
the circles represent the actual accretion rates and they are color-coded
according to the ice volume fraction of the pebbles that are most efficiently
accreted at the given protoplanet mass. Initially the accretion rate is
determined by the geometric cross section and gravitational focusing, with
little contribution from gas drag (green line). When the Bondi radius becomes
larger than the planetesimal, we enter the Bondi regime, where aerodynamic
drag allows efficient capture of pebbles within the Bondi radius (magenta).
Lastly, when the Hill radius becomes larger than the Bondi radius, we enter
the Hill regime, in which the Hill radius becomes the new limiting accretion
radius (orange). In the Bondi regime, the best accreted pebbles are those of
stopping time similar to the time to cross the Bondi radius, which can be
significantly smaller than the biggest pebbles present in the distribution
(Lyra et al., 2023). This is of significant consequence because small bodies
accrete in the Bondi regime (Johansen et al., 2015), opening the possibility
of preferential accretion of small (silicate) pebbles for these bodies.
Indeed, we see in Fig. 3 that the bimodal distribution (Model 2) shows a
window of silicate accretion. Model 2 begins accreting pebbles of $\approx$
50%-50% ices and silicates composition in the geometric regime, but then
accretes nearly 100% silicate pebbles right after the onset of Bondi
accretion. This happens because the transition from the geometric regime to
the Bondi regime is discontinous (Ormel, 2017). In our model, the best
accreted pebble at the geometric regime is of radius $\approx$3 mm, but it
abruptly passes to $\approx$0.5 mm after the onset of Bondi acccretion. The
window of silicate accretion starts to narrow at higher protoplanet masses,
that accrete larger pebbles, of higher ice volume fraction, finally accreting
mostly ices in the Hill regime. Model 1 (power law) never accretes silicates
significantly, as the pure silicate pebbles are too small.

The mass accretion rates are integrated numerically choosing a timestep
$\Delta t$ such that the mass accreted per timestep ($\dot{M}\Delta t$) is no
greater than a fraction $C$ of the planetesimal mass $M_{p}$,

$\Delta t=C\ \frac{M_{p}}{\dot{M}}.$ (28)

We set the value of $C$ at 0.01, found empirically to be a good compromise
between stability and performance. At each time interval, we calculate the new
mass and density of the planetesimal by taking a weighted average of the mass
and density acquired by pebble accretion. Concurrently with pebble accretion,
we reduce the gas and dust density exponentially in time, with an e-folding
time of 2.5 Myr. We terminate the simulation after 4 e-folding times (10 Myr),
which we quote as “the disk lifetime”. We note that at 5-10 Myr there is still
some gas, yet much less than at $t=0$. The pebble accretion rates are impacted
accordingly. To have significant mass accretion rates nearing 10 Myr is
unlikely in our model (because pebbles of all sizes are essentially
decoupled), even though the calculations go until this time.

We highlight and address here the inconsistency of using a two-species pebble
model for streaming instability, while using a continuous size distribution
for the pebble accretion calculation. While inconsistent, this was done
because our goal with the hydroydynamical simulation was to provide a proof-
of-concept that the first planetesimals would be mostly icy, binning the
pebbles in two species, one rocky and one icy, was judged a satisfactory first
order approximation (but see Schaffer et al., 2018; Krapp et al., 2019;
Paardekooper et al., 2020, 2021; Schaffer et al., 2021; Zhu & Yang, 2021; Yang
& Zhu, 2021, for the impact of introducing multiple species in the development
of the streaming instability). For the higher-mass objects, most of the mass
is accreted during the pebble accretion phase.

###  2.6 Porosity Evolution

Finally, we must consider a compaction model to grasp the density evolution of
Kuiper belt objects. We examine here the primary mechanisms by which
planetesimals undergo removal of porosity. The first is compaction through
gravitational pressure. As a planetesimal grows, gravity is continuously
pulling all the material in the planetesimal toward the center. After some
critical mass, around when the central pressure is greater than 10 MPa (Yasui
& Arakawa, 2009; Bierson & Nimmo, 2019), the pull of gravity towards the
center becomes strong enough so that the material yields, and compaction
ensues, removing porosity.

The other mechanism involves the removal of porosity through heating. In the
early solar nebula, there is a non-negligible amount of the short-lived
radioactive isotope, 26Al (initial abundance 26Al/27Al = $5\times{10^{-5}}$
MacPherson et al., 1995; Davidsson et al., 2016). The decay of this isotope,
if present in large quantities in the protoplanet, would act to essentially
melt away porosity from these bodies. Bierson & Nimmo (2019) explore this
effect in detail, along with gravitational compaction. We do not solve the
full system of coupled non-linear partial differential equations; instead, we
use the fact that regardless of the rock-to-ice mass fraction, Bierson & Nimmo
(2019) find that the objects become fully compact at a radius of $R\approx
1500$ km. Thus, we assume the following simplified porosity parametrization

$\phi(R)=\phi_{0}\
\theta(R_{1}-R)+\xi_{R_{1},R_{2}}(R)\,\log_{10}\left(\frac{R_{2}}{R}\right)$
(29)

where $R$ is the radius of the KBO, $R_{2}=1500$ km, and
$R_{1}=R_{2}/10^{\phi_{0}}$. Here $\theta$ is the Heaviside step function

$\theta(x)\equiv\left\\{\begin{array}[]{lr}1,&\text{if }x>0\\\ 1/2,&\text{if
}x=0\\\ 0,&\text{if }x<0\end{array}\right.$ (30)

and $\xi$ is the boxcar function

$\xi_{x_{1},x_{2}}(x)\equiv\theta(x-x_{1})-\theta(x-x_{2}).$ (31)

Figure 4: Left: Graph of Eq. (29), the parametrization we adopt for porosity.
Right: Densities, for constant ice mass fraction $F_{i}$, calculated with Eq.
(29) and Eq. (32) are shown in solid lines. The dashed lines of same color
show the density of a fully compact body of same ice mass fraction.

Eq. (29) is plotted in Fig. 4, left panel. It states that the porosity is
constant at the porosity of formation $\phi_{0}$, up to a radius $R_{1}$,
beyond which the porosity is removed, logarithimically, reaching zero at
$R_{2}$. Choosing $R_{2}=1500$ km and $\phi_{0}=0.5$ results in $R_{1}\approx
474\,$km. Although admittedly crude, this approximation is sufficient for our
purposes. We apply this model concurrently with pebble accretion, using the
newly calculated mass and density to obtain a radius, and then using Eq. (29)
to determine the porosity fraction. To aid the interpretation of our results,
we first isolate the effect the porosity function would have on bodies of
constant composition. The mass fraction of a constituent $j$ (ice or silicate)
is defined as $F_{j}\equiv m_{j}/M_{p}=f_{j}\rho_{j}/\rho$ where $m_{j}$ is
the mass of the constituent in the body, $M_{p}=m_{i}+m_{s}$ is the total
mass, and $f_{j}$ is the constituent’s volume fraction. Given
$f_{v}\equiv\phi$, Eq. (26) can be solved for $\rho$ in terms of the ice mass
fraction $F_{i}$

$\rho=\frac{\left(1-\phi\right)\rho_{i}\rho_{s}}{\rho_{i}+F_{i}\left(\rho_{s}-\rho_{i}\right)}$
(32)

Curves of constant $F_{i}$ calculated according to Eq. (32) are shown in Fig.
4, right panel. The full lines are porous bodies with porosity given by Eq.
(29); dashed horizontal lines mark fully compact (zero porosity) bodies.

##  3 Results

###  3.1 Planetesimal Formation

We run a hydrodynamic simulation where we consider two pebble species, ices
and silicates, with Stokes number $\rm{St}=0.5$ and $\rm{St}=5\times 10^{-3}$,
respectively. We run the simulation until planetesimal formation saturates
(i.e., roughly an orbit goes by without pebbles collapsing into
planetesimals). This condition occurs roughly after five orbits, or 1500
years, considering an orbit at 45 AU is $T=2\pi/\varOmega_{0}\approx 300\text{
years}$.

This is visualized in Fig. 5, where the left panel shows the planetesimal
formation rate as a function of time, and the right panel shows the maximum
grain density within a mesh cell, also as a function of time. The left panel
shows that between the second and fourth orbits, hundreds of planetesimals
were formed, with a peak of 70 per time interval $\varOmega_{0}^{-1}$
($\approx$ 50 yrs). The last planetesimal was formed right after 5 orbits (the
small spike before the rate last goes to zero). The right panel shows that
between the second and fourth orbit, there appears to be at least one mesh
cell with pebble density achieving Roche density, almost continuously. After
about five orbits, the maximum pebble density within a mesh cell steadily
decreases, suggesting that pebble clumps will not reach the Roche density
again. The density within a mesh cell is not allowed to exceed the Roche
density, because once a mesh cell has achieved $\rho_{R}$, the particles
within the accretion radius at formation (set to the mesh spacing $\Delta$)
are removed and replaced by a sink particle.

During the time elapsed, a total of 408 planetesimals are formed, 276 of which
are accreted by other planetesimals, leaving 132 planetesimals by the end of
the simulation.

Figure 5: Left: Planetesimal formation rate (number of planetesimals $N$
formed per time interval $\varOmega_{0}^{-1}$), as a function of time, in the
box. There is a rapid burst of planetesimal formation between two and four
orbits and planetesimal formation comes to a halt after about five orbits.
Right: Time evolution of maximum dust density in a mesh cell. The blue line
corresponds to the Roche density, which if exceeded, particles in that mesh
cell are removed and replaced by a sink particle. The initial rise in density
from zero to one orbit marks the sedimentation of solids and concentration by
the instability, after which the Roche density is achieved for several orbits
but then drops after roughly five orbits.

Figure 6: The vertically integrated (top) and azimuthally integrated (bottom)
dust density after six orbits, when planetesimal formation has saturated. The
left panels correspond to the sum of ice and silicate densities, while the
middle and right panels are the ice and silicate densities. Dots indicate
formed planetestimals and circles in the top panels show the Hill radii of
these planetesimals. The Hill radii appear close to ice overdensities,
consistent with our detailed measurements of high ice mass fractions.

With icy pebbles being less susceptible to the drag force compared to
silicates, they experience lower levels of turbulent diffusion, which provides
support against stellar gravity. The result is that ices form a thinner,
denser mid-plane layer compared to silicates (Dubrulle et al., 1995; Youdin &
Lithwick, 2007) and are therefore more likely to form clumps that collapse
into planetesimals. This is better demonstrated in Fig. 6, where the top row
of panels shows the integrated column density along the z-axis, and the bottom
row the volume density averaged over the y-axis. The first column in each row
shows the combined density of ices and silicates, while in the second and
third columns we disaggregate the pebbles into ices and silicates,
respectively. The circled objects are planetesimals and the size of the
circles represents the Hill radius of the encapsulated planetesimal.

We see in the column density plots that the filamentary structure associated
with the streaming instability is apparent only in the ice plots, whereas the
silicate plot shows a smoother distribution. This is because the silicates are
too tightly coupled to the gas and do not drift as much as the ices, being
thus less prone to the streaming instability (e.g. Yang & Zhu, 2021). The
azimuthal average plots show that the silicates have a height of $H_{d}\approx
0.1$ AU (true bare silicates in stronger turbulence would have a taller scale
height), whereas ices have a denser layer that is more favorable to the
formation of planetesimals.

To help distinguish between these two processes (vertical settling vs
streaming), we analyze the time-evolution of the ice-to-silicate ratio. The
settling time for a grain of a given Stokes number is (Youdin, 2010)

$t_{\rm settle}=\frac{1+2\mathrm{St}^{2}}{\varOmega\mathrm{St}},$ (33)

Figure 7: Evolution of the ice-to-silicate ratio in the vertical plane. The
bright yellow is saturated. The ices settle in about 0.5 orbit, reaching an
ice-to-silicate ratio of about 10 in the midplane. The filamentary structure
of the streaming instability is seen after $\sim$1 orbit.

Figure 8: Left: Ice-to-silicate ratio in the midplane at $t=0.5$ orbits,
approximately the time when the ices settle. Right: Ice-to-silicate ratio in
the midplane at $t=2.0$ orbits, approximately the time of formation of the
first planetesimals. The increase in ice-to-silicate due to settling leads to
a factor 10 enhancement. Further enhancement is due to the streaming
instability.

which yields $\sim$0.3 and 30 orbits, for $\mathrm{St}=0.5$ and
$\mathrm{St}=5\times 10^{-3}$, respectively. We started the simulation with
both dust species at Gaussian stratifications of 0.01H, which is the
equilibrium scale height for the silicate dust grains for $\alpha\sim
10^{-6}$. As a result, the silicate grains do not evolve significantly
vertically. The ice pebbles settle very fast. We show in Fig. 7 the time
evolution of the ice-to-silicate ratio (defined as $\rho_{\rm ice}/\rho_{\rm
sil}$, where $\rho_{\rm ice}$ and $\rho_{\rm sil}$ are the volume densities of
ices and silicates, respectively) in the vertical plane. It mimics the
evolution of the ice pebbles. Indeed we see that the settling time is
$\sim$0.5 orbit, as expected. At 1–2 orbits, the streaming instability
develops and saturates.

We show in Fig. 8 the midplane snapshots of the ice/silicate ratio at $t=0.5$
and 2 orbits. These are the times at $t_{\rm settle}$ for the ices, and when
the first planetesimal form. The left panel shows the ratio resulting from
settling alone. As we see it, it looks homogeneous, with ice-to-silicate ratio
about 10. The right panel shows the filamentary structure expected from
streaming instability, showing the difference in ice-to-silicate ratio between
the filaments and the voids.

Figure 9: Ice mass fraction of the planetesimals formed at the end of the
streaming instability simulations. The composition is not uniform; a trend
emerges where the higher mass planetesimals are formed with higher ice mass
fraction.

Finally, we show in Fig. 9 the mass and ice mass fraction $F_{i}$ of the
planetesimals, at the end of the simulation. The masses range from $\approx
1.5\times{10^{-4}}$ to $\approx 3\times{10^{-3}}$ $M_{\rm Pluto}$, with a
trend that the lower mass planetesimals are more silicate-rich, and the larger
ones more ice-rich. The trend does not seem to be a numerical artifact, a
point we will go back to in Sect. 4.2. We estimate in Appendix B whether this
ice mass fraction would lead to melting from heating from 26Al.

###  3.2 Integrating Pebble Accretion

We take the distribution of planetesimals we found in the previous section,
and feed it into a polydisperse pebble accretion integrator (Lyra et al.,
2023). The streaming instability model was calculated at 45 AU, yet at that
location the pebble accretion rates are too slow, and we do not expect
planetesimals formed at that distance to become protoplanets, given the
existence of the “cold classical” KBOs population (Brown, 2001; Kavelaars et
al., 2008; Petit et al., 2011). In fact, in the context of the “Nice” model
(Tsiganis et al., 2005; Morbidelli et al., 2005; Gomes et al., 2005, see Sect.
4.5) Neptune’s current location at 30 AU constrains that all KBOs aside from
the cold classicals formed up to 30 AU, otherwise Neptune would have continued
its migration further out (Stern et al., 2018, and references therein).
Therefore, we vary the distance at which we perform the pebble accretion
integration, from 10 to 30 AU, in intervals of 5 AU. We also complete the
initial mass function of planetesimals up to $10^{-2}$ $M_{\rm Pluto}\
$(following the composition trend found in the ices and silicates streaming
intability model), given that this is the approximate mass for the onset of
pebble accretion in the Bondi regime (Lyra et al., 2023).

We show in the upper row of Fig. 10 the result for the power law distribution
(model 1), at 20 AU. The panels are a time series, showing the mass vs density
evolution, and color-coded the ice mass fraction of the protoplanets. The
actual KBO data from Fig. 1 are overplotted as green stars. We see that this
model favors ices too strongly to cause any significant changes in density,
even with the removal of porosity through compaction. The silicate pebbles in
this model are simply too small to be accreted efficiently. Pluto mass is
achieved around 4 Myr, but consisting of almost completely pure water ice.

The results change considerably when using the bimodal density distribution
(model 2, middle row of Fig. 10). In the figure, we see that the growth rate
matches the mass-density trend, providing an excellent fit to the
intermediate-mass range objects from $10^{-2}$ (2002 UX25) to $10^{-1}$
$M_{\rm Pluto}\ $(Charon), and also the higher-mass objects, Haumea, Pluto,
and Triton. Yet, the model fails to reproduce the larger density of Eris.

Motivated to reach the density of Eris, we try a 3rd model, which consists of
accretion of pure silicate pebbles. The results are shown in the lower panels
of Fig. 10. Starting from the icy planetesimals produced by the streaming
instability, the density increases monotonically with mass as only silicates
are accreted. This model also provides a good fit to the observed mass-density
trend for intermediate-mass KBOs, up to $10^{-1}$ $M_{\rm Pluto}\ $(evidencing
that model 2 was accreting silicates in this mass range, as we will expand on
in Sect. 4.1). However, beyond this threshold, unlike the bimodal
distribution, the objects have an ice mass fraction that is close to zero. The
model largely overestimates the densities of Haumea, Pluto, and Triton. Yet,
it does reproduce the high density of Eris.

Figure 10: Upper row: Mass, density, and ice mass fraction evolution of
protoplanets with pebble accretion of the power-law density distribution
(model 1, Eqs. (21) and (22)). Actual KBO data is overplotted as green stars.
The best accreted pebbles at each mass are still too icy to lead to
significant density increase through silicate accretion. The most massive
planets formed within 5 Myr are of Pluto mass, but their composition is almost
pure water ice.
Middle row: Same as upper row, but for the bimodal density distribution (model
2, Eqs. (23) and (24)). As the intermmediate-mass objects ($10^{-2}-10^{-1}$
$M_{\rm Pluto}$) accrete silicate pebbles, the mass-density trend is
reproduced. The flattening at higher mass is also reproduced, matching Haumea,
Pluto, and Triton. The model does not reproduce the density of Eris.
Lower row: Same as the other rows, but for constant (silicate) pebble density
(model 3). The intermmediate-mass objects ($10^{-2}-10^{-1}$ $M_{\rm Pluto}$)
is similar to model 2, but the model overshoots the density of the higher mass
objects Haumea, Pluto, and Triton. Eris, however, is reproduced.

####  3.2.1 The effect of distance

We explore now the parameter space of distance in the pebble accretion model,
taking Model 2 as fiducial and varying radii from 10 AU to 30 AU, at 5 AU
intervals. We illustrate the results in Fig. 11. This figure shows that we are
able to best reproduce the mass ranges of KBOs if they formed between 15 AU
and 22 AU. At 10 AU, we see that within 500 thousand years, masses attained
through pebble accretion are a factor of 10 larger than the mass of Pluto.
Considering the lifetime of the solar nebula is of order 10 million years
(i.e., a few e-folding times), it is expected that the masses at this
heliocentric distance will ultimately be several orders of magnitude larger
than the mass of Pluto. On the other hand, considering the accretion rates at
30 AU, we see little to no growth from pebble accretion. After 10 million
years, planetesimals barely grow to $0.2\,M_{\rm{Pluto}}$, suggesting that it
is unlikely that Pluto and the rest of the dwarf planets formed at or beyond
30 AU. A favorable location is between 15 AU and 22 AU, where Pluto’s mass is
reached within 10 million years. We conclude that for planetesimals formed
beyond this range pebble accretion rates would be far too low to achieve Pluto
mass within the assumed lifetime of the solar nebula (as also found by
Lambrechts et al. 2014 for monodisperse pebble accretion). Conversely, if the
planetesimals formed any closer to the Sun, then high accretion rates will
produce planetesimals much larger than Pluto, contrasting the masses seen in
the Kuiper belt.

##  4 Discussion

###  4.1 The Importance of the Silicate Window

To understand these results, we show in Fig. 12 the comparison of the
accretion rates at 10 (diamonds), 20 (stars), and 30 (circles) AU, for the
bimodal distribution. We see that at $\sim 10^{-2}$ $M_{\rm Pluto}$, roughly
the mass where Kuiper belt objects begin to show an increase in density, the
accretion rate at 10 AU is roughly an order of magnitude larger than the
accretion rate at 20 AU, and about two orders of magnitude larger than the
accretion rate at 30 AU. Furthermore, we find that the window of enhanced
silicate accretion provided by the Bondi regime (see Sect. 2.5) shifts to
higher masses with increasing orbital distance. The result is that at 10 AU,
planetesimals with mass $M=10^{-2}M_{\rm{Pluto}}$ have already missed the
window of silicate accretion, resulting in the lower densities seen in the
left-most panel of Fig. 11. Conversely, at 30 AU, the window of silicate
accretion extends beyond $10^{-1}$ $M_{\rm Pluto}$, which could facilitate
this model’s ability to reproduce the density of Eris, but due to the low
accretion rates, these masses are never achieved.

The figure also shows why both model 2 and model 3 reproduce the mass-density
trend for intermmediate masses. This is because at $10^{-2}$ $M_{\rm Pluto}\
$at 20 AU, model 2 is accreting almost pure silicates, like model 3. At
$10^{-1}$ $M_{\rm Pluto}\ $the models diverge significantly as the larger
pebbles are icy in model 2.

Figure 11: The results from applying the polydisperse pebble accretion model
with the bimodal distribution of pebbles across various heliocentric
distances. At 10 AU (upper left) not only do we underestimate the densities of
the KBOs, but the masses produced are much larger than those seen in the
Kuiper belt. At just 500 thousand years, planetesimals are already an order of
magnitude larger than Pluto. From 15 AU to 22 AU, we fit the density trend
well and reach Pluto mass in 1.8 million years, 5.7 million years, and 8.5
million years, respectively. Beyond 22 AU, accretion rates are so low that we
can not reach Pluto’s mass within the assumed lifetime of the disk (10 million
years).

Figure 12: Comparison of accretion rates for the bimodal distribution at
different radii. The window at which silicate accretion is most effective
shifts to higher masses as we move outwards in the disk. This provides further
evidence for 20 AU being a favorable location, as silicate accretion is most
effective around $10^{-2}M_{\rm Pluto}$, right when the densities of Kuiper
belt objects begin to increase.

###  4.2 The planetesimals formed

We perform a streaming instability simulation at 45 AU (a location where no
significant pebble accretion is expected), producing the population of
planetesimals seen in Fig. 9. The span in mass is about an order of magnitude,
from $\approx 1.5\times{10^{-4}}$ to $\approx 3\times{10^{-3}}$ $M_{\rm
Pluto}$. The trend seen in that figure is that smaller planetesimals are more
silicate-rich than the larger ones. This is not a numerical artifact, as Bondi
accretion is too slow at this mass to significantly affect the masses over the
timespan of the hydrodynamical simulation. Also, the Hill radii of these
planetesimals are resolved. For 45 AU and at this resolution, the mass at
which the Hill radius $R_{H}$ is equal to the length $\Delta$ of a mesh cell
is

$\displaystyle M_{p}^{\left(R_{H}=\Delta\right)}$ $\displaystyle=$
$\displaystyle 3\,M_{\sun}\left(\frac{\Delta}{r}\right)^{3}$
$\displaystyle\approx$ $\displaystyle 8.6\times{10^{-5}}\,M_{\rm Pluto}$

or about half the mass of the smallest planetesimal formed. The trend seen in
Fig. 9 is likely physical. We note that this is consistent with the fact that
some small objects are high in density, such as (66652) Borasisi (diameter
$\approx$ 160 km, mean density $\approx$ 2.1 g cm-3 Trujillo et al., 2000;
Noll et al., 2004). The effect of radiogenic heating is a significant function
of ice mass fraction, because 26Al is only present in dust. Naively, we would
expect that the smaller objects would be less affected by radiogenic heating
(bulk heating versus area cooling results in $\propto R$ dependency), but
Davidsson et al. (2016) note that small objects have poor thermal conduction
than larger objects, trapping heat and melting the interior if the ice mass
fraction is low. The melting removes the porosity, making these objects high
density. In our model, significant dispersion in ice mass fraction exists for
streaming instability products less massive than $10^{-3}M_{\rm Pluto}$, so
some objects (the more ice-rich) should avoid melting whereas some (the more
rock-rich) should undergo melting. The scatter for the low mass objects is
consistent with Fig. 1. All streaming instability products more massive than
$10^{-3}M_{\rm Pluto}$ are ice-rich and should avoid melting.

###  4.3 CO ice line

We take this planetesimal population, and feed them into a polydisperse pebble
accretion calculation, at different locations, using three different models
for a pebble size-composition relation, as shown in Fig. 10. We find that the
density and masses of large KBOs are best reproduced if we place their
formation between 15 and 22 AU. Interestingly enough, this distance roughly
coincides with the probable location, at 20 AU, of the CO snowline.

A best location scenario coinciding with the general location of the CO
iceline opens the possibility that CO ice could also be relevant for the
composition of the pebbles and hence, the densities of KBOs. While CO is a
hypervolatile, CO reacts with water to produce methanol (Grundy et al., 2020),
which is refractory. Having a density of $0.64$ g cm-3 at 20 K, methanol ice
could also explain the low densities of the small KBOs, if a significant
fraction of their bulk is of that material. Indeed, near an ice line pebble
growth is boosted through deposition and nucleation of vapor onto bare
silicate and already ice-covered pebbles (see Ros et al., 2019, in the context
of water ice). This possibility is consistent with the lack of observed water
ice spectral features on small KBOs (Barucci et al., 2011; Grundy et al.,
2020).

The CO ice line is also favorable to the formation of the ice giants, not only
owing to the excess of carbon and lack of nitrogen found in the ice giants
(Ali-Dib et al., 2014) but also because ice giants are thought to contain
$>10\%$ methane by mass, which would occur if the ice giants formed between
the CO and N2 ice lines (Dodson-Robinson & Bodenheimer, 2010). Our results are
also in excellent agreement with the constrains that Pluto and the KBOs
(except the cold classicals) formed closer in, as previously discussed in
Sect. 3.2 (see also Malhotra, 1993; Stern et al., 2018; Canup et al., 2021).

###  4.4 No special timing

Our model also has the advantage of not needing to invoke a special formation
time for the Kuiper belt objects, contrasting with the model from Bierson &
Nimmo (2019). This is because silicates, or more specifically the short-lived
radioactive isotope 26Al with a half-life of $7\times 10^{5}$ yr, are not
incorporated in the initial formation of KBOs. Instead, they are incorporated
through Bondi accretion in the later stages of the evolution of KBOs, when a
large fraction of 26Al has already been depleted. Lastly, since we were able
to reproduce the high density of Eris through the pebble distribution model of
pure silicates, our model could suggest that Eris formed from a streaming
instability filament (Lorek & Johansen, 2022) with a large fraction of rock,
or that Eris formed at a later time when volatiles in the disk were lost
through drift and photoevaporation (although that would potentially not leave
enouh time for pebble accretion to operate). Alternatively, Eris could have
formed with the density predicted by model 2, and subsequently lost some of
its ice mantle through collisional evolution (Barr & Schwamb, 2016).

###  4.5 Connection to the Nice Model

Our model is in agreement with the “Nice” model scenario (Tsiganis et al.,
2005; Morbidelli et al., 2005; Gomes et al., 2005; Emel’yanenko, 2022). In
this picture, after the dissipation of the solar nebula, the giant planets are
initially in a more compact configuration than currently, and their masses
decreasing monotonically as heliocentric distance increases (i.e., Uranus and
Neptune swapped). Jupiter is placed in the vicinity of its current location at
5 AU, Uranus around 11 - 17 AU, and Saturn and Neptune between these two
limits. The giant planets are also assumed to have been in near circular and
co-planar orbits.

A belt of protoplanets exists just beyond the orbit of the outermost giant
planet (in this case Uranus), and objects in the inner edge of the belt
interact with the planet, being scattered inward or outward, exchanging
angular momentum (Fernandez & Ip, 1984). The net gain of angular momentum of
Uranus would be zero in this scenario, but the presence of Jupiter breaks the
symmetry. A protoplanet scattered outward will likely return to interact with
the planet again, but a protoplanet scattered inward has a high chance of
interacting with Jupiter and getting ejected out of the solar system. The
inner disk is thus a better sink of angular momentum than the outer disk, and
the net result is that Uranus migrates outward, whereas Jupiter migrates
inward. Adding Neptune and Saturn does not alter the conclusion; these planets
also migrate outward as they scatter protoplanets toward Jupiter, than in turn
ejects them. Once Jupiter and Saturn cross their 2:1 mean orbital resonance,
dynamical instability ensues, throwing the ice giants into the primordial
belt, and implanting a small fraction (0.01%-0.1%) of the objects into the
present-day Kuiper belt (except the cold-classicals, that likely formed in
situ). Further interaction with the belt dynamically cooled the orbits of the
giant planets post-instability. The original belt extended at most up to 30
AU, the current orbit of Neptune (otherwise Neptune would have migrated
further outward).

Our results are in excellent agreement with this model. While planetesimals
can form by streaming instability at large heliocentric distances, explaining
the cold classicals, pebble accretion is only efficient up to $\approx$25 AU
(Johansen et al., 2015). Beyond this distances, planetesimals do not grow up
to Triton or Pluto sizes, even by polydisperse pebble accretion Lyra et al.
(2023), as we have demonstrated. Scattering small planetesimals is not
efficient to drive further migration, so Neptune’s current location coincides
with where pebble accretion stops being efficient at growing large objects.

##  5 Limitations and Future Work

The model presented is an exploratory proof-of-concept model, and as such has
a number of simplications.

First and foremost, we argued in the introduction that compositional
differences in the pebbles should be expected, because of photodesorption of
ices off small dust grains. The compositional difference can only be
maintained if the photodesorption rate and coagulation rate are similar, which
a growth time estimate indeed shows they are. The growth time of dust grains,
up to a factor of order unity, is $t_{\rm grow}\sim\varepsilon^{-1}T_{\rm
orb}$ (Birnstiel et al., 2012; Lambrechts & Johansen, 2014; Lorek & Johansen,
2022), where $T_{\rm orb}$ is the orbital period. The disk starts with
$\varepsilon\sim 10^{-2}$ of dust, of interstellar (sub-micron) size. Assuming
an orbital period $T_{\rm orb}\sim 100$ years yields $t_{\rm grow}\sim 10^{4}$
years. The dust thus coagulates very quickly into pebbles. The steady state,
however, is top-heavy, leaving about $\varepsilon\sim 10^{-4}$ of sub-micron
dust, most of the mass residing in larger pebbles. Under these conditions,
after most small grains are consumed, the growth time then jumps to $\sim$1
Myr, which is similar to the photodesorption rate.

Despite this assurance, the argument remains mostly qualitative, and thus the
conversion between radius and composition we used is arbitrary. A future model
should include coagulation, drift, fragmentation, turbulence, UV and cosmic
ray desorption, calculating the grain size vs composition relation from first
principles. Such a model would also have the benefit of constraining the
transition in the bimodal model, which we arbitrarily placed at 1 mm.

Secondly, it is an inconsistency that we used a two-species pebble model for
the streaming instability, whereas the pebble accretion integrator used a
continuum distribution. Ideally, the two calculations should be consistent,
with the streaming instabiliy model also using a continuum of pebble sizes, of
mixed composition. The computational cost of three-dimensional hydrodynamical
streaming instability simulations also motivated this simplification.

Since we found that the best location for formation is near the CO ice line,
our pebble accretion model could benefit from evolving the CO abundance,
matching chemical evolutionary expectations, and the depletion of volatiles
from heating, radial drift, and disk winds.

It is also an idealization that our model assumes an infinite supply of
pebbles whereas, realistically, the pebbles drift and are eventually lost if
there is no re-supply. Interestingly, this would perhaps lead to a
compositional variation in time, because the larger pebble drift faster than
the smaller ones. At earlier times, when planetesimals first form, ice allows
for the formation of large pebbles that drift, rapidly depleting the ice.
Pebble accretion, conversely, is a slower process, so it would mostly occur
after ice depletion in the outer disk, feeding from the remaining smaller,
more silicate-rich, dust grains. Pebble drift would be specially important for
the evolution of CO-coated pebbles, because once they drift inward to the CO
ice line, the CO gas evaporates and is vulnerable to disk winds,
photoevaporation, or photodissociation.

Our model would also benefit from a more involved porosity evolution model
that takes into consideration the abundance of radioactive elements like 26Al
and the impact that radioactive heating has on the porosity of KBOs.

##  6 Conclusions

We set out to reproduce the density trend of Kuiper belt objects, where small
objects have densities less than water ice and large Kuiper belt objects can
reach densities of $2.5\text{ g cm}^{-3}$. We run a hydrodynamic simulation to
model planetesimal formation via the streaming instability of ices and
silicates, where silicates are small grains with short friction times, and
ices are large grains with long friction times. We find that planetesimals
formed under these conditions are icy and are low in mass $M_{p}\lesssim
10^{-2}M_{\text{Pluto}}$, effectively reproducing the densities observed in
the low mass Kuiper belt objects. We also find a correlation between ice mass
fraction and planetesimal mass, albeit with significant dispersion.

We then use these planetesimals formed by the streaming instability as
starting masses for the recently derived polydisperse pebble accretion model
(Lyra et al., 2023). We consider three different population of pebbles: the
first is a power law distribution where the density of pebbles follow the
power law in Eqs. (21) and (22), with the largest grain being pure ice and the
smallest grain being pure silicate. The second is a bimodal distribution,
given by Eqs. (23) and (24), where grains smaller than $1$ mm are pure
silicates, and larger grains follow a power law, with the largest grain being
pure ice. The last model assumes all pebbles are silicates.

We find the power law model does not reproduce the densities for high-mass
KBOs. The model with purely silicate pebbles is able to reproduce the high
density of Eris, yet fails to reproduce the densities of Haumea, Pluto, and
Triton. The bimodal distribution, however, is capable of reproducing the
densities of the dwarf planets but underestimates the density of Eris. Thus we
speculate that Eris could have formed from a rock-rich filament or that it
formed later in the solar system history when volatiles were lost through
radial drift or photoevaporation. Nevertheless, it is conceivable that Eris
formed with the density predicted our bimodal distribution, and subsequently
lost ice through collisions, which we do not model.

We find that the masses of KBOs are best reproduced between 15 and 22 AU.
Beyond this range, accretion rates are far too low to achieve dwarf-planet
mass by the end of the disk lifetime. Conversely, inwards of 15 AU, accretion
rates are too high, resulting in masses that are orders of magnitudes larger
than Triton and Pluto.

Our solution avoids the timing problem that KBOs formed too early would melt
and become compact, due to the energy released by 26Al. In our model, the
first planetesimals are icy, and 26Al is mostly incorporated in the long phase
of silicate pebble accretion, when most 26Al has already decayed. While the
specific results on location of formation and final mass are dependent on the
disk model adopted, the premise and conclusions of the work, namely the need
to separate the silicates from ices and then preferentially accrete silicates,
and that Bondi accretion and ice desorption from small grains are a way to
accomplish these, are independent of the particular disk model.

Our results lend further credibility to the streaming instability as a
planetesimal formation mechanism and to pebble accretion as a mechanism by
which planetesimals grow into larger bodies. This model also provides yet
another challenge for the previously held idea of planetesimal accretion as
the main driver for growth. Growth through binary planetesimal accretion would
result in planetesimals with similar densities regardless of size, with a
maximum increase of a factor of two due to gravitational compaction. We show
that multi-species streaming instability could result in planetesimals of
nearly uniform composition and that a polydisperse pebble accretion model can
have a significant impact on the final composition of a planetary body.

MHC, WL, DS, JBS, OMU, CCY, and ANY acknowledge support from the NASA
Theoretical and Computational Astrophysical Networks (TCAN) via grant
80NSSC21K0497. MHC and WL are further supported by grant #80NSSC22K1419 from
the NASA Emerging Worlds program, DS by NSF via grant AST-2007422, and CCY
acknowledges the support from NASA via the Astrophysics Theory Program (grant
#80NSSC21K0141) and the Emerging Worlds program (grant #80NSSC20K0347 and
#80NSSC23K0653). We acknowledge suggestions and constructive criticism from
Anders Johansen, Zsolt Sándor, Octavio Guilera, and Paul Estrada. The
simulations presented in this paper utilized the Stampede cluster of the Texas
Advanced Computing Center (TACC) at The University of Texas at Austin, through
XSEDE/Access grant TG-AST140014, and the Discovery cluster at New Mexico State
University (Trecakov & Von Wolff, 2021). This work utilized resources from the
New Mexico State University High Performance Computing Group, which is
directly supported by the National Science Foundation (OAC-2019000), the
Student Technology Advisory Committee, and New Mexico State University and
benefits from inclusion in various grants (DoD ARO-W911NF1810454; NSF EPSCoR
OIA-1757207; Partnership for the Advancement of Cancer Research, supported in
part by NCI grants U54 CA132383 (NMSU)). This research was made possible by
the open-source projects jupyter (Kluyver et al., 2016), IPython (Perez &
Granger, 2007), matplotlib (Hunter, 2007; Caswell et al., 2020), NumPy (Harris
et al., 2020), SymPy (Meurer et al., 2017), and AstroPy (Astropy Collaboration
et al., 2013).

##  Appendix A Drag force backreaction with multiple species

Eqs. (3) and (4) are solved for every particle $k$ in the simulation. The drag
force accelerating particle $k$ is

${\boldsymbol{f}}_{d,k}=-\frac{{\boldsymbol{v}}_{k}-\overline{{\boldsymbol{u}}}({\boldsymbol{x}}_{k})}{\tau_{k}}$
(A1)

where $\tau$, the friction time, is the timescale at which the relative
velocity between the gas and particle $k$ is reduced by a factor of $e$. For
grains smaller than the mean free path of the gas (Epstein drag, Epstein,
1924), the appropriate regime for our investigation (e.g. Johansen et al.,
2014), this timescale is related to the grain size by

$\tau_{k}=\frac{a_{k}\ \rho_{\bullet,k}}{v_{\rm
th}\overline{\rho_{g}}({\boldsymbol{x}}_{k})}.$ (A2)

Here $a$ is the grain radius, $\rho_{\bullet}$ is the grain internal density
and $v_{\rm th}=\sqrt{8/\pi}c_{s}$ is the thermal speed of the molecules. The
friction time is often non-dimensionalized by defining the Stokes number
$\mathrm{St}$,

$\mathrm{St}_{k}\equiv\varOmega_{0}\tau_{k}.$ (A3)

The quantities $\overline{{\boldsymbol{u}}}({\boldsymbol{x}}_{k})$ and
$\overline{\rho_{g}}({\boldsymbol{x}}_{k})$ refer to the gas velocity
${\boldsymbol{u}}({\boldsymbol{x}})$ and density $\rho_{g}({\boldsymbol{x}})$,
defined on the Eulerian mesh positions ${\boldsymbol{x}}$, interpolated to the
position ${\boldsymbol{x}}_{k}$ of particle $k$. The interpolation for a
quantity $\psi$ is done as described in Youdin & Johansen (2007)

$\overline{\psi}({\boldsymbol{x}}_{k})=\sum
W\left({\boldsymbol{x}}_{k}-{\boldsymbol{x}}\right)\psi({\boldsymbol{x}})$
(A4)

where the sum is done over mesh positions. The mesh weight kernel $W$ is
chosen as the Triangular Shaped Cloud algorithm (Hockney & Eastwood, 1988),
where the nearest 3 cells to the particle in each direction (27 cells in
total) contribute to the interpolation. Incidentally, the pebble density
needed for Eq. (5) is calculated on the mesh cells according to

$\rho_{d}({\boldsymbol{x}})=\frac{1}{V({\boldsymbol{x}})}\sum_{k}W({\boldsymbol{x}}_{k}-{\boldsymbol{x}})m_{k},$
(A5)

and once the selfgravity potential $\varPhi({\boldsymbol{x}})$ is calculated
and the gradient ${\boldsymbol{\nabla}}{\varPhi}({\boldsymbol{x}})$ taken on
the mesh, the interpolation
$\overline{{\boldsymbol{\nabla}}{\varPhi}}({\boldsymbol{x}}_{k})$ to the
position of each particle $k$ is done according to Eq. (A4) and added to Eq.
(4).

The drag force backreaction ${\boldsymbol{f}}_{g}$ onto a mesh cell centered
at ${\boldsymbol{x}}$ is calculated in a momentum-conserving way, via Newton’s
3rd law (Youdin & Johansen, 2007)

${\boldsymbol{f}}_{g}({\boldsymbol{x}})=-\frac{1}{\rho_{g}({\boldsymbol{x}})V({\boldsymbol{x}})}\sum_{k}W\left({\boldsymbol{x}}_{k}-{\boldsymbol{x}}\right)m_{k}\,f_{d,k}$
(A6)

where $V({\boldsymbol{x}})$ is the mesh cell volume and $m_{k}$ the mass of
particle $k$. In this formulation, the dragforce backreaction for multiple
particle species is straightforward, as the sum collects the individual
contribution of each particle $k$.

##  Appendix B Heating from decay of 26Al

Figure 13: Body temperature $T_{1}$ resulting from heating from decay of 26Al,
as a function of mass fraction of silicates, $F_{s}$.

We estimate here whether the small bodies produced in the streaming
instability model would melt as a result of heating from 26Al. Considering the
model from Castillo-Rogez et al. (2009), the volumetric heating rate
$\mathcal{H}$ is

$\mathcal{H}=\rho\,F_{s}\,{\rm[^{26}Al]}_{0}\,\mathcal{H}_{0}\,e^{-\lambda t}$
(B1)

Where $\rho$ is the protoplanet’s density, $F_{s}$ is the silicate mass
fraction, ${\rm[^{26}Al]}_{0}$ is the initial isotopic abundance of 26Al in
ordinary chondrites, $\mathcal{H}_{0}$ is the heat production rate per mass,
and $\lambda=\ln(2)/t_{1/2}$ is the decay constant of 26Al; $t_{1/2}$ is the
half-life. We find the energy by integrating in time and multiplying by the
volume $V$

$\displaystyle Q(t)$ $\displaystyle=$ $\displaystyle
V\int_{0}^{t}\mathcal{H}(t^{\prime})dt^{\prime}$ (B2) $\displaystyle=$
$\displaystyle
M_{p}\,F_{s}\,{\rm[^{26}Al]}_{0}\,\mathcal{H}_{0}\,\lambda^{-1}\,\left(1-e^{-\lambda
t}\right)$ (B3)

where we substituted the protoplanet’s mass $M_{p}=\rho V$. If we assume that
all this energy is trapped within the body, we can find the temperature
$T_{1}$ achieved by setting $t\rightarrow\infty$, and equating with the heat
capacity equation

$Q=M_{p}\,c_{p}\,\Delta T$ (B4)

where $c_{p}$ is the heat capacity at constant pressure. Then solving for
$\Delta T$

$\Delta
T=F_{s}\,{\rm[^{26}Al]}_{0}\,\mathcal{H}_{0}\,\lambda^{-1}\,c_{p}^{-1}$ (B5)

Substituting the values $c_{p}=2108\ {\rm J\,kg^{-1}\,K^{-1}}$ (Marquet,
2015), ${\rm[^{26}Al]}_{0}=\left(8.6\times 10^{-8},5.37\times 10^{-7}\right)$
(Parhi & Prialnik, 2023), $\mathcal{H}_{0}=0.355\ {\rm W/kg}$ (Castillo-Rogez
et al., 2009), $t_{1/2}=0.717\ {\rm Myr}$, and $T_{0}=30\,{\rm K}$, we find
the $F_{s}$ vs $T_{1}$ (the final temperature) relationship shown in Fig. 13.
The blue curve is for the lower value of ${\rm[^{26}Al]}_{0}$, the red curve
for the higher value. The light grey dotted line is for heating up to 250 K;
the dark grey dashed line for heating up to 200 K. It seems that according to
this model, with no energy release, a temperature of 200 K (250 K) is achieved
at a mere $F_{s}$=0.06 (0.08) for the higher 26Al abundance. The lower
abundance of 26Al reaches 200 K for $F_{s}$=0.36, or 250 K at $F_{s}$=0.47.

In this model, the higher abundance would not allow for porosity retention,
but the lower abundance accomodates the $\gtrsim 65\%$ ice mass fraction seen
in Fig. 9 for the products of streaming instability. We highlight that these
are likely overestimates for the heating rate because, as noted by Parhi &
Prialnik (2023), part of the heat is used to sublimate the hypervolatiles CO
and CH4, which considerably lowers the attained final temperature. We note
that this is consistent with the very low densities of the small KBOs, without
necessitating unrealistic porosities.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192
  * Ali-Dib et al. (2014) Ali-Dib, M., Mousis, O., Petit, J.-M., & Lunine, J. I. 2014, ApJ, 785, 125
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3
  * Andersson et al. (2006) Andersson, S., Al-Halabi, A., Kroes, G.-J., & van Dishoeck, E. F. 2006, J. Chem. Phys., 124, 064715
  * Andersson & van Dishoeck (2008) Andersson, S., & van Dishoeck, E. F. 2008, A&A, 491, 907
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Baer et al. (2011) Baer, J., Chesley, S. R., & Matson, R. D. 2011, AJ, 141, 143
  * Bai & Stone (2010) Bai, X.-N., & Stone, J. M. 2010, ApJL, 722, L220
  * Bai & Stone (2011) —. 2011, ApJ, 736, 144
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
  * Barr & Schwamb (2016) Barr, A. C., & Schwamb, M. E. 2016, MNRAS, 460, 1542
  * Barucci et al. (2011) Barucci, M. A., Alvarez-Candal, A., Merlin, F., et al. 2011, Icarus, 214, 297
  * Bergin et al. (1995) Bergin, E. A., Langer, W. D., & Goldsmith, P. F. 1995, ApJ, 441, 222
  * Bierson & Nimmo (2019) Bierson, C. J., & Nimmo, F. 2019, Icarus, 326, 10
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  * Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
  * Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
  * Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859
  * Brown (2001) Brown, M. E. 2001, AJ, 121, 2804
  * Brown (2012) —. 2012, Annual Review of Earth and Planetary Sciences, 40, 467
  * Brown (2013) —. 2013, ApJL, 778, L34
  * Brown & Butler (2017) Brown, M. E., & Butler, B. J. 2017, AJ, 154, 19
  * Canup et al. (2021) Canup, R. M., Kratter, K. M., & Neveu, M. 2021, in The Pluto System After New Horizons, ed. S. A. Stern, J. M. Moore, W. M. Grundy, L. A. Young, & R. P. Binzel, 475–506
  * Carrasco-González et al. (2019) Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Castillo-Rogez et al. (2009) Castillo-Rogez, J., Johnson, T. V., Lee, M. H., et al. 2009, Icarus, 204, 658
  * Caswell et al. (2020) Caswell, T. A., Droettboom, M., Lee, A., et al. 2020, matplotlib/matplotlib: REL: v3.3.1, doi:10.5281/zenodo.3984190
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
  * Ciesla (2014) Ciesla, F. J. 2014, ApJL, 784, L1
  * Davidsson et al. (2016) Davidsson, B. J. R., Sierks, H., Güttler, C., et al. 2016, A&A, 592, A63
  * Dodson-Robinson & Bodenheimer (2010) Dodson-Robinson, S. E., & Bodenheimer, P. 2010, Icarus, 207, 491
  * Drażkowska & Dullemond (2014) Drażkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Dullemond & Dominik (2005) Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971
  * Emel’yanenko (2022) Emel’yanenko, V. V. 2022, A&A, 662, L4
  * Epstein (1924) Epstein, P. S. 1924, Physical Review, 23, 710
  * Fernandez & Ip (1984) Fernandez, J. A., & Ip, W. H. 1984, Icarus, 58, 109
  * Geretshauser et al. (2010) Geretshauser, R. J., Speith, R., Güttler, C., Krause, M., & Blum, J. 2010, A&A, 513, A58
  * Gomes et al. (2005) Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466
  * Grigorieva et al. (2007) Grigorieva, A., Thébault, P., Artymowicz, P., & Brandeker, A. 2007, A&A, 475, 755
  * Grundy et al. (2019) Grundy, W. M., Noll, K. S., Buie, M. W., et al. 2019, Icarus, 334, 30
  * Grundy et al. (2007) Grundy, W. M., Stansberry, J. A., Noll, K. S., et al. 2007, Icarus, 191, 286
  * Grundy et al. (2020) Grundy, W. M., Bird, M. K., Britt, D. T., et al. 2020, Science, 367, aay3705
  * Gundlach & Blum (2015) Gundlach, B., & Blum, J. 2015, ApJ, 798, 34
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  * Güttler et al. (2009) Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Harrison & Schoen (1967) Harrison, H., & Schoen, R. I. 1967, Science, 157, 1175
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  * Hockney & Eastwood (1988) Hockney, R. W., & Eastwood, J. W. 1988, Computer simulation using particles
  * Hogerheijde et al. (2011) Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Johansen et al. (2014) Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547
  * Johansen & Klahr (2005) Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353
  * Johansen & Lacerda (2010) Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475
  * Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109
  * 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. 2007, ApJ, 662, 627
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJL, 704, L75
  * Kavelaars et al. (2008) Kavelaars, J., Jones, L., Gladman, B., Parker, J. W., & Petit, J. M. 2008, in The Solar System Beyond Neptune, M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, and A. Morbidelli (eds.), University of Arizona Press, Tucson, 592 pp., p.59-69, ed. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson, 59
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Scmidt (IOS Press), 87–90
  * Krapp et al. (2019) Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJL, 878, L30
  * Krijt et al. (2016) Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, 833, 285
  * Krijt et al. (2015) Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83
  * Krijt et al. (2018) Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78
  * Lambrechts & Johansen (2012) Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32
  * Lambrechts & Johansen (2014) —. 2014, A&A, 572, A107
  * Lambrechts et al. (2014) Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35
  * Lesur et al. (2022) Lesur, G., Ercolano, B., Flock, M., et al. 2022, arXiv e-prints, arXiv:2203.09821
  * Li & Youdin (2021) Li, R., & Youdin, A. N. 2021, ApJ, 919, 107
  * Li et al. (2019) Li, R., Youdin, A. N., & Simon, J. B. 2019, ApJ, 885, 69
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17
  * Lisse et al. (2021) Lisse, C. M., Young, L. A., Cruikshank, D. P., et al. 2021, Icarus, 356, 114072
  * Lorek & Johansen (2022) Lorek, S., & Johansen, A. 2022, arXiv e-prints, arXiv:2208.01902
  * Lyra et al. (2023) Lyra, W., Johansen, A., Cañas, M. H., & Yang, C.-C. 2023, arXiv e-prints, arXiv:2301.03825
  * Lyra et al. (2008a) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008a, A&A, 491, L41
  * Lyra et al. (2008b) —. 2008b, A&A, 479, 883
  * Lyra et al. (2009a) —. 2009a, A&A, 493, 1125
  * Lyra et al. (2009b) Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009b, A&A, 497, 869
  * Lyra & Lin (2013) Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17
  * Lyra et al. (2021) Lyra, W., Youdin, A. N., & Johansen, A. 2021, Icarus, 356, 113831
  * MacPherson et al. (1995) MacPherson, G. J., Davis, A. M., & Zinner, E. K. 1995, Meteoritics, 30, 365
  * Malhotra (1993) Malhotra, R. 1993, Nature, 365, 819
  * Manara et al. (2018) Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3
  * Marquet (2015) Marquet, P. 2015, Quarterly Journal of the Royal Meteorological Society, 141, 67
  * McKinnon et al. (2020) McKinnon, W. B., Richardson, D. C., Marohnic, J. C., et al. 2020, Science, 367, aay6620
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103
  * Morbidelli et al. (2005) Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Nature, 435, 462
  * Morgado et al. (2023) Morgado, B. E., Sicardy, B., Braga-Ribas, F., et al. 2023, Nature, 614, 239
  * Musiolik (2021) Musiolik, G. 2021, MNRAS, 506, 5153
  * Musiolik & Wurm (2019) Musiolik, G., & Wurm, G. 2019, ApJ, 873, 58
  * Nakagawa et al. (1981) Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517
  * Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Nesvorný et al. (2019) Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nature Astronomy, 349
  * Nesvorný et al. (2010) Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785
  * Noll et al. (2004) Noll, K. S., Stephens, D. C., Grundy, W. M., & Griffin, I. 2004, Icarus, 172, 402
  * Norris et al. (1983) Norris, T. L., Gancarz, A. J., Rokop, D. J., & Thomas, K. W. 1983, Lunar and Planetary Science Conference Proceedings, 88, B331
  * Öberg et al. (2009) Öberg, K. I., Linnartz, H., Visser, R., & van Dishoeck, E. F. 2009, ApJ, 693, 1209
  * Ormel (2017) Ormel, C. W. 2017, in Astrophysics and Space Science Library, Vol. 445, Astrophysics and Space Science Library, ed. M. Pessah & O. Gressel, 197
  * Ormel & Klahr (2010) Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43
  * Ortiz et al. (2017) Ortiz, J. L., Santos-Sanz, P., Sicardy, B., et al. 2017, Nature, 550, 219
  * Paardekooper et al. (2020) Paardekooper, S.-J., McNally, C. P., & Lovascio, F. 2020, MNRAS, 499, 4223
  * Paardekooper et al. (2021) —. 2021, MNRAS, 502, 1579
  * Parhi & Prialnik (2023) Parhi, A., & Prialnik, D. 2023, MNRAS, 522, 2081
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration, Brandenburg, A., Johansen, A., et al. 2021, The Journal of Open Source Software, 6, 2807
  * Pereira et al. (2023) Pereira, C. L., Sicardy, B., Morgado, B. E., et al. 2023, A&A, 673, L4
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science Engineering, 9, 21
  * Pérez et al. (2015) Pérez, L. M., Chandler, C. J., Isella, A., et al. 2015, ApJ, 813, 41
  * Petit et al. (2011) Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131
  * Planes et al. (2020) Planes, M. B., Millán, E. N., Urbassek, H. M., & Bringa, E. M. 2020, MNRAS, 492, 1937
  * Planes et al. (2021) —. 2021, MNRAS, 503, 1717
  * Potapov et al. (2019) Potapov, A., Jäger, C., & Henning, T. 2019, ApJ, 880, 12
  * Powell et al. (2022) Powell, D., Gao, P., Murray-Clay, R., & Zhang, X. 2022, Nature Astronomy, arXiv:2208.13806
  * Rawlings (2022) Rawlings, J. M. C. 2022, MNRAS, 517, 3804
  * Ribas et al. (2014) Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54
  * Ros et al. (2019) Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65
  * Safronov (1960) Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979
  * Safronov (1972) —. 1972, Evolution of the protoplanetary cloud and formation of the earth and planets. (Jerusalem (Israel): Israel Program for Scientific Translations, Keter Publishing House)
  * Sai et al. (2023) Sai, J., Yen, H.-W., Ohashi, N., et al. 2023, arXiv e-prints, arXiv:2307.08952
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69
  * Schaffer et al. (2021) Schaffer, N., Johansen, A., & Lambrechts, M. 2021, A&A, 653, A14
  * Schaffer et al. (2018) Schaffer, N., Yang, C.-C., & Johansen, A. 2018, A&A, 618, A75
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Silsbee et al. (2021) Silsbee, K., Caselli, P., & Ivlev, A. V. 2021, MNRAS, 507, 6205
  * Simon et al. (2017) Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJL, 847, L12
  * Simon et al. (2018) Simon, J. B., Bai, X.-N., Flaherty, K. M., & Hughes, A. M. 2018, ApJ, 865, 10
  * Simon et al. (2016) Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169
  * Sipilä et al. (2021) Sipilä, O., Silsbee, K., & Caselli, P. 2021, ApJ, 922, 126
  * Squire & Hopkins (2020) Squire, J., & Hopkins, P. F. 2020, MNRAS, 498, 1239
  * Stansberry et al. (2006) Stansberry, J. A., Grundy, W. M., Margot, J. L., et al. 2006, ApJ, 643, 556
  * Stansberry et al. (2012) Stansberry, J. A., Grundy, W. M., Mueller, M., et al. 2012, Icarus, 219, 676
  * Stern et al. (2018) Stern, S. A., Grundy, W. M., McKinnon, W. B., Weaver, H. A., & Young, L. A. 2018, ARA&A, 56, 357
  * Stuber & Wolf (2022) Stuber, T. A., & Wolf, S. 2022, A&A, 658, A121
  * Tazzari et al. (2016) Tazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130
  * Tominaga et al. (2021) Tominaga, R. T., Inutsuka, S.-i., & Kobayashi, H. 2021, ApJ, 923, 34
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
  * Trecakov & Von Wolff (2021) Trecakov, S., & Von Wolff, N. 2021, in Practice and Experience in Advanced Research Computing, PEARC ’21 (New York, NY, USA: Association for Computing Machinery)
  * Trujillo et al. (2000) Trujillo, C. A., Luu, J. X., Jewitt, D. C., et al. 2000, Minor Planet Electronic Circulars, 2000-O12
  * Tsiganis et al. (2005) Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459
  * Wada et al. (2009) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
  * Westley et al. (1995a) Westley, M. S., Baragiola, R. A., Johnson, R. E., & Baratta, G. A. 1995a, Nature, 373, 405
  * Westley et al. (1995b) —. 1995b, Planet. Space Sci., 43, 1311
  * Yamato et al. (2023) Yamato, Y., Aikawa, Y., Ohashi, N., et al. 2023, ApJ, 951, 11
  * Yang & Johansen (2014) Yang, C.-C., & Johansen, A. 2014, ApJ, 792, 86
  * Yang et al. (2017) Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Yang & Zhu (2021) Yang, C.-C., & Zhu, Z. 2021, MNRAS, 508, 5538
  * Yasui & Arakawa (2009) Yasui, M., & Arakawa, M. 2009, Journal of Geophysical Research (Planets), 114, E09004
  * Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, ApJ, 662, 613
  * Youdin (2010) Youdin, A. N. 2010, in EAS Publications Series, Vol. 41, EAS Publications Series, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange, 187–207
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Youdin & Shu (2002) Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494
  * Zhu & Yang (2021) Zhu, Z., & Yang, C.-C. 2021, MNRAS, 501, 467
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57