# An upper limit on late accretion and water delivery in the TRAPPIST-1 exoplanet system

Sean N. Raymond, Andre Izidoro, Emeline Bolmont, Caroline Dorn, Franck Selsis, Martin Turbet, Eric Agol, Patrick Barth, Ludmila Carone, Rajdeep Dasgupta, et al.

► To cite this version:

Sean N. Raymond, Andre Izidoro, Emeline Bolmont, Caroline Dorn, Franck Selsis, et al. An upper limit on late accretion and water delivery in the TRAPPIST-1 exoplanet system. Nature Astronomy, 2022, 6 (1), pp.80-88. 10.1038/s41550-021-01518-6. hal-03879051

HAL Id: hal-03879051

https://hal.science/hal-03879051v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# An upper limit on late accretion and water delivery in the TRAPPIST-1 exoplanet system

Sean N. Raymond<fnref n="1" />✉, Andre Izidoro², Emeline Bolmont³, Caroline Dorn⁴, Franck Selsis<fnref n="1" />, Martin Turbet⁵, Eric Agol⁵, Patrick Barth⁶,⁷,⁸, Ludmila Carone⁹, Rajdeep Dasgupta<fnref n="1" />,², Michael Gillon¹⁰ and Simon L. Grimm¹¹

The TRAPPIST-1 system contains seven roughly Earth-sized planets locked in a multiresonant orbital configuration<fnref n="1" />,², which has enabled precise measurements of the planets' masses and constrained their compositions³. Here we use the system's fragile orbital structure to place robust upper limits on the planets' bombardment histories. We use *N*-body simulations to show how perturbations from additional objects can break the multiresonant configuration by either triggering dynamical instability or simply removing the planets from resonance. The planets cannot have interacted with more than ~5% of one Earth mass (*M*<sub>⊕</sub>) in planetesimals—or a single rogue planet more massive than Earth's Moon—without disrupting their resonant orbital structure. This implies an upper limit of 10<sup>-4</sup> *M*<sub>⊕</sub> to 10<sup>-2</sup> *M*<sub>⊕</sub> of late accretion on each planet since the dispersal of the system's gaseous disk. This is comparable to (or less than) the late accretion on Earth after the Moon-forming impact⁴,⁵, and demonstrates that the growth of the TRAPPIST-1 planets was complete in just a few million years, roughly an order of magnitude faster than that of the Earth⁶,⁷. Our results imply that any large water reservoirs on the TRAPPIST-1 planets must have been incorporated during their formation in the gaseous disk.

The seven-planet TRAPPIST-1 system is characterized by mean motion resonances between neighbouring planets and three-body (Laplace) resonances between triplets of adjacent planets (Methods and Extended Data Fig. 1)<fnref n="1" />,². Here we use the multiresonant structure of the TRAPPIST-1 system to constrain the planets' bombardment histories after their formation. We use *N*-body simulations to determine the perturbations from additional objects in the system required to break the TRAPPIST-1 resonant chain. Given the survival of the resonant chain from its formation to the present day, we place an upper limit on the integrated mass with which the system has interacted and the impacts that may have taken place.

Starting from a configuration of the system that is consistent with current constraints and stable for billion-year timescales, we invoked gravitational perturbations from additional objects under two limiting assumptions. First, we assumed that a single rogue planetary embryo from the outer disk was scattered inwards towards the planets. The planets undergo stochastic exchanges of orbital energy with the rogue planet during close gravitational encounters,

which can excite the planets' orbits, disrupt the system's multiresonant structure and induce instability. We ran 1,000 simulations varying the mass of a rogue planet in the range 0.001–1 Earth mass (*M*<sub>⊕</sub>). Second, we assumed that a swarm of leftover planetesimals from the outer disk was scattered inwards towards the planets. The planets undergo an exchange of orbital angular momentum with the planetesimals, which causes their orbits to spread out⁸. We ran 140 simulations varying the total mass in 1,000 equal-mass planetesimals between 0.01 *M*<sub>⊕</sub> and 3 *M*<sub>⊕</sub>. Planetesimals either represent leftovers whose orbits were excited by the planets' migration⁹ or objects scattered inwards towards the known planets, perhaps by a more distant planet¹⁰,¹¹. Although other potential sources of impacts exist¹⁰, we focused on inward-scattered rogue objects because they provide analogues to late accretion on Earth⁴,⁵ and have the highest impact probabilities¹⁰ (Methods and Extended Data Table 1). Each system was integrated using the hybrid Mercury integrator¹² for 10 Myr or until it became unstable (Methods, Extended Data Table 2 and Extended Data Figs. 2 and 3). This is two to three orders of magnitude shorter than the lifetime of the system¹³, meaning that our upper limits are conservative, as some systems may become unstable or break their resonances on longer timescales. We performed additional calculations to verify that our results are robust to the initial semimajor axis and periastron distribution of rogue objects and to our choice of the 'best-fit' configurations of the system (Methods, Extended Data Tables 1 and 2, Extended Data Figs. 4 and 5 and Supplementary Table 1)³. We also calculated that star-planet tidal effects did not shift the outer planets' orbits enough to affect our results (Methods and Extended Data Fig. 6).

The TRAPPIST-1 system is dynamically fragile. The planets maintained their multiresonant structure in most simulations containing a rogue planet less massive than ~0.005 *M*<sub>⊕</sub>, but more massive rogue planets were disruptive (Fig. 1). The majority of systems with rogue planet masses between 0.005 *M*<sub>⊕</sub> and 0.02 *M*<sub>⊕</sub> disrupted their resonances but remained stable for the 10 Myr duration of the simulation. More massive rogue planets generally triggered dynamical instability, leading to collisions between planets. The simulation with the most massive rogue planet that maintained its stable resonant structure was just over one lunar mass (~0.01 *M*<sub>⊕</sub>), making this the most massive possible giant impact that could have taken place on any of the TRAPPIST-1 planets at any time since the dispersal of

<footnote n="1">¹Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, Pessac, France. ²Department of Earth, Environmental and Planetary Sciences, Rice University, Houston, TX, USA. ³Observatoire de Genève, Université de Genève, Sauverny, Switzerland. ⁴University of Zurich, Institute of Computational Sciences, Zurich, Switzerland. ⁵Department of Astronomy, University of Washington, Seattle, WA, USA. ⁶Centre for Exoplanet Science, University of St Andrews, St Andrews, UK. ⁷SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, UK. ⁸School of Earth and Environmental Sciences, University of St Andrews, St Andrews, UK. ⁹Max Planck Institute for Astronomy, Heidelberg, Germany. ¹⁰Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Liège, Belgium. ¹¹Center for Space and Habitability, University of Bern, Bern, Switzerland. ¹²e-mail: rayray.sean@gmail.com</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Response of the TRAPPIST-1 system when interacting with rogue planets. Left: Simulations with one rogue planet. Right: simulations with 1,000 rogue planetesimals. Blue circles represent simulations that are consistent with the observed system: the planets' orbits remained stable for the 10 Myr integration and retained their multiresonant configuration. The red symbols denote systems that were stable for 10 Myr (shifted slightly downwards for visibility) but whose resonant structure was disrupted, as quantified using Laplace resonant angles for the outer triplets of planets (Methods). The black circles represent systems that underwent dynamical instabilities, leading to planet-planet collisions. The starting orbits of the rogue planets or planetesimals were chosen assuming that they were scattered inwards from just exterior to the known planets. The results of simulations with different configurations of the TRAPPIST-1 planets produced near-identical results (Methods).</figcaption></figure>

the star's gaseous disk. The three outer planets were most likely to undergo a giant impact, with relative probabilities of 11% (planet f), 28% (planet g) and 31% (planet h). If such an impact occurred, it would have been the only late giant impact on any of the planets.

A population of rogue planetesimals also easily disrupts the TRAPPIST-1 system (Fig. 1). The outer planets are the first to interact with the planetesimals, which are for the most part scattered inwards. The planets' orbits widen and spread out due to the back-reaction of planetesimal scattering⁸. When the planets scatter a total of more than ~0.05 $M_⊕$ in planetesimals, the system's multiresonant structure is disrupted as the angles that characterize mean motion and three-body resonances transition from libration to circulation. We used the evolution of $\phi_5$ (the angle characterizing the observed Laplace resonance between planets f, g and h)² to determine which systems were still in resonance (Methods). After leaving their stabilized resonant state¹⁴, most systems that interacted with more than ~0.1 $M_⊕$ in planetesimals became dynamically unstable and underwent planet-planet collisions. The instability timescale was shorter in systems with more massive swarms of rogue planetesimals, and it is likely that many of the systems from Fig. 1 that left resonance but survived would become dynamically unstable on timescales longer than our 10 Myr integrations.

The maximum integrated bombardment fluxes on the planets are measured in thousandths of an Earth mass (Table 1). These were calculated by imposing a total planetesimal mass of 0.05 $M_⊕$—the maximum mass for which the resonant structure was preserved (Fig. 1)—and measuring the relative impact rates on the planets in simulations that maintained their resonant structure. The spreading of the planets' orbits resulted from the cumulative effect of planetesimal perturbations such that the impacts could have taken place at any time during the system's evolution. Given the short collisional timescales of close-in planetesimals¹⁰, we expect most collisions to happen shortly after the dissipation of the gaseous disk. The impact speeds in our simulations were consistent with previous studies for our chosen impactor population¹⁰,¹¹, with median velocities of 30 km s⁻¹ for planet b decreasing to 13 km s⁻¹ for planet h (Methods and Extended Data Fig. 7).

Our results demonstrate that the TRAPPIST-1 planets grew an order of magnitude faster than Earth (Fig. 2). Cosmochemical constraints⁶,⁷ and accretion models¹⁵ indicate that a large portion of Earth's growth took place after the dissipation of the Sun's gaseous disk. Mars's growth was close to complete by the time the disk dissipated¹⁶, and Earth's constituent embryos were probably roughly

<table><caption>**Table 1** | Upper limits on the total mass in planetesimals that collided with each of the TRAPPIST-1 planets and the maximum amount of water delivered by those impacts</caption><thead><tr><th>Planet</th><th>Orbital radius (au)</th><th>Maximum bombardment mass (M<sub>⊕</sub>)</th><th>Maximum water delivered (Earth oceans)</th></tr></thead><tbody><tr><td>b</td><td>0.0115</td><td>0.00038</td><td>0.15</td></tr><tr><td>c</td><td>0.0158</td><td>0.0015</td><td>0.64</td></tr><tr><td>d</td><td>0.0223</td><td>0.0016</td><td>0.68</td></tr><tr><td>e</td><td>0.0293</td><td>0.0035</td><td>1.54</td></tr><tr><td>f</td><td>0.0385</td><td>0.008</td><td>3.41</td></tr><tr><td>g</td><td>0.0469</td><td>0.018</td><td>7.62</td></tr><tr><td>h</td><td>0.0620</td><td>0.012</td><td>6.16</td></tr></tbody></table>

