Planetesimal Accretion at Short Orbital Periods

We begin our analysis by considering a disk of equal-mass planetesimals with radius r_pl, mass m_pl, and surface density Σ_pl. The collision rate in the vicinity of an orbit defined by Keplerian frequency Ω can be written as nΓv, where n = Σ_plΩ/2m_pl v (where we have assumed that the scale height of the planetesimal disk goes as 2v/Ω). Γ describes the effective collision cross section, and v is the typical encounter velocity between planetesimals. For a swarm of planetesimals on randomly oriented orbits, v is typically taken to be the rms velocity, ${\left\langle {v}^{2}\right\rangle }^{1/2}$, which can be related to the eccentricity and inclination distribution (e, i) in the following way (Lissauer & Stewart 1993):

$$\begin{eqnarray}&&{\left\langle {v}^{2}\right\rangle }^{1/2}=\left(\displaystyle \frac{5}{4}{\left\langle {e}^{2}\right\rangle }^{1/2}+{\left\langle {i}^{2}\right\rangle }^{1/2}\right){v}_{k}.\end{eqnarray} \tag{ 1 }$$

The dynamical interactions between growing planetesimals can be somewhat simplified by scaling the orbital elements of the bodies by the Hill radius,

$$\begin{eqnarray}&&{r}_{h}={\left(\displaystyle \frac{{m}_{\mathrm{pl}}}{3{M}_{* }}\right)}^{1/3},\end{eqnarray} \tag{ 2 }$$

where M_* is the mass of the central star. The Hill radius of a body describes the size scale over which the gravity of the growing planetesimal dominates over the gravity of the star. Using Equation (2), the eccentricity, inclination, and separation between orbiting bodies can be defined as

$$\begin{eqnarray}&&{e}_{h}=\displaystyle \frac{{ae}}{{r}_{h}},{i}_{h}=\displaystyle \frac{{ai}}{{r}_{h}},\tilde{b}=\displaystyle \frac{{a}_{2}-{a}_{1}}{{r}_{h}}.\end{eqnarray} \tag{ 3 }$$

Using this formalism, e_h and i_h describe the radial and vertical excursions of an orbiting body in units of its Hill radius. For e_h > 1, the random velocity dispersion dominates over the shearing motion across a separation of 1r_h and encounters can be treated with a two-body formalism.

Assuming that every collision results in a perfect merger, the growth rate of a planetesimal is given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{M}\displaystyle \frac{{dM}}{{dt}}=\displaystyle \frac{{\rm{\Sigma }}{\rm{\Omega }}}{2{m}_{\mathrm{pl}}}{\rm{\Gamma }}.\end{eqnarray} \tag{ 4 }$$

In the case where the collision cross section depends only on the physical size of the planetesimals, the growth scales sublinearly with mass and the mass distribution is expected to evolve in an "orderly" fashion, in which mass ratios between bodies tend toward unity. However, bodies larger than ∼100 km in size are expected to exert a significant gravitational force on each other during encounters, and the collision cross section depends on both the size of the bodies and their encounter velocities. In this case,

$$\begin{eqnarray}&&{\rm{\Gamma }}=\pi {r}_{\mathrm{pl}}^{2}\left(1+{v}_{\mathrm{esc}}^{2}/{v}^{2}\right)\end{eqnarray} \tag{ 5 }$$

(Safronov 1969), where v_esc is the escape velocity from the two bodies at the point of contact.

In the limit that v_esc ≫ v, it can be shown that dM/dt ∝ M^4/3, which implies a runaway scenario where growth accelerates with mass. This mode of growth was confirmed with N-body simulations by Kokubo & Ida (1996) and appears necessary to construct protoplanets within the lifetime of a protoplanetary disk (Lissauer 1987), although one should note that pebble accretion (Lambrechts & Johansen 2012, 2014; Bitsch et al. 2015) is a viable alternative scenario. Due to the velocity dependence of the gravitational focusing effect, it is not clear how ubiquitous this mode of growth is. In particular, encounter velocities at short orbital periods will be rather large (because v ∼ v_k ) and the v_esc ≫ v condition may not always be satisfied. The effect that a dynamically hot disk has on runway growth will be examined in detail in Section 4.

An important feature is missing from the model described above, which limits its applicability at late times. Gravitational stirring, which converts Keplerian shear into random motion, raises the typical encounter velocity between planetesimals over time (Weidenschilling 1989; Ida 1990) and diminishes the effectiveness of gravitational focusing. As the mass spectrum of the system evolves away from uniformity, these velocity differences become even more pronounced. As the system evolves, it tends toward a state of energy equipartition where v ∼ m^1/2. For a system of equal-mass bodies in which encounters are driven by random motions rather than Keplerian shear (dispersion dominated), the timescale for gravitational stirring is described by the two-body relaxation time (Ida & Makino 1993),

$$\begin{eqnarray}&&{t}_{\mathrm{relax}}=\displaystyle \frac{{v}^{3}}{4\pi {{nG}}^{2}{{m}_{\mathrm{pl}}}^{2}\mathrm{ln}{\rm{\Lambda }}},\end{eqnarray} \tag{ 6 }$$

where $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm, typically taken to be ≈3 for a planetesimal disk (Ida 1990; Stewart & Ida 2000). Despite the fact that the behavior of gravitational stirring is described well by a two-body formula, Ida & Makino (1993) found that the stirring in a planetesimal disk is actually driven by close encounters, which requires a three-body formalism. As we will show in Section 4, gravitational stirring effectively shuts off when the Hill sphere of a body becomes comparable to its physical size. In this case, close encounters tend to result in collisions, and the main pathway for energy exchange between planetesimals and growing protoplanets is unable to operate.

