© The Authors 2025

1 Introduction

Multi-planet systems resemble our own Solar System the most, thereby making them the natural starting point in the search of planets where habitable environments could form in a similar way to the evolution of Earth. However, studies of the Solar System have shown us how the formation and evolution of multi-planet systems are shaped in numerous and often subtle ways by the multifaceted interactions between their components and their surrounding environment. These aspects include the interplay of young planets and the gas and dust of their native protoplanetary disc to that between planets and minor bodies in mature planetary systems. The appearance of the giant planets in the outer Solar System, in particular, has been shown to dynamically excite the planetary bodies within their wide gravitational reach and to be responsible for multiple processes controlling the transfer of volatiles and organics from the outer to the inner Solar System and ultimately to Earth (e.g. Turrini et al. 2011; Turrini & Svetsov 2014; Raymond & Izidoro 2017; Ronnet et al. 2018; Pirani et al. 2019).

Giant planet-driven dynamical transport and impact contamination are processes observed being still active in the Solar System, as exemplified by the flux of Jupiter-family comets presently reaching the inner Solar System (e.g. Morbidelli 2008; Dones et al. 2015; Galiazzo et al. 2016) and impacting Jupiter itself (Taylor et al. 2004; Turrini et al. 2015; Hueso et al. 2018). These processes were significantly more intense at the time of the solar nebula due to the larger population of small bodies existing at the time (Turrini et al. 2011; Raymond & Izidoro 2017; Pirani et al. 2019). The resulting higher impact frequencies between planetesimals allow for the systematic collisional implantation of volatiles and organics in the rocky bodies existing in the inner Solar System with efficiencies comparable to those required by the Earth’s water abundance (Rubin et al. 2005; Turrini & Svetsov 2014; Turrini et al. 2018a; Trigo-Rodríguez et al. 2019), as well as allowing for the widespread fragmentation, compaction and compositional mixing processes recorded by brecciated meteorites (Bischoff et al. 2006, 2018; Trigo-Rodriguez et al. 2006). In parallel, the higher impact frequencies between the excited planetesimals and the giant planet that cause their excitation has been shown to result in the alteration of the atmospheric compositions of the giant planets themselves (see Turrini et al. 2015; Ronnet et al. 2018; Shibata & Helled 2022 for the Solar System and Turrini et al. 2015; Seligman et al. 2022; Sainsbury-Martinez & Walsh 2024 for the generalisation to the case of exoplanetary systems). Finally, carbonaceous planetesimals that do not undergo collisions during their dynamical transport will see their orbits circularised by the aerodynamic drag with the nebular gas and will become part of the population of minor bodies in the inner Solar System (Raymond & Izidoro 2017; O’Brien et al. 2018; Pirani et al. 2019) where they can be accreted by the growing terrestrial planets (Chambers 2001; O’Brien et al. 2006; Alexander 2017).

While the dominant mechanism responsible for the primordial inwards transport of volatiles and organics in the Solar System remains debated – whether it was controlled by Jupiter’s mass growth alone or mainly sculpted by its migration, see Turrini et al. 2018a; O’Brien et al. 2018 and references therein – a growing body of observational evidence supports the role of dynamical transport in shaping the present-day characteristics of the Solar System. In particular, due to the volatile nature of ammonia the presence of ammoniated minerals on Ceres (de Sanctis et al. 2015) and other large main belt asteroids (Kurokawa et al. 2022b) has been interpreted as an indication that they either originated in colder regions of the solar nebula or that accreted material was transported inwards from such outer orbital regions. The presence of dark comets among the near-Earth objects originating from the inner main belt (Taylor et al. 2024) is also consistent with the implantation of cometary material from the outer Solar System.

In parallel, the presence of cometary and water-rich xenoliths in carbonaceous chondrite meteorites (Nittler et al. 2019; Kebukawa et al. 2020) and of apatite minerals in the eucritic meteorites originating from the basaltic surface of Vesta (Sarafian et al. 2013, 2014), one of the oldest bodies known in the Solar System (Bizzarro et al. 2005; Schiller et al. 2011; Consolmagno et al. 2015), as well as in angrite basaltic meteorites (Sarafian et al. 2017), support the notion that this transport process was already active in the first few million years of the life of the Solar System, after the formation of the first generation of planetesimals.

The study of these processes in the case of exoplanetary systems is made even more complex by the fact that many of these systems show indications that their original architectures were altered by phases of dynamical instability and chaotic evolution triggered by planet-planet scattering (Laskar & Petit 2017; Zinzi & Turrini 2017; Turrini et al. 2020; Gajdoš & Vaňko 2023). Even in the case of young planetary systems unaltered by chaos (e.g. Damasso et al. 2024; Mantovan et al. 2024a,b), the degeneracy in the initial conditions caused by the current lack of constraints on the disc-driven migration of their planets hinders these kind of studies. In recent years, ALMA observations have been providing an increasingly detailed view of protoplanetary discs, revealing features in many of them that are suggestive of the presence of forming planets at large distances from their host star (e.g. ALMA Partnership 2015; Isella et al. 2016; Fedele et al. 2017; Long et al. 2018; Andrews et al. 2018). In fact, such observations have recently brought to light the first direct detection of a young exoplanet embedded in its native disc (Keppler et al. 2018; Müller et al. 2018). In parallel, population studies of protoplanetary discs across different star forming regions have revealed a possible resurgence of their dust abundance at ages of 2–3 Myr (Testi et al. 2022). This suggests that many of them could be host to both giant planets and massive planetesimal discs whose dynamical and collisional interactions replenish the disc dust through more intense versions of the collisional cascades that shape debris discs (Bernabò et al. 2022).

HD 163296 is estimated to be 5–6 Myr old (Isella et al. 2016; Wichittanakom et al. 2020) and its protoplanetary disc is currently one of the most well–characterised ones, making it an ideal case study to investigate the dynamical transport of the astrobiologically important volatile elements oxygen (O), carbon (C) and nitrogen (N) outside the boundaries of our Solar System. HD 163296 has been proposed to host at least four giant planets on wide orbits (Isella et al. 2016; Pinte et al. 2018; Zhang et al. 2018; Teague et al. 2018, 2019; Izquierdo et al. 2022; Alarcón et al. 2022; Garrido-Deutelmoser et al. 2023) whose estimated masses, however, fall below the threshold for their direct detections (Guidi et al. 2018; Mesa et al. 2019). Furthermore, ALMA observations showed an unexpectedly high dust abundance in the disc region inwards of the giant planets, where dust transport from the outer disc should be hindered by the dynamical barrier created by the planets themselves (Isella et al. 2016; Liu et al. 2018; Guidi et al. 2022). The planet formation study by Turrini et al. (2019) revealed that both the presence and location of this enhanced dust region are naturally explained by the collisional erosion of planetesimals if HD 163296’s disc hosts an extended and massive population of exocomets dynamically excited by giant planets. The same study highlighted how part of the exocometary population could penetrate the innermost 10–20 au of the protoplanetary disc (although the amount of delivered mass was not quantified) and how a large fraction of the excited planetesimals could be ejected into interstellar space (Turrini et al. 2019). However, the latter estimate was based on a fixed ejection distance linked to the extension of the disc gas without distinguishing between planetesimals ejected into interstellar space or scattered into an extended disc similarly to what is observed in the trans-neptunian region of the Solar System (e.g. Morbidelli 2008).

In this work, we revisit the HD 169236 system so as to explore in more detail the process of dynamical excitation and transport of planetesimals therein, taking advantage of the updated information on the disc and stellar characteristics (Isella et al. 2018; Booth et al. 2019; Kama et al. 2020; Wichittanakom et al. 2020), on the planetary architecture (Pinte et al. 2018; Guidi et al. 2018; Zhang et al. 2018; Teague et al. 2019; Mesa et al. 2019; Izquierdo et al. 2022; Alarcón et al. 2022; Garrido-Deutelmoser et al. 2023) and on the distance of HD 163296 from the Sun (Bailer-Jones et al. 2018, 2021; Gaia Collaboration 2020). Our main goal is twofold: (i) we want to assess the impact of the dynamical excitation on the planetesimals for the transport of the volatile elements C, O, and N, and its implications for the disc itself and the planets it contains (both proposed and hypothetical ones); and (ii) we aim to quantify the efficiency of this class of wide planetary systems in creating interstellar objects like the two that recently crossed the Solar System, 1I/Oumuamua and 2I/Borisov (see Moro-Martín 2018; Bailer-Jones et al. 2020; Pfalzner et al. 2021).

The rest of the paper is organised as follows. Section 2 reports the updated observational constraints on the HD 163296 system and the methods we used to simulate its formation and early evolution. Section 3 presents the results in terms of dynamical transport and ejection of planetesimals and the collisional contamination of possible planets embedded in the disc. In Section 4 we discuss the implications of our results beyond the aspects we directly simulated, while in Section 5 we summarise the conclusions of this work.

