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

Kuiper Belt objects (KBOs) 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 icier. 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.

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 Pluto masses) have a nearconstant density of 0.5 g cm −3 and larger KBOs have an increasingly larger density that reaches ≈ 2.5 g cm −3 (Stansberry et al. 2006;Brown & Butler 2017;Grundy et al. 2007;Brown 2012Brown , 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) × 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 ≈1 g cm −3 ) is the lowest density material of substantial abundance in the composition of the small objects, the bulk density of 0.5 g cm −3 implies a porosity of at least 50% for them.The larger objects (e.g.Eris, Pluto, and Triton) should be less porous even in the unlikely case of 100% rocky composition.Interestingly, Baer et al. (2011) find a similar trend in the asteroid belt, with the largest asteroids showing porosities less than 20%, with an abrupt change for diameters under 300 km where a range of porosities from 0-70% is seen.
Gravitational compaction at constant composition was explored by Bierson & Nimmo (2019).Assuming a rock mass fraction of 70%, these authors successfully match 15 of 18 KBOs within two standard deviations allowed by their model (or 11 within one standard deviation).While promising, this scenario poses a timing problem requiring KBOs to be formed at least four million years after the formation of the Calcium-Aluminum-rich Inclusions (CAIs).The reason for this time constraint is that heat from the decay of the shortlived radioactive isotope 26 Al, 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(Nesvorný et al. , 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 con-taining 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. 2016Krijt et al. , 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-µm size (Krijt et al. 2018).Yet, the optical depth at a given height decreases with time as the sub-µm grains responsible for the opacity coagulate into larger grains that settle toward the midplane.As most of the solid masses is converted into larger grains, the mass remaining in sub-µm grain species decreases significantly.When the dust-to-gas ratio of these sub-µm grain species decreases to 10 −4 , typical of Class I-II disks, the gas becomes optically thin to UV as close to the midplane as 1.5 scale heights, a layer where even low levels of turbulence would suffice to loft sub-mm sized grains into.Bulk motions due to the VSI would make it even easier to get these particles to the UV-irradiated layer, and cosmic ray desorption (Silsbee et al. 2021;Sipilä et al. 2021;Rawlings 2022) would further enhance the loss of ice coat-ing 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 ą żkowska & Dullemond (2014); Carrera et al. (2015); Yang et al. (2017), andLi &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 −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.  the shearing-box approximation (Brandenburg et al. 1995;Hawley et al. 1995), centered at an arbitrary position r 0 and orbiting at the corresponding Keplerian frequency Ω 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.

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 (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 selfgravity potential, G is the gravitational constant, and (x,y,z) are the local Cartesian coordinates of the shearing box.The terms f d and f g are the drag force and its backreaction, respectively, explained in Appendix A.
The quantity ∆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), A table of symbols is provided in Table 1.We have defined the operator where u 0 = − (3/2) Ω 0 x is the linearized Keplerian shear flow.The pressure is related to the sound speed c s by the equation of state which closes the system.Here γ = 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 where ρ 0 is the midplane gas density and is the gas scale height.

Box Domain and Resolution
We choose the length of our box such that the scales of the streaming instability, λ SI ≈ rη, are well contained within the box.Given Eq. ( 10), we can rewrite this length in terms of the scale height as λ SI ≈ hH, where the disk aspect ratio h ≡ H/r 0 .The sound speed is, for constant adiabatic index γ and mean molecular weight µ, solely dependent on temperature We substitute the values γ = 1 for the adiabatic index of an isothermal flow, µ = 2.3 for the mean molecular weight (i.e., 5 H 2 for every 2 He), k B and m H are the Boltzmann constant and atomic mass unit, respectively.
The disk temperature T is found using the temperature relation of Chiang & Goldreich (1997).In this model at r = 45 AU, the disk temperature is T ≈ 30 K, which results in a sound speed of c s ≈ 0.3 km s −1 .The Keplerian frequency is 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 −10 s −1 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 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) 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) 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 α-viscosity around α ≈ 10 −5 .2For 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 6 mesh points (the size of a stencil), and to this end, we have set the resolution of our simulation to be N x = N y = N z = 256 resulting in a points-per-scale-height value of H/∆ = 1280.