Kokubo & Ida (1998) showed that the runaway growth process described above is actually self-limiting. As the runaway bodies develop, they become increasingly effective at dynamically heating the remaining planetesimals, which diminishes the gravitational focusing cross sections and throttles the growth rate. Around the time that the mass of the runaway bodies exceeds the mass of the planetesimals by a factor of ∼50–100 (Ida & Makino 1993), a phase of less vigorous "oligarchic" growth commences, in which the largest bodies continue to accrete planetesimals at similar rates, independent of mass.

Regardless of the mechanism that eventually limits the growth of the planetary embryos, a maximum estimate for the masses produced during this stage of accretion can be obtained using the initial solid surface density profile. A growing protoplanet is expected to accrete material within a distance $b=\tilde{b}{r}_{h}$ of its orbit. The total mass of planetesimals within this distance is then $2\pi a\left(2\tilde{b}{r}_{h}\right){{\rm{\Sigma }}}_{\mathrm{pl}}$, where a is the semimajor axis of the growing protoplanet. The "isolation mass" of the protoplanet can then be obtained by setting the protoplanet mass equal to the total mass of planetesimals within accretionary reach such that

$$\begin{eqnarray}&&{M}_{\mathrm{iso}}=4\pi {a}^{2}\tilde{b}{\left(\displaystyle \frac{{M}_{\mathrm{iso}}}{3{M}_{* }}\right)}^{1/3}{{\rm{\Sigma }}}_{\mathrm{pl}}.\end{eqnarray} \tag{ 7 }$$

Solving for M_iso gives

$$\begin{eqnarray}&&{M}_{\mathrm{iso}}={\left[\displaystyle \frac{{\left(2\pi {a}^{2}{{\rm{\Sigma }}}_{\mathrm{pl}}\tilde{b}\right)}^{3}}{3{M}_{* }}\right]}^{1/2}.\end{eqnarray} \tag{ 8 }$$

For bodies on circular, noninclined orbits, $\tilde{b}=2\sqrt{3}$ is the smallest orbital separation that produces a non-negative Jacobi energy and permits a close encounter (Nakazawa & Ida 1988). This value of $\tilde{b}$ is typically used to calculate the final isolation mass of a protoplanet because oligarchic growth tends to maintain near-circular orbits for the growing protoplanets.

The picture described above relies upon a crucial assumption, which is that the mass distribution evolves slowly enough for gravitational stirring to maintain energy equipartition. In other words, the relaxation timescale must remain short relative to the collision timescale. For typical conditions near the terrestrial region of the Solar System, this timescale condition is satisfied. Due to the steep dependence of the relaxation time on encounter velocity, however, this condition can easily be violated at shorter orbital periods.

In Figure 1, we show the ratio between the relaxation and collision timescale for a population of equal-mass planetesimals as a function of orbital period. Here, the encounter velocity is described by Equation (1). For simplicity, we assume that ${\left\langle {e}^{2}\right\rangle }^{1/2}=2{\left\langle {i}^{2}\right\rangle }^{1/2}$ (Ida et al. 1993) and that the eccentricity dispersion is constant with orbital period. The coloring indicates the ratio between t_relax and t_coll. The dashed lines denote the eccentricity at which the random encounter velocity (calculated according to Equation (1)) is equal to the mutual escape velocity of the bodies. This is shown for planetesimals with an internal density of 3 g cm⁻³ ranging from 10 to 200 km in size.

Zoom In Zoom Out Reset image size

A physically realistic value of e_h for a population of planetesimals is going to depend on the structure of the gaseous disk (which we will address in Section 5). However, the eccentricity dispersion will typically increase over time until ${\left\langle {v}^{2}\right\rangle }^{1/2}={v}_{\mathrm{esc}}$, and this is often used to construct the initial conditions (e.g., Barnes et al. 2009). Therefore, the curves in this figure should be interpreted as upper limits.

The timescale criterion for oligarchic growth is only satisfied in regions where the disk is sufficiently dynamically cold and the orbital period is sufficiently long. In Sections 4 and 5, we will explore the behavior and outcome of planetesimal accretion in regions where this criterion is not satisfied.

In the formalism described above, the mass and velocity distribution of the bodies are both a function of time. Due to the interdependence of these quantities, it is not clear whether the ratio between the relaxation and collision timescales will remain constant as the oligarchs develop. In the case of the Solar System, t_relax ≪ t_coll likely continued to remain true, otherwise runaway growth would have consumed all of the small bodies and there would be nothing left to populate the asteroid or Kuiper Belt. In the case where t_coll ≪ t_relax, however, it is not clear how the system might evolve.

An insight into the expected behavior in this regime can be gained by defining a dimensionless parameter, α, which is the ratio between the physical size of a body and its Hill radius, r_h :

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{r}_{\mathrm{pl}}}{{r}_{h}}=\displaystyle \frac{1}{a}{\left(\displaystyle \frac{9{M}_{\star }}{4\pi {\rho }_{\mathrm{pl}}}\right)}^{1/3},\end{eqnarray} \tag{ 9 }$$

where a is the semimajor axis of the body and ρ_pl is its bulk density. Assuming that the bulk density stays constant as bodies collide and grow, and that no large-scale migration occurs, the scaling of both r_pl and r_h as ${m}_{\mathrm{pl}}^{1/3}$ means that α will be constant with time. For a composition of ice and rock, α is small for any presently populated region of the Solar System (α ∼ 10⁻² near Earth and α ∼ 10⁻⁴ in the Kuiper Belt). As one moves closer to the Sun, α becomes larger than 1, which implies that the Hill sphere of a body becomes smaller than its physical size. As an additional note, the Roche limit of the central body and the distance at which α = 1 are equivalent for a rigid spherical body. After applying a hydrostatic correction to the Roche limit, a_Roche is equivalent to 0.6 times this distance. This accretion mode should therefore be relevant for planetary ring systems, which is a topic that deserves further study using high-resolution N-body techniques.

