Cometary impactors on the TRAPPIST-1 planets can destroy all planetary atmospheres and rebuild secondary atmospheres on planets f, g, h

The TRAPPIST-1 system is unique in that it has a chain of seven terrestrial Earth-like planets located close to or in its habitable zone. In this paper, we study the effect of potential cometary impacts on the TRAPPIST-1 planets and how they would affect the primordial atmospheres of these planets. We consider both atmospheric mass loss and volatile delivery with a view to assessing whether any sort of life has a chance to develop. We ran N-body simulations to investigate the orbital evolution of potential impacting comets, to determine which planets are more likely to be impacted and the distributions of impact velocities. We consider three scenarios that could potentially throw comets into the inner region (i.e within 0.1au where the seven planets are located) from an (as yet undetected) outer belt similar to the Kuiper belt or an Oort cloud: Planet scattering, the Kozai-Lidov mechanism and Galactic tides. For the different scenarios, we quantify, for each planet, how much atmospheric mass is lost and what mass of volatiles can be delivered over the age of the system depending on the mass scattered out of the outer belt. We find that the resulting high velocity impacts can easily destroy the primordial atmospheres of all seven planets, even if the mass scattered from the outer belt is as low as that of the Kuiper belt. However, we find that the atmospheres of the outermost planets f, g and h can also easily be replenished with cometary volatiles (e.g. $\sim$ an Earth ocean mass of water could be delivered). These scenarios would thus imply that the atmospheres of these outermost planets could be more massive than those of the innermost planets, and have volatiles-enriched composition.


INTRODUCTION
The nearby (d=12pc) M8V ultra-cool dwarf star TRAPPIST-1 (2MASS J23062928-0502285) is now known to be surrounded by at least seven terrestrial-like planets Luger et al. 2017). This old (7.6±2.2 Gyr, Burgasser & Mamajek 2017), close-by, multi-planetary system may offer one of our best chances to study the chemistry, and structure of terrestrial planet atmospheres outside our Solar System (de Wit et al. 2016;Morley et al. 2017). Moreover, several of the TRAPPIST-1 E-mail: qkral@ast.cam.ac.uk planets (most likely planets e, f and g, Gillon et al. 2017) lie within the liquid water habitable zone (HZ, e.g. O'Malley-James & Kaltenegger 2017). However, the presence of liquid water and possible life strongly depends on the atmospheric content of these planets, the presence of oceans, the vegetation coverage, etc. (e.g. Wolf 2017; Alberti et al. 2017;Ehlmann et al. 2016;Carone et al. 2016;Godolt et al. 2016).
This system being very close-by, we may soon be able to start characterising the atmospheres of the seven planets with new telescopes such as JWST ( search for tracers of life. Such detailed spectral characterisation may eventually allow us to infer the presence of biological activity via the detection of gases such as ozone (Barstow & Irwin 2016), or via the spectral signatures of pigmented micro-organisms (Poch et al. 2017). Regardless, such observations will inform on the atmospheric compositions of these planets that is necessary to study the possibility that life may develop.
For now, little is known about the atmospheres of these seven planets. The two innermost planets b and c have been observed using transmission spectroscopy (de Wit et al. 2016). This showed that the combined spectrum of both planets (obtained when transiting at the same time) is featureless, which favours atmospheres that are tenous (composed of a variety of chemical species), not hydrogendominated, dominated by aerosols or non-existent. Similar conclusions have been made for planets d, e, f (and potentially g) by de Wit et al. (2018). Also, from the combined measurement of planet radii (transit) and masses (transit timing variations), the derived planets' densities show that TRAPPIST-1 b, d, f, g, and h may require envelopes of volatiles in the form of thick atmospheres, oceans, or ice (Grimm et al. 2018). We thus do not know yet information that would be important for considering the habitability of the planets such as whether these planets' atmospheres are primordial or created later, for instance by cometary impacts, although current observations suggest that current atmospheres may not be primordial due to a lack of hydrogen signatures in the observed spectra (de Wit et al. 2018).
Previous theoretical studies of the atmospheric composition of the TRAPPIST-1 planets have shown that they may vary with time and be affected by the early evolution of the star. Indeed, ultra-cool dwarfs such as TRAPPIST-1 take up to 1Gyr to cool down (Baraffe et al. 2015) and reach the main-sequence after the planets formed. This means that planets that are today in the HZ would have undergone a very hot pre-main-sequence era (with potentially a runaway greenhouse phase) and may have lost all (or part) of their initial water content (Bolmont et al. 2017). Moreover, Bourrier et al. (2017) find that the total XUV emission from the star might be strong enough to entirely strip the primordial atmospheres of the planets over a few Gyr. One could then expect that a few of the TRAPPIST-1 planets are devoid of atmospheres, or left with a gas layer too tenuous for life to persist over long timescales (Roettenbacher & Kane 2017).
Here we consider another process that can strongly influence the atmospheres, both positively and negatively for life: exocomets. Impacting exocomets can influence planetary atmospheres in multiple ways: a) they can provide an energy source that depletes primordial atmospheres. b) They might also deliver volatiles that subsequently replenish a secondary atmosphere (i.e., dry, depleted atmospheres from impacts or XUV irradiation could be replenished by later impacts, and surviving primordial atmospheres could see their elemental abundances significantly transformed via exocomet impacts). c) Impacting exocomets may also act as catalysts for the development of life. Indeed, these impacts may initiate a cascade of chemical reactions, some of which can produce the necessary precursors to nucleobases on these planets (Saladino et  For now, there is no evidence of exocomets in the TRAPPIST-1 system, however, this does not mean they are not present and part of the motivation of this work is to determine if evidence for such a population may be imprinted on the planets' atmospheres. Many stars have large outer reservoirs of planetesimals that produce a detectable infrared excess due to collisional production of dust (Wyatt 2008;Eiroa et al. 2013). Detections of CO gas in several systems are used to infer that these planetesimals are icy with a composition that is similar to Solar System comets (e.g. Kral et al. 2016;Marino et al. 2016;Matrà et al. 2017a). These planetesimal belts are harder to detect around low mass stars such as TRAPPIST-1 due to their low luminosity but this does not mean they are not present (Plavchan et al. 2009;Theissen & West 2014). Some stars also have evidence that comets from these outer regions are being scattered into the inner regions. For example, CO detected at 20au in η Corvi is inferred to originate in the sublimation of such an exocomet population ). In addition, high-velocity metallic gas absorption lines in some systems (Montgomery & Welsh 2012;Kiefer et al. 2014;Eiroa et al. 2016) are inferred to originate in very eccentric comets passing very close to their host star (called falling evaporating bodies, e.g. Beust et al. 1990). Thus, it is not unreasonable that TRAPPIST-1 has (or indeed may have had) comets at some level.
In this study, we hypothesize that such comets exist in the TRAPPIST-1 system and use previous studies that looked at the effect of impacts onto planetary atmospheres (e.g. de Niem et al. 2012;Schlichting et al. 2015), and especially hydrodynamical simulations (Shuvalov 2009) to derive some constraints on the TRAPPIST-1 planets' atmospheres in the presence of impacting comets. We start by estimating the possible mass of a planetesimal belt that could have survived around TRAPPIST-1. In Sec. 3, we then study the dynamics of comets in the TRAPPIST-1 system that come close to the seven planets, i.e within 0.1au. Notably, we look into which planet will receive most impacts, at which velocity and derive the timescales on which impacts happen. In Sec. 4, we describe three plausible scenarios that can potentially scatter many exocomets over the lifetime of the system. In Sec. 5, we show the results of our model, i.e how much atmospheric mass is removed from the primordial atmospheres of the seven planets by a continuous series of impacts and evaluate whether those impacts increase or reduce the amount of volatiles in the planets' atmospheres, and what kind of atmosphere each planet is likely to end up with. We then discuss our results in terms of their implications for the development of life in Sec. 6 before concluding.

THE POSSIBLE PRESENCE OF A DISC AROUND TRAPPIST-1
This paper is based on the potential presence of a yet undetected debris disc around TRAPPIST-1. To consider what this debris disc might look like, we construct a minimum mass extrasolar nebula for the TRAPPIST-1 system similar to Hayashi (1981), or Chiang & Laughlin (2013) who used 1925 extrasolar planets to constrain the minimum surface densities at different distances from the star assuming planets formed in situ.
10 −2 10 −1 10 0 10 1 10 2 Distance to the star (au) Surface density (kg/m 2 ) data Fit Extrapolation After 7Gyr Figure 1. Surface density in the TRAPPIST-1 system assuming a minimum mass extrasolar nebula and extrapolating to tens of au to obtain a plausible mass that would be left in a potential, yet undetected, belt. In red, we show the predicted profile after 7Gyr of collisional evolution.
To get a surface density for each planet, we take the planet mass and divide it by the area of the annulus around the planet. For planets c to g, we define the annulus as being between the two midpoints to the neighbouring planets. For planets b and h, we work out the half width using the centres between planets b and c and between planets g and h and multiply that width by two. This gives the following surface density (in solids) after fitting the data (see Fig. 1) where r is the distance to the star. Our fit of Σ provides values a factor 4 smaller than Chiang & Laughlin (2013) at 1au (who used a large sample of Kepler planets around earlier-type stars) but steeper in r and very close to the fit by Gaidos (2017) who did it specifically for M-dwarf Kepler planets. It is less than a factor 2 from the minimum mass solar nebula (MMSN) in solids for terrestrial planets at 1au (Hayashi 1981). The H2O iceline during planetesimal formation is estimated to have been close to ∼0.1au in the TRAPPIST-1 system (Ormel et al. 2017). It could as well have been slightly closer-in (by a factor ∼2) based on the (still not-well constrained) gradient of water compositions of the 7 planets (Unterborn et al. 2018). We assumed that beyond 0.1au, the solid {rock+ice} surface density is a factor 4 higher following Hayashi (1981). We can now extrapolate the mass that may be present at several au and potentially form a disc of planetesimals rather than planets.
A planetesimal belt at a radius r with dr/r ∼ 0.5 would have a mass of ∼ 12.6(r/1au) 0.03 M⊕. The collisional lifetime of the biggest planetesimals in such a belt is given by 1 tc = 1.4 × 10 −3 r 13/3 (dr/r)DcQ 5/6 D e −5/3 M −4/3 /Mtot 1 We note that this formula can be used when the largest bodies yr (Wyatt 2008). This gives tc ∼ 4 × 10 3 (r/1au) 4.3 yr by assuming typical values (as in Wyatt 2008;Kral et al. 2017c, i.e e = 0.05, Q D = 500J/kg, Dc=100km). In other words, a belt at 1au would be significantly depleted after 7Gyr (the age of the system) of collisional evolution and we expect any belt this close in to have been significantly depleted. However, a belt at >28au could survive over 7Gyr. At shorter radii, the mass that remains after collisional evolution for 7Gyr would be expected to have a radial profile that scales ∝ r 7/3 (Kennedy & Wyatt 2010) as shown by the red dotted line in Fig. 1. While this formula depends on many uncertain parameters, it shows that we expect any potential surviving belt to be located at 10au. Using the extrapolation in Eq. 1, we expect such a leftover belt between 10 and 50au to have a mass of ∼20M⊕, which is compatible with the predicted large initial mass of the protoplanetary disc around TRAPPIST-1 required to have formed the seven planets (Haworth et al. 2018). While this is at least two orders of magnitude more massive than the Kuiper belt (Fraser & Kavelaars 2009;Vitense et al. 2010), note that the Kuiper belt is thought to have formed much more massive, with a solid mass of 20-40M⊕ compatible with the MMSN (e.g. Weidenschilling 1977;Hayashi 1981;Nesvorný & Morbidelli 2012, but see Shannon et al. 2016 for a dissenting view). The left-over belt is not expected to extend much farther than 50-100au because protoplanetary discs around low-mass stars are less extended than around T-Tauri stars (Hendler et al. 2017). One caveat to this estimate is that our approach is only accurate for an in situ formation of the seven planets. For planets that formed further out close to the water iceline as suggested by Ormel et al. (2017), the surface density would go down by a factor 10 at most and so would the belt mass leading to an estimate of 2M⊕.
The only observation of TRAPPIST-1 in the infrared is by WISE at 22µm (Patel et al. 2017), which shows no signs of infrared excess. However, any belt that is warm enough to emit at 22µm would have to be inside 10au and so, as noted above, would be expected to be collisionally depleted. The only region where significant mass is expected to remain at 7Gyr is beyond 10au, where such a belt would be < 15K (assuming a black body) and therefore its emission would peak at λ > 340µm. This WISE observation is thus not constraining and observations at longer wavelengths are required to constrain such a cold belt, for instance using the ALMA interferometer.

DYNAMICS OF IMPACTS FOR COMETS COMING FROM AN OUTER BELT
There are many possible origins for the comets that may impact planets b to h. Rather than studying the details of the specific evolution for each scenario, we will start by assuming that very eccentric comets are produced and we will study their dynamics and look at their interactions with the seven planets. This framework is therefore general as soon from the collisional cascade have a large enough collision velocity that they can fragment after an impact. Depending on the level of stirring, using this formula for radii 50au is therefore not accurate and only gives a lower limit on the timescale tc.
as eccentric comets are produced and will be tied to specific scenarios (planet scattering, Kozai-Lidov mechanism or Galactic tides) in Sec. 4. The pericentre q of the eccentric comets we model can reach a few tenths to hundredths of au where they can collide with one of the seven detected planets around TRAPPIST-1. The apocentre Q can vary from a few au (for comets that originate in close-in belts) to > 100s of au for comets coming from very cold outer belts or exo-Oort clouds. We perform N-body simulations of these very eccentric orbits assuming that the evolution is dominated by perturbations from the known TRAPPIST-1 planets to understand how their fate depends on the comet's orbital parameters q (pericentre) and Q (apocentre). That is, for each of these different comet families (i.e for a given set of {q,Q}) we determine the fraction that is accreted onto the different planets and the fraction that is ejected. We also compute impact velocities for each family of comets, which are used in Sec. 5 to assess if cometary impacts are able to destroy planetary atmospheres and if delivery of volatiles from these comets is possible.

N-body simulations of impacts with the seven planets
The N-body simulations are run with REBOUND (Rein & Liu 2012) with the Hermes integrator which combines the IAS15 integrator for close encounters within a few Hill radii of the planets (Rein & Spiegel 2015) and uses the WHFast integrator otherwise (Rein & Tamayo 2015). The simulations include the seven planets orbiting around the central star TRAPPIST-1 (see Tab. 1 for the parameters used). We use a timestep of 5% of planet b's orbital timescale. We assumed zero eccentricities for the planets as the 2σ upper limits are low (< 0.09 as implied by tidal forces and orbital stability, Gillon et al. 2017;Tamayo et al. 2017;Quarles et al. 2017). The planets gravitationally interact with each other, but their orbits do not evolve significantly over the course of all our simulations. We start each simulation with 2000 test particles that all have a similar pericentre and apocentre {q,Q} spread in a narrow range defined by a grid (see Fig. 2). We run the simulations until all test particles have either been ejected from the system (i.e., if their positions go beyond 100 times the initial comet's apocentre) or accreted onto the planets or the star. We note however that almost no particles collide with the central star. This is because for high-eccentricity orbits the pericentres will be locked and for low-eccentricity orbits, we notice that there are very few scattering events that could potentially send the comets onto the star. Rather, the particles tend to be accreted or ejected by the planet close to their pericentres. We assumed zero inclination, which we expect to be unrealistic but leads to much faster simulations and can be scaled a posteriori to give results for a comet-like inclination distribution (see subsection 3.5). Running the simulations assuming a zero inclination angle was necessary to allow the simulations to be performed in a reasonable timescale (i.e., not exceeding two months). The whole set of 900 simulations took ∼2 months on 20 CPUs, whereas inclined comets would have taken almost two years to compute. This is because we ran each simulation until there are no particles left. As the time to accrete/eject particles is much smaller in the zero inclination case, we gain a factor greater than 10 in overall computational time. We note that of the results we derive in this section, the probability map as well as the accretion/ejection timescales are affected (in a quantifiable way) by a change in inclination but not the impact velocities (subsection 3.4).
We ran a grid of 900 N-body simulations for a wide range of {q,Q} values, with 90 logarithmically-spaced bins in pericentre covering 10 −3 au < q < 10 −1 au and 10 logarithmically-spaced bins in apocentre covering 10 −1 au < Q < 10 2 au, which form the grid seen in Fig. 2. The grid is defined by the pericentres and apocentres at the start of the simulations. The TRAPPIST-1 planets are located between 0.01 and 0.06au (white vertical lines in Fig. 2) so that the chosen range of pericentres is large enough to follow what happens when the comets' orbits cross those of the planets.
3.2 Probability to impact the different planets or to be ejected for comet-like orbits Fig. 2 shows a map of the probability to impact the different planets (each inset is for a given planet, planet b to h from left to right), while Fig. 3 shows the probability to be ejected for each given {q,Q} of our parameter space. Some of the large scale features in these figures can be readily understood.
For example, the extended black regions in Fig. 2 at large pericentres are because in order for a comet to collide with a given planet, the comet's pericentre must be smaller than the planet's semi-major axis a pla . Since the pericentre and apocentre of comets do not evolve significantly from their starting values, this means that the region of the parameter space with q > a pla appears in black. Comets with such pericentres collide with the more distant planets.
Another large scale feature is that the probability to impact one of the planets is higher for smaller cometary apocentres. This can be explained by looking at Fig. 3 which shows that the ejection rate goes up with increasing Q, noting that the sum of the impact probabilities over the seven planets and the ejection probability equals 1.
The increased ejection probability seen in Fig. 3 with Q (for all pericentres) is because the comet's energy (∼ GM /Q) is lower for these larger apocentres and so a comet is ejected by a smaller kick when passing by a planet. The biggest kick in velocity that the comet can receive from a planet (without colliding onto it) is roughly equal to vesc , where vesc is the escape velocity of the planet. The resulting increase in the comet's orbital energy can be enough to unbind the comet if v orb vesc > GM /Q, where the comet's orbital velocity v orb close to a planet is which leads to the apocentre value Q where ejections start becoming dominant The mass of planet h is not well constrained (Luger et al. 2017) and as its radius is similar to that of planet d, we assume the same mass as planet d.  where M pla , R pla and a pla are the planet mass, radius and semi-major axis, respectively. This calculation explains why Fig. 3 shows that for Q 1au, ejection is the more likely outcome.
Another feature in Fig. 2 is that for pericentres inside planet b, the accretion probability is higher for planets closer to the star. In fact, the accretion probability decreases as a −1 pla from planet b to h for a fixed {q, Q} in this regime. This can be explained by the different accretion timescales of each planet as showed in Sec. 3.3.
Finally, another noticeable feature is that the highest probabilities of impacts (the narrow yellow regions) are for comet orbits that have pericentres close to but slightly smaller than the positions of the planets. For instance, on the planet d inset, we see that the yellow region is concentrated in a narrow region of the parameter space between 0.015−0.021au (planets c and d positions). This can be readily explained because comets with such a pericentre cannot collide with planets b, c so increasing the rate of collisions with planets d, e, f, g, h. The most extreme case is for comets that have pericentres just below planet h, thus ensuring that they can only collide with planet h and explaining the very narrow yellow region in the planet h inset.