2 Observational constraints and numerical methods

We simulated the formation of the giant planets embedded in the protoplanetary disc of HD 163296 and their impact on the surrounding planetesimal disc with the N-body code Mercury-Arχes (Turrini et al. 2019, 2021). The set-up of the n-body simulations builds on the one described in Turrini et al. (2019) accounting for the current constraints on the disc gas mass from Booth et al. (2019) and Kama et al. (2020) and the inclusion of the outermost planetary companion proposed by (Pinte et al. 2018). Specifically, we considered two different disc gas masses and two sets of planetary masses for a total of four combinations of planetary and disc parameters (see Table 1). The stellar mass and age of HD 163296 were always set to 1.95 M_⊙ and 6 Myr following Wichittanakom et al. (2020). Since these stellar, disc and planetary data were computed using the Gaia DR2 distance of HD 163296 from Bailer-Jones et al. (2018), throughout this study we will continue using the DR2 distance of 101.5 pc instead of the more recent EDR3 distance of 100.5 pc from Gaia Collaboration (2020) and Bailer-Jones et al. (2021) for consistency. This choice has no impact on the results of the study as the changes in the masses and the distances within HD 163296’s system would be limited to about 1% (i.e. significantly smaller than the observational uncertainties on the relevant values).

For the disc mass we considered a ‘high disc mass’ (HDM) scenario based on the gas mass estimated by Booth et al. (2019) with ALMA using ¹³C¹⁷O, where the total gas mass is 0.215 M_⊙. We further used a ‘low disc mass’ (LDM) scenario based on the gas mass estimates of Isella et al. (2016) and Kama et al. (2020) using ALMA and Herschel observations, respectively, where the total gas mass was set to 0.05 M_⊙. The protoplanetary disc was modelled adopting the gas surface density of $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{r_0}\right)}^{\gamma }{e}^{[-{\left(\frac{r}{r_0}\right)}^{(2-\gamma )}]}$$(1)

where γ = 0.8 (Isella et al. 2016), r is the radial distance from the central object while the characteristic radius, r₀ = 137.3 au, is obtained by scaling the original characteristic radius of r₀ = 165 au from Isella et al. (2016) by the DR2 stellar distance from Bailer-Jones et al. (2018). The gas surface density Σ₀ is 19.4 g cm⁻² in the HDM scenario and 4.5 g cm⁻² in the LDM scenario. The disc gas mass is assumed to be in steady state and does not decline over time across the simulations, a reasonable approximation given that the mass of gas accreted by the giant planets amounts at most to 8% of the disc mass (i.e. the combination of smallest disc mass and highest planetary masses) and the observed high mass-loss rates due to molecular wind have been interpreted as the onset of the disc dispersal phase (Klaassen et al. 2013). The disc temperature profile on the midplane is $$T(r)=T_0{r}^{-0.5}$$(2)

where T₀ = 220 K after scaling the original temperature profile from Isella et al. (2016) by the DR2 stellar distance (Bailer-Jones et al. 2018).

For the planets, we adopted a conservative approach and focussed on the three giant planets suggested to be responsible for the dust gaps between 40 and 140 au (Isella et al. 2016; Liu et al. 2018; Zhang et al. 2018; Teague et al. 2018, 2019; Izquierdo et al. 2022; Alarcón et al. 2022), along with the fourth outermost giant planet whose presence was proposed based on the observations of the gas dynamics (Pinte et al. 2018). We did not consider the hypothesised presence of an additional giant planet responsible for the tentative detection of a gap at about 10 au (Isella et al. 2018; Zhang et al. 2018), although we did explore the implications of our results for such putative planet (in Sect. 3.3). We also do not consider the proposed scenario invoking the presence of two giant planets instead of one in the dust gas at 50 au and a possible multi-resonant orbital architecture between them and the giant planets at 86 and 137.7 au (Garrido-Deutelmoser et al. 2023). For the planetary masses, we again considered a ‘low planetary mass’ (LPM) scenario, where the masses of the three inner giant planets are from Liu et al. (2018) and that of the outermost giant planet comes from Pinte et al. (2018), along with a ‘high planetary mass’ (HPM) scenario, where the masses of the three inner giant planets are from Teague et al. (2018, 2019) and that of the outermost giant planet is from Pinte et al. (2018). In the LPM scenario the masses originally estimated by Liu et al. (2018) based on the radial extensions of the gaps in the disc dust distribution have been rescaled to the DR2 stellar distance measured by Gaia (Bailer-Jones et al. 2018). These scenarios for the planetary masses encompass the range of values proposed by the different studies (Isella et al. 2016; Liu et al. 2018; Zhang et al. 2018; Teague et al. 2018, 2019; Izquierdo et al. 2022; Alarcón et al. 2022) and are consistent with the upper limits posed by the non-detection of the planets by direct imaging (Guidi et al. 2018; Mesa et al. 2019; Huélamo et al. 2022). The disc and planetary parameters adopted in the four scenarios are summarised in Table 1.

The n-body simulations modelled the dynamical evolution of a disc of planetesimals embedded within HD 163296’s gaseous disc under the effects of the mass growth of the four forming giant planets, of the aerodynamic drag by the disc gas and of the disc gravity. Following Turrini et al. (2019, 2021), planetesimals were modelled as test particles possessing inertial mass and no gravitational mass, meaning that the test particles are affected by the disc gas but do not influence each other nor the giant planets.

The inertial mass is computed assuming a characteristic diameter of 100 km (see Klahr & Schreiber 2016; Johansen & Lambrechts 2017) and bulk density of 1 g/cm³ (see Turrini et al. 2019, 2021). The damping effects of aerodynamic gas drag on the planetesimals were simulated following the treatment from Brasser et al. (2007) with updated drag coefficients from Nagasawa et al. (2019) accounting for both the Mach and Reynolds numbers of the planetesimals. The exciting effects of the disc gravity were simulated based on the analytical treatment for axisymmetric discs by Ward (1981) following Marzari (2018) and Nagasawa et al. (2019). We refer to Turrini et al. (2021) for additional details on this step.

The formation of the giant planets was modelled over two growth phases using the parametric approach from Turrini et al. (2011, 2019). The first phase accounts for their core growth and subsequent capture of an expanded atmosphere (e.g. Bitsch et al. 2015; Johansen et al. 2019; D’Angelo et al. 2021). The planetary mass evolves as $$M_p(t)=M_0+\left(\frac{e}{e-1}\right)(M_1-M_0)(1-e^{-t/{\tau }_p})$$(3)

where M₀ = 0.01 M_⊕ is the initial mass of the core, while M₁ is the final cumulative mass of the core and its expanded atmosphere at the end of the first growth phase (see Turrini et al. 2021 for further discussion), e is the Euler number, t is the time, and τ_p is the duration of the first growth phase. The second phase of mass growth accounts for the runaway gas accretion of the two giant planets, where their mass evolves as $$M_p(t)=M_1+(M_2-M_1)(1-e^{-(t-{\tau }_p)/{\tau }_g})$$(4)

with M₂ being the final mass of the giant planets and τ_g the e-folding time of the runaway gas accretion process.

As in Turrini et al. (2019), we model the formation of HD 163296’s four giant planets as occurring in situ, that is, without them undergoing orbital migration while growing in mass. This simplifying assumption is motivated by the fact that dynamical excitation of the planetesimal disc is dominated by the rapid mass growth and migration of the giant planets during their runaway gas accretion phase (Turrini et al. 2011, 2012, 2019; Raymond & Izidoro 2017; Pirani et al. 2019) and recent studies indicate that the migration of giant planets during such phase is limited (Tanaka et al. 2020; Paardekooper et al. 2023). The inclusion of the migration of the planetary core result in slightly earlier onsets of the dynamical excitation (Pirani et al. 2019).

The adopted values of M₂ for the four planets are those of their estimated masses in the HPM and LPM scenarios. For the three innermost giant planets we adopted values of M₁ = 30 M_⊕ (Lissauer et al. 2009; D’Angelo et al. 2021), τ_p = 1 Myr (Mulders et al. 2021; Bernabò et al. 2022; Lichtenberg et al. 2023), and τ_g = 0.1 Myr (Lissauer et al. 2009; D’Angelo et al. 2021) to simulate their growth by core accretion. The timescale of core formation, τ_p, was chosen to fit the temporal interval in the disc lifetime that is most favourable to the formation of giant planets (Savvidou & Bitsch 2023). For the outermost planet, we set M₁ = 0.3 M_⊕, τ_p = 0.01 Myr, and τ_g = 0.01 Myr to simulate its faster formation by disc instability (D’Angelo et al. 2010; Helled et al. 2014).

