Introduction
Recent advances in computational technology have had major impacts on the studies of the formation and evolution of planetary systems. Fast and powerful computers have enabled high-resolution simulations, and sophisticated codes have allowed collisions to be modeled more accurately (e.g., Burger et al., Reference Burger, Bazsó and Schäfer2020; Reinhardt et al., Reference Reinhardt, Meier, Stadel, Otegi and Helled2022). Collectively, these developments have elevated the field to a new level where now the majority of the underlying assumptions that have primarily been developed to circumvent computational complexities can be removed, and almost any physical processes relevant to planet formation can be directly included in the simulations. The latter is a powerful capability that allows for examining the consequences of different physical effects thereby developing a pathway to building a realistic formation model from grounds up. On that front, it is of utmost importance to determine the requirements that any model of planet formation needs to satisfy in order for its results to be realistic.
In the past few years, the idea of removing simplifying assumptions, combined with requiring models to meet criteria for producing realistic outcomes raised a demand from some members of the community to justify its worth and significance. It is important to note that any action to respond to this demand needs to not only identify the requirements for a realistic model, but also to demonstrate, for each criterion, the scientific rational that would justify its consideration. To do so, the foundations of planet formation need to be revisited, and all their initial conditions and underlying assumptions need to be examined. The goal of this article is to address the above in the context of planetesimals runaway growth.
One of the most outstanding assumptions, that despite its unrealistic nature, has been used repeatedly for almost three decades, is the assumption of “inflated planetesimals.” As we will explain in the next section, to avoid long integration times associated with the collisional growth of small bodies, in almost all simulations of planetesimal growth, the rate of the collisions of planetesimals is artificially enhanced by increasing their initial sizes, forcing their growth to happen in shorter times. Table 1 shows an up to date list of the studies that have used inflated planetesimals along with their regions of simulation, integration resolution (number of planetesimals), and their radius-enlargement factor
$\left( f \right)$
.
While this unphysical approach does not affect the mechanics of the occurrence of runaway growth, and may even be useful in demonstrating some of the underlying physics of the process, it will have profound effects on the formation of planetary embryos and the subsequent evolution of the system. For instance, in our solar system, the conventional wisdom is that terrestrial planets formed through a two-stage process: the runaway growth, during which larger planetesimals swept up smaller ones and grew to Moon- to Mars-sized planetary embryos, followed by giant impacts among the latter bodiesFootnote 1 . Because runaway growth is a natural consequence of the dynamical evolution of planetesimals, when studying the formation of terrestrial planets, it is customary to, a priori, assume its outcome and use the latter as the initial condition for simulating the state of giant impacts. However, during this time (i.e., runaway growth), the dynamics of the system is strongly affected by the perturbation of its rapidly growing bodies. As the growth of these objects has been artificially enhanced (by increasing the sizes of their seed planetesimals), their perturbation becomes stronger than it would have been, had their planetesimals not been inflated. The latter deviates the system from its natural evolution within its natural timescale, causing it to take a dynamical path that is dictated to it by the perturbation of its artificially hyper-sized bodies. This deviation of the system from its natural evolution in the runaway growth phase propagates to the phase of giant impacts, rendering the outcome of the simulations of the late stage of terrestrial planet formation unreliable. For instance, because the accretion zone of the artificially enlarged runaway growth bodies increases with their large masses, the planetary embryos that are formed by the growth of these objects become larger than they would have been without the radius inflation. This larger masses of planetary embryos profoundly affects the masses and orbital architecture of the final terrestrial planets as it increases the spacing among the embryos (an important factor in assembling protoplanetary disks at the beginning of the giant-impacts simulations) and decreases the effect of dynamical friction by scattering a large number of smaller bodies out of the system.
Table 1. An up to date list of all
$N$
-body simulations of planetesimals growth showing the region of the simulation, integration resolution (number of planetesimals), and the radius-enlargement factor
$\left( f \right)$
To ensure that the final system is free from these biases, simulations have to be carried out using uninflated planetesimals. The results of such simulations, combined with the mechanics of the process will have implications that point to and/or manifest themselves as requirements for obtaining realistic outcomes. To better understand this process, its worthiness, and the significance of its implications, it is necessary to develop a deep understanding of the history of modeling planetesimal growth and more importantly, the thought-process through which the field has progressed. In this article, we aim to achieve the latter by demonstrating how the requirements for a realistic model emerge organically from the natural evolution of the system. To do so, we begin by presenting an in-depth review of the history of planetesimal growth, and continue by carrying out new simulations and analyzing the implications of the results in connection with the results in the literature. Because, traditionally, planetesimal growth has been studied as an introduction to terrestrial planet formation, we follow the same approach and focus our review and simulations within the same context. However, our analyses are completely general and can be applied to the formation of the cores of the giant planets as well.
That the collision among small planetesimals may result in the rapid formation of much larger bodies was first noted by Greenberg et al. (Reference Greenberg, Wacker, Hartmann and Chapman1978). Using the particle-in-a-box approach, these authors found that in a swarm of km-sized planetesimals, a few locally large bodies may dominate the process and rapidly grow to a few hundred km in size. Prior to that study, Safronov (Reference Safronov1969) had shown that planetesimal growth would proceed in an orderly manner meaning that large planetesimals gradually grow by sweeping up smaller ones while maintaining their mass ratios. Orderly growth was also reported in later studies by Nakagawa et al. (Reference Nakagawa, Hayashi and Nakazawa1983), Ohtsuki et al. (Reference Ohtsuki, Nakagawa and Nakazawa1988) and Hayakawa et al. (Reference Hayakawa, Mizutani, Kawakami and Takagi1989) who advanced the work of Safronov (Reference Safronov1969) by including the effect of gas drag.
The viability of runaway growth was first examined by Wetherill and Stewart (Reference Wetherill and Stewart1989). Because the underlying motivation behind their study was to understand the formation of Earth, these authors considered a swarm of planetesimals between 0.98 and 1.02 AU, and following Safronov (Reference Safronov1962, Reference Safronov1969), distributed them according to the mass distribution function
$${{{dN( m )}} \over{dm}} {\propto {e^{ - m/{m_0}}}},$$
(similar distribution function had also been used by Greenberg et al., Reference Greenberg, Wacker, Hartmann and Chapman1978 and Nakagawa et al., Reference Nakagawa, Hayashi and Nakazawa1983). In this equation,
$N\left( m \right)$
is the number density of planetesimals with mass
$m$
, and
${m_0}$
is the mean mass of the planetesimals distribution at the time
$t=0$
. Wetherill and Stewart (Reference Wetherill and Stewart1989) showed that, under the conditions considered by Safronov (Reference Safronov1969) and Nakagawa et al. (Reference Nakagawa, Hayashi and Nakazawa1983) (that is, when the mutual interactions between planetesimals are neglected), the results obtained by those authors were correct. However, when simulations include planetesimals interactions, the equipartition of energy (Stewart and Kaula, Reference Stewart and Kaula1980), resulted from, for instance, dynamical friction, viscous stirring, the existence of large bodies, low-velocity gravitational focusing, gas drag, and fragmentation, will cause a few large bodies to dominate the growth and, in a relatively short time, reach sizes of the order of several hundreds to a few thousands times those of the initial planetesimals. Wetherill and Stewart (Reference Wetherill and Stewart1989) also carried out simulations similar to those in Greenberg et al. (Reference Greenberg, Wacker, Hartmann and Chapman1978). They showed that unlike the results reported by those authors, their simulations show no signs and indications of runaway growth. Wetherill and Stewart (Reference Wetherill and Stewart1989) attributed the appearance of runaway growth in the simulations by Greenberg et al. (Reference Greenberg, Wacker, Hartmann and Chapman1978) to the inaccuracy of the numerical method used by those authors (Wetherill, Reference Wetherill1990).
In principle, gravitational focusing, combined with the equipartition of energy is the main reason that runaway growth occurs (Ohtsuki et al., Reference Ohtsuki, Ida, Nakagawa, Nakazawa, Levy, Lunine and Tucson1993; Kokubo and Ida, Reference Kokubo and Ida1996). In a system of planetesimals with different sizes, the equipartition of energy manifests itself in the form of dynamical friction. That is, the interaction between a large planetesimal and a population of smaller ones reduces the relative velocity of the large body by lowering its eccentricity and inclination (Ida, Reference Ida1990; Ida and Makino, Reference Ida and Makino1992a, b, Reference Ida and Makino1993). This low relative velocity enhances the collisional growth of the large planetesimal, causing it to dominate the growth process in its surrounding.
Although, as a proof of concept, simulations by Wetherill and Stewart (Reference Wetherill and Stewart1989) successfully demonstrated the appearance and progression of runaway growth, their underlying assumptions limited the applicability and generalization of their results. For instance, these authors did not include the effect of the gravity of the Sun on the collision of objects and assumed, a priori, that the spatial velocities of planetesimals would follow a triaxial Gaussian distribution. These limitations motivated other researchers to revisit runaway growth and carry out new simulations (Ida, Reference Ida1990; Ohtsuki and Ida, Reference Ohtsuki and Ida1990; Barge and Pellat, Reference Barge and Pellat1991; Spaute et al., Reference Spaute, Weidenschilling, Davis and Marzari1991; Kolvoord and Greenberg, Reference Kolvoord and Greenberg1992; Wetherill and Stewart, Reference Wetherill and Stewart1993). However, all these simulations, as well as those of Wetherill and Stewart (Reference Wetherill and Stewart1989) still suffered from one fundamental inconsistency: In all these studies, growth was modeled through a statistical approach using the coagulation equation in the context of kinetic theory. This approach requires a uniform and homogeneous distribution of planetesimals both at the start and during a simulation. The issue is that in reality, the perturbation of the runaway growth bodies disturbs the distribution of planetesimals in their surrounding and breaks its uniformity (Ida and Makino, Reference Ida and Makino1993; Tanaka and Nakazawa, Reference Tanaka and Nakazawa1994).Footnote 2 In all above studies, either these non-uniformities were ignored, or the distribution of planetesimals was held uniform, artificially.
It is fundamentally important to note that the most realistic approach to studying planetesimal growth is simulating their collisions by direct integration of their orbits. Aarseth et al. (Reference Aarseth, Lin and Palmer1993) were first to use this approach in the context of planetesimal growth.Footnote 3 In a study on examining the predictions of their earlier work where Palmer et al. (Reference Palmer, Lin and Aarseth1993) had presented an analytical treatment for the dynamical relaxation of a disk of planetesimals, these authors integrated the orbits of 100 equal-mass and equal-size small bodies uniformly distributed around 1 AU. We recall that triggering runaway growth requires non-uniformities in the size distribution of planetesimals, meaning that one or a few of these bodies must have larger sizes than those in their surrounding. Aarseth et al. (Reference Aarseth, Lin and Palmer1993) considered a disk of equal-mass and equal-size planetesimals (which they referred to as the least favorite environment for the runaway growth) to examine whether and how runaway growth would appear and proceed. They showed that when allowing collisions with low relative velocities to result in coalescence (larger relative velocities were considered to result in bouncing), some planetesimals do in fact show signs of the runaway growth and increase their sizes to several times their original radii. However, because of the low resolution and short time of integrations (which were due to the limitations in their computational capabilities), the simulations by these authors did not show the sharp increase that is expected to be observed in the sizes of the runaway growth bodies. For that reason, their results stayed rather unnoticed.
The first (relatively) high-resolution simulations of runaway growth were carried out by Kokubo and Ida (Reference Kokubo and Ida1996) and were later expanded by (Kokubo and Ida, Reference Kokubo and Ida1998, Reference Kokubo and Ida2000, Reference Kokubo and Ida2002, see Kokubo and Ida, Reference Kokubo and Ida2012 for a review). These authors simulated the growth of 3000 equal-mass planetesimals randomly distributed between 0.98 and 1.02 AU by integrating their orbits and allowing them to collide with one another. The initial eccentricities and inclinations of planetesimals were chosen following the distributions used in Ida and Makino (Reference Ida and Makino1992a). To avoid computational complexities due to the increase in the number of the bodies resulted from the addition of collisional fragments, these authors considered growth by the way of perfect-merging.
As mentioned above, the onset of runaway growth requires one or a few planetesimals to be substantially larger than others. Because in the system considered by Kokubo and Ida (Reference Kokubo and Ida1996, Reference Kokubo and Ida1998, Reference Kokubo and Ida2002) all planetesimals had equal sizes, it would take a long time for their population to naturally produce a few large bodies on its own. To speed up the integrations, these authors artificially increased the sizes of planetesimals by 5–10 folds at the beginning of simulations. Kokubo and Ida (Reference Kokubo and Ida1996) argued that, based on their approximated analytical method, this artificial increase in the planetesimals sizes affects only the accretion time scale and plays no role on whether the runaway growth occurs or not. The results by these authors showed that as planetesimals collide and grow, the locally large bodies undergo runaway growth in their surrounding, increasing their sizes by two orders of magnitude. The perturbation from these bodies disturbs the orbits of planetesimals in their vicinity causing non-uniformities to develop in planetesimals distribution. Kokubo and Ida (Reference Kokubo and Ida1996) showed that the runaway growth appears only when simulations are 3-dimensional, and results in the formation of planetary embryos within the first 20,000 years (similar timescale has also been reported by Kokubo and Ida, Reference Kokubo and Ida1998). When co-planar systems are considered, growth through perfect-merging does not result in the runaway growth. Instead, similar to the results of Safronov (Reference Safronov1969) and Nakagawa et al. (Reference Nakagawa, Hayashi and Nakazawa1983), it proceeds in an orderly manner.
The development of computational models capable of resolving collisions, combined with the advances in computer technology prompted researchers to improve the direct simulations of runaway growth by (i) increasing the resolution of
$N$
-body integrations, (ii) expanding the spatial range of simulations, and (iii) including a certain number of collisional fragments. For instance, Richardson et al. (Reference Richardson, Quinn, Stadel and Lake2000) simulated the collisional growth of
${10^6}$
uninflated planetesimals via perfect-merging in the region of 0.8–3.8 AU considering the perturbation of the present-day outer planets. These authors showed that, although their integrations were so slow that it reached only 1000 years after 200 wall-clock hours, at which point they stopped the integrations, some of the planetesimals still reached the state of runaway growth. However, due to the short time of integration, the largest planetesimals did not reach the protoplanetary sizes and grew to only nine times their initial radii. To demonstrate the reliability of their results, Richardson et al. (Reference Richardson, Quinn, Stadel and Lake2000) reproduced the results of Kokubo and Ida (Reference Kokubo and Ida1998) by integrating the orbits of 4000 equal-mass planetesimals within the annulus of 0.96 to 1.04 AU. These authors inflated the initial planetesimals by 6 fold and showed that, in agreement with the results of Kokubo and Ida (Reference Kokubo and Ida1998), their integrations were able to produce planetary embryos within the first 20,000 years. Kokubo and Ida (Reference Kokubo and Ida2002) also presented high-resolution integrations of planetesimal growth where they integrated the orbits of
$\left( {5 - 10} \right) \times {10^4}$
planetesimals distributed from 0.5 to 1.5 AU. These authors too increased the initial sizes of planetesimals by 6–10 folds and showed that when the growth is modeled via perfect-merging, the first set of protoplanetary bodies with masses of 100 times larger than the initial planetesimals appears in 50,000 yearsFootnote
4
.
High-resolution simulations of planetesimal growth were also carried out by Morishima et al. (Reference Morishima, Schmidt, Stadel and Moore2008) and Barnes et al. (Reference Barnes, Quinn, Lissauer and Richardson2009). To examine the effect of remnant planetesimals on the orbital eccentricities and growth of planetary embryos, Morishima et al. (Reference Morishima, Schmidt, Stadel and Moore2008) integrated the orbits of up to 5000 equal-mass planetesimals in two annuli centered at 1 AU, one with a width of 0.3 AU and one with a width of 0.5 AU. These authors inflated the initial sizes of planetesimals by a factor of 4.3 and showed that when growth is modeled through perfect-merging, their systems reach the state of runaway growth at ~104 years. Barnes et al. (Reference Barnes, Quinn, Lissauer and Richardson2009) considered
${10^5}$
km-sized planetesimals at 0.4 AU, and integrated their orbits for 89 to 135 years without inflating their initial sizes. These authors demonstrated that, while some of the bodies did in fact grow by two orders of magnitude, runaway growth did not start as these bodies did not detach themselves from the rest of the population. We believe this lack of detachment (and with the same token, the non-appearance of runaway growth) is due to the short time of integrations by these authors.
The first
$N$
-body integrations of planetesimal growth that included collisional fragments were carried out by Leinhardt and Richardson (Reference Leinhardt and Richardson2005) and were later followed by Leinhardt et al. (Reference Leinhardt, Richardson, Lufkin and Haseltine2009).Footnote
5
These authors simulated the growth of
${10^4}$
planetesimals (initially inflated in size by 6 folds) in a 1 AU-wide area around the orbit of Earth and included collisional fragments derived from the collision catalogs of Leinhardt et al. (Reference Leinhardt, Richardson and Quinn2000) and Leinhardt and Richardson (Reference Leinhardt and Richardson2002). Following those works, Bonsor et al. (Reference Bonsor, Leinhardt, Carter, Elliott, Walter and Stewart2015), Leinhardt et al. (Reference Leinhardt, Robinson, Carter and Lines2015) and Carter et al. (Reference Carter, Leinhardt, Elliott, Walter and Stewart2015) carried out similar simulations where, increasing the initial sizes of planetesimals by 6 folds, these authors integrated the orbits of
${10^4} - {10^5}$
bodies in the region between 0.5 and 1.5 AU. The growth in these studies were simulated using the accretion equation of Leinhardt and Richardson (Reference Leinhardt and Richardson2005), and collisional fragments were included using the analytic predictions of Leinhardt and Stewart (Reference Leinhardt and Stewart2012) and Stewart and Leinhardt (Reference Stewart and Leinhardt2012). The results of all these studies broadly matched those of previous simulations demonstrating that the inclusion of fragments does in fact affect the evolution of planetesimals and runaway growth, although the effect may be minor.
Most recently, Clement et al. (Reference Clement, Kaib and Chambers2020) studied the collision and evolution of planetesimalsFootnote
6
. These authors considered five narrow annuli between 0.5 and 3 AU, and integrated the orbits of 5000 uninflated planetesimals in each annulus. Using an
$N$
-body integrator within a GPU environment, these authors showed that, when the effect of gas drag is mimicked by including the analytic gas-disk model of Morishima et al. (Reference Morishima, Stadel and Moore2010), runaway growth appears within the first 80,000 years. As we will demonstrate later, given the uninflated sizes of the planetesimals, this time of the appearance of runaway growth is too short, even with the effect of gas drag included (gas drag will enhance the rate of growth).
Using a clever approach, Clement et al. (Reference Clement, Kaib and Chambers2020) expanded their annuli integrations to the entire region of 0.48–1.65 AU (as one large and continuous area) and showed that, when the effects of giant planets are included, while as expected, the perturbation of these planets mainly affects the outer region of the asteroid belt, it is still large enough to have a moderate effect on the planetesimals orbits in the inner regions as well (see Haghighipour and Scott, Reference Haghighipour and Scott2012, for more details on the effects of giant planets on planetesimals and planetary embryos in the terrestrial planet region). We will return to this topic and the work of these authors in more details in section “Implications of the results.”
At the time of this writing, the highest resolution simulations of planetesimal growth were those of Wallace and Quinn (Reference Wallace and Quinn2019). These authors integrated a set of low-resolution (4000) and a set of high-resolution (
${10^6}$
) planetesimals in a narrow annulus between 0.94 and 1.04 AU for 20,000 years. In order to demonstrate the reliability of their integrations by comparing their results with previous studies, these authors set the initial eccentricities and inclinations of their bodies similar to those in Ida and Makino (Reference Ida and Makino1992a), and increased the initial planetesimal sizes by 6 folds. The 20,000 year integration time was chosen as the shortest time that would allow the authors to compare their results with those of Kokubo and Ida (Reference Kokubo and Ida1998) while their extremely slow high-resolution simulations would finish in a reasonable time. Results showed that, in agreement with the work of Kokubo and Ida (Reference Kokubo and Ida1998), simulations did in fact reach the state of runaway growth within ~2000 years. They also demonstrated that the mean-motion resonances between growing planetary embryos (i.e., Moon- to Mars-size bodies) and planetesimals play an important role in the collisional growth of these bodies, resulting in the appearance of a bump in their final mass distribution (see the bottom-left panel of Figure 3 in Wallace and Quinn, Reference Wallace and Quinn2019, and the bottom-left panel in Figure 4 of this study).
Before we continue, we would like to note that the studies of planetesimal accretion and runaway growth have not been limited to the direct integration approach. Although
$N$
-body integrations present the most precise treatment of the growth and dynamical evolution of planetesimals, because realistic scenarios require integrating billions of bodies, which have been and still are beyond the capability of computational resources, many scientists studied planetesimal growth using statistical methods. For instance, Kenyon and Luu (1998, Reference Kenyon and Luu1999a, b), Kenyon and Bromley (Reference Kenyon and Bromley2001, Reference Kenyon and Bromley2002a, b), Ohtsuki et al. (Reference Ohtsuki, Stewart and Ida2002) and Kenyon and Bromley (Reference Kenyon and Bromley2004a, b) used the particle-in-a-box approach of Safronov (Reference Safronov1969) and studied planetesimal growth at the inner and outer parts of the solar system, as well as in extrasolar planets. To overcome the breakdown of this approach due to the appearance of non-uniformities, Weidenschilling et al. (Reference Weidenschilling, Spaute, Davis, Marzari and Ohtsuki1997), Kenyon and Bromley (Reference Kenyon and Bromley2006) and Bromley and Kenyon (Reference Bromley and Kenyon2006, Reference Bromley and Kenyon2011) developed hybrid methods by combining
$N$
-body integrations with an improved version of the coagulation scheme. Collectively, results of these studies showed that, contrary to direct
$N$
-body integrations, runaway growth appears within 10-100 Myr, with this time increasing with the distance from the central star.
Ormel et al. (Reference Ormel, Dullemond and Spaans2010a) and Ormel and Okuzumi (Reference Ormel and Okuzumi2013) also studied runaway growth in the context of giant planet formation. Ormel et al. (Reference Ormel, Dullemond and Spaans2010a) developed a statistical approach that would use a full Monte Carlo coagulation-fragmentation scheme, and showed that their approach preserves the individual nature of particles while it treats a large number of them statistically. In a later study, Ormel et al. (Reference Ormel, Dullemond and Spaans2010b) used the results of this statistical approach and refined the conditions for transition from runaway growth to oligarchic growth. These authors showed that a state consistent with the runaway growth appears in the mass distribution of planetesimals at about
${10^5}$
years. Ormel and Okuzumi (Reference Ormel and Okuzumi2013) developed a semi-analytical model through which they studied the effect of dead-zones in turbulent disks on the onset of runaway growth. These authors found that in the context of giant planet formation, the minimum size of a planetesimal to begin the runaway growth is 100 km or larger, and that the growth of such bodies will overrun the life time of the nebular gas.
On the analytical fronts, Kobayashi et al. (Reference Kobayashi, Tanaka, Krivov and Inaba2010) studied planetesimal growth and developed a formula for the radius of runaway growth bodies. These authors defined their analytical radius in such a way that when used in studying planetesimal growth, the results would reproduce those of previous numerical studies, in particular those of Inaba et al. (Reference Inaba, Tanaka, Nakazawa, Wetherill and Kokubo2001). Results of the numerical integrations by these authors demonstrated that, when considering growth via perfect-merging, their defined version of runaway growth appears in
$\left( {1 - 4} \right){10^5}$
years. Using the statistical code of Inaba et al. (Reference Inaba, Tanaka, Ohtsuki and Nakazawa1999, Reference Inaba, Tanaka, Nakazawa, Wetherill and Kokubo2001), Kobayashi et al. (Reference Kobayashi, Tanaka, Krivov and Inaba2010) also showed that their analytical results match those obtained from statistical analysis. The analytical and numerical approaches of Kobayashi et al. (Reference Kobayashi, Tanaka, Krivov and Inaba2010) have been used in the study of debris disks (Kobayashi and Löhne, Reference Kobayashi and Löhne2014), dynamical evolution of planetesimals in gaseous disks (Kobayashi, Reference Kobayashi2015), runaway growth in turbulent disks (Kobayashi et al., Reference Kobayashi, Tanaka and Okuzumi2016), and the growth of gas-giant planets (Kobayashi and Tanaka, Reference Kobayashi and Tanaka2018). Recently Walsh and Levison (Reference Walsh and Levison2019) studied the growth of planetesimals to terrestrial planets where they used tracer particles to simulate the growth of planetesimals and planetary embryos.
Although the above-mentioned indirect approaches, as well as those prior to the work of Aarseth et al. (Reference Aarseth, Lin and Palmer1993) have greatly contributed to our understanding of planetesimal growth, they are not discussed in this study. As mentioned earlier, direct integration is the most realistic approach to planetesimal growth, and for that reason, it allows for studying the underlying physical processes that actually affect the growth and evolution of planetesimals. We, therefore, maintain our focus on direct integrations, using real-size planetesimals, and continue by examining the implications of the results for developing a realistic and self-consistent approach to modeling planet formation.
As mentioned before, the reason for inflating planetesimals is to shorten the integration time. Although it has been known that this approach is not realistic, it has been argued that because, in general, enlarging planetesimals radii does not disturb the effectiveness of gravitational focusing, it will not affect the occurrence of runaway growth. In other words, as long as the study is qualitative, or for understanding the underlying physics of planetesimal growth, or for the purpose of proving a concept, the error associated with inflating planetesimals may not be relevant. However, when the study is quantitative, and when the purpose is to develop formation models that are capable of making realistic predictions, the approach needs to be realistic. In such cases, unrealistic initial conditions will not allow simulations to portray the natural evolution of the system. For instance, as shown in Figure 1 of Bonsor et al. (Reference Bonsor, Leinhardt, Carter, Elliott, Walter and Stewart2015), the time of growth increases by almost two orders of magnitude from
$5 \times {10^4}$
years to 1.8 Myr for the appearance of the runaway growth, and from
$2 \times {10^5}$
years to 7.2 Myr for the formation of planetary embryos when the initial set-up is changed from using 6-fold inflated planetesimals to objects with no artificially expanded radii. As shown by Levison and Agnor (Reference Levison and Agnor2003), Haghighipour and Scott (Reference Haghighipour and Scott2012), and Haghighipour and Winter (Reference Haghighipour and Winter2016), this time falls within the timescale that the gravitational perturbation of (growing) Jupiter disturbs the protoplanetary disk, especially within the time that the secular resonance of Jupiter begin to appear and disturb the region around 1 AU.
As the first step toward developing a realistic planet formation model, it is fundamentally important to revisit runaway growth and examine its consequences using initial conditions that allow for the natural evolution of the system. The rest of this article has been devoted to this purpose. As mentioned earlier, our approach is to revisit runaway growth by redoing some of the simulations and reviewing the results by combining them with those in the literature. When carrying out simulations, we maintain focus on principle concepts and take a conservative approach: We consider a system consisting of equal-size and equal-mass planetesimals. It is understood that such distribution of planetesimals is not entirely natural (a natural distribution will include planetesimals of different sizes and masses). However, quoting Aarseth et al. (Reference Aarseth, Lin and Palmer1993), as “the least favorable environment for runaway growth,” such a system will give us a less biased view into the process of the formation of planetary embryos. Also, to avoid complications due to the proper handling of collisions, we consider growth through perfect-merging. Although unrealistic, this scenario presents the most efficient mode of growth meaning that, any conditions and requirements that appear when using perfect-merging will only be enhanced when the collisions between planetesimals are resolved realistically. We refer the reader to Haghighipour and Maindl (Reference Haghighipour and Maindl2022) for a detailed analysis of the errors due to perfect-merging, and the remedy for them. Finally, we do not enlarge any planetesimal at the beginning of the integrations and to maintain focus on the underlying physics of runaway growth and its timescale, we do not include additional physical processes such as the perturbation of the growing giant planet(s) (these effects are included in a subsequent study), the damping effect of the nebular gas drag, the increase in the number of bodies due to fragmentation, and the dynamical friction due to the collisional debris (the dynamical friction due to the bi-modal mass of the disk is implicitly included).
Revisiting planetesimal growth
As mentioned earlier, to identify the requirements for a realistic model, the foundations of planet formation have to be revisited. Because in the context of planetesimal growth, many of the fundamental studies are more than 25 years old, and because those studies were carried out in contexts other than identifying requirements for realistic models, we revisited this topic by performing those fundamental simulations, but without inflating planetesimals. In the rest of this article, we will present the results of these simulations in the context of all previous investigations and discuss their implications both in connection to one another and in connection with previous studies.
As in many of the previous
$N$
-body simulations of planetesimal growth, we followed Kokubo and Ida (Reference Kokubo and Ida1996) and randomly distributed approximately 3000 equal-mass planetesimals around 1 AU. Several studies, such as those cited here, have demonstrated that a population of 3000–4000 planetesimals would be sufficient to reveal and study the underlying physics of planetesimal growth without overloading the integrations. We set the mass of each planetesimal to
$m = {10^{23}}$
g and its density to 2 g cm
${{\rm{\;}}^{\!\!- 3}}$
. The total mass of the disk was 0.05 Earth-masses.
In their simulations, Kokubo and Ida (Reference Kokubo and Ida1996) distributed planetesimals between 0.98 and 1.02 AU. These authors argued that, based on the results of Wetherill and Stewart (Reference Wetherill and Stewart1989), this width is large enough for the planetesimals to stay within its boundaries for the duration of the integrations. To examine the effect of the spatial distribution of planetesimals on their growth, we followed similar approach and considered two sets of simulations, one with an annulus from 0.96 to 1.04 AU (hereafter, set A) and one with an annulus extending from 0.98 to 1.02 AU (hereafter, set B). In each set, we ran five different simulations.
As mentioned earlier, runaway growth is triggered when simulations are carried out in 3D. We, therefore, considered a Gaussian distributionFootnote
7
for the eccentricities and inclinations of planetesimals with a zero mean and a dispersion of
$\lt{e}^2{\gt}^{1/2} = 2 {\lt}i^2\gt^{1/2} = 2h$
. Here
$h = {R_h}/a$
is the reduced Hill’s radius (Ida and Makino, Reference Ida and Makino1992a) and
$${R_H} = {\left( {{m}\over{{3{M_{\rm{\odot}}}}}} \right)^{1/3}}a$$
is the Hill radius of a planetesimal. In equation (2),
$m$
and
$a$
are the planetesimal’s mass and semimajor axis, respectively, and
${M_{\rm{\odot}}}$
is the mass of the Sun. The argument of the pericenter and the mean anomaly of each planetesimal were randomly chosen using a uniform distribution, and the longitudes of their ascending nodes were set to zero.
We integrated each system for
$5 \times {10^5}$
years using the hybrid routine in the
$N$
-body integration package MERCURY (Chambers, Reference Chambers1999). The time-steps of integrations were set to 6 days.Footnote
8
As stated above, integrations did not include the artificial inflation of planetesimals. We considered two objects with masses
${m_1}$
and
${m_2}$
on a possible collision course when they were closer than three times their mutual Hill’s radius
$\left( {{R_{{H_{12}}}}} \right)$
,
$${R_{{H_{12}}}} = {\left( {{{{m_1} + {m_2}}}\over{{3{M_{\rm{\odot}}}}}} \right)^{1/3}}\left( {{{{a_1} + {a_2}}}\over{2}} \right){\rm{\;}},$$
and, in cases when the objects collided, we simulated collisions as perfect-merging. In all integrations, the total energy and angular momentum were conserved with a relative error of
${10^{ - 10}} - {10^{ - 11}}$
.
Results and comparison with previous studies
Table 2 shows the masses and growth times of the three largest embryos at the end of each integration. As shown by the column “Final Mass,”Footnote 9 the masses of these bodies vary between 20 and 142 times the mass of the initial planetesimals. That means, not only did runaway growth occur in all our simulations, in most systems, the runaway bodies continued their growth even up to the boundary of the oligarchic regime.Footnote 10
Table 2. The final mass and the times of the growth of the three largest bodies (the black circles in Figures 1 and 2) in all simulations. The final mass is the mass of the body at the end of the 500,000 years of integration and is given in terms of the planetesimals’ initial mass (
${m_{{\rm{min}}}} = {10^{23}}$
g). Set A refers to initial distribution of planetesimals from 0.96 to 1.04 AU, and set B corresponds to an initial planetesimals distribution of 0.98 to 1.02 AU
Figures 1 and 2 show the snapshots of the evolution of a sample of our systems. In Figure 1, we show systems A1 and B4 as examples of those in which the mass of the largest embryo is more than 100 times the initial planetesimals, and Figure 2 shows systems A2 and B5 as samples of the rest of the systems. The black circles in each figure show the growth and orbital evolution of the three largest bodies at the end of the simulations. See the figures captions for more detail.
Figure 1. Snapshots of the evolution of systems A1 (left) and B4 (right) where the mass of the largest body is more than 100 times the initial planetesimals. Each object is represented by a red circle with its radius proportional to its mass. Blue circles represent bodies with masses at least 20 times their initial masses. The black circles in the bottom panel show the largest three bodies at the end of the integrations. Black circles in prior panels show the same bodies as they grow in time. Note the spreading of the disk and the dynamical diffusion of planetesimals (both red and blue circles) to outside the disk’s initial boundaries during its evolution.
Figure 2. Same as Figure 1, showing the evolution of systems A2 (left) and B5 (right) as samples of the rest of the systems. Note the spreading of the disk and the dynamical diffusion of the planetesimals (red and blue circles) to outside the disk’s initial boundaries during its evolution.
Runaway growth
As shown by Figures 1 and 2, runaway growth has occurred in these (and by the same token, in all our) simulations. We show the latter using two approaches. First, we use the definition of runaway growth as presented by Kokubo and Ida (Reference Kokubo and Ida1996). We will then explain the inconsistencies in using this definition and demonstrate the occurrence of runaway growth by using a more accurate criterion based on the analysis of the distribution of masses after a runaway growth system has reached relaxation (i.e., when its mass distribution follows a power-law; see equation 4 and Figures 4 and 5).
Figure 3. Graphs of the growth of the large bodies in the systems of Figures 1 and 2. Top panels show systems A1 (left) and B4 (right), and bottom panels show systems A2 (left) and B5 (right). The red, blue and green curves correspond to the three largest bodies (the black circles in Figures 1 and 2) and the curves in gray show the growth of bodies with a mass larger than 20 times the initial planetesimals
$\left( {{m_{{\rm{min}}}}} \right)$
. The black curve in each panel shows the variation of the mean mass of the system without including the largest body.
Figure 4. Graphs of the evolution of the mass distribution in systems A1 (left) and B4 (right). Each point represents the number of bodies with that mass. Because growth is through perfect-merging, each mass is a multiple of the mass of the initial planetesimals
$\left( {{m_{{\rm{min}}}}} \right)$
. The bottom panels show the fits to the final mass distribution and their corresponding slopes (see section “Mass distribution” and Table 3 for more details).
Kokubo and Ida (Reference Kokubo and Ida1996) defined runaway growth as when (1) the largest object (i.e., the runaway growth body) grows locally (i.e., within its accretion feeding zone) faster than the second largest body, and (2) that the ratio of the mass of the runaway growth body to the mean mass of the rest of the system increases with time. We show in Figure 3 the growth of all large bodies of Figures 1 and 2, as well as the time-variation of the mean mass of each system when the largest body is not included. As shown here, all large bodies grow much faster than the mean mass of their respective systems, fulfilling the above criteria for runaway growth. We have observed similar trend in all our simulations.
Figure 5. Graphs of the evolution of the mass distribution in systems A2 (left) and B5 (right). Each point represents the number of bodies with that mass. Because growth is through perfect-merging, each mass is a multiple of the mass of the initial planetesimals
$\left( {{m_{{\rm{min}}}}} \right)$
. The bottom panels show the fits to the final mass distribution and their corresponding slopes (see section “Mass distribution” and Table 3 for more details).
Although the above definition of runaway growth can explain the formation of our runaway bodies, it suffers from a fundamental inconsistency: The term “locally” in this definition refers to the local feeding zone of an object and requires this zone to be static. However, as an object grows, its position changes with time and so does its local feeding zone. This can be seen in Figures 1 and 2 where the semimajor axes, and therefore, the spatial locations of the growing bodies (e.g., the black circles) change during the integration. This change of the orbit and mass causes the location and size of the feeding zones of the bodies to change as well, making it practically impossible to determine the feeding zone and the second largest body to which the above definition of runaway growth applies. As noted by Kokubo and Ida (Reference Kokubo and Ida2000), and as explained below, a more accurate way of demonstrating whether runaway growth has occurred would be to study the mass distribution of the system when it has reached relaxation.
Mass distribution
To study the mass distribution, we note that, as shown by Figure 3, at the end of an integration, the mean mass of the rest of the system only doubles. In other words, most of the remaining mass is still in the form of the initial planetesimals (as we explain below, this is another characteristic of a system in which runaway growth has occurred). We demonstrate the latter in Figures 4 and 5 where we show the time evolution of the distribution of mass (that is, the number of objects for different values of their masses) in systems of Figures 1 and 2. Each point in these figures corresponds to the number of planetesimals with the same mass. Note that, because growth has been modeled through perfect-merging and because, initially, all planetesimals had equal masses, the growth of an object appears as the multiples of its initial mass.
As shown by Figures 4 and 5, there is a strong correlation between the number of bodies and their masses up to 10 times the initial mass of individual planetesimals (i.e., up to
${10^{24}}$
g.) The mass distribution in this interval seems to follow a power law of the form
$$N\left( m \right) = {N_1}{\rm{\;}}{m^{ - \alpha }},$$
where
${N_1}$
is the number of the planetesimals remaining from the initial population. As shown by Kokubo and Ida (Reference Kokubo and Ida1996, Reference Kokubo and Ida1998, Reference Kokubo and Ida2000, and references therein), in systems where runaway growth occurs,
$\alpha \gt 2$
. We show the latter in Table 3 where we have listed the values of
${N_1}$
and
$\alpha $
for all our systems at the end of their integrations. As shown here,
$\alpha $
stays larger than 2 for the entire time during integrations confirming that runaway growth occurred in all our simulations.
Table 3. Values of
${N_1}$
and
$\alpha $
for power-law fitting (equation 4) to the mass distributions at the end of all simulations. The values for
${N_1}$
were rounded to the nearest integer
Because equation (4) shows the distribution of small masses, if, while the objects grow (for instance, through perfect-merging), there is no re-supply of small bodies, the value of
$\alpha $
will decrease with time. Kokubo and Ida (Reference Kokubo and Ida1996) showed that in their simulations (where the planetesimals had been inflated),
$\alpha $
decreased from 2.6 to 2.4 while maintaining an averaged value of
$\sim 2.5 \pm 0.4$
. Kokubo and Ida (Reference Kokubo and Ida2000) reported that at the end of their simulations (where planetesimals had not been inflated), the final value of
$\alpha $
varied between 2.2 and 1.9. In our simulations,
$\alpha $
exhibited similar behavior. For example, in the systems of Figure 4,
$\alpha=2.57$
in the top-left panel and 2.6 in the top-right panel, and in the systems of Figure 5, these values are 2.8 and 2.6 for the top left and right panels, respectively. During the course of the integrations, these values dropped to those in Table 3, which are in stark agreement with the values reported by Kokubo and Ida (Reference Kokubo and Ida2000).
It is important to note that the mass distribution portrayed by Figures 4 and 5, and formulated by equation (4), while the direct consequence of the growth via perfect-merging, is independent of radius inflation. Because when a planetesimal is inflated, it is only its radius that is enhanced (and its mass is not changed), the planetesimal’s bulk density will be reduced by the same factor as its radius enhancement. Since the mass of the perfect-merging body is equal to the sum of the masses of the two colliding objects, the decrease in the bulk density of the impactors cancels out the effect of their size inflation, leaving their mass distribution intact.
However, as mentioned earlier, perfect-merging is unrealistic. In reality, collisions result in break-up, fragmentation, and shattering. Each of these events affects the number, as well as the mass and size distributions of the bodies, differently. The latter implies that, although the appearance of the runaway growth is independent of the modes of collision and growth (as its occurrence is the result of the appearance of non-uniformities in the planetesimals population), the time of its appearance, its duration, and the mass and size distributions of the resulted bodies will be different from those obtained from perfect-merging scenario. Unfortunately, at present, it is not possible to determine mass distribution realistically as computational facilities do not have the necessary capabilities to resolve collisions and growth among tens of millions of objects, accurately.
Remaining small bodies
An inspection of Figures 4 and 5 indicates that during the evolution of a system, a large number of small bodies (in the form of both initial planetesimals and objects a few times more massive) stay in the system. Figure 6 shows this by demonstrating the decrease in the number of objects in each system. We also show the number of the remaining bodies at each time of the integration on the upper left corner of the panels in Figures 1 and 2. As shown here, the number of small bodies at the earlier stage of integrations, for instance at
$t=125,000$
, where the value of
$\alpha $
is between 2.8 and 2.5, is even larger indicating that when runaway growth is dominant, most of the mass is still in the form of small objects (see the graph of the mean mass in Figure 3).
Figure 6. Graphs of the number of the bodies in each system of Figures 1 and 2 during the evolution of the system. Top panels show systems A1 (left) and B4 (right), and bottom panels show systems A2 (left) and B5 (right). Note that in systems A1 and B4 where the largest body is 100 times more massive than the initial planetesimals (as well as in simulation B5 where the largest body is 92 times more massive), the number of bodies drop by over 60% whereas in simulation A2 where the largest body is 47 times the initial planetesimals, the number of bodies drop by only 50%.
These remaining small bodies play a vital role in the subsequent evolution of the system. The collective effect of these objects increases the efficacy of collisional growth by damping the orbital eccentricities and inclinations of growing bodies through dynamical friction. We have demonstrate the latter in Figures 7 and 8 where we show the RMS values of the eccentricities (filled circles) and inclinations (open circles) of the bodies of the systems of Figures 1 and 2. As shown by the top panels, at the early stage of the evolution, when objects have just started to grow, their eccentricities and inclinations are strongly damped. This trend is maintained throughout the integration by small bodies ranging from 1 to 10 times the mass of the initial planetesimals (shown by
${m_{{\rm{min}}}}$
). As larger bodies appear, their perturbation causes the eccentricities and inclinations of all objects to increase. The latter can be seen for the value of the mass larger than
$10{\rm{\;}}{m_{{\rm{min}}}}$
and is more pronounced in Figure 7 which corresponds to systems A1 and B4. In these systems, embryos are larger than the systems of Figure 8 (A2 and B5, see Table 2) with the largest embryos being 1.5–3 times more massive, causing the eccentricities and inclinations of smaller bodies to raise to slightly higher values.
Figure 7. Graphs of the evolution of the RMS values of the eccentricities (filled circles) and inclinations (open circles) of planetesimals in systems A1 (left) and B4 (right). Planetesimals masses have been binned as in Figures 3 and 4. Note the damping of the eccentricity and inclination up to the point when the objects reach
$10{\rm{\;}}{m_{{\rm{min}}}}$
. As shown by Figures 4 and 5, this values of mass marks the onset of the runaway growth when the large bodies decouple from the rest of the planetesimals. At this stage, the perturbation of the larger bodies disturbs the dynamics of the planetesimals and their mutual interactions cause their eccentricities and inclinations to increase. As the objects grows, their perturbing effects becomes stronger to the extent that the dynamical friction due to the remaining population of small planetesimals can hardly damp their eccentricities and inclinations to lower values.
Figure 8. Graphs of the evolution of the RMS values of the eccentricities (filled circles) and inclinations (open circles) of planetesimals in systems A2 (left) and B5 (right). Planetesimals masses have been binned as in Figures 3 and 4. Note the damping of the eccentricity and inclination up to the point when the objects reach
$10{\rm{\;}}{m_{{\rm{min}}}}$
. As shown by Figures 4 and 5, this values of mass marks the onset of the runaway growth when the large bodies decouple from the rest of the planetesimals. At this stage, the perturbation of the larger bodies disturbs the dynamics of the planetesimals and their mutual interactions cause their eccentricities and inclinations to increase. As the objects grows, their perturbing effects becomes stronger to the extent that the dynamical friction due to the remaining population of small planetesimals can hardly damp their eccentricities and inclinations to lower values.
Similar results have been reported by other researchers, though with some differences. For instance, at the end of the simulations by Kokubo and Ida (Reference Kokubo and Ida1996), 65% of the initial planetesimals were accreted by the runaway growth bodies. Kokubo and Ida (Reference Kokubo and Ida2000) report this number at 56%. These authors carried out simulations for 200,000 years. In our simulations, the number of accreted bodies ranged between 48% and 57%. We believe that the larger number of accreted bodies in the work of Kokubo and Ida (Reference Kokubo and Ida1996) is due to the inflating of planetesimals at the beginning of simulation, and is, therefore, unrealistic. When realistic sizes are used, as in our simulations and those of Kokubo and Ida (Reference Kokubo and Ida2000), the collision cross-section is not artificially enhanced and as a result, the number of accreted bodies drops. The larger value of this number in the simulations of Kokubo and Ida (Reference Kokubo and Ida2000) can be attributed to the effect of nebular gas drag. As explained in the Appendix, the drag force of the nebula operates as an additional factor in damping the eccentricity which subsequently results in lowering the relative velocities of colliding bodies.
We would like to emphasize that because accretion is a function of time, it would, generally, be more meaningful to compare the rates of accretion in different simulations at the same times. However, in this study, that is unnecessary as our goal in the above paragraph is to merely demonstrate that the high rate of accretion in simulations with inflated planetesimals is in fact artificial and, therefore, such high accretion rates and any results obtained from them must not be considered as the natural evolution of the system.
An interesting result depicted by Figures 7 and 8 is the manner by which the small material at the end of the simulations reach their final eccentricities and inclinations. As shown here, in systems of Figure 7, the bottom panels show an almost constant distribution for the objects in the 1–10 mass range. However, in the systems of Figure 8, where the final embryos are smaller than those of Figure 7, objects in the same mass range continue their linear trend toward lower values of eccentricity and inclination. Note that the trend in the system on the right in Figure 8 (system A5) is weaker as in this system embryos are larger than those in system A2 shown by the left column. It is important to note that in these simulations, the dynamical evolution of the large bodies cannot be studied statistically (meaning, their evolution cannot be fitted or used to identify a trend) as the number of these objects is too small for the statistical analysis of their dynamics to produce meaningful results.
Implications of the results
In the last section, we presented the fundamental characteristics of a system of planetesimals within a narrow annulus in a protoplanetary disk. We showed what it means when the runaway growth occurs, and how the system evolves and responds to it. While most of the results can be found in previously published articles, we presented them here collectively, so that the reader could see the connection between them at one place. We also demonstrated how these results were obtained by carrying out relevant simulations. In this section, we follow the same approach and present the implications of the results. We discuss how specific initial conditions manifest themselves in the final outcome, and what the characteristics of the system suggest for carrying out realistic simulations (i.e., simulation with no ad hoc, unrealistic, and/or simplifying assumptions).
In general, the results presented in the previous section have four implications for the late stage of terrestrial planet formation. They show that (1) when planetesimals are not inflated, planetary embryos do not reach the high masses that simulations with inflated planetesimals suggest. That means, (2) in actual systems, the evolution of the protoplanetary disk and the process of planetesimal growth are subject to the perturbation of embryos smaller than those obtained from these simulations. They also imply that (3) the growth of embryos takes much longer to the extent that while they go through the runaway and oligarchic modes, their dynamical evolution and, therefore, their isolation masses and orbital architecture may be affected by the growing giant planets. Finally, results show that while simulating planetesimal growth using isolated distributions (e.g., a planetesimals annulus) provides an excellent approach to understanding the mechanics of the process and its underlying physics, (4) in order to develop realistic initial conditions that can be used in simulating the last stage of terrestrial planet formation, simulations of planetesimal growth have to be carried out for the entire of a planetesimal disk at the same time.
In the following we explain each of these implications in more detail.
The effect of inflating planetesimals
A comparison between the results of the last section and those in which planetesimals were initially inflated immediately demonstrates that our runaway growth bodies are, in general, much smaller than in those simulations. For instance, in the simulations of Kokubo and Ida (Reference Kokubo and Ida1996) and Richardson et al. (Reference Richardson, Quinn, Stadel and Lake2000), after only 20,000 years of integrations, the runaway growth bodies reached the masses of 300-400 time the mass of the initial planetesimals
$\left( {{m_{{\rm{min}}}}} \right)$
. In the works of Leinhardt and Richardson (Reference Leinhardt and Richardson2005) and Wallace and Quinn (Reference Wallace and Quinn2019), these values rose up to 1000–1500
${m_{{\rm{min}}}}$
during the same time. In our simulations, however, after 500,000 years of integrations, the masses of our runaway growth bodies reached only 20–142
${m_{{\rm{min}}}}$
with our largest embryos having a mass between 107
${m_{{\rm{min}}}}$
and 142
${m_{{\rm{min}}}}$
.
Similar results have also been reported by Kokubo and Ida (Reference Kokubo and Ida2000) and Clement et al. (Reference Clement, Kaib and Chambers2020). These authors did not initially inflate the planetesimals, and also included the effect of gas drag. In the simulations of Kokubo and Ida (Reference Kokubo and Ida2000), it took 200,000 years for the largest planetesimals to reach a mass of 200
${m_{{\rm{min}}}}$
, and in the work of Clement et al. (Reference Clement, Kaib and Chambers2020), the largest body in their annulus centered at 1 AU reached 550 times its initial mass after 100,000 years. It is important to note that the shorter time of growth in the work of Kokubo and Ida (Reference Kokubo and Ida2000) is due to the effect of gas drag, and that in the simulations of Clement et al. (Reference Clement, Kaib and Chambers2020), the time of growth is even shorter is due to the fact that these authors used almost twice as many planetesimals in their 1 AU centered annulus as those in ours and in the simulations of Kokubo and Ida (Reference Kokubo and Ida2000).
These small sizes of runaway growth bodies in simulations with uninflated planetesimals, combined with the value of the exponent
$\alpha $
in each of our systems (see Table 3), and the fact that in such simulations, integrations had to be carried out for 5–25 times longer than those with inflated planetesimals, strongly implies that in simulations where planetesimals are not inflated, the dynamical evolution of the system would follow a path in which it would naturally take much longer for planetary embryos to reach their isolation masses and for the growth to transition from runaway
$\left( {\alpha \ge 2} \right)$
to oligarchic mode
$(\alpha \lt 2)$
. For instance, an examination of Table 3 indicates that while systems A1, B4 and B5 are at the verge of transitioning to oligarchic growth, other systems still require more time, some of them (e.g., systems A2, A3, A5) much longer than our integrations time of 500,000 years.
The fact that in systems with inflated planetesimals, large bodies form in much shorter times indicates that from the early stages, the dynamical evolution of these systems, as well as the growth and orbital architecture of their bodies are heavily influenced by the perturbation of their unrealistically large embryos. We note that the large sizes of runaway growth bodies in simulations with inflated planetesimals is an expected result that has roots in the fact that increasing the initial sizes of planetesimals enlarges their collision cross sections and increases the rate of their growth causing very large bodies to form at short times. As a result, soon after the start of the simulations, these systems evolve along a path that is dictated to them by the perturbation of their unrealistically large bodies. This unreal evolution propagates to all stages of planetesimal growth and embryo formation making the spatial and mass distributions of these bodies unreal as well. The latter will affect the final stage of terrestrial planet formation causing the system to produce planetary masses and orbital architecture that are inherently contaminated by the unnatural evolution of their original systems.
We would like to note that, in general, a better comparison between the results of two simulations would be when the ratio of the number of their remaining bodies to the number of their initial bodies are similar. For instance, the number of remaining bodies at the end of the integrations of Kokubo and Ida (Reference Kokubo and Ida1996)
$(\sim 35\%)$
is about
$6.7{\rm{\% }}$
to
$17{\rm{\% }}$
smaller than the final bodies in our simulations implying that the results would have been better compared if our integrations had been continued till the number of the remaining bodies in our simulations would reach
$35{\rm{\% }}$
. However, for the purpose of this article, such extension of integrations would be unnecessary. As shown by Table 2, given the rate of the planetesimal growth in our system, it is very unlikely that the masses of our runaway growth bodies would have reached 300–400 initial planetesimals even if we had continued the integrations. Even if the effect of nebular gas drag had been included, given the similarity between our initial planetesimals population and those of Kokubo and Ida (Reference Kokubo and Ida2000) who considered gas drag, our planetesimals would have at most reached the maximum mass of 200
${m_{{\rm{min}}}}$
as obtained by these authors. That, however, would have had no significant qualitative impact on the implications of the results. In other words, while extending integrations to longer times would have allowed for a quantitatively more accurate comparison with previous works, the implications of the results, as stated in this section and summarized in section “Summary and conclusions,” stay intact.
The effect of the spatial distribution of planetesimals
As shown by Table 2, at the end of the integrations, with the exception of system A1, the masses of the three largest bodies in the systems of set B are much larger than those in set A. This is an expected result that is the consequence of the initial distribution of planetesimals. As mentioned in section “Results and comparison with previous studies,” in the simulations of set A, planetesimals were initially distributed over a range of 0.08 AU whereas in set B, same number of planetesimals were scattered over 0.04 AU. The latter concentration of planetesimals in a smaller region naturally increased the rate of their collisions, and because collisions were taken to be perfectly inelastic, it resulted in their growth to objects larger than in set A. This can also be seen in the work of Kokubo and Ida (Reference Kokubo and Ida2000) where the initial distribution of planetesimals was over a smaller region (0.02 AU) and as a result (and also because of the effect of gas drag), their largest body was even larger than those in set B (200 times the mass of the initial planetesimals).Footnote 11 In contrast, when the same number of planetesimals were distributed in a larger region (e.g., set A), the surface density of their annulus was smaller, which, because initial planetesimals were equal-mass, resulted in their number density to be smaller as well. A smaller number density resulted in a smaller number of collisions, and therefore, less massive bodies (see, e.g., the explanation about over-populated annuli as used in Clement et al., Reference Clement, Kaib and Chambers2020, in the second paragraph of section “The effect of inflating planetesimals”).
Expansion of initially localized distributions
It is important to note that during the evolution of a system, any local distribution of planetesimals will gradually expand. Figures 1 and 2 show that by the end of the integrations, the range of the planetesimals’ semimajor axes has expanded by about 0.04–0.06 AU. Similar expansions can also be found in the works of Kokubo and Ida (Reference Kokubo and Ida1996), Richardson et al. (Reference Richardson, Quinn, Stadel and Lake2000), Morishima et al. (Reference Morishima, Schmidt, Stadel and Moore2008), and Wallace and Quinn (Reference Wallace and Quinn2019). For instance, at the end of the simulations in Kokubo and Ida (Reference Kokubo and Ida1996), the semimajor axes of planetesimals expanded by 0.06 AU.
In general, this expansion is caused by (1) diffusion due to the mutual interaction among planetesimals (as in our simulations), and (2) in cases where the drag effect of the nebular gas is included, by gas drag which causes planetesimals to migrate toward the central star (known as gas drag-induced migration). When the initial distribution of planetesimals is localized, this expansion reduces planetesimals surface density and after a while, impedes further formation and growth of large bodies.
While diffusion and gas drag-induced migration are real physical processes that occur naturally during the dynamical evolution of planetesimals, the delay and/or the disruption of the growth is an artificial effect that appears because of the decrease in the surface density which itself is a consequence of using localized distributions. In a real system, where the initial extent of the planetesimals disk is over a much larger region, when a local distribution expands, its neighboring distributions also expand. As a result, while a local distribution is losing bodies due to migration and diffusion, planetesimals from its neighboring distributions migrate/diffuse into it, compensating for the ones that were lost. The latter prevents the decrease in the surface density of the disk and allows the growth to continue. It is important to note that because the mechanics of the expansion is different for disks with different surface densities, and because that is what determines how neighboring regions spill into each other, the efficiency of this mechanism is heavily governed by the disk’s initial surface density profile (see e.g., Izidoro et al., Reference Izidoro, Raymond, Morbidelli and Winter2015; Haghighipour and Winter Reference Haghighipour and Winter2016)
The fact that in a real disk, planetesimals from one local distribution enter their neighboring units prompted Kokubo and Ida (Reference Kokubo and Ida2000) and Clement et al. (Reference Clement, Kaib and Chambers2020) to account for the reduction of the surface density by introducing a boundary condition in which after a planetesimal has left one boundary of a local distribution, that planetesimal is removed from the integration and a new one with the same mass is added at the opposite boundary. In the simulations of the last section, we chose not to adopt this approach because while on the surface, this boundary condition remedies the loss of material and reduction of surface density, it suffers from the following shortcoming: The effect of an object artificially added to one boundary is not identical to the effect of a body that diffuses out of the system in the opposite side. Also, and more importantly, in the simulations presented here, our main intention has been to reveal the issues associated with using isolated distributions, and that reliable results cannot be obtained using this type of systems. Finally, as we explain in section “Simulating the entire disk at the same time”, unless integrations are expanded to the entire disk at the same time, any workaround, especially using a collection of isolated distributions will still suffer from above issues and internal inconsistencies.
Simulating the entire disk at the same time
That the decrease in the disk surface density, and the subsequent impeding of the growth are artifacts of using a localized distribution strongly implies that simulations with localized mass distributions, or a collection of localized mass distributions (i.e., rings), do not produce reliable results, and in a realistic simulation, the entire disk of planetesimals needs to be integrated at the same time.
At the time of this writing, the only attempt in extending integrations to the entire of a disk was due to Clement et al. (Reference Clement, Kaib and Chambers2020). As mentioned in the Introduction, these authors simulated planetesimal growth in five separate rings between 0.5 and 3 AU (see Table 1 for the locations of their planetesimal annuli). To extend their isolated simulations to a larger region (0.48–1.65 AU), the authors made the assumption that any ring-like distribution of planetesimals between each two of their main annuli, when integrated, would, in general, exhibit the same dynamical behavior as those in their five main rings, and produce similar results. With that assumption, the authors divided the region of 0.48–1.65 AU into small rings and interpolated the results of their simulations into these new areas.
It is important to note that the simulations of planetesimal growth are extremely stochastic. That means, they are unpredictable and cannot be regulated. This lack of predictability implies that in principle, it is not possible to ensure that the dynamical evolution of planetesimals in different distributions, even if geometrically similar, would be identical, or carry certain predictable features. The interpolation approach used by Clement et al. (Reference Clement, Kaib and Chambers2020), although undoubtedly clever and most likely useful in demonstrating some proofs of concepts, is subject to this unpredictability meaning that, its results, may not be the true representation of the dynamical evolution of planetesimals, if used in lieu of the direct integration of the entire system. In other words, to obtain a realistic image of the dynamical evolution of planetesimals, the entire disk has to be integrated at the same time.
A note on recent efforts on simulating planetesimal growth in ring-like distributions
Despite the above-mentioned effects of the expansion of isolated distributions on the growth and dynamics of planetesimal, in the past few years, a series of articles have promoted planetesimal growth in rings as a way of explaining the formation and orbital architecture of terrestrial planets in our solar system (Ogihara et al., Reference Ogihara, Kokubo, Suzuki and Morbidelli2018; Broz et al., 2022; Izidoro et al., Reference Izidoro, Dasgupta, Raymond, Deienno, Bitsch and Isella2022; Morbidelli et al., Reference Morbidelli, Baillié, Batygin, Charnoz, Guillot, Rubie and Kleine2022), as well as the formation of super-Earths in extrasolar planets (Batygin and Morbidelli, Reference Batygin and Morbidelli2023). In these studies, converging radial migration of solids in the gaseous component of the disk has been presented as the way of accumulating these bodies in ring-like distributions. This radial migration has been attributed to the appearance of pressure enhanced regions due to the interaction of stellar wind with MRI active/inactive regions of a disk (Ogihara et al., Reference Ogihara, Kokubo, Suzuki and Morbidelli2018), the mere assumption of a pressure enhancement near the CO snowline, water snowline, and silicate sublimation line (Izidoro et al., Reference Izidoro, Dasgupta, Raymond, Deienno, Bitsch and Isella2022; Morbidelli et al., Reference Morbidelli, Baillié, Batygin, Charnoz, Guillot, Rubie and Kleine2022; Batygin and Morbidelli, Reference Batygin and Morbidelli2023), and to the Lindblad, co-rotation, and heating torques applied to a protoplanet due to its interaction with the gas. These studies suggest that the accumulated bodies in such regions can grow through collision as well as accretion of pebbles, and form larger bodies including terrestrial planets and super-Earths. Recently Kambara and Kokubo (Reference Kambara and Kokubo2025) have shown that it is also possible for runaway growth bodies to undergo oligarchic growth while in such isolated distributions, and as their distributions expand.
We would like to note that while the underlying physics of converging gas drag-induced migration in the vicinity of pressure enhanced regions is solid and fully supported by the physics of solid-gas interaction (we refer the reader to Haghighipour and Boss, Reference Haghighipour and Boss2003a, b, for the full theory of dynamical evolution of solid objects in the vicinity of gas pressure/density enhancements), and although the appearance of Lindblad and co-rotation torques can indeed affect the orbits of protoplanetary bodies, caution must be taken in using results of these simulations as these models include assumptions that have been tailored to facilitate the goals of their studies. For instance, in the simulations of Ogihara et al. (Reference Ogihara, Kokubo, Suzuki and Morbidelli2018), the accumulation of planetesimals at the pressure enhanced regions and their subsequent growth occur only if the planetesimals are km-size and no larger than 10 km. The process fails when the objects are larger. In contrast, in the work of Broz et al. (2022), torques are more efficient in accumulating bodies when the objects are large, as the interaction of these bodies with the gaseous disk needs to be strong so that the resonance density waves and their resulting spiral arms can produce strong Lindblad and co-rotation torques. Finally, in the works of Izidoro et al. (Reference Izidoro, Dasgupta, Raymond, Deienno, Bitsch and Isella2022), Morbidelli et al. (Reference Morbidelli, Baillié, Batygin, Charnoz, Guillot, Rubie and Kleine2022) and Batygin and Morbidelli (Reference Batygin and Morbidelli2023), the pressure enhancements do not naturally appear, but they have merely been assumed, and their associated gas drag-induced migration has been included in the equations of motion, analytically.
The effect of the long time of the runaway growth
As shown by Figures 1, 2 and Table 2, for any given mass, the time of growth in our simulations is considerably long. For instance, after almost
$\left( {4.5 - 5} \right) \times {10^5}$
years, the five largest embryos in all our simulations grew to only 92–142 times the initial planetesimals
$\left( {{m_{{\rm{min}}}}} \right)$
, and the three largest bodies in each system grew only slightly larger than
$100{\rm{\;}}{m_{{\rm{min}}}}$
. This long formation time is a characteristic of realistic simulations and is the main reason that many researchers opted for inflating planetesimals. Kokubo and Ida (Reference Kokubo and Ida1996) showed that when the gravitational focusing is effective, an increase in the initial sizes of planetesimals by a factor
$f$
reduces the time of the formation of the runaway growth bodies by a factor between
$f$
and
${f^2}$
. These authors inflated planetesimals by a factor
$f=5$
and their runaway growth bodies reached the mass of 400 times the initial planetesimals (almost 4 times the masses of our final bodies) in ~20,000 years, a time that is
${5^2}=25$
times shorter than the growth time of the five largest embryos in our simulations.
Similar long formation times have also been reported by Kokubo and Ida (Reference Kokubo and Ida2000). Recall that these authors did not inflate the initial planetesimals. The results of the simulations by these authors showed a runaway growth body with a mass 200 times the initial planetesimals forming within
$2 \times {10^5}$
years. Although this time is shorter than the formation time in our simulations, it is still 10 times longer than the time of the formation of runaway growth bodies in simulations of Kokubo and Ida (Reference Kokubo and Ida1996) and Richardson et al. (Reference Richardson, Quinn, Stadel and Lake2000). Other simulations such as those by Kortenkamp et al. (Reference Kortenkamp, Wetherill and Inaba2001), Carter et al. (Reference Carter, Leinhardt, Elliott, Walter and Stewart2015) and Wallace and Quinn (Reference Wallace and Quinn2019) also hinted at long formation times. For instance, Carter et al. (Reference Carter, Leinhardt, Elliott, Walter and Stewart2015) and Wallace and Quinn (Reference Wallace and Quinn2019) simulated planetesimal growth using inflated planetesimals and stated that in order to produce similar results using uninflated planetesimals, they would have had to continue integrations to at least 750,000 years or more likely to a few million years.
Given that runaway growth is an inevitable consequence of the stochastic nature of planetesimals collisions, the timescales presented in our simulations, supported by their agreements with theoretical predictions and results of similar studies, strongly suggest that, in real disks, the time of the formation of runaway growth bodies is most likely of the order of a few to several hundred thousand years and maybe even longer. These long timescales, fall within the time of the growth of the cores of giant planets and the process of their gas-accretion. As shown by Haghighipour and Scott (Reference Haghighipour and Scott2012), in our solar system, the growth of giant planets strongly affects the dynamics of small bodies interior to the orbit of Jupiter. It is, therefore, imperative that the simulations of planetesimal growth in the inner region of the Solar System protoplanetary disk to be carried out simultaneously and in concert with the simulations of the growth of giant planets in the outer regions. Carter et al. (Reference Carter, Leinhardt, Elliott, Walter and Stewart2015) have presented such calculations in the context of the Grand Tack model. However, their study suffers from two shortcomings. First, these authors used inflated planetesimals, and second, as shown by Morbidelli et al. (Reference Morbidelli, Nesvorny, Laurenz, Marchi, Rubie, Elkins-Tanton, Wieczorek and Jacobson2018), Grand Tack model is inconsistent with the compositional properties of the moon and terrestrial planets, implying that the problem is still open and simulations needs to be carried out with giant planets forming in or in the vicinity of their current orbits.
Summary and conclusions
In this article, we presented the requirements that simulations of planetesimal growth need to fulfill in order to portray a more realistic image of the early stages of planet formation, and to produce physically meaningful and reliable results. To place the analyses in the right context, and to maintain focus on the basics, we started by an in-depth review of the history of the field where we discussed different approaches to modeling the growth of planetesimals, as well as their advantages and shortcomings. We focused our analysis on direct,
$N$
-body integrations as these integrations present the most precise approach to planetesimal growth. We identified the artificially increasing planetesimals’ sizes at the start of simulations, a practice that has been used by many researcher for the past three decades for the mere purpose of reducing the integration times (see Table 1), as the most unrealistic aspect of these integrations. We showed how inflating planetesimals affects the results, and how its unrealistic outcome extends to the subsequent phases of oligarchic growth and giant impacts.
With the goal of identifying the requirements for a realistic approach, we discussed the results of the simulations of planetesimal growth by carrying out such simulations in a localized environment without inflating their initial sizes. We identified the above requirements and demonstrated how they would emerge from the natural evolution of the system by analyzing the implications of the results in connection with those of similar simulations in the literature. To maintain focus on the underlying physics of runaway growth and the dynamical characteristics of the system, we considered collisions to be perfectly inelastic and did not include additional physical processes such as gas drag, fragmentation, and the effect of debris. It is understood that neglecting some of these processes, in particular, gas drag, may not allow the results to be fully realistic. However, comparisons with previous works indicate that the lack of these processes has only moderate quantitative effects on the final results, and do not change the nature of their implied constraints and requirements, qualitatively.
The following presents a summary of the analyses and their implied requirements for any realistic model.
Requirement 1 – Using uninflated planetesimals: For similar integration times, the masses of the runaway growth bodies in simulations with uninflated planetesimals are considerably smaller than those with inflated planetesimals. For instance, as demonstrated in sections “Results and comparison with previous studies” and “Implications of the results,” the largest body reaches a mass of 142 times that of the initial planetesimals whereas in simulations where planetesimals were initially inflated, this value varies between 300 and 1500. In the latter simulations, the perturbation of these highly exaggerated bodies affect the dynamics of their surrounding in an unrealistic way. Because when using inflated planetesimals, growth happens in much shorter time, these unreal perturbations affect the evolution of the disk from the beginning, causing the formation, growth, and orbital evolution of other bodies to be unrealistically perturbed as well. Using real-size planetesimals prevents all this by allowing the system to evolve naturally and not be forced through an unrealistic dynamical path.
Requirement 2 – Including the effect of growing giant planets: In systems with uninflated planetesimals, the formation of large planetary embryos and the onset of the oligarchic growth takes several hundred thousands of years. For instance, in the simulations presented here, after 500,000 years of integrations, the runaway growth was still in progress and the largest body reached a mass of only 142 times the mass of the initial planetesimals. In contrast, in simulations with inflated planetesimals (e.g., Kokubo and Ida, Reference Kokubo and Ida1996; Wallace and Quinn, Reference Wallace and Quinn2019), bodies 3.5–120 times larger than ours formed during an integration time that was 25 times shorter. That means, to produce systems that would show transition into the oligarchic growth, we had to continue the integrations much longer than 500,000 years. Similar long times for the growth of embryos were also reported by Kokubo and Ida (Reference Kokubo and Ida2002) where, using an analytical model, the authors show that the growth timescale of embryos at 1 AU is of the order of 1 Myr and it may reach 10 Myr at 5 AU for different radial profiles of the disk surface density function. As shown by Haghighipour and Scott (Reference Haghighipour and Scott2012), such integration times are comparable to the time that a growing giant planet in the orbit of Jupiter and Saturn manifests itself by perturbing the dynamical architecture of the asteroid belt. The latter strongly implies that in realistic simulations of planetesimal growth, in addition to using uninflated planetesimals, integrations need to be carried out for long times and include the perturbation of growing giant planets as well.
Requirement 3 – Simulating the entire disk at the same time: During an integration, the spatial distribution of planetesimals expands. In the systems of set A, where the initial distribution had a localized width of 0.08 AU, expansion was as large as 0.04 AU. In set B, initial distribution was more concentrated and had a width of 0.04 AU. In these systems, the expansion was about 0.06 AU. This spatial expansion reduces the surface density of planetesimals which causes their growth to slow down or even be hindered. In a real disk, where the distribution of planetesimals expands over a much larger area, the drop in the local number/surface density of planetesimals is complemented by the expansion of adjacent local distributions and, as a result, the overall surface density of the disk remains unchanged. This means, while simulating planetesimal growth in a localized distribution can be helpful in unraveling the underlying physics of the processes, and while some remedies such as the interpolation scheme presented by Clement et al. (Reference Clement, Kaib and Chambers2020) can help with proofs of concepts, in real systems, simulations need to be carried out for the entire disk at the same time.
Requirement 4 – Including dynamical friction due to the population of small bodies: In general, the eccentricities and inclinations of all large bodies stay low due to the dynamical friction with the remaining small planetesimals. In our simulations, the largest bodies maintained an eccentricity smaller than 0.01 (as a point of comparison, in the simulations of Kokubo and Ida, Reference Kokubo and Ida1996, Reference Kokubo and Ida2000, this values ranged between 0.001 and 0.005). This is an important result that has significant implications for the formation of the final planets and their orbital architecture. It indicates that, although during an integration, the orbits of individual planetesimals are affected by the perturbation of the growing bodies, and although their numbers may drop by 40–60%, still at any time during an integration, their population remains large enough for the dynamical friction to manifest itself and be effective. In other words, small bodies such as small planetesimals (and in the modern simulations of planet formation, debris and small fragments created by collisions, see e.g., Crespi et al., Reference Crespi, Dobbs-Dixon, Georgakarakos, Haghighipour, Maindl, Schäfer and Winter2021) are essential parts of planet formation without which simulations will not produce reliable results.
Requirement 5 – Including nebular gas drag: Although, in order to maintain focus on the underlying physics of runaway growth, we did not consider the drag force of the nebula, we would like to note that any realistic simulation of planetesimal growth needs to include the effect of the gas drag as well. Some authors have stated that the inclusion of gas drag would only moderately change the results (e.g., Wallace and Quinn, Reference Wallace and Quinn2019). However, simulations presented here demonstrate that when this effect is not included, results will be noticeably different. For instance, the three largest bodies in the simulations of Section “Results and comparison with previous studies” reached the masses of 119–142 times the initial planetesimals in 500,000 years whereas in the simulations of Kokubo and Ida (Reference Kokubo and Ida2000), where the authors used real-size planetesimals and gas drag was included, the largest body reached 200 times the initial planetesimals in less than half the time of the integrations presented here (
$\sim 200,000$
years). This shorter time of the appearance of large bodies strongly implies that the subsequent dynamical evolution of the protoplanetary disk in systems with gas drag will also be considerably different from those that do not include this effect. We refer the reader to the Appendix for more details on the significance of this effect.
In conclusion, our analyses demonstrate that accurate simulations of terrestrial planet formation need to start from a disk of planetesimals extended over a large area (e.g., from 0.5 au to at least 4.5 au), use uninflated bodies (i.e., bodies with their actual sizes), and a mass distributions that includes objects of different sizes and masses. These simulations need to be carried out simultaneously with the growth of giant planets so that the gravitational perturbation of these bodies are taken into account. They also need to resolve collisions accurately, by using, for instance, SPH (Smoothed Particle Hydrodynamics) simulations as in the works of Burger et al. (Reference Burger, Bazsó and Schäfer2020) and Reinhardt et al. (Reference Reinhardt, Meier, Stadel, Otegi and Helled2022), include dynamical friction due to the background small bodies and the impact debris, and include the nebular gas drag as these two effects will contribute to the growth of planetesimals by damping their eccentricities and, therefore, their impact velocities.
Acknowledgments
We would like to express our deepest gratitude to Dr. Patryk Sofia Lykawka for his in-depth review of our manuscript and excellent recommendations that have greatly improved our paper. We are also thankful to the La Plata Institute for Astrophysics (IALP) for extensive use of their computational facilities and to the Information Technology division of the Institute for Astronomy at the University of Hawaii-Manoa for maintaining computational resources that were used in this project.
Financial support
LAD acknowledges financial support from the Argentine National Institute for Advancement in Science and Technology (ANPCyT) through grants PICT 2016-2635 and PICT 201-0505, the National University at La Plata, Argentina, (UNLP) through grant PID G144, and from The Department of Astronomical & Geophysical Science at the UNLP (FCAGLPUNLP). NH acknowledges support from NASA XRP through grant numbers 80NSSC18K0519, 80NSSC21K1050, and 80NSSC23K0270 and NSF grant AST-2109285.
Competing interests
None.
Appendix. The effect of nebular gas
Because the growth of planetesimals begins at the early stage of the evolution of the protoplanetary disk, due to their small sizes, the dynamics, and consequently, growth of these objects are strongly affected by their interaction with the nebula through gas drag. For a spherical body with a radius
${R_{\rm{p}}}$
, the resistive force of the gas can be formulated as (Landau and Lifshitz, Reference Landau and Lifshitz1959; Adachi et al., Reference Adachi, Hayashi and Nakazawa1976; Haghighipour and Boss, Reference Haghighipour and Boss2003a)
$${F_{{\rm{drag}}}} = - {\rm{\;}}{{1}\over{2}}{C_{\rm{D}}}\pi R_{\rm{p}}^2{\rm{\;}}{\rho _{\rm{g}}}\left( {{r_{\rm{p}}}} \right){\rm{\;}}v_{{\rm{rel}}}^2{\rm{\;}}.$$
In this equation,
${r_{\rm{p}}}$
is the radial distance of the planetesimal to the central star,
${\rho _{\rm{g}}}\left( {{r_{\rm{p}}}} \right)$
is the density of the gas at the location of the planetesimal, and
${v_{{\rm{rel}}}}$
represents the velocity of the planetesimal relative to the local gas velocity. The quantity
${C_{\rm{D}}}$
in equation (A1) is the drag coefficient whose value depends on the size of the object and properties of the gas. When applied to the simulations of planetesimal growth, it is customary to take
${C_{\rm{D}}}=2$
. For more details on this quantity and its values, we refer the reader to section 3.1 in Adachi et al. (Reference Adachi, Hayashi and Nakazawa1976) and section 2.1 in Haghighipour and Boss (Reference Haghighipour and Boss2003a).
Studies of the effect of gas drag on the growth of planetesimals can be found in the fundamental works of Adachi et al. (Reference Adachi, Hayashi and Nakazawa1976), Ohtsuki et al. (Reference Ohtsuki, Ida, Nakagawa, Nakazawa, Levy, Lunine and Tucson1993), Kokubo and Ida (Reference Kokubo and Ida2000), Thebault et al. (Reference Thébault, Marzari, Scholl, Turrini and Barbieri2004) and Haghighipour (Reference Haghighipour2005). Collectively, these studies have demonstrated that gas drag has a positive effect on the rate of planetesimal growth as it reduces their impact velocities by damping their orbital eccentricities and inclinations. For instance, as shown by Thebault et al. (Reference Thébault, Marzari, Scholl, Turrini and Barbieri2004), gas drag aligns the periastrons of the orbits of planetesimals reducing their mutual impact velocities, and thereby enhancing their collisional growth especially among objects of the same size. The above works have also shown that the damping effect of gas drag is smaller for larger objects meaning that at the early stages of planetesimal growth, this effect has been more prominent.
The evidence to the fact that gas drag enhances collisional growth can be found in the results of the study by Kokubo and Ida (Reference Kokubo and Ida2000). These authors presented the first numerical simulations in which uninflated planetesimals grow through perfect-merging while subject to the drag force of the nebula. Using the same initial set-up as those presented here (section “Revisiting planetesimal growth”), these authors found that, after 200,000 years of integration, their largest runaway growth body reached the mass of 200 times the initial planetesimals and after 500,000 years, this value rose up to several hundred. By comparison, in simulations of section “Results and comparison with previous studies,” the largest object is 142 times more massive than initial planetesimals and it took 470,000 years to form. Also, as shown by these authors, after 200,000 years of integration, the RMS values of the eccentricities and inclinations of planetesimals show a relatively flat distribution. In systems studied here, however, while after 250,000 years of integration, the RMS values of these quantities show, in general, a similar trend, their distribution is not flat and contains scattered values as well.
The comparison between the results presented here and those of Kokubo and Ida (Reference Kokubo and Ida2000) shows that while, as expected, results of the simulations with and without gas drag are qualitatively similar, some subtle quantitative difference exist. These differences strongly imply that in order to be able to build a realistic and quantitatively accurate model of terrestrial planet formation, the simulations of planetesimal growth need to include the effect of gas drag as well.
Open access
