This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

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 μ\mum-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-μ\mum 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).

Refer to caption
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 (≲2×10−2\lesssim 2\times 10^{-2} Pluto masses) have a near-constant density of 0.50.5 g cm-3 and larger KBOs have an increasingly larger density that reaches ≈2.5​ g cm−3\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 (1.20±0.05)×1024\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 ≈\approx1 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-μ\mum size (Krijt et al., 2018). Yet, the optical depth at a given height decreases with time as the sub-μ\mum 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-μ\mum grain species decreases significantly. When the dust-to-gas ratio of these sub-μ\mum grain species decreases to 10−410^{-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 Dr\każ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​ g cm−32.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
tt time r0r_{0} arbitrary origin of shearing box
x,y,zx,y,z Cartesian coordinates pp Eq. (8) gas pressure
ε\varepsilon ρd/ρg\rho_{d}/\rho_{g} dust-to-gas density ratio csc_{s} Eq. (11) the sound speed
ZZ ∫ρd​𝑑V∫ρg​𝑑V\frac{\int\rho_{d}dV}{\int\rho_{g}dV} dust-to-gas mass ratio d​VdV d​x​d​y​d​zdx\,dy\,dz volume infinitesimal
vKv_{K} Ω0​r0\varOmega_{0}r_{0} Keplerian velocity Δ​v\Delta v η​vK\eta v_{K} gas velocity reduction
γ\gamma adiabatic index τ\tau Eq. (A2) particle friction time
vthv_{\rm th} 8π​cs\sqrt{\frac{8}{\pi}}c_{s} thermal velocity ρ∙\rho_{\bullet} pebble material density
kBk_{B} Boltzmann constant 𝒗{\boldsymbol{v}} (vx,vy,vz)(v_{x},v_{y},v_{z}) pebble velocity
TT temperature λSI\lambda_{\rm SI} η​r\eta r streaming instability length scale
mHm_{H} atomic mass unit HdH_{d} Eq. (13) pebble scale height
ΩK\varOmega_{K} Eq. (12) Keplerian frequency NN Nx,Ny,NzN_{x},N_{y},N_{z} box resolution
MM stellar mass G~\widetilde{G} Eq. (18) dimensionless gravitational constant
δ\delta dimensionless pebble diffusion ρd\rho_{d} volume density of pebbles
Δ\Delta mesh spacing aa pebble radius
ϕ\phi porosity St\rm{St} Ω0​τ\varOmega_{0}\tau Stokes number of pebbles
RR protoplanet radius Φ\varPhi self-gravitational potential
ρs\rho_{s} silicate density μ\mu mean molecular weight
ρi\rho_{i} ice density LL LxL_{x}=LyL_{y}=LzL_{z} length of shearing box sides
HH cs/Ω0c_{s}/\varOmega_{0} gas scale height α\alpha Eq. (14) dimensionless gas viscosity
𝒖{\boldsymbol{u}} (ux,uy,uz)(u_{x},u_{y},u_{z}) gas velocity in Cartesian coordinates NpN_{p} number of numerical particles
Ω0\varOmega_{0} Keplerian frequency at r0r_{0} QQ Eq. (19) Toomre QQ parameter
ρg\rho_{g} gas density ρR\rho_{R} Eq. (20) Roche density
hh H/r0H/r_{0} disk aspect ratio M˙\dot{M} mass accretion rate
η\eta Eq. (6) dimensionless pressure support rr astrocentric distance
u0u_{0} −32​Ω0​x-\frac{3}{2}\varOmega_{0}x Keplerian shear flow MpM_{p} protoplanet mass
θ\theta Eq. (30) Heaviside step function ξ\xi Eq. (31) boxcar function
GG gravitational constant Π\Pi Δ​v/cs\Delta v/c_{s} dimensionless gas velocity reduction
VV protoplanet volume mjm_{j} mass of protoplanet in constituent jj
ρ\rho protoplanet density VjV_{j} mj/ρjm_{j}/\rho_{j} volume occupied by constituent jj
FjF_{j} mj/Mpm_{j}/M_{p} mass fraction of constituent jj fjf_{j} Vj/VV_{j}/V volume fraction of constituent jj
tsettlet_{\rm settle} Eq. (33) pebble settling time [26Al] isotopic abundance of 26Al
ℋ\mathcal{H} heat production rate per mass cpc_{p} heat capacity of water ice
t1/2t_{1/2} half-life of 26Al λ\lambda ln⁡(2)/t1/2\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 r0r_{0} and orbiting at the corresponding Keplerian frequency Ω0\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

𝒟​ρg𝒟​t\displaystyle\frac{{\mathcal{D}}\rho_{g}}{{\mathcal{D}}t} =\displaystyle= −ρg​∇⋅𝒖,\displaystyle-\rho_{g}{\boldsymbol{\nabla}}\cdot{\boldsymbol{u}}, (1)
𝒟​𝒖𝒟​t\displaystyle\frac{{\mathcal{D}}{\boldsymbol{u}}}{{\mathcal{D}}t} =\displaystyle= −1ρg​∇p−2​Ω0​𝒛^×𝒖+32​Ω0​ux​𝒚^−Ω02​𝒛+𝒇g\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)
+2​Ω0​Δ​v​𝒙^,\displaystyle+2\varOmega_{0}\Delta v\hat{{\boldsymbol{x}}},
d​𝒙d​t\displaystyle\frac{d{\boldsymbol{x}}}{dt} =\displaystyle= 𝒗−32​Ω​x​𝒚^,\displaystyle{\boldsymbol{v}}-\frac{3}{2}\varOmega x\hat{{\boldsymbol{y}}}, (3)
d​𝒗d​t\displaystyle\frac{d{\boldsymbol{v}}}{dt} =\displaystyle= −2​Ω0​𝒛^×𝒗+32​Ω0​vx​𝒚^−Ω02​𝒛+𝒇d−∇Φ,\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)
∇2Φ\displaystyle\nabla^{2}{\varPhi} =\displaystyle= 4​π​G​ρd.\displaystyle 4\pi G\rho_{d}. (5)