During their gas accretion phase all four giant planets form gaps in the disc gas whose widths are modelled as W_gap = C · R_H (Isella et al. 2016; Marzari 2018), where the numerical proportionality factor C = 4 is from Isella et al. (2016) and Marzari (2018) and R_H is the relevant planetary Hill’s radius. The gas density Σ_gap(r) inside the gap evolves over time with respect to the local unperturbed gas density Σ(r) as $${\mathrm{\Sigma }}_{gap}(r)=\mathrm{\Sigma }(r)\cdot e^{[-(t-{\tau }_p)/{\tau }_g]}$$(5)

(Turrini et al. 2021). When planetesimals cross the gaps the effects of the gas on their orbital evolution are computed using the local gas density Σ_gap.

We set the spatial density of planetesimals in the N-body simulations to 1000 particles/au, with the inner edge of the planetesimal disc at 1 au and the outer edge at 200 au (see Turrini et al. 2021 for the discussion of the choice of the disc inner edge). The orbital regions corresponding to the feeding zones of the giant planets (see D’Angelo et al. 2010 and references therein) are not populated by planetesimals to account for the local solids accreted in their cores or trapped within their collapsing gas. The resulting planetesimal discs are populated by about 10⁵ test particles. To accurately reproduce the orbital evolution of planetesimals close to the inner edge of the disc, the time step is set smaller than 5% of the orbital period at 1 au (e.g. Rein & Tamayo 2015) by adopting a value of 15 days. Planetesimals are considered ejected from the protoplanetary disc into interstellar space once their orbits become unbound, namely when their eccentricity becomes equal or greater than one. We recorded snapshots of the evolving architecture of the planetary system every 10⁴ years, keeping track of the original formation region of the planetesimals. The recorded eccentricity, e, and inclination, i, values allow us to quantify the dynamical excitation of the planetesimals as $\sqrt{e^2+{\sin }^2(i)}$ (Petit et al. 2001).

To quantify how the solid mass is redistributed across the protoplanetary disc by the appearance of the four giant planets, we treated each test particle in the N-body simulations as a swarm of real planetesimals. Since constraining the original mass of gas and solids of HD 163296 is a challenging task (Turrini et al. 2019; Booth et al. 2019; Mulders et al. 2021), particularly in light of the high mass loss (Klaassen et al. 2013) and gas accretion rate (Wichittanakom et al. 2020) presently experienced by this protoplanetary disc, we adopted the following simplified approach. The cumulative mass of each swarm was computed integrating the adopted disc gas density profile over a ring 0.1 au wide centred on the initial orbit of the impacting particle, and multiplying the resulting gas mass by the local solid-to-gas ratio. The solid-to-gas ratio, in turn, is a function of the disc metallicity and local disc midplane temperature and is described by the simplified radial condensation profile from Turrini et al. (2023).

The disc metallicity was set to 1.4% (Asplund et al. 2009) under the assumption of solar metallicity for consistency with the characterisation of HD 163296’s stellar parameters by Wichittanakom et al. (2020). The solid-to-gas ratio is 0.5 times the disc metallicity for planetesimals formed at temperatures between 1200 and 140 K, i.e. between the condensation of silicates and that of water. The solid-to-gas ratio grows to 0.75 times the disc metallicity for planetesimals formed at temperatures between 140 and 30 K, namely between the snowlines of water and carbon monoxide, and reaches 0.9 times the disc metallicity for planetesimals formed at temperatures below 30 K, namely, those beyond the carbon monoxide snowline.

The impacts on the proposed four giant planets forming in HD 163296’s disc rings were recorded directly during the n-body simulations. Then the formation region of the impacting particles was used to quantify the mass of the swarm of planetesimals they represent based on the methods described above. To investigate the collisional implications of the dynamical transport of planetesimals to the inner regions of the protoplanetary disc, we also computed the impact probabilities and velocities among the particles that originate within the inner unperturbed planetesimal disc and those originating from the outer excited planetesimal disc using the well-tested statistical collisional methods developed for the study of the asteroid belt (Wetherill 1967; Farinella & Davis 1992; O’Brien & Sykes 2011). The transition between the two planetesimal populations is set at 40 au based on the results of the n-body simulations (see Sects. 3.1 and 3.2). The resulting values are used to estimate the accretional fluxes of impactors on possible undiscovered inner planets hidden within HD 163296’s gas as discussed in Sect. 3.3.

Fig. 1 Dynamical excitation of the planetesimal disc at 0.75, 2 and 5 Myr in the four simulated scenarios. Dynamical excitation accounts for the contributions of both eccentricity and inclination and is defined as $\sqrt{e^2+{\sin }^2(i)}$ (Petit et al. 2001). The colour bar indicates the region in the disc from which the planetesimals originated (in units of au). The dashed grey line at 40 au indicate the division between inner and outer disc. The rectangles in the rightmost panels mark the “fountain” feature we refer to in the main text, i.e. the damping of the orbital eccentricities of the excited planetesimals by aerodynamic drag that creates the implanted planetesimals. The light blue filled circles indicate the location of the four simulated planets in the disc. First row: high disc mass – high planetary mass. Second row: high disc mass – low planetary mass. Third row: low disc – high planetary mass. Fourth row: low disc mass – low planetary mass.

3 Results

In all simulated systems (Figure 1) we can distinguish two planetesimal populations and two temporal phases concerning their dynamical evolution. One population is that of the dynamically excited planetesimals that formed in the disc region where the giant planets reside. The other population is that of the unperturbed (or limitedly excited) planetesimals that formed inwards of the giant planets. As shown by Figure 1 and introduced in Sect. 2, the transition between the two populations occurs at about 40 au. In the following we refer to the first 40 au as ‘inner disc’ and to the region beyond 40 au as ‘outer disc’.

The two populations remain well separated during the first temporal phase that extends until all four giant planets reach their current masses, namely, before 1 Myr in our simulations. In this phase the most dynamically excited region is the one between the third and fourth planets, mainly due to the assumed rapid formation of the outermost planet by disc instability (see Section 2). During this phase the excited planetesimals mostly cross the orbital region of the giant planets and only limitedly penetrate within the inner disc. In this phase we can see a difference in the level of global dynamical excitation between high mass and low mass discs. Specifically, in the HDM scenarios, the stronger effect of the disc gravity adds to the planetary perturbations producing a higher level of dynamical excitation with respect to the LDM cases (see the planetesimals in the rings between the planets at t=0.75 Myr in Fig. 1).

The second temporal phase starts once the giant planets complete their growth. Planetesimals in the outer disc reach very high levels of dynamical excitation (producing eccentricity values above 0.4, see also Fig. 2) and are able to systematically penetrate the innermost 40 au where they start being significantly affected by the gas drag (see Sect. 3.2). During this second phase, therefore, we observe the systematic mixing between the two planetesimal populations within the inner disc (see Sect. 3.1). In the LDM scenarios the effect of gas drag remains limited and only the most excited planetesimals, moving at high relative velocity with respect to the gas, see their excitation damped and have their orbits circularised. In the HDM scenarios, instead, planetesimals that are injected inside the orbit of the second planet (i.e. inwards of 80 au) are affected by a stronger gas drag and see their eccentricity significantly damped over time (the ‘fountain’ feature inwards of 40 au at 5 Myr identified by the box in Fig. 1). By 5 Myr, a fraction of these damped planetesimals has been stably implanted within the inner disc, in some cases reaching inwards of 10 au (see Figs. 1 and 2).

In the rest of this work we identify the excited planetesimals that formed in the outer disc, as being ‘implanted’ in the inner disc when their pericentres are in the inner disc while their orbits underwent circularisation due to the gas drag to the point that their apocentres became less than 40 au (see Fig. 2), and, therefore, became decoupled from the outer disc. This process causes the implanted planetesimals to mix with the original unperturbed population of the inner disc region (see Sect. 3.1 for details). We use the term ‘injected’ to identify the excited planetesimals that formed in the outer disc, have pericentres in the inner disc, and remain on eccentric orbits with apocentres in the outer disc (see Fig. 2), namely where gas drag has not been able to circularise the orbits. These planetesimals can contribute to the compositional evolution of the inner disc through collisional processes (see Sect. 3.2).

The dynamical excitation process scatters a marked population of planetesimals in an extended disc beyond the orbit of planet e and causes the ‘ejection into interstellar space’ of one–tenth to one–fourth of the initial planetesimal disc. The three inner giant planets accrete in total between 4 and 5% of the planetesimals initially embedded in the disc in the HDM scenarios, while in the LDM scenarios their accretion efficiency drops to 1–2% (see Sect. 3.4). The outermost giant planet, due to its faster formation and larger final mass, quickly becomes more efficient in scattering planetesimals rather than accreting them and undergoes no impacts during the simulations.

Fig. 2 Pericentre vs. apocentre plot of the planetesimal disc at the end of our simulations at 5 Myr, where we can clearly distinguish the excited and unperturbed planetesimal populations. Excited planetesimals from the outer disc can penetrate the inner disc down to a few au from the host star. Planetesimals are colour-coded based on their formation region, the colour scale is truncated at 100 au to better resolve the inner disc. The vertical and horizontal dashed lines at 40 au split the disc in the inner and outer region.