Simulation Parameters
We set the dust-to-gas mass ratio (metallicity) to be Z = 0.03, slightly above Solar, which is known to trigger strong clumping by the streaming instability (Carrera et al. 2015;Yang et al. 2017;Li & Youdin 2021).The bulk mass of solids is evenly split into N p = 1.536 × 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 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 −3 , representing silicates.
For the pebble drift, we use the dimensionless parameter defined by Bai & Stone (2010) 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 and is related to the usual Toomre Q parameter (Safronov 1960;Toomre 1964) by The constants on the right-hand side of Eq. ( 5) are non-dimensionalized 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 ρ 0 = Ω 0 = 1, then 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, 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.

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(Planes et al. , 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).

Pebble Accretion
Once gravitationally bound objects are produced, they continue to grow by accreting pebbles (Lyra et al. 2008a(Lyra et al. , 2009a,b;,b;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 with 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 µm and a max = 1 cm.The second model is our physically motivated bimodal distribution with We choose for this model a 0 = 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 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 f i + ρ s f s (25) 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 which is the quantity plotted in the bottom row of Fig. 2. We use ρ i = 1 g cm −3 and ρ 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 The accretion rates for the three dust distributions at t = 0.The orange lines are the accretion rates in the Hill regime, the magenta lines are the accretion rates in the Bondi regime, and the green lines are the accretion rates in the geometric/focusing regime.The scatter points are color-coded corresponding to the ice volume fraction of the pebble of most efficient accretion at the given protoplanet mass.(Left) Accretion rates of the power law distribution, which accrete no silicates except for a small window in the transition from geometric to Bondi, where roughly 20% of the material being accreted is silicate.(Middle) The bimodal distribution begins accreting roughly 50% ice and 50% silicates in the geometric regime, but then accretes nearly 100% silicates upon entering the Bondi regime.As a planetesimal grows through Bondi accretion, it begins to favor larger pebbles, thereby increasing the ice volume fraction that is being accreted, finally accreting mostly ices in the Hill regime.(Right) Accretion rates for the constant distribution, where all pebbles being accreted have a 0% ice volume fraction.
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 ≈ 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 ≈3 mm, but it abruptly passes to ≈0.5 mm after the onset of Bondi acccretion.The window of silicate accretion starts to narrow at higher protoplanet masses, that accrete larger pebbles, of higher ice volume fraction, finally accreting mostly ices in the Hill regime.Model 1 (power law) never accretes silicates significantly, as the pure silicate pebbles are too small.The mass accretion rates are integrated numerically choosing a timestep ∆t such that the mass accreted per timestep ( Ṁ∆t) is no greater than a fraction C of the planetesimal mass M p , We set the value of C at 0.01, found empirically to be a good compromise between stability and performance.
At each time interval, we calculate the new mass and density of the planetesimal by taking a weighted average of the mass and density acquired by pebble accretion.Concurrently with pebble accretion, we reduce the gas and dust density exponentially in time, with an e-folding time of 2.5 Myr.We terminate the simulation after 4 e-folding times (10 Myr), which we quote as "the disk lifetime".We note that at 5-10 Myr there is still some gas, yet much less than at t = 0.The pebble accretion rates are impacted accordingly.To have significant mass accretion rates nearing 10 Myr is unlikely in our model (because pebbles of all sizes are essentially decoupled), even though the calculations go until this time.
We highlight and address here the inconsistency of using a two-species pebble model for streaming instability, while using a continuous size distribution for the pebble accretion calculation.While inconsistent, this was done because our goal with the hydroydynamical simulation was to provide a proof-of-concept that the first planetesimals would be mostly icy, binning the pebbles in two species, one rocky and one icy, was judged a satisfactory first order approximation (but see Schaffer et al. 2018;Krapp et al. 2019;Paardekooper et al. 2020Paardekooper et al. , 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.

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, 26 Al (initial abundance 26 Al/ 27 Al = 5 × 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 ≈ 1500 km.Thus, we assume the following simplified porosity parametrization where R is the radius of the KBO, R 2 = 1500 km, and Here θ is the Heaviside step function and ξ is the boxcar function Eq. ( 29) is plotted in Fig. 4, left panel.It states that the porosity is constant at the porosity of formation ϕ 0 , up to a radius R 1 , beyond which the porosity is removed, logarithimically, reaching zero at R 2 .Choosing R 2 = 1500 km and ϕ 0 = 0.5 results in R 1 ≈ 474 km.Although admittedly crude, this approximation is sufficient for our purposes.We apply this model concurrently with pebble accretion, using the newly calculated mass and density to obtain a radius, and then using Eq. ( 29) to determine the porosity fraction.To aid the interpretation of our results, we first isolate the effect the porosity function would have on bodies of constant composition.The mass fraction of a constituent j (ice or silicate) is defined as F j ≡ 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 ≡ ϕ, Eq. ( 26) can be solved for ρ in terms of the ice mass fraction Curves of constant F i calculated according to Eq. ( 32) are shown in Fig. 4, right panel.The full lines are porous bodies with porosity given by Eq. ( 29); dashed horizontal lines mark fully compact (zero porosity) bodies.

Planetesimal Formation
We run a hydrodynamic simulation where we consider two pebble species, ices and silicates, with Stokes number St = 0.5 and St = 5 × 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 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 Ω −1 0 (≈ 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 , the particles within the accretion radius at formation (set to the mesh spacing ∆) 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.
With icy pebbles being less susceptible to the drag force compared to silicates, they experience lower levels of turbulent diffusion, which provides support against stellar gravity.The result is that ices form a thinner, denser mid-plane layer compared to silicates (Dubrulle et al. 1995;Youdin & Lithwick 2007) and are therefore more likely to form clumps that collapse into planetesimals.This is better demonstrated in Fig. 6, where the top row of panels shows the integrated column density along the z-axis, and the bottom row the volume density averaged over the y-axis.The first column in each row shows the combined density of ices and silicates, while in the second and third columns we disaggregate the pebbles into ices and silicates, respectively.The circled objects are planetesimals and the size of the circles represents the Hill radius of the encapsulated planetesimal.
We see in the column density plots that the filamentary structure associated with the streaming instability is apparent only in the ice plots, whereas the silicate plot shows a smoother distribution.This is because the silicates are too tightly coupled to the gas and do not drift as much as the ices, being thus less prone to the streaming instability (e.g.Yang & Zhu 2021).The azimuthal average plots show that the silicates have a height of H d ≈ 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 timeevolution of the ice-to-silicate ratio.The settling time for a grain of a given Stokes number is (Youdin 2010)  29), the parametrization we adopt for porosity.Right: Densities, for constant ice mass fraction F i , calculated with Eq. ( 29) and Eq. ( 32) are shown in solid lines.The dashed lines of same color show the density of a fully compact body of same ice mass fraction.), 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.
which yields ∼0.3 and 30 orbits, for St = 0.5 and St = 5 × 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 .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 , 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 orbit, as expected.At 1-2 orbits, the streaming instability develops and saturates.We show in Fig. 8 the midplane snapshots of the ice/silicate ratio at t = 0.5 and 2 orbits.These are the times at t 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.
Finally, we show in Fig. 9 the mass and ice mass fraction F i of the planetesimals, at the end of the simulation.The masses range from ≈ 1.5 × 10 −4 to ≈ 3 × 10 −3 M 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 26 Al.

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;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.Morbidelli et al. 2005;Gomes et al. 2005, see Sect.4.5) Neptune's current location at 30 AU constrains that all KBOs aside from the cold classicals formed up to 30 AU, otherwise Neptune would have continued its migration further out (Stern et al. 2018, and references therein).Therefore, we vary the distance at which we perform the pebble accretion integration, from 10 to 30 AU, in intervals of 5 AU.We also complete the initial mass function of planetesimals up to 10 −2 M Pluto (following the composition trend found in the ices and silicates streaming intability model), given that this is the approximate mass for the onset of pebble accretion in the Bondi regime (Lyra et al. 2023).
We show in the upper row of Fig. 10 the result for the power law distribution (model 1), at 20 AU.The panels are a time series, showing the mass vs density evolution, and color-coded the ice mass fraction of the protoplanets.The actual KBO data from Fig. 1 are overplotted as green stars.We see that this model favors ices too strongly to cause any significant changes in density, even with the removal of porosity through compaction.The silicate pebbles in this model are simply too small to be accreted efficiently.Pluto mass is achieved around 4 Myr, but consisting of almost completely pure water ice.
The results change considerably when using the bimodal density distribution (model 2, middle row of Fig. 10).In the figure, we see that the growth rate matches the mass-density trend, providing an excellent fit to the intermediate-mass range objects from 10 −2 (2002 UX 25 ) to 10 −1 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 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 den- sity increases monotonically with mass as only silicates are accreted.This model also provides a good fit to the observed mass-density trend for intermediate-mass KBOs, up to 10 −1 M 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.

The effect of distance
We explore now the parameter space of distance in the pebble accretion model, taking Model 2 as fiducial and varying radii from 10 AU to 30 AU, at 5 AU intervals.We illustrate the results in Fig. 11.This figure shows that we are able to best reproduce the mass ranges of KBOs if they formed between 15 AU and 22 AU.At 10 AU, we see that within 500 thousand years, masses attained through pebble accretion are a factor of 10 larger than the mass of Pluto.Considering the lifetime of the solar nebula is of order 10 million years (i.e., a few e-folding times), it is expected that the masses at this heliocentric distance will ultimately be several orders of magnitude larger than the mass of Pluto.On the other hand, considering the accretion rates at 30 AU, we see little to no growth from pebble accretion.After 10 million years, planetesimals barely grow to 0.2 M 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.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 M Pluto , roughly the mass where Kuiper belt objects begin to show an increase in density, the accretion rate at 10 AU is roughly an order of magnitude larger than the accretion rate at 20 AU, and about two orders of magnitude larger than the accretion rate at 30 AU.Furthermore, we find that the window of enhanced silicate accretion provided by the Bondi regime (see Sect. 2.5) shifts to higher masses with increasing orbital distance.The result is that at 10 AU, planetesimals with mass M = 10 −2 M Pluto have already missed the window of silicate accretion, resulting in the lower densities seen in the left-most panel of Fig. 11.Conversely, at 30 AU, the window of silicate accretion extends beyond 10 −1 M Pluto , which could facilitate this model's ability to reproduce the density of Eris, but due to the low accretion rates, these masses are never achieved.
The figure also shows why both model 2 and model 3 reproduce the mass-density trend for intermmediate masses.This is because at 10 −2 M Pluto at 20 AU, model 2 is accreting almost pure silicates, like model 3.At 10 −1 M Pluto the models diverge significantly as the larger pebbles are icy in model 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 to ≈ 3 × 10 −3 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 timespan of the hydrodynamical simulation.Also, the Hill radii of these planetesimals are resolved.For 45 AU  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.and at this resolution, the mass at which the Hill radius R H is equal to the length ∆ of a mesh cell is 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 ≈ 160 km, mean density ≈ 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 26 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 poor thermal conduction than larger objects, trapping heat and melting the interior if the ice mass fraction is low.The melting removes the porosity, making these objects high density.
In our model, significant dispersion in ice mass fraction exists for streaming instability products less massive than 10 −3 M Pluto , so some objects (the more ice-rich) should avoid melting whereas some (the more rockrich) should undergo melting.The scatter for the low mass objects is consistent with Fig. 1.All streaming instability products more massive than 10 −3 M Pluto are ice-rich and should avoid melting.

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 pro-  21) and ( 22)).Actual KBO data is overplotted as green stars.The best accreted pebbles at each mass are still too icy to lead to significant density increase through silicate accretion.The most massive planets formed within 5 Myr are of Pluto mass, but their composition is almost pure water ice.
Middle row: Same as upper row, but for the bimodal density distribution (model 2, Eqs. ( 23) and ( 24)).As the intermmediate-mass objects (10 −2 − 10 −1 M Pluto ) accrete silicate pebbles, the mass-density trend is reproduced.The flattening at higher mass is also reproduced, matching Haumea, Pluto, and Triton.The model does not reproduce the density of Eris.
Lower row: Same as the other rows, but for constant (silicate) pebble density (model 3).The intermmediate-mass objects (10 −2 − 10 −1 M 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.
Having a density of 0.64 g cm −3 at 20 K, methanol ice could also explain the low densities of the small KBOs, if a significant fraction of their bulk is of that material.Indeed, near an ice line pebble growth is boosted through deposition and nucleation of vapor onto bare silicate and already ice-covered pebbles (see Ros et al. 2019, in the context of water ice).This possibility is consistent with the lack of observed water ice spectral features on small KBOs (Barucci et al. 2011;Grundy et al. 2020).The CO ice line is also favorable to the formation of the ice giants, not only owing to the excess of carbon and lack of nitrogen found in the ice giants (Ali-Dib et al. 2014) but also because ice giants are thought to contain > 10% methane by mass, which would occur if the ice giants formed between the CO and N 2 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).

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 26 Al with a half-life of 7 × 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 evolu- tion of KBOs, when a large fraction of 26 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 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).

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, throw- Ice Volume Fraction (%) Figure 12.Comparison of accretion rates for the bimodal distribution at different radii.The window at which silicate accretion is most effective shifts to higher masses as we move outwards in the disk.This provides further evidence for 20 AU being a favorable location, as silicate accretion is most effective around 10 −2 M Pluto , right when the densities of Kuiper belt objects begin to increase.
ing 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 postinstability.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 ≈25 AU (Johansen et al. 2015).Beyond this distances, planetesimals do not grow up to Triton or Pluto sizes, even by polydisperse pebble accretion Lyra et al. (2023), as we have demonstrated.Scattering small planetesimals is not efficient to drive further migration, so Neptune's current location coincides with where pebble accretion stops being efficient at growing large objects.