Here, ρg\rho_{g} is the gas density, tt is time, 𝒖{\boldsymbol{u}} is the gas velocity, pp is the gas pressure, 𝒙{\boldsymbol{x}} and 𝒗{\boldsymbol{v}} are the position and velocity vectors of the pebbles, respectively, ρd\rho_{d} is the volume density of pebbles, Φ\varPhi is the selfgravity potential, GG is the gravitational constant, and (xx,yy,zz) are the local Cartesian coordinates of the shearing box. The terms fdf_{d} and fgf_{g} are the drag force and its backreaction, respectively, explained in Appendix A.

The quantity Δ​v=η​vK\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),

η≡−12​ρg​ΩK2​r​∂p∂r.\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

𝒟𝒟​t=∂∂t+u0​∂∂y+𝒖⋅∇,\frac{{\mathcal{D}}\ }{{\mathcal{D}}t}=\frac{\partial}{\partial t}+u_{0}\frac{\partial}{\partial y}+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}, (7)

where u0=−(3/2)​Ω0​xu_{0}=-\left(3/2\right)\varOmega_{0}x is the linearized Keplerian shear flow. The pressure is related to the sound speed csc_{s} by the equation of state

p=ρg​cs2/γ,p=\rho_{g}c_{s}^{2}/\gamma, (8)

which closes the system. Here γ=1\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

ρg​(z)=ρ0​exp⁡(−z22​H2),\rho_{g}\left(z\right)=\rho_{0}\exp\left(-\frac{z^{2}}{2H^{2}}\right), (9)

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

H≡csΩ0H\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, λSI≈r​η\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 λSI≈h​H\lambda_{\rm SI}\approx hH, where the disk aspect ratio h≡H/r0h\equiv H/r_{0}. The sound speed is, for constant adiabatic index γ\gamma and mean molecular weight μ\mu, solely dependent on temperature

cs2=γ​kB​Tμ​mH.c_{s}^{2}=\frac{\gamma k_{B}T}{\mu m_{H}}. (11)

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

ΩK2≡G​Mr3,\varOmega^{2}_{K}\equiv\frac{GM}{r^{3}}, (12)

where MM is the mass of the host star, and rr is the astrocentric distance. At 45 AU for a solar mass star, the Keplerian frequency is ΩK=6.596×10−10​ s−1\varOmega_{K}=6.596\times 10^{-10}\text{ s}^{-1} resulting in a scale height of H≈3.3​ AUH\approx 3.3\text{\,AU}. Thus, we have the characteristic length of the streaming instability being η​r≈0.07​H\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.2​HL=0.2H (0.66​ AU0.66\text{\,AU} for H=3.3​ AUH=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 λSI\lambda_{\rm SI}, the required resolution of our mesh is primarily constrained by the dust scale height, HdH_{d}, which is determined as a balance between vertical stellar gravity and turbulent diffusion

Hd=H​(1+Stδ)−1/2,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)