Accretion/ejection timescales
It is important to consider the timescale on which particles are accreted (or ejected) in the simulations, because we will later be considering how these outcomes compete with other processes that may be acting to modify the particles' orbits (such as the processes that brought them onto comet-like orbits in the first place).
In Fig. 4, we plot the loss timescale t loss which is the timescale for half of the 2000 test particles to be lost from the simulation (through accretion or ejection) as a function of the apocentre Q. Since there is little dependence on the pericentre of the comets' orbits (because the comet velocity is almost independent of q), this shows that for Q 1au, the loss timescale scales as Q 3/2 and as Q 1/2 for larger Q. For Q 1au, the loss of particles is dominated by accretion onto the planets. For a 2D geometry, the rate of collisions between a given comet and planet is proportional to R col = nσσ2Dv rel , where nσ is the fraction of the comet's orbit per unit cross-section spent in a region dr around the planet's orbit, v rel is the relative velocity at encounter and σ2D is the collisional cross section. Considering the fraction of the orbit spent in an annulus dr around the planet's orbit we find that nσ (per cross-sectional area) is ∝ Q −3/2 a −1/2 pla . In practice, the velocity at encounter v rel is the same as  The white colour is used to show the part of the parameter space that have pericentres too far from the planets to either collide or be ejected. The black colour is for orbits that collide and eject particles but with a very low ejection probability 10 −2 .
10 −1 10 0 10 1 10 2 Q (au)   . Accretion timescale tacc = R −1 col as a function of the semi-major axes of the planets a pla for an i = 0 • (crosses) and a realistic comet-like inclination distribution (filled dots). We show the analytically predicted tacc (blue line) for an inclined distribution of comets from the i = 0 • case (see subsection 3.5). Here, tacc is plotted for planet b for Q ∼ 1au and it scales as Q 3/2 . the impact velocity vimp, and we show in Sec. 3.4 that vimp is close to the comet's velocity (see Eq. 2), which is large enough for gravitational focusing to be ignored such that σ2D = 2R pla . Therefore, we find that R col ∝ Q −3/2 a −1 pla , so that the accretion timescale is tacc ∝ R −1 col ∝ Q 3/2 a pla , explaining why the loss timescale scales as Q 3/2 for small Q. It also shows that the accretion timescale scales as a pla as shown in Fig. 5, where we plot tacc for planet i by computing t loss /pi, where pi is the probability to be accreted on planet i (that we have for every {q, Q} cell in Fig. 2). This also explains why the accretion probability (∝ t −1 acc ) decreases as a −1 pla from planet b to h as noted in Sec. 3.2. For Q 1au, the loss is dominated by ejections. In that case, the cross section σej used to calculate the rate of ejection is proportional to the impact parameter bej at which encounters are just strong enough to cause ejection. The kick ∆v that the comet receives from a planet after a close encounter scales with 1/b, and for the ejection to happen vcom∆v > GM /Q (see Sec. 3.2). This means that for a flat geometry σej ∝ Q and so tej ∝ (nσσejv rel ) −1 ∝ Q 1/2 , explaining the dependencies.

Impact velocities for the different planets
An important parameter to determine the effects of a cometary impact onto the atmosphere of a planet is the impact velocity.
In Fig. 6, we show histograms of impact velocities for the different planets. We computed the impact velocities for each simulation (i.e. for specific pericentres and apocentres), but find that the distributions of impact velocity do not depend significantly on the comet's pericentre. To get Fig. 6, we therefore average the vimp distributions over the pericentres in the grid, assuming that comet orbits are uniform in log q (keeping a fixed apocentre). Averaging in this way results in more accurate histograms of impact velocities. To do so, the impact velocities from the different simulations are weighted by the probability to impact the different planets (using Fig. 2). Furthermore, Fig. 7 shows that the medians of the impact velocity distributions for each planet also do not depend significantly on apocentre. Thus while Fig. 6 shows the distributions for an apocentre of ∼1au, these distributions are also representative of that of a large range of apocentres.
We see that the impact velocity distribution is peaked at a different location for each planet from ∼ 15 to ∼110km/s from planet h to b. A much smaller secondary peak can also be seen for each planet. This is because there are two extreme types of impacts. Collisions can occur when the comet is on a near radial orbit approaching or receding from pericentre. They may also occur when the planet and comet velocities are parallel (i.e when the comet encounters the planet near its pericentre). As shown by Eq. 2, the comet velocity at impact is ∼ 2GM /a pla (for a pla a), which is thus always higher than the planet's Keplerian velocity of GM /a pla (which varies from ∼35km/s for the farthest to ∼80km/s for the closest planet). Therefore, we find that the impact velocity distributions should peak at GM /a pla ( √ 2 − 1) (33, 29, 24, 21, 18, 16, 14km/s, for planet b to h) for parallel orbits at impact and would be maximal for radial encounters at 3GM /a pla . We note that the impact velocities are much greater than the escape velocities of the planets (∼10km/s) and therefore gravitational focusing is not important.
Thus the high velocity peaks correspond to comets colliding on radial orbits and the low velocity peaks to comets falling on the planets at their pericentres (i.e., parallel collision). By looking at Fig. 2, we see that for planet h, the highest impact probability region (the yellow region) is very narrow and restricted to comets whose pericentres are close to planet h's position so that most collisions are going to be parallel. This explains why the low velocity peak is higher for this planet. For planet b, however, the yellow region is large and not peaked close to planet b's semi-major axis. Therefore, most impacts will happen with comets on nearly radial orbits and the high velocity peak is therefore higher than the low velocity peak. Histograms for the other planets can be understood following the same procedure.
The non-dependence of impact velocities on apocentres shown in Fig. 7 also derives from the velocity at impact, which, as shown by Eq. 2, only depends on a pla and not a. We notice that the median velocities are close to the Keplerian velocities of the corresponding planets (Fig. 7).

