Habitable Planet Formation around Low-mass Stars: Rapid Accretion, Rapid Debris Removal, and the Essential Contribution of External Giants

For the purposes of our current investigation, we are chiefly interested in the total masses and orbital distributions of planetesimals and collisional fragments that would be generated in (or stranded by) the ultimate epoch of giant impacts that complete the formation of our systems. The second part of our analysis involves studying the dynamical evolution of this debris in three multi-planet, M-dwarf-hosted systems (TRAPPIST-1, Proxima Centauri, and TOI-700: Gillon et al. 2016; Anglada-Escudé et al. 2016; Gilbert et al. 2020). In all cases, we utilize the current published orbits and masses for the systems' planets as inputs for our numerical simulations (Anglada-Escudé et al. 2016; Gillon et al. 2017; Lienhard et al. 2020; Damasso et al. 2020; Kervella et al. 2020; Suárez Mascareño et al. 2020; Agol et al. 2021). TRAPPIST-1 is an 0.0898 M_⊙ M8V dwarf (Lienhard et al. 2020) with three planets in the HZ (e, f, and g). Proxima Centauri is an 0.12 M_⊙ M5.5V dwarf (Kervella et al. 2017) with the ∼1.6 M_⊕ planet Proxima Centauri b situated in the HZ. TOI-700 is an 0.41 M_⊙ M2V dwarf (Gilbert et al. 2020) hosting three planets, the last of which (d) is potentially habitable and has a mass of ∼1.72 M_⊕ (Rodriguez et al. 2020). In the case of Proxima Centauri, we elect to include the suspected small, ∼0.3 M_⊕ interior planet reported by Suárez Mascareño et al. (2020). We chose these systems because they encompass a range of stellar masses and each possess an approximately Earth-sized, potentially rocky planet in the HZ.

We derive leftover planetesimal populations directly from our simulations sets that incorporate planetesimals (see Table 1 and the left panel of Figure 1), and perform an additional batch of ∼10⁵ simulations of the final few giant impacts in our systems to derive collisional fragment distributions. This additional suite of computations utilizes a modified version of the Mercury integration package described in Chambers (2013) that incorporates collisional parameter space mapped in Stewart & Leinhardt (2012) to track hit-and-run and fragmenting impacts. Aside from the fragmentation algorithm, the integration parameters (i.e., time step, ejection radius, etc.) are identical to those utilized in Section 2.1. To prevent the calculation from becoming intractable, the user must establish a minimum fragment mass (MFM; 0.005 M_⊕ in our simulations based on previous experience with the code; Quintana et al. 2016; Clement et al. 2019b). When the algorithm detects a fragmenting collision, the total remnant mass is divided into a number of equal-mass fragments with M > M_MFM that are ejected in uniformly spaced directions in the collisional plane at v ≃ 1.05v_esc. If the total eroded mass is less than the MFM, the collision is considered totally accretionary. While this treatment of fragmentation is obviously somewhat contrived (see Section 4 for a discussion of some of these caveats), a first-order dissection of the distinctive populations of collisionally generated debris and unaltered planetesimals is important for our study as the fragment material might posses enhanced volatile and light element inventories if it is ejected from near the surface of partially or fully differentiated proto-planets. Moreover, since we utilize these simulations to infer the orbital distribution of debris given a supposed size frequency distribution (SFD), our particular selection of MFM is not carried forward into our study of late bombardment, and is therefore only important for keeping our experiments computationally feasible.

Zoom In Zoom Out Reset image size

To generate initial conditions for these simulations, we extract the state of our planet-formation models at the point where N_emb = 10 (15 for simulations that do not include planetesimals), remove all small bodies, and apply small changes (δ_a ≲ 0.001) to the embryos' orbits to generate unique evolutions. These systems are then integrated for 20 Myr. Through this process, we generate a population of 1285 surviving collisional fragments (right panel of Figure 1).