The magnitude of α controls the relative importance of gravitational scattering and collisions in driving the evolution of the planetesimal disk. When α is small, most close encounters will result in a gravitational interaction, moving the system toward a state of relaxation. If, however, the Hill sphere is largely filled by the body itself, these same encounters will instead drive evolution of the mass spectrum. Because α stays constant with mass, the boundary in the disk where collisions or gravitational encounters dominate will stay static with time.

We also introduce a second dimensionless quantity, which relates the physical size of the bodies to the velocity state of the system:

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{r}_{\mathrm{pl}}}{{r}_{g}},\end{eqnarray} \tag{ 10 }$$

where r_g = Gm_pl/v² is the gravitational radius of a body (see Equation (4.1) of Ida 1990). Encounters between bodies inside of a distance of r_g result in significant deflections of their trajectories. It should be noted that the gravitational focusing enhancement factor ${v}^{2}/{v}_{\mathrm{esc}}^{2}$ is equal to 1 for β = 1. In the case where r_g is smaller than the size of a planetesimal, the gravitational focusing enhancement factor will be between 0 and 1, and the collision cross section is mostly set by the geometric value. For very large values of r_g (β ≪ 1), the effective collision cross section is almost entirely set by gravitational scattering.

These scaling considerations motivate the range of parameters we choose for the numerical experiments presented in the next section, where we aim to understand where and when runaway and oligarchic growth can operate. In Figure 2, we show a schematic that relates α and β to the geometry of a two-body encounter. For large values of α, the Hill radius of a body (dashed circle) becomes comparable to its physical radius (solid circle). As β increases, the trajectory of a body undergoing a flyby is less affected by the encounter.

Zoom In Zoom Out Reset image size

Motivated by the qualitative dependence of accretion on α, we next investigate whether this highly efficient, non-oligarchic growth should be expected to operate near the innermost regions of a typical planet-forming disk. Given that N-body simulations of short-period terrestrial planet formation typically begin with a chain of neatly spaced, isolation-mass (see Kokubo & Ida 2000 Equation (20)) protoplanets, it is pertinent to determine whether the high-α growth mode we revealed in the previous section invalidates this choice of initial conditions.

Given the dearth of short-period terrestrial planets observed around M stars (e.g., TRAPPIST-1 (Gillon et al. 2016, 2017; Agol et al. 2021)), we chose to model the evolution of a series of wide planetesimal disks, which span from 1 to 100 days in orbital period, orbiting a late-type M star of mass 0.08 M_⊙. For a population of planetesimals with a bulk density of 3 g cm⁻³, this orbital period range corresponds to α ∈ (0.7, 0.05). By simultaneously modeling a broad range of orbital periods, we can determine the critical value of α that divides these two modes of accretion, and also explore how the oligarchic/non-oligarchic accretion boundary affects the resulting distribution of protoplanets.

Four wide-disk simulations are run in total (see Table 1). In each case, the solid surface density follows a power-law profile

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=10\,{\rm{g}}\,{\mathrm{cm}}^{-2}\times A\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right){\left(\displaystyle \frac{r}{1\,\mathrm{au}}\right)}^{-\delta },\end{eqnarray} \tag{ 11 }$$

where M_* is the mass of the central star, 10 g cm⁻² is the surface density of the minimum-mass solar nebula (Hayashi 1981, MMSN) at 1 au, A is an enhancement factor, and δ is the power-law index. In the first case (fdHi), we model a disk that follows a MMSN power-law slope (δ = 1.5), with the overall normalization enhanced by a factor of 100. This choice of normalization for the solid surface density profile appears necessary in order to reproduce many observed short-period terrestrial worlds in situ (Hansen & Murray 2012). Many recent planetesimal formation models invoke streaming instabilities, which first require sufficiently large amounts of solid material concentrating at preferential locations in the disk. This can occur via zonal flows (Johansen et al. 2009; Simon et al. 2012), vortices (Klahr & Bodenheimer 2003), or through mechanisms that produce a pressure bump in the gas disk, such as an ionization (Lyra et al. 2008) or condensation front (Brauer et al. 2008; Drkażkowska et al. 2013), or even the perturbations from an existing planet (Shibaike & Alibert 2020). Drkażkowska & Dullemond (2018) showed that evolution of the snow-line boundary can cause planetesimal formation over a significant area of the disk, producing mass concentrations at least as large as the ones we use here. We argue that this is a particularly attractive mechanism for widespread planetesimal formation around low-mass stars, as their extreme pre-main-sequence evolution is particularly conducive to significant movement of the condensation fronts (Baraffe et al. 2015).

Additionally, we vary the power-law index (fdHiShallow, fdHiSteep) and overall normalization (fdLo) of Σ(r). Although there is no way to reliably measure the uncertainty on the MMSN power-law slope, Chiang & Laughlin (2013) applied a similar analysis to a sample of Kepler multiplanet systems and found a variation of ∼0.2. Submillimeter observations of the outer regions of cold protoplanetary disks find that δ can be as low as 0.5 (Mundy et al. 2000; Andrews et al. 2009, 2010). Therefore, we vary δ by a value of 1.0 relative to the MMSN value for the fdHiShallow and fdHiSteep simulations.