Simulations for realistic inclinations
The simulations assumed comets with zero inclination. To check how our results change for different inclinations, we ran a set of 30 additional simulations (spread across the {q, Q} parameter space) with more realistic comet-like inclinations. The chosen inclination distribution follows a Rayleigh distribution peaking at 10 degrees, i.e close to the distribution of JFC 2 comets Di Sisto et al. 2009). We find that the loss timescale (see Fig. 4, filled dots) and the timescale for accretion onto the different planets are affected (see Fig. 5, filled dots), but that the impact velocities are unaffected. The difference in tacc between the inclined and flat cases in Fig. 5 can be explained by generalising the analytics in subsection 3.3. The ratio of the rates at which comets collide with a planet is expected to be (nvσ3D)/(nσσ2D) = πR pla /(4Imaxa pla ), where nv is now number per unit volume in the vicinity of the planet and Imax the median inclination of the comets in the 3D case. We plot this analytical prediction (blue line) together with some numerical simulations for a distribution of inclinations (filled dots) in Fig. 5. A similar comparison shows that the dependence in Q remains the same for tacc for the two types of simulations. We thus conclude that the zero inclination simulation collisional rates can be scaled to account for the inclined case.
To recover the probability map shown in Fig. 2 for the case of an inclined distribution, we also need to rescale the loss timescale t loss shown in Fig. 4 because the probability pi to be accreted on a given planet i is equal to t loss /tacc. The results for t loss from the inclined numerical simulations are shown in Fig. 4 (filled dots). We can predict the change in ejection timescale (which dominates t loss at Q 0.2au) for the 3D case from the 2D simulations in the same way as for the timescales for accretion onto the planets. This prediction is that ejection timescale should be longer by 0.8(M /M pla )(a pla Imax)/Q. This is reasonably accurate but 10 −1 10 0 10 1 10 2 Apocentre Q (au) we prefer to use the numerical ratio of t loss for the inclined and zero-inclination cases which is best fit by a power law equal to 63 (Imax/10 • ) Q −0.61 . Using these different scalings, we calculate the new probability map (see Fig. 8) and use that the sum of the probabilities to be accreted onto each of the planets and to be ejected equals 1 to compute a new ejection map (see Fig. 9). Comparing the predictions from our scalings to the different results from the inclined distribution simulations, we find that we are accurate within a factor 2.

DIFFERENT SCENARIOS TO MAKE ECCENTRIC COMETS
We have studied the dynamics of highly eccentric comets in the presence of the seven TRAPPIST-1 planets in the previous section. Here, we consider three different scenarios that can send comets from the outer regions of the TRAPPIST-1 system onto such eccentric orbits (see Fig. 10); 1) A planetesimal disc is perturbed by a nearby planet and comets are scattered inwards by this single planet or through a chain of planets (similar to comets scattered in our Solar System, e.g. Bonsor et al. 2012;Marino et al. 2018). 2) A distant companion to TRAPPIST-1 forces comets in a Kuiper belt-like disc to undergo Kozai-Lidov oscillations (e.g. Nesvold et al. 2016), which can bring the comets to very close pericentres. 3) Galatic tides perturb a far away exo-Oort-cloud and send comets to decreasing pericentres. We assume that the evolution of comets' orbits in these three scenarios can be approximated as an evolution in which their apocentres Q remain constant and their pericentres q decrease at a constant rateq. This approximation allows us to use the results of Sec. 3 to consider the outcome for comets scattered into the inner regions without having to consider the detailed dynamics of the comets' origin. The simplified dynamics allows us to study a wide range of different possible scenarios. Owing to this simplification, the results are expected to give order of magnitude correct estimates, which is justified by the uncertainties on the presence of a belt in this system and its yet unknown properties. We explore expectations for the differentq values for each of these three scenarios below.

Impacts from comets scattered by a single or a series of planets
In our Solar System, comets from the Kuiper belt are thrown into the inner Solar System thanks to a series of scattering by different planets. This planet scattering scenario has been invoked multiple times (Nesvorný et al. 2010;Booth et al. 2009;Bonsor et al. 2014) to try to explain the presence of hot dust around many stars (see Kral et al. 2017a, for a review). More recently, Marino et al. (2018) studied the effect of scattering by a chain of planets for a large parameter space so as to understand which planetary systems are more suited to create large hot dust levels or maximise impacts on the chain of planets. They ran simulations for 1Gyr for different chains of planets with semi-major axes ranging from 1 to 50au and planet masses ranging from a few to 100M⊕. In our case, we consider that interactions with the innermost planet of the chain dominate the comets' dynamical evolution as they reach to the very small pericentres considered here. However, we also have to consider that some fraction of comets would have been ejected or accreted by the other planets before reaching the innermost planet of the chain. Marino et al. (2018) show that for a wide range of planet chain architectures, fin = 1 − 7% of comets originating in an outer belt end up reaching the inner system. For a closein belt similar to the debated (see MacGregor et al. 2018) belt recently invoked around the M-dwarf Proxima Centauri (Anglada et al. 2017), a single planet at 1au could be enough to scatter comets into the inner regions (but we note that we showed in Sec. 2 that such a close-in belt is not likely to have survived around TRAPPIST-1 unless it formed recently). In that case, no comets are lost on the way through the chain and fin = 1. We only consider the case of a planet coplanar to the 7 planets as this system seems well-aligned. A non-coplanar configuration would lead to an increased inclination distribution, which effect could be quantified using the analytics from Sec. 3.5. Conservation of the Tisserand parameter 3 (Murray & Dermott 1999) means that comets being scattered by an innermost planet that is on a circular orbit can only reach down to a minimum pericentre (see Bonsor et al. 2012). Thus to scatter comets to small enough pericentres to reach the TRAPPIST-1 planets' locations, we consider that the innermost planet must be on an eccentric orbit, since in that case there is no minimum pericentre constraint (see also Frewen & Hansen 2014). We also consider that this planet   should not be too massive, so as not to eject the comets before they can reach the innermost parts of the system.
Guided by the results of Frewen & Hansen (2014) we consider 1 and 10M⊕ planets with a 0.4 eccentricity orbiting at 1au to be representative of the kind of planets that are able to put comets on orbits that are capable of colliding with the seven known TRAPPIST-1 planets. We note that such planets are not massive enough and not close enough to gravitationally disturb the orbits of the seven currently known planets as can be checked directly from Read & Wyatt (2016), so that the system of the seven inner planets stays stable even in the presence of such an additional planet (see also Quarles et al. 2017). We also note that these planet masses agree with current mass upper limits by Boss et al.  (2017) (i.e < 4.6MJup within a 1 yr period, and < 1.6MJup within a 5 yr period). While an eccentricity of 0.4 is above the median eccentricity found for Earth to Super-Earth mass planets, such eccentric planets are observed. We also note that our scenario would still work for lower eccentricities as described in Frewen & Hansen (2014) but theq value would vary.
Moreover, we find that the outermost planet interact-ing with the belt and causing the scattering could migrate outwards if it has a mass 10M⊕ (see Eq. 58 in Ormel et al. 2012, where we used the surface density of the potential surviving belt shown in Fig. 1) and stall if it is more massive. Such a migration is beneficial to sending more comets inwards as shown in Bonsor et al. (2014) because more time is available for the scattering process and it can access more material to scatter from. However, for a planet mass 0.1M⊕, the scattering would not be efficient anymore as shown by Eq. 52 of Ormel et al. (2012). Too massive an outer planet would also prevent material from being scattered inwards as it would be more likely ejected but disc evolution models find that having another Jupiter in the TRAPPIST-1 system is unlikely (Haworth et al. 2018). Marino et al. (2018) show that for planet masses 100M⊕, 2% of the scattered comets still reach the inner region. More massive planets such as a Jupiter would more likely eject most of the material (see Wyatt et al. 2017).
We ran N-body simulations to follow the evolution of test particles initially randomly located in the chaotic zone (where resonances overlap, as classically defined in Wisdom 1980) of such an eccentric planet. The planet is located at 1au and simulations are run until the planet runs out of material to scatter among the initial 5000 particles in the chaotic zone. Fig. 11 shows the evolution of the distribution of pericentres of particles that have their apocentres at the 10M⊕ planet location, which decreases steadily with time. Quantifying the rate of this decrease by looking at the evolution of the median of the distribution, we find thaṫ qP ∼ 5 × 10 −5 au/yr, over many orbits for the 10M⊕ case. Running another simulation for an Earth mass planet with similar eccentricity at the same location, we find thatqP ∼ 10 −5 au/yr. While the path of an individual comet could be somewhat stochastic through the parameter space (jumping in individual scattering events), the effect for an ensemble of comets is that of a slow inward migration of q. Therefore, we model this population, and how it is depleted due to interactions with the seven inner planets, by assuming that comets have an apocentre Q that is fixed at the position of the innermost planet of the chain (1au), and considering the various depletion pathways as the comets cross the parameter space in Figs. 2 and 3 at a constant rateq. That rate depends on the mass of the planet, so we keep this as a free parameter, noting that Eq. 4 gives realistic values. That rate has a strong influence on the outcome (see Sec. 5).