α≡⟨ρg​δ​vx​δ​vy⟩p.\alpha\equiv\frac{\langle\rho_{g}\,\delta v_{x}\,\delta v_{y}\rangle}{p}. (14)

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

δ=α​(1+St2)−1.\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 α≈10−5\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 St=0.5\mathrm{St}=0.5, this results in a dust scale height of Hd≈H/250H_{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 Nx=Ny=Nz=256N_{x}=N_{y}=N_{z}=256 resulting in a points-per-scale-height value of H/Δ=1280H/\Delta=1280.

2.3 Simulation Parameters

We set the dust-to-gas mass ratio (metallicity) to be Z=0.03Z=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 Np=1.536×107N_{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

St\displaystyle\mathrm{St} ≈\displaystyle\approx 0.1​(a1​cm)​(ρ∙1​g​cm−3)​(T35​K)−1/2​(ρg10−13​g​cm−3)−1\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 St=0.5\mathrm{St}=0.5, representing ices, and the other half a Stokes number of St=5×10−3\mathrm{St}=5\times 10^{-3}, representing silicates.

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

Π≡Δ​vcs,\Pi\equiv\frac{\Delta\,v}{c_{s}}, (17)

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

G~≡4​π​G​ρ0Ω02,\widetilde{G}\equiv\frac{4\pi G\rho_{0}}{\varOmega_{0}^{2}}, (18)

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

Q=8π​G~−1.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​π​G=0.14\pi G=0.1, consequently setting the gravitational constant in code units to be G=0.1/4​πG=0.1/4\pi. Furthermore, with ρ0=Ω0=1\rho_{0}=\varOmega_{0}=1, then G~=0.1\widetilde{G}=0.1 and Q≈16Q\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,

ρR≡9​Ω024​π​G=9​ρ0​G~−1.\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

Refer to caption
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 μ\mum 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 ρ∙=3\rho_{\bullet}=3 g cm-3.
Refer to caption
Figure 3: The accretion rates for the three dust distributions at t=0t=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

ρ∙​(a)=ρi​(aamax)−q1,\rho_{\bullet}(a)=\rho_{i}\left(\frac{a}{a_{\rm max}}\right)^{-q_{1}}, (21)

with

q1≡log⁡(ρs/ρi)log⁡(amax/amin)>0q_{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 ρs\rho_{s}, and ρ∙\rho_{\bullet} decreases in internal density with increasing grain size such that the largest particle has a density ρi\rho_{i}. We choose for this model amin=1​μa_{\rm min}=1\,\mum and amax=1a_{\rm max}=1\,cm.

The second model is our physically motivated bimodal distribution

ρ∙​(a)={ρs,if ​a<a0ρi​(aamax)−q2,if ​a≥a0\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

q2≡log⁡(ρs/ρi)log⁡(amax/a0)>0.q_{2}\equiv\frac{\log(\rho_{s}/\rho_{i})}{\log(a_{\rm max}/a_{0})}>0. (24)

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

The last model assumes that all pebbles are silicates, with internal density ρs\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

ρ=ρi​fi+ρs​fs\displaystyle\rho=\rho_{i}f_{i}+\rho_{s}f_{s} (25)
fi+fs+fv=1\displaystyle f_{i}+f_{s}+f_{v}=1 (26)

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

fi=1−(ρ∙−ρiρs−ρi)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 ρi=1\rho_{i}=1 g cm-3 and ρs=3\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 ≈\approx3 mm, but it abruptly passes to ≈\approx0.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 Δ​t\Delta t such that the mass accreted per timestep (M˙​Δ​t\dot{M}\Delta t) is no greater than a fraction CC of the planetesimal mass MpM_{p},

Δ​t=C​MpM˙.\Delta t=C\ \frac{M_{p}}{\dot{M}}. (28)

We set the value of CC 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=0t=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×10−55\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≈1500R\approx 1500 km. Thus, we assume the following simplified porosity parametrization

ϕ​(R)=ϕ0​θ​(R1−R)+ξR1,R2​(R)​log10⁡(R2R)\phi(R)=\phi_{0}\ \theta(R_{1}-R)+\xi_{R_{1},R_{2}}(R)\,\log_{10}\left(\frac{R_{2}}{R}\right) (29)

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

θ​(x)≡{1,if ​x>01/2,if ​x=00,if ​x<0\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

ξx1,x2​(x)≡θ​(x−x1)−θ​(x−x2).\xi_{x_{1},x_{2}}(x)\equiv\theta(x-x_{1})-\theta(x-x_{2}). (31)
Refer to caption
Refer to caption
Figure 4: Left: Graph of Eq. (29), the parametrization we adopt for porosity. Right: Densities, for constant ice mass fraction FiF_{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 ϕ0\phi_{0}, up to a radius R1R_{1}, beyond which the porosity is removed, logarithimically, reaching zero at R2R_{2}. Choosing R2=1500R_{2}=1500 km and ϕ0=0.5\phi_{0}=0.5 results in R1≈474R_{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 jj (ice or silicate) is defined as Fj≡mj/Mp=fj​ρj/ρF_{j}\equiv m_{j}/M_{p}=f_{j}\rho_{j}/\rho where mjm_{j} is the mass of the constituent in the body, Mp=mi+msM_{p}=m_{i}+m_{s} is the total mass, and fjf_{j} is the constituent’s volume fraction. Given fv≡ϕf_{v}\equiv\phi, Eq. (26) can be solved for ρ\rho in terms of the ice mass fraction FiF_{i}

ρ=(1−ϕ)​ρi​ρsρi+Fi​(ρs−ρ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 FiF_{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 St=0.5\rm{St}=0.5 and St=5×10−3\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​π/Ω0≈300​ yearsT=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 Ω0−1\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 ρR\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.

Refer to caption
Figure 5: Left: Planetesimal formation rate (number of planetesimals NN formed per time interval Ω0−1\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.
Refer to caption
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 Hd≈0.1H_{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)

tsettle=1+2​S​t2Ω​St,t_{\rm settle}=\frac{1+2\mathrm{St}^{2}}{\varOmega\mathrm{St}}, (33)
Refer to caption
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 ∼\sim1 orbit.
Refer to caption
Refer to caption
Figure 8: Left: Ice-to-silicate ratio in the midplane at t=0.5t=0.5 orbits, approximately the time when the ices settle. Right: Ice-to-silicate ratio in the midplane at t=2.0t=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 ∼\sim0.3 and 30 orbits, for St=0.5\mathrm{St}=0.5 and St=5×10−3\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 α∼10−6\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 ρice/ρsil\rho_{\rm ice}/\rho_{\rm sil}, where ρice\rho_{\rm ice} and ρsil\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 ∼\sim0.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.5t=0.5 and 2 orbits. These are the times at tsettlet_{\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.

Refer to caption
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 FiF_{i} of the planetesimals, at the end of the simulation. The masses range from ≈1.5×10−4\approx 1.5\times{10^{-4}} to ≈3×10−3\approx 3\times{10^{-3}} MPlutoM_{\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−210^{-2} MPlutoM_{\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−210^{-2} (2002 UX25) to 10−110^{-1} MPlutoM_{\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−110^{-1} MPlutoM_{\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.

Refer to caption
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−110^{-2}-10^{-1} MPlutoM_{\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−110^{-2}-10^{-1} MPlutoM_{\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​MPluto0.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 ∼10−2\sim 10^{-2} MPlutoM_{\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​MPlutoM=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−110^{-1} MPlutoM_{\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−210^{-2} MPlutoM_{\rm Pluto}\ at 20 AU, model 2 is accreting almost pure silicates, like model 3. At 10−110^{-1} MPlutoM_{\rm Pluto}\ the models diverge significantly as the larger pebbles are icy in model 2.

Refer to caption
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).
Refer to caption
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​MPluto10^{-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 ≈1.5×10−4\approx 1.5\times{10^{-4}} to ≈3×10−3\approx 3\times{10^{-3}} MPlutoM_{\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 RHR_{H} is equal to the length Δ\Delta of a mesh cell is

Mp(RH=Δ)\displaystyle M_{p}^{\left(R_{H}=\Delta\right)} =\displaystyle= 3​M☉​(Δr)3\displaystyle 3\,M_{\sun}\left(\frac{\Delta}{r}\right)^{3}
≈\displaystyle\approx 8.6×10−5​MPluto\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 ∝R\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​MPluto10^{-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​MPluto10^{-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.640.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%>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×1057\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 ≈\approx25 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 tgrow∼ε−1​Torbt_{\rm grow}\sim\varepsilon^{-1}T_{\rm orb} (Birnstiel et al., 2012; Lambrechts & Johansen, 2014; Lorek & Johansen, 2022), where TorbT_{\rm orb} is the orbital period. The disk starts with ε∼10−2\varepsilon\sim 10^{-2} of dust, of interstellar (sub-micron) size. Assuming an orbital period Torb∼100T_{\rm orb}\sim 100 years yields tgrow∼104t_{\rm grow}\sim 10^{4} years. The dust thus coagulates very quickly into pebbles. The steady state, however, is top-heavy, leaving about ε∼10−4\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 ∼\sim1 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​ g cm−32.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 Mp≲10−2​MPlutoM_{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 11 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 kk in the simulation. The drag force accelerating particle kk is

𝒇d,k=−𝒗k−𝒖¯​(𝒙k)τk{\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 kk is reduced by a factor of ee. 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

τk=ak​ρ∙,kvth​ρg¯​(𝒙k).\tau_{k}=\frac{a_{k}\ \rho_{\bullet,k}}{v_{\rm th}\overline{\rho_{g}}({\boldsymbol{x}}_{k})}. (A2)

Here aa is the grain radius, ρ∙\rho_{\bullet} is the grain internal density and vth=8/π​csv_{\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 St\mathrm{St},

Stk≡Ω0​τk.\mathrm{St}_{k}\equiv\varOmega_{0}\tau_{k}. (A3)

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

ψ¯​(𝒙k)=∑W​(𝒙k−𝒙)​ψ​(𝒙)\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 WW 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

ρd​(𝒙)=1V​(𝒙)​∑kW​(𝒙k−𝒙)​mk,\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 ∇Φ¯​(𝒙k)\overline{{\boldsymbol{\nabla}}{\varPhi}}({\boldsymbol{x}}_{k}) to the position of each particle kk is done according to Eq. (A4) and added to Eq. (4).

The drag force backreaction 𝒇g{\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)

𝒇g​(𝒙)=−1ρg​(𝒙)​V​(𝒙)​∑kW​(𝒙k−𝒙)​mk​fd,k{\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​(𝒙)V({\boldsymbol{x}}) is the mesh cell volume and mkm_{k} the mass of particle kk. In this formulation, the dragforce backreaction for multiple particle species is straightforward, as the sum collects the individual contribution of each particle kk.

Appendix B Heating from decay of 26Al

Refer to caption
Figure 13: Body temperature T1T_{1} resulting from heating from decay of 26Al, as a function of mass fraction of silicates, FsF_{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

ℋ=ρFs[26Al]0ℋ0e−λ​t\mathcal{H}=\rho\,F_{s}\,{\rm[^{26}Al]}_{0}\,\mathcal{H}_{0}\,e^{-\lambda t} (B1)

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

Q​(t)\displaystyle Q(t) =\displaystyle= V​∫0tℋ​(t′)​𝑑t′\displaystyle V\int_{0}^{t}\mathcal{H}(t^{\prime})dt^{\prime} (B2)
=\displaystyle= MpFs[26Al]0ℋ0λ−1(1−e−λ​t)\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 Mp=ρ​VM_{p}=\rho V. If we assume that all this energy is trapped within the body, we can find the temperature T1T_{1} achieved by setting t→∞t\rightarrow\infty, and equating with the heat capacity equation

Q=Mp​cp​Δ​TQ=M_{p}\,c_{p}\,\Delta T (B4)

where cpc_{p} is the heat capacity at constant pressure. Then solving for Δ​T\Delta T

ΔT=Fs[26Al]0ℋ0λ−1cp−1\Delta T=F_{s}\,{\rm[^{26}Al]}_{0}\,\mathcal{H}_{0}\,\lambda^{-1}\,c_{p}^{-1} (B5)

Substituting the values cp=2108​J​kg−1​K−1c_{p}=2108\ {\rm J\,kg^{-1}\,K^{-1}} (Marquet, 2015), [26Al]0=(8.6×10−8,5.37×10−7){\rm[^{26}Al]}_{0}=\left(8.6\times 10^{-8},5.37\times 10^{-7}\right) (Parhi & Prialnik, 2023), ℋ0=0.355​W/kg\mathcal{H}_{0}=0.355\ {\rm W/kg} (Castillo-Rogez et al., 2009), t1/2=0.717​Myrt_{1/2}=0.717\ {\rm Myr}, and T0=30​KT_{0}=30\,{\rm K}, we find the FsF_{s} vs T1T_{1} (the final temperature) relationship shown in Fig. 13. The blue curve is for the lower value of [26Al]0{\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 FsF_{s}=0.06 (0.08) for the higher 26Al abundance. The lower abundance of 26Al reaches 200 K for FsF_{s}=0.36, or 250 K at FsF_{s}=0.47.

In this model, the higher abundance would not allow for porosity retention, but the lower abundance accomodates the ≳65%\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
  • Dr\każkowska & Dullemond (2014) Dr\każ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