In all cases, the eccentricities and inclinations of the bodies are randomly drawn from a Rayleigh distribution, with ${\left\langle {e}^{2}\right\rangle }^{1/2}=2{\left\langle {i}^{2}\right\rangle }^{1/2}={e}_{\mathrm{eq}}$ (Ida & Makino 1993). Following Kokubo & Ida (1998), the value of e_eq is chosen such that the timescales for viscous stirring and aerodynamic gas drag on the planetesimals are in equilibrium. Although this approach assumes that these two mechanisms are in balance, there is nothing preventing planetesimal accretion from getting underway before the disk is sufficiently hot to be limited by gas drag. However, as we showed in the previous section, the initial dynamical state of the planetesimals does not seem to affect the outcome of accretion, so it is safe to assume that the resulting distribution of protoplanets would remain unchanged had we started with a colder disk. The viscous stirring timescale is given by Ida & Makino (1993) as

$$\begin{eqnarray}&&{t}_{{vs}}=\displaystyle \frac{\left\langle {e}^{2}\right\rangle }{d\left\langle {e}^{2}\right\rangle /{dt}}\approx \displaystyle \frac{1}{40}{\left(\displaystyle \frac{{{\rm{\Omega }}}^{2}{a}^{3}}{2{{Gm}}_{\mathrm{pl}}}\right)}^{2}\displaystyle \frac{4{m}_{\mathrm{pl}}\langle {e}^{2}{\rangle }^{2}}{{\rm{\Sigma }}{a}^{2}{\rm{\Omega }}},\end{eqnarray} \tag{ 12 }$$

where Ω, a, and e are the orbital frequencies, semimajor axes, and eccentricities of the individual planetesimals, respectively. In the Stokes regime, where the mean-free path of the gas is much smaller than the solid particles, the gas can be treated as a fluid and the drag timescale is given by Adachi et al. (1976) as

$$\begin{eqnarray}&&{t}_{s}=\displaystyle \frac{2{m}_{\mathrm{pl}}}{{C}_{D}\pi {r}_{\mathrm{pl}}^{2}{\rho }_{g}{v}_{g}},\end{eqnarray} \tag{ 13 }$$

where C_D is a drag coefficient of order unity, ρ_g is the local gas volume density, and v_g is the headwind velocity of the gas experienced by the planetesimal. The local gas volume density is given by

$$\begin{eqnarray}&&{\rho }_{g}=\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{\sqrt{2\pi }{h}_{g}}\exp \left[-{z}^{2}/\left(2{h}_{g}^{2}\right)\right],\end{eqnarray} \tag{ 14 }$$

where Σ_g is the gas surface density (taken to be 240 times the solid surface density (Hayashi 1981)), h_g = c_s /Ω is the local gas scale height, and z is the height above the disk midplane. The sound speed profile is given by ${c}_{s}=\sqrt{{k}_{{\rm{B}}}T(r)/\left(\mu {m}_{{\rm{H}}}\right)}$, where k_B is Boltzmann's constant, $T(r)={T}_{0}{\left(r/1\mathrm{au}\right)}^{-q}$, μ = 2.34, and m_h is the mass of a hydrogen atom. For a protoplanetary disk around a typical M star, T₀ = 148 K and q = 0.58 (Andrews & Williams 2005).

Finally, the headwind velocity of the gas, due to the fact that the gas disk is pressure supported, is given by

$$\begin{eqnarray}&&{v}_{g}={v}_{k}\left[1-\sqrt{{{qc}}_{s}^{2}/{v}_{k}^{2}}\right],\end{eqnarray} \tag{ 15 }$$

where v_k is the local Keplerian velocity (see Equation (4.30) of Armitage 2020). As in Section 4, the argument of perihelion ω, longitude of ascending node Ω, and mean anomaly M for the planetesimals are drawn from a uniform distribution ∈[0, 2π).

One should note that this choice for the gas disk profile almost certainly does not capture the wide range of possibilities in real planet-forming disks. On one hand, a larger initial gas surface density could act to completely remove solids via radial drift, rendering in situ accretion of solids impossible. On the other hand, a more tenuous gas disk might render aerodynamic drag forces completely unimportant. In this case, the random velocity of the initial planetesimals should be close to their mutual escape velocity. As we showed in Section 4, the initial dynamical state of the solids seems to have a very minimal effect on the final outcome of planetesimal accretion. In a similar vein to Hansen & Murray (2012), we choose to use an MMSN-like profile for the gas disk and instead vary the solid surface density profile to capture the range of mechanisms that might have acted to facilitate planetesimal formation in the first place.

In addition to the mutual gravitational forces, a Stokes drag force due to the gas disk is applied to each particle, following the prescription described in Section 2.2.1 of Morishima et al. (2010). For the initial mass planetesimals, the respective Stokes numbers (St = t_s Ω) at the inner and outer disk edges are roughly 2 × 10⁵ and 10⁷. For t_s ≫ 1, bodies are decoupled from the gas and are only weakly affected by it. Because t_s ∼ m^1/3, the Stokes number grows as planetesimal accretion proceeds, and the drag force plays an increasingly minor role. Although the aerodynamic gas drag is not expected to significantly alter the final protoplanet distribution, we include its effects here in order to be self-consistent with the initial conditions, which were constructed by balancing the effects of viscous stirring with gas drag.

In the case of the fdHi simulation, there are nearly one million particles, whose orbital periods vary by two orders of magnitude. Because the interaction timescales near the inner edge of the disk are exceedingly short, a fixed time step size would required a prohibitively large number of steps to follow planetesimal growth throughout the entire disk. For this reason, we use a multitiered time-stepping scheme, in which each particle is placed onto the nearest power-of-two time step based on its most recently calculated gravitational acceleration. This scheme is used in almost all works using ChaNGa, and it is common among large-scale simulation codes.