Impacts from comets undergoing Kozai-Lidov oscillations due to an outer companion
The incidence of binaries around M-dwarfs is around 27% (Janson et al. 2012). Comets located at tens of hundreds of au (either in a disc or a more spherical Oort-cloud like distribution) could be perturbed by such distant companion stars. If the mutual inclination i0 of some comets with this companion is greater than 39.23 degrees, the so-called Kozai-Lidov cycle can start and the mutual inclination starts decreasing while the eccentricity of the comet increases to reach a maximum (Kozai 1962;Lidov 1962  For the case of a circular outer companion, the maximum eccentricity reached by the comets is given by 1 − (5/3) cos 2 (i0) (Innanen et al. 1997). This means that to reach a pericentre q < q1 = 10 −2 au (to be able to reach the seven planets), the initial mutual inclination should be greater than i0 > arccos 3 5 (1 − (1 − q1/a) 2 ) , where a is the semi-major axes of the comets. If the belt is really closein, i.e a = 0.1au, this corresponds to i0 > 70.3 degrees, and at 100au, it gives i0 > 89.4 degrees (i.e., an almost perpendicular orbit is necessary in this latter case). We note that this inclination is between the perturber and comets and the latter can be inclined compared to the planets. Therefore, finely tuned companions would be needed to send comets to the right location.
However, for an eccentric outer companion, the comets' eccentricity can reach values arbitrarily close to 1 (e.g. Lithwick & Naoz 2011;Teyssandier et al. 2013). While the periastron precession (due to GR) can dominate the dynamics (Liu et al. 2015) and stop the Kozai mechanism from working when the eccentricity comes close to 1, this occurs only for pericentres interior to the known TRAPPIST-1 planets (see Eq. 50 and Fig. 6 of Liu et al. 2015).
The Kozai oscillations will occur even if the perturbing companion is very distant and/or not very massive, only the timescale to achieve the eccentricity change will be longer in that case. Assuming the initial eccentricities e0 of comets in a disc are small then the timescale TK to reach the maximum eccentricity given an eccentric perturber is (Antognini 2015; Naoz 2016) where ac, ec, Mc and Pc are, respectively, the semi-major axis, eccentricity, mass and orbital period of the companion and P is the orbital period of the comet being perturbed (i.e P = 2π a 3 /(GM )). The parameter oct = ec(a/ac)/(1 − e 2 c ) quantifies the relative size of the octupole term of the Hamiltonian compared to the quadrupole term and is not equal to zero for an eccentric perturber (the timescale for a circular orbit can be found by substituting oct ∼ 259, see Antognini 2015). Then, we can determine an order of magnitude for the rate of pericentre evolution given byqK Therefore, considering a belt at 100au perturbed by an eccentric companion of mass 0.01M at 150au 4 , we find thaṫ qK ∼ 2 × 10 −4 au/yr. While a much farther companion could decrease that value and a farther exo-Kuiper belt would in-creaseqK , we consider in Sec. 5.1 how evolution at typical q might affect the planetary atmospheres of TRAPPIST-1 planets.
We have also checked that for such a configuration the Kozai mechanism cannot be suppressed by the precession induced by unknown planets in the system. Imagining the worst case scenario of the presence of a Jupiter-mass planet at 10au in this system, the Kozai dynamics remains dominated by the outer companion if the belt is located further than ∼30au (Petrovich & Muñoz 2017), which is assumed here. We also looked at the effect of the seven known TRAPPIST-1 planets on the Kozai mechanism. The bodies that can reach these planets must be very eccentric and to take that into account properly, we model the effect of these planets as being an effective J2 (quadrupole moment) and check whether the precession rate due to J2, i.e. ωJ 2 is able to counteract the precession due to Kozai. Using Eq. 35 in Fabrycky & Tremaine (2007), we find that the effective J2 of the 7 TRAPPIST-1 planets starts contributing and reduce the maximum Kozai eccentricity for where a7 is the semi-major axis of the outermost 4 We note that such a low mass companion at large distances is not yet ruled out ( TRAPPIST-1 planet, and e is the eccentricity of the comet, which for a belt of semi-major axis a should be 1 − 0.01/a to be able to reach the innermost planet or 1 − 0.06/a to reach the outermost one. Using Batygin & Morbidelli (2017), we find that the J2 due to the 7 TRAPPIST-1 planets would be ∼ 2 × 10 −5 . Therefore, from Eq. 7, we estimate that for a 0.01M at 150au, the belt of planetesimals should be at 10au to be able to reach the outermost planet or at 70au to reach the innermost one.
We acknowledge that the change in inclination while undergoing Kozai oscillations is not taken into account in our previous general simulations shown in Sec. 3. However, depending on the exact inclination of the companion compared to the belt, we can quantify using the equations given in Sec. 3.5 how it will affect the probability to be accreted onto the planets, which scales as I −1 max .

Impacts from Oort-cloud comets perturbed by Galactic tides
TRAPPIST-1 may have an Oort cloud, either because comets were captured from their neighbouring stars' belts at the cluster stage , or because comets were scattered out by its planetary system (Tremaine 1993). In our Solar System, Duncan et al. (1987) propose that leftover comets between Uranus and Neptune would be thrown onto more extended orbits by the two planets until they reach a semi-major axis of ∼ 5000au where Galactic tides change their angular momentum, therefore moving their periastron from reach of the planets. While planets that are efficient at forming Oort clouds need to have the right ranges of mass and semi-major axis, which does not include the known TRAPPIST-1 planets (e.g. Wyatt et al. 2017), other (as yet unseen) planets in the system could have scattered material into such an Oort cloud. Moreover, an Oort-cloud forming planet does not necessarily need to be at this exact location now as it could have migrated. The same mechanism, i.e. Galactic tides, which increased angular momentum of leftover comets pumped up by Uranus and Neptune and thus detaching the comet orbits from the planets can also decrease angular momentum and bring back an outer Oort-cloud comet to the planetary system. For an Oort cloud, the Galactic tidal force (due to the Galactic disc potential) will slowly make the comets lose angular momentum resulting in a slow drift inwards of pericentre (because e increases, a is constant) at a rateqG (e.g. Heisler & Tremaine 1986;Matese & Whitman 1992;Veras & Evans 2013). The eccentricity reaches a maximum that is given by Breiter & Ratajczak (2005), which is greater for comets perpendicular to the orbital plane. It is usually assumed that when a comet reaches a few au, it is lost from the Oort cloud due to planetary perturbations (e.g. Heisler & Tremaine 1986;Fouchard et al. 2006).
The value ofqG can be estimated from the mean square change in angular momentum per orbit ∆J 2 = 1.2 × 10 −29 ρ 2 0 a 7 /M (in au 4 /yr 2 , Eq. A4 of Wyatt et al. 2017), with ρ0 the stellar mass density in units of 0.1M /pc 3 (local stellar mass density), a in au and M in M . Thus, sincė We note that thisqG value is very small but it varies strongly with a.
Therefore, for the case of the TRAPPIST-1 planets, it means that if the location of the Oort cloud is closer than a few 10 3 au, the time for moving the bulk of Oort-cloud bodies down to small pericentres close to the planet positions (i.e., < 0.06au) would be greater than the age of the system. However, for an Oort-cloud at 10 5 au, it only takes ∼5Myr to reach the inner region but an origin at such a large distance becomes unlikely given that such comets would have been stripped by passing stars (Tremaine 1993). Given the age and low-mass of TRAPPIST-1, comets with a semi-major axis beyond 2000au should be strongly depleted by passing stars but some may still remain. We also note that the presence of massive Jupiter-like planets in the outer regions of the TRAPPIST-1 system would have strong effects on the dynamics of the system (e.g. Kaib & Quinn 2009) and our prescription would need to be revised if this type of planet is discovered.

For a generalq
Here we use the results from sections 3.1 and 3.5 to determine how much material will be accreted on the different planets depending on how fast comets move inwards, which is assumed to be set by a constant rate of change of pericentreq. For a givenq and apocentre Q, each (pericentre) q cell of our parameter space is progressively crossed as the comet moves inwards to smaller q values. Taking into account the timescales shown in Fig. 4 and scalings from Sec. 3.5 (as we consider realistic inclined comets), we can use Figs. 2 and 3 (or their counterparts Figs. 8 and 9) to work out the fraction of comets that are accreted onto the different planets or ejected along the way. Hereafter we only consider the inclined case using the results from Sec. 3.5. Fig. 12 shows the fraction accreted on the different planets facc for Q = 1 and 100au. For smallq (i.e < 10 −5 au/yr for the Q = 1au case) most comets end up on planets g (yellow curve) and h (brown), while for largeq each planet gets a fraction of comets accreted. We also show the fraction ejected as grey lines for Q = 1 and 100au, which is close to 1 for very smallq and decreases for larger values, meaning that comets can go past the planets without being ejected nor accreted onto the planets forq > 10 −5 au/yr. These comets may end up on the star or collisionally deplete before reaching it.
We see that the fraction of comets accreted 5 facc onto the different planets varies significantly from 0.05 to < 10 −8 for 10 −6 <q < 1au/yr. For smallq, the fraction accreted is dominated by planets h and g because q decreases so slowly that these outermost planets catch all impacting comets before they reach further in. On the other hand, for largeq, the comets cannot efficiently accrete on the planets (as the loss timescale is long compared to q/q, see Fig. 4) and end up at small radii (where they either accrete onto the star or deplete collisionally). In between these two regimes, each planet accretes a fraction of the scattered comets. The fraction of comets accreted facc is also higher for smaller apocentres, as expected from Fig. 2. We find that for largeq, facc ∝q −1 Q −1 . The fraction accreted by the different planets vary by one order of magnitude in this regime (with b and h representing the extremes). This is due to both the difference in collisional cross sections and positions (since tacc ∝ a 2 pla , see Fig. 5). In subsection 5.1.2, we assess the outcome for the specific values ofq that have been derived in Sec. 4 for the different scenarios.

Forq derived from specific scenarios
In Fig. 12, we see that for a planet scattering scenario (both for a single planet or a chain) in which Q = 1au andq = 10 −5 au/yr (i.e. corresponding to a 1M⊕ planet in Sec. 4.1), we end up in the regime where a fraction of comets is accreted onto each planet. The fraction accreted is rather high in this case (between 0.01 and 0.03 for planets b to g) because the probability to be accreted for comets with Q = 1au is rather high (as expected from Fig. 3). This fraction accreted is valid for a single eccentric planet scattering material from a close-in belt at ∼ 1au similar to the debated belt potentially found around Proxima Cen (Anglada et al. 2017). However, comets coming from tens of au belts would have to be scattered through a planet chain before making it to the innermost planet of the chain and some will be lost on the way. Here we assume that fin ∼ 5% of the comets will make it to the innermost planet (see Sec. 4.1), which reduces the fraction accreted on the different planets from the initial reservoir (the Kuiper-belt like disc) to ∼ 5×10 −4 .
For the Kozai scenario, we consider that Q represents the disc location from which the comets are perturbed by an outer companion and we take a typical distance of 100au as being representative. Considering the typicalqK value derived (2 × 10 −4 au/yr), Fig. 12 shows that we are in the second regime where a fraction of comets (∼ 10 −5 ) is accreted onto each of the seven planets. This is close to two orders of magnitude smaller than the chain of planets case. Here, we do not have to reduce the number of comets that arrives on to the seven planets (i.e fin = 1) as Kozai oscillations operate directly from the outer belt.
For the Galactic tide scenario, from Fig. 12, we evaluate that an Oort cloud at a few 10 4 au (that has a fast enoughqG to send comets to small pericentres within a fraction of the age of the system), i.e with large apocentres, the probability to be accreted is always 10 −6 (for all the plottedq, i.e more than order of magnitude smaller than for the Kozai mechanism). Therefore, the fraction of of the accreted comet's material will stay in the atmosphere as we will show later. comets accreted is very low, which will not have any impacts on the atmospheres. Therefore, we rule out Galactic tides as being an efficient mechanism 6 to modify the atmospheres of the TRAPPIST-1 planets. We also note that the same forces driving particles with high pericentres to low pericentres could also drive back these low pericentre orbits to high values, sometimes before they had time to reach the 7 inner planets (e.g. Emel'Yanenko et al. 2007;Rickman et al. 2008), and thus makes this scenario even more unlikely. We are, thus, left with two plausible mechanisms to throw comets on the seven planets, namely, scattering by planets and Kozai oscillations due to an outer companion. 6 We note that we could fine tune the position of the Oort cloud to be in a narrow range in between 10 4 and 10 5 au to maximise the fraction accreted while allowing enough time for the comets to reach the planetary system but this would always result in an order of magnitude less efficient mechanism than Kozai. Moreover, it is not likely that an Oort cloud around a low mass star such as TRAPPIST-1 forms farther out than in our Solar System ) as required here for maximising Galactic tide effects. And as shown in Sec. 4.3, such distant belts should be depleted owing to passing stars.

