A Solution for the Density Dichotomy Problem of Kuiper Belt Objects with Multispecies Streaming Instability and Pebble Accretion

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 self-gravity of the gas (explained a posteriori in Section 2.3) but consider the gravity of the dust grains. Under this assumption, the equations of motion for the system are

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal D }{\rho }_{g}}{{ \mathcal D }t}=-{\rho }_{g}{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{ \mathcal D }{\boldsymbol{u}}}{{ \mathcal D }t} & = & -\displaystyle \frac{1}{{\rho }_{g}}{\boldsymbol{\nabla }}p-2{{\rm{\Omega }}}_{0}\hat{{\boldsymbol{z}}}\times {\boldsymbol{u}}+\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{0}{u}_{x}\hat{{\boldsymbol{y}}}-{{\rm{\Omega }}}_{0}^{2}{\boldsymbol{z}}+{{\boldsymbol{f}}}_{g}\\ & & +2{{\rm{\Omega }}}_{0}{\rm{\Delta }}v\hat{{\boldsymbol{x}}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\boldsymbol{v}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=-2{{\rm{\Omega }}}_{0}\hat{{\boldsymbol{z}}}\times {\boldsymbol{v}}+\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{0}{v}_{x}\hat{{\boldsymbol{y}}}-{{\rm{\Omega }}}_{0}^{2}{\boldsymbol{z}}+{{\boldsymbol{f}}}_{d}-{\boldsymbol{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G{\rho }_{d}.\end{eqnarray} \tag{ 5 }$$

Here ρ_g is the gas density, t is time, u is the gas velocity, p is the gas pressure, x and v are the position and velocity vectors of the pebbles, respectively, ρ_d is the volume density of pebbles, Φ is the self-gravity potential, G is the gravitational constant, and (x,y,z) are the local Cartesian coordinates of the shearing box. The terms f_d and f_g are the drag force and its back-reaction, respectively, explained in Appendix A.

The quantity Δv = η v_K is the orbital velocity reduction of the gas due to the pressure support, and η is related to the global-scale radial pressure gradient (Nakagawa et al. 1986),

$$\begin{eqnarray}&&\eta \equiv -\displaystyle \frac{1}{2{\rho }_{g}{{\rm{\Omega }}}_{{\rm{K}}}^{2}r}\displaystyle \frac{\partial p}{\partial r}.\end{eqnarray} \tag{ 6 }$$

A table of symbols is provided in Table 1. We have defined the operator

$$\begin{eqnarray}&&\frac{{ \mathcal D }=}{{ \mathcal D }t}\frac{\partial }{\partial t}+{u}_{0}\frac{\partial }{\partial y}+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }},\end{eqnarray} \tag{ 7 }$$

where ${u}_{0}=-\left(3/2\right){{\rm{\Omega }}}_{0}x$ is the linearized Keplerian shear flow. The pressure is related to the sound speed c_s by the equation of state

$$\begin{eqnarray}&&p={\rho }_{g}{c}_{s}^{2}/\gamma ,\end{eqnarray} \tag{ 8 }$$

which closes the system. Here γ = 1 is the adiabatic index for an isothermal gas.

Table 1. Symbols Used in This Work

Symbol	Definition	Description	Symbol	Definition	Description
t		time	r0		arbitrary origin of shearing box
x, y, z		Cartesian coordinates	p	Equation (8)	gas pressure
ε	ρd /ρg	dust-to-gas density ratio	cs	Equation (11)	sound speed
Z	$\tfrac{\int {\rho }_{d}{dV}}{\int {\rho }_{g}{dV}}$	dust-to-gas mass ratio	dV	dx dy dz	volume infinitesimal
vK	Ω0 r0	Keplerian velocity	Δv	η vK	gas velocity reduction
γ		adiabatic index	τ	Equation (A2)	particle friction time
vth	$\sqrt{\tfrac{8}{\pi }}{c}_{s}$	thermal velocity	ρ•		pebble material density
kB		Boltzmann constant	v	(vx , vy , vz )	pebble velocity
T		temperature	λSI	η r	streaming instability length scale
mH		atomic mass unit	Hd	Equation (13)	pebble scale height
ΩK	Equation(12)	Keplerian frequency	N	Nx , Ny , Nz	box resolution
M		stellar mass	$\widetilde{G}$	Equation (18)	dimensionless gravitational constant
δ		dimensionless pebble diffusion	ρd		volume density of pebbles
Δ		mesh spacing	a		pebble radius
ϕ		porosity	St	Ω0 τ	Stokes number of pebbles
R		protoplanet radius	Φ		self-gravitational potential
ρs		silicate density	μ		mean molecular weight
ρi		ice density	L	Lx = Ly = Lz	length of shearing box sides
H	cs /Ω0	gas scale height	α	Equation (14)	dimensionless gas viscosity
u	(ux , uy , uz )	gas velocity in Cartesian coordinates	Np		number of numerical particles
Ω0		Keplerian frequency at r0	Q	Equation (19)	Toomre Q parameter
ρg		gas density	ρR	Equation (20)	Roche density
h	H/r0	disk aspect ratio	$\dot{M}$		mass accretion rate
η	Equation (6)	dimensionless pressure support	r		astrocentric distance
u0	$-\tfrac{3}{2}{{\rm{\Omega }}}_{0}x$	Keplerian shear flow	Mp		protoplanet mass
θ	Equation (30)	Heaviside step function	ξ	Equation (31)	boxcar function
G		gravitational constant	Π	Δv/cs	dimensionless gas velocity reduction
V		protoplanet volume	mj		mass of protoplanet in constituent j
ρ		protoplanet density	Vj	mj /ρj	volume occupied by constituent j
Fj	mj /Mp	mass fraction of constituent j	fj	Vj /V	volume fraction of constituent j
tsettle	Equation (33)	pebble settling time	[26Al]		isotopic abundance of 26Al
${ \mathcal H }$		heat production rate per mass	cp		heat capacity of water ice
t1/2		half-life of 26Al	λ	$\mathrm{ln}(2)/{t}_{1/2}$	decay constant of 26Al

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&{\rho }_{g}\left(z\right)={\rho }_{0}\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}^{2}}\right),\end{eqnarray} \tag{ 9 }$$