This more efficient scheme introduces two issues, however. First, momentum is not completely conserved when bodies switch time-step tiers. The error introduced becomes particularly severe for a particle on an eccentric orbit, whose perihelion and aphelion distances straddle a time-step boundary. For a large collection of particles, this problem manifests itself as the development of a V-shaped gap in the a–e plane, centered on the boundary itself. To correct this problem, we introduce a slightly modified time-stepping criterion, which is based on the expected gravitational acceleration of the particle at pericenter. Only in the case of a close encounter with another planetesimal (in which the acceleration is no longer dominated by the star) is the time step allowed to reduce based on the original instantaneous criterion.

A second issue is introduced when two particles on different time steps undergo a collision. As in the previous case, momentum is not completely conserved, because the most recent "kick" steps did not happen simultaneously for these bodies. Early in the simulation, we find that this problem tends to trigger runaway growth at the time-step boundaries first. This issue carries itself forward through the embryo-formation phase, and protoplanets tend to form at the boundaries. To correct this issue, we ignore collisions between bodies on different time steps early in the simulation. We find that preventing multi-time-step collisions until after the maximum mass grows by a factor of 10 prevents any artifacts from developing at the time-step boundaries, while also minimizing the number of "skipped" collisions. In the case of the fdHi simulation, only about 20 collisions out of an eventual 900,000 are ignored. To verify that this time-stepping scheme does not affect protoplanet growth, we tested an annulus of growing planetesimals with both fixed steps and our two-phase variable time-stepping scheme. The results of these tests are shown in Appendix.

The timescales for embryo formation depend on the chosen surface density profile, along with the local orbital timescale. Protoplanets form first at the inner edge of the disk, where the dynamical timescales are short. Growth proceeds in an inside-out fashion, with the outermost regions of the disk completing the protoplanet assembly phase last (as an example, see Figure 1 of Kokubo & Ida (2002)). This radial timescale dependence is not typically accounted for in planet formation simulations (a notable exception being Emsenhuber et al. (2021a, 2021b)), and it appears to be an important component to forming realistic Solar System analogs (Clement et al. 2020). As with the narrow-annulus simulations, we stop the integration once the masses of protoplanets in the outermost region of the disk reach a steady value. In Table 1, we summarize the outcomes of the four "full-disk" cases.

We show the final state of the "fdHi" simulation in Figure 6. In the top panel, the initial (contours) and final (points) state of the simulation are shown in the orbital period–eccentricity plane. The size of the points indicates the relative mass of the bodies. In the bottom panel, the final masses of the bodies (in units of feeding zone size $\tilde{b}$) are shown as a function of orbital period. The y-values in the bottom panel of Figure 6 are calculated by solving Equation (8) for $\tilde{b}$, and inputting the initial surface density and final particle mass into the expression. In other words, $\tilde{b}$ is describing the size of the annulus that must be cut out of the planetesimal disk in order to produce a protoplanet of the current mass. By plotting the derived value of $\tilde{b}$ as a function of orbital period, differences in the dynamical interactions at different locations of the disk are made more clearly visible. The feeding zone size $\tilde{b}=2\sqrt{3}$ permitted by bodies on circular, noninclined orbits (Nakazawa & Ida 1988) is shown by the horizontal dashed line. In typical oligarchic growth simulations (Kokubo & Ida 1998), protoplanets tend to space themselves apart by $\tilde{b}=10$, although it should be noted that they do not consume all of the planetesimals within this distance.

Zoom In Zoom Out Reset image size

A qualitative shift in the protoplanet and planetesimal distribution is visible inside of ∼60 days. Interior to this location, there are very few remaining planetesimals and the embryos formed have larger feeding zones. Exterior to the boundary, the residual planetesimal population is much more visible, and protoplanets more closely follow the $\tilde{b}=2\sqrt{3}$ line. This suggests that the transition between the low-α and high-α accretion modes seen in Section 4 is happening near this location.

In Section 4, we postulated that the increased importance of inelastic damping in the inner, non-oligarchic growth region of the disk should lower the overall eccentricity of the protoplanets there. This behavior is not immediately apparent in the top panel of Figure 6. In fact, the opposite appears to be true. There are, however, a couple of factors in the wide-disk simulations that could make this extra dynamical cooling mechanism difficult to see. First, the initial eccentricity distributions of the inner and outer disk are different because of the dependence of the viscous stirring and gas-drag timescales on orbital period. The mean eccentricity at the outer disk edge is 4× larger than at the inner disk edge. Additionally, the protoplanet formation timescales for the inner and outer disk are vastly different, making a comparison between these regions at the same moment in time somewhat inappropriate. A quick back-of-the-envelope calculation yields ${\left\langle {e}^{2}\right\rangle }^{1/2}=0.05$ for a population of ∼1M_⊕ bodies with a random velocity dispersion equal to their mutual escape velocity. It is therefore likely the case that the innermost protoplanets have had ample time to self-stir.

To ensure that the boundary seen around 60 days in orbital period is not simply a transient product of the inside-out growth throughout the disk, we examine the time evolution of σ/Σ, which compares the planetesimal and total solid surface density at multiple orbital periods. In Figure 7, the value of σ/Σ is plotted as a function of time in 10 orbital period bins, each with a width of 10 days. To determine whether the evolution of the planetesimal surface density behaves self-similarly across the disk, we normalize the time values in each bin by the local accretion timescale at the beginning of the simulation, which is given by

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}={\left(n{\rm{\Gamma }}v\right)}^{-1}={\left(\displaystyle \frac{{{\rm{\Sigma }}}_{0}{\rm{\Omega }}}{2{m}_{\mathrm{pl}}}{\rm{\Gamma }}\right)}^{-1},\end{eqnarray} \tag{ 16 }$$