3.1 Implantation in the inner disc

Once the system enters its second temporal phase of dynamical excitation, the orbital evolution of the planetesimals injected in the inner disc becomes dominated by gas drag and their architecture is shaped by its damping effects, as exemplified by the presence of the ‘fountain’ features that can be seen in Figure 1. The orbital circularisation caused by this damping process permanently implants the excited planetesimals into the inner disc and mixes them with the local unperturbed planetesimal population, similarly to what has been suggested for the asteroid belt in the Solar System (Turrini & Svetsov 2014; de Sanctis et al. 2015; Raymond & Izidoro 2017; Pirani et al. 2019).

To assess the efficiency and the potential impact of the inwards transport and implantation of planetesimals and astrobiologically important elements, we quantify the ratio between local planetesimals and implanted planetesimals in each 10 auwide ring within the inner disc region (see Table 2). The planetesimals are assumed to be implanted in the rings within whose boundaries their pericentres are located at the end of the simulations. As shown in Figure 1, a number of implanted planetesimals are still characterised by marked eccentricities (e > 0.1), meaning that they can spend part of their orbits outside the boundaries of the ring. Given that HD 163296 still hosts its protoplanetary disc after the timespan covered by our simulations, our approach is equivalent to assuming that the gas drag over the subsequent life of the disc will circularise the orbits of the implanted planetesimals at their pericentres.

As shown in Figure 2 and Table 2, implantation is most effective in the region between 10 and 30 au, spanning two rings that are each 10 au wide, in the case of the HDM scenarios. In the LDM scenarios the peak efficiency is reached between 10 and 20 au. The majority of implanted planetesimals originate beyond the CO line at 60 au with a small fraction originating instead in the region between 40 and 60 au, namely, between the CH₄ and CO snowlines (see Fig. 3). About 50–60% of the planetesimals originating beyond the CO snowline actually come from beyond the N₂ snowline, enriching the local ring composition in nitrogen alongside carbon and oxygen (see Fig. 3).

In the HDM-HPM scenario the giant planets can implant planetesimals between 10–20 au and 20–30 au with similar efficiencies (Table 2), while their less massive counterparts in the HDM-LPM scenario are three times less efficient in implanting planetesimals between 10 and 20 au than between 20 and 30 au (Table 2). The implantation efficiency in the outermost (30–40 au) and innermost (0–10 au) rings proves always negligible in the HDM scenarios. In the LDM scenarios the most enriched ring is that between 10 and 20 au, which sees its mass increased by 3 and 1% in the HPM and LPM scenarios, respectively. Differently from the HDM case, however, in the LDM-HPM scenario also the innermost ring (0–10 au) sees its mass increased by 1.3% thanks to the interplay between the stronger planetary perturbations and the weaker gas drag that allows planetesimals to penetrate deeper in the inner disc. The comparison of the implantation efficiencies in the different orbital regions as a function of the disc and planetary masses is shown in Fig. 4. Specifically, the solid lines in Fig. 4 show the ratios between the implantation efficiencies for the two different sets of planetary masses while keeping the disc mass constant, while the dashed lines show the ratios for the two different disc masses while keeping the planetary masses constant. As can be immediately seen from the slopes of the curves in Fig. 4, increasing the disc mass enhances the implantation efficiency in the part of the inner disc closer to the giant planets, while increasing the planetary masses enhances the implantation efficiency closer to the host star.

The enrichment in mass due to implantation, while moderate even in the most efficiently implanted ring, is comparable to that estimated to have occurred in the inner Solar System and the asteroid belt following the formation of Jupiter and Saturn (1–10% Raymond & Izidoro 2017; Ronnet et al. 2018) while affecting an orbital region an order of magnitude larger. In other words, HD 163296’s protoplanetary disc is characterised by higher degrees of dynamical excitation and compositional remixing than the Solar System. The higher abundance of N in the implanted planetesimals (Öberg & Bergin 2021) makes the enrichment in this element much more significant than what the implanted mass suggests.

More generally, due to the important roles of N₂ and CO as carriers of nitrogen and carbon (Öberg & Bergin 2021), the enrichment in these elements of the inner disc region (0–40 au) proves always larger than the one in H₂O and volatile elements estimated as having occurred in the asteroid belt and the inner Solar System due to cometary impactors (0.1–0.01%, see Turrini & Svetsov 2014; Sarafian et al. 2017; Turrini et al. 2018a, and references therein). Finally, the mass implanted between 10 and 30 au in the HDM scenario and between 0 and 20 au in the LDM scenario is comparable in magnitude or greater than that estimated for the ‘late veneer’ that affected the planetary bodies in the asteroid belt and the inner Solar System (see Day et al. 2012, 2016; Turrini et al. 2018a, and references therein).

Fig. 3 Same as Fig. 2 with the excited planetesimals being colour–coded according to the compositional region they originate from. This compositional gradient is based on HD 163296’s temperature profile and on the condensation temperatures as ices of the main volatile molecules.

Fig. 4 Relative efficiency between the HPM and LPM scenarios (solid curves) and between the HDM and LDM scenarios (dashed curves) in implanting planetesimals in the inner disc. The orange colour identifies the high-mass disc when comparing the efficiencies of the HPM and LPM scenarios, and the high-mass planets when comparing the efficiencies of the HDM and LDM scenarios. The blue curves identify the LDM and LPM scenarios, respectively. High relative efficiencies in the curves can still be associated with limited absolute efficiencies (given in Table 2). Increasing the disc mass enhances the implantation efficiency in the part of the inner disc closer to the giant planets, while increasing the planetary masses enhances the implantation efficiency closer to the host star.

3.2 Injection in the inner disc

Alongside the effect of dynamical implantation we also evaluated the efficiency of dynamical injection and collisional delivery of planetesimals to the inner regions of the disc, similarly to the case of the cometary flux in the inner solar system (Turrini & Svetsov 2014; Turrini et al. 2018a). Since the dynamical injection process depends only on the orbital eccentricity of the excited planetesimals, it involves significantly more bodies than the implantation process and is particularly efficient with respect to the highest planetary masses considered here (see Table 2).

While lower disc masses result in decreased damping efficiency of gas drag and should allow for the injection of a larger population of high eccentricity planetesimals, we observe (as shown in Table 2) this behaviour only with the HPM cases. For the LPM cases the HDM scenario is always characterised by higher injection efficiencies than the respective LDM counterpart due to the effects of the stronger disc gravity. The exciting effects of the local disc gas increase the number of planetesimals that cross the orbital region of the giant planets and get their eccentricity increased during close encounters (see Fig. 1 at 0.75 Myr).

The injection efficiency peaks in the outermost ring of the inner disc, i.e. between 30 and 40 au, and decreases monotonically when moving to rings closer to the star. The injection efficiency decreases a factor of two faster for the LPM scenarios than for the HPM scenarios. Specifically, in the LPM scenarios it drops by a factor of 10–15 moving from the ring between 30 and 40 au to the ring between 10 and 20 au, while over the same orbital range the efficiency of the HPM scenarios drops only by a factor of 4–6. The orbital region between 20 and 40 au is crossed by a mass flux of planetesimals representing 10–30% of the local mass of planetesimals.

As shown in Fig. 2, excited planetesimals can reach the disc region inwards of (or very close to) 10 au in all scenarios. In the HDM-HPM scenario, the mass of injected planetesimals reaches up to 1% of the local planetesimal population, while in all other scenarios, the injected mass is always less than this value (see Table 2). Notwithstanding, the efficiency of the injection process in HD 163296 proves comparable to the efficiency characterising the scattering process in the inner Solar System triggered by the formation of Jupiter and Saturn (Raymond & Izidoro 2017; Ronnet et al. 2018), while affecting an orbital region that is wider by an order of magnitude.

3.3 Collisional enrichment and contamination of possible inner planets

As introduced in Sect. 2, from the analysis of the observations from the DSHARP ALMA survey (Andrews et al. 2018) Zhang et al. (2018) and Isella et al. (2018) report evidence for an additional gap at 10 au within HD 163296’s protoplanetary disc, not seen in the previous works on the system (e.g. Isella et al. 2016). In particular, Zhang et al. (2018) postulate that such a ring can be the result of a new planet with mass ranging between 0.35–0.7 M_J embedded in the disc. While we do not account for this possible planet in our simulations, in light of the estimated efficiency of the injection process this result prompted us to look into the collisional enrichment and contamination effects of the planetesimal flux crossing the inner disc on possible planetary bodies existing within 40 au, following Turrini & Svetsov (2014) and Turrini et al. (2015, 2018a). Specifically, following Turrini et al. (2015) we aim both to quantify the increase in the abundance of the volatile elements (enrichment) and to assess whether this enrichment can significantly alter or mask the compositional signatures of the formation process (contamination). To calculate the numbers of potential impacts on these assumed targets, we focussed on the orbital architectures of the planetesimal discs in the final snapshot of our simulations and computed the intrinsic impact probability (see Sect. 2 and Wetherill 1967; Farinella & Davis 1992) between all possible target-impactor pairs.

