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

Sean N. Raymond1 Andre Izidoro2 Emeline Bolmont3 Caroline Dorn4 Franck Selsis1
Martin Turbet3 Eric Agol5 Patrick Barth6,7,8 Ludmila Carone9 Rajdeep Dasgupta2
Michael Gillon10 Simon L. Grimm11

###### Abstract

The Trappist-1 system contains seven roughly Earth-sized planets locked in a
multi-resonant orbital configuration,[1, 2] which has enabled precise
measurements of the planets’ masses and constrained their compositions.[3]
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 multi-resonant
configuration by either triggering dynamical instability or simply removing
the planets from resonance. The planets cannot have interacted with more than
$\mathbf{\sim 5\%}$ of an Earth mass ($\mathbf{M_{\oplus}}$) 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
$\mathbf{10^{-4}}$ to $\mathbf{10^{-2}M_{\oplus}}$ 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[4, 5],
and demonstrates that the Trappist-1 planets’ growth was complete in just a
few million years, roughly an order of magnitude faster than Earth’s.[6, 7]
Our results imply that any large water reservoirs on the Trappist-1 planets
must have been incorporated during their formation in the gaseous disk.

Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée
Geoffroy Saint-Hilaire, 33615 Pessac, France

Department of Earth, Environmental and Planetary Sciences, Rice University,
6100 Main Street, MS 126, Houston, TX 77005, USA

Observatoire de Genève, Université de Genève, Chemin Pegasi 51, 1290,
Sauverny, Switzerland

University of Zurich, Institute of Computational Sciences, Winterthurerstrasse
190, CH-8057, Zurich, Switzerland

Department of Astronomy, University of Washington, Seattle, WA 98195-1580, USA

Centre for Exoplanet Science, University of St Andrews, North Haugh, St
Andrews, KY16 9SS, UK

SUPA, School of Physics & Astronomy, University of St Andrews, North Haugh, St
Andrews, KY16 9SS, UK

School of Earth & Environmental Sciences, University of St Andrews, Irvine
Building, St Andrews, KY16 9AL, UK

Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

Space sciences, Technologies and Astrophysics Research (STAR) Institute,
Universite de Liege, 4031 Liege, Belgium

Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6,
3012 Bern, Switzerland

The seven-planet Trappist-1 system is characterized by mean motion resonances
between neighboring planets and 3-body (Laplace) resonances between triplets
of adjacent planets (see Methods sections 1 and 3 and Extended Data Fig.
1).[1, 2] Here we use the multi-resonant 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 Gyr timescales, we invoke gravitational
perturbations from additional objects under two limiting assumptions. First,
we assume that a single rogue planetary embryo from the outer disk was
scattered inward toward 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 multi-resonant
structure and induce instability. We ran 1000 simulations varying the mass of
a rogue planet in the range $0.001-1{\rm\,M_{\oplus}}$. Second, we assume that
a swarm of leftover planetesimals from the outer disk was scattered inward
toward the planets. The planets undergo an exchange of orbital angular
momentum with the planetesimals, which causes their orbits to spread out.[8]
We ran 140 simulations varying the total mass in 1000 equal-mass planetesimals
between $0.01{\rm\,M_{\oplus}}$ and $3{\rm\,M_{\oplus}}$. Planetesimals either
represent leftovers whose orbits were excited by the planets’ migration[9] or
objects scattered inward toward the known planets, perhaps by a more distant
planet.[10, 11] While other potential sources of impacts exist[10], we focused
on inward-scattered rogue objects because they provide analogs to late
accretion on Earth[4, 5] and have the highest impact probabilities[10] (see
Methods Section 4 and Extended Data Table 2). Each system was integrated using
the hybrid Mercury integrator[12] for 10 million years or until it became
unstable (see Methods Section 2, Extended Data Table 1, and Extended Data
Figs. 2 and 3). This is 2-3 orders of magnitude shorter than the lifetime of
the system,[13] meaning that our upper limits are conservative, as some
systems may go 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 (see Methods Sections 4
and 5, Extended Data Tables 1 and 2, Supplementary Table 1, and Extended Data
Figs. 4 and 5).[3] We also calculated that star-planet tidal effects did not
shift the outer planets’ orbits enough to affect our results (see Methods
Section 7 and Extended Data Fig. 7).

The Trappist-1 system is dynamically fragile. While the planets maintained
their multi-resonant structure in most simulations containing a rogue planet
less massive than $\sim 0.005{\rm\,M_{\oplus}}$, more massive rogue planets
were disruptive (Fig. 1). The majority of systems with rogue planet masses
between $0.005$ and $0.02{\rm\,M_{\oplus}}$ disrupted their resonances but
remained stable for the 10 million year 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 a lunar mass
($\sim 0.01{\rm\,M_{\oplus}}$), 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 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.