LIMITATIONS AND FUTURE WORK
The model presented is an exploratory proof-ofconcept 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 coagula-tion rate are similar, which a growth time estimate indeed shows they are.The growth time of dust grains, up to a factor of order unity, is t grow ∼ ε −1 T orb (Birnstiel et al. 2012;Lambrechts & Johansen 2014;Lorek & Johansen 2022), where T orb is the orbital period.The disk starts with ε ∼ 10 −2 of dust, of interstellar (submicron) size.Assuming an orbital period T orb ∼ 100 years yields t grow ∼ 10 4 years.The dust thus coagulates very quickly into pebbles.The steady state, however, is top-heavy, leaving about ε ∼ 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 ∼1 Myr, which is similar to the photodesorption rate.
Despite this assurance, the argument remains mostly qualitative, and thus the conversion between radius and composition we used is arbitrary.A future model should include coagulation, drift, fragmentation, turbulence, UV and cosmic ray desorption, calculating the grain size vs composition relation from first principles.Such a model would also have the benefit of constraining the transition in the bimodal model, which we arbitrarily placed at 1 mm.
Secondly, it is an inconsistency that we used a twospecies 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 resupply.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 26 Al and the impact that radioactive heating has on the porosity of KBOs.

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 −3 .We run a hydrodynamic simulation to model planetesimal formation via the streaming instability of ices and silicates, where silicates are small grains with short friction times, and ices are large grains with long friction times.We find that planetesimals formed under these conditions are icy and are low in mass M p ≲ 10 −2 M Pluto , effectively reproducing the densities observed in the low mass Kuiper belt objects.We also find a correlation between ice mass fraction and planetesimal mass, albeit with significant dispersion.
We then use these planetesimals formed by the streaming instability as starting masses for the recently derived polydisperse pebble accretion model (Lyra et al. 2023).We consider three different population of pebbles: the first is a power law distribution where the density of pebbles follow the power law in Eqs. ( 21) and ( 22), with the largest grain being pure ice and the smallest grain being pure silicate.The second is a bimodal distribution, given by Eqs. ( 23) and ( 24), where grains smaller than 1 mm are pure silicates, and larger grains follow a power law, with the largest grain being pure ice.The last model assumes all pebbles are silicates.
We find the power law model does not reproduce the densities for high-mass KBOs.The model with purely silicate pebbles is able to reproduce the high density of Eris, yet fails to reproduce the densities of Haumea, Pluto, and Triton.The bimodal distribution, however, is capable of reproducing the densities of the dwarf planets but underestimates the density of Eris.Thus we speculate that Eris could have formed from a rockrich 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 26 Al.In our model, the first planetesimals are icy, and 26 Al is mostly incorporated in the long phase of silicate pebble accretion, when most 26 Al 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 multispecies streaming instability could result in planetesimals of nearly uniform composition and that a poly- St k ≡ Ω 0 τ k . (A3) The quantities u(x k ) and ρ g (x k ) refer to the gas velocity u(x) and density ρ g (x), defined on the Eulerian mesh positions x, interpolated to the position x k of particle k.The interpolation for a quantity ψ is done as described in Youdin & Johansen ( 2007) where the sum is done over mesh positions.The mesh weight kernel W is chosen as the Triangular Shaped Cloud algorithm (Hockney & Eastwood 1988), where the nearest 3 cells to the particle in each direction (27 cells in total) contribute to the interpolation.Incidentally, the pebble density needed for Eq. ( 5) is calculated on the mesh cells according to and once the selfgravity potential Φ(x) is calculated and the gradient ∇Φ(x) taken on the mesh, the interpolation ∇Φ(x k ) to the position of each particle k is done according to Eq. (A4) and added to Eq. ( 4).The drag force backreaction f g onto a mesh cell centered at x is calculated in a momentum-conserving way, via Newton's 3rd law (Youdin & Johansen 2007) where V(x) is the mesh cell volume and m k the mass of particle k.In this formulation, the dragforce backreaction for multiple particle species is straightforward, as the sum collects the individual contribution of each particle k.