where ρ₀ is the midplane gas density and

$$\begin{eqnarray}&&H\equiv \displaystyle \frac{{c}_{s}}{{{\rm{\Omega }}}_{0}}\end{eqnarray} \tag{ 10 }$$

is the gas scale height.

We choose the length of our box such that the scales of the streaming instability, λ_SI ≈ r η, are well contained within the box. Given Equation (10), we can rewrite this length in terms of the scale height as λ_SI ≈ hH, where the disk aspect ratio h ≡ H/r₀. The sound speed is, for constant adiabatic index γ and mean molecular weight μ, solely dependent on temperature,

$$\begin{eqnarray}&&{c}_{s}^{2}=\displaystyle \frac{\gamma {k}_{{\rm{B}}}T}{\mu {m}_{{\rm{H}}}}.\end{eqnarray} \tag{ 11 }$$

We substitute the values γ = 1 for the adiabatic index of an isothermal flow, μ = 2.3 for the mean molecular weight (i.e., five H₂ for every two He), and k_B and m_H are the Boltzmann constant and atomic mass unit, respectively. The disk temperature T is found using the temperature relation of Chiang & Goldreich (1997). In this model at r = 45 au, the disk temperature is T ≈ 30 K, which results in a sound speed of c_s ≈ 0.3 km s⁻¹. Similarly, we find the Keplerian frequency

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{K}}}^{2}\equiv \displaystyle \frac{{GM}}{{r}^{3}},\end{eqnarray} \tag{ 12 }$$