Figure 1: Response of the Trappist-1 system when interacting with a single
rogue planet (left panel) or 1000 rogue planetesimals (right panel). Blue
circles represent simulations that are consistent with the observed system:
the planets’ orbits remained stable for the 10 Myr integration and retained
their multi-resonant configuration. The red symbols denote systems that were
stable for 10 Myr (shifted slightly downward for visibility) but whose
resonant structure was disrupted, as quantified using Laplace resonant angles
for the outer triplets of planets (see 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 inward from just exterior to the
known planets. The results of simulations with different configurations of the
Trappist-1 planets produced near-identical results (see Methods Section 5).

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 inward. The planets’ orbits widen and
spread out due to the back-reaction of planetesimal scattering.[8] When the
planets scatter a total of more than $\sim 0.05{\rm\,M_{\oplus}}$ in
planetesimals the system’s multi-resonant structure is disrupted, as the
angles that characterize mean-motion and 3-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[2] –
to determine which systems were still in resonance (see Methods Section 3).
After leaving their stabilized resonant state,[14] most systems that
interacted with more than $\sim 0.1{\rm\,M_{\oplus}}$ 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{\rm\,M_{\oplus}}$ – 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,[10] 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,[10, 11] with median velocities of 30 km/s for
planet b decreasing to 13 km/s for planet h (see Methods Section 6 and
Extended Data Fig. 6).

Planet | Orbital radius | Maximum bombardment | Maximum water delivered
---|---|---|---
| (AU) | mass (${\rm\,M_{\oplus}}$) | (Earth oceans)
b | 0.0115 | 0.00038 | 0.15
c | 0.0158 | 0.0015 | 0.64
d | 0.0223 | 0.0016 | 0.68
e | 0.0293 | 0.0035 | 1.54
f | 0.0385 | 0.008 | 3.41
g | 0.0469 | 0.018 | 7.62
h | 0.0620 | 0.012 | 6.16
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 (assuming perfect mergers). Here we have adopted a maximum value
of $0.05{\rm\,M_{\oplus}}$ in scattered rogue planetesimals (see Fig. 1) and
assumed that they contain 10% water by mass, typical of the wettest
carbonaceous chondrite meteorites.[29] One “ocean” of water is equivalent to
the mass of the water on Earth’s surface, $2.3\times
10^{-4}{\rm\,M_{\oplus}}$.

Our results demonstrate that the Trappist-1 planets grew an order of magnitude
faster than Earth (Fig. 2). Cosmochemical constraints[6, 7] and accretion
models[15] indicate that a large portion of Earth’s growth took place after
the dissipation of the Sun’s gaseous disk. Mars’ growth was close to complete
by the time the disk dissipated[16], and Earth’s constituent embryos were
likely $\sim$Mars-sized at that time.[15] Earth’s final, Moon-forming impact
took place $\sim 50-100$ million years later.[7] The abundance of highly
siderophile elements (HSEs) in the mantle and crust suggest that Earth
accreted $\sim 0.003-0.007{\rm\,M_{\oplus}}$ in planetesimals in the 4.45 Gyr
between the Moon-forming impact and the present day;[4, 5] 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 significantly lower than on Earth
(Table 1). The four outer Trappist-1 planets could not have accreted more than
$\sim 1\%$ of their masses during this interval, comparable to upper limits on
Earth’s late accretion.[17] 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.[18, 19, 20] Migration appears to be essential in producing an
orbital configuration with a long dynamical lifetime.[14] The starting orbital
distance of migrating planets is unconstrained, although models expect large
embryos to form fastest just beyond the snow line.[18, 21] While the low
occurrence rate of multi-resonant systems suggest that most resonant chains
undergo dynamical instabilities after the dispersal of the disk,[22, 9]
Trappist-1 represents a rare surviving resonant chain.

Figure 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.[6, 7, 16] After the dispersal of the gaseous
disk, Earth’s growth was likely characterized by giant impacts[15] (including
the Moon-forming impact), which are not reflected in the Hf-W and U-Pb model
curves.[6] We used a simulation by Izidoro et al[9] that formed a close-in,
long-term stable multi-resonant chain of planets. The growth curves of seven
planets with masses between 0.4 and 1.5 ${\rm\,M_{\oplus}}$ – characterized by
an early phase of pebble accretion followed by giant impacts – are included as
analogs for the accretion of the Trappist-1 planets (see Methods Section 9).
Each planet’s growth curve is shown in a shade of blue, with darker colors
representing more distant orbits. While this is only one illustrative example
among many,[20, 21] the Trappist-1 system’s growth must have completed by the
end of the gaseous disk phase, whose duration is constrained by observations
of the occurrence rate of hot dust around stars in clusters of different
ages.[84]

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[23] and atmospheric[24] models constrained by the latest
planetary mass and radius estimates.[3] Given that inferred water contents are
subject to interior degeneracies,[25, 3, 26, 27, 28] we estimated water mass
fractions assuming five different rocky interior models (see 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 through h) are more likely to have
significant water contents (consistent with previous models[28, 3]). Yet while
our mean water contents are $\sim$1-10% wt% for the outer planets, 2-$\sigma$
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\sigma$ level.

Figure 3: Late water delivery to the Trappist-1 planets. The thick and thin
error bars indicate the one and two-sigma confidence intervals (encompassing
68% and 95% probability), respectively, for estimates of the bulk water
content of each planet, for each of the five models discussed in the text
(note that certain models do not apply to the planets b, c and d; see Methods
Section 8). 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,[85] which swept
inward in time as the central star evolved to the main sequence.[33, 34]
Earth’s surface water budget of one “ocean” is shown as the dashed horizontal
line.

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 wettest carbonaceous chondrite
meteorites.[29] 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).[10, 30] Taking the upper limits at face value,
late accretion can deliver enough water for planets e through h to each have
bulk water contents similar to Earth’s value of $\sim 10^{-3}$ by mass.[31]
However, late accretion falls an order of magnitude short of delivering enough
water to explain the mean water contents from our models (Fig. 3). If future
work confirms water contents at the wt% 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.[21, 9, 18, 26, 20] To acquire significant water
on planets e-h requires long-range orbital migration. Some constituent pieces
of each planet must have originated beyond the snow line, which was more
distant during the star’s brighter early phases.[32, 26] Two viable
possibilities are (i) that the cores of the planets originated just past the
snow line and migrated inward[18], or (ii) that the snow line swept inward,
allowing rocky cores to efficiently accrete icy pebbles and migrate inward.[9]

Each of the planets may have lost a significant amount of water after
accretion. High-velocity planetesimal impacts may cause water (and
atmospheric) loss.[30, 10] 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.[10] Given the radial dependence of impact speed
(see Methods Sections 5 and 6), the outer planets are more likely to
accumulate volatiles.[10] In addition, as an ultracool star Trappist-1 evolved
very slowly onto the main sequence; the habitable zone swept inward such that
even the planets in the present-day habitable zone spent time with much higher
equilibrium temperatures.[33, 34] 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).[35, 36] Models of
the early evolution of the Trappist-1 planets that also include the magma
ocean phase[37, 27] find that each planet may have lost a few to 20 Earth
oceans of water.[27] Calculations that make assumptions designed to maximize
water loss still find an upper limit of 25 oceans lost in the interval before
planets e through h reached the habitable zone.[36, 27] This is an order of
magnitude lower than the water contents allowed by observationally-constrained
internal structure modeling (Fig. 3) for each planet, compatible with a
relatively water-rich formation scenario.

Our work highlights the exceptional value of multi-resonant planetary systems,
only a handful of which are currently well-characterized, notably TOI-178[38]
and Kepler-223.[39] Not only are these systems ideally-suited for mass
measurement via transit-timing variation analysis,[3] their dynamical
fragility allows us to strongly constrain their bombardment histories. Our
results demonstrate that the Trappist-1 system satisfies the International
Astronomical Union’s (IAU) criteria for planethood. IAU resolution B5 (see
https://www.iau.org/static/resolutions/Resolution_GA26-5-6.pdf) 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 neighborhood around its orbit. Assuming criteria 1 and 2 to be satisfied
using observations of the central star[1] and internal structure modeling of
the planets[25, 24], our simulations show that the inner Trappist-1 system has
been cleared of rogue material to a stringent limit.

## 1\. 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 Trappist-1 planets’
sizes are between 0.75 and $1.15{\rm\,R_{\oplus}}$, and their masses have been
measured with a precision of 3-6% using transit-timing variation analysis.[3]
The system’s two innermost planets were discovered in 2016[40] and the next
five in 2017.[1] The central star is an ultracool dwarf star of spectral type
M8 located 12 pc from the Sun.[40, 1] A number of analyses have constrained
the star’s physical characteristics.[41, 42, 13, 43] It has a bolometric
luminosity of roughly $5.5\times 10^{-4}L_{\odot}$[43], a mass of $\sim
0.09M_{\odot}$ and a size of $\sim 0.12R_{\odot}$. Its age has proven hard to
pin down. While one study suggest a relatively old age of $\sim 7.6$ Gyr[41],
a recent study[13] found analog 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 multi-resonant 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.[2]
Extended Data Fig. 1 shows the evolution of a number of resonant angles in a
100 million year 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.[44, 45, 46, 9, 20]
Papaloizou et al[19] proposed that each of these planet pairs was in fact
captured into 3:2 or 4:3 resonance; i.e. first-order resonances. This scenario
requires fast circulation of the longitudes of periastron for planets b and c,
which requires damping of their free eccentricities, which naturally happens
due to tides. This remains a plausible scenario to be confirmed with future
studies combining orbital and tidal dynamics.

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

## 2\. Simulations

Our simulations started from a configuration of the Trappist-1 system derived
from the most comprehensive analysis of the system to date.[3] We integrated
the system on its own and verified that it is dynamically stable for at least
1 billion years. Table 2 lists the orbital parameters of our starting system
(as well as those of three other best-fit systems that we tested; see Methods
Section 5). While the system is two-dimensional, we gave non-zero inclinations
to rogue embryos and planetesimals to ensure that the systems as a whole were
3D.

Planet | Mass | Radius | Semimajor | Eccentricity | Longitude of | Mean
---|---|---|---|---|---|---
| (${\rm\,M_{\oplus}}$) | (${\rm\,R_{\oplus}}$) | Axis $a$ (AU) | $e$ | periastron $\varpi$ (∘) | Anomaly $M$ (∘)
b | 1.3925 | 1.1174 | 0.011551 | 0.002344 | 253.61247 | 105.78489
c | 1.2943 | 1.0967 | 0.015820 | 0.001224 | 132.62793 | 54.89836
d | 0.3958 | 0.7880 | 0.02229 | 0.005045 | 202.45580 | 171.39157
e | 0.6824 | 0.9200 | 0.02930 | 0.007013 | 52.42997 | 30.97582
f | 1.0634 | 1.0448 | 0.038551 | 0.008298 | 170.04247 | 247.44087
g | 1.3464 | 1.1294 | 0.046896 | 0.003760 | 355.97714 | 87.27858
h | 0.3198 | 0.7552 | 0.061963 | 0.003571 | 172.18673 | 118.58431
Extended Data Table 1. Starting orbital configuration of the Trappist-1 system
that were used in our fiducial simulations (referred to as ’set 1’ in Methods
Section 5). 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 2-dimensional. The system’s exact configuration was
drawn from the best-fit posteriors of Agol et al.[3]

The initial conditions of our simulations were built on the assumption that
one or more rogue objects was scattered inward toward the Trappist-1 planets.
Rogue objects may be thought of as analogs to the leftover planetesimals that
are thought to have impacted the Earth and other terrestrial planets during
the late accretion phase.[4, 5, 47, 48, 49, 17] If rogue objects were
scattered inward 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[50]), which is poorly constrained from observations.[51, 52]
In addition to the Trappist-1 planets listed in Table 2 we therefore included
either one rogue planet with a randomly-chosen mass between $0.001$ and
$1{\rm\,M_{\oplus}}$, or 1000 equal-mass rogue planetesimals whose masses
totaled between $0.01$ and $3{\rm\,M_{\oplus}}$. The initial semi-major 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 for which the outermost planets are encountered first.[50] We performed
additional sets of simulations in which the semi-major axes of rogue
planetesimals were randomly chosen between either 0.2-3 au or 0.7-1 au (see
Methods Section 5). The initial inclinations of rogue objects were randomly
assigned between zero and ten degrees, which matches the inclination
distribution of small leftover bodies in simulations that form close-in,
Trappist-1-like systems from pebble accretion and migration;[9] we inherently
assume this distribution to apply regardless of the exact migration
scenario.[22, 53, 54, 55, 56, 46, 18, 57, 20, 19, 21, 58] In simulations with
1000 planetesimals, the planetesimals felt the planets’ gravity but not each
others’. Even though all planetesimals were included from the start of the
simulation, the absence of mutual gravity allows us to roughly mimic the
interaction of individual rogue planetesimals with the planets and estimate
their cumulative effects.