While the distributions of fragment orbital elements are relatively similar in our various simulation sets investigating different stellar and disk masses, perturbations from an external, Neptune-mass planet mildly heat the eccentricities and inclinations of the small remnant bodies compared to the simulations that do not incorporate a giant planet model. We also find no differences between the distribution of fragments produced in systems that originated in 200 embryos compared with those beginning with 20 embryos (Table 1). However, it is clear from Figure 1 that surviving fragments in the inner component of the disk are much rarer than surviving planetesimals. We find this to be the result of high rates of accretion by the central body for fragments generated at low heliocentric distances (by virtue of the low perihelia distances they are endowed with when ejected by our code; see also Clement et al. 2019a). As these fragments would likely be removed rapidly by Yarkovsky drift (see Section 4), we do not view this as a serious caveat on our results. However, it should be noted that the exact physics underlying the process of collisional fragmentation during planetary-scale collisions is still not fully understood. For these reasons, we elect to treat our fragment and planetesimal populations separately in the majority of our bombardment experiments.

The largest concentration of surviving debris material in all of our simulations lies in the vicinity of the most distant rocky planets, including a sizable number of planetesimals and fragments that are scattered into the "exo-asteroid belt" region exterior to our outer terrestrial disk boundary of 0.5 au (e.g., Izidoro et al. 2016; Raymond & Izidoro 2017). On average, these belt-analogs possess total masses of 0.034 M_⊕ in scattered planetesimals and 0.012 M_⊕ in collisional fragments. Simulations incorporating a Neptune-mass planet at 1.0 au have systematically less massive belts (an average of just ∼0.0064 M_⊕ in scattered planetesimals) as a result of scattering events with the massive planet, direct accretion of planetesimals, and the presence of overlapping first-order mean motion resonances near 1.0 au (e.g., Wisdom 1980; Deck et al. 2013). Conversely, integrations considering more massive central stars strand a larger fraction of material in the terrestrial and asteroid belt regions (an average of 0.076 M_⊕ and 0.11 M_⊕ worth of leftover planetesimals in our 0.5 M_⊙ and 0.6 M_⊙ simulations, respectively). We attribute this result to the relatively larger orbital velocities in these simulations providing stronger gravitational scattering events that disfavor accretion.

The average total leftover fragment (∼0.04 M_⊕) and planetesimal (∼0.06 M_⊕) mass in our simulations (combining debris in the vicinity of the terrestrial planets with exo-asteroid belt material) is remarkably similar to that in similar simulations of terrestrial planet formation in the solar system (presumably several lunar masses to provide a sufficient population of late impactors to source the late veneer, e.g., Raymond et al. 2013; Brasser et al. 2016; Clement et al. 2019b). Systems considering an external, Neptune-mass body tend to possess lower total leftover masses in the terrestrial and exo-belt regions (for additional data, see Table 3). This is likely attributable to resonant perturbations from the planet exciting the velocity dispersion of the small-body population, thereby providing increased high-speed collisions that disfavor perfect accretion (e.g., Stewart & Leinhardt 2012). However, the range of outcomes in terms of total debris mass overlap significantly in each of our simulation batches.

We experiment with two different approaches for deriving late bombardment curves in our chosen systems of M-dwarf hosted planets. In both cases, we adjust the semimajor axes of all our debris particles by linearly scaling the orbital periods of our initial inner and outer disk boundaries (0.01 and 0.5 au) to the inner and outermost planet in each system. In our first experiment, we simply integrate each debris particle ⁵ as a semi-active body in each system of interest for 1 Gyr (in practicality, we perform 52 separate simulations, each incorporating ∼50 different debris particles). These simulations leverage the same Mercury integrator and settings as our planet-formation simulations described in Section 2.1. Through this process, we generate chronological impact curves for each planet that we can calibrate by varying both the presumed initial SFD and each impactor's assumed initial composition (which we relate to its initial origin in our planet-formation simulations).

To validate our above methodology, we perform a second set of similar bombardment experiments (studying our collisional fragment distributions) that incorporate much higher particle resolution by leveraging the GENGA GPU integration package (a GPU-parallelized version of the Mercury code described in Grimm & Stadel 2014). Aside from the change in integrator, the set-up for these simulations is identical to that of the models described above. To take advantage of GENGA's superior performance, we include 10,000 debris particles in these simulations. The orbits of these new objects are interpolated by sampling from the original debris populations' distributions of semimajor axes, eccentricities and inclinations. The masses of the particles are assigned in a manner to mimic the SFD of the modern asteroid belt down to 5 km objects.