We considered as targets the orbits and positions of all planetesimals originating within 40 au as a proxy for planetary–sized objects, while we considered as impactors all the planetesimals originating beyond 40 au whose pericentres reach inwards of 40 au. We focussed on six impact scenarios based on templates from both the Solar System¹ and exoplanets² spanning the mass range from small planetary embryos to fully-formed massive planets. Specifically, we assumed targets with lunar-like (0.01 M_⊕), Earth-like (1 M_⊕), Neptune-like (17 M_⊕) and Jupiter-like (318 M_⊕) masses as well as targets modelled after the super-Earth GJ 486b (M_t=3 M_⊕, R_t=1.34 R_⊕, v_esc=16.7 km/s) and the mini-Neptune TOI-1260b (M_t=8.56 M_⊕, R_t=2.41 R_⊕, v_esc=21.1 km/s), based on the classification from Fulton et al. (2017). We focussed our evaluation on a temporal window of 1 Myr, mamely the difference between the timespan covered by our simulations and the estimated age of HD 163296 (Wichittanakom et al. 2020), to account for the possible slower or later formation of the putative inner planets.

To compute the impact fluxes from the intrinsic impact probabilities we used the collisional formula (see O’Brien & Sykes 2011, and references therein): $$N_i=P_i\times {(R_{t,eff}+R_{pl})}^2\times \mathrm{\Delta }t\times N_{pl}$$(6)

where N_i is the number of impacts on the target body, P_i is the intrinsic impact probability (see Sect. 2 and Wetherill 1967; Farinella & Davis 1992, the value includes the π term of the cross section), R_pl is the projectile radius in km (here fixed at 50 km), Δt the temporal interval considered (here 1 Myr), N_pl is the number of planetesimals that are injected into the inner disc, and R_{t,e f f} is the target’s effective radius in km defined as $$R_{t,eff}^2=R_t^2\times (1+{\left(\frac{{\nu }_{esc}}{{\nu }_{imp}}\right)}^2)$$(7)

where R_t is the target’s physical radius, ν_esc and ν_imp are the escape and impact velocity respectively. The tethe injecrm $1+{\left(\frac{{\nu }_{esc}}{{\nu }_{imp}}\right)}^2$ is the gravitational focussing factor (Safronov 1972).

The resulting numbers of impacts were multiplied by the inertial mass of the colliding planetesimals and normalised depending on the mass of the target body. In the case of the mini-Neptune-like, Neptune-like and Jupiter-like targets, the mass flux of impactors is normalised to the mass of heavy elements contained in an atmospheric shell of solar metallicity accounting for 20% of 1.5 the target mass following Turrini et al. (2015). The heavy elements in these atmospheric shells amount to about 0.02 M_⊕ for the mini-Neptune targets, 0.05 M_⊕ for the Neptune-like targets and 0.9 M_⊕ for the Jupiter-like targets. These values are used to compute the amount of heavy elements required to reach an atmospheric enrichment that is three times solar as, as in the case of Jupiter (Atreya et al. 2018), which was used to compare the enrichments produced by collisional contamination.

For the lunar-sized, Earth-sized and, super-Earth-sized targets, the cumulative mass flux of planetesimals is instead normalised to the whole target mass. As reference value to compare the magnitude of the resulting late accretion with Solar System analogues, we adopt the water mass fraction estimated for the Earth (5 × 10⁻⁴ planetary masses, Morbidelli et al. 2000). We need to emphasise, however, that Earth’s water enrichment does not translate into an unequivocal late accretion mass value due to the uncertainty on the source of Earth’s water. Specifically, Earth’s water enrichment requires the accretion of a mass of planetesimals twice as large if these planetesimals are cometary in nature and have half their mass as water ice (e.g. Turrini et al. 2018b, and references therein). If the source of Earth’s water was instead linked to carbonaceous asteroids, whose water content amounts to about 10% in mass (e.g. Turrini et al. 2018b, and references therein), the mass of planetesimals required to produce Earth’s water abundance would be ten times as large. Because of this uncertainty and the lack of detailed information on the chemical setup of HD 163296’s disc, Earth’s water content is adopted as a convenient reference but the outcomes of the comparison should not be over-interpreted.

The effects of the collisional contamination by impacting planetesimals are summarised in Fig. 5³. In the HDM scenarios lunar-like, Earth-like and super-Earth targets can become enriched in ices from the outer disc to levels comparable to Earth’s water enrichment; whereas in the LDM scenarios, the enrichment of this class of targets is always smaller than this reference value by at least a factor of a few (see left panels in Fig. 5). The highest enrichments are reached by target planets beyond 10–15 au in the HPM scenarios and beyond 15–20 au in the LPM scenarios, with the magnitude of the enrichment increasing proportionally to the gravitational cross-section of the target planet (e.g. super-Earths experience higher enrichments than Earth-like and lunar-like targets). In these orbital regions, the amount of outer disc planetesimals accreted by these target planets falls between the cases of cometary and carbonaceous asteroidal impactors as the source of Earth’s water discussed above (Turrini et al. 2018b).

Also, in the case of the three more massive targets (right panels in Fig. 5), the larger gravitational cross-section of the Jupiter-like targets allows for a greater enrichment than their Neptune-like and mini-Neptune counterparts. Specifically, in the HDM scenarios the atmospheric enrichment of the Jupiter-like targets is always super-solar beyond 10 au and can grow to about 9x (i.e. Saturn’s atmospheric enrichment, Atreya et al. 2018) between about 15 and 30 au as well as around 40 au. The atmospheres of Neptune-like planets can become three to four times more enriched in heavy elements with respect to the solar composition (i.e. as enriched as Jupiter) between 10 and 30 au in the HPM cases and between 20 and 30 au in the LPM cases (see Fig. 5), as well as around 40 au in both cases. The atmospheres of mini-Neptune targets are the least enriched yet they range between two and three times solar beyond 15 au in the HPM case and 20 au in the LPM case. In the LDM scenarios, the atmospheres of these massive planets are characterised by mostly solar abundances (see Fig. 5), with only Jupiter-like target showing limited super-solar enrichments.

These results highlight how, in the HDM scenarios, the atmospheric enrichment of possible inner giant planets can be so marked that the resulting collisional contamination can significantly alter or overwrite the original compositional signatures of their formation process. In the LDM scenarios, on the other hand, the compositional signatures left by the planet formation process on giant planets appear to be preserved.

Fig. 5 Enrichment of possible undiscovered planets in the inner disc by the impacts of planetesimals that formed in the outer disc for all four simulation scenarios (with blue we denote the HDM-HPM scenario, with green the HDM-LPM one, with orange the LDM-HPM one and with yellow the LDM-LPM one). Top row: normalised mass fraction of accreted material onto lunar-like, Earth-like and super-Earth targets, with the horizontal dashed line marking the water mass fraction on Earth (5×10−4, Morbidelli et al. 2000) shown from left to right. Bottom row: atmospheric enrichment of mini-Neptune, Neptune-like an Jupiter-like targets assuming a solar-like metallicity atmosphere, with the horizontal dashed line marking Jupiter’s 3× solar atmospheric metallicity (Atreya et al. 2018), shown from left to right.

3.4 Collisional enrichment and contamination of the known giant planets

In parallel to the enrichment and contamination of possible undiscovered inner planets, we investigated the effects the excited planetesimals have on the very giant planets embedded in the gaps of HD 163296’s disc that are at the origin of their dynamical excitation. Across our simulations, the first three planets (b, c and d) are hit by a flux of impactors ranging from a few hundred to a couple of thousand test particles. The impacting particles originated at distances spanning the locations of the H₂S, CH₄, CO, and N₂ snowlines, located at semimajor distances of about 19, 47, 61, and 71 au for the reconstructed temperature profile of HD 163296’s disc midplane (see Fig. 3), and enrich the forming giant planets with heavy elements.

The resulting compositional enrichment is significant in the HDM scenarios, where the total enrichment can reach values as high as 5 times the solar metallicity, with as much as 11 M_⊕ of heavy elements being added to a single planet. In the LDM scenarios, on the other hand, none of the planets accrete more that the equivalent of an Earth mass of planetesimals; this means that the enrichment effect is negligible when assuming a solar metallicity for HD 163296’s host star and its disc. Planet e suffers no impact in all scenarios, as its faster formation timescale and larger mass make it more efficient in scattering the planetesimals that undergo close encounters with it than in accreting them.