B. HEATING FROM DECAY OF 26 AL
We estimate here whether the small bodies produced in the streaming instability model would melt as a result of heating from 26 Al.Considering the model from Castillo-Rogez et al. (2009), the volumetric heating rate H is Where ρ is the protoplanet's density, F s is the silicate mass fraction, [ 26 Al] 0 is the initial isotopic abundance of 26 Al in ordinary chondrites, H 0 is the heat production rate per mass, and λ = ln(2)/t 1/2 is the decay constant of 26 Al; t 1/2 is the half-life.We find the energy by integrating in time and multiplying by the volume V
2. THE MODELWe use the PENCIL CODE 1 (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

Figure 2 .
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 µm and mm are pure silicates, and pebbles above this size are increasingly icier with pebbles at one centimeter being pure ice.(Right) The constant distribution, in which pebbles of all grain sizes have a 0% ice volume fraction and a density of ρ • = 3 g cm −3 .

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

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

Figure 7 .
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 ∼1 orbit.
The Importance of the Silicate Window

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

Figure 9 .
Figure9.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.

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

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

Figure 13 .
Figure 13.Body temperature T 1 resulting from heating from decay of 26 Al, as a function of mass fraction of silicates, F s .

Table 1 .
Symbols used in this work.
x , u y , u z ) gas velocity in Cartesian coordinates