The relative effect of different impactor sizes on the atmospheres of the different planets
In the previous subsection, we analysed the fraction of comets accreted facc by each planet. However, we want to quantify the effect of these impacts on the atmospheres of the seven planets. For example, we have seen in Fig. 7 that impact velocities are much higher for planet d compared to further planets, so even if the fraction accreted is the same as that of the more distant planets in the planet scattering scenario, the effect on atmospheric mass loss may still be more important. Here, we quantify the atmospheric mass loss, and projectile mass accreted in the atmosphere (relative to impactor masses), i.e., the impactor mass that does not escape the atmosphere after impact, for the different planets and for different impactor sizes.
We use the numerical study of the effect of impacts on atmospheres by Shuvalov (2009) to derive some conclusions for the TRAPPIST-1 planets. We first present the set of equations from Shuvalov (2009) that we use to derive the atmospheric mass loss and projectile mass accreted after a given impact. The outcome depends on the dimensionless variable η (Shuvalov 2009)  It shows how much atmospheric mass is lost after a given comet is thrown in, taking into account that the fraction of comets that hit the different planets is not equal to 1 as already seen in Fig. 12. Right: Accreted projectile mass to impactor mass ratio (M impacc /M imp )facc as a function of D. It shows how much projectile mass is accreted after a given comet is thrown in. The different lines are for Q = 1 (solid lines) and 100au (dashed lines) withq = 10 −5 au/yr and for Q = 1au andq = 10 −1 au/yr (dotted line).
where D is the impactor diameter, H the atmosphere scale height (H = kT /(µmH g) for an isothermal atmosphere with g = GMp/R 2 p ) and ρt, ρpr, ρatm0 are the densities of the target (planet), projectile (exocomet), and atmosphere at the surface, respectively. We assumed ρt = 5000 (terrestrial planet-like), ρpr = 1200 (comet-like), ρatm0 = 1.2kg/m 3 , µ = 28.97 (we assume an Earth-like atmosphere for now) and T is taken to be the equilibrium temperature of the planets (assuming a null Bond albedo as calculated in Gillon et al. 2016). We note that recent observations suggest that some of the TRAPPIST-1 planet densities may be slighlty lower because of the potential presence of ice layers. Grimm et al. (2018) find that water mass fractions < 5% can largely explain the observed mass-radius relationship of the less dense planets. Therefore, we can expect densities that are 10s of percent lower than assumed here, which would translate as a small uncertainty on η, which is however much lower than the uncertainties on the dynamics (see Sec. 4), and is thus not considered here in details. Vimp is the impact velocity and Vesc = 2GM pla /R pla is the escape velocity for the different planets. We have seen in Sec. 3.4 that Vimp is much greater than Vesc, which simplifies the previous and following equations for most cases.
To get meaningful results, we compare the atmospheric mass loss M atmloss to the impactor mass (of size D) that makes it to the inner regions and that is accreted on to the planet. Therefore, using the previous notations, we are interested in (M atmloss /Mimp)finfacc where we recall that fin is the proportion of comets that are scattered from an outer belt and make it to the inner regions and facc is the accreted fraction onto a given planet. This ratio can therefore be understood as the atmospheric mass that is removed by one comet scattered from an outer belt, where only "a fraction" of the comet makes it to the inner regions and a fraction of that is accreted onto a specific planet. In Fig. 13 (left), we plot (M atmloss /Mimp)facc keeping in mind that we should multiply that value by fin (if it is different than 1, see Table 2) to get the real value of accreted comets that make it to the inner regions. Fig. 13 (left) shows (M atmloss /Mimp)facc as a function of impactor diameter D for Q = 1 and 100au (anḋ q ∼ 10 −5 au/yr) and for a higherq (10 −1 au/yr) and Q = 1au, using the impact velocity distributions shown in Fig. 6. The overall shape of the curves in Fig. 13 (left) is explained in Shuvalov (2009). Impactors of size a few kms are the most harmful at removing atmospheric mass. Impactors smaller than 100m do not create large impact plumes and cannot accelerate large atmospheric mass to high latitudes. For impactors larger than a few 10s of kms, atmospheric erosion continues to grow very slowly but the mass an impact removes cannot be greater than the total local atmospheric mass available. Therefore, for large impactors the atmospheric mass removed per increasing impactor mass becomes smaller.
The most harmful impactor size shifts along the x-axis for the different planets mainly because of the change in impactor velocity and the different properties of the planets through H (the atmosphere scale height) with the relative scalings given in Eq. 9. The variations along the y-axis are mainly due to the different fraction accreted facc for each different planet (see Fig. 12) and the different impact velocities (see Fig. 6) and scale as shown in Eq. 10. For example, we see that even though planets d, e, f, g accrete at the same level (see Fig. 12), the atmospheric mass loss is greater for planet d because impact velocities are higher for the closer in planets (see Fig. 6).
The effect of increasing Q from 1 to 100au (solid to dashed lines) is to shift all the lines down by a factor 100 because facc decreases by a factor 100. Changingq from 10 −5 (solid) to 10 −1 au/yr produces a shift downwards of four orders of magnitude since facc decreases by a factor 10 4 between these two cases. Fig. 13 (left) can therefore be used to work out the relative effectiveness of comets at removing mass from the atmosphere of each planet for any givenq and Q, even though we show the results for only two different Q (i.e., it is a general plot, not tied to a specific scenario from Sec. 4, and only facc ∝q −1 Q −1 changes for different values ofq and Q, making it easy to compute results for differenṫ q and Q).
The simulations of Shuvalov (2009) also showed that the projectile mass accreted per impactor is given by where χ pr = min{1, 0.07(ρt/ρpr)(Vimp/Vesc)(log 10 η − 1)}. Similarly to atmospheric mass loss, Fig. 13 (right) shows the accreted projectile mass per comet (Mimpacc/Mimp)facc as a function of impactor diameter D for Q = 1 and 100au (andq ∼ 10 −5 au/yr) and for a higherq (10 −1 au/yr) and Q = 1au. The shape of the curves in Fig. 13 (right) is already known from Shuvalov (2009). The ejecta from impacting bodies that are 1km does not have enough energy to escape after impact and is stranded in the atmosphere (though some material may condense on the planet surface at a later point, see Sec. 6.3). For more massive bodies, the ejecta after impact is increasingly more energetic until the airless limit is reached (i.e., when atmospheric drag can be neglected before the after-impact plume expansion) where all the projectile material escapes. This cut-off happens for bodies larger than a few km.
In Fig. 13 (right), the variations along the x-axis (e.g. of the cut-off position) are due to different impact velocities (for instance a larger planetesimal can deliver material onto planet h because impacts happen at lower velocities) and it can also vary with the planets' properties through H and the atmospheric density (assumed constant for now) with the scalings given by Eq. 9. Planets g and h can therefore get volatiles delivered from larger comets than further in planets. The variations along the y-axis are mainly due to the fraction of comets accreted onto the planets and the different impact velocities and scale as depicted by Eq. 11.
The effect of increasing Q from 1 to 100au (solid to dashed lines) or increasingq from 10 −5 to 10 −1 au/yr is the same as explained when describing Fig. 13 (left), i.e., due to the change in facc. This plot is therefore also general and can be used to compute the outcome of an impact for any values ofq and Q, and is not tied to any of the specific scenarios explained in Sec. 4.
The volatile mass that ends up in the atmospheres is a fraction f vol of the mass delivered. We assume that volatiles are delivered to the atmospheres in proportion to their fraction of the mass of the parent body. For a comet-like body, we assume a rock-to-ice mass ratio of 4 based on recent measurements in the 67P comet (Rotundi et al. 2015), i.e 20% of ice by mass. For an asteroid-like body, the water mass fraction is lower and is found to vary between 10 −3 and 0.1 (Abe et al. 2000). We will assume an intermediate value of 1% for asteroid-like bodies 7 , which is typical of ordinary chondrites in our Solar System (but we note that carbonaceous chondrites can reach 10% of water by mass, Raymond et al. 2004). This gives us two extreme volatile delivery scenarios to consider with our model.
The CO or H2O content of exocomets can be probed for the most massive belts and are found to be similar to Solar System comets (e.g. Kral et al. 2016;Marino et al. 2016;Matrà et al. 2017a). The potential to detect gas in debris disc systems will improve with new missions (see Kral et al. 2017b) and the assumptions used in this study could be refined with future estimates of the volatile content of exocomets in the TRAPPIST-1 system to get a better handle on the final atmospheric composition.

5.3
The integrated effect of these impacts over the age of the system

Total incoming mass over the system's age
We now work out the effect of impacts on the TRAPPIST-1 planets over the age of the system and more specifically, how much atmospheric mass is lost and how much projectile/volatile mass is accreted for a given total incoming mass of comets. To do so, we assume a typical N (D) ∝ D γ size distribution with γ = −3.5 for the comets that are expelled from the belt (e.g. Dohnanyi 1969;Thébault & Augereau 2007) up to a maximum size of 10km 8 . Indeed, integrating over the assumed size distribution for the total atmospheric mass loss (or accreted material) shows that > 10km impactors are unimportant (as already concluded by Schlichting et al. 2015) because M atmloss ∝ D −2 for large bodies as seen from Fig. 13 (left), which decreases faster than the gain in mass of these larger bodies (∝ D 0.5 ). Very massive giant impacts (e.g. Kral et al. 2015) of bodies with radius > 1000km (i.e Pluto-sized or greater) can have a devastating effect on the atmosphere of a planet (Schlichting et al. 2015), which is not modelled in Shuvalov (2009), but these impacts are rare and thus neglected in this study. We consider an incoming mass of comets Minc that reaches and can potentially hit the TRAPPIST-1 planets after a mass Msca of comets has been scattered from this outer belt over the system's age. Taking into account the efficiency to reach inner regions, Minc = Mscafin (see Fig. 14).
The integrated amount of mass scattered from a belt Msca over the system's age can be evaluated. In Sec. 2, we predicted that a potential planetesimal belt of 20M⊕ could potentially have survived around TRAPPIST-1 at tens of 7 We assume that the bulk of the volatile mass is in water so that this value is representative of the total volatile mass, though a lower limit. 8 We note that the size distribution of the Kuiper belt for the largest bodies is complicated and best-fitted by two shallow power laws and a knee or a divot between the two (Lawler et al. 2018), which would imply that most of the cross section would be in the biggest bodies. This is not representative of what is observed in general for the debris disc population, for which a -3.5 slope all the way through the largest bodies is able to explain the observations. Table 2. Table describing the parameters used for the different scenarios we tested. We list the rate of change of pericentreq, the apocentre Q, the mean fraction accreted facc on each planet, the fraction of comets that makes it to the inner regions f in , the mass fraction of volatiles on the exocomets/exoasteroids f vol , the minimum scattered mass M scadestroy to destroy all 7 primordial atmospheres (Msca = M inc /f in ), the mass of delivered volatiles M volmin and water M watmin (assuming Solar-System comet-like compositions) for a belt scattering at the low scattering rate of the current Kuiper belt (i.e M inc ∼ 10 −2 M ⊕ f in ) for each of the planets f, g, h, and M volT , M watT for a belt of 20M ⊕ (close to the expected mass for a potential leftover belt around TRAPPIST-1, see section 2) scattering 5% (i.e M inc ∼ 1M ⊕ f in ) of its mass over 7Gyr. For the case of exoasteroids, f vol = 0.01, and M volmin as well as M volT should be divided by 20, and M watmin , M watT by 10. Meo means 1 Earth ocean (i.e 2.5 × 10 −4 M ⊕ ).
Scenariosq au. By the action of a nearby planet, many planetesimals may have been scattered inwards over the lifetime of the system. Assuming that 5% of the belt mass is scattered over 7 Gyrs (using results by Marino et al. 2018), we get that Msca ∼ 1M⊕ leading to Minc ∼ fin M⊕. In our Solar System, ∼ 0.27 comet/yr leave the Kuiper belt towards the inner regions . The typical mass of comets in 's study is ∼ 4 × 10 13 kg so that the rate of scattered incoming comets isṀsca ∼ 2 × 10 −3 M⊕/Gyr. Therefore, a similar Kuiper belt around TRAPPIST-1 would give Msca ∼ 10 −2 M⊕ over 7Gyr leading to Minc ∼ 10 −2 fin M⊕.
However, the Kuiper belt is thought to have been a lot more massive in its youth (e.g. Levison et al. 2011, and see Sec. 2) and in general, debris discs that are observed can have fractional luminosities of up to 10 4 greater than this low-mass belt (Wyatt 2008), which is an indicator of them being more massive. We note that the Kuiper belt is so light (∼ 0.1M⊕, Fraser & Kavelaars 2009;Vitense et al. 2010) that current instruments could not even detect it around another star (Vitense et al. 2012;Kral et al. 2017b). From an MMSN-like calculation, the initial Kuiper belt mass may have been of several 10s of Earth masses (Hayashi 1981;Morbidelli et al. 2003), meaning that Msca could have been of the order of a few 10M⊕ owing to the depletion of the belt to reach its current mass. In other words, we expect 10 −2 fin M⊕ Minc 30fin M⊕. (12)