Dynamically speaking, our numerically generated systems are qualitatively similar to known multi-planet transiting systems, and the solar system's terrestrial worlds. We quantify the degree of dynamical excitation in our systems by calculating the angular momentum deficit (AMD) of all planets with M > 0.05 M_⊕ (Laskar 1997):

$$\begin{eqnarray}&&\mathrm{AMD}=\displaystyle \frac{{\sum }_{i}{m}_{i}\sqrt{{a}_{i}}[1-\sqrt{(1-{e}_{i}^{2})}\cos {i}_{i}]}{{\sum }_{i}{m}_{i}\sqrt{{a}_{i}}}\end{eqnarray} \tag{ 3 }$$

The AMD of a system essentially gauges the degree to which the collection of orbits deviate from that of a perfectly circular, co-planar system. Figure 3 plots the cumulative distribution of AMDs in our various simulation suites testing different values of M_Disk and giant planet models. While the inclusion of an external massive planet (on a circular orbit) has a negligible effect on the resulting planetary system's dynamical excitation, the difference in AMDs between our simulation sets testing M_Disk = 3.0 and 6.0 M_⊕ demonstrates the trade-offs of forming HZ planets in situ from a massive disk. While a high initial surface density of solids is essential for growing planets of the mass of the Earth, such concentrated regions of proto-planets naturally produce systems of fully evolved planets on more eccentric and inclined orbits. Thus, the HZ planets themselves would likely undergo larger magnitude secular oscillations (most consequentially in obliquity and perihelia, e.g., Berger 1978; Laskar et al. 1993) that change their climates over kiloyear to megayear timescales and would also be more susceptible late planetary collisions and ejections (e.g., Laskar & Gastineau 2009; Batygin et al. 2015). While these phenomena might present challenges for the emergence of life, it is important to note that models of terrestrial planet formation in the solar system incorporating numerical methodologies akin to that presented here systematically struggle to replicate the low eccentricities of Earth and Venus (e.g., Chambers & Wetherill 1998; Raymond et al. 2006; Clement et al. 2021). Moreover, it is still unclear whether short-period rocky planets around M-dwarfs universally possess low eccentricities and inclinations (as expected for a resonant system like TRAPPIST-1: Grimm et al. 2018), or more moderate values capable of driving large obliquity and perihelia cycles (e.g., Shan & Li 2018; Quarles et al. 2020).

Zoom In Zoom Out Reset image size

Figure 4 plots the cumulative distribution of neighboring-planet period ratios in our various simulation sets against those of all M-dwarf-hosted multi-planet systems. Here, we only consider exoplanets with a < 0.5 au and R < 2.0 R_⊕ to more accurately compare with our modeled scenario. It is clear that our simulations considering disk masses of 6.0 M_⊕ provide a reasonable match to the observed period ratios of exoplanets, however, it is important to note that these observations are biased due to potentially undetected planets. In particular, undetected small planets (e.g., Faridani et al. 2021) similar to those formed in our M_Disk = 3.0 M_⊕ simulations would likely shift the observed curve left by virtue of potentially inhabiting more compact orbital architectures. Nevertheless, Figure 4 depicts a reasonable match to the observed pileup of observed systems with orbital period ratios near, but not exactly within, the 3:2 MMR. This perceived under-abundance of resonant planets compared to predictions from convergent migration models (namely applied to Super-Earths orbiting Sun-like stars, e.g., Ogihara et al. 2015, 2018) has also been invoked to suggest that such planets acquired their orbital architectures via delayed episodes of dynamical instability (Izidoro et al. 2017, 2021) after forming in quasi-stable resonant configurations (e.g., Volk & Gladman 2015). While such a formation model is certainly plausible for explaining chains of resonant and non-resonant short-period planets around M-dwarfs like TRAPPIST-1 and Kepler-186, our simulations indicate that in situ formation is also capable of reconciling the period ratios of observed planets.

Zoom In Zoom Out Reset image size

Our simulated terrestrial disks experience moderate to strong radial mixing of material that systematically boosts the water-ice contents of HZ planets. Specifically, ≳83° of HZ planets in each of our simulation sets where the snow line is located interior to our outer disk boundary accrete at least one object from beyond the line. To quantify this enhancement, we crudely assume that all planetesimals and embryos originating inside of 0.8a_snow possess initial WMFs of 10⁻⁵, those between 0.8 and 1.0 a_snow are endowed with 0.1% of their mass in the form of water, and exterior disk objects begin with 10% of their mass in the form of water-ice. We discuss the limitations of these assumptions in Section 2.1. Figure 5 plots the cumulative distribution of final WMF for all HZ planets formed in our our M_* = 0.1, 0.2, 0.3, and 0.4 M_⊙ simulation batches (note that the snow line is beyond the edge of our initial terrestrial disk in our investigations of more massive M-dwarfs). The boosted WMFs in simulations investigating smaller stars thus demonstrate the ability of HZ planets to accrete material from the more distant regions of their initial planet-forming disks (i.e., WMFs are not boosted simply as a result of accretion from the edge of the snow line). Indeed, we note multiple incidences of HZ planets in our 0.1 M_⊙ batch with 0.025 < a < 0.057 au accreting planetesimals initialized at semimajor axes beyond 0.3 au. Overall, the radial mixing in our simulations is similar to that noted in classic models of terrestrial planet formation in the solar system employing similar disk setups to our models (e.g., Raymond et al. 2004, 2007a; O'Brien et al. 2018).

Zoom In Zoom Out Reset image size

The large pileup of Earth analogs with water contents of ∼1% corresponds to objects that accrete exactly one embryo from beyond the snow line. Assuming a nominal ∼1.0 M_⊕ planet in our batch of 6.0 M_⊕ disks (or, equivalently, a ∼0.5 M_⊕ planet in our set of 3.0 M_⊕ disks), an impact from a water-ice-rich embryo would deliver approximately 14% of the Earth analogs' ultimate mass. If the impactor possessed a water-ice content of 10%, the resulting HZ planet's WMF would be approximately 10⁻². Thus, the final water contents themselves are essentially an artifact of our initial disk's presumed compositional structure, and the more important result is the observed strong radial mixing.

While our numerical estimates of each HZ planets' final WMF are obviously contrived, it is clear that Earth analogs formed in situ around M-dwarfs might be endowed with initially high volatile inventories when compared to the modern Earth and Mars. Thus, it might be possible for such planets to replenish atmospheres desiccated during their host star's lengthy cooling epochs via volcanic outgassing of water-rich mantle material (e.g., Moore & Cowan 2020). While the presence of an external giant planet slightly curtails this incorporation of water-rich material, the effect is fairly minor. Moreover, it is important to note that the more dynamically compact configurations significantly bolster the chances of an HZ planet accreting material from beyond the snow line. While a planetesimal in the solar nebula located at 2.8 au must attain e ≃ 0.66 (an unlikely occurrence in models of terrestrial planet formation; Clement et al. 2019a) to be directly accreted by the young Earth without undergoing a series of fortuitous scattering events, a similarly situated icy planetesimal around an 0.1 M_⊕ M-dwarf at 0.085 au could be accreted by an HZ planet at 0.06 au if it achieves an eccentricity of around 0.3. While this analysis is clearly complicated by the relative location of the snow line during the enhanced luminosity epoch of these star's early lives, this might not pose a significant issue if planetesimals and embryos form further out in the disk before migrating inwards.

We develop a semi-analytic model to analyze whether our extended simulations of debris evolution in TRAPPIST-1, Proxima Centauri, and TOI-700 support the occurrence of late, large, wet impacts on the star's respective HZ planets. The raw data (cumulative fraction of total system debris accreted by each planet of interest) is plotted for our Mercury CPU and GENGA GPU simulations in Figure 6. For the purposes of our subsequent analysis, we focus on the TRAPPIST-1 planets in the conservative HZ (e and f; Kopparapu et al. 2013), though we find that TRAPPIST-1g (situated in the optimistic HZ) accretes a similar number of remnant particles at a similar rate. We also compare our measured bombardment rates to those for the young Earth extrapolated from a similar study of planet formation in the solar system presented in Clement et al. (2019b). Notably, the computations in that paper utilized the same code and fragmentation algorithm employed here. Specifically, we integrate each remnant planetesimal and collisional fragment from all "instability" simulations that finish with the Jupiter-Saturn period ratio less than 2.8 in the presence of the modern solar system (these model the formation of the terrestrial planets occurring in conjunction with the so-called Nice Model instability, and have been demonstrated success when measured against a number of observational constraints: Tsiganis et al. 2005; Clement et al. 2018; Deienno et al. 2018; Mojzsis et al. 2019; Brasser et al. 2020; Nesvorný et al. 2021).

Zoom In Zoom Out Reset image size

The right panel of Figure 6 compares the results of our high-resolution GPU simulations (dashed lines) with our co-added CPU simulations. While the agreement between the respective curves for each planet is reasonable, it is clear that our GPU simulations tend to produce fewer late impacts than our CPU runs (the largest difference can be seen in the yellow curves for TOI-700). We attribute this result to the more realistic size distribution assumed in our GPU simulations (in our CPU simulations, we assign all debris particles equal masses), however, we do not test this hypothesis with additional simulations. Nevertheless, it is important to note that, by focusing on the bombardment curves derived from our co-added CPU simulations in our subsequent analyses, we are essentially assuming a best-case scenario in order to place upper limits on the late bombardment flux in our systems.

From our raw impact chronologies (Figure 6), we fit exponential decay curves and utilize their respective probability distribution functions as inputs for our analytic experiments. Examples of these curves are plotted in Figure 7. We also generate a series of similar curves for eight separate radial origin bins to determine the probability of accreting an object from a particular initial location in the disk at any given time. As expected given the shorter orbital periods of HZ planets around low-mass stars, debris is removed significantly faster in our systems of interest than in the solar system model. In this manner, the planets around TOI-700 (a 0.41 M_⊙ star) experience delayed episodes of bombardment more similar to that of the young Earth than those of our other systems' HZ planets. We also find that our planetesimal populations are removed far quicker than the collisional fragments. We attribute this dissimilarity to the differences in orbital distributions of the debris plotted in Figure 1. Our fragment distributions span a radial range that extends from the exo-asteroid belt region inwards toward the outer edge of the HZ. Conversely, our remnant planetesimal population possess two components: (1) un-accreted objects that survive the giant impact phase on cold orbits in between planets at small semimajor axes and (2) implanted asteroids on relatively stable orbits with large semimajor axes outside of the planetary regime. While the first population is accreted rather easily in the first few tens of kiloyears of our bombardment simulations, the second population is dynamically detached from the orbits of TOI-700d, TRAPPIST-1h, and Proxima Centauri b Thus, in the absence of perturbations from distant giant planets and non-Keplerian accelerations (such as Yarkovsky drift, see discussion in Section 4 and Vokrouhlicky 1998; Bottke et al. 2001; Dencs & Regály 2019; Raymond et al. 2021), the HZ planets in our planetesimal bombardment simulations are unable to access the majority of the asteroid belt material. While the Proxima Centauri system does possess a super-Earth at 1.5 au (outside of this hypothetical asteroidal reservoir), we find it to be too distant and too small to perturb asteroids onto orbits that intersect the HZ.

Zoom In Zoom Out Reset image size

With the rate of material delivery from each initial radial region of our system's presumed primordial terrestrial-forming disk defined, we need only input an SFD for the debris field in order to generate hypothetical bombardment chronologies. While this analysis assumes that these planets formed in a manner similar to that modeled in our terrestrial formation simulations (Section 2.2), we argue in more detail in Section 4 that the generation and stranding of debris should occur in these systems regardless of where the planets form in the disk. Moreover, since our simulations that do not include planetesimals generate collisional fragments at equal rates to our integrations incorporating a sea of such small bodies, it seems reasonable that Earth-size planets implanted in the HZ from the more distant regions of the disk would still be surrounded by debris fields if they experience a few giant impacts with other proto-planets en route to achieving their final orbital architectures.

While the precise SFD of remnant objects is difficult to constrain in the solar system (as it is tied to the size-dependent efficiency of planetesimal formation, e.g., Morbidelli et al. 2009; Johansen et al. 2015), the SFD of the modern asteroid belt serves as a useful proxy for our current purposes (e.g., Sinclair et al. 2020). However, the SFD of asteroids that do not belong to collisional families (Delbo' et al. 2017, 2019) seems to suggest that D ∼ 100 km objects were far more prevalent in the primordial belt than in the modern one (thus potentially evidencing the efficient formation of larger, D ∼ 100-1,000 km planetesimals). Therefore, we also experimented with both debris populations composed entirely of 100 km objects and one that combines the main belt's SFD with an additional 100 km component and found that the different presumed distributions only slightly altered our results.

Figure 8 plots hypothetical remnant planetesimal impact chronologies for our four HZ, Earth-size planets of interest. For the scenario of secondary atmospheric development via small-body delivery to be successful, we require our models produce at least one impact of an object with D ≥ 100 km that originated outside the snow line at t ≥ 500 Myr. Our size requirement is based on an impactor of roughly carbonaceous chondritic composition (and a WMF of ∼0.10 as discussed in Section 2.1) possessing roughly one Earth atmosphere worth of volatiles (∼5 ⨯ 10²⁰ kg). Our time limit is roughly related to the time after star formation that the outer region of an 0.1 M_⊙ star's main sequence HZ is habitable (e.g., Bolmont et al. 2017). While our reference solar system models from Clement et al. (2019b) typically yield several D ≳ 100 km impactors that originated outside of the snow line at t ≳ 500 Myr, such events are nonexistent in our M-dwarf simulations. ⁶ Indeed, the tail-end phase of bombardment is mostly complete in each of our systems of M-dwarf-hosted planets around t ≃ 100 Myr. While collisions persist over gigayear timescales in our systems, our model finds that such events are statistically far more likely to consist of small projectiles that originated inside of the snow line. We also looked at the possibility that Lidov–Kozai (Lidov 1962; Kozai 1962) cycles of high-inclination exo-belt objects might lengthen dynamical lifetimes and prolong bombardment, however the characteristic timescale for such oscillations in our investigated systems is quite short:

$$\begin{eqnarray}&&{\tau }_{\mathrm{Kozai}}\simeq \displaystyle \frac{{P}_{2}^{2}}{{P}_{1}}{(1-{e}_{2}^{2})}^{3/2}.\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

To demonstrate a scenario where significant volatile-rich bombardment persists on ∼500 Myr timescales in our systems, we relax the requirement that impactors originate outside of the snow line when analyzing collisional fragment impacts in our systems. While one might loosely justify this assumption by speculating that fragments might possess enhanced volatile contents by virtue of being ejected from near the surfaces of differentiated proto-planets, we caution the reader that this analysis is the least conservative of our study, and is meant to serve as more of a limiting case than a proof of concept. Figure 9 plots our derived fragment impact chronologies in the same manner as Figure 8, and depicts appreciable bombardment approaching gigayear timescales in our systems (particularly on the TRAPPIST planets). It is important to note that, even under these favorable assumptions, large late impacts are not guaranteed to result in net atmospheric enhancement. Depending on the geometry of the impact (e.g., Svetsov 2007; Shuvalov 2009), such events can result in no volatile accretion, and even erode a significant fraction of the planets' existing atmosphere. Thus, we conclude that, provided these systems do not host distant massive planets capable of dynamically perturbing debris onto HZ-crossing orbits over lengthy timescales (discussed further in the subsequent section: Quintana & Lissauer 2014; Dencs & Regály 2019), the outlook for delayed, impact-driven atmospheric reconstitution of M-dwarf-hosted planets is bleak.

Zoom In Zoom Out Reset image size