where we have assumed the local number density of particles n is set by the surface density and the local scale height of the disk (see Section 2). The effective collision cross section is set by gravitational focusing and is given by Equation (5).

Zoom In Zoom Out Reset image size

The color of the curves indicate the orbital period bin that is being measured. From about 40 to 100 days in orbital period, the planetesimal surface density follows a similar trajectory as accretion proceeds. Interior to about 40 days, σ actually decays more slowly. In other words, growth is actually fueled less vigorously by planetesimals in this region. This highlights the fact that accretion proceeds in a qualitatively different way in the inner disk. For the outer disk, gravitational focusing tends to facilitate collisions between protoplanets and preferentially smaller bodies. At short period, however, all close encounters result in a collision, regardless of mass. In a rather counterintuitive fashion, planetesimals in the inner disk actually persist for longer. In Section 5.5, we examine the assembly history of the embryos and show that there is much less of a preference for planetesimal–embryo collisions at short period as well.

In the inner disk, this value asymptotes to zero as the planetesimal population entirely depletes. In the outer disk, dynamical friction between the embryos and planetesimals eventually throttles subsequent accretion and leaves ∼10% or more of the mass surface density as planetesimals. It should be noted that, in a typical oligarchic growth scenario, where protoplanets space themselves apart by 10 r_h and settle onto circular orbits (giving $\tilde{b}=2\sqrt{3}$), roughly 30% ($2\sqrt{3}/10\simeq 0.3$) of the planetesimals should remain out of reach of the protoplanets.

Next, we investigate how the resulting planetesimal and protoplanet distribution changes as we vary the initial solid surface density profile. The final orbital period–eccentricity state of the particles in the fdHi, fdHiSteep, fdHiShallow, and fdLo simulations are shown in Figure 8, with point sizes indicating the relative masses of the bodies. In all cases, the inner disk is largely depleted of planetesimals, while the outer disk contains a bimodal population of planetesimals and embryos, with a clear separation in eccentricity between the two. Despite having significantly different masses, the semimajor axis–eccentricity distributions of the planetary embryos formed in all simulations are remarkably similar. This is likely due to the fact that inelastic collisions play a more significant role where the solid surface density is highest, which offsets the fact that the initial bodies started off in a dynamically hotter state (due to the increased effectiveness of viscous stirring). The only exception to this is the fdLo simulation, where the resulting eccentricities are a couple orders of magnitude smaller. Inelastic damping likely plays an even more significant role here, due to the much larger masses of the initial planetesimals.

Zoom In Zoom Out Reset image size

In Figure 9, we plot the masses of the resulting protoplanets and planetesimals in all four simulations in units of $\tilde{b}$ (see Equation (8)). To make the trend between $\tilde{b}$ and orbital period more clear, we ran four more versions of the fdHi, fdHiSteep, and fdHiShallow simulations using different random number seeds and included these in the figure as well. As mentioned previously, $\tilde{b}=2\sqrt{3}$ (indicated by the horizontal dashed line) is the feeding zone width that a body on a circular orbit will have. In all four simulations, the feeding zone width exceeds the minimum value in the inner disk and approaches $2\sqrt{3}$ around ∼40–60 days. The orbital periods at which this transition occurs are quite similar between simulations, despite the vastly different solid surface density profiles used. This indicates that the boundary between accretion modes is driven entirely by the local value of α, and it also supports our conclusion that planetesimal accretion is largely complete everywhere in the disk.

Zoom In Zoom Out Reset image size

Further insight regarding the difference between the short- versus long-period accretion modes can be gained by looking at the growth history of the planetary embryos. Because all collisions are directly resolved by the N-body code, a lineage can be traced between each planetary embryo and the initial planetesimals. For the fdHi, fdHiSteep, and fdHiShallow simulations, protoplanets gain a factor of ∼10⁶ in mass relative to the initial planetesimals. For the fdLo simulation, this growth factor is nearly a thousand times smaller, which produces rather shallow and noisy collision histories. For this reason, we choose to exclude the fdLo simulation from our analysis in this section.