Atmospheric mass loss
The total atmospheric mass loss for a given planet over the system's age is where N (D) is the number of bodies in each impactor diameter bin D that make it to the inner regions. Fig. 15 shows M totatmloss /Minc, i.e the total atmospheric mass loss compared to the incoming mass Minc of comets injected into the inner regions over the lifetime of the star. Once again, this figure is general (and can be used for anyq and Q) and is not tied to a specific scenario (only the black vertical lines are scenario dependent). We show the atmospheric mass removed for specific values of apocentres Q = 1 and 100au but values for other Q can also be estimated (as M totatmloss ∝ facc ∝ Q −1 ). Atmospheric mass loss remains lower for planets g and h because impacts happen at lower velocities (see also Fig. 13 left). The mean total atmospheric mass loss for the seven planets can be approximated as where we note that this ratio is accurate for planets d, e and f but can be a factor 10 more or less for a specific planet (e.g. 10 times higher for planet b and 10 times lower for planet h), and Fig. 15 should be used to get more accurate values.
To assess whether the impact process is capable of destroying an entire primordial atmosphere, we first estimate the primordial atmospheric masses of the different planets. These primordial atmopsheric masses are not known and so for reference we asssume an Earth-like composition and density. Computing the scale height for each planet (as in Eq. 9) and assuming an isothermal atmosphere of temperature T (the equilibrium temperature of the planets), we integrate over the height of the planet atmospheres to get their masses Matm = 4πρatm0H(R 2 pla +2HR pla +2H 2 ). This gives primordial atmospheric masses of 2, 0.9, 0.7, 0.9, 1.1, 0.7, 0.4 × 10 −6 M⊕ for planets b to h. This is shown on Fig. 15 as horizontal lines, where this mass has been divided by 1M⊕ to show the effect of an incoming mass of 1M⊕.
Therefore, a primordial Earth-like density atmosphere on the TRAPPIST-1 planets could be destroyed if 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 q (au/yr) Minc > 5 × 10 −4 q 10 −5 au/yr For the specific physical scenarios from Sec. 4 (see vertical black lines on Fig. 15 for the planet scattering and Kozai scenarios), Table 2 shows the minimum scattered mass needed M scadestroy from an outer belt (the minimum incoming mass would be finM scadestroy ) to destroy the primordial atmospheres of the seven planets. For example, for the planet scattering scenario with a single Earth-mass planet at 1au (i.eq ∼ 10 −5 au/yr, Q = 1au and fin = 1, see Table 2), using Eq. 15 we find that Minc 5 × 10 −4 M⊕ can destroy the primordial atmospheres of the seven planets. This corresponds to a belt that is being depleted for 7Gyr at a rate ten times lower than that at which the current Kuiper belt is being depleted. If the comets had to be passed in through a planetary system before reaching the planet at 1au, the inefficiency in the inward scattering process results in an additional factor fin = 0.05. This means that even with this factor, the current Kuiper belt scattering rate is enough to destroy the atmospheres of the seven TRAPPIST-1 planets.
For the Kozai scenarioq values are higher (q ∼ 2 × 10 −4 au/yr) and Q is at larger distances (100au), meaning that interactions with planets are much more likely to result in ejections rather than accretions (see Fig. 3). We find that Minc > 1M⊕ is needed to destroy the primordial atmospheres, i.e two orders of magnitude larger than in the planet chain case. For a 1M⊕ incoming mass (i.e. 100 times the current Kuiper-belt like incoming mass rate), Fig. 15 shows that the atmospheric mass loss is ∼ 2 × 10 −5 M⊕ for planet b and a factor 10 less for planets c, d, e, and f, and about another factor 5-10 less for planets g, h (all of which are higher than the primordial Earth-like atmospheric masses assumed here except for planets g and h that are a factor 2 too small).
Given that the exo-Kuiper belts detected around F, G, K stars are much more massive than the Kuiper belt, and that the possible belt mass we derive for the TRAPPIST-1 belt in Sec. 2 is ∼ 20M⊕), the scattering may be even higher than assumed here (i.e., up to a factor of a few 10 3 the Kuiper-belt incoming mass), and we conclude that if a scattering belt is around TRAPPIST-1, the primordial atmospheres would not survive impacts over the system's lifetime for both a planet scattering and a Kozai scenario.

Water mass loss
In Table 3, we also quantify the maximum water mass loss for the single and planet chain scenarios. The water mass loss MwatLossT is given for each planet for a belt of 20M⊕, which is close to the expected mass for a potential leftover belt around TRAPPIST-1 (see section 2) scattering 5% (i.e., Minc ∼ 1M⊕ fin) of its mass over 7Gyr. For the planet chain scenario, the planets can lose up to 4, 1.2, 0.8, 0.6, 0.4, 0.12, 0.06 Meo (Earth ocean mass), for b, c, d, e, f, g, h, respectively, and 20 times more for the single planet case. These values can be compared to the water mass loss from hydrodynamic escape due to XUV irradiation during the runaway greenhouse phase, for which they found upper limits of (Bourrier et al. 2017), 80, 40, 2.9, 1.5, 0.9, 0.4, 0.1 Meo, for b, c, d, e, f, g, h, respectively. These values are, however, to be taken as strict upper limits because it is uncertain that hydrogen can reach the very top layers at the base of the hydrodynamic wind, which is needed for it to escape (Bolmont et al. 2017). Also, this hydrodynamic escape works well to eject hydrogen but other atoms are difficult to drag along (Bolmont et al. 2017). For the impact case, not only hydrogen would escape but the whole fluid in the ejected plume. Bearing these caveats in mind, we can now compare the water mass loss from hydrodynamic escape to the impact scenario. For the planet chain case, the water mass loss due to impacts seems to be less efficient than hydrodynamic escape for planets b and c and both scenarios are within a factor of a few for the other planets. For the most optimistic case of the single planet case, impacts could produce the same water loss as hydrodynamic escape for planets b and c and be an order of magnitude higher for planets d to h.

Delivery of volatiles
We now evaluate the total mass of material and volatiles that can be delivered from the impactors over the system's lifetime. We derive the total accreted projectile mass Mtotimpacc by integrating the mass accreted per impactor (Fig. 13 right) over the assumed size distribution. This accreted mass is assumed to be deposited in the planets' atmospheres. Fig. 16 shows Mtotimpacc/Minc, the total accreted projectile mass compared to mass of comets injected into the inner regions over the lifetime of the star. The overall shape is similar to Fig. 15, but note that planet h is far better for delivery of mass into its atmosphere than having its atmosphere depleted because impacts are at lower velocities (which means material from larger planetesimals can be accreted, see Fig. 13 right). This means that the mass delivered on planet h (and planets with similar impact velocities) may be greater than that lost after each impact.
To quantify this, Fig. 17 shows the ratio of the accreted projectile mass and atmospheric mass lost, which does not significantly depend onq, instead only depending on the size distribution of comets and slightly on Q. Thus, this ratio is plotted as a function of the slope in the size distribution γ for two different values of Q (1 and 100au). For γ = −3.5, this ratio is greater than one for planets g and h and close to 1 for planet f but lower for the other planets. This means that even if all of the accreted mass ends up in the atmosphere, the total atmospheric mass must be decreasing for the inner planets and can only increase for the outer three planets if all mass ends up in the atmospheres. Regardless, all planets will have their atmospheres enriched by the planetesimals' composition and the situation is similar for all Q.
Consider now the fraction of this delivered projectile mass that will be in volatiles, i.e M totvolacc = Mtotimpaccf vol , which could be delivered to planets from comets. To assess the amount of volatiles that are delivered to the planets we consider two types of material that impact on to these planets (presented in Sec. 5.2); 1) Cometary-like material with 20% of ice by mass (f vol = 0.2), and 2) Asteroid-like material with ∼1% of volatiles by mass (f vol = 0.01).
One important question is whether the icy material will have dissappeared through sublimation before impacting the planets. Marboeuf et al. (2016) show that a 1km comet survives sublimation for ∼ 560 orbits around a 0.1L star. TRAPPIST-1 is 200 times less luminous and so the comets will survive much longer. Extrapolating Marboeuf et al. (2016)'s formula to TRAPPIST-1 luminosity, we get that a 1km comet passing at small pericentres (0.1au) would need more than 10 5 orbits to sublimate. As it only takes a few 100s of orbits for the comets to be accreted on the planets (see Fig. 4), we assume that most of the icy content of the comets will not have sublimated and so will be available to be delivered at impact. We note that during theq evolution, the sublimation will start happening only in the very last phase, i.e when the pericentre is close to the planets already (because for larger pericentres the mass loss from comet sublimation is very slow and the timescale of evolution of the pericentre is much faster, Marboeuf et al. 2016). Thus, the impact timescale of 100s of orbits is a good indicator of the number of orbits before impact during which sublimation could happen. Therefore, we estimate the mean of the total volatile mass delivered on each of the seven planets as which is, for all planets, within a factor 3 of that from Fig. 16. For all planets, we can also estimate the incoming mass needed to deliver more volatiles than the primordial atmospheric mass where we assumed primordial atmospheres of Earth-like densities. Thus, for the planet scattering scenario, we find that only a small incoming mass is needed to deliver enough volatiles to potentially replenish an atmosphere with an Earth-like density (e.g., Minc > 3 × 10 −3 M⊕ for comet-like bodies scattered from an outer belt to a 1M⊕ planet at 1 au). The incoming mass needed for the Kozai scenario is larger, 5M⊕, but not implausible to reach as shown by Eq. 12.
From Fig. 17, we have shown that only planets g and h (and possibly f) would be able to retain the largest part of the delivered volatiles. This means that for the planetscattering and Kozai scenarios, the new atmospheric compositions of planets f, g and h could be entirely set by the comet volatile content, which would replenish the atmospheres over the system's age. However, the absolute level of the volatile Table 3. Amount of water lost due to impacts for the planet scenario (single and chain). The water mass loss M watLossT is given for each planet for a belt of 20M ⊕ (close to the expected mass for a potential leftover belt around TRAPPIST-1, see section 2) scattering 5% (i.e M inc ∼ 1M ⊕ f in ) of its mass over 7Gyr. Meo means 1 Earth ocean (i.e 2.5 × 10 −4 M ⊕ ).
Single planet (10M ⊕ ) 5 × 10 − content that will remain in the atmosphere is difficult to constrain as some fraction of the volatile mass will be ejected by later impacts or end up on the planet's surface and some other sources of volatiles could be present (see Sec. 6.2). We can, however, estimate the amount of volatiles that will survive after each impact assuming that a fraction fr of the accreted material remains in the atmosphere rather than condensing on the planet. Therefore, after a given impact frMtotimpacc of material will be added to the atmosphere and the next impact could remove a maximum of M totatmloss from this added material. Assuming that fr = 1 (if impacts are frequent enough, e.g. LHB-like, material does not have time to condense back on the surface), we compute the fraction of volatiles that would accumulate from subsequent impacts in Fig. 18. We note that some additional volatiles could be added by degassing of the planets' interiors but that fr may also be smaller so that the exact volatile mass that can accumulate depends on complex physics that cannot be modelled in this paper. We see that indeed, only planets f, g, and h have positive values (i.e. they gain volatiles over time) and therefore appear 9 in Fig. 18 showing M vol /M incvol , where M incvol = Mincf vol is the incoming mass of volatiles. We can also derive a general formula as a function ofq and Q (similar to Eq. 16) that gives the mass of volatiles that can accumulate M vol rather than the total volatile mass delivered. We do that in Sec. 6.1.2 and give the temporal evolution (assuming a constant rate of impact) of the build up of the secondary atmospheres that are created for planets f, g, and h (see Eq. 21). We thus conclude that the atmospheres of planets f, g for planets g and h (and for fr < 0.8 for planet f) so that no secondary atmospheres would accumulate in this case, but this neglects outgassing which would add more volatiles and would make it harder to not build up secondary atmospheres on these three planets.
and h might be more massive than that of the innermost planets of the TRAPPIST-1 system if cometary bombardement has happened, and that a fraction of their composition should reflect the cometary abundances in this system. We note that the build-up of secondary atmospheres for planets f, g and h is mainly allowed by the impact velocities that are low enough on these outermost planets to both reduce the atmospheric mass loss after each impact and allow to deliver more volatiles (from larger bodies).