where M is the mass of the host star and r is the astrocentric distance. At 45 au for a solar-mass star, the Keplerian frequency is Ω_K = 6.596 × 10⁻¹⁰ s⁻¹, resulting in a scale height of H ≈ 3.3 au. Thus, we have the characteristic length of the streaming instability being η r ≈ 0.07H, placing an upper limit for the distance between mesh cells, Δ. Even for low resolutions of 16 mesh cells along each axis, a cubic box of sides L = 0.2H (0.66 au for H = 3.3 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, the required resolution of our mesh is primarily constrained by the dust scale height, H_d , which is determined as a balance between vertical stellar gravity and turbulent diffusion

$$\begin{eqnarray}&&{H}_{d}=H{\left(1+\displaystyle \frac{\mathrm{St}}{\delta }\right)}^{-1/2},\end{eqnarray} \tag{ 13 }$$

where δ 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 α viscosity (Shakura & Sunyaev 1973)

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{\langle {\rho }_{g}\,\delta {v}_{x}\,\delta {v}_{y}\rangle }{p}.\end{eqnarray} \tag{ 14 }$$

Here δ v_i is the turbulent deviation from the mean flow in the direction i. The α-viscosity parameter is related to δ by (Youdin & Lithwick 2007)

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

Equation (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. In addition, Equation (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 α-viscosity around α ≈ 10⁻⁵. ⁹ For Stokes number of St = 0.5, this results in a dust scale height of H_d ≈ H/250. We want to resolve this layer with at least six mesh points (the size of a stencil), and to this end, we have set the resolution of our simulation to be N_x = N_y = N_z = 256, resulting in a points-per-scale-height value of H/Δ = 1280.

We set the dust-to-gas mass ratio (metallicity) to be Z = 0.03, slightly above solar, which is known to trigger strong clumping by the streaming instability (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021). The bulk mass of solids is evenly split into N_p = 1.536 × 10⁷ numerical superparticles; 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 centimeter and submillimeter size scales, respectively. Given

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{St}\approx 0.1\left(\displaystyle \frac{a}{1\ \mathrm{cm}}\right)\left(\displaystyle \frac{{\rho }_{\bullet }}{1\ {\rm{g}}\ {\mathrm{cm}}^{-3}}\right)\\ {\left(\displaystyle \frac{T}{35\ {\rm{K}}}\right)}^{-1/2}{\left(\displaystyle \frac{{\rho }_{g}}{{10}^{-13}\ {\rm{g}}\ {\mathrm{cm}}^{-3}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 16 }$$

we assign half the particles a Stokes number of St = 0.5, representing ices, and the other half a Stokes number of St = 5 × 10⁻³, representing silicates.

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

$$\begin{eqnarray}&&{\rm{\Pi }}\equiv \displaystyle \frac{{\rm{\Delta }}\,v}{{c}_{s}},\end{eqnarray} \tag{ 17 }$$

and set it to Π = 0.05 (note also that Πc_s = η v_K). The relative strength of self-gravity to tidal shear is determined by another dimensionless parameter,

$$\begin{eqnarray}&&\widetilde{G}\equiv \displaystyle \frac{4\pi G{\rho }_{0}}{{{\rm{\Omega }}}_{0}^{2}},\end{eqnarray} \tag{ 18 }$$

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

$$\begin{eqnarray}&&Q=\sqrt{\displaystyle \frac{8}{\pi }}{\widetilde{G}}^{-1}.\end{eqnarray} \tag{ 19 }$$

The constants on the right-hand side of Equation (5) are nondimensionalized and in code units are set to be 4π G = 0.1, consequently setting the gravitational constant in code units to be G = 0.1/4π. Furthermore, with ρ₀ = Ω₀ = 1, then $\widetilde{G}=0.1$ and Q ≈ 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,

$$\begin{eqnarray}&&{\rho }_{{\rm{R}}}\equiv \frac{9{{\rm{\Omega }}}_{0}^{2}}{4\pi G}=9{\rho }_{0}{\widetilde{G}}^{-1}.\end{eqnarray} \tag{ 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.

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 centimeter 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). In addition, Musiolik (2021) reports higher surface energies for irradiated ice owing 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).

Once gravitationally bound objects are produced, they continue to grow by accreting pebbles (Lyra et al. 2008b, 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

$$\begin{eqnarray}&&{\rho }_{\bullet }(a)={\rho }_{i}{\left(\displaystyle \frac{a}{{a}_{\max }}\right)}^{-{q}_{1}},\end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray}&&{q}_{1}\equiv \displaystyle \frac{\mathrm{log}({\rho }_{s}/{\rho }_{i})}{\mathrm{log}({a}_{\max }/{a}_{\min })}\gt 0,\end{eqnarray} \tag{ 22 }$$

so that the smallest grain has a density ρ_s and ρ_• decreases in internal density with increasing grain size such that the largest particle has a density ρ_i . We choose for this model ${a}_{\min }=1\,\mu {\rm{m}}$ and ${a}_{\max }=1\,$ cm.

The second model is our physically motivated bimodal distribution

$$\begin{eqnarray}{\rho }_{\bullet }(a)=\left\{\begin{array}{ll}{\rho }_{s}, & \mathrm{if}\ a\lt {a}_{0}\\ {\rho }_{i}{\left(\displaystyle \frac{a}{{a}_{\max }}\right)}^{-{q}_{2}}, & \mathrm{if}\ a\geqslant {a}_{0}\end{array}\right.,\end{eqnarray} \tag{ 23 }$$

with

$$\begin{eqnarray}&&{q}_{2}\equiv \displaystyle \frac{\mathrm{log}({\rho }_{s}/{\rho }_{i})}{\mathrm{log}({a}_{\max }/{a}_{0})}\gt 0.\end{eqnarray} \tag{ 24 }$$

We choose for this model a₀ = 1 mm, that is, particles between the sizes 1 μm ≤ a ≤ 1 mm are assigned internal density ρ_s . Pebbles beyond 1 mm in size then follow a power law such that the largest pebble has internal density ρ_i .

The last model assumes that all pebbles are silicates, with internal density ρ_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 Figure 2, where each column represents the different models. The top 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

$$\begin{eqnarray}&&\rho ={\rho }_{i}{f}_{i}+{\rho }_{s}{f}_{s}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{f}_{i}+{f}_{s}+{f}_{v}=1,\end{eqnarray} \tag{ 26 }$$

where f_j ≡ V_j /V is the volume fraction of component j, V_j is the volume occupied by the component, and V is the total volume. Here f_v is the fractional volume of empty space, or "porosity." For compact pebbles, f_v = 0 (and hence ρ = ρ_•). Solving for the ice volume fraction of a pebble,

$$\begin{eqnarray}&&{f}_{i}=1-\left(\displaystyle \frac{{\rho }_{\bullet }-{\rho }_{i}}{{\rho }_{s}-{\rho }_{i}}\right),\end{eqnarray} \tag{ 27 }$$

which is the quantity plotted in the bottom row of Figure 2. We use ρ_i = 1 g cm⁻³ and ρ_s = 3 g cm⁻³.

Zoom In Zoom Out Reset image size

The accretion rates of each model at 20 au are illustrated in Figure 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 Figure 3 that the bimodal distribution (model 2) shows a window of silicate accretion. Model 2 begins accreting pebbles of ≈50%/50% ice/silicate 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 discontinuous (Ormel 2017). In our model, the best accreted pebble at the geometric regime is of radius ≈3 mm, but it abruptly passes to ≈0.5 mm after the onset of Bondi accretion. The window of silicate accretion starts to narrow at higher protoplanetary masses that are able to 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\rm{\Delta }}t=C\ \displaystyle \frac{{M}_{p}}{\dot{M}}.\end{eqnarray} \tag{ 28 }$$

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

We highlight and address here the inconsistency of using a two-species pebble model for streaming instability, while using a continuous size distribution for the pebble accretion calculation. While inconsistent, this was done because our goal with the hydrodynamical simulation was to provide a proof of concept that the first planetesimals would be mostly icy. Thus, binning the pebbles in two species—one rocky and one icy—was judged a satisfactory first-order approximation (but see Schaffer et al. 2018, 2021; Krapp et al. 2019; Paardekooper et al. 2020, 2021; Yang & Zhu 2021; Zhu & Yang 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.

Finally, we must consider a compaction model to grasp the density evolution of KBOs. 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 toward the center becomes strong enough 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 nonnegligible amount of the short-lived radioactive isotope ²⁶Al (initial abundance ²⁶Al/²⁷Al = 5 × 10⁻⁵; 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 nonlinear 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 ≈ 1500 km. Thus, we assume the following simplified porosity parameterization:

$$\begin{eqnarray}&&\phi (R)={\phi }_{0}\ \theta ({R}_{1}-R)+{\xi }_{{R}_{1},{R}_{2}}(R)\,{\mathrm{log}}_{10}\left(\displaystyle \frac{{R}_{2}}{R}\right),\end{eqnarray} \tag{ 29 }$$

àwhere R is the radius of the KBO, R₂ = 1500 km, and ${R}_{1}={R}_{2}/{10}^{{\phi }_{0}}$. Here θ is the Heaviside step function,

$$\begin{eqnarray}\theta (x)\equiv \left\{\begin{array}{ll}1, & \mathrm{if}\ x\gt 0\\ 1/2, & \mathrm{if}\ x=0\\ 0, & \mathrm{if}\ x\lt 0\end{array}\right.,\end{eqnarray} \tag{ 30 }$$

and ξ is the boxcar function,

$$\begin{eqnarray}&&{\xi }_{{x}_{1},{x}_{2}}(x)\equiv \theta (x-{x}_{1})-\theta (x-{x}_{2}).\end{eqnarray} \tag{ 31 }$$

Equation (29) is plotted in the left panel of Figure 4. It states that the porosity is constant at the porosity of formation ϕ₀, up to a radius R₁, beyond which the porosity is removed, logarithmically, reaching zero at R₂. Choosing R₂ = 1500 km and ϕ₀ = 0.5 results in R₁ ≈ 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 Equation (29) to determine the porosity fraction. To aid the interpretation of our results, we first isolate the effect the porosity function would have on bodies of constant composition. The mass fraction of a constituent j (ice or silicate) is defined as F_j ≡ m_j /M_p = f_j ρ_j /ρ, where m_j is the mass of the constituent in the body, M_p = m_i + m_s is the total mass, and f_j is the constituent's volume fraction. Given f_v ≡ ϕ, Equation (26) can be solved for ρ in terms of the ice mass fraction F_i ,

$$\begin{eqnarray}&&\rho =\displaystyle \frac{\left(1-\phi \right){\rho }_{i}{\rho }_{s}}{{\rho }_{i}+{F}_{i}\left({\rho }_{s}-{\rho }_{i}\right)}.\end{eqnarray} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

Curves of constant F_i calculated according to Equation (32) are shown in the right panel of Figure 4. The solid lines are porous bodies with porosity given by Equation (29); dashed horizontal lines mark fully compact (zero-porosity) bodies.

We run a hydrodynamic simulation where we consider two pebble species, ices and silicates, with Stokes number St = 0.5 and St = 5 × 10⁻³, 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 yr, considering that an orbit at 45 au is T = 2π/Ω₀ ≈ 300 yr.

This is visualized in Figure 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 ${{\rm{\Omega }}}_{0}^{-1}$ (≈50 yr). The last planetesimal was formed right after five orbits (the small spike before the rate last goes to zero). The right panel shows that between the second and fourth orbits 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, the particles within the accretion radius at formation (set to the mesh spacing Δ) are removed and replaced by a sink particle.

Zoom In Zoom Out Reset image size

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.

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 midplane 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 Figure 6, where the top row of panels shows the integrated column density along the z-axis and the bottom row of panels shows 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.

Zoom In Zoom Out Reset image size

We see in the column density plots that the filamentary structure associated with the streaming instability is apparent only in the ice plots, whereas the silicate plot shows a smoother distribution. This is because the silicates are too tightly coupled to the gas and do not drift as much as the ices, being thus less prone to the streaming instability (e.g., Yang & Zhu 2021). The azimuthal average plots show that the silicates have a height of H_d ≈ 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 versus 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)

$$\begin{eqnarray}&&{t}_{\mathrm{settle}}=\displaystyle \frac{1+2{\mathrm{St}}^{2}}{{\rm{\Omega }}\mathrm{St}},\end{eqnarray} \tag{ 33 }$$

which yields ∼0.3 and 30 orbits for St = 0.5 and St = 5 × 10⁻³, 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⁻⁶. As a result, the silicate grains do not evolve significantly vertically. The ice pebbles settle very fast. We show in Figure 7 the time evolution of the ice-to-silicate ratio (defined as ρ_ice/ρ_sil, where ρ_ice and ρ_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 ∼0.5 orbits, as expected. At one to two orbits, the streaming instability develops and saturates.

Zoom In Zoom Out Reset image size

We show in Figure 8 the midplane snapshots of the ice/silicate ratio at t = 0.5 and 2 orbits. These are the times at t_settle for the ices and when the first planetesimals form. The left panel shows the ratio resulting from settling alone. As we see it, it looks homogeneous, with an ice-to-silicate ratio of 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.

Zoom In Zoom Out Reset image size

Finally, we show in Figure 9 the mass and ice mass fraction F_i of the planetesimals, at the end of the simulation. The masses range from ≈1.5 × 10⁻⁴ M_Pluto to ≈3 × 10⁻³ M_Pluto, with a trend that the lower-mass planetesimals are more silicate-rich and the larger ones are more ice-rich. The trend does not seem to be a numerical artifact, a point we will go back to in Section 4.2. We estimate in Appendix B whether this ice mass fraction would lead to melting from heating from ²⁶Al.

Zoom In Zoom Out Reset image size

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" KBO population (Brown 2001; Kavelaars et al. 2008; Petit et al. 2011). In fact, in the context of the "Nice" model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005; see Section 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 farther 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⁻² M_Pluto (following the composition trend found in the ice and silicate streaming instability 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 top row of Figure 10 the result for the power-law distribution (model 1), at 20 au. The panels are a time series, showing the mass versus density evolution, and with the ice mass fraction of the protoplanets color-coded. The actual KBO data from Figure 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.

Zoom In Zoom Out Reset image size

The results change considerably when using the bimodal density distribution (model 2, middle row of Figure 10). In the figure, we see that the growth rate matches the mass–density trend, providing an excellent fit to the intermediate-mass objects in the range from 10⁻² M_Pluto (2002 UX₂₅) to 10⁻¹ M_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 third model, which consists of accretion of pure silicate pebbles. The results are shown in the bottom panels of Figure 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⁻¹ M_Pluto (evidencing that model 2 was accreting silicates in this mass range, as we will expand on in Section 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. However, it does reproduce the high density of Eris.

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 to 30 au, at 5 au intervals. We illustrate the results in Figure 11. This figure shows that we are able to best reproduce the mass ranges of KBOs if they formed between 15 and 22 au. At 10 au, we see that within 500,000 yr masses attained through pebble accretion are a factor of 10 larger than the mass of Pluto. Considering that 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.2M_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 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.

Zoom In Zoom Out Reset image size

To understand these results, we show in Figure 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⁻² M_Pluto, roughly the mass where KBOs 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 Section 2.5) shifts to higher masses with increasing orbital distance. The result is that at 10 au planetesimals with mass M = 10⁻² M_Pluto have already missed the window of silicate accretion, resulting in the lower densities seen in the left panel of Figure 11. Conversely, at 30 au, the window of silicate accretion extends beyond 10⁻¹ M_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.

Zoom In Zoom Out Reset image size

The figure also shows why both model 2 and model 3 reproduce the mass–density trend for intermediate masses. This is because at 10⁻² M_Pluto at 20 au model 2 is accreting almost pure silicates, like model 3. At 10⁻¹ M_Pluto the models diverge significantly, as the larger pebbles are icy in model 2.

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 Figure 9. The span in mass is about an order of magnitude, from ≈1.5 × 10⁻⁴ M_Pluto to ≈3 × 10⁻³ M_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 time span of the hydrodynamical simulation. In addition, the Hill radii of these planetesimals are resolved. For 45 au and at this resolution, the mass at which the Hill radius R_H is equal to the length Δ of a mesh cell is

$$\begin{eqnarray}\begin{array}{rcl}{M}_{p}^{\left({R}_{{\rm{H}}}={\rm{\Delta }}\right)} & = & 3\,{M}_{\odot }{\left(\displaystyle \frac{{\rm{\Delta }}}{r}\right)}^{3}\\ & \approx & 8.6\times {10}^{-5}\,{M}_{\mathrm{Pluto}},\end{array}\end{eqnarray} \tag{ 34 }$$

or about half the mass of the smallest planetesimal formed. The trend seen in Figure 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 ≈ 160 km, mean density ≈ 2.1 g cm⁻³; Trujillo et al. 2000; Noll et al. 2004). The effect of radiogenic heating is a significant function of ice mass fraction because ²⁶Al 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 dependency), but Davidsson et al. (2016) note that small objects have poorer 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⁻³ M_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 Figure 1. All streaming instability products more massive than 10⁻³ M_Pluto are ice-rich and should avoid melting.

We take this planetesimal population and feed it into a polydisperse pebble accretion calculation, at different locations, using three different models for a pebble size–composition relation, as shown in Figure 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 ice line opens the possibility that CO ice could also be relevant for the composition of the pebbles and hence the densities of KBOs. While CO is a hypervolatile, CO reacts with water to produce methanol (Grundy et al. 2020), which is refractory. Having a density of 0.64 g cm⁻³ 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 due to the excess of carbon and lack of nitrogen found in the ice giants (Ali-Dib et al. 2014) but also because ice giants are thought to contain >10% methane by mass, which would occur if the ice giants formed between the CO and N₂ 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 Section 3.2 (see also Malhotra 1993; Stern et al. 2018; Canup et al. 2021).

Our model also has the advantage of not needing to invoke a special formation time for the KBOs, contrasting with the model from Bierson & Nimmo (2019). This is because silicates, or more specifically the short-lived radioactive isotope ²⁶Al with a half-life of 7 × 10⁵ 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 ²⁶Al 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 enough 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).

Our model is in agreement with the "Nice" model scenario (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis 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 decrease 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 coplanar 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, which 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, which 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 farther 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 ≈25 au (Johansen et al. 2015). Beyond this distance, 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.