Perfect mergers were assumed in calculations. Here we have adopted a maximum value of 0.05 $M_⊕$ in scattered rogue planetesimals (Fig. 1) and assumed that they contain 10% water by mass, typical of the wettest carbonaceous chondrite meteorites¹⁰. One 'ocean' of water is equivalent to the mass of the water on Earth's surface (2.3×10⁻⁴ $M_⊕$).

Mars-sized at that time¹⁵. Earth's final, Moon-forming impact took place ~50–100 Myr later⁷. The abundance of highly siderophile elements in the mantle and crust suggest that Earth accreted ~0.003–0.007 $M_⊕$ in planetesimals in the 4.45 Gyr between the Moon-forming impact and the present day⁴,⁵; this is generally called late accretion. Our dynamically derived upper limits on late accretion on the TRAPPIST-1 planets are more stringent than those for the Earth because the time window extends all the way back to the end of the gaseous disk phase. Nonetheless, we constrain late accretion on the three innermost TRAPPIST-1 planets to have been substantially lower than on Earth (Table 1). The four outer TRAPPIST-1 planets could not have accreted more than ~1% of their masses during this interval, comparable to upper limits on Earth's late accretion¹⁷. The TRAPPIST-1 planets must therefore have been fully assembled during the few-million-year-long gaseous disk phase. This is in line with successful formation models, which invoke gas-driven orbital migration¹⁸⁻²⁰. Migration seems to be essential to produce an orbital configuration with a long dynamical lifetime¹⁴. The starting orbital distance of migrating planets is unconstrained, although models expect large embryos to form fastest just beyond the snow line¹⁸,²¹. The low occurrence rate of multiresonant systems
<figure><img src="image_4.png" /></figure>

**Fig. 2 |** A possible growth history for the TRAPPIST-1 planets compared with constraints for Solar System objects. Characteristic growth timescales for Vesta, Mars and Earth (and for the Moon-forming impact) come from analyses of different radioactive systems<sup>6,7,16</sup>. After the dispersal of the gaseous disk, Earth's growth was probably characterized by giant impacts<sup>15</sup> (including the Moon-forming impact) that are not reflected in the Hf-W and U-Pb model curves<sup>5</sup>. We used a simulation by Izidoro et al.<sup>9</sup> that formed a close-in, long-term stable multiresonant chain of planets. The growth curves of seven planets with masses between 0.4 M<sub>⊕</sub> and 1.5 M<sub>⊕</sub>—characterized by an early phase of pebble accretion followed by giant impacts—are included as analogues for the accretion of the TRAPPIST-1 planets (Methods). Each planet's growth curve is shown in a shade of blue, with darker colours representing more distant orbits. Although this is only one illustrative example among many<sup>20,21</sup>, the growth of the TRAPPIST-1 system must have completed by the end of the gaseous disk phase, the duration of which is constrained by observations of the occurrence rate of hot dust around stars in clusters of different ages<sup>85</sup>.

suggests that most resonant chains undergo dynamical instabilities after the dispersal of the disk<sup>9,22</sup>; TRAPPIST-1 therefore represents a rare surviving resonant chain.

To investigate the role of late water delivery to the TRAPPIST-1 planets, we performed new calculations to estimate the planets' water contents. We used interior structure<sup>23</sup> and atmospheric<sup>24</sup> models constrained by the latest planetary mass and radius estimates<sup>3</sup>. Given that inferred water contents are subject to interior degeneracies<sup>3,25-28</sup>, we estimated water mass fractions assuming five different rocky interior models (Methods). Across all of our models we found that planets b, c and d are likely to be volatile-depleted, whereas the outer planets (e–h) are more likely to have significant water contents (consistent with previous models<sup>3,28</sup>). Yet although our mean water contents were ~1–10 wt% for the outer planets, 2σ error bars reach down to zero in at least one model for each planet (Fig. 3). We therefore cannot be certain of a water-rich composition for any of the planets. Planet g is our strongest candidate for being water-rich, with just one model reaching zero water content at the 2σ level.

Late accretion may have delivered the full present-day water budgets of planets b, c and d (Fig. 3). These upper limits assume that late impactors contained 10% water by mass, typical of the wet-test carbonaceous chondrite meteorites<sup>29</sup>. They also assume perfect merging and no water loss, although only a fraction of material is likely to be accreted (as little as 1–10% for the highest-velocity impacts)<sup>10,30</sup>. Taking the upper limits at face value, late accretion can deliver enough water for planets e–h to each have bulk water contents similar to Earth's value of ~10<sup>-3</sup> by mass<sup>31</sup>. However, late accretion falls an order of magnitude short of delivering enough water to

<figure><img src="image_5.png" /></figure>

**Fig. 3 |** Late water delivery to the TRAPPIST-1 planets. The error bars indicate the 1σ (thick bars) and 2σ (thin bars) confidence intervals (encompassing 68% and 95% probabilities) for estimates of the bulk water content of each planet for each of the five models discussed in Methods (note that certain models do not apply to the planets b, c and d). Each circle is proportional to the planet's size, and its vertical position represents an upper limit on the water delivered during late accretion. The grey vertical lines show snapshots of the location of the inner edge of the habitable zone calculated with 3D simulations<sup>86</sup>, which swept inwards in time as the central star evolved to the main sequence<sup>33,34</sup>. Earth's surface water budget of one ocean is shown as the dashed horizontal line.

explain the mean water contents from our models (Fig. 3). If future work confirms water contents at the weight per cent level or higher for these planets, then a separate source of water would be required beyond late accretion. The most natural source is water incorporated during the planets' accretion, as predicted by formation models<sup>9,18,20,21,26</sup>. Long-range orbital migration would be required for planets e–h to acquire significant amounts of water. Some constituent pieces of each planet must have originated beyond the snow line, which was more distant during the star's brighter early phases<sup>26,32</sup>. Two viable possibilities are (1) that the cores of the planets originated just past the snow line and migrated inwards<sup>18</sup> or (2) that the snow line swept inwards, allowing rocky cores to efficiently accrete icy pebbles and migrate inwards<sup>9</sup>.

Each of the planets may have lost a significant amount of water after accretion. High-velocity planetesimal impacts may cause water (and atmospheric) loss<sup>10,30</sup>. On long timescales, impacts may cause a net water gain or loss on each planet depending on the balance between stripping impacts and lower-velocity impacts<sup>10</sup>. Given the radial dependence of impact speed (Methods), the outer planets are more likely to accumulate volatiles<sup>10</sup>. In addition, as an ultra-cool star TRAPPIST-1 evolved very slowly onto the main sequence; the habitable zone swept inwards such that even the planets in the present-day habitable zone spent time with much higher equilibrium temperatures<sup>33,34</sup>. During this time, water in the planets' atmospheres would have been vulnerable to loss driven by ultraviolet radiation from the central star (especially for the inner planets)<sup>35,36</sup>. Models of the early evolution of the TRAPPIST-1 planets that also include the magma ocean phase<sup>27,37</sup> find that each planet may have lost a few to 20 Earth oceans of water<sup>27</sup>. Calculations that make assumptions designed to maximize water loss still find an upper limit of 25 Earth oceans lost in the interval before planets e–h reached the habitable zone<sup>27,36</sup>. This is an order of magnitude lower than the water contents allowed by observationally constrained internal structure modelling (Fig. 3) for each planet, compatible with a relatively water-rich formation scenario.
Our work highlights the exceptional value of multiresonant planetary systems, only a handful of which are well-characterized, notably TOI-178³⁸ and Kepler-223³⁹. Not only are these systems ideally suited for mass measurement via transit-timing variation analysis³, their dynamical fragility allows us to tightly constrain their bombardment histories. Our results demonstrate that the TRAPPIST-1 system satisfies the International Astronomical Union's criteria for planethood. The International Astronomical Union's resolution B5⁴⁰ defines a planet as an object that (1) is in orbit around the Sun, (2) is massive enough to assume hydrostatic equilibrium in a nearly round shape and (3) has cleared the neighbourhood around its orbit. Assuming criteria 1 and 2 to be satisfied using observations of the central star¹ and internal structure modelling of the planets²⁴,²⁵, our simulations show that the inner TRAPPIST-1 system has been cleared of rogue material to a stringent limit.

## Methods

**Current state of knowledge of the TRAPPIST-1 central star and the system's resonant structure.** TRAPPIST-1 has become iconic among exoplanet systems. The sizes of TRAPPIST-1 planets are between 0.75 and 1.15 Earth radii ($R_⊕$), and their masses have been measured with a precision of 3–6% using transit-timing variation analysis³. The system's two innermost planets were discovered in 2016⁴¹ and the next five in 2017¹. The central star is an ultracool dwarf star of spectral type M8 located 12 pc from the Sun¹,⁴¹. A number of analyses have constrained the star's physical characteristics¹³,⁴²⁻⁴⁴. It has a bolometric luminosity of roughly $5.5 \times 10^{-4} L_⊕$ (where $L_⊕$ is the luminosity of the Sun; ref. ⁴⁴), a mass of ~0.09 solar masses ($M_⊕$) and a size of ~0.12 solar radii ($R_⊕$). Its age has proven hard to pin down. While one study suggest a relatively old age of ~7.6 Gyr (ref. ⁴⁵), a recent study¹³ found analogue stars that are as young as 0.5 Gyr and as old as 10 Gyr. While this remains an interesting astronomical detective story, even a 'young' age of 0.5 Gyr is far longer than our N-body simulations. As discussed in the main text, our upper limits on late accretion are conservative, as an older system age implies a longer timescale for instability to occur (among the systems that were stable for our 10 Myr simulations).

TRAPPIST-1's multiresonant structure is one of its most remarkable attributes. Judging from orbital period ratios, the sequence of orbital resonances from planets b:c outward is: 8:5, 5:3, 3:2, 3:2, 4:3, 3:2 (ref. ²). Extended Data Fig. 1 shows the evolution of a number of resonant angles in a 100 Myr integration of our fiducial best-fit TRAPPIST-1 system. The innermost pair of resonances are rarely seen in simulations of resonant capture, which tend to produce first-order resonances⁹,²⁰,⁴⁵⁻⁴⁷. Papaloizou et al.¹⁵ proposed that each of these planet pairs was in fact captured into 3:2 or 4:3 resonance; that is, first-order resonances. This scenario requires fast circulation of the longitudes of periastron for planets b and c via damping of their free eccentricities, which happens naturally due to tides. This remains a plausible scenario to be confirmed with future studies combining orbital and tidal dynamics.