Since the planetesimal impact rate is not constant during the growth process of the giant planets but peaks at the beginning of the runaway gas accretion process (Podolak et al. 2020; Turrini et al. 2021), the previous bulk enrichments may not be reflected in the atmospheric composition of the planets unless their envelopes are well mixed. We therefore took advantage of the information on the times of the impacts recorded by our simulations to quantify the atmospheric contamination due to late accretion, to explore the scenario of inhomogeneous planetary envelopes. We focussed on the contamination affecting the outermost 20% of the planetary mass, assumed to represent the molecular gas shell (Turrini et al. 2015); hence, we examined the impacts occurring after 0.16 Myr since the onset of the runaway gas accretion (or, equivalently, 1.16 Myr since the beginning of the growth process in our simulations).

We once again find that in the HDM scenarios the atmosphere of the giant planets can be enriched by as much as 2–3 M_⊕ of heavy elements, resulting in their super-stellar metallicities. In the LDM scenarios the atmospheric contamination by planetesimal impacts always results in the delivery of less than 1 M_⊕ of heavy elements and therefore produces negligible effects. In Table 3 we summarise the above results for both the bulk envelope and atmospheric enrichment cases.

The above results, together with those discussed in Sect. 3.3, highlight how the interplay between the formation process of multi-planet systems and the protoplanetary discs in which they are embedded can play a marked role in altering the primordial compositional signatures of the native environments in the planetary atmospheres, even in absence of orbital migration, and how collisional contamination is an important process in shaping multi-planet systems born in massive discs.

3.5 Ejection as interstellar objects

As discussed in Section 2, we considered as ejected from HD 163296’s system all planetesimals whose orbits become unbound to the host star, namely, their eccentricity is equal to 1 or greater. We also estimated the fraction of planetesimals that are scattered on high eccentricity orbits with semimajor axis beyond 260 au, namely, beyond the orbit of the outermost planet. As shown in Table 4, we find that significant fractions of the planetesimals originally present within HD 163296’s protoplanetary disc are ejected into interstellar space, confirming the findings of Turrini et al. (2019). About half as many planetesimals as those ejected into interstellar space are scattered into an extended disc beyond the orbit of the fourth giant planets (see Table 4), similar to the so-called scattered disc in the trans-Neptunian region of the Solar System (e.g. Morbidelli 2008).

The key factor in controlling the ejection efficiency is the mass of the planets (see Table 4): in the HPM scenarios about one fourth of the planetesimals initially present in the disc are ejected, while in the LPM scenarios this value drops to about one tenth. The disc mass has a smaller impact on the ejection efficiency driven by the exciting effect of the disc gravity, which proves dominant with respect to the damping effect of gas drag.

Specifically, the constructive interplay between the planetary perturbations and the disc gravity leads to slightly higher ejection efficiencies in the high disc mass scenarios (see Table 4). The impact of the disc gravity is more pronounced in the LPM scenarios, where the fractional increase in ejection efficiency when moving from the low disc mass to the high disc mass is four times larger than in the case of the HPM scenarios (18% instead of 4%).

In Figure 6 we show how the ejection efficiency and the source regions of the ejected planetesimals vary as a function of the disc and planetary masses. The ejection efficiency is expressed as the fraction of the planetesimals originally contained into each 10 au-wide ring that is ejected into interstellar space. The regions experiencing the highest ejection rates for all combinations of disc and planetary masses are those between the second and third and between the third and fourth giant planets. The ring of planetesimals between the first and the second planets is significantly less depleted in the HPM scenarios than in the LPM ones (see Figure 6).

Under the assumptions adopted for the computation of the mass of the planetesimal discs discussed in Section 2, we find that the LDM scenarios can eject into the interstellar space between 17 and 37 M_⊕. In the HDM, on the other hand, we see a factor of 5 increase in the ejected mass, that ranges between 88 and 168 M_⊕. Note, however, that these values provide only general indications of the efficiency of the process in producing interstellar objects, as we do not know the original mass and distribution of solids nor the planetesimal formation efficiency across the disc. The comparison between Figs. 3 and 6 points to the majority of ejected planetesimals being N-rich with an almost stellar composition with the remaining ones being rich in C and O, making them as rich in volatiles and organics or even richer than the most pristine comets observed in the Solar System (see Mumma & Charnley 2011; Altwegg et al. 2019, and references therein).

Fig. 6 Planetesimal ejection efficiency across the outer disc in the four disc-planets scenarios (from left to right: HDM-HPM, HDM-LPM, LDL-HPM and LDL-LPM). The ejection efficiency is expressed as the fraction of planetesimals originally present in each 10 au-wide ring that end up being ejected from HD 163296’s disc.

4 Discussion

In this section we first discuss the caveats on our results due to the model assumptions, followed by the physical processes and effects that where not modelled directly, but still extend the impact of the dynamical and collisional processes we simulated for the compositional setup of HD 163296’s protoplanetary disc and its planets.

4.1 Model assumptions and caveats

In this work, we assumed a planetesimal disc whose extension matches that of the observed dust disc and characterised by uniform planetesimal formation efficiency. We further assumed that the planetesimal disc was already present at the beginning of the giant planet formation process. Finally, we did not account for the enrichment of the disc gas in heavy elements due to the sublimation of ices from the inwardly drifting pebbles possibly supporting the core growth process.

While the assumption of matching extensions of the planetesimal and dust discs is compatible with constraints from the Solar System (Kretke et al. 2012), the radial dependence of the formation efficiency of planetesimals across the disc is currently largely unknown. Markedly lower planetesimal formation efficiencies in the outer disc regions with respect to the inner disc regions would result in limited or negligible impact of the inwards transport and enrichment processes. It would also result in a sparser scattered disc and a smaller production of interstellar objects. However, the large dust abundances observed inwards of the inner giant planet (Isella et al. 2016; Guidi et al. 2022) and their proposed origin as second-generation dust produced by high-velocity planetesimal collisions argue against widely different planetesimal formation efficiencies between inner and outer disc regions.

The presence of significant amounts of planetesimals in the disc at the onset of giant planet formation is consistent with recent results on the early conversion of dust into planetesimals during the Class 0/I stages (Cridland et al. 2022). As we discuss in Sect. 4.4, our results remain qualitatively valid also in case of a prolonged planetesimal formation process throughout the giant planet formation process. The main impact on our results in such a scenario is limited to the temporal evolution of the transport and enrichment processes.

The effects of the pebble flux on the disc gas metallicity and the planetary compositions would mainly influence our results on the enrichment of the giant planets present in the disc. Specifically, the enhanced gas metallicity produced by the sublimation of the ices from the drifting pebbles would translate in an higher metallicity of the envelopes of the giant planets (which we assumed to be characterised by solar metallicity) and lower enrichments by planetesimal accretion. However, for typical disc parameters (Rosotti 2023) the enrichment of the disc gas mainly occurs in a limited radial region inwards of the snowlines (Booth & Ilee 2019). This means that this process would be relevant only for the innermost giant planet at 50 au and possible undiscovered giant planets close to the snowlines in the inner disc (see Fig. 3). In the case of solid planets, the only impact on our results would be the reduction of the enrichment by a factor equal to the ratio between the mass in pebbles and the total mass of solids (pebbles and planetesimals) in the disc.

4.2 Interplay between dynamical transport and geophysical evolution of the planetesimals

In modelling the composition of the planetesimals and discussing the implications of their dynamical transport for the protoplanetary disc, we implicitly assumed that all planetesimals preserved their original composition, set by their formation regions in the disc midplane, and their initial budget of volatiles and organics. The meteoritic data, however, reveal that the Solar System experienced the formation of multiple generations of planetesimals and that the earliest generations experienced extensive melting and volatile loss during their differentiation and geophysical evolution (e.g. Scott 2007; Lichtenberg et al. 2023, and references therein). Furthermore, meteoritic data also seem to indicate that the earliest generations of planetesimals formed in the inner Solar System, with planetesimal formation in the outer Solar System having being delayed by about 1 Myr (Lichtenberg et al. 2023).

In parallel, observational data on the evolution of the dust abundances in discs across the Class 0/I and Class II phases (Tychoniec et al. 2020; Mulders et al. 2021; Bernabò et al. 2022) and planetesimal formation studies (Cridland et al. 2022) are pushing the beginning of the planet formation process closer to the beginning of the star formation process. The combination of these indications suggests that the planetesimals present in HD 163296’s inner disc may have experienced various degrees of volatile and organics loss with respect to their counterparts in the outer disc. In such a scenario, the impact of the dynamical injection and implantation of outer disc planetesimals for enriching HD 163296’s inner disc in N, C, and O would be significantly enhanced with respect to the processes discussed in Sect. 3.

4.3 Impact of dynamical transport on HD 163296’s habitable zone

Carbon, oxygen, and nitrogen are three of the most important heavy elements representing the building blocks of life (e.g. Herbort et al. 2024). Because of their different volatility these elements will be present with different abundances in the planetesimals populating HD 163296’s inner disc. Specifically, carbon will be more depleted than oxygen and nitrogen will be more depleted than carbon, with nitrogen being significantly depleted in solids up to the N₂ snowline (Öberg & Wordsworth 2019; Öberg & Bergin 2021; Turrini et al. 2023). As we explained in Sects. 3.1 and 3.2, the inwards transport of outer disc planetesimals enriches the inner disc in all three elements.

This enrichment in astrobiologically important elements can play a significant role in shaping HD 163296’s habitable zone, where the probability that planetesimals experienced major volatile loss during their geophysical evolution is high (see Sect. 4.2 and Lichtenberg et al. 2023 and references therein). A simple calculation of HD 163296’s habitable zone (Kasting et al. 1993; Kopparapu et al. 2013) using the current luminosity of the star (16 L_⊙; Millan-Gabet et al. 2016) gives a ring spanning from approximately 3.8–5.5 au. Considering the physical properties of a Main Sequence star with the same mass and metallicity (solar) as HD 163296, the habitability zone can expand to around 9 au (see Figure 6 from Danchi & Lopez 2013).

In two of our simulated scenarios (HDM–HPM and LDM–HPM), we find that as much as 1.3% of the mass contained in the innermost ring (10–0.1 au) is made up of planetesimals rich in N, C and O that formed at distances greater than 40 au. Even if we assume the total loss of volatiles by the planetesimals originally formed in the habitable zone, the enrichment in volatile elements in this region proves greater than that estimated for the Earth and the asteroid belt in terms of water (a few 10⁻⁴, see Morbidelli et al. 2000; Turrini & Svetsov 2014; Turrini et al. 2018a, and references therein). Based on the study of the Earth, the enrichment in N is particularly important in light of its low solubility into magma oceans (Kurokawa et al. 2022a). Specifically, Kurokawa et al. (2022a) argue that if the Earth’s current N budget is primordial and inherited from its formation process, Earth’s initial abundance of N should have been an order of magnitude higher. The delivery of N-rich planetesimals to the habitable zone could therefore support the formation of N-dominated atmospheres like that of the Earth.

Finally, as shown in Fig. 3 the inwards transport of planetesimals can also affect planetesimals enriched in sulphur by the condensation of H₂S, but data from meteorites (Lodders 2010; Palme et al. 2014) and protoplanetary discs (Kama et al. 2019; Rivière-Marichalar et al. 2020, 2021; Le Gal et al. 2021) suggest that most S is incorporated into refractory materials close to the host star. As a result, the enrichment in S due to the dynamical transport of planetesimals is plausibly limited.

4.4 Extended impact flux and dynamical transport in the case of a pebble-rich disc

While the focus of this study is on the planetesimals, the process of dynamical transport we discuss here applies also to the case of dust-dominated or pebble-dominated scenarios for HD 163296’s disc. Specifically, when the growing planets reach the so-called pebble isolation mass, dense rings of pebbles accumulate at the pressure bumps outside each planet orbit (Morbidelli & Nesvorny 2012; Lambrechts & Johansen 2014). These pebble-dense rings provide favourable sites for forming planetesimals by streaming instability (Johansen & Lambrechts 2017) and can continue to grow in mass until an outer planet reaches its pebble isolation mass, blocking the incoming flux of pebbles, or the gas begins to dissipate typically within 10–20 Myr (Hernández et al. 2007; Fedele et al. 2010; Ribas et al. 2014).

Once planetesimals form in these rings, they get dynamically excited and undergo the same orbital evolution of the dynamical tracers in our simulations (Eriksson et al. 2021). This means that if HD 163296’s disc has been characterised by a continuous process of planetesimal formation in its outer regions (is also proposed also for the outer Solar System, Lichtenberg et al. 2023; Neumann et al. 2024), the dynamical injection of the planetesimals and the collisional contamination of the planets embedded in the disc can evenly spread across the whole life of the system – instead of reaching their peak intensity when the giant planets undergo their runaway gas accretion. This scenario likely applies to the two outermost rings, as the dust production by planetesimal impacts in these orbital regions does not appear to be efficient enough to explain the local dust solely as second-generation dust (Turrini et al. 2019).

Furthermore, if pebbles are still present in the system when the gaseous disc dissipates an additional phase of dynamical evolution occurs. Once the gas has dissipated, the pebbles do not feel the pressure bump and the damping effects of the gas any longer. As a consequence, the pebble rings are quickly excited by the giant planets and the pebbles start to behave dynamically as planetesimals. An additional flux of ice-rich pebbles is expected, which would reach the inner regions of the planetary system like discussed in Sect. 3. This third wave of contaminating bodies is expected to extend in time the impacting flux on existing planets, increasing their budget of carbon, oxygen, and nitrogen.

Differently from planetesimals, however, due to the small sizes of the pebbles the ices and organics they contain will undergo prompt sublimation due to the higher temperatures they experience during their transport inwards. The capability of this pebble flux to transport volatiles and organics to the inner disc will therefore depend on the ratio between the time they need to reach thermal and chemical equilibrium with the surrounding environment and their flight time before being accreted. This late flux of pebbles will likely be more effective in contaminating the observed giant planets, which will be closer to their source regions, than any possible inner planet.

4.5 Interplay between dynamical transport and pebble accretion in the inner disc

As discussed by Turrini et al. (2019), the dynamical excitation and transport of the planetesimals will result in their super-sonic motion with respect to the gas, which, in turn will result in their heating and thermal ablation (Tanaka et al. 2013). In their study focussed of the Solar System and the disc of HL Tau, Eriksson et al. (2021) argued that the ablation of scattered planetesimals may replenish the pebble population in the inner disc regions, supporting the process of planetary growth by pebble accretion by local planetary bodies (Johansen & Lambrechts 2017).

In parallel, Turrini et al. (2019) showed how the collisional cascade triggered by excited planetesimals can produce tens of Earth masses of dust and explain the abundant dust population observed inwards of HD 163296’s first giant planet (Isella et al. 2016). The study of the ejecta size distribution after Hayabusa 2’s impact experiment on asteroid Ryugu point to the ejected dust particles ranging in size from one mm to several decimetres, with a characteristic size of about one cm (Wada et al. 2021).

Collisional debris and second-generation dust sharing a similar size distribution will behave dynamically (i.e. as pebbles) and also support the planetary growth process by pebble accretion (Turrini et al. 2023) in the inner disc of HD 163296, overcoming the obstacle posed to the inwards drift of pebbles from beyond 40 au by the barrier effect of the four giant planets. As a result, the dynamical excitation and inwards transport of the outer disc planetesimals by the giant planets can support the formation of inner undiscovered planets (as those discussed in Sect. 3.3) and, more generally, extend the duration of the pebble accretion and streaming instability processes in the inner disc.

5 Conclusions

Over the past decade, a growing body of study has argued that the appearance of giant planets in protoplanetary discs naturally triggers an intense phase of dynamical excitation of the surrounding planetesimal disc, resulting in the dynamical and collisional transport of volatile elements across different orbital regions and planetary bodies (Turrini et al. 2011, 2015, 2018a,b, 2019; Turrini & Svetsov 2014; Raymond & Izidoro 2017; Ronnet et al. 2018; Pirani et al. 2019; Bernabò et al. 2022; Seligman et al. 2022; Shibata & Helled 2022; Sainsbury-Martinez & Walsh 2024). In parallel, observational and laboratory studies have been providing more and more evidence supporting the importance of these processes in shaping the early history of the Solar System (Sarafian et al. 2013, 2014, 2017; de Sanctis et al. 2015; Bischoff et al. 2018; Trigo-Rodríguez et al. 2019; Nittler et al. 2019; Kebukawa et al. 2020; Kurokawa et al. 2022a; Taylor et al. 2024).

In this work we use the HD 163296 system as our case study, as it is one of the best characterised protoplanetary discs (Isella et al. 2016, 2018; Booth et al. 2019; Kama et al. 2020; Wichittanakom et al. 2020) suggested to be host to four or more giant planets (Isella et al. 2016; Liu et al. 2018; Pinte et al. 2018; Guidi et al. 2018; Zhang et al. 2018; Teague et al. 2018, 2019; Mesa et al. 2019; Izquierdo et al. 2022; Alarcón et al. 2022; Garrido-Deutelmoser et al. 2023) and a dynamically and collisionally excited planetesimal disc (Turrini et al. 2019). Our goal is to investigate the impact of the dynamical and collisional transport processes in infant planetary systems with giant planets on wide orbits, such as those revealed by ALMA surveys of protoplanetary discs (e.g. ALMA Partnership 2015; Isella et al. 2016; Fedele et al. 2017; Long et al. 2018; Andrews et al. 2018).

We explored four end-member scenarios that span the possible combinations of the estimated disc and planetary masses. These scenarios have allowed us to assess the impact of these physical parameters on the dynamical evolution of the system as a whole. We find that the appearance of the giant planets in the host disc always results in a large-scale dynamical excitation of the planetesimals. The intensity of the dynamical excitation process has been shown to be controlled by the masses of the giant planets yet we find that the disc mass has smaller but non-negligible effects, with the excitation effect of the disc gas proving more important than that of the damping of gas drag in the outer disc regions. The dynamical excitation of the planetesimal disc results in the creation of a massive population of exo-comets that transport volatile elements to the inner disc regions, enriching existing planets and their atmospheres in these astrobiologically important elements, and producing a large population of interstellar planetary objects.

In particular, we find that:

1. The giant planets in the outer disc deliver planetesimals to the inner disc both through dynamical injection and implantation. Injected planetesimals periodically cross the inner disc on eccentric orbits, while implanted planetesimals become decoupled from the outer disc to permanently reside in the inner disc. HD 163296’s giant planets between 50 and 260 au are capable of delivering large amounts of planetesimals inwards of 40 au.

2. Dynamical implantation and injection influence different but overlapping regions of the inner disc. Implantation is more efficient between 10 and 30 au, while injection is more efficient between 20 and 40 au. Implantation can increase the local mass of the affected ring up to ≈10%, while injection can temporarily increase it up to ≈50%. These values are comparable or higher than those characterising the same process in the Solar System (Turrini & Svetsov 2014; Raymond & Izidoro 2017; Ronnet et al. 2018) but the affected area is an order of magnitude larger.

3. Higher planetary masses result in excited planetesimals penetrating deeper into the inner disc and allow for their implantation in the innermost 10 au. Higher disc masses enhance the excitation in the outer disc but the stronger gas drag they exert circularises the orbits farther away from the star, favouring the region between 20 and 30 au in the case of HD 163296.

4. The dynamical transport process results in the collisional enrichment of possible planets populating the inner disc in the volatile elements C, O, and N (as well as possibly, S). For higher disc masses, solid planets spanning from lunar to super-Earth masses and residing beyond 10 au are enriched to higher or comparable levels than the Earth; whereas for lower disc masses, the enrichment is on average an order of magnitude smaller. Giant planets ranging from sub-Neptunian to Jovian masses and residing in the same orbital region can have their atmospheres enriched from two to nine times the stellar metallicity proportionally to their mass.

5. The dynamical transport process can affect also the very giant planets that trigger it. For the higher disc masses, planetesimal accretion significantly enhances the bulk metallicity of planets b and c and moderately that of planet d. The atmospheric metallicity of these planets is increased moderately to scarcely, respectively. For the lower disc masses, the impact of planetesimal accretion on the planetary metallicity is generally negligible. Due to our assumptions on its faster formation, planet e is never affected by planetesimal accretion.

6. The dynamical excitation of the planetesimals scatters a significant fraction of them (up to 40%) beyond the orbit of the outermost planet. In all simulated cases, about one-third of the scattered planetesimals remain bound to HD 163296 forming an extended and excited planetesimal belt as in the case of the Solar System. The remaining two-thirds are instead ejected into interstellar space, demonstrating that systems such as HD 163296 are efficient factories of interstellar objects.

Building on the above findings, we can draw the following global implications for the evolution of protoplanetary discs hosting young massive planets:

1. Depending on the specific architecture of the system, the dynamical and collisional transport process can influence the composition of already formed planets, including those in the habitable zone under the right conditions. In the case of solid planets, this can have a major impact on their budget of volatiles and astrobiological material. In the case of giant planets, the same processes can alter or overwrite their original atmospheric compositional signatures.

2. The interplay between dynamical transport, thermal ablation and collisions can revert significant fraction of the injected and implanted planetesimals back into dust and pebbles. This can prolong the temporal window available to form planets by pebble accretion, bypassing the barrier effect of existing massive planets to the inwards drift of the pebbles from the outer disc.

3. The dependence of the collisional enrichment on the disc mass can, in principle, put further constraints on the disc mass. Specifically, signatures of planetesimal impacts, such as the presence of dust or refractory elements in the atmospheres of giant planets, would point to high disc masses, while their absence would favour lower disc masses.

Overall the picture depicted by our results indicates that studies aimed at constraining the formation history of planets through their atmospheric characterisation need to account for the global architecture of the host planetary system. This is important because the presence of massive planets could imply that the present-day atmospheric signatures are not a direct reflection of the primordial ones.

Acknowledgements

We thank the reviewer for their comments that helped give more clarity to the manuscript. This work is supported by the Fondazione ICSC, Spoke 3 “Astrophysics and Cosmos Observations”, National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19)” - Next Generation EU (NGEU). D.P. and P.S. acknowledge the support from the Istituto Nazionale di Oceanografia e Geofisica Sperimentale (OGS) and CINECA through the program “HPC-TRES (High Performance Computing Training and Research for Earth Sciences)” award numbers 2022-05 and 2022-02. The authors acknowledge the support from the ASI-INAF grant no. 2021-5-HH.0 plus addenda no. 2021-5-HH.1-2022 and 2021-5-HH.2-2024 and grant no. 2016-23-H.0 plus addendum no. 2016-23-H.2-2021 and from the PRIN INAF 2019 PLATEA and HOT-ATMOS and the INAF Main Stream project, CUP: C54I19000700005. This work was partly supported by the Italian Ministero dell’Istruzione, Università e Ricerca through the grant Progetti Premiali 2012-iALMA (CUP C52I13000140001). This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska–Curie grant agreement No. 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). JMT-R acknowledges financial support from the project PID2021-128062NB-I00 funded by MCIN/AEI/10.13039/501100011033. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made extensive use of NASA’s Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC21M00561. The authors would like to acknowledge the computational support from John Scige Liu and the Genesis cluster at INAF-IAPS.

References

All Tables

All Figures

Fig. 1 Dynamical excitation of the planetesimal disc at 0.75, 2 and 5 Myr in the four simulated scenarios. Dynamical excitation accounts for the contributions of both eccentricity and inclination and is defined as $\sqrt{e^2+{\sin }^2(i)}$ (Petit et al. 2001). The colour bar indicates the region in the disc from which the planetesimals originated (in units of au). The dashed grey line at 40 au indicate the division between inner and outer disc. The rectangles in the rightmost panels mark the “fountain” feature we refer to in the main text, i.e. the damping of the orbital eccentricities of the excited planetesimals by aerodynamic drag that creates the implanted planetesimals. The light blue filled circles indicate the location of the four simulated planets in the disc. First row: high disc mass – high planetary mass. Second row: high disc mass – low planetary mass. Third row: low disc – high planetary mass. Fourth row: low disc mass – low planetary mass.

In the text

Fig. 2 Pericentre vs. apocentre plot of the planetesimal disc at the end of our simulations at 5 Myr, where we can clearly distinguish the excited and unperturbed planetesimal populations. Excited planetesimals from the outer disc can penetrate the inner disc down to a few au from the host star. Planetesimals are colour-coded based on their formation region, the colour scale is truncated at 100 au to better resolve the inner disc. The vertical and horizontal dashed lines at 40 au split the disc in the inner and outer region.

In the text

	Fig. 3 Same as Fig. 2 with the excited planetesimals being colour–coded according to the compositional region they originate from. This compositional gradient is based on HD 163296’s temperature profile and on the condensation temperatures as ices of the main volatile molecules.
In the text

Fig. 4 Relative efficiency between the HPM and LPM scenarios (solid curves) and between the HDM and LDM scenarios (dashed curves) in implanting planetesimals in the inner disc. The orange colour identifies the high-mass disc when comparing the efficiencies of the HPM and LPM scenarios, and the high-mass planets when comparing the efficiencies of the HDM and LDM scenarios. The blue curves identify the LDM and LPM scenarios, respectively. High relative efficiencies in the curves can still be associated with limited absolute efficiencies (given in Table 2). Increasing the disc mass enhances the implantation efficiency in the part of the inner disc closer to the giant planets, while increasing the planetary masses enhances the implantation efficiency closer to the host star.

In the text

Fig. 5 Enrichment of possible undiscovered planets in the inner disc by the impacts of planetesimals that formed in the outer disc for all four simulation scenarios (with blue we denote the HDM-HPM scenario, with green the HDM-LPM one, with orange the LDM-HPM one and with yellow the LDM-LPM one). Top row: normalised mass fraction of accreted material onto lunar-like, Earth-like and super-Earth targets, with the horizontal dashed line marking the water mass fraction on Earth (5×10−4, Morbidelli et al. 2000) shown from left to right. Bottom row: atmospheric enrichment of mini-Neptune, Neptune-like an Jupiter-like targets assuming a solar-like metallicity atmosphere, with the horizontal dashed line marking Jupiter’s 3× solar atmospheric metallicity (Atreya et al. 2018), shown from left to right.

In the text

	Fig. 6 Planetesimal ejection efficiency across the outer disc in the four disc-planets scenarios (from left to right: HDM-HPM, HDM-LPM, LDL-HPM and LDL-LPM). The ejection efficiency is expressed as the fraction of planetesimals originally present in each 10 au-wide ring that end up being ejected from HD 163296’s disc.
In the text