© ESO 2020

1 Introduction

The first exoplanet discovery (Mayor & Queloz 1995) around a main-sequence star started a new chapter in the field of planetary sciences. More than 4000 confirmed exoplanets have been detected as of May 2020 and the number doubles every ~27 months¹. The variety of the discovered exoplanets shows that planetary systems are very common (planets are probably as common as stars) and diverse (Perryman 2018). There are exoplanets with masses ranging from the Moon’s mass to hundreds of Jupiter masses, with orbital periods ranging from hours to thousands of years and with a great variety in mass densities from that of H₂ gas to that of solid iron (Winn & Fabrycky 2015). A large fraction of the discovered exoplanets do not have analogues in the Solar System; this is the case for super-Earths, with masses between about 1 and 10 Earth masses, and hot Jupiters, giant planets orbiting very close to their host stars (Raymond et al. 2019), for example. Importantly, though, there appears to be a population of terrestrial planets at distances to the star that may allow liquid water to exist on the surface, possibly allowing life as we know it to form.

Apart from looking for analogues of the Solar System planets, a search for analogues of small body populations is ongoing. Recent observations(Acke et al. 2012; Booth et al. 2016; Close 2010; Lagrange et al. 2010; Matthews et al. 2014; Moór et al. 2013; Su et al. 2013; Welsh & Montgomery 2013) have found analogues of the main asteroid belt (MAB) and the Kuiper belt (KB) around ~20% of the nearest stars (Ren et al. 2019). Within our own Solar System, the small body belts (MAB and KB) play an important role in deliveringvolatiles to planets (e.g. Frantseva et al. 2018, 2020).

Delivery is one of the two major processes that has shaped the atmospheres of the terrestrial planets in our Solar System, with outgassing being the other (de Pater & Lissauer 2015; Smallwood et al. 2018; Chen et al. 2019). In the MAB, certain types of asteroids, C-types, together with comets and small bodies from the KB are known to contain significant amounts of water and organic compounds. The migration of these small bodies occasionally results in impacts with planets, which in turn leads to volatile enrichment of the planets. The existence of asteroid belts in other planetary systems implies that the same volatile delivery mechanisms might be at play around other stars. This way, volatile material, in particular water and organic compounds, can be delivered to exoplanets. Studying these delivery mechanisms is of astrobiological interest and will lead to better understanding of planetary formation, evolution, and habitability (Raymond et al. 2019). Enrichment of exoplanet atmospheres in refractory materials is of less astrobiological relevance, but may become observable using the Mid-InfraRed Instrument aboard the James Webb Space Telescope (JWST-MIRI) due to be launched in 2021. Silicate features in particular feature prominently in MIRI’s wavelength range (e.g. Rieke et al. 2015). The atmospheres of exoplanets beyond the snow line would not be expected to contain significant amounts of refractory materials, so any detection thereof would be diagnostic of impact-driven enrichment.

In this paper we focus on the exoplanetary system HR 8799, which is known to host both a warm and a cold debris belt, respective analogues of the Main Asteroid Belt and the Kuiper Belt. Modelling such a planetary system helps in understanding the interaction between planets and planetary debris.

HR 8799 is a nearby young A5V star with mass ≈1.56 M_⊙ (Gray & Kaye 1999; Gray et al. 2003; Baines et al. 2012; Goździewski & Migaszewski 2014). The exact age of the star is somewhat uncertain between 30 and 60 Myr (Marois et al. 2010; Zuckerman et al. 2011). HR 8799 is the host star of four giant planets HR 8799 e, HR 8799 d, HR 8799 c, and HR 8799 b detected via direct imaging (Marois et al. 2008, 2010). HR 8799 e, HR 8799 d and HR 8799 c have masses of approximately 9 M_jup, and HR 8799 b of 7 M_jup. For comparison, the giant planets in the Solar System have masses of 1 M_jup, 0.3 M_jup, 0.05 M_jup and 0.04 M_jup, see Fig. 1. HR 8799 e, HR 8799 d, HR 8799 c and HR 8799 b orbit their host star at distances of 15–69 AU, see Table 1 for more details. The Solar System giant planets orbit closer to the host star at distances of 5–30 AU.

The orbits of the HR 8799 planets are not well characterised observationally, yet. The available astrometry spans only ≃12% of the innermost planet’s orbit and ≃3% of the outermost planet’s. A variety of orbital configurations is consistent with the data, most of which are dynamically unstable; however, a few stable configurations were found (Goździewski & Migaszewski 2014, 2018; Götberg et al. 2016; Wang et al. 2018). The best long-term stable solution assumes that the four planets are locked in a protective mean motion resonance 1:2:4:8 (Goździewski & Migaszewski 2014, 2018). The higher planetary masses relative to the star as well as the larger eccentricities make planet-planet interactions more important in HR 8799 relative to the Solar System. This leads to many unstable orbital configurations consistent with the data (see above), and also requires extra care in numerical modelling. Detection of double Laplace resonances or higher-order resonant chains, such as proposed for the HR 8799 planets or is known for the Galilean satellites (Sinclair 1975; Peale & Lee 2002), is demanding (Michtchenko & Ferraz-Mello 2001; Gallardo et al. 2016). A good understanding of the orbital architecture of the system is essential for studying the system’s debris belts that have been observed with the Herschel Space Observatory, the Atacama Large Millimetre/submillimetre Array (ALMA) and the Spitzer Space Telescope (Su et al. 2009; Su & Rieke 2014; Matthews et al. 2014; Booth et al. 2016).

The system contains at least two distinct debris belts: an inner belt and an outer belt. The outer belt was discovered in imaging with the Multi-Band Imaging Photometer for Spitzer (Spitzer-MIPS) and is spatially resolved, while the inner belt is too small to be spatially resolved, but was inferred from the observed infrared (IR) excess in the spectral energy distribution (SED) using the Spitzer Infrared Spectrograph (Sptizer-IRS) data (Su et al. 2009). The inner belt lies between ~6 AU and ~15 AU, and the outer belt extends from ~90 AU to ~310 AU. The inner belt lower dust mass limit was estimated to be 1.1 × 10⁻⁶ M_⊕ and the outer belt’s is 0.12 M_⊕ (Su et al. 2009). These estimated belt masses are dust-only masses, while large objects (if they exist) would by far dominate the mass budget. The structure of HR 8799 resembles the structure of our own Solar System but on a larger scale as shown in Fig. 1. The Solar System up to the Kuiper belt can be fitted inside the orbit of the outermost planet of the HR 8799 system. Both systems contain an inner warm belt close to the host star, within ~6 and ~15 AU for HR 8799 and within ~2 and ~4 AU for the Solar System. Furthermore, both systems contain an outer cold belt, within ~90 to ~310 AU for HR 8799 and ~35 to ~100 AU for the Solar System. Between the two belts each of the systems have four giant planets. Remarkably, HR 8799 and the Solar System show the same overall structure although the HR 8799 system is much larger than the Solar System and its planets are much more massive (see Fig. 1). We note that terrestrial planets may exist inside the warm belt of HR 8799; those would not be detectable using current telescopes.

The interaction between the giant planets of HR 8799 and the debris belts have been studied by Contro et al. (2015, 2016); Read et al. (2018); Goździewski & Migaszewski (2018). Numerical simulations by Contro et al. (2016) indicate that the inner belt is structured, with gaps at the locations of mean motion resonances with the innermost planet similar to the Kirkwood gaps in the MAB. Moreover, the simulations suggest that the belt is located between 6 and 8 AU and collisions within the belt occur at velocities of the order of 1.2 km s⁻¹ or less, lower than the collisional velocities in the MAB, which are of order 5 km s⁻¹. We note that Contro et al. (2016) suggest that the outer edge of the inner belt has to be located at 8 AU, and not 15 AU as suggested by past observations by Spitzer (Su et al. 2009). The simulations demonstrate that due to the influence of the innermost planet it is not possible to keep a population of small bodies between 8 and 15 AU. Similarly, Read et al. (2018) performed N-body simulations to study theeffect of the giant planets on the outer belt and demonstrated that the belt has a similar structure as the inner belt, with gaps at the locations of the mean motion resonances with the outermost planet. Read et al. (2018) show that adding a hypothetical planet with mass 0.1−1 M_jup outside HR 8799 b, the outermost confirmed planet, pushes the belt outwards, providing a better fit to ALMA observations (Booth et al. 2016; Read et al. 2018; Goździewski & Migaszewski 2018). Such a planet cannot be detected with current telescopes.

In our own Solar System the interaction between the belts and the planets leads to scattering of small bodies within the system and beyond. Sometimes this scattering leads to impacts with the planets. For the terrestrial planets impacts can appreciably enrich the surface in various materials including volatiles. For the giant planets refractory enrichment is even observable, as in the case of the Shoemaker-Levy 9 impacts in 1994 (Atreya et al. 1999; Harrington et al. 2004; Fletcher et al. 2010). Here, we test for the first time how strong this effect is in the HR 8799 system. We perform several sets of dynamical N-body simulations with 1 600 000 test particles representing exoasteroids and exocomets. These simulations result in minor body impact rates with the giant planets that are then converted to volatile delivery rates. We also perform a set of N-body simulations of dust to analyse a potential mass transport from the inner and outer dust belts onto the HR 8799 planets.

In Sect. 2 we describe our N-body simulations. Section 3 presents the results of the numerical simulations, checks for consistency with previous results, and determines impact rates. In Sect. 4 we convert the latter to volatile delivery rates to the planets. Discussion and conclusions are presented in Sects. 5 and 6.

Fig. 1 Orbital structure of the Solar System (left) and the HR 8799 system (right). The inner belts are represented in red and the outer belts in blue. Planet masses are indicated by the size of the dots; each plot has a separate log scale. We note that the orbital eccentricities of the giant planets in the HR 8799 system are several times larger than the eccentricities of the Solar System giants, which lead to more impacts between planets and small bodies.

2 Numerical simulations

We have performed numerical simulations of the dynamical evolution of dust and minor bodies in the inner and outer belts in the exoplanetary system HR 8799 to study the volatile and refractory influx onto the four known giant exoplanets.

2.1 Minor bodies

To model the minor body populations of the inner and outer belts, we used the N-body integrator MERCURIUS (Rein et al. 2019) from the software package REBOUND (Rein & Liu 2012). We set up a model of the motion of one gravitationally dominant object (HR 8799; M_⋆ = 1.56 M_⊙) and N massive objects (in our case:N = 4 representingthe four planets in the system) under the influence of their mutual gravity over Myr timescales forward in time. Simulations focusing on the inner and outer belts were done separately, as described in Sects. 2.1.1 and 2.1.2. The inner belt is represented by 650 000 test particles and the outer belt by 1 450 000. Our first run of the outer belt simulation contained 650 000 test particles and one of the planets received 0 impacts. To overcome small number statistics we added more test particles to the outer belt simulation. To achieve a 3σ result while adopting $\sqrt{N}$, where N is the number of impacts, for the noisewe need N = 10 as a minimum. Therefore we require the total number of impactors on each planet to be always well above 10. Increasing the number of test particles resulted in total impacts on each planet to be more than 10 and helped us to achieve a 3σ result.

In our simulation, the test particles move passively under the influence of the combined gravitational potential of the star and the planets. The gravitational effect of the test particles on one another and on the massive objects is neglected. Test particles are removed from the simulation once they collide with a planet or with the star. Test particles are considered ejected from the planetary system and discarded when they exceed a belt-dependent heliocentric distance: 1000 AU for the inner-belt simulations and 10 000 AU for the outer-belt simulations following Contro et al. (2016) and Read et al. (2018), respectively. Timestep values were set separately for the inner and outer belts, see Sects. 2.1.1 and 2.1.2 for a detailed description. The simulation results are recorded by taking snapshots of the simulations at predefined times (6, 30, 60 and 70 Myr).

Table 1 presents the planetary initial conditions that have been adopted for our simulations. The masses and orbital parameters were taken from Goździewski & Migaszewski (2014, their Table 1). The inclination, i, of the four planets is the best-fitting inclination of coplanar orbits to the sky plane, while the inclination of the HR 8799 equator to the sky plane is ~ 23 degrees as follows from the statistical analysis of the rotational speed of A5 stars. For the longitude of the ascending node, Ω, Goździewski & Migaszewski (2014) assumed the same value for all four planets. The planetary radii (0.0005327 AU, 0.0005327 AU, 0.0005327 AU and 0.000538629 AU) were calculated using the “Jovian Worlds” mass-to-radius relation from Chen & Kipping (2017): $$\frac{R}{R_{\oplus }}=17.74{\left(\frac{M}{M_{\oplus }}\right)}^{-0.044},$$(1)

which is valid for 0.414 M_jup < M < 83.79 M_jup. The relation represents a population of planets with masses larger than a sufficient mass for gravitational self-compression that starts reversing the growth of the planet (therefore the exponent is negative). Radii are needed to model impacts: a test particle is considered to have impacted the star or a planet once it ventures within its radius. For the radius of the star, we adopt a value of 1.440 ± 0.006 R_⊙ (Baines et al. 2012) and 60 Myr for its age (Marois et al. 2010).

2.1.1 Inner belt

To model the inner belt of the HR 8799 we followed the initial conditions described in Contro et al. (2016). We used 650 000 test particles (500 000 in Contro et al. 2016) with semi-major axes between 1 and 10 AU.

We adopt 1 AU as the inner edge value of the simulations following Contro et al. (2016). The outer boundary, 10 AU, follows the estimated observational value as in Marois et al. (2010). The eccentricities of the particles are distributed uniformly between 0 and 0.1 and the inclinations are set between 0° and 5°. The remaining, angular, orbital elements are set to random numbers between 0° and 360°. The simulations were performed for 70 Myr forward in time. The timestep of the simulations is set to 7 days in order to resolve the orbits of the innermost test particles orbiting the star at 1 AU with an orbital period of 292 days.

2.1.2 Outer belt

The outer belt simulations were based on the initial conditions from Read et al. (2018). We created 1 450 000 test particles (50 000 in Read et al. 2018) with semi-major axes between 69.1 and 429 AU. The inner boundary of the outer belt is set to the semi-major axis of the outermost planet HR 8799 b. We note that ALMA observations (Booth et al. 2016) suggest that the inner edge of the disc has to be further out at 145 AU, which might be explained by an additional fifth planet. However, we stick to the initial conditions described in Read et al. (2018) for consistency. The outer edge of the disc follows from the ALMA observations by Booth et al. (2016). Eccentricities were set between 0 and 0.05 and inclinations between 0° and 2.86° (corresponding to 0.05 rad). The remaining, angular, orbital elements are randomly distributed between 0° and 360° as for the inner belt. The simulations were performed forward in time for 70 Myr. The simulation’s timestep was set to 0.48 yr in order to resolve the orbit of the innermost planet HR 8799 e, which has an orbital period of 48 yr, and any close encounters that would occur between the planet and test particles.

2.2 Dust

We model dust originating from the inner and outer belts. In this manuscript, we call all particles smaller than minor bodies collectively dust, even though, by definition, dust larger than 30 micrometres in diameter should be called meteoroids².

For the source body populations we use the orbital element distributions of the small bodies in the inner and outer belts at time 70 Myr from our N-body simulation, see Sects. 2.1.1 and 2.1.2. This is because the dust populations in HR 8799 based on our simulations are not dynamically stable for more than 40 Myr and the parent body configuration is quite similar between 30 and 70 Myr. For the inner belt we released 875 000 particles with diameters between 1.5 and 4.5 μm, uniformly distributed over 7 bins. For the outer belt we released 17 500 particles with diameters between 10 and 1000 μm, also uniformly distributed over 7 bins. These size distributions are chosen based on the Spitzer, Herschel and ALMA observations of dust grains in the inner and outer belts (Su et al. 2009; Matthews et al. 2014; Booth et al. 2016). The difference between the number of particles in both belts is caused by the difference in the computational power needed to complete the dynamical evolution of the dust particles: 50 Myr for the outer belt versus several Myr for the inner belt. We assume a dust bulk density of ρ = 2 g cm⁻³.

All particles were numerically integrated under the gravitational influence of four planets, the host star, and the effects of radiation pressure and Poynting-Robertson (PR) drag assuming a stellar luminosity of L_⋆ = 4.92 L_⊙ (Gray & Kaye 1999) using the SWIFT_RMVS_3 numerical integrator (Levison & Duncan 1994). Particles were tracked until they vanished in one of the sinks: too close to the host star (R <0.05 AU), too far from the host star (R > 1000 AU), or impact on one of the four planets. No particle in our simulation was stable for more than 70 Myr, that is, the age of the system. Particle orbital elements were recorded every 1000 yr, where each recorded set of orbital elements is assumed to be a unique representation of meteoroids. This in other words means that we assume that the dust population is in a steady-state and dust is generated every 1000 yr to feed the dynamically evolving dust cloud. This results in a simulation containing 875 000 and 15 000 dynamical pathways for the inner and outer belt dust populations respectively, that is sampled by millions of unique dust orbital configurations.

3 Simulation results: orbital structure, impact rates

3.1 Minor bodies

3.1.1 Inner belt

The snapshot at 6 Myr shows that around that time the inner belt develops structure due to the interaction between the test particles and the planets, see Fig. 2. The interaction with the innermost planet HR 8799 e shapes the belt in a similar way Jupiter shapes our own Main Asteroid Belt. There are broad gaps at the locations of the mean motion resonances with the giant planets and there are many dynamically excited objects with large eccentricities and inclinations. At the end of the simulation, after 70 Myr, 0.143% of test particles were discarded due to collisions with the star and 31.224% of test particles were found to be ejected from the system. As shown in Table 2, 0.248% of test particles collided with the innermost planet HR 8799 e, 0.014% with HR 8799 d, 0.007% with HR 8799 c, and 0.002% with HR 8799 b. As seen from Fig. 3, the number of impacts peaks prominently in the first 1 Myr. We attribute this to test particles that were initialised on highly unstable orbits. After the first ~1 Myr, impacts continue at a roughly continuous rate, see Fig. 4; this is referred to as steady state in the following. The steady stateimpact rates is what we are interested in for the purposes of this project.

The initial conditions of the inner belt simulations followed the example of Contro et al. (2016). As can be seen in our Fig. 2 (see Fig. 3 in Contro et al. 2016), we reproduce the resulting structure of the inner belt from their simulations. At the end of the simulations there are almost no particles left beyond 8 AU. At the location of the 3:1 mean motion resonance with the innermost planet, at ~ 7.4 AU, a broad gapis formed as well as a smaller gap at the location of the 4:1 mean motion resonance (~ 6 AU). Unlike Contro et al. (2016) we resolve and analyse impacts between planets and test particles, see Sect. 4.

Fig. 2 Time evolution of 650 000 test particles representing minor bodies in the inner belt for a period of 70 Myr. All four planets are located outside of the plot. The innermost planet HR 8799 e has a semi-major axis of 15.4 AU.

Fig. 3 Fraction of test particles representing minor bodies originating from the inner belt that impacted with the planets, the star or exceeded a user-provided heliocentric distance, 1000 AU, during 70 Myr of the simulations.

Fig. 4 Same as Fig. 3, but only after the first 1 Myr of the simulations, when the system reached steady state.

3.1.2 Outer belt

In analogy with the inner belt, the structure of the outer belt is clearly visible by the time of the first snapshot at 6 Myr due to the interaction between the test particles and the planets as shown in Fig. 5. After 70 Myr, 0.0143% of test particles collided with the star and 9.4581% of test particles were ejected from the system. As seen in Table 2, the outermost planet HR 8799 b suffered impacts from 0.0410% of test particles, HR 8799 c from 0.0077%, HR 8799 d from 0.0050% and the innermost HR 8799 e from 0.0037%. Figure 6 shows that most impacts occurred in the first 1 Myr. The majority of the impactors is produced by the test particles initialised on highly unstable orbits. Steady state is reached around 1 Myr. We are interested in those impacts that occur after steady state was reached, see Fig. 7. The structure of the outer belt in our simulations is the same as in previous work (compare to Fig. 2 in Read et al. 2018). The range of orbits within the initial inner edge of the outer belt, 69.1 AU, and 100 AU is depleted by the end of the simulations. At the location of the 2:1 mean motion resonance with the outermost planet, a broad gap is formed. Also the simulations produce a few HR 8799 b “Trojans” orbiting at the same semi-major axis as the planet. Unlike Read et al. (2018) we resolve and analyse impacts between planets and test particles, see Sect. 4.

Fig. 5 Same as Fig. 2, but for the outer belt. 70 Myr evolution of 1 450 000 test particles representing minor bodies. The outermost planet HR 8799 b has a semi-major axis of 69.1 AU.

3.1.3 Impact rates

The number of test particles representing minor bodies impacting the planets from both belts are shown in Table 2. We only consider those impacts that occurred after the system reached steady state, 1–70 Myr. The innermost planet receives 7.2 and 0.06 impacts Myr⁻¹ from the inner and the outer belts. Planet d receives 0.6 and 0.1 impacts Myr⁻¹, while planet c 0.3 and 0.09 impacts Myr⁻¹. The outermost planet receives 0.07and 0.3 impacts Myr⁻¹.

The total number of test particles, N_tp, in the inner belt is 544 659 and 1 347 075 in the outer belt. We note that these values represent the total number of test particles in either belt after the first 1 Myr when the system reached steady state and differ from the 650 000 test particles and 1 450 000 per inner and outer belts initialised at the start of the simulations.

3.2 Dust

In our numerical simulations of two distinctive dust populations originating in the inner and outer belts, we recorded the dust dynamical evolution in the HR 8799 system. Since our simulations do not include effects of dust-dust collisions that are estimated in the post-processing, we obtained the complete dynamical history of each simulated particle. These simulations resulted in 150 million unique records of orbital elements for the outer belt dust population, and 188 million records for the inner belt. Due to their significantly different source regions and size range, we describe each population separately in the following paragraphs.

The dynamical pathways of the outer belt dust (10 μm < D < 1000 μm) are mostly driven by a slow migration into outer mean-motion resonances with HR 8799 b, influence of secular resonances and planetary close-encounters that lead to scattering of dust into outer space (R > 10 000 AU). Only less than 5% of dust is able to migrate below R < 3 AU assuming that there are no dust-dust collisions. The average collision probability P of particles migrating via PR drag from the outer belt is P ≈ 10⁻¹⁰ per year for HR 8799 b, and decreases by a factor of 2 for each subsequent planet. We calculate P adopting the collision probability between orbiting objects from Kessler (1981). In reality, the mass of the outer dust ring is 18 million times higher than the mass of the Zodiacal Cloud (Nesvorný et al. 2011) 4 × 10¹⁹ g vs. 0.12 M_⊕ = 7.2 × 10²⁶ g, and therefore destructive dust-dust collisions happen on orders of magnitude shorter timescales than the PR drag dynamical timescales, as suggested in Wyatt (2005) for young debris discs like HR 8799.

To check the effects of dust-dust collisions in the HR 8799 system, let us assume that the HR 8799 outer belt dust has a similar shape as the Zodiacal Cloud, and use the collisional balance estimates from Grun et al. (1985). We multiply the density of the cloud by a factor of 18 million and scale it to a distance of 100 AU. This results in PR drag driven migration, where only the smallest D = 10 μm particles can reach HR 8799 b (a = 69.1 AU) with an impact probability P = 10⁻¹³ per year. For any dust larger than D = 10 μm the impact probability even with the outermost planets drops quickly to zero for D < 21.5 μm; P < 10⁻³⁰ per year as shown in the Fig. 8 panel B. From this we conclude, that shortly after the outer belt is shaped by the gravitational interaction with planets (Fig. 3, 6 Myr), the dust contribution from the outer belt is negligible due to the overwhelming dust-dust collisional rate. This is also independently confirmed in Wyatt (2005), who concludes that for all resolved debris discs with ages < 100 Myr the effects of PR drag are insignificant. The source populations of dust are far behind the outermost planets and thus particles cannot reach planet crossing orbits before their collisional demise. The dust fragments are either blown out of the system or subject to more collisions. The PR timescale for a D = 10 μm particle starting at a = 100 AU on a circular orbit to reach planet HR 8799 b is 7.2 Myr (Wyatt 2005), for our assumed bulk density of ρ = 2 g cm⁻³. We can thus effectively disregard any possible delivery of material from the outer dust belt on any of the four planets.

The situation is vastly different for the inner dust population due to the particle size range, 1.5 μm < D < 4.5 μm. These sizes are very close to the blow outregime, β > 0.5 for circular orbits of their parent bodies, where β ≈ (1.15 × 10⁻³∕ρD)(L_⋆∕L_⊙)(M_⊙∕M_⋆) (Wyatt et al. 1999), where ρ and D are in SI units. The critical diameter when β = 0.5 is D_⋆ = 3.85 μm for ρ = 2 kg m⁻³. Since our parent bodies have non-zero eccentricity, particles smaller than D_⋆ = 3.85 μm that are ejected close to aphelion can retain bound orbits. Particles with D < 3.0 μm have almost zero impact probability with all four planets (P < 10⁻¹³ per year as shown in the Fig. 8 panel A), thus we conclude they do not deliver any significant portion of material to any planet. Dust particles with 3.5 < D < 4.5 μm have the highest impact probability (P ~ 10⁻⁸ per year) with the innermost planet, and the impact probability decreases by a factor of 2 for each subsequent planet further from HR 8799. The total mass of dust in the inner belt is M_dustInner = 1.1 × 10⁻⁶ M_⊕, so at the impact rate we calculated here, the HR 8799 planets would in 100 Myr accrete < 10⁻⁶ M_⊕, which is 9 orders of magnitude less than their estimated mass.

To summarise, our numerical integrations show that both the inner and outer dust populations do not significantly contribute to material delivery onto any of the four known planets orbiting HR 8799. This is unlikely to change until the rate of collisions/optical depth of the outer belt drops significantly and allows inward migration of dust via PR drag.

Additionally, it is important to highlight that the dust size ranges, 1.5−4.5 μm for the inner belt and 10−1000 μm for the outer belt, have been chosen based on the existing observations of the system. Such discontinuous size ranges between these two populations can be due to the observational bias because Spitzer is sensitive to small grains and ALMA to large ones. Since we do not have observational constraints on the distribution of the larger grains in the inner belt and the smaller grains in the outer belt, we do not try to add these size ranges in our simulations.

Fig. 6 Fraction of planet impactors representing minor bodies originating from the outer belt after 10, 20, 30,..., 70 Myr, respectively.

Fig. 7 Fraction of planet impactors representing minor bodies which happened after the first 1 Myr of the simulations, after the system reached steady state.

Fig. 8 Panel A: average collisional probability with each HR 8799 planet for 7 different particle diameters originating in the inner belt. Each colour corresponds to a different planet, where the dashed lines are assuming no particle-particle collisions, while solid lines are including collisions (see the main text for details). The units are particle impacts onto the entire planet surface per year. Panel B: same as panel A but now for the outer belt particles and 6 different particle diameters. In order to show the difference between scenarios with and without particle-particle collisions, we added a constant value of P = 10−30 to all data points in Panel B. Except for D = 30 μm impacts ontoplanet HR 8799 b all particles with D ≥ 30 μm have P ≪ 10−30 when particle-particle collisions are included.

4 Delivery rates

In Sect. 3 we have derived impact rates. Here, we use these to estimate the corresponding delivery rates of volatile (water and organics) and refractory (metals and silicates) material to the four planets. In Sect. 3.2 we have demonstrated that the dust populations of both belts do not significantly contribute to the delivery of material to the giant planets, therefore the delivery rates calculated in this section refer only to the delivery by minor bodies.

The volatile delivery rate R_vol and the refractory delivery rate R_refr, for each planet, can be expressed as (see Schwarz et al. 2018, Eq. (3)): $$R_{\text{vol}}=\frac{M_{\text{belt}}\times f_{\text{vol}}\times N_{\text{imp}}}{N_{\text{tp}}{\tau }_{\text{sim}}},$$(2) $$R_{\text{refr}}=\frac{M_{\text{belt}}\times f_{\text{refr}}\times N_{\text{imp}}}{N_{\text{tp}}{\tau }_{\text{sim}}},$$(3)

where M_belt is the total mass of the belt, f_vol and f_refr are the mass fraction of volatiles and refractories therein, N_imp is the number of impacts on a planet from the corresponding belt, N_tp is the number of test particles in the corresponding belt after the first 1 Myr of the simulation and τ_sim is the simulation time. The simulation time τ_sim is 69 Myr and not 70 Myr since we do not consider the impacts that happened in the first 1 Myr of the simulations. To estimate values of M_belt, f_vol and f_refr for the inner belt we use the Main Asteroid belt as a proxy and the Kuiper belt for the outer belt.

The values of N_imp, N_tp, τ_sim have been determined in Sect. 3.1.3. In this section, we estimate M_belt, f_vol and f_refr.

4.1 Belt masses

The total dust mass of either belt was determined by dust modelling based on infrared observations. According to Su et al. (2009) the total mass of the 1.5−4.5 μm sized dust grains in the inner belt is M_dustInner = 1.1 × 10⁻⁶M_⊕ and the mass of the 10−1000 μm sized dust grains in theouter belt is M_dustOuter = 0.12 M_⊕. Geiler et al. (2019) created a disc model that would explain Herschel and ALMA observations of the outer belt. Their “preferred” model predicts M_plOuter = 134 M_⊕ as the total mass of the outer belt, assuming a size-frequency distribution of minor bodies up to a maximum diameter of 100 km. We adopt this value for our further calculations.

For the total minor body mass of the inner belt, however, there are no published estimates that we are aware of. We estimate the total belt mass (minor bodies only) based on the dust mass determined by Su et al. (2009). In doing so, we assume a power law describing the differential size-frequency distribution n(a) (diameter a) of a colliding steady state population of bodies: $$n(a)\propto a^{-q},$$(4)

where the power-law index q = 3.5 (Dohnanyi 1969; Bottke et al. 2015). The total mass of the inner belt for minor bodies can be calculated by integrating over the size distribution from the minimum to maximum dust sizes (from a_minD = 1.5 to a_maxD = 4.5 μm as observed by Su et al. 2009) and minor body sizes (a_minD = 1 m to a_maxD = 1000 km, based on the size distribution of the Main Asteroid Belt): $$m(<a_{\text{max}})=\frac{\pi }{6}{\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}}n(a)\rho a^3\text{d}a\propto \sqrt{a}{|}_{a_{\text{min}}}^{a_{\text{max}}}.$$(5)

Normalising the total inner belt mass for minor bodies by the total inner belt mass for dust particles gives: $$\frac{M_{\text{plInner}}}{M_{\text{dustInner}}}=\frac{\sqrt{a_{\text{maxA}}}-\sqrt{a_{\text{minA}}}}{\sqrt{a_{\text{maxD}}}-\sqrt{a_{\text{minD}}}}.$$(6)

In doing so, we found the total mass of the inner belt to be M_plInner = 1.23 M_⊕. For comparison, the estimated dust mass of the inner Solar System is 7 × 10⁻⁶ M_⊕ (Nesvorný et al. 2011) and the estimated mass of the Main Asteroid Belt is 0.0004 M_⊕ (de Pater & Lissauer 2015). The planetesimal to dust mass ratio for the inner belt of the HR 8799 is 10⁶ compared to 10³ for the Solar System.

4.2 Volatile content

While it is reasonable to assume that the debris belts in HR 8799 contain some volatile material, no observational constraints are available. We estimate the volatile content f_vol of the inner and outer debris belt, basing ourselves on knowledge from the Solar System.

Solar System asteroids are classified in several taxonomic types. Only members of one such type, the so-called C-type asteroids, show appreciable amounts of volatiles. About 33% of all asteroids belong to this type (DeMeo & Carry 2013). The water content of the C-type asteroids is ≈ 10% by mass and their carbon content is ≈2% by mass (National Research Council 2007; Sephton et al. 2002; Sephton 2014). Therefore, we adopt 0.33 × (0.10 + 0.02) = 0.04 as the f_vol value for the inner belt of the HR 8799 system.

In contrast, the Kuiper Belt objects all are formed beyond the snow line, which means that on average all those objects are expected to be volatile rich. However, the ice fraction measured in KBOs spans the full range from 0 to 1 (Brown 2012). Since the range of the volatile fraction is rather wide we adopted an average value of 0.5 as the f_vol for the HR 8799 outer belt. The same value was adopted by Ciesla et al. (2015) and Schwarz et al. (2018) in their studies of volatile delivery.

After inserting all necessary values into Eq. (2) we calculate the volatile delivery rates, as shown in Table 3 from the minor bodies of the inner and the outer belts to the four giant planets. The outer belt delivers an order of magnitude more volatiles to the planets despite the smaller total number of impacts. Especially for the outer belt our simulations yield few impacts resulting in large Poisson noise, up to 50%, leading to large uncertainties in delivery rates.

4.3 Refractory content

Like for volatiles, we base our estimates of refractory-material content on our knowledge of the Solar System. In particular, we assume that any material not counted as volatile is refractory; in practice, the refractory component is likely to be dominated by silicates and metals. We therefore adopt a refractory content f_refr of 0.96 for the inner belt and 0.5 for the outer belt, respectively. The resulting delivery rates are presented in Table 3.

5 Discussion

As Table 3 shows, volatiles and refractories are delivered from both belts to all four planets by minor bodies. The inner belt delivers to the planets 7.3 × 10⁻⁷ M_⊕ of volatile material per Myr and 1.7 × 10⁻⁵ M_⊕ of refractory material per Myr. The outer belt delivers 2.7 × 10⁻⁵ M_⊕ of volatiles per Myr and 2.7 × 10⁻⁵ M_⊕ of refractories per Myr.

In terms of the impact rates the inner belt dominates the outer belt except for planet HR 8799 b. For example, planet HR 8799 d experiences 0.6 impacts Myr⁻¹ from the inner belt and 0.1 impacts Myr⁻¹ from the outer belt. However, the outer belt is much more massive than the inner belt, two orders of magnitude, and its volatile fraction is an order of magnitude higher than for the inner belt. As a result, the volatile delivery to the planets is dominated by the outer belt. On the other hand, the refractory fraction of the inner belt is almost twice as large as for the outer belt, which leads to similar refractory delivery rates from both belts.

Over the course of 69 Myr the four giant planets receive between 2 ± 1 × 10⁻⁴ M_⊕ and 9 ± 2 × 10⁻⁴ M_⊕ of volatile material and between3 ± 1 × 10⁻⁴ M_⊕ and 13 ± 2 × 10⁻⁴ M_⊕ of refractorymaterial. The uncertainty is dominated by the Poisson noise in the simulated number of impacts, and by the planetesimal masses of the belts. These total volatile and refractory fluxes are small compared to the total planet mass, which is 10⁶ −10⁷ times larger. Such an amount of volatile and/or refractory infall may be detectable in the upper layers of the planetary atmospheres. For example, gases with mixing ratios of ~10⁻⁶ and even lower are often detectable, depending on species (de Pater & Lissauer 2015).

Since the four giants HR 8799 e, d, c, b are beyond the snow line (and presumably formed there), we expect them to be born volatile-rich.Any future detection of volatiles would therefore not necessarily imply delivery through impacts. Silicates or other refractory material would be more diagnostic in this regard. Past observations of Jupiter after the impact of comet Shoemaker-Levy 9 and smaller objects may be used as a proxy for post impact enrichment observations of the giant planets in the HR 8799 system (Atreya et al. 1999; Fletcher et al. 2010).

The HR 8799 system may contain terrestrial planets, which presumably formed dry within the snow line. Volatile delivery from the belts may be of astrobiological importance for those. The same is true, probably, for possible terrestrial planets around other stars with minor body belt analogues. Given the lower mass of terrestrial planets (an Earth analogue would be about 2800 less massive than HR 8799 e), the relative contribution of impactor material to the overall composition of the planet and its atmosphere could be much higher. Volatiles derived by impacts could be enough to explain an Earth-like atmosphere mass: Earth’s atmosphere contributes about 10⁻⁶ of Earth’s total mass. We caution, however, that we did not (yet) model the impactor flux on terrestrial planets in the HR 8799 system.

We also analysed a potential transport of material from the dust belts to the four giant planets. For the outer belt model only meteoroids with D = 10 μm deliver material to HR 8799 planets with a non-zero probability due to their short PR timescales. We can extrapolate this behaviour to particles smaller than that, which are migrating even faster due to their larger β values, and thus are less susceptible to destruction in mutual meteoroid collisions. However, too small particles, D < 3 μm, are immediately blown out of the system, thus only a narrow size range D = 3−10 μm remains. This fact greatly diminishes any possibility of any considerable mass transport from the outer belt to planets via micron sized dust. Particles with larger diameters (D > 1000 μm) migrate inwards slower than D < 1000 μm, are destroyed in dust-dust collisions, and should not deliver material to HR 8799 planets.

Inner belt particles with D > 4.5 μm are released on stable orbits. These orbits are migrating inwards due to PR drag. Due to the effects of the radiation pressure the initial aphelion distance of released particles can extend beyond the orbits of the HR 8799 planets, however particles with D > 10.2 μm, released from a = 8 AU and e = 0.1 bodies have aphelion distances Q < 15.4 AU, meaning that orbits not crossing any of HR 8799 planets. It is possible to increase the semimajor axes of the inner belt dust via close encounters or secular resonances, however the dust-dust collisions are expected to act on much shorter timescales, and we do not believe that this outward particle migration is efficient enough to effectively deliver a significant portion of material onto the HR 8799 planets. Moreover, the current constraints give a very low mass of the inner belt dust (M_dustInner = 1.1 × 10⁻⁶M_⊕), which further decreases the possibility of material transport from the inner belt via dust particles.

Our results show that the delivery process in the HR 8799 is dominated by minor bodies, while in the current Solar System delivery processes are dominated by dust populations. This follows from the fact that HR 8799 is a young non depleted system thus the minor body influx on planets dominates, while Solar System is 4.5 billion yr old, and most of the minor bodies are in stable reservoirs that do not intersect any planet.

Collisions between minor bodies in the outer and inner belts of the HR 8799 system produce dust, analogous to the zodiacaldust in the Solar system, that may be observable with the upcoming HOSTS (Hunt for Observable Signatures of Terrestrial Systems) survey (Ertel et al. 2018a,b). Possible observation of the exozodiacal dust around HR 8799 may provide insights on the existence of terrestrial planets in the system and on the nature of the belts.

6 Conclusions

We performed numerical simulations of the inner and outer debris belts in the exoplanetary system HR 8799 to calculate for the first time delivery rates from the belts to the four planets. Delivery rates were separated into a volatile and a refractory component, respectively, where we base ourselves on Solar System knowledge.

After the first 1 Myr the system reaches steady state. Also, in the first several Myr the HR 8799 belts develop orbital structure with gaps caused by the mean motion resonances with the planets (confirming previous work by Contro et al. 2016 and Read et al. 2018). All four planets experience minor body impacts from the inner and outer belts. The innermost planet HR 8799 e is affected the most by the objects from the inner belt. In turn, the outermost planet HR 8799 b experiences the most impacts from the outer belt. Our numerical integrations of the dust populations in the inner and outer belts show that there is no significant contribution from these populations to the materialdelivery to the giant planets. The outer belt delivers an order of magnitude more of the volatile material to the planets. The refractory delivery rates from the inner and the outer belt are similar, within the error bars.

We expect the four giant planets to be born volatile-rich, so volatile delivery through impacts is probably insignificant in comparison. The enrichment in refractory material, however, may well be significant. JWST-MIRI observations targeting silicate features could be especially diagnostic and should be considered.

Acknowledgements

We thank an anonymous referee for providing us with helpful comments that allowed us to improve and clarify the manuscript. We are thankful to Hanno Rein for valuable help with REBOUND, Jonathan Horner for input on the inner belt, Matthew Read and Mark Wyatt for input on the outer belt. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. We would like to thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Peregrine high performance computing cluster. PP’s work was supported by NASA’s ISFM award.

References

All Tables

All Figures

Fig. 1 Orbital structure of the Solar System (left) and the HR 8799 system (right). The inner belts are represented in red and the outer belts in blue. Planet masses are indicated by the size of the dots; each plot has a separate log scale. We note that the orbital eccentricities of the giant planets in the HR 8799 system are several times larger than the eccentricities of the Solar System giants, which lead to more impacts between planets and small bodies.

In the text

	Fig. 2 Time evolution of 650 000 test particles representing minor bodies in the inner belt for a period of 70 Myr. All four planets are located outside of the plot. The innermost planet HR 8799 e has a semi-major axis of 15.4 AU.
In the text

	Fig. 3 Fraction of test particles representing minor bodies originating from the inner belt that impacted with the planets, the star or exceeded a user-provided heliocentric distance, 1000 AU, during 70 Myr of the simulations.
In the text

	Fig. 4 Same as Fig. 3, but only after the first 1 Myr of the simulations, when the system reached steady state.
In the text

	Fig. 5 Same as Fig. 2, but for the outer belt. 70 Myr evolution of 1 450 000 test particles representing minor bodies. The outermost planet HR 8799 b has a semi-major axis of 69.1 AU.
In the text

	Fig. 6 Fraction of planet impactors representing minor bodies originating from the outer belt after 10, 20, 30,..., 70 Myr, respectively.
In the text

	Fig. 7 Fraction of planet impactors representing minor bodies which happened after the first 1 Myr of the simulations, after the system reached steady state.
In the text

Fig. 8 Panel A: average collisional probability with each HR 8799 planet for 7 different particle diameters originating in the inner belt. Each colour corresponds to a different planet, where the dashed lines are assuming no particle-particle collisions, while solid lines are including collisions (see the main text for details). The units are particle impacts onto the entire planet surface per year. Panel B: same as panel A but now for the outer belt particles and 6 different particle diameters. In order to show the difference between scenarios with and without particle-particle collisions, we added a constant value of P = 10−30 to all data points in Panel B. Except for D = 30 μm impacts ontoplanet HR 8799 b all particles with D ≥ 30 μm have P ≪ 10−30 when particle-particle collisions are included.

In the text