Planet Formation by Gas-assisted Accretion of Small Solids

We compute the accretion efficiency of small solids, with radii 1 cm ≤ R s ≤ 10 m, on planets embedded in gaseous disks. Planets have masses 3 ≤ M p ≤ 20 Earth masses (M ⊕) and orbit within 10 au of a solar mass star. Disk thermodynamics is modeled via 3D radiation-hydrodynamics calculations that typically resolve the planetary envelopes. Both icy and rocky solids are considered, explicitly modeling their thermodynamic evolution. The maximum efficiencies of 1 ≤ R s ≤ 100 cm particles are generally ≲10%, whereas 10 m solids tend to accrete efficiently or be segregated beyond the planet’s orbit. A simplified approach is applied to compute the accretion efficiency of small cores, with masses M p ≤ 1 M ⊕ and without envelopes, for which efficiencies are approximately proportional to Mp2/3 . The mass flux of solids, estimated from unperturbed drag-induced drift velocities, provides typical accretion rates dM p /dt ≲ 10−5 M ⊕ yr−1. In representative disk models with an initial gas-to-dust mass ratio of 70–100 and total mass of 0.05–0.06 M ⊙, the solids’ accretion falls below 10−6 M ⊕ yr−1 after 1–1.5 Myr. The derived accretion rates, as functions of time and planet mass, are applied to formation calculations that compute dust opacity self-consistently with the delivery of solids to the envelope. Assuming dust-to-solid coagulation times of ≈0.3 Myr and disk lifetimes of ≈3.5 Myr, heavy-element inventories in the range 3–7 M ⊕ require that ≈90–150 M ⊕ of solids cross the planet’s orbit. The formation calculations encompass a variety of outcomes, from planets a few times M ⊕, predominantly composed of heavy elements, to giant planets. The peak luminosities during the epoch of the solids’ accretion range from ≈10−7 to ≈10−6 L ⊙.


Introduction
Classic theories of planet formation assume that initial growth proceeds by collisions with planetesimals, bodies ranging in size from a fraction to hundreds of kilometers (e.g., Lissauer 1993).Relict populations of these bodies are found in the asteroid and Kuiper belts (e.g., Morbidelli & Nesvorný 2020;Bourdelle de Micas et al. 2022).The path to formation of planetesimals around stars remains mysterious, but conventional wisdom generally assumes that they are assembled from aggregates of much smaller particles, which are collected through some hydrodynamical process (in the circumstellar gas) and eventually clump together under their own gravity (see, e.g., Weiss et al. 2023, and references therein).These smaller particles, whose presence can be inferred from observations of protoplanetary disks (e.g., Natta et al. 2007;Andrews 2020), would form through coagulation of primordial dust grains entrained in the gas.
This scenario for planetesimal formation has led to the suggestion that swarms of "small particles" orbiting a star in a protoplanetary disk may directly supply solid material to emerging planetary embryos and planets, determining their heavyelement inventories (Ormel & Klahr 2010;Lambrechts & Johansen 2012).Over the past decade, such delivery mode of solids has been referred to as "pebble accretion".However, despite the terminology, these solids are unrelated to clastic rocks and are not required to fulfill the geological definition (Gary et al. 1972).In fact, they are not characterized by size or composition, but rather by drag interaction with circumstellar gas, corresponding to Stokes numbers (see below) in the range between ∼ 10 −3 and ≈ 1 (Drazkowska et al. 2023).These numbers quantify the coupling timescale, in units of the orbital time (i.e., the inverse of the Keplerian frequency), between particle and gas dynamics.The Stokes number is proportional to the radius and material density of a particle, and inversely proportional to the gas density.Therefore, according to current terminology, astrophysical pebbles orbiting in a solar-nebula type of disk, between ≈ 1 and ≈ 10 au, can range from submillimeter silicate grains (i.e., astrophysical dust, D' Alessio et al. 2001) to meter-sized ice blocks.At larger distances, as the gas density lowers, smaller particles would be considered as astrophysical pebbles (and the opposite would occur closer to the star).An overview of the typical outcomes from pebble accretion calculations can be found in, e.g., Johansen & Lambrechts (2017), Drazkowska et al. (2023), and references therein.
Herein, we build first-principles numerical models of gas and solids thermodynamics.Since Stokes numbers cannot be constrained in the calculations (because they depend on primitive variables), we consider specific particle sizes (from 1 cm to 10 m) and compositions (SiO 2 and H 2 O), and we will refer to these particles simply as "small solids".Overlap with the domain of astrophysical pebbles varies according to local gas conditions.
The largest difference between planetesimal and small-solid accretion scenarios rests on their dynamics as they evolve in circumstellar gas.On the one hand, drag forces only provide a correction to the Keplerian orbits of the large bodies (Whipple 1973).The drag-induced drifting timescale of 1 km planetesimals is ∼ 10 6 orbital periods (in the 1-10 au region), hence radial mobility due to aerodynamic drag can be largely ignored during planet formation.The accretion rate on a planet is determined by the local surface density of solids and by a characteristic cross-section for collisions, which can be enhanced by gravitational focusing and atmospheric drag (see, e.g., Pollack et al. 1996;Inaba & Ikoma 2003).On the other hand, drag forces strongly affect the motion of small solids (Weidenschilling 1977).Radial drift occurs on much shorter timescales than those of planetesimals.Hence, accretion is a non-local process.The global transport of solids through the disk determines the availability of the supply and the gas thermodynamics at and below the orbital length-scale determines the probability of accretion.The former cannot be evaluated without a global model of the circumstellar gas.The latter requires modeling tidal interactions between the planet and the disk (e.g., Morbidelli & Nesvorny 2012) and gas dynamics in proximity of the planet (e.g., Picogna et al. 2018;Popovas et al. 2018).Therefore, any outcome of a planet formation calculation based on small-solid accretion is bound to depend on both local and remote gas thermodynamics.The applicability of generic formulations of accretion rates, found in the literature, is obviously limited by underlying assumptions (see, e.g., Lambrechts & Johansen 2012;Drazkowska et al. 2023).
Planetesimal accretion proceeds as long as there is material in a region, the "feeding zone", spanning a few Hill radii on either side of the planet's orbit.Once this zone is depleted, accretion is much reduced and limited to the resupply rate into the region (e.g., Pollack et al. 1996;D'Angelo et al. 2014).Small-solid accretion can proceed as long as there is transport of material across the planet's orbit, which can be hindered by the overall depletion of solids, dissipating gas, and disk-planet tidal interactions.
In terms of interior compositions, the heavy-element inventory provided by the two accretion scenarios may differ, since solids are accreted locally in one case and (effectively) remotely in the other.However, if planetesimals form out of small solids, which have drifted considerable distances prior to clumping together, planetesimals too can have (largely or to some extent) non-local compositions.Hence, diversity in compositions may not discriminate between the two scenarios.Differences in heavy-element stratification may also be difficult to ascertain.Clearly, small and large solids may coexist (e.g., Kessler & Alibert 2023), although mass partitioning would depend on the mechanism that converts one population into the other (see, e.g., Drazkowska et al. 2023;Weiss et al. 2023, and references therein).
Herein, we present direct calculations of small-solid accretion, based on three-dimensional radiation-hydrodynamics calculations of disks with embedded planets that resolve the planetary envelope.Models also account for the thermal evolution of the particles.We focus on planet masses between 3 and 20 M ⊕ .These calculations provide the accretion efficiency (or probability) of the particles on the planets at different orbital locations (1 to 10 au) and disk conditions.The accretion efficiency is the fraction of the local accretion rate of solids in-tercepted by the planet.The results are complemented with those obtained from a simpler approach, adopted to compute the accretion efficiency of bare cores (i.e., without a substantial gaseous envelope) less massive than one Earth mass.The global transport of small solids, computed from drag-induced drift velocities in unperturbed disks, is combined with the accretion efficiencies to provide the accretion rates of the planets.These rates are then applied to actual formation and structure calculations, which compute dust opacity in the planet's envelope self-consistently with the delivery of solids.These latter calculations are intended as demonstrations of the framework provided herein and do not target any specific planet.
In the following, the orbital dynamics of small solids in unperturbed and perturbed disks is described in Section 2. The radiation-hydrodynamics models and accretion efficiencies of envelope-bearing planets are presented in Section 3; the accretion efficiencies of small, bare-core planets are discussed in Section 4. The formation calculations are presented in Section 5. Finally, our conclusions are summarized in Section 6.

Dynamics in Unperturbed Disks
Small solids orbiting a star in a smooth unperturbed disk may experience strong drag forces exerted by the gas, and thus undergo rapid orbital decay.Assuming that the solids' mass is locally much smaller than the gas mass, the time required for the velocity of a solid particle, u s , to converge toward the gas velocity, u g , is where a D is the drag acceleration.The time τ s is referred to as the "stopping time" of the particle (Whipple 1973;Weidenschilling 1977).Indicating with ρ g and ρ s respectively the gas and solid density, R s the particle radius, and a s its orbital distance, Equation (1) becomes (Whipple 1973;Weidenschilling 1977) The drag coefficient, C D , is a function of the thermodynamic properties of both gas and solids.In the free molecular flow regime (i.e., for very small particles), The rotation rate of the gas, Ω g , is generally different from the Keplerian rate, Ω K , and depends on the radial gradient of the gas pressure, P g , typically a function of gas density and temperature, T g .By requiring conservation of the gas linear momentum in the disk's radial direction, r, it can be approximated as Φ being the gravitational potential in the disk.Neglecting the self-gravity of gas and solids, ∂Φ/∂r = rΩ 2 K and Equation (3) applied to the midplane reduces to where H g is the disk's pressure scale-height at the midplane.The quantity η depends on the normalized gradients ∂ ln ρ g /∂ ln r and ∂ ln T g /∂ ln r (see, e.g., Takeuchi & Lin 2002;Tanaka et al. 2002), and is of order unity.Hereafter, η is included in the ratio H g /a s .For typical protoplanetary disks, (Ω K − Ω g )/Ω K ≈ (H g /a s ) 2 /2 is on the order of a percent or less (e.g., D'Alessio et al. 1998).
The midplane (non-transient) velocity of a particle subject to drag can be approximated to (e.g., Chiang & Youdin 2010, and references therein) (5) φ s being the azimuthal angle of the particle.Equations ( 5) and ( 6) assume a negligible radial component of u g .Quantity τ s Ω K is referred to as the Stokes number of the particle.For 2 τ s and the radial displacement converges to that of the gas.
Using the coefficient C D reported in D'Angelo & Podolak (2015, hereafter, DP15), which depends on |u g − u s | and hence on τ s , Equation ( 2) can be inverted to provide a function τ s = τ s (R s ), plotted in Figure 1 (top panel).Through Equation (5), one can then recover the timescale for drag-induced orbital decay, τ D (middle panel).For reference, C D is plotted in the bottom panel.The results in Figure 1 are based on disk thermodynamic conditions relevant to the radiationhydrodynamics calculations described in Section 3, also detailed in Table 1.Under such conditions, τ s Ω K ≈ 1 when R s is a few times 10 to ≈ 100 cm.Smaller particles may satisfy the condition τ s Ω K ≈ 1 in a low density disk.The timescale τ D has a minimum of a few hundred orbital periods or less, and diverges as R s → 0 (τ s Ω K → 0) because the drift velocity | ȧs | ∝ Ω 2 K τ s a s approaches zero.In this limit, however, ȧs converges to the radial velocity of the gas, u r g , which cannot be neglected.In such cases, a corrected version of Equation ( 5) is (see, e.g., Takeuchi & Lin 2002) ȧs In the classic theory of accretion disks (e.g., Lynden-Bell & Pringle 1974;Pringle 1981), in which ν g and Σ g are respectively the kinematic viscosity and surface density of the gas.Writing the kinematic viscosity as ν g = α g H 2 g Ω K (Shakura & Sunyaev 1973), if ∂ ln (ν g Σ g )/∂ ln r ≪ 1 then the bracket in Equation ( 8) is ≈ 1, and the correction term in Equation ( 7) is negligible as long as τ s Ω K ≫ α g .Indicating with Σ s the solids' surface density, the mass accretion rate of solids through an unperturbed disk is dM s /dt = The different lines refer to different disk locations and/or gas densities (as indicated in g cm −3 ) and temperatures (see also Table 1).
−2πΣ s a s ȧs .Applying Equation ( 7), the accretion rate becomes Neglecting the term in u r g , maximum accretion occurs for At this rate, a mass of solids equal to πa 2 s Σ s would be transported through the radius r = a s in a timescale (a s /H g ) 2 /Ω K , or ∼ 10 3 years at r = 10 au.In the decoupled regime, In the well-coupled regime, if τ s Ω K α g , the term in u r g cannot be neglected and Equation (9) implies possible stalling or outward transport of solids in expanding disk regions (where u r g > 0, Lynden-Bell & Pringle 1974) or in transition regions with steep density gradients (see Equation ( 8)).
Equation ( 9) is plotted in Figure 2, for a radially constant value of Σ s = 1 g cm −2 and u r g ≈ −(3/2)ν g /r, in which ν g is constant and corresponds to α g ≈ 0.005 at 1 au (H g ∝ r 9/7 ).The effect of u r g can be seen when τ s Ω K α g (α g is roughly proportional to 1/r).The curves in the figure should only be interpreted as the value at r for the imposed Σ s .In fact, if ∂ Ṁs /∂r = 0, mass conservation requires that where Λ includes source/sink terms to account for, e.g., the ongoing coagulation of solids from dust or fragmentation and erosion.
A solution of Equation ( 10) is plotted in Figure 3, requiring that τ s Ω K = 1 in Equation ( 9).The initial conditions are based on the nebula model of Chiang & Youdin (2010), with a total gas mass of 0.055 M ⊙ (Σ g ∝ r −3/2 and gas temperature T g ∝ r −3/7 ) and a total solids' mass of ≈ 183 M ⊕ (for an initial gas-to-solids mass ratio of 100).The conversion of dust into the larger particles is assumed to have occurred at time t = 0, with no further growth or fragmentation (i.e., Λ = 0).Over the first ≈ 10 5 years, the disk loses ≈ 120 M ⊕ of solids to the star, with the remainder removed during the next ≈ 1.2 × 10 5 years.For comparison, the condition τ s Ω K = 0.1 removes ≈ 90 M ⊕ in ≈ 3 × 10 5 years and the entire reservoir of solids within ≈ 10 6 years, whereas the condition τ s Ω K = 0.01 causes a depletion of ≈ 50 M ⊕ in ≈ 10 6 years.Clearly, these numbers are determined by the initial mass and distribution of solids.Provided that H g remains roughly constant in time, the condition τ s Ω K = constant allows one to ignore the disk's gas evolution.However, due to the changing thermodynamic properties of the gas, this condition does not characterize any single solid's size, but it rather corresponds to a changing R s , in space and time.
The evolution of solids of a particular size can be obtained by solving Equation (10) for a fixed R s .In this case, τ s Ω K is a function of both t and r, and Equation (2) must be inverted at any time and distance in order to self-consistently determine ȧs and τ s .The evolution of the disk's gas cannot be neglected in this case.Thus, it is assumed that gas is depleted via accretion on the star and through photo-evaporation by hard radiation from the star.The gas kinematic viscosity is taken as ν g ∝ 1/Σ g (α g ≈ 10 −3 at 1 au), so that d(ν g Σ g )/dr = 0 and Equation (8) reduces to u r g = −(3/2)ν g /r.For simplicity, the time evolution of Σ g is obtained by rescaling the initial condition (the same as in Figure 3) according to the current gas mass (∂ ln Σ g /∂ ln r is constant in time).The disk's gas is entirely depleted in ≈ 3.5 Myr (e.g., Weiss et al. 2021).
Figure 4 illustrates the evolution for rocky solids of 1 (left), 10 (center), and 100 cm (right) radius.The initial inventory of solids, ≈ 183 M ⊕ , is entirely removed within ≈ 9 × 10 5 years for the case shown in the left panels.The reason is that maximal accretion (τ s Ω K ≈ 1) is initially sustained at distances of several tens of astronomical units, rapidly transferring large amounts of solids toward the star (see the top left and middle left panels).Consequently, two-thirds of the initial mass is removed in ≈ 3 × 10 5 years.For particles with R s = 10 and 100 cm, respectively, ≈ 60 and ≈ 110 M ⊕ are left at 1 Myr, but 1 M ⊕ remains inside 20 au.At t ≈ 5 × 10 5 years, dM s /dt inside 100 au ranges from a few to several 10 −5 M ⊕ yr −1 , dropping below 10 −5 M ⊕ yr −1 after 1 Myr.
An experiment conducted with R s = 10 m rocks indicates that over 75% of the initial solids' mass survives after 3 Myr (because τ s Ω K ≫ 1), but only an amount of ≈ 5 M ⊕ orbits within 10 au, where dM s /dt 10 −5 M ⊕ yr −1 after 10 5 years.Beyond 20 au, dM s /dt never exceeds ≈ 10 −5 M ⊕ yr −1 during the disk's evolution.
The bottom panels of Figure 4 indicate that for 0.1 < r < 100 au, τ s Ω K can vary by up to three orders of magnitude for R s = 1 and 10 cm particles and by over an order of magnitude for R s = 1 m (and 10 m) bodies.The largest particles also display a significant increase of τ s , at any given r, as the disk's gas dissipates.Clearly, the condition τ s Ω K = constant is only indicative of a particle size for a given gas state.In these examples, τ s Ω K = 1 corresponds to R s ≈ 1 cm particles between r ≈ 50 and ≈ 80 au, and to R s ≈ 10 cm particles between r ≈ 10 and ≈ 20 au.Larger particles (in the ≈ 10-100 cm range) would satisfy this condition inside 10 au.
Figure 5 shows dM s /dt versus time at r = 1 (top), 5 (middle), and 10 au (bottom), for the calculations presented in Figure 4 (solid lines).For the first few 10 5 years, dM s /dt can be 10 −4 M ⊕ yr −1 , but it becomes 10 −5 M ⊕ yr −1 by ≈ 1 Myr.The same calculations were executed applying the source term Λ(t, r) = 0.9 Σ s (0, r)/τ c in Equation (10), with a "coagulation time" τ c = 3 × 10 5 yr (Voelkel et al. 2020) and Λ = 0 for t > τ c .At time t = 0, 10% of the initial solids' mass is already in the form of particles of radius R s .These latter results are illustrated in Figure 5 as dashed lines.Compared to the solid curves, peaks in dM s /dt are lower and broader, but differences become small after ≈ 0.5 Myr.The results presented in Figures 3-5 are approximate but representative of the evolution of the solids' density and accretion rates in an unperturbed disk.

Dynamics in Perturbed Disks
In the presence of a perturbing body, e.g., a planet, Equation (9) can be interpreted as the rate at which solids can be supplied to the proximity of its orbit, thus determining the maximum rate of the solids' accretion on the planet.Assuming that the thickness of the disk of solids, H s , is infinitesimal, then the accretion rate on a small planet is essentially two-dimensional and can be approximated by in which M p and a p are the planet mass and orbital radius, respectively.This expression assumes that the mass flux of solids in the radial direction can be approximated as azimuthindependent, so the fraction of the integrated flux intercepted by the planet is given by the ratio of the planet's "effective" size to the orbit circumference.This approximation is useful when working in one-dimensional, radial disks.Corrections due to azimuthal variations of the solids' transport can be accounted for in the definition of the planet's effective radius for accretion.In fact, the radius R eff can be much larger than the planet's physical radius, R p (the two-dimensional assumption implies that R eff > H s ).Applying a formalism developed for the accretion of planetesimals in the two-body problem framework (see, e.g., Lissauer 1987), the effective radius for the capture of solids can be written as , where u esc = 2GM p /R p is the escape velocity from the planet surface and u rel is the relative velocity (in magnitude) between the planet and the solids.This relation originates from the conservation of relative energy and angular momentum of the approaching particle.If τ s Ω K ≫ α g , the relative velocity can be obtained directly from Equations ( 5) and ( 6), assuming that the planet moves on a circular orbit.
The efficiency, or probability, of the accretion of solids onto a planet is defined as where Ṁs is the accretion rate of solids toward the planet, in the proximity of its orbit (and Ṁp ≤ Ṁs by assumption).Inserting the effective radius given above in Equation ( 11), we find in which M ⋆ is the stellar mass.For τ s Ω K 1, Equation ( 14) scales with planet mass and stopping time as the accretion probability estimated by Kary et al. (1993) in the "nongravitating" planet limit, , where R H is the Hill radius of the planet.Maximum values of the large-scale transport of solids, i.e., dM s /dt in Equation ( 9), occur when For a Moon-sized planetary embryo orbiting at 5 au, the efficiency would be ≈ 10 −4 .Therefore, the values of dM s /dt reported in Figure 5 would result in dM p /dt < 10 −7 M ⊕ yr −1 , and in a mass-doubling timescale exceeding 10 5 years.
Alternatively, the effective radius for the capture of solids can be estimated by equating the gravitational and relative kinetic energies (Lambrechts & Johansen 2012), a condition required for escape, resulting in R eff ≈ 2GM p /u 2 rel .By using Equation ( 12), this approximation yields an efficiency In this case, a Moon-sized embryo orbiting at 5 au would have ζ ≈ 10 −2 at maximum values of dM s /dt (τ s Ω K ≈ 1), hence dM p /dt 10 −5 M ⊕ yr −1 , much larger than the previous assessment.However, both approximations for the effective radius neglect stellar gravity, implying that in neither case can The accretion efficiency of a Moon-mass body would therefore be limited to ≈ 10 −3 .The remainder of the difference between the estimates provided by Equations ( 14) and ( 15) can be resolved by applying a more appropriate evaluation of u rel close to the planet, as explained below.The trajectory of the smallest particle is initiated as the others, but a smaller portion is plotted.The gray lines represent streamlines of the perturbed gas.The unperturbed disk properties are as in Figure 1 at 5 au (see also Table 1).
If R eff H s , the two-dimensional approximation for the accretion of solids breaks down and the right-hand side of Equation (11), hence the efficiency ζ (see Equation ( 13)), must be multiplied by (π/4)(R eff /H s ) to account for the vertical distribution of the solids.In the small-particle regime, the thickness H s can be affected by turbulent stirring.If the kinematic viscosity ν g is written in terms of α g , then (Dubrulle et al. 1995) where the numerical factor depends on the nature of the turbulence in the gas (and may be somewhat different from 2, see the discussion in Dubrulle et al. 1995).Typical values of the turbulence parameter α g ∼ 10 −4 -10 −2 would result in H s /H g ∼ 0.01-0.1 for τ s Ω K ≈ 1 (see also Lambrechts & Johansen 2012).Since the gas temperatures in planet-forming regions correspond to H g /a p ∼ 0.01-0.1 (e.g., D 'Alessio et al. 1998), the relative thickness of the solid layer would be H s /a p ∼ 10 −4 -10 −2 .Thus, in the two-dimensional accretion regime (i.e., R eff > H s ), ζ > (1/π)H s /a p , and thus 10 −5 -10 −3 , for τ s Ω K ≈ 1.
Equations ( 14) and ( 15) neglect the gravitational perturbations exerted by the planetary body, which alter the dynamics of both solids and gas, i.e., u rel is different from the unperturbed expression derived from Equations ( 5) and ( 6).An approximation to the perturbed rotation rate of the gas is still given by Equation (3), in which Φ must comprise the contributions due to the planet, including non-inertial terms (if applicable), i.e., where r and r p are the generic and planet position vectors.Neglecting contributions from gas pressure variations along the azimuthal direction around the star, the perturbed radial component of the gas velocity is (Ogilvie & Lubow 2006) in which φ is the azimuthal angle.The radial velocity arising from global transport, Equation ( 8), should be added to the right-hand side of Equation ( 18), if relevant.Given the importance of gas drag, corrections to the gas flow introduced by Equations ( 3), ( 17) and ( 18) can have non-trivial effects on the dynamics of small solids.The perturbed gas velocity can be used to integrate the trajectories of solid particles through the disk (see, e.g., DP15).Some examples are shown in Figure 6, in a frame co-rotating with the planet, for R s = 1 (green line), 10 (orange line), and 100 cm (purple line) rocky particle drifting through the gas perturbed by 5 (top) and 15 M ⊕ (bottom) planets on circular orbits.The gray lines represent perturbed gas streamlines.The particle motion is altered by the planet's gravity and by drag forces due to the perturbed gas velocity field, invalidating the approximations applied in deriving Equations ( 5) and (6).Those equations remain valid only sufficiently far from the orbit of the planet.Note that this approach does not account for planet-induced tidal perturbations on the gas density, and hence pressure, which can also affect particle dynamics, as discussed later, and are expected to become increasingly important as M p grows.
Figure 7 illustrates results from similar experiments, at lower masses: M p = 0.01 (top), 0.1 (middle), and 1 M ⊕ (bottom) at a p = 5 au.The particle radius is R s = 1 (dashed lines) and 30 cm (solid lines), corresponding respectively to τ s Ω K ≈ 0.01 and ≈ 1 for the disk conditions at 5 au in Figure 1.Gas streamlines are also plotted, along with the Roche lobe contours in the disk midplane.The calculation of the Roche lobe contours neglects modifications induced by gas drag (Murray 1994).The numerical results indicate that the estimate of R eff used in Equation ( 14) may be a reasonable approximation in these cases.Although the accretion stream entering the Roche lobe is asymmetric with respect to the planet (see the orange trajectories), an average size of R eff can be estimated as a few tenths of R H .An outcome similar to that of the R s = 30 cm particles is obtained from the trajectory of the R s = 100 cm particles (not shown).In these experiments, the particle trajectories are integrated until they impact the planet, assumed to be a condensed object.As argued below, an atmosphere could extend at most over some fraction of the Bondi Figure 7. Trajectories of Rs = 1 (dashed lines) and 30 cm (solid lines) rocky particles in the proximity of 0.01 (top), 0.1 (middle), and 1 M⊕ (bottom) planets, orbiting at ap = 5 au on circular orbits.The trajectories are represented in a frame co-rotating with the planet.The gray lines represent gas streamlines.The planet's Roche lobe is also shown.The unperturbed disk properties are as in Figure 1 at 5 au (τsΩK ≈ 1 for Rs = 30 cm and ≈ 10 −2 for Rs = 1 cm).radius, 0.2 R H in the M p = 1 M ⊕ case (and much less in the other cases), a length R eff .
It must be pointed out that in the presence of a perturbing body, the unperturbed relative velocity provided by Equation (12), which is independent of the distance r from the perturbing object and of its mass M p , is appropriate for use only far away from the perturber.However, both Equations ( 14) and ( 15) neglect three-body effects, and therefore the relative velocity should be sampled where the gravity field of the planetary body dominates, that is, where gravitational perturbations on u rel are non-negligible.If τ s Ω K ≪ 1, the perturbed relative velocity can be estimated from Equations (3), ( 17) and ( 18), and it depends on M p , r, and the direction of approach.At a distance r ≈ R H , its magnitude is roughly ∝ R H and, for Moon-mass bodies, it is several times as large as the unperturbed velocity predicted by Equation ( 12).For finite (non-zero) values of the particle stopping time, analytical derivations of the perturbed velocity are more involved, but it can be easily evaluated through numerical integration of the trajectories (as in Figure 7).For τ s Ω K ≈ 1, such experiments indicate that particles approaching a Moon-mass body can achieve relative velocities many times (up to a factor of ten) as large as the right-hand side of Equation ( 12) at distances r R H . Applying this correction, the efficiencies predicted by Equations ( 14) and (15) become similar.

Simple Estimates of Accretion Efficiencies
The approach used for the trajectory integration in Figures 6  and 7 is applied to conduct experiments on the accretion of particles, with radii of 1, 10, 100, and 1000 cm, on M p = 3 and 5 M ⊕ planets orbiting at a p = 1, 5, and 10 au.These experiments provide a statistical estimate of the accretion efficiency ζ for the simplified gas flow given by Equations ( 3), (17), and ( 18).The gas density is axisymmetric around the star and proportional to a power of r.A summary of the results is displayed in Figure 8. Histograms show the distributions of the closest approach distance along the particles' trajectories for the M p = 3 M ⊕ planet at a p = 1 au (top panel) and for the M p = 5 planets at a p = 5 and 10 au (middle and bottom panels, respectively; see also the figure's caption).The gas densities are those quoted in Figure 1 (see also Table 1), according to the planet's orbital radius (at a p = 1 au, Σ g ≈ 1100 g cm −2 ).All trajectories begin as random circles at radial distances beyond the 2:3 mean-motion resonance with the planet, out to a distance somewhat beyond the 1:2 mean-motion resonance (in the orbital plane of the planet).These experiments assume equal masses of icy and rocky particles, in each size bin, at a p = 5 and 10 au, and only rocks at a p = 1 au.Accretion is occurs if a particle enters the planetary envelope, r < R p and R p is taken as the Bondi radius, Equation ( 19).For a solid to be permanently captured, its relative velocity in the envelope must not exceed the escape velocity, u esc .These experiments do not account for such requirement, which is generally irrelevant for the small particles considered here (as confirmed by the calculations discussed in Section 3).
The closest approach distance from the planet in Figure 8 is in units of the Hill radius, R H .For a given particle radius, the histogram bins are normalized to indicate the fraction of the initial mass of the solids with radius R s considered in the experiment.The bins are then overlaid in order of increasing solid size.The vertical dashed line indicates the distance r = R B , the assumed envelope radius.Consistent with the results of Figure 7, the histograms show that particles of all sizes can enter the planet's Hill sphere without necessarily being captured.The pile-up in the bin interior to or overlapping R p is caused by particles flagged as accreted.Since each color represents a normalized mass, the innermost bin in each panel provides a measure of the efficiency ζ, estimated here as the accreted mass divided by the mass of the solids approaching and interacting with the planet.
Particles with an initial τ s Ω K ∼ 1 have R s ≈ 10-100 cm.For both planet masses, the efficiency ranges from a few to several percent.For 1 cm particles (τ s Ω K 10 −2 initially), ζ ranges from a few to ≈ 15%.The largest, 1000 cm particles (τ s Ω K > 10 2 initially) can accrete efficiently: ζ > 50% in with Ṁs ∼ 10 −5 M ⊕ yr −1 over the first 1 Myr, would result in Ṁp of order 10 −7 -10 −6 M ⊕ yr −1 .During the first 1 Myr, R s = 1000 cm solids also have Ṁs ∼ 10 −5 M ⊕ yr −1 , which would generate values of Ṁp up to ∼ 10 −5 M ⊕ yr −1 .It is important to notice that high values of ζ for a particular solid size imply an efficient removal of these solids, which then cannot significantly contribute to the growth of other planets orbiting closer to the star.
A number of effects are neglected in these experiments.Gas density perturbations driven by the planet's tidal field can impact the drag forces exerted on the particles, possibly causing a size-dependent segregation effect.Additionally, the complex three-dimensional flow circulation in the proximity of the planet can alter the balance of forces and change the outcome of the capture process.Moreover, for the icy solids, variations of the gas temperature close to the planet can promote ablation and vaporize some fraction of the local mass of solids.These effects, which require more realistic disk-planet interaction models and more sophisticated calculations of the solids' evolution in disks, are included in the simulations discussed below.

Gas Thermodynamics
During the relatively rapid motion of small particles through the disk, the gas component of the disk can be assumed to be in a quasi-steady state, so that its thermodynamic fields, e.g., velocity u g , temperature T g , and density ρ g , remain nearly constant in time when described in a reference frame co-rotating with an embedded planet on a circular orbit.This assumption is based on the fact that global disk evolution is mainly driven by viscous transport during the stages considered here.This study uses the thermodynamic steady-state fields (u g , ρ g , T g ) from the models of D'Angelo & Bodenheimer (2013, hereafter, DB13), who performed global three-dimensional (3D) radiation-hydrodynamics (RHD) calculations of planets embedded in disks.They modeled planets of mass M p = 5, 10, and 15 M ⊕ on circular orbits at a p = 5 and 10 au from a solarmass star.In all cases, the envelope mass was a small fraction of M p .Table 1 lists some properties of these models.
Briefly, DB13 adopted a spherical polar discretization of the disk, with coordinates {r, θ, φ}, solving the Navier-Stokes equations for a compressible and viscous fluid in a reference frame rotating about the star at the same angular velocity as the planet's orbital frequency.They applied a numerical method that is effectively second-order accurate in both space and time (see, e.g., Boss & Myhill 1992).The disks extended in radius from 0.5 a p to 2 a p , about 15 • in the meridional direction and 2π radians in azimuth around the star.The energy equation accounted for the advection of gas and radiation energy, for the work done by gas and radiation pressure, and for viscous dissipation, radiation transport, and energy released by accreted solids.Radiative transfer was approximated via fluxlimited diffusion (Levermore & Pomraning 1981;Castor 2007) and solved by means of an implicit algorithm second-order accurate in both space and time (an approach recently validated a Ratio measured for models at 5 and 10 au, assumed at 1 au.
b Midplane gas temperature and density at r = ap, averaged around the star.
by Bailey et al. 2023).The gas equation of state was one of a mixture of hydrogen and helium in solar proportions that accounted for the ionization of atomic species, the dissociation of molecular hydrogen, and for the translational, rotational, and vibrational energy states of H 2 .
Once a thermodynamic quasi-steady state has been achieved, as in the models of DB13, the disk evolution proceeds on the viscous timescale, ∼ r 2 /ν g .The calculations adopted a turbulent viscosity parameter α g ≈ 10 −3 , corresponding to a viscous timescale of several 10 4 orbital periods at r = a p .This level of turbulence is in accordance with observational estimates (e.g., Flaherty et al. 2018).Toward the end of a disk's lifetime, when the gas density is very low, evolution is mainly driven by stellar photo-evaporation and the gas is cleared inside-out (e.g., Gorti et al. 2009;Ercolano & Pascucci 2017).
For the purposes of this study, an important aspect of DB13's models is that they are global (disks extend many astronomical units in radius) and resolve the actual envelopes of the planetary cores.They did so by applying hierarchies of nested grids with finest resolutions comparable to the condensed core radii.In fact, a detailed comparison was performed between 3D envelopes and one-dimensional, spherically symmetric envelopes obtained from planet structure and evolution calculations.The high resolution of the models allows the trajectories of the solids to be integrated until they enter the actual planetary envelopes.
Models with planets at a p = 5 and 10 au are complemented by additional RHD calculations of M p = 3 and 20 M ⊕ planets orbiting at a p = 1 au from a 1 M ⊙ star.In these cases, however, the gas density is varied to mimic different epochs of the disk's evolution.Values of Σ g at r = a p vary from ≈ 20 to ≈ 1100 g cm −2 (see Table 1).The steady-state gas density and temperature in the disk midplanes of a p = 1 au models are illustrated in Figure 9. Contrary to the calculations of DB13, these models do not resolve in detail possible gaseous envelopes bound to the cores.Therefore, the planetary cores are assumed to bear an envelope whose radius is the smaller of R B and R H .Under this assumption, R p in these six models ranges from two to ten times the linear resolution of the grid.A close-up around the planets is illustrated in Figure 10.As a reference, R H /a p ≈ 0.015 and 0.027 in the top and bottom panels, respectively.
When the envelope volume is determined by an energy balance rather than by gravity, R p is comparable to the Bondi radius (see Table 1).This distance is found by equating the thermal velocity of the gas around the planetary core and its escape velocity from the core, where G, m H and k B are standard physical constants, µ g is the mean molecular weight of the gas mixture, and T g is the average gas temperature at the disk midplane at the radial distance r = a p .Equation ( 19) is valid when the energy of the gas is dominated by thermal energy.However, in the 20 M ⊕ planet models, R B is similar (warmest disk model) or exceeds R H . Therefore, in those three cases, it is assumed that R p = 0.7R H , the mean volume radius of the Roche lobe (Eggleton 1983).The planetary radii in the models, compared to R B and R H , are reported in Table 1 (see also Kuwahara & Kurokawa 2024, for a comparison).The approach of simulating the evolution of the solids in the gas steady-state fields is dictated by the computational cost of these RHD calculations, which can hardly be executed along with the thermodynamic evolution of solids for the required timescales.Nonetheless, the validity of this approach is tested in Appendix A, where results from two RHD calculations are compared.In the first calculation, a population of small solids evolves in the 3D steady-state fields (u g , ρ g , T g ) of a disk, as described in this section.In the second, the same population evolves together with the disk's gas, starting from the steadystate fields used in the first calculation.As discussed in Appendix A, the outcomes of the two approaches are statistically consistent, with relative differences in accretion rates 10% and typically hovering around a few percent.
Figure 9. Steady-state gas density (color scale) and temperature (contour lines) at the disk's midplane, perturbed by a (circular orbit) 3 M⊕ planet located at r/ap = 1 and φ = φp.The density is in units of grams per cubic centimeter and the temperature is in kelvin degrees.The azimuthally averaged surface density at r = ap is Σg ≈ 1100 (top), ≈ 220 (middle), and ≈ 20 g cm −2 (bottom).See also Table 1.

Thermodynamics of Solids
The physical model for the thermodynamic evolution of the solids is described in DP15.Briefly, the equations of motion of a solid particle are written in terms of its linear and absolute angular momenta per unit mass.The particle is subjected to gravitational forces by the star and the planet, to non-inertial forces, to aerodynamic drag, and to the effect of mass loss due to ablation.The particle temperature, T s , is determined through an energy equation that takes into account the work done by gas drag, the energy absorbed from the radiation field of the ambient gas, and black body emission from the particle's surface.The energy equation also accounts for the energy removed/supplied during phase transitions of the parti- cle's outer layers.Details on the numerical methods and tests can be found in DP15.
Two types of material are considered here: ice (H 2 O) and rock (SiO 2 ).Particles can lose mass by ablation.In this study, the saturated vapor pressure, P v , of SiO 2 is upgraded to that reported by Melosh (2007): where the critical pressure is P cr = 1.89 × 10 9 dyne cm −2 , the critical temperature is T cr = 5398 K and A = 10.675270434.
The temperature of the solid T s is defined in DP15.Ablation in the ambient disk's gas is important for ice (or an admixture of ice and rock), but much less for rock owing to low gas temperatures.Ablation of rock becomes significant only inside the planetary envelopes (see Appendix B), when particles are already accreted.

Distributions and Segregation of Solids
Four size bins are initially populated with 10 5 particles each, with initial radii R s = 1, 10, 100, and 1000 cm.The initial Figure 11.Histograms of the particles' semi-major axes (left), eccentricities (center), and inclinations (right) as they approach and cross the planet's orbit (ap = 5 au and Mp = 5 M⊕).The top and bottom rows refer to icy and rocky solids, respectively.The gray histograms include all particles (Rs ≤ 10 m); the red histograms include particles with radii Rs ≤ 100 cm; the green histograms include particles with radii Rs ≤ 10 cm; and the blue histograms include particles with radii Rs ≤ 1 cm.
(circular) orbits of the solids are randomly distributed between a s = 1.35 and 1.65 a p , with orbital inclinations varying between i s = 0 and 1.7 degrees and random longitudes of the ascending node.Separate calculations are carried out for icy and rocky particles.
The normalized stopping time at deployment differs little between cases at a p = 5 and 10 au, ranging from τ s Ω K ∼ 10 −2 (R s ≈ 1 cm) to ∼ 10 2 (R s ≈ 1000 cm).The range is wider at a p = 1 au (see Figure 1).Figure 11 shows the distributions of the orbital properties of the solids approaching and crossing the orbit of an M p = 5 M ⊕ planet at 5 au.The histograms of different colors include particles grouped by a range of sizes, as specified in the caption.Particles of 10 and 100 cm in radius drift inward the fastest (see the left panels).As discussed later, in either case, only a small fraction of these solids is actually accreted.The orbital eccentricity, e s , remains small during the evolution, efficiently damped by gas drag, and so does the orbital inclination (center and right panels, respectively).The inclination of the largest particles and of the smallest icy particles is also damped to small values, i s 0.1 • , by the end of the simulations.Similarly small values for e s and i s , to those in Figure 11, are also found in the other models, with somewhat longer tails in the eccentricity distributions of higher-mass planet models.
As particles radially drift inward, they cross mean-motion resonances with the planet.At those locations, resonant perturbations tend to raise a particle's semi-major axis, whereas gas drag tends to lower it, which may lead to capture in a stable orbit.Weidenschilling & Davis (1985) found that capture is possible for a range of the parameter It can be shown (see Peale 1993) that, in an unperturbed disk, the rate of change of a s of a solid on a near-circular orbit is If K exceeds some critical value, the particle can break through the mean-motion resonance and inward migration continues.Weidenschilling & Davis (1985) quantified this critical value in the limit of small deviations from Keplerian orbits, i.e., for τ s Ω K 2π.(Note that, in Equation ( 21), their original definition is multiplied by a p so as to render the parameter nondimensional.)Kary et al. (1993) performed numerical experiments of particles' capture in exterior resonances, with both planets and smaller planetary embryos, and determined the values of K for which capture is likely.For example, they found that an M p ≈ 0.3 M ⊕ planet orbiting at a p = 5 au in a minimum-mass solar nebula (ρ g ∼ 10 −10 g cm −3 ) can trap R s ∼ 100 m rocky bodies into mean-motion resonances.Since resonant perturbations increase with M p , smaller particles (i.e., larger K) can in principle be trapped by larger planets (see Kary et al. 1993).In fact, for quasi-Keplerian orbits, the critical value of K for resonance capture scales as (Weidenschilling & Davis 1985).
In the two cases shown in Figure 11 (icy and rocky particles), K is large enough for gas drag to overcome the resonant forcing so that the simulated particles can drift inward, inside the planet's orbit.Other models, however, do result in particles segregated in exterior orbits (over the course of the calculation).Some examples are illustrated in Figure 12.The trapped particles in the figure have orbits basically co-planar with the planet, and their eccentricities are typically small, e s H g /a s .The semi-major axis distribution of the 10 m solids (gray histograms) can be very narrow, indicating segregation at locations of mean-motion resonances with the planet.Such are the cases shown in the left panels (4 : 5 and 3 : 4 mean-motion resonances) and in the bottom center panel (6 : 7 Figure 12.Semi-major axis distributions of segregated particles for models with ap = 1 (left), 5 (center), and 10 au (right).In the left panels, ρg represents the unperturbed gas density at r = ap, in units of grams per cubic centimeter.The blue, red, and gray histograms refer to Rs = 1, 100, and 1000 cm particles, respectively.Segregation also occurs in other models not shown.mean-motion resonance), but other instances occur.The equilibrium eccentricities of solids trapped in mean-motion resonances are H g /(2a s ), in accordance with the estimate of Weidenschilling & Davis (1985).
The minimum distance ∆r = r − a p at which particles can be trapped at resonant locations depends on resonance overlap, which may lead to chaotic orbits (Wisdom 1980).Duncan et al. (1989) found that a limit for resonance overlap is given by |a s − a p |/a p ≈ 1.5(M p /M ⋆ ) 2/7 , ranging from ≈ 0.06 to ≈ 0.09 for the cases illustrated in Figure 12.Kary et al. (1993) found an additional limit related to orbital energy dissipation by gas drag and estimated that orbital stability is possible as long as |a s − a p |/a p 1.6[(H g /a p )(M p /M ⋆ )] 2/9 , whose value is comparable to the limit of Duncan et al. (1989) for the conditions found in the RHD calculations.The semi-major axis distributions of 10 m particles trapped in mean-motion resonances satisfy these theoretical limits.
The distributions of smaller particles in Figure 12 tend to be much broader, indicating that their segregation radii are associated to orbital locations where Ω g exceeds its unperturbed value.In the proximity but interior to these radii, drag torques tend to push solids outward.A planetary-mass body exerts a positive tidal torque on the exterior disk, transferring angular momentum and hence augmenting the gas rotation rate relative to that of an unperturbed disk.The torque density per unit disk mass is ∝ (M p /M ⋆ ) 2 (a p /∆r) 4 (Lin & Papaloizou 1986) and peaks for ∆r ≈ H g when R H ≤ H g (D'Angelo & Lubow 2010).When integrated over the disk mass exterior to the planet, one finds the well-known result that the onesided torque exerted by the planet is ∝ (M p /M ⋆ ) 2 (a p /H g ) 3 (see, e.g., Lubow & Ida 2010).Since the torque density function has a peak of finite width, particles well coupled to the gas can attain minimum trapping distances a s − a p ∼ H g (if R H ≤ H g ).Angular momentum transfer beyond |∆r| ≈ 3H g becomes more inefficient (D'Angelo & Lubow 2010), even though tidal perturbations in the gas can extend farther (see Figure 9).Gas density perturbations become increasingly prominent as M p grows, eventually leading to gap formation, which begins when tidal torques exceed viscous torques, i.e., when (M p /M ⋆ ) 2 3πα g (H g /a p ) 5 .In the models presented herein, only the M p = 20 M ⊕ cases at 1 au (marginally) satisfy this condition in the two coldest disks (see Table 1).
The 10 m particles too can stop at equilibrium locations dictated by disk-planet tidal perturbations, as suggested by the broad distributions of a s in, e.g., the top center and bottom right panels of Figure 12.In fact, the values of K are comparable in these two cases, as is the width ∆r of the confinement region.
Segregation would prevent further accretion, isolating the planet from the solids.Lambrechts et al. (2014) argued that isolation can be achieved when M p /M ⊕ ≈ 20 [H g /(0.05a p )] 3 .Dipierro & Laibe (2017) found a comparable limit for τ s Ω K ∼ 1, while Bitsch et al. (2018) estimated a mass larger by ≈ 25% for the conditions realized in the RHD calculations.In the models with a p = 5 and 10 au, these "isolation" masses (≈ 30 M ⊕ and ≈ 50 M ⊕ , respectively) would be much larger than the planet masses modeled herein.Indeed, in none of the models are 10-100 cm particles trapped in exterior orbits.However, at both 5 and 10 au, R s = 1 cm and 10 m particles can be segregated by smaller, 10 M ⊕ planets (see Figure 12).At a p = 1 au, the proposed isolation limit would range from ≈ 8 M ⊕ (coldest, least dense disk) to ≈ 22 M ⊕ (warmest, most dense disk).At least two of the 20 M ⊕ models should thus result in the segregation of solids.At the lowest density (ρ g = 2.0 × 10 −11 g cm −3 ), R s = 1 cm solids are efficiently accreted, while R s ≥ 10 cm solids are segregated.At the intermediate density (ρ g ≈ 10 −10 g cm −3 ), R s ≤ 10 cm particles are efficiently accreted, whereas larger ones are segregated.At the largest density, 1 cm and 10 m solids are segregated.Nonetheless, even 3 M ⊕ planets can segregate solids under appropriate conditions, as shown in the top left panel of Figure 12.  1).Non-accreted solids either are segregated in exterior orbits or drift toward the star.Missing symbols indicate ζ = 0.In the top panel, efficiencies are comparable, at all three values of the gas density, for Rs = 1 and 10 cm particles.
The continued accumulation of solids at segregation sites may enhance collisional comminution, possibly releasing some of the mass in the form of smaller particles, which may drift inward.Coagulation into larger bodies, less coupled to the gas, may also reactivate the accretion of heavy elements on the planet.

Accretion Efficiencies
The efficiency or probability of accretion, ζ, in these models is computed as the ratio of the solids' mass accreted by the planet to the solids' mass interacting with the planet, as they drift inward.Interacting particles can be accreted, transferred to inner orbits, or segregated in the outer disk.Scattering out of the computational domain is possible but rare.
Icy particles can undergo ablation.In the models with a p = 5 au, all particles drifting inward of r ≈ 4 au (T g ≈ 150 K) are consumed to some extent.For a fixed particle temperature T s , the ablation timescale is ∝ R s (DP15).At T s ≈ 150 K, this time is in excess of 10 4 years for 10 m solids, but much shorter for centimeter-sized particles.Since the gas temperature in the proximity of a planet is higher than the average disk temperature at r = a p (DB13), particles can shed mass before being accreted.This process does indeed occur, typically reducing the accreted mass by 15% at a p = 5 au and by 7% at a p = 10 au.The ablated mass is expected to increase in warmer disks and as M p increases.
The accretion efficiencies obtained from the RHD calculations are plotted in Figures 13 and 14.At a p = 1 au, the M p = 3 M ⊕ planet typically accretes rocks with efficiencies 10%, with few exceptions (Figure 13, top panel).The largest, 10 m particles can also be segregated, more often in the M p = 20 M ⊕ planet cases (see the bottom panel).In these latter models, when particles accrete, they do so at high efficiency.
At a p = 5 and 10 au, there are some similarities among the models in terms of accretion efficiency (see Figure 14).Particles of 1 cm in radius are accreted inefficiently (a few percent or less).The missing symbols in the figure indicate zero efficiency, but not necessarily segregation in exterior orbits.Solids of 10 and 100 cm in radius are generally accreted at an efficiency ζ of a few to several percent, with peaks of ≈ 10%.The 10 m particles are accreted efficiently, when accreted, with 0.5 ζ 0.8.However, in a number of cases (M p ≥ 10 M ⊕ ), 10 m icy and rocky solids are segregated in the exterior disk (see also Figure 12).
Overall, the results in Figures 13 and 14 are in general accordance with those from the simplified calculations shown in Figure 8.The values of ζ for 1, 10, and 100 cm solids are comparable.The typically large efficiency for the accretion of 10 m particles (when not segregated) is also in reasonable agreement.
Morbidelli & Nesvorny (2012) estimated small solids' accretion on M p = 1 and 5 M ⊕ planets with a p = 0.8 au.They found efficiencies ranging from a few to several percent for R s ≈ 10-100 cm particles.They also found that 10 m solids have much larger accretion efficiencies, in excess of ≈ 50%.These predictions are generally consistent with those shown in Figures 13 and 14.Although they used an unperturbed surface density twice as large as the largest density applied in the RHD calculations (see Figure 9), the value of τ s Ω K for these particles was comparable, ≈ 0.1-100.Picogna et al. (2018) estimated particle accretion efficiencies from hydrodynamics calculations of laminar and turbulent disks.In models of 5 and 10 M ⊕ planets, they found ζ ≈ 0.2 at τ s Ω K ≈ 0.01 and ≈ 0.4 at τ s Ω K ≈ 100, with a minimum around a few percent at τ s Ω K ∼ 1.Their results are similar to those presented in Figure 14 in the range 0.1 τ s Ω K 10, although their estimate of ζ is somewhat larger for ≈ 1 cm particles and somewhat smaller for ≈ 10 m particles than reported herein.
At an accretion rate of Ṁs ∼ 10 −5 M ⊕ yr −1 (see Figure 5) and ζ ≈ 0.1, a 3 M ⊕ planet at 1 au would double its mass in ∼ 10 6 yr, whereas a 20 M ⊕ planet could grow by up to 1 M ⊕ every 10 5 yr (ζ ≈ 1), if not isolated.At 5 and 10 au, by accreting either icy or rocky R s ≤ 100 cm solids, the accretion rate of a 5 M ⊕ planet would be Ṁp 10 −6 M ⊕ yr −1 , but it could be many times as large if it accreted 10 m particles, resulting in a mass-doubling timescale of several times 10 5 yr.More massive, 10-15 M ⊕ planets would accrete at comparable rates, although the interaction with 10 m solids would typically lead to segregation rather than to accretion.Kary et al. (1993) found that the accretion efficiency is inversely proportional to the radial velocity of a particle, defined by Equation ( 22), in units of Hill radius per synodic period (at ∆r = 2 √ 3R H ), i.e.,

Accretion Efficiency of Bare Cores
hence ζ KLG ∝ 1/K (see Equation ( 21)), for given planet and disk conditions, and proportional to R s /C D for a given material.Their calculations and those presented in Sections 2.2 and 2.3 are most appropriate in the case of small tidal perturbations by the planet, i.e., when If α g ≈ 10 −3 and H g /a p = 0.022(a p /au) 2/7 , this limit would require M p /M ⋆ < 2 × 10 −5 , or ≈ 6 M ⊕ at 5 au from a solarmass star.A set of calculations, similar to those discussed in Section 2.2, were performed to evaluate the accretion efficiency of R s = 1, 10, 100, and 1000 cm rocky/icy solids on small cores, M p ≤ 1 M ⊕ , orbiting at a p = 1, 5, and 10 au from a 1 M ⊙ star.The nebula model is based on that of Chiang & Youdin (2010).The cores are "bare" in the sense that they do not bear a bound atmosphere (of significant mass), an appropriate approximation at such low masses (see, e.g., Kuwahara & Kurokawa 2024).Therefore, these planets are assumed to be composed entirely of condensible materials and solids accrete when they hit the condensed surface of the planet.Although there is no local enhancement of the gas density close to the planet, gas drag operates all the way down to its surface, effectively mimicking a low-density atmosphere.
Results from these calculations are illustrated in Figure 15.The efficiency ζ is plotted as a function of R s , for each core mass (see legend of the top panel).The missing symbols indicate that no collision was detected (ζ = 0).The largest particles tend to accrete more efficiently, although they can become trapped in exterior mean-motion resonances (see the varies by a factor 10).At a p = 1 au, this scaling is restricted to particles with τ s Ω K 0.1.

Formation and Structure Models
The flow rate of 1-1000 cm solids through a disk, Ṁs = Ṁs (t, R s ), obtained from Equation (10) (see, e.g., Figure 5) and the efficiencies illustrated in Figures 13, 14 to formation calculations at 1, 5, and 10 au.Accretion of rocks (a p = 1, 5 au) and ice (a p = 5, 10 au) is considered in separate models.The total mass of the disk is ≈ 0.055 M ⊙ at time t = 0.The initial gas-to-solids mass ratio is 100 and 71.5 (Pollack et al. 1994) in the rocky and icy disk, respectively.Formation starts from Moon-mass embryos, M p = 10 −2 M ⊕ at t = 1000 orbits (at each orbital distance).Assuming mass equipartition among particle sizes, these embryos take ≈ 0.25 to ≈ 0.55 Myr to attain ≈ 1 M ⊕ , when structure calculations begin ( Ṁp = ṀZ at smaller masses).
Planet evolution is modeled as in D' Angelo et al. (2014), by solving the structure equations in spherical symmetry (Kippenhahn et al. 2013), including mass and energy deposited by incoming solids and gas, particle break-up and ablation.Dust opacity in the envelope is computed self-consistently with the sedimentation and coagulation of grains released by accreted gas and solids (Movshovitz et al. 2010).Dissolution and mixing of heavy elements in H-He gas (e.g., Bodenheimer et al. 2018) is neglected and all heavy elements settle to the condensed core.When the planet is in contact with the nebula, its radius R p is evaluated as in Lissauer et al. (2009), and depends on both the Bondi and Hill radius.Mass loss is enabled when the envelope overflows its boundary.Disk-limited accretion of gas is accounted for as in D' Angelo & Bodenheimer (2016).Disk evolution is modeled as outlined in Section 2.1 (see Figure 4).Formation continues until gas dispersal, at t ≈ 3.5 Myr.
The heavy-element and H-He masses of the planets are shown in Figure 16.Models at 1 and 5 au accreting silicates achieve heavy-element masses M Z ≈ 3 and 3.4 M ⊕ , and H-He inventories contributing respectively 12% and 24% of the total masses.Instead, models accreting ice become gas giants, M p = 321 M ⊕ (M Z = 6.9 M ⊕ ) at 5 au and M p = 46 M ⊕ (M Z = 4.8 M ⊕ ) at 10 au.For an initial gas-to-solids mass ratio equal to 100, the heavy-element inventory at this latter distance would be limited to ≈ 1 or 2.5 M ⊕ by accreting rocks or ice, respectively.
The different outcomes at 5 au can be attributed to the different amounts of solids transported across the planet's orbit during formation: ≈ 150 M ⊕ in the case of ice and ≈ 90 M ⊕ in the case of rocks.Had the initial gas-to-solids mass ratio in the former case been 100 (as in the rocky disk), this mass would have been ≈ 107 M ⊕ .The cumulative accretion efficiency, i.e., M Z divided by the amount of solids crossing the orbit of the planet, also contributes, with that of icy particles (≈ 0.046) exceeding by about 20% that of rocky particles.In both models, most of the heavy-element mass (≈ 70%) is delivered by 10 m boulders/ice blocks.Only 10%-15% of M Z is supplied by 10-100 cm solids.
In the models at 1 and 10 au, the total solid mass crossing the planet's orbit is ≈ 103 M ⊕ (of rocks) and 136 M ⊕ (of ice), respectively, corresponding to cumulative accretion efficiencies of 0.03 and 0.04.In these two calculations, about half of M Z is delivered by 100 cm boulders, at 1 au, and by 10 m ice blocks, at 10 au.Note that a 50/50 solid mixture (by mass) of ice and silicates would have a material density closer to that of ice and, therefore, is expected to behave somewhat closer to ice than to silicates in terms of transport and accretion.
The H-He mass at the beginning of the structure calculations, ≈ 10 −5 M ⊕ , is negligible and remains so for ∼ 1 Myr (see Figure 16).In fact, the total accretion rate Ṁp remains close to ṀZ for a large part of the formation epoch (see left panels of Figure 17).The maximum of ṀZ is achieved between M p = 0.5 and 1 M ⊕ when heavy elements are supplied by silicates, and between M p = 1 and 2 M ⊕ when they are supplied by ice.At 1 au, the accretion rate of gas exceeds 50% that of solids relatively late in the evolution, after 1.6 Myr, when ṀZ 4 × 10 −7 M ⊕ yr −1 .The planet accreting rocks at 5 au first reaches this condition around 1.5 Myr ( ṀZ ≈ 10 −6 M ⊕ yr −1 ; see the middle panel), and then again when ṀZ 4×10 −7 M ⊕ yr −1 .Mass loss by overflow is significant in this case during the first ≈ 2 Myr of evolution.In the models resulting in the formation of gas giants, the accretion rate of gas exceeds ṀZ after ≈ 1 Myr, when the heavy-element mass M Z has nearly attained its final value (see Figure 16).
Instead of mass equipartition, if all of the initial mass had been carried by particles of a single size, in general M Z would have been largest by accreting R s = 1-10 m boulders/ice blocks.However, the final heavy-element inventory of a planet is not only determined by the mass distribution of particles, but also by the accretion of gas, since the accretion efficiency ζ is a function of M p .
The luminosity during formation is shown in the right panels of Figure 17.Prior to beginning the structure calculations (M p 1 M ⊕ ), the luminosity is computed as L = GM Z ṀZ /R p , where R p is the radius of the condensed core.The curves are color-coded by the total accretion rate.The first luminosity peak corresponds to the maximum of ṀZ .For the gas giants, a second peak occurs around the maximum of the gas accretion rate.Although the planet accreting silicates at 5 au exhibits a maximum in gas accretion around 2.6 Myr, the corresponding luminosity peak is smaller than the first, because of the larger planet radius.In fact, both lower-mass planets (i.e., those accreting silicates) remain in contact with the disk until it completely disperses.Therefore, their envelopes are extended throughout formation.The gas giants, instead, detach from the disk when t ≈ 1.2 and 2.2 Myr at 5 and 10 au, respectively.They undergo disk-limited accretion thereafter.
The dissolution of heavy elements in the interior of the planets may produce a stratification in composition and possibly affect energy transport and temperatures (e.g., Bodenheimer et al. 2018;Stevenson et al. 2022).Although heavy-element deposition resulting from the accretion of small solids can differ from that resulting from planetesimal accretion, the radial composition of an interior may not provide distinctive information on the original carriers after sustained mixing.Moreover, the deposition of refractory material may not clearly differentiate between the accretion of small and large solids, as deposition patterns in the planet's envelope can overlap (see Appendix B).

Discussion and Conclusions
The accretion of small solids (10 −2 τ s Ω K 10 2 ) on a planet depends on the transport of solids through the disk, dM s /dt, and on the efficiency or probability of accretion, ζ (see Equation ( 25)).The former quantity is non-local, and determined by the global evolution of the disk's gas, the initial distribution of solids, coagulation, and fragmentation.The latter quantity (i.e., the fraction of dM s /dt at a p that is accreted by the planet) is determined by multiple processes, at orbital or smaller scales.Small cores without an envelope of significant mass) allow for relatively straightforward calculations of ζ, including planet-induced velocity perturbations on the gas (see Section 4).However, as the planet grows larger, the computation of ζ becomes increasingly complex.In fact, tidal interactions between the planet and the disk have important effects on the dynamics of the solids.Close-range, 3D gas thermodynamics is also expected to significantly affect ζ.Tidal segregation (and capture in mean-motion resonances) can impede the flux of solids through the planet's orbit, a process that depends on solid size and composition, planet mass, and on the gas thermal state (and, therefore, on disk age and mass).The epoch of accretion is also limited by the duration of the solids' transport across the planet-forming region.Boulder-sized solids tend to accrete more efficiently than smaller solids, unless they become segregated.A disk of initial mass 0.055 M ⊙ in gas and ≈ 260M ⊕ in solids would provide typical accretion rates of heavy elements (in the 1-10 au region) between 10 −6 and 10 −5 M ⊕ yr −1 , for the first 1-1.5 Myr of evolution, and 10 −6 M ⊕ yr −1 thereafter (see Figure 17).The results presented in Figures 4 and 5 are specific to the applied (global) disk conditions, but appropriate solutions can be found by integrating Equations 9 and 10 in relevant gaseous disks.The efficiencies of accretion presented in Figures 13 and 14 are more general, since they only depend on local quantities, and can be used as long as the gas densities and temperatures are comparable to those found in the RHD calculations (see Table 1).The same considerations apply to the efficiencies of small cores with no significant envelope reported in Figure 15.
The formation calculations presented herein apply accretion rates of solids based on simple, but consistent, models of gas and solid disk evolution and first-principles determinations of accretion efficiencies.The structure calculations include detailed computations of grain opacity as the solids travel through the envelope and release material.The outcomes encompass a broad range, from low-mass, heavy-element-dominated planets to H-He gas giants.Heavy-element masses from 3 to 7 M ⊕ require an amount of solids from ≈ 90 to 150 M ⊕ to cross the planet's orbit.The maximum luminosity emitted during the epoch of solids' accretion ranges from ≈ 10 −7 to ≈ 10 −6 L ⊙ .Cumulative efficiencies of accretion in the formation calculations are below 5%.Multiple planets may form if solids can drift inward, past a planet.Heavy elements would be delivered during the epoch when dM s /dt at a given distance is non-negligible.As outer planets remove solids, lesser amounts would be available to inner planets.The higher the efficiency of accretion in a size range, the lower is the availability of solids to planets downstream.In addition, because of segregation, a planet large enough (a few to several times M ⊕ ; see Figure 12) could hinder or terminate the supply of solids (in some size ranges) to inner planets.In summary, predictions of heavy-element delivery (to one or more planets) require accounting for both local and remote processes, over the disk lifetime.Considerations of the evolving size distribution of solids are also necessary, including at segregation sites, where enhanced comminution and/or growth may release solids and promote further accretion.
Figure 18.Azimuthal averages (around the star) versus radial distance of the midplane gas density (left) and midplane gas temperature (right) at the quasi-steady state Si (green) and at the end of the test calculation (red), several hundred planet's orbits later.Curves at the two epochs overlap one another in both panels (for better visibility, the green curves are thicker than the red curves).
Figure 19.Histograms of the particles' semi-major axes (left), eccentricities (center), and inclinations (right) as they approach the planet's orbit.The top panels refer to the model in which the evolution of both the gas and the solids is calculated.In the bottom panels, the gas quasi-steady state Si (defined in the text) is applied and only the evolution of the solids is calculated.The gray histograms include all particles; the red, green, and blue histograms include particles with radii Rs ≤ 100 cm, Rs ≤ 10 cm, and Rs = 1 cm, respectively.planet and do not change significantly over the course of the calculations.The approach is intended to mitigate the computational requirements of those RHD models and render the problem more tractable.Here a test is presented on the validation of this approach.
An RHD calculation of a 3 M ⊕ planet, embedded in a disk at r = a p = 1 au from a 1 M ⊙ star, is performed and evolved until the system attains a state of quasi-equilibrium, hereafter referred to as state S i .The setup of the model is equivalent to that of the models with a p = 1 au discussed in Section 3 (reference density ρ g = 1.5 × 10 −10 g cm −3 ), but at a somewhat lower grid resolution.The planet radius is assumed to be approximately equal to R B .Four size bins, R s = 1, 10, 100, and 1000 cm, are populated with equal numbers of rocky particles.They are deployed outside the planet's orbit by applying the same random process as described in Section 3.3.
The calculation is then continued, from state S i , by evolving both the gas and the solids for several hundred orbits of the planet.During this time, the disk's gas density and temperature (red curves) do not deviate significantly from those of state S i (green curves), as illustrated in Figure 18 (in which curves basically overlap and are difficult to distinguish).In another calculation, the same distribution of solids evolves in the gas fields provided by state S i .
The stopping times of the solids at deployment range from τ s Ω K ≈ 10 −2 (R s = 1 cm) to τ s Ω K ≈ 10 3 (R s = 1000 cm).Intermediate size particles (0.1 τ s Ω K 10) drift toward the planet relatively quickly.For solids in the various size ranges, Figure 19 shows the histograms of the particles' semimajor axes (left), orbital eccentricities (center), and inclinations (right) during the approaching phase (see the figure's caption for further details).The orbits remain nearly circu- lar, although the smallest particles are somewhat more eccentric than the largest ones (see also Figure 11).The inclination with respect to the midplane of the disk is basically zero, except for that of the largest particles, which is damped over a longer timescale (∝ R s , Adachi et al. 1976).In statistical terms, the orbital elements of the solids in the two calculations remain quantitatively similar.By the end of the simulations, the R s = 10 and 100 cm particles have either accreted on the planet or moved toward interior orbits, whereas the smallest and largest particles continue drifting toward the planet.
In order to gauge variations in the accretion rates of solids on the planet, Figure 20 shows the relative difference of the accreted mass versus time, since state S i , between the two calculations.The R s = 10 cm particles intercept the planet's orbit first, followed by the R s = 100 cm particles.The differences in accreted mass are below ≈ 10%, and typically at the level of a few percent on average.The accretion efficiencies ζ are also in close agreement, within ≈ 1%.The fact that the accretion rates of the particles are statistically consistent during close encounters with the planet implies that the particle dynamics is indeed similar in the two simulated scenarios, even during the phases of the evolution preceding accretion.

B. Dissolution of Solids in Planetary Envelopes
Solid bodies moving through an atmosphere warm up and ablate, disseminating heavy elements along their path.This process can affect the composition of planetary atmospheres and interiors.We performed experiments of solids' dissolution in planetary envelopes taken from the giant planet formed at 5 au and presented in Section 5.The envelopes refer to times at which their masses are 0.15 M ⊕ (M p = 6.8 M ⊕ ) and 6.9 M ⊕ (M p = 13.8M ⊕ ).To circumvent segregation issues, solids were deployed around the planet's L 2 point (or close enough to allow for accretion).The spread of the results due to possible differences in entry velocities was not assessed.Outside of the envelope, the disk is modeled as in the planet structure calculation (see Section 5).
The energy budget of a particle includes heating by frictional work and by the ambient thermal field of the envelope, and cooling by black-body radiation at the surface temperature and by latent heat release through phase change.Above the critical temperature of the material, all energy input drives mass loss (see DP15, for details).A body can break-up when subjected to a dynamical pressure larger than the compressive strength of the material, although break-up did not occur in these experiments.
The heavy-element deposition profiles of the rocky (SiO 2 ; top) and icy solids (H 2 O; bottom) are illustrated in Figure 21.The high-mass envelope planet is shown on the right (the shaded areas represent the condensed core region).Upon entry in the envelope, and prior to significant phase change, the surface temperature of R s 10 cm particles tends to attain values similar to gas temperatures.In these cases, significant ablation begins only after envelope temperatures rise above some value, dictated by the vapor pressure curve of the material.When frictional heating exceeds warming by the local thermal field, deposition may begin at lower envelope temperatures.In the top panels of Figure 21, because of the contribution of frictional heating, 1 km planetesimals can disseminate heavy elements over a broad range of depths, including layers in which small solids would ablate.At all sizes considered in these experiments, about 50% of the input silicate mass would be released in the outer 15%-20% (top left) and outer 25%-35% (top right) of the envelope mass.Small icy solids (bottom panels) would be disseminated farther up in the planet, compared to 1 km planetesimals.Nonetheless, for all sizes shown in the figure, 50% of the accreted ice mass would be dispersed in the outer 5-15% of the envelope mass.Thus, in all cases shown in Figure 21, most of the accreted solids would not settle to the core, as assumed in Section 5, but would result in an envelope enriched in heavy elements.

Figure 1 .
Figure 1.Stopping time τs (top) and orbital decay time τD (middle) of rocky (solid lines) and icy (dashed lines) particles of radius Rs, orbiting in a smooth and unperturbed protoplanetary disk.The time τs is computed by inverting Equation (2) and τD = as/| ȧs| (see Equation (5)).The drag coefficient CD (DP15) is also shown (bottom).The different lines refer to different disk locations and/or gas densities (as indicated in g cm −3 ) and temperatures (see also Table1).

Figure 3 .
Figure 3. Top: Solution of Equation (10) by imposing τsΩK = 1 in Equation (9) and the initial gaseous nebula model of Chiang & Youdin (2010) (see the text for details).The times in the legend are in years.The disk's gas scale-height Hg is constant in time.Bottom: Accretion rate of solids corresponding to the evolution of Σs in the top panel.

Figure 4 .
Figure 4. Solutions of Equation (10) for rocky solids of radius Rs, equal to 1 (left), 10 (center), and 100 cm (right).From top to bottom, the panels show Σs, dMs/dt, and τsΩK.The times in the legends are in years.The black curves in the top panels represent the initial distribution of solids.See the text for further details.

Figure 5 .
Figure 5. Solutions of Equation (10) for solids of radius Rs, as indicated, for the accretion rate through the disk, dMs/dt.The panels illustrate the evolution at distance r = 1 (top), 5 (middle), and 10 au (bottom).The solid lines represent results from the disk model in Figure4.The dashed lines indicate results from similar models with a source term, Λ (see the text for details).

Figure 6 .
Figure6.Trajectories in the rotating frame of the planet of Rs = 1 (green), 10 (orange), and 100 cm (purple) particles, orbiting in a disk perturbed by a 5 (top) and a 15 M⊕ (bottom) planet on circular orbits.The trajectory of the smallest particle is initiated as the others, but a smaller portion is plotted.The gray lines represent streamlines of the perturbed gas.The unperturbed disk properties are as in Figure1at 5 au (see also Table1).

Figure 8 .
Figure8.Histograms of the closest approach, in units of the planet's Hill radius, along the trajectories of Rs = 1, 10, 100, and 1000 cm particles (as indicated).Each bin is normalized by the total mass used in the experiment for the given particle size.Bins are stacked in order of increasing Rs, rendered by different colors.The dashed line indicates the planet's Bondi radius.The unperturbed disk properties correspond to those of Figure1, according to ap.The case at 1 au (top panel) uses the highest gas density in Figure1.

Figure 10 .
Figure 10.Close-up of log ρg (color scale) and log Tg (contour lines), at the disk's midplane, for the 3 (top) and 20 M⊕ (bottom) planet located at ap = 1 au.The gas density and temperature are in units of grams per cubic centimeter and kelvin degrees, respectively.The unperturbed surface density at r = ap is Σg ≈ 220 g cm −2 .

Figure 13 .
Figure13.Accretion efficiencies versus particle radius, obtained from the RHD calculations at ap = 1 au, for an Mp = 3 (top) and 20 M⊕ (bottom) planet.The different colors/symbols refer to different values of the reference gas density, as indicated (see also Table1).Non-accreted solids either are segregated in exterior orbits or drift toward the star.Missing symbols indicate ζ = 0.In the top panel, efficiencies are comparable, at all three values of the gas density, for Rs = 1 and 10 cm particles.

Figure 14 .
Figure 14.Accretion efficiencies obtained from the RHD calculations at ap = 5 (top) and 10 au (bottom).Icy and rocky solids are represented on the left and right, respectively.The different colors/symbols refer to different planet masses, as indicated.The different line styles connect symbols of the same type and color.Missing symbols indicate ζ = 0. Non-accreted solids either are segregated on exterior orbits or drift toward the inner disk.Icy particles can also evaporate.
in the top panel).The mass scaling in Figure 15 agrees reasonably well with the predictions of Equations (14) and (23), ζ ∝ M 2/3 p , at a p = 5 and 10 au (ζ/M p 2/3 , and 15 (appropriately interpolated in planet mass), ζ = ζ(M p , R s ), are combined to provide a size-integrated accretion rate,

Figure 15 .
Figure 15.Accretion efficiencies of bare cores (i.e., without an atmosphere), versus particle radius, for rocky (red colors) and icy particles (green colors).The planet's semi-major axis is ap = 1 (top), 5 (middle), and 10 au (bottom), and the mass is indicated in the legend.The disk model is the same as that in Figure 4. Missing symbols indicate ζ = 0.

Figure 16 .
Figure 16.Heavy-element and H-He masses obtained from formation models at ap = 1 (top), 5 (middle), and 10 au (bottom).The heavyelement mass is color-coded by the solids' accretion rate, in units of M⊕ yr −1 .The thinner lines in the middle panel refer to the planet accreting icy particles.

Figure 17 .
Figure17.Accretion rates versus time (left) and luminosity versus mass (right) from formation models at ap = 1 (top), 5 (middle), and 10 au (bottom).The accretion rate of solids is color-coded by the planet mass (in units of M⊕) and the luminosity by the total accretion rate (in units of M⊕ yr −1 ).The thinner lines in the middle panels refer to the planet accreting icy particles.

Figure 20 .
Figure 20.Normalized relative difference in the accretion rate of solids on the planet versus time, since the reference time of state Si (defined in the text), between the two calculations described in the text.Accretion first involves Rs = 10 cm particles and then Rs = 100 cm particles.

Table 1 .
Properties of Planets and Disks in the RHD Models ap [au] Mp/M⊕ RB/ap Rp/RB a Rp/RH a Tg b [K] ρg b [g cm −3 ]