Each simulation was integrated for 10 million years but stopped if an
instability occurred, marked by a collision between planets. We used the
hybrid integrator in the publicly-available Mercury integration package[12]
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 $\sim 0.12$ Solar
radii[13, 43]). In terms of calculating collisions, the Trappist-1 planets
were assigned their measured sizes (see Extended Data Table 1) whereas rogue
planets and planetesimals were assigned sizes assuming a bulk density of 3 g
cm-3. Collisions were treated as inelastic mergers conserving linear momentum.

## 3\. 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.[2] From inside to out, 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\varpi_{1,2}\\\ $ (1)

where $\varpi$ are longitudes of pericenter, $\lambda$ are mean longitudes
(equal to the sum of $\varpi$ 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 Eq. 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 3-body (Laplace)
resonant angles. There are five such angles in the system, one for each
adjacent triplet of neighboring planets. Using subscripts appropriate for the
Trappist-1 planets, these are (see Table 1 from Luger et al[2]):

$\displaystyle\phi_{1}=2\lambda_{b}-5\lambda_{c}+3\lambda_{d}$
$\displaystyle\phi_{2}=\lambda_{c}-3\lambda_{d}+2\lambda_{e}$
$\displaystyle\phi_{3}=2\lambda_{d}-5\lambda_{e}+3\lambda_{f}$
$\displaystyle\phi_{4}=\lambda_{e}-3\lambda_{f}+2\lambda_{g}$
$\displaystyle\phi_{5}=\lambda_{f}-2\lambda_{g}+\lambda_{h}$ (2)

Extended Data Fig. 1 demonstrates the multi-resonant 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
million year integration of the fiducial system (orbital parameters in Table 2
of the main text). 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.[59].

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.

Extended Data Fig. 2 presents three examples of the evolution of resonant
behavior from simulations with 1000 rogue planetesimals. In the simulation
interacting with the lowest mass in planetesimals (a total of
$0.0213{\rm\,M_{\oplus}}$; 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{\rm\,M_{\oplus}}$ (center 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 center panels). However, the 3:2 resonance between
planets g and h was compromised, as the resonant angle $\theta_{2}$ no longer
shows clear libration but rather circulates. However, all three 3-body
resonances were maintained, as shown by the continued libration of the Laplace
angles $\phi_{3}$, $\phi_{4}$, and $\phi_{5}$ (bottom three center panels).
While the exact center 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{\rm\,M_{\oplus}}$ (right panel). In that case, all of the resonant
angles associated with mean-motion and Laplace resonances transition from
libration to circulation within roughly 1 Myr of the start of the integration.
While the system remained stable for the 10 Myr length of the integration, it
is possible that it would be unstable on longer timescales.[14]

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{\rm\,M_{\oplus}}$ (left
panel; resonance maintained), $0.0429{\rm\,M_{\oplus}}$ (middle panel;
resonance barely maintained), and $0.0645{\rm\,M_{\oplus}}$ (right panel;
resonance disrupted).

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{\rm\,M_{\oplus}}$. The widening the outer planets’ orbits is small,
and the changes in the planets’ semimajor axes can 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 $\sim 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.

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{\rm\,M_{\oplus}}$. 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.

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 & Lithwick[60]
(see their equations 5, 11 and 15) we estimate the width of the 3:2 mean
motion resonance between planets g and h to be $\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 $\sim
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:

$L=M_{p}\sqrt{GM_{\star}a}\left(1-cos(i)\sqrt{1-e^{2}}\right),$ (3)

where $G$ is the gravitational constant, and $i$ is the orbital inclination.

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 $\sim 0.005-0.01{\rm\,M_{\oplus}}$ in order
to move planet h by the width of its common 3:2 resonance with planet g.
Indeed, this is roughly the minimum mass in a single rogue planet that
destroys the system’s resonant structure (see Fig. 1 in the main text and
Extended Data Fig. 5). However, a total mass in rogue planetesimals of $\sim
0.05{\rm\,M_{\oplus}}$ 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 rather scatters the
planetesimal inward to interact with planet g. The back-reaction of the
scattering causes planet h’s semimajor axis to shift slightly outward.
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 outward together, as is clearly seen in the example from Extended Data
Fig. 3. This process is not perfectly balanced, and 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 outward before they separate enough to leave the
resonance. For this reason it takes substantially more than
$0.01{\rm\,M_{\oplus}}$ 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 multi-resonant 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 Methods Section 7) 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 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{\rm\,M_{\oplus}}$ in rogue planetesimals during the first 10 Myr of its
history would retain its multi-resonant structure. Given the upper limit of
$0.05{\rm\,M_{\oplus}}$ from Fig. 1 in the main text, this system could not
interact with more than an additional $0.03{\rm\,M_{\oplus}}$ on longer
timescales without losing its resonances.

## 4\. 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[50] and
of leftovers from terrestrial planet formation[48], and from the orbital
distribution of analogous leftover material in simulations of the formation of
close-in systems like Trappist-1.[9] Yet it is worth questioning these
choices, and exploring how robust our conclusions are to other plausible
populations of rogue objects (for a discussion see Kral et al[10]).

We therefore performed two additional sets of simulations. In each additional
set the orbits of the planets were unchanged. Rogue planetesimals were
included with the same inclination distribution (randomly selected from zero
to 10 degrees – derived from simulations of the formation of Trappist-1-like
systems[9]) 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 in rogue
planetesimals evenly chosen in the range $0.01-1{\rm\,M_{\oplus}}$.

  * •

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 in rogue planetesimals
evenly chosen in the range $0.01-1{\rm\,M_{\oplus}}$.

All other integration parameters (e.g., timestep, ejection radius, stellar
mass; see Methods Section 2) were the same as for our fiducial set of
simulations.

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.  Planet | Orbital radius | Max. bombard. | Max. bombard. | Max. bombard. | Max. bombard.
---|---|---|---|---|---
| (AU) | mass (${\rm\,M_{\oplus}}$) | mass (${\rm\,M_{\oplus}}$) | mass (${\rm\,M_{\oplus}}$) | mass (${\rm\,M_{\oplus}}$)
|  | [0.07-0.1 au] | [0.2-0.3 au] | [0.7-1 au] | [0.7-1 au]
|  |  |  |  | [super-ecc.]
b | 0.0115 | 0.00035 | 0.00013 | 0.0001 | 0.0082
c | 0.0158 | 0.0015 | 0.0005 | 0.00023 | 0.0073
d | 0.0223 | 0.0015 | 0.00047 | 0.0003 | 0.0033
e | 0.0293 | 0.0034 | 0.00098 | 0.00087 | 0.0044
f | 0.0385 | 0.0078 | 0.0038 | 0.003 | 0.0063
g | 0.0469 | 0.017 | 0.01 | 0.0083 | 0.0074
h | 0.0620 | 0.014 | 0.013 | 0.011 | 0.0045
Extended Data Table 2. 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{\rm\,M_{\oplus}}$ for $a_{rogue}=0.07-0.1$ au and
$a_{rogue}=0.2-0.3$ au, $0.1{\rm\,M_{\oplus}}$ for $a_{rogue}=0.7-1$ au with
the standard periastron distribution, and $0.018{\rm\,M_{\oplus}}$ for
$a_{rogue}=0.7-1$ au with super-eccentric orbits (see Extended Data Fig. 4.)

Extended Data Fig. 4 compares the outcomes of these new sets of simulations
with our fiducial set. There is no clear difference in the figure between
simulations with $a_{rogue}=0.07-0.1$ au and those with $a_{rogue}=0.2-0.3$ au
(for the same periastron distribution). In both cases a total mass in
planetesimals larger than $\sim 0.05{\rm\,M_{\oplus}}$ universally disrupts
the system’s multi-resonant architecture. In contrast, simulations with
$a_{rogue}=0.7-1$ au the critical mass in planetesimals is increased to $\sim
0.1{\rm\,M_{\oplus}}$.

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. Only considering simulations that retained their multi-resonant
structure, 7.1% of planetesimals in our fiducial simulations (with
$a_{rogue}=0.07-0.1$ 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$ au.

The tendency for increased planetesimal ejected for larger $a_{rogue}$ can be
understood from simple dynamical considerations (see also the discussion in
Section 3.2 of Kral et al[10]). 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.[61] Orbital
energy scales with the apoastron distance as $\sim GM_{\star}/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 Eq. 3 in Kral et al[10]). Our simulations
match this expectation, as rogue objects are preferentially accreted in our
fiducial simulations with $a_{rogue}=0.07-0.1$ 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 2). These
upper limits were calculated by tabulating the impact rates on each of the
planets in simulations that preserved their multi-resonant 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 in the main paper that were derived using our fiducial initial
conditions. While other orbital distributions of rogue planetesimals may imply
somewhat higher values for the total mass in 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 2). In addition,
the higher impact speeds of more distant planetesimals may often lead to
erosion rather than accretion.[10]

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_{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{\rm\,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
modeled by Kral et al.[10] The most striking features are the absence of
unstable orbits and the larger critical mass in rogue planetesimals below
which the systems maintained their multi-resonant 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 2). 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 $\sim 0.005{\rm\,M_{\oplus}}$.[4, 5] For planet e the
upper limit on late accretion was similar in both cases ($\approx
0.003-0.004{\rm\,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{\rm\,M_{\oplus}}$ (for planet b;
Extended Data Table 2). 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.[10]

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 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 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 counter-balanced by the much
increased ejection efficiency and the lower accretion efficiency, as
quantified by Kral et al[10]. 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.[10]

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,[62] although the bulk
water contents of the best-studied comet (67P/Churyumov-Gerasimenko) is
roughly constrained to $\sim 25\%$ or lower (see Fig 5b in Choukroun et
al[63]). 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 (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.[10, 30]

## 5\. 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 data set. 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.[3] Long-term simulations show that each of these three
configurations is stable in isolation for at least 1 billion years.

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

  * •

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

  * •

A set of 50 simulations that included 1000 planetesimals, with a total mass of
$0.01-1{\rm\,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 (see
Methods), utilizing the hybrid Mercury integrator.[12]

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.[3] 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).

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{\rm\,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 a single simulation maintained
its multi-resonant structure with a rogue planet more massive than
$0.01{\rm\,M_{\oplus}}$. From all sets together, only a single simulation
remained in a multi-resonant configuration after interacting with a rogue
planet more massive than $0.01{\rm\,M_{\oplus}}$ and even in that case (in Set
1) the rogue planet was barely more massive than the Moon
($M_{Rogue}=0.0167{\rm\,M_{\oplus}}=1.4M_{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 just $\sim 2\%$ of our simulations. We therefore consider our
upper limit for a single late giant impactor of a lunar mass to be robust.

## 6\. 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 $\sim 5-10$ kyr for each planet,
consistent with previous work.[10] 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_{rogue}=0.07-0.1$ au), as
well as for the simulations with larger initial semimajor axes described above
(see Methods Section 4 and Extended Data Fig. 4).

Extended Data Fig. 6. Distribution of impact speeds of planetesimals in
simulations in which the planets remained in resonance on each of the
Trappist-1 planets.

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.[10, 11] 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_{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_{\star}/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.[10] 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 that on the outer planets. The trend of higher impact speeds for
larger initial rogue semimajor axes $a_{rogue}$ can be explained by the second
term, as rogue objects had fixed initial periastron distances such that larger
$a_{rogue}$ corresponds to larger eccentricity.

Kral et al[10] 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 $\sim 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 due to 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 (see 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[10]). 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.

## 7\. Tidal evolution

One might wonder whether resonances could be re-formed by orbital migration
driven by tidal interactions between the planets and the central star. In
order to obtain a convergent migration with tides, one needs an 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 separately the stellar tide
induced migration for planet g and h using an equilibrium tide formalism[64],
which accounts for the fact that the stellar radius evolves with time (we used
the evolution tracks of Baraffe et al. (2015)[65] for a stellar mass of
$0.09~{}M_{\odot}$). We computed as a function of time the evolution of the
quantity $\Delta a$, which is the evolution of the separation between the two
planets $a_{h}(t)-a_{g}(t)$ with respect to the initial separation
$a_{h}(0)-a_{g}(0)$. No tidal migration would lead to $\Delta a=0$, a
convergent migration would lead to $\Delta a<0$ and a divergent migration to
$\Delta a>0$. Due 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. Besides, the migration timescale highly
depends on the assumed stellar dissipation. Which is why we treat stellar
dissipation as a free parameter and explore a range of $[10^{-3},...,10^{3}]$
times the dissipation of a fully convective star[66, 67]. Note that due to the
lack of a radiative core for low-mass stars, we do not expect the dynamical
tide to play a role as important as for Sun-like stars[68]. 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.

Extended Data Fig. 7.  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_{init}$, where $t_{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.

The results for an initial time of $3$ Myr can been seen in Extended Data Fig.
7. Due 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 in resonance after a few million years,
it never does. 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 in resonance.

Our reasoning assumes that the inner planets do not play a role on the outward
migration of planet 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[69, 70]) 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.[71] If the star spun fast enough during
its early evolution, the corotation radius would have been very close-in and
tidally-induced migration could have been directed outward for all of the
planets (see Bolmont et al[72] for a discussion of the effect of the stellar
spin history on tidal migration). While the calculations presented here do not
account for that possibility, it is something to be addressed in future,
detailed modeling of the tidal evolution of the system.

We have also assumed here that only star-planet tides influence the orbits of
the planets, meaning that we have neglected other processes like planet-planet
tides. However, it has been shown by Hay & Matsuyama (2019)[73] that the
migration planet-planet tides induce is negligible compared to the star-planet
tides. Moreover, planet-planet tides would induce an inward migration, which
would rather 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.

## 8\. Internal structure modeling 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[40, 1] and masses.[74, 3] 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.[75, 26, 76]

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 fit the measured density of the planet with the highest
(uncompressed) density, planet c.

  * •

In model (iii) we assumed bulk molar Fe/Mg of 0.75 and Mg/Si of 1.02 that was
suggested by Unterborn et al[26] 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.[77]

  * •

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

As suggested by Dorn et al.[25], 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
adaption of model (ii) but assumes 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[24] indeed
suggest that the three inner planets may be depleted in water and other
volatiles[3] (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) are both constrained by the stellar proxy[26] that is based on
refractory abundances of F-G-K stars with metallicities similar to Trappist-1.
While model (iii) assumes all of each planet’s iron is found in the core,
model (v) has all the planets’ iron oxidized in their mantles with core-free
interiors, in line with recent estimates of oxygen budgets of extrasolar
rocks.[78] Model (iv) is for reference and assumes solar relative refractory
abundances.

The tested models (i-v) were then used to calculate the possible water mass
fractions (wmf) needed to fit the planetary mass and radius data as well as
stellar irradiation. To accomplish this we calculated the wmf for all
planetary mass-radius pairs of the 104 posterior samples as in Agol et al.[3],
thereby accounting for non-Gaussian distributed data errors. Whenever the
sample’s density is equal or higher than the pure rocky interior, we set the
wmf of the sample to zero.

Interior profiles are calculated similar to previous studies.[23, 3]
Specifically, we use equations of state for pure iron (hcp),[79] mantle
silicates,[76] solid or liquid water[80] and steam water.[81] Following Agol
et al.[3], water is treated as a condensed layer (liquid and/or solid) for
planets less irradiated than the runaway greenhouse insolation threshold
(i.e., Trappist-1 e,f,g,and h) ; water is treated as steam for planets more
irradiated than the runaway greenhouse insolation threshold (i.e., Trappist-1
b,c,d),[82, 24] The atmospheric structure and radiative balance were computed
using the LMD Generic 1-D climate model.[82, 24]

We neglected hydrated mantle rocks as they only affects the planet radius by
less than 2 % at maximum[83], such that water is assumed to only be present 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 lead
marginal changes in planet radii ($<$1%) compared with pure Fe cores while
keeping the same bulk Fe/Si ratios.

Tables containing all of the results of our modeling are included in the
Supplementary Information (Supplementary Tables 2-6).

## 9\. Description of the simulation from Fig. 2 in the main text that formed
our illustrative Trappist-1 ‘analog.’

The Trappist-1 analog planets in Figure 2 from the main text represent a
super-Earth system produced in simulations of Izidoro et al[9]. Their planet
formation model includes the effects of pebble accretion, type-I migration,
eccentricity and inclination damping due to tidal interaction of
protoplanetary objects with the disk. This particular system shown was
produced in the setup of Model-III of Izidoro et al. It started with $\sim 70$
Moon-mass 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 snowline 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 gaseous
disk lifetime a system of 7 planets formed in a multi-resonant configuration
between 0.05 au and 0.24 au. The planets’ masses were 0.77 $M_{\oplus}$, 0.92
$M_{\oplus}$, 1.03 $M_{\oplus}$, 1.42 $M_{\oplus}$, 1.09 $M_{\oplus}$, 0.49
$M_{\oplus}$, and 0.42 $M_{\oplus}$. 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 $\sim 14\%$ water
by mass, the sixth planet was dry, and the outermost planet contained $\sim
20\%$ water by mass. The disk dissipated after 5 Myr and we continued the
simulation in a gas-free environment for additional 50 Myr to verify that the
resonant configuration remained dynamically stable.

The primary goal of Figure 2 in the main text 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 of all,
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 Papaloizou et al[19],
Coleman et al[20] and Lin et al[58] for careful studies of the resonant
structure of simulated Trappist-1 analogs). 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 analog does have outstanding 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.

Code availability
The Mercury N-body integrator[12] 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.
<EMAIL_ADDRESS>upon reasonable request, and the atmospheric code is
based on polynomial fits from Turbet et al (2020).[24]

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.

This paper is dedicated to the memory of our friend and colleague Franck
Hersant.

We thank Sylvie Rousseau and Arnaud 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 Aurelie
Guilbert-Lepoutre and Anastasis Gkotsinas for input on the water contents of
comets. S. N. R. and F. S. thank the CNRS’s PNP program for support. A. I and
R. D. acknowledge NASA grant 80NSSC18K0828 for financial support during
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 CA 1795/3. E.A. and S.N.R. acknowledge support from
from the Virtual Planetary Laboratory Team, a member of the NASA Nexus for
Exoplanet System Science, funded via NASA Astrobiology Program Grant No.
80NSSC18K0829. This project has received funding from the European Union’s
Horizon 2020 research and innovation program under the Marie Sklodowska-Curie
Grant Agreement No. 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 the allocation 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 PZ00P2_174028.

S. N. R. initiated the project after discussions with L. C., P. B. and A. I.,
ran and analyzed the simulations using codes optimized by A. I. and
configurations of the system provided by S. G. and E. A.. C. D. and M T.
performed the interior modeling to obtain constraints on the water contents of
the Trappist-1 planets. E. B. performed the analysis of 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. G.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to S.N.R
(rayray.sean@gmail.com).

## Supplementary Information

## Exact orbital configurations of the Trappist-1 system that were used in our
analysis

As explained in the Methods, we tested four different best-fit orbital
configurations of the system, each drawn from Agol et al.[3] They are listed
in Table 1.

Planet | Mass | Radius | Semimajor | Eccentricity | Longitude of | Mean
---|---|---|---|---|---|---
| (${\rm\,M_{\oplus}}$) | (${\rm\,R_{\oplus}}$) | Axis $a$ (AU) | $e$ | periastron $\varpi$ (∘) | Anomaly $M$ (∘)
Fiducial (Set 1) |  |  |  |  |  |
b | 1.3925 | 1.1174 | 0.011551 | 0.002344 | 253.61247 | 105.78489
c | 1.2943 | 1.0967 | 0.015820 | 0.001224 | 132.62793 | 54.89836
d | 0.3958 | 0.7880 | 0.02229 | 0.005045 | 202.45580 | 171.39157
e | 0.6824 | 0.9200 | 0.02930 | 0.007013 | 52.42997 | 30.97582
f | 1.0634 | 1.0448 | 0.038551 | 0.008298 | 170.04247 | 247.44087
g | 1.3464 | 1.1294 | 0.046896 | 0.003760 | 355.97714 | 87.27858
h | 0.3198 | 0.7552 | 0.061963 | 0.003571 | 172.18673 | 118.58431
Set 2 |  |  |  |  |  |
b | 1.3925 | 1.1174 | 0.011551 | 0.004677 | 181.37652 | 178.45848
c | 1.2896 | 1.0979 | 0.015821 | 0.002755 | 276.77544 | 270.77957
d | 0.3909 | 0.7889 | 0.022289 | 0.007069 | 217.02820 | 156.90568
e | 0.7068 | 0.9210 | 0.029296 | 0.005359 | 44.52192 | 38.91996
f | 1.0477 | 1.0460 | 0.038551 | 0.009050 | 179.50878 | 238.00290
g | 1.3232 | 1.1307 | 0.046896 | 0.003686 | 335.57244 | 107.74799
h | 0.3309 | 0.7561 | 0.061965 | 0.003714 | 188.97120 | 101.87152
Set 3 |  |  |  |  |  |
b | 1.3405 | 1.1174 | 0.011550 | 0.003402 | 194.66082 | 165.01640
c | 1.3062 | 1.0979 | 0.015818 | 0.000802 | 196.57200 | 351.02606
d | 0.3987 | 0.7889 | 0.022290 | 0.006354 | 210.05726 | 163.88358
e | 0.7003 | 0.9210 | 0.029295 | 0.005947 | 51.87353 | 31.58265
f | 1.0777 | 1.0460 | 0.038549 | 0.008950 | 175.24378 | 242.26801
g | 1.3652 | 1.1307 | 0.046897 | 0.003526 | 347.44533 | 95.88219
h | 0.3404 | 0.7561 | 0.061966 | 0.003855 | 172.27644 | 118.53050
Set 4 |  |  |  |  |  |
b | 1.3370 | 1.1174 | 0.011551 | 0.000611 | 167.48106 | 191.83734
c | 1.2811 | 1.0979 | 0.015818 | 0.004240 | 330.78097 | 216.35048
d | 0.3803 | 0.7889 | 0.022288 | 0.005124 | 210.85020 | 162.86095
e | 0.7112 | 0.9210 | 0.029294 | 0.007072 | 45.29093 | 38.04304
f | 1.0464 | 1.0460 | 0.038551 | 0.007802 | 173.00821 | 244.39623
g | 1.3231 | 1.1307 | 0.046892 | 0.004384 | 352.07474 | 91.11644
h | 0.3050 | 0.7561 | 0.061960 | 0.003833 | 182.31057 | 108.43892
Table 2: Starting orbital configurations of the Trappist-1 system that were
used in our simulations. The fiducial system is at the top (and referred to as
Set 1 in our tests). In all systems 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 2-dimensional.
Each system was taken from the best-fit posteriors of Agol et al.[3]

## Interior structure modeling to estimate the water contents of the
Trappist-1 planets

In Tables 2-6 we include the full results of our interior structure modeling
described in the Methods. Each of the five tables contains the derived water
contents (including the means, 1-$\sigma$ and $2-\sigma$ errors) for each
planet.

Planet | Mean water content | Uncertainties ($1\sigma$) | Uncertainties ($2\sigma$)
---|---|---|---
| by mass |  |
b | – | – | –
c | – | – | –
d | – | – | –
e | $6.88\times 10^{-3}$ | ${}^{+1.58\times 10^{-2}}_{-6.88\times 10^{-3}}$ | ${}^{+4.31\times 10^{-2}}_{-6.88\times 10^{-3}}$
f | $2.05\times 10^{-2}$ | ${}^{+1.69\times 10^{-2}}_{-1.20\times 10^{-2}}$ | ${}^{+4.40\times 10^{-2}}_{-2.05\times 10^{-2}}$
g | $3.64\times 10^{-2}$ | ${}^{+1.84\times 10^{-2}}_{-1.42\times 10^{-2}}$ | ${}^{+4.82\times 10^{-2}}_{-2.66\times 10^{-2}}$
h | $3.53\times 10^{-2}$ | ${}^{+4.37\times 10^{-2}}_{-3.02\times 10^{-2}}$ | ${}^{+9.57\times 10^{-2}}_{-3.53\times 10^{-2}}$
Table 3: Water contents from interior model i. Planet | Mean water content | Uncertainties ($1\sigma$) | Uncertainties ($2\sigma$)
---|---|---|---
| by mass |  |
b | – | – | ${}^{+3.11\times 10^{-5}}$
c | – | – | –
d | – | – | –
e | $1.60\times 10^{-2}$ | ${}^{+1.88\times 10^{-2}}_{-1.21\times 10^{-2}}$ | ${}^{+4.66\times 10^{-2}}_{-1.60\times 10^{-2}}$
f | $3.25\times 10^{-2}$ | ${}^{+1.74\times 10^{-2}}_{-1.36\times 10^{-2}}$ | ${}^{+4.50\times 10^{-2}}_{-2.43\times 10^{-2}}$
g | $4.90\times 10^{-2}$ | ${}^{+1.89\times 10^{-2}}_{-1.47\times 10^{-2}}$ | ${}^{+4.94\times 10^{-2}}_{-2.78\times 10^{-2}}$
h | $4.70\times 10^{-2}$ | ${}^{+4.46\times 10^{-2}}_{-3.37\times 10^{-2}}$ | ${}^{+9.65\times 10^{-2}}_{-4.70\times 10^{-2}}$
Table 4: Water contents from interior model ii. Planet | Mean water content | Uncertainties ($1\sigma$) | Uncertainties ($2\sigma$)
---|---|---|---
| by mass |  |
b | – | ${}^{+1.05\times 10^{-5}}$ | ${}^{+6.89\times 10^{-5}}$
c | – | – | ${}^{+2.69\times 10^{-5}}$
d | – | – | –
e | $2.75\times 10^{-2}$ | ${}^{+2.06\times 10^{-2}}_{-1.57\times 10^{-2}}$ | ${}^{+5.03\times 10^{-2}}_{-2.59\times 10^{-2}}$
f | $4.59\times 10^{-2}$ | ${}^{+1.88\times 10^{-2}}_{-1.46\times 10^{-2}}$ | ${}^{+4.84\times 10^{-2}}_{-2.76\times 10^{-2}}$
g | $6.40\times 10^{-2}$ | ${}^{+2.03\times 10^{-2}}_{-1.60\times 10^{-2}}$ | ${}^{+5.31\times 10^{-2}}_{-3.03\times 10^{-2}}$
h | $6.11\times 10^{-2}$ | ${}^{+4.68\times 10^{-2}}_{-3.79\times 10^{-2}}$ | ${}^{+1.01\times 10^{-1}}_{-6.11\times 10^{-2}}$
Table 5: Water contents from interior model iii. Planet | Mean water content | Uncertainties ($1\sigma$) | Uncertainties ($2\sigma$)
---|---|---|---
| by mass |  |
b | – | ${}^{+1.65\times 10^{-5}}$ | ${}^{+1.25\times 10^{-4}}$
c | – | – | ${}^{+5.11\times 10^{-5}}$
d | – | – | –
e | $3.69\times 10^{-2}$ | ${}^{+2.14\times 10^{-2}}_{-1.76\times 10^{-2}}$ | ${}^{+5.17\times 10^{-2}}_{-2.96\times 10^{-2}}$
f | $5.61\times 10^{-2}$ | ${}^{+1.95\times 10^{-2}}_{-1.54\times 10^{-2}}$ | ${}^{+5.00\times 10^{-2}}_{-2.87\times 10^{-2}}$
g | $7.49\times 10^{-2}$ | ${}^{+2.12\times 10^{-2}}_{-1.64\times 10^{-2}}$ | ${}^{+5.56\times 10^{-2}}_{-3.14\times 10^{-2}}$
h | $7.16\times 10^{-2}$ | ${}^{+4.77\times 10^{-2}}_{-4.02\times 10^{-2}}$ | ${}^{+1.03\times 10^{-1}}_{-6.91\times 10^{-2}}$
Table 6: Water contents from interior model iv. Planet | Mean water content | Uncertainties ($1\sigma$) | Uncertainties ($2\sigma$)
---|---|---|---
| by mass |  |
b | – | – |
c | – | – | –
d | – | – | –
e | – | ${}^{+2.55\times 10^{-2}}$ | ${}^{+4.92\times 10^{-2}}$
f | $7.01\times 10^{-3}$ | ${}^{+1.39\times 10^{-2}}_{-7.01\times 10^{-3}}$ | ${}^{+3.91\times 10^{-2}}_{-7.01\times 10^{-3}}$
g | $2.01\times 10^{-2}$ | ${}^{+1.70\times 10^{-2}}_{-1.20\times 10^{-2}}$ | ${}^{+4.50\times 10^{-2}}_{-2.01\times 10^{-2}}$
h | $2.13\times 10^{-2}$ | ${}^{+3.98\times 10^{-2}}_{-2.13\times 10^{-2}}$ | ${}^{+9.04\times 10^{-2}}_{-2.13\times 10^{-2}}$
Table 7: Water contents from interior model v.

## References

  * [1] Gillon, M. _et al._ Seven temperate terrestrial planets around the nearby ultracool dwarf star TRAPPIST-1.  _Nature_ 542, 456–460 (2017).  1703.01424.
  * [2] Luger, R. _et al._ A seven-planet resonant chain in TRAPPIST-1.  _Nature Astronomy_ 1, 0129 (2017).  1703.04166.
  * [3] Agol, E. _et al._ Refining the Transit-timing and Photometric Analysis of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides.  _The Planetary Science Journal_ 2, 1 (2021).  2010.01074.
  * [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– (2007).
  * [5] Walker, R. J.  Highly siderophile elements in the Earth, Moon and Mars: Update and implications for planetary accretion and differentiation.  _Chemie der Erde / 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.  _A &A_ 650, A152 (2021).  1902.08772.
  * [10] Kral, Q. _et al._ Cometary impactors on the TRAPPIST-1 planets can destroy all planetary atmospheres and rebuild secondary atmospheres on planets f, g, and h.  _MNRAS_ 479, 2649–2672 (2018).  1802.05034.
  * [11] Došović, V., Novaković, B., Vukotić, B. & Ćirković, M. M.  Water transport throughout the TRAPPIST-1 system: the role of planetesimals.  _MNRAS_ 499, 4626–4637 (2020).  2007.05366.
  * [12] Chambers, J. E.  A hybrid symplectic integrator that permits close encounters between massive bodies.  _MNRAS_ 304, 793–799 (1999).
  * [13] Gonzales, E. C. _et al._ A Reanalysis of the Fundamental Parameters and Age of TRAPPIST-1.  _ApJ_ 886, 131 (2019).  1909.13859.
  * [14] Tamayo, D., Rein, H., Petrovich, C. & Murray, N.  Convergent Migration Renders TRAPPIST-1 Long-lived.  _ApJ_ 840, L19 (2017).  1704.02957.
  * [15] Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N. & Walsh, K. J.  Building Terrestrial Planets.  _Annual Review of Earth and Planetary Sciences_ 40, 251–275 (2012).  1208.4694.
  * [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.  _Nature Geoscience_ 11, 77–81 (2018).
  * [18] Ormel, C. W., Liu, B. & Schoonenberg, D.  Formation of TRAPPIST-1 and other compact systems.  _A &A_ 604, A1 (2017).  1703.06924.
  * [19] Papaloizou, J. C. B., Szuszkiewicz, E. & Terquem, C.  The TRAPPIST-1 system: orbital evolution, tidal dissipation, formation and habitability.  _MNRAS_ 476, 5032–5056 (2018).  1711.07932.
  * [20] Coleman, G. A. L., Leleu, A., Alibert, Y. & Benz, W.  Pebbles versus planetesimals: the case of Trappist-1.  _A &A_ 631, A7 (2019).  1908.04166.
  * [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.  _MNRAS_ 491, 1998–2009 (2020).  1909.12320.
  * [22] Terquem, C. & Papaloizou, J. C. B.  Migration and the Formation of Systems of Hot Super-Earths and Neptunes.  _ApJ_ 654, 1110–1120 (2007).  arXiv:astro-ph/0609779.
  * [23] Dorn, C. _et al._ A generalized Bayesian inference method for constraining the interiors of super Earths and sub-Neptunes.  _A &A_ 597, A37 (2017).  1609.03908.
  * [24] Turbet, M. _et al._ Revised mass-radius relationships for water-rich rocky planets more irradiated than the runaway greenhouse limit.  _A &A_ 638, A41 (2020).  1911.08878.
  * [25] Dorn, C., Mosegaard, K., Grimm, S. L. & Alibert, Y.  Interior Characterization in Multiplanetary Systems: TRAPPIST-1.  _ApJ_ 865, 20 (2018).  1808.01803.
  * [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.  _Nature Astronomy_ 2, 297–302 (2018).  1706.02689.
  * [27] Barth, P. _et al._ Magma ocean evolution of the trappist-1 planets.  _Astrobiology_ 21 (2021).  URL https://doi.org/10.1089/ast.2020.2277.  PMID: 34314604, %****␣TRAP1_final_complete.tex␣Line␣925␣****https://doi.org/10.1089/ast.2020.2277.
  * [28] Acuña, L. _et al._ Characterisation of the hydrospheres of TRAPPIST-1 planets.  _A &A_ 647, A53 (2021).  2101.08172.
  * [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).  1406.6435.
  * [31] Marty, B.  The origins and concentrations of water, carbon, nitrogen and noble gases on Earth.  _Earth and Planetary Science Letters_ 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.  _ApJ_ 807, 9 (2015).  1505.03516.
  * [33] Ramirez, R. M. & Kaltenegger, L.  The Habitable Zones of Pre-main-sequence Stars.  _ApJ_ 797, L25 (2014).  1412.1764.
  * [34] Tian, F. & Ida, S.  Water contents of Earth-mass planets around M dwarfs.  _Nature Geoscience_ 8, 177–180 (2015).
  * [35] Luger, R. & Barnes, R.  Extreme Water Loss and Abiotic O2Buildup on Planets Throughout the Habitable Zones of M Dwarfs.  _Astrobiology_ 15, 119–143 (2015).  1411.7412.
  * [36] Bolmont, E. _et al._ Water loss from terrestrial planets orbiting ultracool dwarfs: implications for the planets of TRAPPIST-1.  _MNRAS_ 464, 3728–3741 (2017).  1605.00616.
  * [37] Kislyakova, K. G. _et al._ Magma oceans and enhanced volcanism on TRAPPIST-1 planets due to induction heating.  _Nature Astronomy_ 1, 878–885 (2017).  1710.08761.
  * [38] Leleu, A. _et al._ Six transiting planets and a chain of Laplace resonances in TOI-178.  _arXiv e-prints_ arXiv:2101.09260 (2021).  2101.09260.
  * [39] Mills, S. M. _et al._ A resonant chain of four transiting, sub-Neptune planets.  _Nature_ 533, 509–512 (2016).  1612.07376.
  * [40] Gillon, M. _et al._ Temperate Earth-sized planets transiting a nearby ultracool dwarf star.  _Nature_ 533, 221–224 (2016).  1605.07211.
  * [41] Burgasser, A. J. & Mamajek, E. E.  On the Age of the TRAPPIST-1 System.  _ApJ_ 845, 110 (2017).  1706.02018.
  * [42] 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.  _ApJ_ 871, 63 (2019).  1811.06938.
  * [43] Ducrot, E. _et al._ TRAPPIST-1: Global results of the Spitzer Exploration Science Program Red Worlds.  _A &A_ 640, A112 (2020).  2006.13826.
  * [44] Lee, M. H. & Peale, S. J.  Dynamics and Origin of the 2:1 Orbital Resonances of the GJ 876 Planets.  _ApJ_ 567, 596–609 (2002).
  * [45] Cresswell, P., Dirksen, G., Kley, W. & Nelson, R. P.  On the evolution of eccentric and inclined protoplanets embedded in protoplanetary disks.  _A &A_ 473, 329–342 (2007).  0707.2225.
  * [46] Izidoro, A. _et al._ Breaking the chains: hot super-Earth systems from migration and disruption of compact resonant chains.  _MNRAS_ 470, 1750–1770 (2017).  1703.03634.
  * [47] 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– (2010).
  * [48] 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).  1306.4325.
  * [49] Morbidelli, A. & Wood, B. J.  Late Accretion and the Late Veneer.  _Washington DC American Geophysical Union Geophysical Monograph Series_ 212, 71–82 (2015).  1411.4563.
  * [50] Levison, H. F. & Duncan, M. J.  The long-term dynamical behavior of short-period comets.  _Icarus_ 108, 18–36 (1994).
  * [51] Boss, A. P. _et al._ Astrometric Constraints on the Masses of Long-period Gas Giant Planets in the TRAPPIST-1 Planetary System.  _AJ_ 154, 103 (2017).  1708.02200.
  * [52] Jontof-Hutter, D., Truong, V. H., Ford, E. B., Robertson, P. & Terrien, R. C.  Dynamical Constraints on Nontransiting Planets Orbiting TRAPPIST-1.  _AJ_ 155, 239 (2018).
  * [53] Ogihara, M. & Ida, S.  N-Body Simulations of Planetary Accretion Around M Dwarf Stars.  _ApJ_ 699, 824–838 (2009).  0904.4543.
  * [54] McNeil, D. S. & Nelson, R. P.  On the formation of hot Neptunes and super-Earths.  _MNRAS_ 401, 1691–1708 (2010).  0910.5299.
  * [55] 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.  _ApJ_ 719, 810–830 (2010).  1006.2584.
  * [56] Cossou, C., Raymond, S. N., Hersant, F. & Pierens, A.  Hot super-Earths and giant planet cores from different migration histories.  _A &A_ 569, A56 (2014).  1407.6011.
  * [57] Schoonenberg, D., Liu, B., Ormel, C. W. & Dorn, C.  Pebble-driven planet formation for TRAPPIST-1 and other compact systems.  _A &A_ 627, A149 (2019).  1906.00669.
  * [58] Lin, Y.-C., Matsumoto, Y. & Gu, P.-G.  Formation of Multiple-planet Systems in Resonant Chains around M Dwarfs.  _ApJ_ 907, 81 (2021).  2012.03469.
  * [59] Murray, C. D. & Dermott, S. F.  _Solar system dynamics_ (Cambridge university press, 1999).
  * [60] Hadden, S. & Lithwick, Y.  A Criterion for the Onset of Chaos in Systems of Two Eccentric Planets.  _AJ_ 156, 95 (2018).  1803.08510.
  * [61] 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.  _MNRAS_ 464, 3385–3407 (2017).  1610.00714.
  * [62] A’Hearn, M. F.  Comets as Building Blocks.  _ARA &A_ 49, 281–299 (2011).
  * [63] 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).
  * [64] Bolmont, E., Raymond, S. N. & Leconte, J.  Tidal evolution of planets around brown dwarfs.  _A &A_ 535, A94 (2011).  1109.2906.
  * [65] 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.  _A &A_ 577, A42 (2015).  1503.04107.
  * [66] Hansen, B. M. S.  Calibration of Equilibrium Tide Theory for Extrasolar Planet Systems.  _ApJ_ 723, 285–299 (2010).  1009.3027.
  * [67] Bolmont, E., Raymond, S. N., Leconte, J. & Matt, S. P.  Effect of the stellar spin history on the tidal evolution of close-in planets.  _A &A_ 544, A124 (2012).  1207.2127.
  * [68] Bolmont, E. & Mathis, S.  Effect of the rotation and tidal dissipation history of stars on the evolution of close-in planets.  _Celestial Mechanics and Dynamical Astronomy_ 126, 275–296 (2016).  1603.06268.
  * [69] Blanco-Cuaresma, S. & Bolmont, E.  Studying Tidal Effects In Planetary Systems With Posidonius. A N-Body Simulator Written In Rust.  In _EWASS Special Session 4 (2017): Star-planet interactions (EWASS-SS4-2017)_ (2017).  1712.01281.
  * [70] Bolmont, E. _et al._ Impact of tides on the transit-timing fits to the TRAPPIST-1 system.  _A &A_ 635, A117 (2020).  2002.02015.
  * [71] Herbst, W., Eislöffel, J., Mundt, R. & Scholz, A.  The Rotation of Young Low-Mass Stars and Brown Dwarfs.  In Reipurth, B., Jewitt, D. & Keil, K. (eds.) _Protostars and Planets V_ , 297 (2007).  astro-ph/0603673.
  * [72] Bolmont, E., Raymond, S. N., Leconte, J. & Matt, S. P.  Effect of the stellar spin history on the tidal evolution of close-in planets.  _A &A_ 544, A124 (2012).  1207.2127.
  * [73] Hay, H. C. F. C. & Matsuyama, I.  Tides Between the TRAPPIST-1 Planets.  _ApJ_ 875, 22 (2019).  1903.04501.
  * [74] Grimm, S. L. _et al._ The nature of the TRAPPIST-1 exoplanets.  _A &A_ 613, A68 (2018).  1802.01377.
  * [75] Dorn, C. _et al._ Can we constrain the interior structure of rocky exoplanets from mass and radius measurements?  _A &A_ 577, A83 (2015).  1502.03605.
  * [76] Sotin, C., Grasset, O. & Mocquet, A.  Mass–radius curve for extrasolar earth-like planets and ocean planets.  _Icarus_ 191, 337–351 (2007).
  * [77] Lodders, K., Palme, H. & Gail, H.-P.  Abundances of the elements in the solar system.  _arXiv preprint arXiv:0901.1149_ (2009).
  * [78] 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).
  * [79] Hakim, K. _et al._ A new ab initio equation of state of hcp-fe and its implication on the interior structure and mass-radius relations of rocky super-earths.  _Icarus_ 313, 61–78 (2018).
  * [80] Haldemann, J., Alibert, Y., Mordasini, C. & Benz, W.  Aqua: a collection of h2o equations of state for planetary models.  _Astronomy & Astrophysics_ 643, A105 (2020).
  * [81] 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 Publ. Corp., Washington, D. C., 1984).
  * [82] 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.  _A &A_ 628, A12 (2019).  1906.03527.
  * [83] Shah, O., Alibert, Y., Helled, R. & Mezger, K.  Internal water storage capacity of terrestrial planets and the effect of hydration on the mr relation.  _arXiv preprint arXiv:2012.06455_ (2020).
  * [84] Mamajek, E. E.  Initial Conditions of Planet Formation: Lifetimes of Primordial Disks.  In Usuda, T., Tamura, M. & Ishii, M. (eds.) _American Institute of Physics Conference Series_ , vol. 1158 of _American Institute of Physics Conference Series_ , 3–10 (2009).  0906.5011.
  * [85] Turbet, M. _et al._ Modeling climate diversity, tidal dynamics and the fate of volatiles on TRAPPIST-1 planets.  _A &A_ 612, A86 (2018).  1707.06927.