Delivery of water
Water on Solar System comets makes up more than fwat = 50% of the volatiles (Mumma & Charnley 2011). Depending on fwat for exocomets, the amount of water Mwater delivered on the seven planets can be approximated by Minc, (18) where f vol ∼ 0.2 for exocomets (∼0.01 for asteroids) and fwat ∼ 0.5 (∼1 for asteroids). For example, for the single Earth-mass planet scattering scenario (i.eq ∼ 10 −5 au/yr, Q = 1au, and fin = 1), we find that a belt scattering at the same low rate as the current Kuiper-belt would result in the planets accreting ∼ 8 × 10 −3 Earth oceans of water (or 10 times less for asteroid-like bodies), assuming that one Earth ocean equals 1.5×10 21 kg (see Mwatmin in Table 2). We note that for the planet chain case (where fin = 0.05), these values would be a factor 20 smaller and for a larger incoming mass Minc these values could go up by a factor more than 10 3 (see Eq. 12). We find that a belt of 20M⊕ (similar to the plausible belt mass we predict around TRAPPIST-1 in Sec. 2) that would scatter 5% of its mass over 7Gyr (i.e., Minc ∼ 1M⊕ fin) would deliver ∼ 1 Earth ocean of water to the planets for the single planet case and ∼ 0.04 Earth ocean for a planet chain (see MwatT in Table 2).
For the Kozai scenario, we find that between ∼ 10 −5 (pessimistic case with a Kuiper belt scattering rate) and ∼ 0.01 (optimist case with Minc ∼ 20M⊕) Earth oceans of water could be delivered to the planets.
This delivered water will presumably recondense as ice on the surface of planet h (but when the star was younger this planet was in the HZ and water could have been in liquid form for a long period, see Sec. 6.1.2), but for warmer planets such as planets f and g, we expect that a rain cycle would create liquid water on these planets that would then be reinjected into the atmospheres cyclically (see Sec. 6.3). The temporal evolution of the build up of the amount of water in these secondary atmospheres can be obtained from Fig. 18 or from the coming Eq. 21 for planets f, g, and h.

Comparison between timescales of the different processes
The consideration of timescales is important because it constrains the duration over which atmosphere loss/gain occurs 10 5 10 6 10 7 10 8 10 9 Time (yr) compared with other processes which may be taking place, but which are beyond the scope of this manuscript to consider in detail.

Timescale to lose primordial atmospheres from impacts
Assuming a constant rate of scatteringṀsca over 7Gyr, we compute the atmospheric mass lost as a function of time M atmlossc (t) ∼ 2 × 10 −3 fin q 10 −5 au/yr where we note thatṀsca = 0.1M⊕/Gyr corresponds to a belt with a total incoming mass of ∼1M⊕ fin, i.e similar to what would be expected for a 20M⊕ belt scattering 5% of its material over the age of the star. M atmlossc (t) becomes greater than an atmospheric mass of 10 −6 M⊕ for Now, we consider the planet chain scenario (i.e witḣ q = 10 −5 au/yr, Q = 1au, and fin = 0.05) with a scattering rateṀsca = 0.1M⊕/Gyr and look at the temporal evolution of the atmospheric mass loss M atmlossc (t) due to the series of impacts over the system's age as shown by Fig. 19. By comparing to the primordial atmospheric masses of the planets (horizontal lines in Fig. 19), we see that for this scenario, it takes between 10 and 400Myr to destroy the primordial atmospheres of all seven planets (assuming an Earth-like atmospheric density). This is very fast compared to the age of the system. This shows that the timescales over which the primordial atmospheres can be destroyed are much shorter than the age of the system. Therefore, we confirm the previous conclusion (see Sec. 5.3.2) that cometary impacts may have entirely stripped all planets of their primordial atmospheres by 7Gyr, even if the scattering rate is smaller by a factor more than 10 than assumed here (i.e., close to the Kuiper-belt scattering rate level).
6.1.2 Timescale to regenerate secondary atmospheres from impacts for planets f, g, and h We also compute the temporal evolution of the volatiles M vol that are deposited and accumulate after each impact (i.e we take into account that subsequent impacts remove part of the volatiles delivered by the preceding impact as in Fig. 18). For planets g and h, M vol is given by (for anyq and Q) (21) and is a factor 5 smaller for planet f. The amount of water delivered is simply Mwater = fwatM vol . Now, we work out the timescale to replenish the secondary atmospheres of planets g and h in cometary volatiles at the level of a 10 −6 M⊕ atmospheric mass and a factor 5 longer for planet f. The replenishment timescale shows that in most physically motivated cases planets f, g, and h will have had time (over the age of the system) to rebuild secondary atmospheres with masses of at least 10 −6 M⊕, i.e., equal or greater than an Earth-like primordial atmosphere. We note that most of the volatiles delivered by the comets have low condensation temperatures and thus would remain in the atmosphere rather than go on the planet's surface but water could condense as ice on planet h and cycle from the surface to the atmosphere on planets f, g owing to rain (see Sec. 6.3). Therefore, we expect M vol to be a good estimate of the amounts of volatiles that can accumulate for planets f and g and note that up to 50% of the volatiles (to account for water) could transform into ice on planet g and thus reduce M vol by a factor 2 (but this ice could outgas at a later stage because of the planet activity). show that for a 0.08M star, the HZ location moves inwards to its present-day position after ∼ 1Gyr. This means that planet h will be the first to enter the liquid water HZ, which it will do at a point when the closer-in planets are still in a runaway greenhouse state (assuming they have retained any atmospheres). According to the Luger & Barnes (2015) model, planet h crosses into the empirical habitable zone at ∼30Myr. Coupled with our results, this scenario indicates that planet h could have received significant volatile delivery at a point in its history (i.e., between 30Myr and 1Gyr) when liquid water was stable at its surface (Fig. 20). This raises the prospect for an early carbon cycle being established on this planet, stabilising climate through water-rock interaction as is inferred for Earth (Walker et al. 1981).
6.2 Additional sources of volatiles 6.2.1 Volatiles created by vapourised material from the planet's surface during impact The volatile fraction that ends up in the atmospheres of the TRAPPIST-1 planets does not only build up from the impactor material but also from the vapourised material from the planet surface, as was probably the case for the Chicxulub impact that may have released large quantities of gas and dust contributing to the environmental stress that led to the demise of dinosaurs on Earth (Pope et al. 1997). From Okeefe & Ahrens (1977, we can estimate the volume of material Vvap vapourised from a given meteoritic impact (with a volume Vpr). They find that Vvap = 0.4SVpr, where S = (ρpr/ρt)(Vimp/Cp) 2 , using the same notations as in previous sections and Cp being the bulk sound speed of planetary surface, which varies depending on the planet ground composition (Melosh 1989). We assume an Earthlike composition for which Cp ∼ 7km/s. We thus find that the vapour mass Mvap produced for a given impactor of mass However, some of the vapour ejecta will escape and only a fraction will have a low enough velocity to be retained in the atmosphere. Once again, using results from Shuvalov (2009), we get that the maximum ejected fraction of target material after impact is Mtaresc ∼ 0.02Mimp(Vimp/Vesc) 2 . This maximum is reached for bodies that are larger than ∼1km and for smaller bodies, the planet retains almost all of the target material created at impact. Of course, above a certain threshold it means that the whole target mass escapes (as Mtaresc becomes greater than the total atmospheric mass), which is similar to the projectile mass behaviour (where volatiles from bodies larger than ∼10km cannot be retained in the atmosphere). We thus find that Mvap is a good indicator of the vapourised mass that will remain in the atmosphere (as Mtaresc Mvap). From Eq. 11, we notice that the mass delivered from the projectile quickly tends to Mimp for bodies smaller than about 10km. Thus, for planets g and h that have median impact velocities of ∼ 25 and 20km/s, Mvap will be slightly higher but of the same order of magnitude as Mimpacc. This means that some volatiles such as SO2, CO2 or water could also be formed from the vapourised planets' crust (see Pope et al. 1997). However, we note that the typically low concentration of volatiles in planetary basalts that would form the bulk of a crust would not release as many volatiles as for the Chicxulub impact (e.g. Dreibus & Wanke 1987;Saal et al. 2002).

Outgassing on the planets
Degassing may happen early during accretion when forming the planets but this is not a concern in our study as we expect the primordial atmospheres to be totally destroyed. Degassing from tectonic activity may also happen at a later stage that could affect the amount of volatiles in the atmospheres. Another way of producing degassing is from stellar induction heating. A recent study that focused on the effect of this mechanism on the TRAPPIST-1 planets finds that induction heating could create strong outgassing on planets b, c, d that are very close to their host star but it should not affect the outermost planets e, f, g, and h (Kislyakova et al. 2017).
For the plate-tectonic degassing, we take the degassing on Earth as an upper bound because plate tectonics is very active on Earth and may be less efficient/active on other planets 10 . Earth produces ∼ 22km 3 of basaltic magmas each year (Crisp 1984). Given a magma density of 2600kg/m 3 , we estimate a total degassing rate of ∼ 6×10 13 kg/yr. Assuming a typical water content of 0.3wt% and the extreme case of perfectly efficient degassing with no subduction recycling of water to the planet's mantle, we find that an upper bound on the tectonically driven water degassing rate is ∼ 3 × 10 −5 M⊕/Gyr (0.11 Earth oceans per Gyr). Therefore, if the tectonic activity on planets f, g and h were as active as on Earth, degassing of water could occur at a similar rate to the water delivered from impacting comets (see Table 2), thus enhancing the amount of water on planets f, g and h.

Volatiles that are ejected of the atmosphere and reaccreted later
The material that escapes the planetary atmospheres after each impact because they have velocities greater than the escape velocity will end up in an eccentric torus around the star close to the given planet location (e.g. Jackson et al. 2014;Cataldi et al. 2017). The eccentricity will vary depending on the ejection velocity of the material. While we expect that high-velocity ejecta may reach neighbouring planets (e.g. in a Panspermia-manner, Krijt et al. 2017;Lingam & Loeb 2017), most of the material in the torus would interact with the planet it has been ejected from. We note that for an Earth-like planet on a slightly wider orbit than planet h, the escape velocity (of about 10km/s) could become greater than the planet's Keplerian velocity and thus the material would not form a torus but rather be ejected on unbound orbits. The fate of the material in the torus is not straightforward to model. The material could deplete collisionally due to high-velocity collisions in the elliptic torus and be ground down to dust, which would be blown out from the system by stellar wind radiation pressure (Wyatt 2008) and at the same time eject the ices or volatiles present on the grains. While the ejecta is also partly made up of gas, one could also expect that the gas material (at least the fraction that is not blown out by radiation pressure) in the torus will viscously spread (maybe dragging dust with it) and end up on more distant planets. The fate of the material that would be able to interact with a planet for a long enough timescale is to be reaccreted onto the progenitor planet . The exact outcome depends on the exact chemico-physical conditions in the TRAPPIST-1 planets environment, which is not known, and thus goes beyond the scope of this paper.

Composition of the atmospheres at the end of the impact process
Thanks to our model, we are able to retrieve the amount of volatiles that is delivered to the different planets as well as the atmospheric mass removed by a long-term series of impacts. For the outermost planets, we find that the volatiles delivered by impacts may accumulate and be abundant, which could give us a way to constrain the atmospheric composition of planets f, g, h, the former two being in the HZ. However, we need to understand how these delivered volatiles would evolve in their new atmospheres to predict the current atmospheric compositions of these planets. They could chemically react to form new species, condense on the surface as ice and some additional volatiles may be produced as seen in the previous Sec. 6.2. For instance, the delivered water will presumably condense as ice on the surface of the colder planet h (when it finishes being in the HZ, see Sec. 6.1.2) but for warmer planets in the HZ (e.g. planets f and g), a rain cycle could create liquid water on the planets that is then reinjected into the atmospheres cyclically. Volatiles such as CO, CO2, or CH4 have a low condensation temperature and will remain in the atmosphere along with other similar volatiles delivered by the comets. However, when liquid water is on the planet, this can draw the CO2 content down by silicate weathering that fixes CO2 in the planet's surface (forming carbonates) as shown in Siever (1968).
Over longer timescales, these volatiles can further chemically react to form new molecules. However, the exact composition of the delivered volatiles depends on the composition of exocomets in our scenario. The latter has been found to be consistent with the composition of comets in our Solar System (e.g. Matrà et al. 2017a) but there is still a wide range of observed compositions amongst the Solar System's comets (e.g. Mumma & Charnley 2011).
Another complication is that, as discussed in the previous subsection 6.2, volatiles may also be formed from the vapourised planet's crust during impact, from outgassing and even by reaccretion of previously ejected material, which would mix with the volatiles delivered by impacts.
All of these factors (active chemistry, potential additional volatiles, exocomet composition) makes it hard to predict the exact final compositions of the atmospheres after a few Gyr of evolution. An atmosphere model that would make assumptions about what happens without impacts could be fed by our impact predictions to come up with a plausible likely composition, but this goes beyond the scope of the present paper. We note however that these extra sources of volatiles do not change our conclusion that in the presence of a belt scattering comets, the atmospheres of the outermost planets f, g, h should be more massive.