While three-body (Laplace) resonance is confirmed for adjacent triplets of planets in the TRAPPIST-1 system (for example, see fig. 25 in Agol et al.¹), current constraints cannot unequivocally guarantee that two-body resonant angles are in libration, rather than circulation. This is because three-body resonant angles (denoted $\phi$) depend only on mean longitudes, whereas two-body angles ($\theta$) explicitly depend on the longitudes of periastron, which are poorly constrained by the transit times and transit timing variations. We strongly suspect that two-body angles are indeed librating, and below we show empirically that both the two-body and three-body resonances seem to break at roughly the same time as the planets are spread apart by dynamical interaction with rogue planetesimals (Extended Data Fig. 3). Even if the two-body frequencies are wide in terms of resonance, at small eccentricity the longitudes of periastron can undergo fast circulation, causing two-body angles to librate.

**Simulations.** Our simulations started from a configuration of the TRAPPIST-1 system derived from the most comprehensive analysis of the system so far¹. We integrated the system on its own and verified that it is dynamically stable for at least 1 Gyr. Extended Data Table 2 lists the orbital parameters of our starting system (as well as those of three other best-fit systems that we tested; see the 'Testing the effect of our chosen configuration for the Trappist-1 planets' orbits' section). Although the system is two-dimensional, we gave non-zero inclinations to rogue embryos and planetesimals to ensure that the systems as a whole were three-dimensional.

The initial conditions of our simulations were built on the assumption that one or more rogue objects were scattered inwards towards the TRAPPIST-1 planets. Rogue objects may be thought of as analogues of the leftover planetesimals that are thought to have impacted the Earth and other terrestrial planets during the late accretion phase⁴,⁵,³⁷,⁴⁸⁻⁵⁰. If rogue objects were scattered inwards from more distant orbits, they would first interact with the outer planets, and their apoastron distances would correspond to the next-most-distant planet in the system (as

is the case for Jupiter-family comets in the Solar System⁵¹), which is poorly constrained from observations⁵²,⁵³. In addition to the TRAPPIST-1 planets listed in Extended Data Table 2, we therefore included either one rogue planet with a randomly chosen mass between 0.001 and 1 $M_⊕$, or 1,000 equal-mass rogue planetesimals whose masses totalled between 0.01 and 3 $M_⊕$. The initial semimajor axis of each rogue object was randomly chosen between 0.07 and 0.1 au, with a periastron distance randomly selected to lie between 0.038 and 0.062 au (between the orbits of planets f and h), to reflect an inward-scattering origin from which the outermost planets are encountered first³¹. We performed additional sets of simulations in which the semimajor axes of rogue planetesimals were randomly chosen between either 0.2 and 3 au or 0.7 and 1 au (see the 'Testing the effect of our chosen configuration for the Trappist-1 planets' orbits' section). The initial inclinations of rogue objects were randomly assigned between 0 and 10°, which matches the inclination distribution of small leftover bodies in simulations that form close-in, TRAPPIST-1-like systems from pebble accretion and migration⁵; we inherently assumed this distribution to apply regardless of the exact migration scenario¹⁸⁻²²,⁴⁷,⁵⁴⁻⁵⁹. In simulations with 1,000 planetesimals, the planetesimals felt the planets' gravity but not each other's. Even though all planetesimals were included from the start of the simulation, the absence of mutual gravity allowed us to roughly mimic the interactions of individual rogue planetesimals with the planets and estimate their cumulative effects.

Each simulation was integrated for 10 Myr but stopped if an instability occurred, marked by a collision between planets. We used the hybrid integrator in the publicly available Mercury integration package¹² with a time step of 0.1–0.15 days, which tests showed was sufficient for long-term integration without significant numerical error. Objects were considered to be ejected if they strayed more than 50 au from the central star, and were considered to have hit the star if they passed closer than 0.005 au, the rough size of our Sun (although note that TRAPPIST-1 is only ~0.12 $R_⊕$; ref. ³,⁴⁴). In terms of calculating collisions, the TRAPPIST-1 planets were assigned their measured sizes (Extended Data Table 2), whereas rogue planets and planetesimals were assigned sizes assuming a bulk density of 3 g cm⁻³. Collisions were treated as inelastic mergers conserving linear momentum.

**Characterizing and breaking orbital resonances.** We used a variety of resonant angles to characterize the planets' resonances. The TRAPPIST-1 resonant chain involves a series of mean motion and Laplace resonances⁵. From the inside outwards, the planet pairs are in or near 8:5, 5:3, 3:2, 3:2, 4:3 and 3:2 mean motion resonance. In addition, each set of three planets is in Laplace resonance.

For each planet pair in mean motion resonance $p+q:p$, resonant angles $\theta_i$ are of the form

$$ \theta_{1,2} = (p+q)\lambda_1 - p\lambda_2 - q\omega_{1,2} \quad (1) $$

where $\omega$ are longitudes of pericentre, $\lambda$ are mean longitudes (equal to the sum of $\omega$ and the mean anomalies) and subscripts 1 and 2 refer to the inner and outer planet, respectively. Resonant angles measure the angle between the planets at conjunction; if any angle librates, rather than circulates, then the planets are in resonance. Different resonances have different numbers of resonant arguments, involving various permutations of the final terms in equation (1). The order of a resonance is given by $q$.

In our analysis, we calculated the relevant $\theta$ values for each planet pair. We focused on the outermost pair of planets, as we found that their 3:2 resonance was often the first to break. We also calculated three-body (Laplace) resonant angles. There are five such angles in the system, one for each adjacent triplet of neighbouring planets. Using subscripts appropriate for the TRAPPIST-1 planets, these are (see table 1 in Luger et al.²):

$$ \begin{aligned} \phi_1 &= 2\lambda_b - 5\lambda_c + 3\lambda_d \\ \phi_2 &= \lambda_c - 3\lambda_d + 2\lambda_e \\ \phi_3 &= 2\lambda_d - 5\lambda_e + 3\lambda_f \\ \phi_4 &= \lambda_e - 3\lambda_f + 2\lambda_g \\ \phi_5 &= \lambda_f - 2\lambda_g + \lambda_h \end{aligned} \quad (2) $$

Extended Data Fig. 1 demonstrates the multiresonant nature of the initial conditions for the TRAPPIST-1 system that we used in our simulations. It shows the evolution of different mean motion and Laplace resonant angles over a 100-Myr integration of the fiducial system (orbital parameters are in Extended Data Table 2). Each of these angles is in a state of libration, meaning that its evolution is constrained to a limited range of possible values. This indicates a repeated alignment of the planets during their conjunctions, and is the definition of orbital resonance⁶⁰.

Extended Data Fig. 2 presents three examples of the evolution of resonant behaviour from simulations with 1,000 rogue planetesimals. In the simulation interacting with the lowest mass in planetesimals (a total of 0.0213 $M_⊕$; left panel), the planets maintained their resonant configuration with no noticeable difference compared with the fiducial system from Extended Data Fig. 1. When the mass in planetesimals was doubled to 0.0429 $M_⊕$ (centre panel), the resonances were
strained. The mean motion resonances between planets d and e, e and f, and f and g were maintained (top three centre panels). However, the 3:2 resonance between planets g and h was compromised, as the resonant angle $\theta_2$ no longer showed clear libration, but instead circulates. However, all three three-body resonances were maintained, as shown by the continued libration of the Laplace angles $\phi_3$, $\phi_4$ and $\phi_5$ (bottom three centre panels). Although the exact centre of libration underwent shifts during the integration, we consider the resonance to have been maintained to a sufficient degree that, if the system were detected and observed over a relatively short temporal arc, it would appear to be in resonance. Resonances were completely disrupted in the system interacting with just a slightly higher mass in planetesimals, 0.0645 $M_⊕$ (right panel). In that case, all of the resonant angles associated with mean motion and Laplace resonances transitioned from libration to circulation within roughly 1 Myr of the start of the integration. Although the system remained stable for the 10 Myr length of the integration, it is possible that it would be unstable on longer timescales¹⁴.

Extended Data Fig. 3 shows the disruption of the resonances in a simulation in which our fiducial TRAPPIST-1 planets interact with a rogue planetesimal swarm of 0.0967 $M_⊕$. The widening of the outer planets' orbits was small, and the changes in the planets' semimajor axes could be measured in ten-thousandths of an astronomical unit. The resonances were maintained for the first several thousand years. The 3:2 mean motion resonance between planets g and h was disrupted once the relative separation of the two planets increased by roughly 0.0002 au. The 4:3 resonance between planets f and g was likewise destroyed when those planets' relative separation increased by ~0.00015 au, which happened at a later time than the disruption of the 3:2 resonance between planets g and h because of the different rates of relative drift. The $\phi_5$ Laplace angle for planets f/g/h was disrupted at roughly the same time as the 3:2 mean motion resonance between planets g and h, and the $\phi_4$ Laplace angle for planets e/f/g was disrupted on a similar timeframe.

To understand why a given degree of planetesimal scattering disrupts orbital resonances, we must first calculate the width of a resonance, and then determine how much mass a planet must scatter to shift its orbital radius enough to leave the resonance. Following the method of Hadden and Lithwick⁶¹ (see their equations (5), (11) and (15)), we estimated the width of the 3:2 mean motion resonance between planets g and h as $\Delta a \approx 2.4 \times 10^{-4}$ au. This is in good agreement with the evolution of the system shown in Extended Data Fig. 3, in which the 3:2 resonance between planets g and h was indeed disrupted once their relative semimajor axes increased by ~2 $\times 10^{-4}$ au.

How much mass in planetesimals must be scattered to remove a pair of planets from resonance? The interactions between the planets and rogue planetesimals involve the exchange of orbital angular momentum, L. The angular momentum of a planet of mass $M_p$ is defined as:

$$ \text{Angular momentum} = M_p \sqrt{GM_* a} (1 - \cos(i) \sqrt{1 - e^2}), \quad (3) $$

where G is the gravitational constant, i is the orbital inclination and $M_*$ is the stellar mass.

If a planetesimal of mass m transfers all of its angular momentum to a planet of mass $M_p$, that planet's orbit will shift by $\delta a/a \approx 2m/M_p$ (neglecting terms related to eccentricity and inclination, and assuming a common semimajor axis for the planetesimal and planet). Based on this calculation, it should require ~0.005-0.01 $M_⊕$ to move planet h the width of its common 3:2 resonance with planet g. This is roughly the minimum mass in a single rogue planet that destroys the system's resonant structure (Fig. 1 and Extended Data Fig. 5). However, a total mass of rogue planetesimals of ~0.05 $M_⊕$ is needed to break the resonance. The reason for this is that planet h does not orbit its star in isolation. Planet h does not take all of a typical planetesimal's angular momentum, but instead scatters the planetesimal inwards to interact with planet g. The backreaction of the scattering causes planet h's semimajor axis to shift slightly outwards. Subsequently, the same process is repeated as the planetesimal interacts with planet g, causing its semimajor axis to also shift slightly. Thus, the planets migrate outwards together, as is clearly seen in the example from Extended Data Fig. 3. This process is not perfectly balanced; in time the planets' orbits do drift apart, but this mutual drift is significantly slower than each of their individual drift rates. Thus, the planets actually drift several times their resonant width outwards before they separate enough to leave the resonance. For this reason it takes substantially more than 0.01 $M_⊕$ to disrupt the resonance. Hence, the system's resonant structure requires less mass to be disrupted by a single rogue planet than by a swarm of rogue planetesimals.

Our upper limits on the mass with which the TRAPPIST-1 system could have interacted are conservative. Once a system loses its multiresonant structure or goes unstable, it is irreversible: the system cannot regain its resonances (after the end of the gaseous disk phase—tidal evolution is not strong enough; see the 'Tidal evolution' section) and it is forever inconsistent with the present-day system. Yet a fraction of systems that remained in resonant chains for the duration of our 10 Myr simulations could lose their resonances or become unstable on longer timescales. A hypothetical system that loses its resonant structure on a longer timescale could not have interacted with more total mass than our upper limits. By using non-self-interacting particles, our simulations constrain the cumulative effect of a given mass in rogue planetesimals. For example, a system that interacted

with 0.02 $M_⊕$ in rogue planetesimals during the first 10 Myr of its history would retain its multiresonant structure. Given the upper limit of 0.05 $M_⊕$ from Fig. 1, this system could not interact with more than an additional 0.03 $M_⊕$ on longer timescales without losing its resonances.

**Testing the effect of our initial orbital distribution of rogue objects.** Our choices for the orbital distribution of rogue planets and planetesimals were motivated by the well-studied dynamics of Jupiter-family comets²¹ and of leftovers from terrestrial planet formation¹⁶, and from the orbital distribution of analogous leftover material in simulations of the formation of close-in systems like TRAPPIST-1⁹. Yet it is worth questioning these choices, and exploring how robust our conclusions are to other plausible populations of rogue objects (see Kral et al.¹⁰ for a discussion).

We therefore performed two additional sets of simulations in which the orbits of the planets were unchanged. Rogue planetesimals were included with the same inclination distribution (randomly selected from 0–10°; derived from simulations of the formation of TRAPPIST-1-like systems⁹) and periastron distributions (randomly chosen between 0.038 and 0.062 au) as in our fiducial set. The reason for keeping these parameters identical was to isolate the importance of the rogue objects' semimajor axes.

The new sets of simulations were set up as follows:

* A set of 50 simulations in which the semimajor axes of rogue objects were randomly chosen between 0.2 and 0.3 au, roughly three times larger than in our fiducial set. 50 simulations were performed with a total mass of rogue planetesimals evenly chosen in the range 0.01–1 $M_⊕$.

* In the second set, the semimajor axes of rogue objects were randomly chosen between 0.7 and 1 au, a full order of magnitude larger than in our fiducial set. 50 simulations were performed with a total mass of rogue planetesimals evenly chosen in the range 0.01–1 $M_⊕$.

All other integration parameters (time step, ejection radius, stellar mass and so on; see the 'Simulations' section) were the same as for our fiducial set of simulations.

Extended Data Fig. 4 compares the outcomes of these new sets of simulations with our fiducial set. There is no clear difference between simulations with $a_{rogue} = 0.07-0.10$ au and those with $a_{rogue} = 0.2-0.3$ au (for the same periastron distribution). In both cases, a total mass of planetesimals larger than ~0.05 $M_⊕$ universally disrupted the system's multiresonant architecture. In contrast, simulations with $a_{rogue} = 0.7-1.0$ au the critical mass of planetesimals increased to ~0.1 $M_⊕$.

The central difference between the different sets of simulations is the fraction of planetesimals that is ejected from the system, rather than accreted by a planet. When considering only simulations that retained their multiresonant structure, 7.1% of planetesimals in our fiducial simulations (with $a_{rogue} = 0.07-0.10$ au) were scattered beyond 50 au and considered to be ejected from the system. This fraction increased dramatically in the new sets of simulations: 41.6% of planetesimals in simulations with $a_{rogue} = 0.2-0.3$ au were ejected, and 76.5% of planetesimals were ejected in simulations with $a_{rogue} = 0.7-1.0$ au.

The tendency for increased ejection of planetesimals for larger $a_{rogue}$ can be understood from simple dynamical considerations (see also the discussion in section 3.2 of Kral et al.¹⁰). To be ejected, rogue objects must be gravitationally scattered by the planets, and the maximal change in velocity that a planet can impart is the escape speed from its surface²². Orbital energy scales with the apoastron distance as ~$GM_*/Q$, such that a smaller velocity kick is required to eject an object on a more distant orbit. For an Earth-mass, Earth-sized planet, comets can be efficiently unbound for Q larger than roughly 0.2 au (see equation (3) in Kral et al.¹⁰). Our simulations match this expectation, as rogue objects are preferentially accreted in our fiducial simulations with $a_{rogue} = 0.07-0.10$ au, preferentially ejected in simulations with $a_{rogue} = 0.7-1$ au and are roughly evenly split between accretion and ejection in simulations with $a_{rogue} = 0.2-0.3$ au.

Compared with the fiducial set, simulations with larger values of $a_{rogue}$ allow for less late accretion on each planet (Extended Data Table 1). These upper limits were calculated by tabulating the impact rates on each of the planets in simulations that preserved their multiresonant structure, then multiplying by the maximum allowed total planetesimal mass. The difference between the fiducial set and the new sets is largest for the inner planets. This is because the initial periastron distances of planetesimals only crossed those of planets f through h; close encounters with the planets are required for a planetesimal's orbit to shrink its close approach distance. Yet for planetesimals with large $a_{rogue}$ such encounters often lead to ejection such that, compared with planetesimals with small $a_{rogue}$, a smaller fraction of objects enters the innermost parts of the system to potentially collide with the inner planets.

These sets of simulations demonstrate the robustness of the upper limits on late accretion that were derived using our fiducial initial conditions. Although other orbital distributions of rogue planetesimals may imply higher values for the total mass of planetesimals that could have interacted with the planets (Extended Data Fig. 4), the total accreted mass is universally lower (first three columns in Extended Data Table 1). In addition, the higher impact speeds of more distant planetesimals may often lead to erosion, rather than accretion¹⁰.

We performed another set of simulations in which the planetesimals were on even more extreme orbits. The planetesimals' semimajor axes were
randomly chosen in the range $a_{\text{rogue}}=0.7-1$ au, and their periastron distances were chosen following a logarithmic distribution between 0.005 au and 0.063 au; the planetesimals' eccentricities thus spanned from 0.91 to 0.995. As before, we performed 50 simulations varying the total mass in rogue planetesimals between 0.01 and $1 M_{\oplus}$.

Extended Data Fig. 4 (bottom right panel) presents the results of this set of simulations with rogue planetesimals on extreme orbits, similar to those modelled by Kral et al.¹⁰. The most striking features are the absence of unstable orbits and the larger critical mass of rogue planetesimals below which the systems maintained their multiresonant structure. This can be understood as a consequence of the fact that planetesimals on extremely high eccentricities carry very little orbital angular momentum, such that a relatively large mass is needed to move the planets' orbits out of resonance.

Despite the increased mass with which the planets interacted, the amount of late accretion did not change dramatically (Extended Data Table 1). Part of this was a consequence of ejection: as in other sets of simulations with large semimajor axes, a large fraction (77%) of planetesimals was ejected, reducing the reservoir of potential impactors. For planets f, g and h, the upper limits on late accretion in this set of simulations is significantly lower than in our fiducial set of simulations, with maximum late accreted masses comparable to Earth's value of ~0.005 $M_{\oplus}$ (refs. ⁴,⁵). For planet e the upper limit on late accretion was similar in both cases (~0.003–0.004 $M_{\oplus}$). For planets b, c and d, however, the upper limits were significantly higher in simulations with super-eccentric rogue planetesimals, with values up to 0.008 $M_{\oplus}$ (for planet b; Extended Data Table 1). This is a consequence of the much higher number of planetesimals initially placed on orbits crossing the orbits of the inner planets (due to the logarithmic distribution of periastron distances), and the fact that impact rates are higher for the inner planets¹⁰.

One might wonder whether a population of planetesimals at even higher semimajor axis could in principle increase our derived upper limits on late accretion (and perhaps on a longer timescale if the semimajor axes are sufficiently large). We do not think that this is the case. In our extra sets of simulations the upper limits on accretion did not increase for any planet (see Table 1). (The exception is the set of simulations with super-eccentric orbits from Extended Data Table 2 in which the increase for planets b–d comes from a change in the initial periastron distribution of rogue planetesimals and not from an increase in the semimajor axes). Further increasing the semimajor axes of rogue planetesimals would further decrease their specific orbital angular momentum and increase the upper limits on mass with which the planets could have interacted, but this effect would be counterbalanced by the much-increased ejection efficiency and the lower accretion efficiency, as quantified by Kral et al.¹⁰. Of course, the origin of a massive flux of objects crossing the planets' orbits is a separate issue that was also addressed by Kral et al.¹⁰.

Upper limits on late water delivery from rogue planetesimals may be affected if planetesimals originating farther from the star have higher water contents. Comets are generally assumed to be half water by mass⁶³, although the bulk water contents of the best-studied comet (67P/Churyumov-Gerasimenko) is roughly constrained to ~25% or lower (see fig. 5b in Choukroun et al.⁶⁴). If planetesimals from 0.7–1 au have bulk water contents of 25–50%, then the corresponding upper limits could modestly increase relative to our fiducial case—for instance, for planets g and h in the simulations presented in Table 1 and for all of the planets in the simulation with extreme eccentricities (Extended Data Table 2). However, the higher impact velocities of planetesimals originating at 0.7–1 au could reduce the actual delivery rate of water through imperfect accretion and impact volatilization¹⁰,³⁰.

**Testing the effect of our chosen configuration for the TRAPPIST-1 planets' orbits.** One might worry that our results were affected by the exact choice of initial conditions—in other words, by the particular 'best-fit' system that we chose. Of course, this is just one of a population of configurations that are consistent with the current observational dataset. To address this concern, we ran similar batches of simulations with rogue planets and planetesimals in three other best-fit configurations that match all current constraints for the TRAPPIST-1 system³. Long-term simulations showed that each of these three configurations is stable in isolation for at least 1 Gyr.

For each of these three systems (full configurations listed in Supplementary Table 1), we ran two sets of simulations that were very similar to our fiducial set:

* A set of 500 simulations that included a single rogue planet with a mass in the same range as for our fiducial set ($M_{\text{rogue}}=0.001-1 M_{\oplus}$) and the same orbital distribution as described above.

* A set of 50 simulations that included 1,000 planetesimals, with a total mass of 0.01–1 $M_{\oplus}$ and the same orbital distribution as for the fiducial simulations.

The integration method was the same as in the main set of simulations, utilizing the hybrid Mercury integrator¹².

Extended Data Fig. 5 shows the outcome of these sets of simulations, which are virtually identical to our fiducial set (referred to as Set 1). The critical mass in planetesimals above which resonance was destroyed was always in the range 0.04–0.05 $M_{\oplus}$. For the case of a single rogue planet, the results of our additional simulations are again very close to those for our fiducial set. In fact, the results are slightly more restrictive in terms of the maximum mass of a late rogue planet: not one simulation maintained its multiresonant structure with a rogue planet more

massive than 0.01 $M_{\oplus}$. From all sets together, only a single simulation remained in a multiresonant configuration after interacting with a rogue planet more massive than 0.01 $M_{\oplus}$ and even in that case (in Set 1), the rogue planet was barely more massive than the Moon ($M_{\text{Rogue}}=0.0167 M_{\oplus}=1.4 M_{\text{Moon}}$). That simulation had an unusual evolutionary path, as the rogue planet passed farther than 50 au from the star and was considered to be ejected from the system. The rogue planet was ejected in only ~2% of our simulations. We therefore consider our upper limit for a single, late, giant impactor of a lunar mass to be robust.

**Distribution of collision parameters.** We calculated the collisional parameters for each planetesimal collision that took place within simulations in which the planets remained in resonance. Almost all collisions took place within the first 100 kyr of each simulation, with a median collision timescale of ~5–10 kyr for each planet, consistent with previous work¹⁰. Extended Data Fig. 6 shows the distribution of these impacts for our fiducial set of simulations (top panel, with initial orbital semimajor axes for rogue planetesimals of $a_{\text{rogue}}=0.07-0.1$ au), as well as for the simulations with larger initial semimajor axes described above (see the 'Testing the effect of our initial orbital distribution of rogue objects' section and Extended Data Fig. 4).

There are two trends in Extended Data Fig. 6. The first is that the impact speed increases monotonically for planets closer to the star. Among just the fiducial set of simulations (top panel), median impact speeds onto the planets varied from 30 km s⁻¹ for the innermost planet (b) to 13 km s⁻¹ for the outermost (h), which is consistent with results of previous studies for similar impacting populations¹⁰,¹¹. The second trend is that impact speeds are modestly higher for larger initial orbital semimajor axis of the rogue objects. For instance, the median impact velocities on planets g are 17.2, 20.3 and 21.1 km s⁻¹ for increasing values of $a_{\text{rogue}}$ (from the top to bottom panels).

Both of these trends can be explained with simple analytical arguments. If we assume that impacts happen when a rogue object is at periastron, the rogue object's velocity is roughly the local Keplerian velocity $v_K = \sqrt{GM_*/a}$ times $\sqrt{1+e}$, where e is the rogue object's orbital eccentricity. The expected encounter speed between a rogue object at periastron and a planet on a roughly circular orbit at the same orbital radius is therefore $v_K(\sqrt{1+e}-1)$ (assuming parallel velocity vectors), and this value should represent a typical impact speed¹⁰. This is indeed the case in our simulations. The scaling of $v_K$ with orbital radius can explain why, within a given system, impacts on the inner planets occur at higher speed than on the outer planets. The trend of higher impact speeds for larger $a_{\text{rogue}}$ can be explained by the second term, as rogue objects had fixed initial periastron distances such that larger $a_{\text{rogue}}$ corresponds to larger eccentricity.

Kral et al.¹⁰ found a bimodal distribution of impact velocities on the TRAPPIST-1 planets (see their section 3.4). Their lower-velocity impacts came from the same type of impacts seen in our simulations, in which collisions tended to take place when rogue objects were at periastron and their velocities were roughly parallel to the planets' orbits. They also found a separate class of impacts from incoming high-eccentricity objects on roughly radial trajectories for which the velocity vectors were roughly perpendicular. In this case, the expected impact speed should be ~$v_K\sqrt{2+e}$. We did not find any examples of such high-velocity collisions among the tens of thousands of impacts in our simulations; the origin of this discrepancy is almost certainly differing assumptions related to the impacting populations, specifically the initial periastron distributions of the impactors. Nonetheless, we did see a broader distribution in the impact speeds on the inner planets (Extended Data Fig. 6). This may be explained by the fact that the inner planets have higher collision probabilities with objects on a wide range of periastron distances, whereas the outer planets are only likely to collide with rogue objects with periastron distances close to their orbital radius (see fig. 2 in Kral et al.¹⁰). The inner planets—especially planets b and c—are the most susceptible to impacts with non-parallel velocity vectors and therefore to a broader distribution of impact speeds.

**Tidal evolution.** One might wonder whether resonances could be reformed by orbital migration driven by tidal interactions between the planets and the central star. To obtain a convergent migration with tides, one needs outward migration. The tides raised by planets g and h in the star can drive an outward migration due to the fact that the rotation period of TRAPPIST-1 is smaller than the orbital period of the two planets, placing them outside the corotation radius.

To evaluate this possibility, we calculated the stellar-tide-induced migration for planets g and h separately using an equilibrium tide formalism⁶⁵, which accounts for the fact that the stellar radius evolves with time (we used the evolution tracks of Baraffe et al.⁶⁶ for a stellar mass of 0.09 $M_{\odot}$). We computed the evolution of $\Delta a$ as a function of time, which is the evolution of the separation between the two planets $a_b(t)-a_g(t)$ with respect to the initial separation $a_g(0)-a_g(0)$. No tidal migration would lead to $\Delta a=0$, convergent migration would lead to $\Delta a<0$ and divergent migration to $\Delta a>0$. Owing to the fact that the radius of a low-mass star decreases significantly during the pre-main sequence, the strength of the stellar tide also decreases. Therefore, depending on the initial time, we can expect different migration timescales. The migration timescale also greatly depends on the assumed stellar dissipation. This is why we treated stellar dissipation as a free parameter and explored a range of [10⁻³, ..., 10³] times the dissipation of a fully convective star⁶⁷,⁶⁸.
Note that owing to the lack of a radiative core for low-mass stars, we do not expect the dynamical tide to play as important a role as for Sun-like stars⁶⁰. Looking at a very wide range of dissipation is therefore a way to encompass the potential additional effect of inertial waves propagating in the star.

The results for an initial time of 3 Myr can been seen in Extended Data Fig. 7. Owing to the decrease of the stellar radius, most of the convergent migration occurs in the first million years. After that, the evolution timescales become too long and the evolution freezes. Therefore, if the convergent migration was not enough to bring the planets back into resonance after a few million years, it never will. We find that unless the star has a dissipation higher than 80 times the fiducial value, the convergent migration does not exceed $10^{-4}$au and is thus insufficient to bring the planets back into resonance.

Our reasoning assumes that the inner planets do not play a role in the outward migration of planets g and h; however, in reality all planets will be submitted to the stellar-tide-induced migration. Therefore, for completeness, we performed full N-body simulations including tides (using the Posidonius code⁷⁶,⁷⁷) for different values of the stellar and planetary dissipation. We found that for stellar dissipations as low as the fiducial value, the system is destabilized. In particular, as planet b is closer than the corotation radius, it collides with the star before the star has time to spin down. This means that the dissipation needed to drive a sufficient convergent migration of planets g and h is not compatible with the system as we observe it today.

The star's spin rate—and therefore the location of the corotation radius—at early times is highly uncertain²². If the star span fast enough during its early evolution, the corotation radius would have been very close-in and tidally induced migration could have been directed outwards for all of the planets (see Bolmont et al.⁷³ for a discussion of the effect of the stellar spin history on tidal migration). Although the calculations presented here do not account for that possibility, it is something to be addressed in future detailed modelling of the tidal evolution of the system.

We have also addressed here that only star–planet tides influence the orbits of the planets, meaning that we have neglected other processes, such as planet–planet tides. However, it has been shown by Hay and Matsuyama⁴ that the migration planet–planet tides induce is negligible compared with star–planet tides. Moreover, planet–planet tides would induce an inward migration, which would lead to a spreading of the orbits, rather than the convergent migration we are investigating, making the convergent migration of planets g and h even more challenging.

**Internal structure modelling to estimate the planets' water contents.** We estimated the water mass fraction of each of the TRAPPIST-1 planets as constrained by their measured sizes¹,⁴¹ and masses¹,⁷⁵. Each of the models (i)–(v) assumes uniform rocky interiors among all seven planets. This is because within a system, planets are expected to be similar in relative refractory abundances²⁶,⁷⁶,⁷⁷.

The assumptions made in the five models were:

* In model (i) we assumed a rocky interior with a core mass fraction of 23.5% that best fits the measured densities of planets b, c and d.

* In model (ii) we assumed rocky interiors with a core mass fraction of 26.6% that best fitted the measured density of the planet with the highest (uncompressed) density, planet c.

* In model (iii) we assumed bulk molar Fe/Mg = 0.75 and Mg/Si = 1.02 as suggested by Unterborn et al.²⁶ to represent a proxy for the rocky interiors of the TRAPPIST-1 planets.

* In model (iv) we constrained rocky interiors to solar ratios of Fe/Mg = 0.832 and Mg/Si = 0.886 (ref. ⁷⁸).

* In model (v) we use identical abundance constraints as in model (iii) but assume a core-free interior with all iron oxidized in the mantle.

As suggested by Dorn et al.²⁵, the densest planet may be a good proxy for the rocky components of other planets within the same system. This is the foundation of model (ii), for which the rocky interior is constrained by planet c's density (for a core mass fraction of 26.6%). Model (i) is an adoption of model (ii) but assumed that the densities of the three inner, least-volatile-rich planets are representative of the interiors of all seven planets (for a core mass fraction 23.5%). Atmosphere models²⁹ suggest that the three inner planets may be depleted in water and other volatiles¹ (see also Fig. 3), which justifies using the density of these planets as a proxy for the interior composition of the outer planets. Models (iii) and (v) were both constrained by the stellar proxy³⁶, which is based on refractory abundances of F–G–K stars with metallicities similar to TRAPPIST-1. Whereas model (iii) assumed all of each planet's iron is found in the core, model (v) had all the planets' iron oxidized in their mantles with core-free interiors, in line with recent estimates of oxygen budgets of extrasolar rocks²⁹. Model (iv) was for reference and assumed solar relative refractory abundances.

Models (i)–(v) were then used to calculate the possible water mass fractions needed to fit the planetary mass and radius data, as well as stellar irradiation. To accomplish this, we calculated the water mass fractions for all planetary mass–radius pairs of the $10^4$ posterior samples as in Agol et al.¹, thereby accounting for non-Gaussian-distributed data errors. Whenever the sample's density was equal to or higher than the pure rocky interior, we set the water mass fraction of the sample to zero.

Interior profiles were calculated in a similar way to previous studies¹,⁴³. Specifically, we used equations of state for pure iron⁴⁰, mantle silicates²⁷, solid or liquid water⁴¹ and steam water⁴². Following Agol et al.¹, water was treated as a condensed layer (liquid and/or solid) for planets less irradiated than the runaway greenhouse insolation threshold (that is, TRAPPIST-1 e, f, g and h) and as steam for planets more irradiated than the runaway greenhouse insolation threshold (that is, TRAPPIST-1 b, c and d)²⁴,⁸³. The atmospheric structure and radiative balance were computed using the LMD Generic 1D climate model²¹,⁸³.

We neglected hydrated mantle rocks as they only affect the planet radius by less than 2% at maximum⁴⁴, such that water is assumed to be present only in ice layers, oceans and/or an atmosphere. Also, we do not show interiors with light core alloys, as our tested interiors with 70% FeS in the core only led to marginal changes in planetary radii (<1%) compared with pure Fe cores while keeping the same bulk Fe/Si ratios.

Tables containing all of the results of our modelling are included as Supplementary Tables 2–6.

**Description of the simulation that formed our illustrative TRAPPIST-1 analogue.** The TRAPPIST-1 analogue planets in Fig. 2 represent a super-Earth system produced in simulations of Izidoro et al.¹. Their planet formation model includes the effects of pebble accretion, type-I migration, eccentricity and inclination damping due to tidal interactions of protoplanetary objects with the disk. This particular system was produced in the set-up of Model-III of Izidoro et al. It started with ~70 $M_{\text{moon}}$ planetary seeds in the gas disk, which were initially placed on orbital radii between 0.2 au and 2 au and radially separated from each other by 0.025 au. The inner edge of the gas disk was located at 0.1 au. At the beginning of the simulation, the snow line was beyond 2 au. Seeds first grew via pebble accretion (interior to the snow line), then type-I migrated and eventually underwent mutual collisions. As the disk evolved, the snow line moved inwards, eventually sweeping interior to the orbits of the outermost seeds and allowing them to accrete water-ice-rich pebbles. The outermost seeds grew fast enough to regulate the pebble flux to the inner ones. During the 5 Myr lifetime of the gaseous disk, a system of seven planets formed in a multiresonant configuration between 0.05 au and 0.24 au. The planets' masses were 0.77 $M_{\text{⊕}}$, 0.92 $M_{\text{⊕}}$, 1.03 $M_{\text{⊕}}$, 1.42 $M_{\text{⊕}}$, 1.09 $M_{\text{⊕}}$, 0.49 $M_{\text{⊕}}$ and 0.42 $M_{\text{⊕}}$. The period ratio of adjacent planet pairs were approximately 1.25, 1.33, 1.33, 1.33, 1.25 and 2.42 (the two outer planets show libration of certain angles associated with the 12:5 and 5:2 resonances). The four innermost planets formed dry, meaning that they did not accrete water-rich pebbles or collide with any water-carrying seeds during the entire course of the simulation. The fifth planet contained ~14% water by mass, the sixth planet was dry and the outermost planet contained ~20% water by mass. The disk dissipated after 5 Myr and we continued the simulation in a gas-free environment for a further 50 Myr to verify that the resonant configuration remained dynamically stable.

The primary goal of Fig. 2 is to illustrate one possible pathway for the rapid growth of TRAPPIST-1-like systems. This particular system is clearly not a perfect match to the TRAPPIST-1 system. First, the resonant structure is not quite right, as the system contains too many 5:4 and 4:3 mean motion resonances, rather than 3:2 (see Papalozou et al.¹⁹, Coleman et al.²⁰ and Lin et al.⁵⁹ for careful studies of the resonant structure of simulated TRAPPIST-1 analogues). In addition, the two planets with >10% water by mass are at the (upper) edge of the allowed parameter space in the interior structure models presented in Fig. 3. Nonetheless, our chosen TRAPPIST-1 analogue does have close similarities to the TRAPPIST-1 system, such as the high planet multiplicity and resonant state of adjacent planets. The simulation also shows that planet formation models can qualitatively account for the different bulk compositions of the TRAPPIST-1 planets as predicted by interior models.

## Data availability

All simulation data and model data that support the findings of this study or were used to make the plots are available from the corresponding author upon reasonable request.

## Code availability

The Mercury N-body integrator¹² is publicly available and can be downloaded at https://github.com/4xxi/mercury. We made trivial modifications to stop simulations in which two planets collided and to create an additional data dump any time a planetesimal was removed from the simulation. For water content calculations, the interior structure code is available from C.D. (cdorn@physik.uzh.ch) upon reasonable request, and the atmospheric code is based on polynomial fits from Turbet et al.²¹.

Received: 18 March 2021; Accepted: 17 September 2021; Published online: 25 November 2021

## References

1. Gillon, M. et al. Seven temperate terrestrial planets around the nearby ultracool dwarf star TRAPPIST-1. *Nature* **542**, 456–460 (2017).

2. Luger, R. et al. A seven-planet resonant chain in TRAPPIST-1. *Nat. Astron.* **1**, 0129 (2017).
3. Agol, E. et al. Refining the transit-timing and photometric analysis of TRAPPIST-1: masses, radii, densities, dynamics, and ephemerides. *Planet. Sci. J.* **2**, 1 (2021).

4. Day, J. M. D., Pearson, D. G. & Taylor, L. A. Highly siderophile element constraints on accretion and differentiation of the Earth-Moon system. *Science* **315**, 217–219 (2007).

5. Walker, R. J. Highly siderophile elements in the Earth, Moon and Mars: update and implications for planetary accretion and differentiation. *Geochemistry* **69**, 101–125 (2009).

6. Wood, B. J. & Halliday, A. N. Cooling of the Earth and core formation after the giant impact. *Nature* **437**, 1345–1348 (2005).

7. Kleine, T. et al. Hf-W chronology of the accretion and early evolution of asteroids and terrestrial planets. *Geochim. Cosmochim. Acta* **73**, 5150–5188 (2009).

8. Fernandez, J. A. & Ip, W. Some dynamical aspects of the accretion of Uranus and Neptune: the exchange of orbital angular momentum with planetesimals. *Icarus* **58**, 109–120 (1984).

9. Izidoro, A. et al. Formation of planetary systems by pebble accretion and migration. Hot super-Earth systems from breaking compact resonant chains. *Astron. Astrophys.* **650**, A152 (2021).

10. Kral, Q. et al. Cometary impacts on the TRAPPIST-1 planets can destroy all planetary atmospheres and rebuild secondary atmospheres on planets f, g, and h. *Mon. Not. R. Astron. Soc.* **479**, 2649–2672 (2018).

11. Došović, V., Novaković, B., Vukotić, B. & Ćirković, M. M. Water transport throughout the TRAPPIST-1 system: the role of planetesimals. *Mon. Not. R. Astron. Soc.* **499**, 4626–4637 (2020).

12. Chambers, J. E. A hybrid symplectic integrator that permits close encounters between massive bodies. *Mon. Not. R. Astron. Soc.* **304**, 793–799 (1999).

13. Gonzales, E. C. et al. A reanalysis of the fundamental parameters and age of TRAPPIST-1. *Astrophys. J.* **886**, 131 (2019).

14. Tamayo, D., Rein, H., Petrovich, C. & Murray, N. Convergent migration renders TRAPPIST-1 long-lived. *Astrophys. J. Lett.* **840**, L19 (2017).

15. Morbidelli, A., Lunine, J. I., O'Brien, D. P., Raymond, S. N. & Walsh, K. J. Building terrestrial planets. *Annu. Rev. Earth Planet. Sci.* **40**, 251–275 (2012).

16. Dauphas, N. & Pourmand, A. Hf-W-Th evidence for rapid growth of Mars and its status as a planetary embryo. *Nature* **473**, 489–492 (2011).

17. Marchi, S., Canup, R. M. & Walker, R. J. Heterogeneous delivery of silicate and metal to the Earth by large planetesimals. *Nat. Geosci.* **11**, 77–81 (2018).

18. Ormel, C. W., Liu, B. & Schoonenberg, D. Formation of TRAPPIST-1 and other compact systems. *Astron. Astrophys.* **604**, A1 (2017).

19. Papalozou, J. C. B., Szuszkiewicz, E. & Terquem, C. The TRAPPIST-1 system: orbital evolution, tidal dissipation, formation and habitability. *Mon. Not. R. Astron. Soc.* **476**, 5032–5056 (2018).

20. Coleman, G. A. L., Leleu, A., Alibert, Y. & Benz, W. Pebbles versus planetesimals: the case of Trappist-1. *Astron. Astrophys.* **631**, A7 (2019).

21. Miguel, Y., Cridland, A., Ormel, C. W., Fortney, J. J. & Ida, S. Diverse outcomes of planet formation and composition around low-mass stars and brown dwarfs. *Mon. Not. R. Astron. Soc.* **491**, 1998–2009 (2020).

22. Terquem, C. & Papalozou, J. C. B. Migration and the formation of systems of hot super-Earths and Neptunes. *Astrophys. J.* **654**, 1110–1120 (2007).

23. Dorn, C. et al. A generalized Bayesian inference method for constraining the interiors of super Earths and sub-Neptunes. *Astron. Astrophys.* **597**, A37 (2017).

24. Turbet, M. et al. Revised mass-radius relationships for water-rich rocky planets more irradiated than the runaway greenhouse limit. *Astron. Astrophys.* **638**, A41 (2020).

25. Dorn, C., Mosegaard, K., Grimm, S. L. & Alibert, Y. Interior characterization in multiplanetary systems: TRAPPIST-1. *Astrophys. J.* **865**, 20 (2018).

26. Unterborn, C. T., Desch, S. J., Hinkel, N. R. & Lorenzo, A. Inward migration of the TRAPPIST-1 planets as inferred from their water-rich compositions. *Nat. Astron.* **2**, 297–302 (2018).

27. Barth, P. et al. Magma ocean evolution of the TRAPPIST-1 planets. *Astrobiology* https://doi.org/10.1089/ast.2020.2277 (2021).

28. Acuña, L. et al. Characterisation of the hydrospheres of TRAPPIST-1 planets. *Astron. Astrophys.* **647**, A53 (2021).

29. Alexander, C. M. O., McKeegan, K. D. & Altwegg, K. Water reservoirs in small planetary bodies: meteorites, asteroids, and comets. *Space Sci. Rev.* **214**, 36 (2018).

30. Schlichting, H. E., Sari, R. & Yalinewich, A. Atmospheric mass loss during planet formation: the importance of planetesimal impacts. *Icarus* **247**, 81–94 (2015).

31. Marty, B. The origins and concentrations of water, carbon, nitrogen and noble gases on Earth. *Earth Planet. Sci. Lett.* **313**, 56–66 (2012).

32. Mulders, G. D., Ciesla, F. J., Min, M. & Pascucci, I. The snow line in viscous disks around low-mass stars: implications for water delivery to terrestrial planets in the habitable zone. *Astrophys. J.* **807**, 9 (2015).

33. Ramirez, R. M. & Kaltengegger, L. The habitable zones of pre-main-sequence stars. *Astrophys. J.* **797**, L25 (2014).

34. Tian, F. & Ida, S. Water contents of Earth-mass planets around M dwarfs. *Nat. Geosci.* **8**, 177–180 (2015).

35. Luger, R. & Barnes, R. Extreme water loss and abiotic O₂ buildup on planets throughout the habitable zones of M dwarfs. *Astrobiology* **15**, 119–143 (2015).

36. Bolmont, E. et al. Water loss from terrestrial planets orbiting ultracool dwarfs: implications for the planets of TRAPPIST-1. *Mon. Not. R. Astron. Soc.* **464**, 3728–3741 (2017).

37. Kislyakova, K. G. et al. Magma oceans and enhanced volcanism on TRAPPIST-1 planets due to induction heating. *Nat. Astron.* **1**, 878–885 (2017).

38. Leleu, A. et al. Six transiting planets and a chain of Laplace resonances in TOI-178. *Astron. Astrophys.* **649**, A26 (2021).

39. Mills, S. M. et al. A resonant chain of four transiting, sub-Neptune planets. *Nature* **533**, 509–512 (2016).

40. IAU Resolution B5: Definition of a Planet in the Solar System (IAU, 2006).

41. Gillon, M. et al. Temperate Earth-sized planets transiting a nearby ultracool dwarf star. *Nature* **533**, 221–224 (2016).

42. Burgasser, A. J. & Mamajek, E. E. On the age of the TRAPPIST-1 system. *Astrophys. J.* **845**, 110 (2017).

43. Mann, A. W. et al. How to constrain your M dwarf. II. The mass-luminosity-metallicity relation from 0.075 to 0.70 solar masses. *Astrophys. J.* **871**, 63 (2019).

44. Ducrot, E. et al. TRAPPIST-1: global results of the Spitzer Exploration Science Program Red Worlds. *Astron. Astrophys.* **640**, A112 (2020).

45. Lee, M. H. & Peale, S. J. Dynamics and origin of the 2:1 orbital resonances of the GJ 876 planets. *Astrophys. J.* **567**, 596–609 (2002).

46. Cresswell, P., Dirksen, G., Kley, W. & Nelson, R. P. On the evolution of eccentric and inclined protoplanets embedded in protoplanetary disks. *Astron. Astrophys.* **473**, 329–342 (2007).

47. Izidoro, A. et al. Breaking the chains: hot super-Earth systems from migration and disruption of compact resonant chains. *Mon. Not. R. Astron. Soc.* **470**, 1750–1770 (2017).

48. Bottke, W. F., Walker, R. J., Day, J. M. D., Nesvorny, D. & Elkins-Tanton, L. Stochastic late accretion to Earth, the Moon, and Mars. *Science* **330**, 1527–1530 (2010).

49. Raymond, S. N., Schlichting, H. E., Hersant, F. & Selsis, F. Dynamical and collisional constraints on a stochastic late veneer on the terrestrial planets. *Icarus* **226**, 671–681 (2013).

50. Morbidelli, A. & Wood, B. J. in *The Early Earth: Accretion and Differentiation* (eds Badro, J. & Walter, M.) 71–82 (Geophysical Monograph Series, Vol. 212, John Wiley & Sons, 2015).

51. Levison, H. F. & Duncan, M. J. The long-term dynamical behavior of short-period comets. *Icarus* **108**, 18–36 (1994).

52. Boss, A. P. et al. Astrometric constraints on the masses of long-period gas giant planets in the TRAPPIST-1 planetary system. *Astron. J.* **154**, 103 (2017).

53. Jontof-Hutter, D., Truong, V. H., Ford, E. B., Robertson, P. & Terrien, R. C. Dynamical constraints on nontransiting planets orbiting TRAPPIST-1. *Astron. J.* **155**, 239 (2018).

54. Ogihara, M. & Ida, S. N-body simulations of planetary accretion around M dwarf stars. *Astrophys. J.* **699**, 824–838 (2009).

55. McNeil, D. S. & Nelson, R. P. On the formation of hot Neptunes and super-Earths. *Mon. Not. R. Astron. Soc.* **401**, 1691–1708 (2010).

56. Ida, S. & Lin, D. N. C. Toward a deterministic model of planetary formation. VI. Dynamical interaction and coagulation of multiple rocky embryos and super-Earth systems around solar-type stars. *Astrophys. J.* **719**, 810–830 (2010).

57. Cossou, C., Raymond, S. N., Hersant, F. & Pierens, A. Hot super-Earths and giant planet cores from different migration histories. *Astron. Astrophys.* **569**, A56 (2014).

58. Schoonenberg, D., Liu, B., Ormel, C. W. & Dorn, C. Pebble-driven planet formation for TRAPPIST-1 and other compact systems. *Astron. Astrophys.* **627**, A149 (2019).

59. Lin, Y.-C., Matsumoto, Y. & Gu, P.-G. Formation of multiple-planet systems in resonant chains around M dwarfs. *Astrophys. J.* **907**, 81 (2021).

60. Murray, C. D. & Dermott, S. F. *Solar System Dynamics* (Cambridge Univ. Press, 1999).

61. Hadden, S. & Lithwick, Y. A criterion for the onset of chaos in systems of two eccentric planets. *Astron. J.* **156**, 95 (2018).

62. Wyatt, M. C., Bonsor, A., Jackson, A. P., Marino, S. & Shannon, A. How to design a planetary system for different scattering outcomes: giant impact sweet spot, maximizing exocomets, scattered discs. *Mon. Not. R. Astron. Soc.* **464**, 3385–3407 (2017).

63. A'Hearn, M. F. Comets as building blocks. *Annu. Rev. Astron. Astrophys.* **49**, 281–299 (2011).

64. Choukroun, M. et al. Dust-to-gas and refractory-to-ice mass ratios of comet 67P/Churyumov-Gerasimenko from Rosetta observations. *Space Sci. Rev.* **216**, 44 (2020).

65. Bolmont, E., Raymond, S. N. & Leconte, J. Tidal evolution of planets around brown dwarfs. *Astron. Astrophys.* **535**, A94 (2011).
66. Baraffe, I., Homeier, D., Allard, F. & Chabrier, G. New evolutionary models for pre-main sequence and main sequence low-mass stars down to the hydrogen-burning limit. *Astron. Astrophys.* **577**, A42 (2015).

67. Hansen, B. M. S. Calibration of equilibrium tide theory for extrasolar planet systems. *Astrophys. J.* **723**, 285–299 (2010).

68. Bolmont, E., Raymond, S. N., Leconte, J. & Matt, S. P. Effect of the stellar spin history on the tidal evolution of close-in planets. *Astron. Astrophys.* **544**, A124 (2012).

69. Bolmont, E. & Mathis, S. Effect of the rotation and tidal dissipation history of stars on the evolution of close-in planets. *Celest. Mech. Dynam. Astron.* **126**, 275–296 (2016).

70. Blanco-Cuaresma, S. & Bolmont, E. Studying tidal effects in planetary systems with Posidonius. A N-body simulator written in Rust. Zenodo https://doi.org/10.5281/zenodo.1095095 (2017).

71. Bolmont, E. et al. Impact of tides on the transit-timing fits to the TRAPPIST-1 system. *Astron. Astrophys.* **635**, A117 (2020).

72. Herbst, W., Eislöffel, J., Mundt, R. & Scholz, A. in *Protostars and Planets V* (eds Reipurth, B. et al.) 297–311 (Univ. Arizona Press, 2007).

73. Bolmont, E., Raymond, S. N., Leconte, J. & Matt, S. P. Effect of the stellar spin history on the tidal evolution of close-in planets. *Astron. Astrophys.* **544**, A124 (2012).

74. Hay, H. C. F. C. & Matsuyama, I. Tides between the TRAPPIST-1 planets. *Astrophys. J.* **875**, 22 (2019).

75. Grimm, S. L. et al. The nature of the TRAPPIST-1 exoplanets. *Astron. Astrophys.* **613**, A68 (2018).

76. Dorn, C. et al. Can we constrain the interior structure of rocky exoplanets from mass and radius measurements? *Astron. Astrophys.* **577**, A83 (2015).

77. Sotin, C., Grasset, O. & Mocquet, A. Mass-radius curve for extrasolar Earth-like planets and ocean planets. *Icarus* **191**, 337–351 (2007).

78. Lodders, K., Palme, H. & Gail, H.-P. Abundances of the elements in the solar system. Preprint at https://arxiv.org/abs/0901.1149 (2009).

79. Doyle, A. E., Young, E. D., Klein, B., Zuckerman, B. & Schlichting, H. E. Oxygen fugacities of extrasolar rocks: evidence for an Earth-like geochemistry of exoplanets. *Science* **366**, 356–359 (2019).

80. Hakim, K. et al. A new ab initio equation of state of hcp-Fe and its implication for the interior structure and mass-radius relations of rocky super-Earths. *Icarus* **313**, 61–78 (2018).

81. Haldemann, J., Alibert, Y., Mordasini, C. & Benz, W. Aqua: a collection of H₂O equations of state for planetary models. *Astron. Astrophys.* **643**, A105 (2020).

82. Haar, L., Gallagher, J. & Kell, G. NBS/NRC Steam Tables Thermodynamic and Transport Properties and Computer Programs for Vapor and Liquid States of Water in SI Units (Hemisphere, 1984).

83. Turbet, M., Ehrenreich, D., Lovis, C., Bolmont, E. & Fauchez, T. The runaway greenhouse radius inflation effect. An observational diagnostic to probe water on Earth-sized planets and test the habitable zone concept. *Astron. Astrophys.* **628**, A12 (2019).

84. Shah, O., Alibert, Y., Helled, R. & Mezger, K. Internal water storage capacity of terrestrial planets and the effect of hydration on the M-R relation. Preprint at https://arxiv.org/abs/2012.06455 (2020).

85. Mamajek, E. E. Initial conditions of planet formation: lifetimes of primordial disks. *AIP Conf. Proc.* **1158**, 3–10 (2009).

86. Turbet, M. et al. Modeling climate diversity, tidal dynamics and the fate of volatiles on TRAPPIST-1 planets. *Astron. Astrophys.* **612**, A86 (2018).

## Acknowledgements

We thank S. Rousseau and A. Collioud for managing the machines on which our simulations were run. S.N.R. is grateful to the ECLIPSE research team at the Laboratoire d'Astrophysique de Bordeaux and the 'Zooming In On Rocky Planet Formation' ISSI Team for helpful discussions, and to A. Guilbert-Lepoutre and A. Gkotsinas for input on the water contents of comets. S.N.R. and F.S. thank the CNRS's PNP programme for support. A.I. and R.D. acknowledge NASA grant number 80NSSC18K0828 for financial support during the preparation and submission of this work. A.I. also thanks the Brazilian Federal Agency for Support and Evaluation of Graduate Education (CAPES), in the scope of the Program CAPES-PrInt, process number 88887.310463/2018-00, International Cooperation Project number 3266 and CNPq (313998/2018-3). P.B. acknowledges a St Leonard's Interdisciplinary Doctoral Scholarship from the University of St Andrews. L.C. acknowledges support from the DFG Priority Programme SP1833 grant number CA 1795/3. E.A. and S.N.R. acknowledge support from the Virtual Planetary Laboratory Team, a member of the NASA Nexus for Exoplanet System Science, funded via NASA Astrobiology Program grant number 80NSSC18K0829. This project has received funding from the European Union's Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement number 832738/ESCAPE. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. M.T. and E.B. acknowledge the financial support of the SNSF. This work was performed using the High-Performance Computing (HPC) resources of Centre Informatique National de l'Enseignement Supérieur (CINES) under allocation number A0080110391 made by Grand Équipement National de Calcul Intensif (GENCI). M.T. thanks the Gruber Foundation for its support to this research. C.D. acknowledges support from the Swiss National Science Foundation under grant number PZ00P2_174028.

## Author contributions

S.N.R. initiated the project after discussions with L.C., P.B. and A.I., and ran and analysed the simulations using codes optimized by A.I. and configurations of the system provided by S.L.G. and E.A. C.D. and M.T. performed the interior modelling to obtain constraints on the water contents of the TRAPPIST-1 planets. E.B. performed the analysis of the tidal effects. S.N.R. wrote the paper with feedback from A.I., E.B., C.D., F.S., M.T., E.A., P.B., L.C., R.D., M.G. and S.L.G.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at https://doi.org/10.1038/s41550-021-01518-6.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-021-01518-6.

**Correspondence and requests for materials** should be addressed to Sean N. Raymond.

**Peer review information** *Nature Astronomy* thanks Shigeru Ida and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2021
<figure><img src="image_6.png" /><figcaption>Trappist-1 system, no extra objects</figcaption></figure>

**Extended Data Fig. 1** | Evolution of resonant angles for our fiducial Trappist-1 system with no added objects. The top four panels each show one angle that characterizes the mean motion resonance between a pair of neighboring planets, as defined in Equations 1 and 2. The bottom three panels show angles that characterize the 3-body Laplace resonances between triplets of planets, as defined in Eq. 3.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 2 | Evolution of the resonant angles for three simulations that included our fiducial Trappist-1 system and 1000 rogue planetesimals.**

Each system was stable for the full 10 Myr simulation. The resonant angles that are shown are the same ones presented in Extended Data Fig. 1. The total mass in planetesimals for each system was 0.0213 $M_⊕$ (left panel; resonance maintained), 0.0429 $M_⊕$ (middle panel; resonance barely maintained), and 0.0645 $M_⊕$ (right panel; resonance disrupted).
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 3** | Disruption of the resonances in the outer parts of the Trappist-1 system. In this simulation, the fiducial system of planets interacted with 1000 planetesimals containing a total of 0.0967 M<sub>⊕</sub>. The top panel shows the increase in the semimajor axis of the orbits of planets f, g and h. The other panels show the angles associated with mean-motion and Laplace resonances between the planets.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 4** | Effect of the semimajor axis distribution of rogue planetesimals. The top left panel shows our the subset of our fiducial batch of simulations with 1000 rogue planetesimals with a total mass in rogue planetesimals between 0.01 and 1 Earth mass. The top right and bottom left panels show the sets of simulations with $a_{rogue} = 0.2 - 0.3$ au and $a_{rogue} = 0.7 - 1$ au, and the same periastron distribution as in the fiducial set (with orbits initially crossing those of only the outer planets). The bottom right panel shows a set of simulations in which the rogue planetesimals were on extremely eccentric orbits, with semimajor axes $a_{rogue} = 0.7 - 1$ au and periastron distance $q_{rogue}$ randomly chosen to sample a logarithmic distribution in the range 0.005 – 0.06 au. As in Fig. 1 in the main text, blue filled circles indicate systems that remained in resonant chains, red diamonds those that remained stable for 10 Myr but lost their resonant structure, and black empty circles those that underwent dynamical instabilities.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

Extended Data Fig. 5 | See next page for caption.
**Extended Data Fig. 5 | Effect of the assumed configuration of the Trappist-1 system.** The left-hand panels show the outcomes of simulations containing a single rogue planet, for four different ‘best-fit’ configurations of the Trappist-1 system that are consistent with current observations³. Likewise, the right-hand panels show the outcomes of simulations containing 1000 rogue planetesimals. As in Fig. 1 in the main text, systems that were stable for 10 Myr but in which resonances were broken (red diamonds) are shifted downward slightly for clarity. The top panels (‘Set 1’) represent our fiducial simulations (note that we have only used the subset of fiducial simulations with 1000 rogue planetesimals that fall in the identical mass range as the simulations for the other sets).
<figure><img src="image_23.png" /></figure>

**Extended Data Fig. 6** | Effect of star-planet tidal interactions on the outer planets in the Trappist-1 system. Left panel: Evolution of the stellar radius with time (given as $t - t_{\text{init}}$, where $t_{\text{init}} = 3$ Myr). Right panel: map showing the values of $\log(\Delta|a|)$ as a function of the stellar dissipation and time. The stellar dissipation is given with respect to a fiducial value corresponding to a fully convective star or brown dwarf.
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

**Extended Data Fig. 7** | Distribution of impact speeds of planetesimals in simulations in which the planets remained in resonance on each of the Trappist-1 planets.
<table><caption>**Extended Data Table 1** | Upper limits on the total mass in planetesimals that could have collided with each of the Trappist-1 planets (assuming perfect mergers), for four sets of initial orbital distributions for the rogue planetesimals. The first three columns share the same periastron distribution, with starting semimajor axes of 0.07-0.1 au (our fiducial set), 0.2-0.3 au, and 0.7-1 au. The fourth column is from the set in which rogue planetesimals have starting semimajor axes of 0.7-1 au and super-eccentric orbits, with starting periastron distances of 0.005-0.063 au. For each set we used the maximum total initial planetesimal mass that preserved the resonant structure of the system: 0.05 M<sub>⊕</sub> for a<sub>rogue</sub> = 0.07-0.1 au and a<sub>rogue</sub> = 0.2-0.3 au, 0.1 M<sub>⊕</sub> for a<sub>rogue</sub> = 0.7-1 au with the standard periastron distribution, and 0.018 M<sub>⊕</sub> for a<sub>rogue</sub> = 0.7-1 au with super-eccentric orbits (see Extended Data Fig. 4.)</caption><thead><tr><th>Planet</th><th>Orbital radius<br>(AU)</th><th>Max. bombard.<br>mass (M<sub>⊕</sub>)<br>[0.07-0.1 au]</th><th>Max. bombard.<br>mass (M<sub>⊕</sub>)<br>[0.2-0.3 au]</th><th>Max. bombard.<br>mass (M<sub>⊕</sub>)<br>[0.7-1 au]</th><th>Max. bombard.<br>mass (M<sub>⊕</sub>)<br>[0.7-1 au]<br>[super-ecc.]</th></tr></thead><tbody><tr><td>b</td><td>0.0115</td><td>0.00035</td><td>0.00013</td><td>0.0001</td><td>0.0082</td></tr><tr><td>c</td><td>0.0158</td><td>0.0015</td><td>0.0005</td><td>0.00023</td><td>0.0073</td></tr><tr><td>d</td><td>0.0223</td><td>0.0015</td><td>0.00047</td><td>0.0003</td><td>0.0033</td></tr><tr><td>e</td><td>0.0293</td><td>0.0034</td><td>0.00098</td><td>0.00087</td><td>0.0044</td></tr><tr><td>f</td><td>0.0385</td><td>0.0078</td><td>0.0038</td><td>0.003</td><td>0.0063</td></tr><tr><td>g</td><td>0.0469</td><td>0.017</td><td>0.01</td><td>0.0083</td><td>0.0074</td></tr><tr><td>h</td><td>0.0620</td><td>0.014</td><td>0.013</td><td>0.011</td><td>0.0045</td></tr></tbody></table>
<table><caption>**Extended Data Table 2** | Starting orbital configuration of the Trappist-1 system that were used in our fiducial simulations. These are referred to as 'set 1' in the Methods. The planets' orbits have zero inclinations, but in our simulations rogue planets and planetesimals were given non-zero inclinations to ensure that the systems were not 2D. The system's exact configuration was drawn from the best-fit posteriors of Agol et al.³</caption><thead><tr><th>Planet</th><th>Mass (M<sub>⊕</sub>)</th><th>Radius (R<sub>⊕</sub>)</th><th>Semimajor Axis <i>a</i> (AU)</th><th>Eccentricity <i>e</i></th><th>Longitude of periastron ϖ (°)</th><th>Mean Anomaly <i>M</i> (°)</th></tr></thead><tbody><tr><td>b</td><td>1.3925</td><td>1.1174</td><td>0.011551</td><td>0.002344</td><td>253.61247</td><td>105.78489</td></tr><tr><td>c</td><td>1.2943</td><td>1.0967</td><td>0.015820</td><td>0.001224</td><td>132.62793</td><td>54.89836</td></tr><tr><td>d</td><td>0.3958</td><td>0.7880</td><td>0.02229</td><td>0.005045</td><td>202.45580</td><td>171.39157</td></tr><tr><td>e</td><td>0.6824</td><td>0.9200</td><td>0.02930</td><td>0.007013</td><td>52.42997</td><td>30.97582</td></tr><tr><td>f</td><td>1.0634</td><td>1.0448</td><td>0.038551</td><td>0.008298</td><td>170.04247</td><td>247.44087</td></tr><tr><td>g</td><td>1.3464</td><td>1.1294</td><td>0.046896</td><td>0.003760</td><td>355.97714</td><td>87.27858</td></tr><tr><td>h</td><td>0.3198</td><td>0.7552</td><td>0.061963</td><td>0.003571</td><td>172.18673</td><td>118.58431</td></tr></tbody></table>