We begin by investigating the "smoothness" of the accretion events that give rise to each embryo. Drawing from a common technique used for cosmological simulations of galaxy formation, we divide growth events for a given body into "major" and "minor" mergers (Kauffmann & White 1993; Murali et al. 2002; L'Huillier et al. 2012). Here, we define minor events as any collision involving an initial mass planetesimal, while major events consist of a merger between any two larger bodies. In Figure 10, we retrieve the collision events for all bodies in all five iterations of the fdHi, fdHiSteep, and fdHiShallow simulations, and we plot the total mass fraction attained through minor merger (smooth accretion) events as a function of the final mass of the body. Here, we define a minor merger to be any collision involving a planetesimal with m < 100m₀. The color of the points indicates the final orbital period of the body. Beyond ∼20–30 days in orbital period, minor mergers make up a significant fraction of the final mass of a body. Interior to this, the smooth accretion fraction drops significantly and the mass contribution to minor mergers can vary by over an order of magnitude.

Zoom In Zoom Out Reset image size

The variation in smooth accretion fraction with mass for the short-period bodies suggests that the planetesimal and embryo populations interact differently than those in the outer disk. Exterior to the accretion mode boundary, the growing embryos continue to accrete planetesimals while avoiding each other as they near their final mass. Inside the boundary, however, any and all bodies collide with each other, and the occasional embryo–embryo collision tends to dominate growth and drive down the smooth accretion fraction. Gravitational scattering between embryos and planetesimals is a key ingredient for orbital repulsion (Kokubo & Ida 1998), and so a lack of gravitational scattering in the inner disk should prevent the embryos from settling onto neatly spaced, isolated orbits. As we showed in Figure 9, the embryos in the inner disk appear to reach well beyond the typical feeding zone size predicted by an oligarchic growth model. Figure 10 suggests that the extra mass here comes from mergers with the other nearby embryos.

Another line of evidence pointing to a lack of gravitational scattering and orbital repulsion in the inner disk can be seen in Figure 11. Here, we measure the initial orbital period distribution of bodies used to construct each embryo and calculate its mean $\left\langle {P}_{\mathrm{acc}}\right\rangle$. Point sizes indicate the relative final masses of bodies. As in the previous figure, we have included data from all five versions of the fdHi, fdHiSteep, and fdHiShallow simulations. For each body, this quantity is then compared with its final orbital period. In this figure, bodies that closely follow the gray dashed line still reside close to their initial feeding zones. On the other hand, bodies further from the dashed line must have experienced a strong gravitational scattering event during or after their accretion has completed. For all of the high surface density simulations, the accretion zones and present positions of the embryos appear to diverge beyond ∼20 days in orbital period.

Zoom In Zoom Out Reset image size

Coupled with the strong decrease in smooth accretion fraction for bodies in this region (seen in Figure 10), it appears as if the relative importance of collisions and gravitational scattering begins to shift around ∼20 days. For the shortest-period bodies, growth events are sudden and stochastic, often involving collisions between bodies of comparable mass. For longer-period bodies, a significant amount of growth is driven by accretion of smaller planetesimals. We postulate that this qualitative difference is driven by the role that embryo–embryo close encounters play in the inner and outer disk. In the inner disk, these encounters tend to result in a merger, which drives down the smooth accretion fraction. In the outer disk, these encounters tend to result in a scattering event that moves bodies away from their initial feeding zones. We find that the accretion zone shapes for the longer-period bodies are also much more smooth and unimodal, which suggests that scattering tends to occur after accretion has largely completed.

In all cases shown so far, the boundary between the oligarchic growth and the highly efficient short-period accretion region lies between 40 and 70 days in orbital period. As discussed in Section 2.2, the mode of accretion is set entirely by the local value of α, which scales with both distance from the star and the bulk density of the planetesimals (see Equation (9)). Because we chose to artificially inflate the collision cross section of the particles in our simulations by a factor of f, the bulk densities of the particles are reduced, and the accretion boundary is shifted outward. However, the scaling relation between α and ρ ($\alpha \sim {\rho }_{\mathrm{pl}}^{-1/3}$) can be used to predict where this accretion boundary should lie in a disk with realistic-sized planetesimals. The simulations presented in this paper use a collision cross-section enhancement factor of 6, which moves the boundary outward in orbital period by a factor of approximately 15. (For a fixed value of α, Equation (9) gives a ∼ r_pl ∼ f and therefore P_orbit ∼ f^−3/2.) One would therefore expect the accretion boundary to lie between 3 and 5 days in orbital period for 3 g cm⁻³ bodies.

Although a simulation with f = 1 is not computationally tractable, we can test whether the accretion boundary moves in the way we expect by modestly changing the value of f. In Figure 12, we compare the fdHi simulation to a nearly identical run using f = 4. In the top panel, we show the feeding zone width required for each particle to attain its present mass. As in Figure 9, we indicate the feeding zone size expected for oligarchic growth with a horizontal dashed line. In the bottom panel, the value of α as a function of orbital period is shown for 3 g cm⁻³ bodies with artificial radius enhancements of f = 1, 4, and 6. The horizontal dashed line indicates the empirical value of alpha below which the accretion mode switches to oligarchic. Comparing the top and bottom panels, the intersection of the feeding zone width seen in our simulations and the feeding zone width predicted by oligarchic growth matches well with the orbital period at which α ∼ 0.1 for both values of f. Also shown by the shaded region are the expected α values for realistic-sized bodies with ρ_pl between 1 and 10 g cm⁻³. Although the removal of the cross-section enhancement greatly reduces the size of the non-oligarchic region, it still should be expected to cover a portion of the disk where planetesimals might be expected to form (Mulders et al. 2018) for a wide range of ρ_pl.

Zoom In Zoom Out Reset image size

For the simulations presented in this work, every collision results in a perfect merger between pairs of bodies, with no loss of mass or energy. Although simpler and less computationally expensive to model, allowing every collision to produce a perfect merger might result in overly efficient growth, in particular in the innermost region of the disk where the encounter velocities are largest. Given that we have just shown that a distinctly non-oligarchic growth mode emerges in the inner disk when the collision timescale is short relative to the gravitational scattering timescale, one might be concerned that a more realistic collision model would act to lengthen the growth timescale enough for this condition to no longer be true. In the outer regions of the disk where oligarchic growth still operates, more realistic collision models have been shown to simply lengthen the timescale for planetary embryos to form (Wetherill & Stewart 1993; Leinhardt & Richardson 2005). A proper way to handle this would be to allow for a range of collision outcomes, based on a semianalytic model (see Leinhardt & Stewart 2012). However, resolving collisional debris, or even prolonging growth by forcing high-velocity pairs of bodies to bounce, is too expensive to model, even with ChaNGa.

To test whether a more restrictive collision model should alter the growth mode of the inner disk, we ran a smaller scale test using a more restrictive collision model. In this case, a collision can result in one of two outcomes: if the impact velocity is smaller than the mutual escape velocity of the colliding particles, defined as

$$\begin{eqnarray}&&{v}_{\mathrm{mut},\mathrm{esc}}=\sqrt{\displaystyle \frac{2G({m}_{1}+{m}_{2})}{{r}_{1}+{r}_{2}}},\end{eqnarray} \tag{ 17 }$$

where m₁, m₂ and r₁, r₂ are the masses and radii of colliding particles 1 and 2, then the bodies merge. For impact velocities larger than v_mut,esc, no mass is transferred, and the bodies undergo a completely elastic bounce. Because the accretion outcome is all or nothing, this model should restrict growth more than a partial accretion model (Leinhardt & Stewart 2012). Below, we will show that the bounce–merge model does not meaningfully affect the outcome of the inner disk's planetesimal accretion phase, and so a more realistic partial accretion model should do the same.

To compare the outcome of the two collision models, we have chosen to use the initial conditions from the fdLo simulation, but have truncated the disk beyond 3 days in orbital period. This offsets the increased computational cost of the more restrictive collision model, while still allowing the disk to evolve in the region where mergers would be most difficult to achieve. For the initial conditions we have chosen, the typical encounter velocity (defined by ${v}_{\mathrm{enc}}={\left\langle {e}^{2}\right\rangle }^{1/2}{v}_{k}$, where v_k is the local Keplerian velocity) is about 25% larger than v_mut,esc. Because the encounter velocities follow a Gaussian distribution, there should still be a small subset of collisions that still meet the merger criteria to occur early on. In addition, v_mut,esc becomes larger as the bodies grow, and the merger criteria should become easier to meet as the system evolves. For these reasons, one would expect the inhibition of growth due to the more restrictive collision model to be temporary.

In Figure 13, we compare the outcomes of the simulations, one with mergers only (shown in orange) and one with the bounce–merge model (shown in green). The blue points in the top panel show the initial conditions used for both cases. Although the bounce–merge simulation takes much longer to reach the same phase of evolution, the resulting orbital properties are indistinguishable from the merger-only case. Performing a Kolmogorov–Smirnov test on the two mass distributions yields a p-value of 2 × 10⁻⁵, which tells us that the two mass distributions are quite firmly statistically different. If we remove the initial mass planetesimals, a KS test yields a p-value of 0.1, which suggests that the distributions are statistically similar. Because the remaining planetesimals only make up about 0.1% of the total mass of the disk, we conclude that the embryo populations are nearly indistinguishable, while the bounce–merge model produces a small amount of residual planetesimals.

Zoom In Zoom Out Reset image size

To investigate the differences in growth between the two collision models early on, we show the time evolution of the ratio between the maximum and mean mass in Figure 14. In both cases, this ratio first increases, which indicates that runaway growth still operates, regardless of the collision model used. In the bounce–merge case, the mass ratio peaks at a higher value, while also undergoing a longer runaway growth phase. This suggests that the mass distribution becomes much less unimodal during this growth process, but as Figure 13 shows, this does not affect the resulting embryos or allow for a residual planetesimal population.

Zoom In Zoom Out Reset image size

As a final note, Childs & Steffen (2022) found that a more realistic collision model also enhanced radial mixing in their simulations. Upon calculating the planetesimal accretion zones using the same method as was done to produce Figure 11, we find that the embryos in the bounce–merge simulation annulus do have modestly wider accretion zones than those produced in the merger-only simulation.

To date, there have been no other studies of planetesimal accretion with such a large value of α. Typically, it is assumed that α ≪ 1 (e.g., Lithwick 2014), which is certainly true for material at and beyond the Earth's orbit. However, a value of α = 1 corresponds to the Roche limit of a three-body system, and so one might wonder this high-α accretion mode might be relevant for a circumplanetary accretion. There is a collection of previous works that use N-body methods to examine in situ satellitesimal accretion (Ida et al. 1997, 2020; Kokubo et al. 2000; Richardson et al. 2000), although some of these simulations involve a complex interaction between spiral density waves formed inside of the Roche limit and the material exterior to it, making the dynamics driving accretion distinctly nonlocal, in contrast to what we have presented in this work. Ida et al. (1997) were able to form 1–2 large moons just exterior to the Roche limit, depending on the extent of the disk with very little satellitesimal material left over. The widest disk they modeled extended out to α = 0.5. Qualitatively, this result is very similar to the short-period planetesimal accretion mode observed in our simulations. Ida et al. (2020) modeled a much wider satellitesimal disk, which extends out to about α ≈ 0.05. Inside to the α = 0.1 accretion boundary (which lies near 15 planetary radii in Figure 1 of Ida et al. 2020), bodies grow beyond the isolation mass, while the opposite is true on the other side of the boundary. In addition, a residual population of satellitesimals is still present beyond the boundary, which suggests that oligarchic growth is indeed operating only on the far side.

Presently, the implications that this non-oligarchic accretion mode has for the formation of short-period terrestrial planets, and whether the accretion boundary would leave any lasting imprint on the final orbital architecture, remain unclear. The extreme efficiency of planetesimal accretion at the inner edge of the disk suggests that no residual populations of small bodies should be expected to exist here. A crucial point that our results do highlight is that the initial conditions used for most late-stage planet formation simulations are overly simplistic. Clement et al. (2020) recently simulated planetesimal accretion in a disk extending from the orbit of Mercury to the asteroid belt and found that the disk never reaches a state in which equally spaced, isolation-mass embryos are present everywhere simultaneously. Instead, different annuli reach a "giant impact" phase at different times, preventing the onset of a global instability throughout the entire disk, as is common in classic terrestrial planet formation models (Chambers & Wetherill 2001; Raymond et al. 2009).

To connect these accretion modes to the final orbital architecture, and to ultimately determine what implications an in situ formation model has for the growth of STIPs, we will evolve the final simulation snapshots presented here with a hybrid-symplectic integrator for Myr timescales. The final distribution of planets formed, along with composition predictions generated by applying cosmochemical models to our initial planetesimal distributions and propagating compositions through the collision trees, will be examined in a follow-up paper.