Impacts in very dense Venus-like atmospheres
We note that our model is not valid for very massive atmospheres. If the atmospheres of planets f, g, h become massive enough (Venus-like, i.e., 200bars) due to impacts (or if the primordial atmospheres were Venus-like), 1-10km impactors do not create craters anymore but rather decelerate and get fragmented before touching the ground and create big aerial bursts that are very effective at removing atmospheric mass (Shuvalov et al. 2014). The amount of accreted projectile material is also very high (close to 100%) for these aerial burst type of impacts (Shuvalov et al. 2014). Therefore, for very dense atmospheres, we expect an increased delivery of volatiles from the impactors and less from the vapourised crust. We also expect that Venus-like primordial atmospheres would still be destroyed, since those impacts are more effective at removing mass, therefore not changing our conclusions.

Implications for life on these planets
One of the prime motives in searching for planets orbiting very low mass stars is to study the chemical composition of their atmospheres, and discover whether they contain large quantities of gas of a likely biological origin (e.g. Seager et al. 2016). Here, we consider the implications of our results concerning impacts towards creating the first forms of life.
Many elements can affect the emergence of life, most of which currently remain unconstrained empirically. We chose to apply our study to the TRAPPIST-1 system because its seven planets mark an important milestone. In addition to the multiple advantages of having a very low-mass host star for atmospheric characterisation (e.g. He et al. 2017), these seven worlds allow us to compare each to one another. All seven have followed a similar history in terms of UV irradiation for instance (modulo their distance to the star). Here we have tried to quantify whether all planets would receive a similar impact history, which may be important to kick start life as explained further.
UV irradiation has often been seen as prejudicial to habitability. Its main disadvantages are: 1) to photodissociate water molecules, of which the hydrogen is then lost to the space, depleting its oceans (e.g. Bourrier et al. 2017), and 2) to break complex molecules on the surface, and affect replication (e.g. O'Malley-James & Kaltenegger 2017). The situation is particularly sensitive for planets orbiting very low-mass stars like TRAPPIST-1, since these spend a long time contracting onto the main-sequence, in a 1 Gyr stage of particularly heightened far UV activity (e.g. Rugheimer et al. 2015).
However, these issues might be mitigated by several effects: 1) Ocean loss depends on the initial water reservoir (e.g. Ribas et al. 2016;Bolmont et al. 2017), and the TRAPPIST-1 planets might have been initially rich in water, having possibly assembled beyond the snow-line (Alibert & Benz 2017;Ormel et al. 2017) and/or accreted water at a later stage owing to impacts (as shown in this study); 2) UV photons do not penetrate water well, and organisms can protect themselves under a few metres of water (e.g. Estrela & Valio 2017); 3) UV irradiation accelerates mutations, leading to Darwinian evolution; 4) the non-illuminated side of a tidally synchronised planet is protected; and 5) UV irradiation, impacts, and a hard surface might be required to kick-start life (abiogenesis).
The literature contains much debate on many of the points above, except on the very last one, which we describe in more detail here as it is related to the outcome of this paper.
Recent advances in biochemistry (summarised in Sutherland 2017) have shown a prebiotic chemical path leading from hydrogen cyanide (HCN) to formaldehyde (CH2O), a known precursor to ribonucleotides (the building block to biologically relevant molecules such as ATP, RNA and DNA), amino acids (required for proteins) and lipids (Patel et al. 2015). Hydrogen cyanide, the initial molecule needed to inititate the process, can be produced in the plasma created when impactors enter in contact with an atmosphere (Ferus et al. 2015). In the presence of UV radiation, hydrogen cyanide can then react with other compounds that can be found concentrated on a planetary surface to create the building blocks of life. The impactor itself may have another role to play, which is to excavate underground material, and reveal chemically interesting strata (Patel et al. 2015), thereby acting as a chemical reactor.
We show in this paper that, if a belt scattering comets is present in the system, numerous impacts with different energies will happen throughout the history of the TRAPPIST-1 planets. From these impacts, we expect to create a subsequent amount of HCN in the impactor plasma (Ferus et al. 2015). We also note that as HCN is found in comets (e.g. Mumma & Charnley 2011), it may also be present on the potential exocomets of TRAPPIST-1 and be delivered along with the other volatiles (e.g. see Matrà et al. 2017b). We also emphasise that if the planets are tidally locked, it does not affect the emergence of life in this scenario as we predict that about half of the impacts would happen on the night side and the other half on the day side so that the UV photons from the star necessary for reactions to happen will be able to play their role.
Thus, our scenario offers the seed to create the first building blocks of life and more detailed modelling is needed to quantify how many ribonucleotides, amino acids and lipids could be created from the impact properties (e.g. impact velocities, rate of impacts) we predict. This is beyond the scope of this paper but should give birth to new interesting studies in the near future. Panspermia may also be viable to transport some potential life forms to other planets, which can enhance the probability of life spreading in the system (Krijt et al. 2017;Lingam & Loeb 2017).
To conclude, we cannot be certain yet that such a path is where biology originated, however, it provides a different narrative, one that requires UV irradiation, impacts and a limited amount of water. Ultraviolet, in this context, becomes beneficial by removing excess liquid water and transforming hydrogen cyanide into formaldehyde, whereas impacts would bring in energy to create hydrogen cyanide, and replenish the planet in volatiles such as water, much like what happened in the LHB (e.g. Court & Sephton 2014;Nesvorný et al. 2017) after a desiccating moon-forming impact (e.g. Canup 2014).

CONCLUSION
In this paper, we have studied the effects of impacts on the seven TRAPPIST-1 planets in terms of atmospheric mass loss and delivery of volatiles and water. We derive general results for any scenario where the comet pericentres slowly migrate inwards at a rateq. We also specifically test three scenarios for the delivery of comets from an outer belt to the inner planets (located within 0.1au): 1) Planet scattering by a single or a chain of planets, 2) An outer companion forcing Kozai oscillations on comets leading them to small pericentres, 3) Galactic tides on an exo-Oort cloud. We model these three scenarios by a steadily decreasing pericentre (constanṫ q) that is quantified in Sec. 4 for each of the scenarios. The results can be summed up as follows: • We find that applying a minimum mass TRAPPIST-1 nebula approach lead to a surface density Σ ∼ 122 (r/1au) −1.97 kg/m 2 . We show that a potential belt around TRAPPIST-1 could not survive within 10au because of collisional erosion (if it was created at the end of the protoplanetary disc phase). Assuming that such a belt is between 10 and 50au, and extrapolating the derived minimum surface density, we infer that this belt would have mass of at least 20M⊕ and may be observable in the far-IR or sub-mm with ALMA (see Sec. 2).
• We ran a suite of N-body simulations to understand the dynamics of comets that impact onto the seven different planets. We find the impact and ejection probabilities for each comet's orbit (see Figs. 2 and 3). We also provide the accretion timescales for these different comet families (see Fig. 4). We analytically explain the main dependencies for these probabilities and timescales.
• We give the impact velocity distributions for each planet, and we find that they typically have double-peaked profiles (see Fig. 6). The median impact velocity for planet b is close to 100km/s, whilst for planet h, it is close to 20km/s (see Fig. 7). These impact velocities are always much above the escape velocities of the planets and gravitational focusing is not important.
• We find that the fraction of comets accreted on each planet depends on the decreasing rate of pericentres (q) and apocentre Q (scaling asq −1 Q −1 ). We find two regimes, for smallq, most of the impacts end up on planets g and h and for higherq, each planet gets a fraction of comets accreted (see Fig. 12).
• The atmospheric removal is dominated by comets of a few km in diameter (see Fig. 13 left).
• The delivery of volatiles is only possible for comets 3km in size (see Fig 13 right). For bigger comets, the projectile material escapes and no delivery is possible.
• We find that the higher impact velocities for the innermost planets lead to a higher atmospheric removal rate for a given cometary impact rate and a lower amount of volatile delivered.
• In general, we find that if the incoming mass of comets that reach the inner regions Minc > 5 × 10 −4 q 10 −5 au/yr Q 1au M⊕, the primordial atmospheres of the seven planets would be totally destroyed (see Fig. 15), i.e a belt with a low scattering rate similar to the current Kuiper belt is enough to destroy all primordial planetary atmospheres.
• We quantify the amount of water lost owing to impacts and find that it is similar (possibly higher) to the amount of water lost through hydrodynamic escape (see Sec. 5.3.3 and Table. 3).
• As for the delivery of volatiles to the comets (see Fig. 16), we find that planets g and h (and most likely f) may retain volatiles from the impacting comets in their atmospheres and the conclusion holds for any size distribution of incoming comets between -3 and -4 (see Fig. 17).
• We thus predict that if the planets were hit by comets, the atmospheres of planets f, g, and h would be more massive, which could be checked by future missions in the next decade.
• We also show that for an incoming mass of comets Minc > 5 × 10 −4 f −1 vol q 10 −5 au/yr Q 1au M⊕ (where f vol is the volatile fraction on solids), the volatile mass delivered by comets is greater than Earth-like atmospheric masses (assuming Earth-like densities for the 7 planets).
• We provide a prescription for the amount of water or volatiles that can accumulate as a function of time (see Eq. 21) that could be used to feed an atmospheric model to check the actual composition of atmospheres dominated by the delivery of comets.
• We find that a large quantity of volatiles may have been delivered to planet h while it was still in the liquid water habitable zone.
• We find that a planet chain that would scatter comets from an exo-Kuiper belt or an outer companion that would force Kozai oscillations on a comet belt are two plausible mechanisms to throw an important number of comets on the seven planets over the system's lifetime (see Secs. 4.1 and 4.2).
• On the other hand, we rule out a potential Oort-cloud around TRAPPIST-1 as being a significant source of impacting comets (see Sec. 5.1.2).
• For the planet-scattering scenario, we find that even a belt with a low scattering rate similar to the current Kuiper-belt is enough to destroy typical Earth-like primordial atmospheres for the seven planets. Taking into account that typically observed debris belts are much more massive than the Kuiper belt, we find that the Kozai (slightly less efficient) scenario can also strip primordial atmospheres even if the impact process only lasts a fraction of the system's age.
• As for the volatile delivery, we find that for the planetscattering scenario, planets f, g, and h can get (more than) an Earth ocean mass of water (and other volatiles) delivered, which can accumulate impact after impact. We find that the primordial atmospheres are gradually replaced by cometary material and may lead to subsequent build up of new secondary atmospheres with exocomet-like compositions. These new secondary atmospheres may become more massive than the initial primordial atmospheres.
• Table 2 summarises the results for the different scenarios as for the minimum scattered (incoming) mass needed to destroy the primordial atmospheres and the volatile/water masses that can be delivered onto each planet.
• We also discuss the implications of impacts to create the building blocks of life. We detail new emerging pathways that can lead to life showing that UV irradiation, impacts and a hard planetary surface might be enough to kick start biological reactions and form ATP, RNA, DNA, amino acids and lipids that are essential to life (see Sec. 6.5).
In brief, we find that the primordial atmospheres of the seven planets orbiting around TRAPPIST-1 would not survive over the lifetime of the system if a belt scattering comets at a similar low rate than the Kuiper belt (or faster) were around TRAPPIST-1. According to our calculations based on applying a minimum mass extrasolar nebula approach for the TRAPPIST-1 system, we expect a potential 20M⊕ belt may have survived around TRAPPIST-1 that would be observable with ALMA. We also show that a large fraction of the delivered cometary volatiles remains in the atmospheres of the outermost planets f, g and h, which gradually replace their primordial atmospheres. We predict that the new secondary atmospheres of planets f, g and h may be more massive than that of the innermost planets (which may soon be checkable/observable with the JWST) and their composition might be dominated by the composition of exocomets in this system (i.e., impacts leave an imprint). We also predict that more than an Earth ocean mass of water could be delivered to planets f, g, and h owing to impacts that may be in liquid form on planets f and g.

ACKNOWLEDGMENTS
This paper is dedicated to Mila. We thank the two referees for comments that greatly improved the quality of the paper. QK and MCW acknowledge funding from STFC via the Institute of Astronomy, Cambridge Consolidated Grant. QK thanks J. Teyssandier for interesting discussions about the Kozai mechanism. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound.