© The Authors 2023

1 Introduction

High-resolution observations of dust in protoplanetary discs by the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed complex features in the form of gaps and rings, arcs, and spirals, commonly referred to as protoplanetary disc substructures (Bae et al. 2022). Since then, many studies have focussed on these disc substructures and their relation to disc evolution and planet formation. Theoretical explanations that have been proposed to explain these features include, among others, pebble growth near condensation zones (Zhang et al. 2015), zonal flows (Flock et al. 2015), self-induced dust traps (Gonzalez et al. 2017), gravitational instabilities (Takahashi & Inutsuka 2016), or the interaction between the disc and the embedded growing planets (e.g. Pinilla et al. 2012; Long et al. 2018).

In some cases, there is additional evidence, derived from gas kinematics, for the presence of growing planets in the observed gaps of protoplanetary discs (Müller et al. 2018; Keppler et al. 2018; Benisty et al. 2021). By measuring the rotation curves of CO emission, Teague et al. (2018) showed that the perturbations of the gas surface density in the HD 163296 disc are likely caused by the presence of two Jupiter-mass planets. Similar work by Pinte et al. (2018) also found that a two Jupiter mass planet in the HD 163296 disc located at 260 AU can reproduce the observed CO emission curve.

For these reasons, many authors have investigated the possibility that the observed substructures in discs are generally caused by growing planets, which would allow indirect information to be derived on the early stages of the planet-formation process. We note that the DSHARP project has detected disc substructures across the range of 5–150 AU (Huang et al. 2018), while the inner regions of protoplanetary discs, that is at distances where fully formed exoplanets are more readily detected, remain difficult to resolve (Andrews et al. 2016). Indeed, there are strong observational biases to keep in mind, with limitations in terms of orbital separations as well as disc and stellar masses (Bae et al. 2022). Still, if the observed rings and gaps are indeed caused by the presence of planets, their existence implies that they should continue to grow and migrate, potentially into regions of the parameter space that are accessible via radial velocity (RV) surveys. These putative planets would then constitute the progenitors of the RV exoplanet sample observed today.

Several studies have been performed in order to investigate the origin of gaps and rings in protoplanetary discs by examining planet-disc and planet-planet interactions and relate current models with exoplanet observations. Lodato et al. (2019) studied 48 gaps in the context of planet-disc interactions. Assuming the same fixed disc properties and disc lifetime for all the investigated gaps, they evolved and migrated the planets for 3–5 Myr. After 3–5 Myr, they found that the final distribution of planets in a mass-semi-major axis plot is consistent with the observed exoplanet distribution, in particular for cold Jupiters. A similar study, again assuming one planet per disc, was carried out by Ndugu et al. (2019). They studied the evolution of planets above pebble isolation mass for a few million years and they compared their resulting population with observations. In contrast to Lodato et al. (2019), Ndugu et al. (2019) find that their planet population is not in agreement with giant planet observations. The main differences arise from the different viscosity values used to model the disc: the low viscosity of the study by Ndugu et al. (2019) prevents the planets from migrating into orbital separations where they could be observed via RV, while the higher viscosity in the study by Lodato et al. (2019) results in efficient inward migration into locations where the giants are observed today.

The first work also including the planet-planet interaction was done by Müller-Horn et al. (2022). Instead of covering a range of observed disc substructures in a statistical sense but treated individually, they focussed on the HD 163296 disc, where they assumed the presence of three interacting planets in observed gaps at 48, 83, and 137 AU; they evolved their systems for a total of 100 Myr, including a few million years lifetime of the gas disc. Their findings regarding the eccentricity distribution showed that in that system the interaction between planets did not result in large eccentricities to match the observed giant planet distribution. They also note that higher viscosity values are more favourable in order to reproduce a planet population that is consistent with current observations, in line with the single body simulations (Ndugu et al. 2019; Lodato et al. 2019).

Our work aims to also include planet-planet interactions while keeping a more general approach without focussing on a specific target. Motivated by the gaps and rings of the observed protoplanetary discs, we run N-body simulations with three, four, five, and seven planets. The initial locations and masses of the planets are informed by the positions of the gaps in the observed discs and the local pebble isolation mass (the mass at which a planet is able to create a sufficiently strong pressure gradient to trap dust particles outside its orbit) and we study how these systems would continue to grow, migrate and gravitation-ally interact for a total of 100 Myr (including a few million years of gas disc evolution). Our goal is to then compare our resulting systems to the observed exoplanet distributions to see if these forming planets could indeed represent the progenitors of the observed exoplanets.

The paper is organised as follows. In Sect. 2, we describe our N -body code FLINSTONE and discuss the different parameters we used for our disc model as well as the initial conditions of our simulations (Sect 2.4). In Sect. 3, we present the results of our simulated systems with three, four, five, and seven planets in Sects. 3.1–3.3 respectively. In Sect. 4, we discuss our findings for the different planetary systems and we finally summarise in Sect. 5 the main outcomes.

2 Method

For our simulations we used the N-body code FLINTSTONE, which includes pebble and gas accretion, gas-driven planetary migration and eccentricity/inclination damping and disc evolution (Izidoro et al. 2017; Bitsch et al. 2019). The integrations were performed by the original N-body integrator that is based on the Mercury hybrid symplectic integrator (Chambers 1999). Collisions between planets were considered as inelastic merging events. In the following we describe the different parts of the model and our setup.

2.1 Disc model

The gas surface density profile Σ_g of the disc is radius and time dependent and can be described as $$\begin{array}{l}{\text{Σ}}_{\text{g}}\left(r,t\right)={\text{Σ}}_0\cdot {\left(\frac{r}{r_0}\right)}^{-1}\exp \left[-\left(\frac{r}{r_c}\right)\right]\exp \left[-\left(\frac{t}{{\tau }_{\text{visc}}}\right)\right]\hfill \\ \{\begin{array}{ll}1,\hfill & t<t_{\text{diss,ph.ev.}}\hfill \\ \exp \left[-\left(\frac{t-t_{\text{diss,ph.ev.}}}{{\tau }_{\text{diss,ph.ev.}}}\right)\right]\hfill & t>t_{\text{diss,ph.ev.}}\hfill \end{array}\hfill \end{array}$$(1)

Here Σ₀ is the gas surface density at r = 1 AU. The surface density profile starts to decrease exponentially with radius at the characteristic value of r_c = 200 AU, giving a total disc mass of 0.074 M_⊙. Moreover, the surface density also decreases exponentially on a timescale of τ_visc = 3 Myr. To take the dissipation process during the final stages of the disc lifetime into account we assumed that the surface density profile decays exponentially with a timescale of τ_diss,ph.ev. = 35 Kyr from t_diss,ph.ev. = t_disc − 100 Kyr until the end of the lifetime of the disc. We summarized all the disc parameters in Table 1.

For the temperature profile of the disc, under the assumption of a locally isothermal state, we chose power law profile $$T\left(r\right)=T_0\cdot {\left(\frac{r}{r_0}\right)}^{-{\beta }_{\text{T}}},$$(2)

where T₀ = 150 K and ß_T is a power law index. We chose ß_T = 3/7 to achieve a flaring disc structure in the outer regions, as expected in discs dominated by stellar irradiation (Chiang & Goldreich 1997; Dullemond & Dominik 2004; Bitsch et al. 2013, 2015a): $$h\left(r\right)=\frac{H}{r}\propto {\left(\frac{r}{r_0}\right)}^{2/7}.$$(3)

Using the above parameters, we obtain a disc that is Toomre stable with Q parameter $Q=\frac{c_s\kappa }{\pi G{\text{Σ}}_{\text{g}}}>$ 10 for all the radii.

2.2 Gas accretion

We are interested in studying the evolution of planets that can open gaps in the pebble distribution in the disc, as we observe with ALMA. These planets need to be at least of the order of the pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018), the mass at which the planet opens a partial gap in the disc resulting in a pressure bump exterior to its orbit, where the pebbles accumulate. We use here the pebble isolation mass prescription from Bitsch et al. (2018): $$M_{\text{iso}}=25f_{\text{fit}}{M}_{\text{Earth}}$$(4)

with $$f_{\text{fit}}={\left(\frac{h}{0.05}\right)}^3\left[0.34{\left(\frac{\log \left(0.001\right)}{\log \alpha }\right)}^4+0.66\right]\left(1-\frac{\frac{\partial lnP}{\partial lnr}+2.5}{6}\right)$$(5)

where α is the viscosity parameter. We used the α-disc model of Shakura & Sunyaev (1973) for which the α parameter is related to the turbulent viscosity as $\upsilon =\alpha H^2{\text{Ω}}_{\text{K}}.\frac{\partial \text{In\hspace{0.17em}}P}{\partial \text{In\hspace{0.17em}}r}$ is the radial pressure gradient.

Once planets have reached the pebble isolation mass they cannot accrete pebbles anymore because the pressure bump that the planets create outside of their orbit stops further pebbles from drifting inwards; thus, they can start to accrete gas into their envelope, because the heating from infalling pebbles stops. The first stage of gas accretion is regulated by the thermodynamics of the envelope, which must cool and contract in order for the planet to accrete more gas. For the contraction rate of the planetary envelope we use the same formula as in (Bitsch et al. 2015b) $$\begin{array}{l}{\dot{M}}_{\text{gas}}=0.000175f^{-2}{\left(\frac{{\kappa }_{\text{env}}}{1{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}{\left(\frac{{\rho }_{\text{c}}}{5.5{\text{cm}}^{-3}\text{g}}\right)}^{-1/6}{\left(\frac{M_{\text{c}}}{M_{\text{E}}}\right)}^{11/3}\hfill \\ \times {\left(\frac{M_{\text{env}}}{M_{\text{E}}}\right)}^{-1}{\left(\frac{T}{81\text{K}}\right)}^{-0.5}\frac{M_{\text{E}}}{\text{Myr}}\hfill \end{array}$$(6)

where f is a factor to change the accretion rate and is set to f = 0.2 as in Piso & Youdin (2014) while for the opacity in the planetary envelope κ_env we used the same value as in Bitsch et al. (2019) κ_env = 0.05 cm² g⁻¹. The core density is ρ_c = 5.5 g cm⁻³. If M_core ≃ M_env the gas accretion phase regulated by cooling and contracting in hydrostatic equilibrium ends and rapid runaway gas accretion starts (M_core < M_env). For rapid gas accretion we follow Machida et al. (2010), who calculated the gas accretion rate onto a planetary core using 3D hydrodynamical simulations. They found that the effective gas accretion rate is given by Ṁ_gas = min {Ṁ_{gas, low}, Ṁ_{gas, high}} where the two different gas accretion limits are $$\begin{array}{l}{\dot{M}}_{\text{gas,low}}=0.83{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{g}}{H}^2{\left(\frac{r_H}{H}\right)}^{\frac{9}{2}}\hfill \\ {\dot{M}}_{\text{gas,high}}=0.14{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{g}}{H}^2\hfill \end{array}$$(7)

where r_H is the Hill radius of the planet r_H = r · (m_pl/3 M_⋆)^1/3, Ω_K is the Keplerian frequency and Σ_g is the surface density. Because the gas accretion rate of a planet depends on what the disc can provide, in our implementation of gas accretion we follow D’Angelo & Lubow (2008), where they found that the gas accretion is limited to 80% of the accretion rate of the disc, which is $${\dot{M}}_{\text{disc}}=3\pi \alpha h^2{r}^2{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{g}}.$$(8)

2.3 Orbital migration

Planets interact with their surrounding disc material by exchanging energy and angular momentum. The exchange of angular momentum can be explained by the wave torques in the disc (Goldreich & Tremaine 1980). Low-mass planets (planets around a few Earth masses), migrate in the type-I regime (Paardekooper et al. 2011) where the total torque is given as $${\text{Γ=Γ}}_{\text{L}}{\text{Δ}}_{\text{L}}+{\text{Γ}}_{\text{C}}{\text{Δ}}_{\text{C}}$$(9)

where Г_L is the Lindblad torque and Г_C is the corotation torque while Δ_L and Δ_C are the rescaling functions of the Lindbland torque and corotation torque respectively (Izidoro et al. 2017) to account the reduction of the total torque according to the orbital eccentricity and inclination of the planet. The nominal torques for zero eccentricity are calculated according to Paardekooper et al. (2011). This migration prescription is consistent with 3D simulations (Bitsch & Kley 2011). The same setup was used in Bitsch & Kley (2010), who studied the eccentricity evolution of low mass planets and find that the damping rates in this case are consistent with the isothermal case.

When the planets are more massive, they open deep gaps in the protoplanetary discs, and they start to migrate in the type-II regime. This happens when (Crida et al. 2006) $$𝒫=\frac{3}{4}\frac{H}{r_{\text{H}}}+\frac{50}{qℛ}\le 1.$$(10)

Here q is the star to planet mass ratio m_pl/M_⋆ and ℛ is the Reynolds number r²Ω_K/υ. For a smooth transition from type-I to type-II migration we follow Kanagawa et al. (2018), who relate the type-II migration timescale with the type-I as $${\tau }_{\text{mig},\text{II}}=\frac{{\text{Σ}}_{\text{un},\text{p}}}{{\text{Σ}}_{\min }}{\tau }_{\text{mig},\text{I}}$$(11)

where Σ_unp is the surface density of the unperturbed disc while Σ_min corresponds to the surface density at the bottom of the gap. They are related as (Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2015) $$\frac{{\text{Σ}}_{\text{un},\text{p}}}{{\text{Σ}}_{\min }}=1+0.04\text{\hspace{0.17em}K}$$(12)

with $$K={\left(\frac{M_{\text{pl}}}{M_{\star }}\right)}^2{\left(\frac{H}{r}\right)}^{-5}{\alpha }^{-1}.$$(13)

To account for the damping of the orbital eccentricity and inclination due to the gravitational interaction between the disc and the planets, we follow Cresswell & Nelson (2006, 2008) who found that the timescale of the eccentricity damping is best described by $$\begin{array}{l}{\tau }_{\text{e}}=\frac{{\tau }_{\text{wave}}}{0.780}[1-0.14{\left(\frac{e}{h}\right)}^2+0.66{\left(\frac{e}{h}\right)}^3\hfill \\ +0.18\left(\frac{e}{h}\right){\left(\frac{i}{h}\right)}^2]\hfill \end{array}$$(14)

while the timescale for the inclination damping is given by $$\begin{array}{l}{\tau }_{\text{i}}=\frac{{\tau }_{\text{wave}}}{0.544}[1-0.30{\left(\frac{i}{h}\right)}^2+0.24{\left(\frac{i}{h}\right)}^3\hfill \\ +0.14{\left(\frac{e}{h}\right)}^2\left(\frac{i}{h}\right)]\hfill \end{array}$$(15)

where τ_wave is the characteristic time of the orbital evolution (Tanaka & Ward 2004) $${\tau }_{\text{wave}}=\left(\frac{M_{\star }}{M_{\text{pl}}}\right)\left(\frac{M_{\star }}{{\text{Σ}}_{\text{g},\text{pl}}{\alpha }_{\text{pl}}^2}\right)h^4{\text{Ω}}_{\text{K},\text{pl}}^{-1}.$$(16)

These timescales for eccentricity and inclination damping are valid for low mass planets that undergo type-I migration. For planets that open gaps and start to migrate in the type-II regime, we follow the classical K-damping prescription (Lee & Peale 2002). The damping rates of eccentricity and inclination for massive planets that migrate in type-II regime are $$\begin{array}{ll}\dot{e}/e=-K_{\text{damp}.}\left|\dot{a}/a\right|,\hfill & \dot{i}/i=-K_{\text{damp}.}\left|\dot{a}/a\right|\hfill \end{array}.$$(17)

Here we use K_damp = 5. Higher K_damp values result in higher damping rates. For a better understanding of how different K_damp values affect the formation of giant planets, we refer to Bitsch et al. (2020) where they investigated the influence of different damping rates on the eccentricity and inclination of growing and migrating giant planets.

Fig. 1 Power law fit between the position of the gaps observed by ALMA 1.25 mm continuum images (DSHARP; Huang et al. 2018) and the distance between them.

2.4 Initial conditions

We modelled the evolution in the case where we have three, four, five, and seven putative planets. For the position of these planets in each case, we took the annular substructures of the 18 disc systems from the Disc Substructures at High Angular Resolution Project (DSHARP; Huang et al. 2018) and fit a power law between the position of these substructures and the distance between them. The result of the fit is shown in Fig. 1. Considering an inner planet located at 10 AU, this yields initial reference locations at 10, 25, 47, 72, 97, 118, 135, 147, and 156 AU. To account for the uncertainties in the initial position of the planets, we ran simulations for two different configurations, one where the planets are located in a wide configuration and another where they are in a compact configuration. The wide configuration for the three planet case harbors initial positions at 10, 47, and 97 AU while the compact configuration has initial planetary positions at 25, 47, and 72 AU. For the four planet case the wide configuration gives starting positions at 10, 47, 97, and 135 AU and for the compact configuration at 25, 47, 72, and 97 AU. To further address any uncertainties, we also add an extra 5% of deviation for the individual starting positions of our planets. For this planet distribution, we assume that the planets have just reached the pebble isolation mass (Eq. (4)) and they start to accrete gas. We chose to test two distinct simulation setups for their initial mass. The first case is where we considered that the planets have already accreted 10% of the gas on top of the pebble isolation mass. As second case we have a mass distribution chosen randomly between 1 and 5 times of the pebble isolation mass plus 10% of gas. Finally, the disc lifetime t_disc were randomly selected from the range of [3–10] Myr (Table 1) although we evolve the systems until 100 Myr to study also their evolution after the dissipation of the disc.

As in Bitsch et al. (2019) the initial inclinations and eccentricities of the planets were distributed uniformly with inclinations between 0.01°–0.05° and eccentricities between 0.001 – 0.01. The other orbital angles were randomly selected between 0°–360°. We used three different α-viscosities, α = 10⁻³, 3.16×10⁻⁴, and 10⁻⁴. Higher values of viscosity result in higher migration and accretion rates. Higher α values result in a shallower gap because of Eq. (13), and the higher accretion rate originates from Eq. (8). For each α value, we ran 50 simulations with different initial conditions for the positions of the planets. Figure 2 shows a tree diagram illustrating all the different sets of our simulations.

3 Results

In this section, we present the outcomes of our simulations regarding the distribution of mass, semi-major axis, and eccentricity for different planetary systems. We organise our results by considering the number of planets in each system. We first show results for the three-planet system in Sect. 3.1 and next we present our results for the four-planet system in Sect. 3.2. We then investigate systems consisting of five and seven planets in Sect. 3.3.

3.1 Three planet system

Figure 3 shows the final locations and masses of the simulated three-planet systems on a mass vs. semi-major axis plot after 100 Myr of integration time (left) as well as their eccentricity distribution (right). This plot includes both configurations for the position of the planets and also both configurations for their initial mass. Orange dots refer to the initial conditions where the planets are located in a wide configuration whereas the brown dots refer to the more compact configuration. The different colours represent simulations of different α viscosity values. The horizontal line segments give an indication of how eccentric our simulated planets are as they extend from the perihelion a(1−e) to aphelion a(1 + e). To provide a clearer indication of the final eccentricity, we have included the average eccentricity values for the three sets of different α values. These values are represented by coloured numbers.

In general the planets experience gas accretion and migration towards the central star, as evident by the increased mass and inward migration. As we mentioned earlier, we expect to have higher migration and accretion rates for the higher values of the α-viscosity 10⁻³. This is exactly what we see in the left panel of Fig. 3. In Table 2, we show the outcomes of our simulations for the different α viscosity values. For higher values of the α viscosity the accretion of gas and the orbital migration is more efficient. On the other hand, simulations with the lowest value α = 10⁻⁴ exhibit relatively slower rates of migration and accretion, while for the intermediate α value, we also find intermediate planetary masses. These results are consistent with those of Müller-Horn et al. (2022), where they investigated a specific target, the disc of HD 163296, by modelling the evolution of the three planet candidates.

To compare our results with observations, in Fig. 3 we also plotted the current observed exoplanets (grey dots). The data we used were obtained from the NASA Exoplanet Archive and were taken on 3 of February in 2023¹. The vertical black dashed line represents the current RV detection limit at 5 AU. This limit corresponds to the time limit of how long observations would need to be done, so it is a mass-independent limit which is valid for the mass range considered here (Wright & Gaudi 2013). Only a small fraction of the synthesised planets migrate into the inner disc and none interior to 1 AU. This already hints to the fact that planets such as those assumed to be originating from gaps in protoplanetary discs might not be the sole source of the currently observed giant planet population. In addition, there is a considerable number of planets that lie beyond the RV detection limit of 5 AU. When these detection limitations are removed, we will be able to confirm or exclude the presence of such wide-orbit giants.

In the right panel of Fig. 3, we present the cumulative distribution of the eccentricities, separately for each α-viscosity value, as well as for the combined dataset (blue line). We also show the cumulative distribution of eccentricity for the observed exoplanet sample as a dashed black line. The dataset is limited to include only those simulated planets with semi-major axis of 1 ≤ a [AU] ≤ 5 and masses within the range of 0.5 ≤ m [M_Jup] ≤ 13. These values are marked by the black box in the left panel of Fig. 3. We also limit the observed exoplanet sample to include only systems for which the mass of the central star is between 0.7 ≤ m [M_⊙] ≤ 1.3 in order to have a fair comparison with our simulated systems where we have assumed the central star mass of m = 1 M_⊙. It is clear that the eccentricity distribution of our simulated planets is too low (last column of Table 2) compared to the eccentricity distribution of the observed giant planets (mean value of 0.302).

Fig. 2 All the different setups we used for our simulations. These setups were tested for three α values (see Table 1).

Fig. 3 Left panel: comparison between the observed giant exoplanet population (grey dots) and all the simulations for the three-planet systems. The initial values of these planets are denoted as orange for the case where the planets located in a wide configuration while brown dots denote planets located in a compact configuration, where the small numbers indicate the position of the first, second, and third planet. Purple, green, and red data points symbolise the simulations for the different α, 10−3, 3.16 × 10−4, and 10−4, respectively, whereas the coloured numbers shows the mean value and the standard deviation of the eccentricity for each case for all the simulated three- planet systems. The horizontal lines refers to the perihelion and aphelion of the planet a(1 ± e). The vertical black dashed line represents the current RV detection limit at 5 AU. Right panel: cumulative frequency of the eccentricity for planets up to 5 AU and mass between 0.5 ≤ mass(MJup) ≤ 13. We also note that these plots include both cases for the initial mass (fixed and random value).

Fig. 4 Similar to Fig. 3, but for systems with initially four planets.

3.2 Four planet system

Figure 4 shows the mass vs. semi major axis plot of the simulated four-planet systems (left) and their eccentricity distribution (right) after 100 Myr of evolution. The plots are analogous to Fig. 3. We also note that for the initial mass in these systems we only tested the case where m = M_iso + 10% gas. If we increase the mass any further, the systems would become unstable on a very short time scale, see Sect. 4.2 for further discussion on this issue.

In the left panel of Fig. 4, we see that the planets experience a similar evolutionary pattern to the three-planet system. The planets undergo gas accretion and migration towards the central star, and the efficiency of these processes increases as the α viscosity increases. We present the final mass and eccentricity in Table 2 for the different α viscosity values as well as for the combined dataset.

The right panel of Fig. 4, shows the cumulative distribution of eccentricities for our simulated four-planet systems after 100 Myr of integration time. We compared to the same set of planets regarding their mass, distance and stellar host mass as for the three-planet case. For the lowest α viscosity value, α = 104⁻⁴, the eccentricities are too low, whereas for higher viscosities for which the rates of gas accretion and migration are more effective, the final eccentricities are a little bit larger. However, our simulated eccentricities are clearly too low to match the observations. The precise values of the eccentricities for the limited dataset are presented in the last column of Table 2.

Compared to the three-planet system, the distribution of mass and semi-major axis of our simulated planets are very similar to each other. On the other hand, the eccentricity distribution for systems consist of four planets is slightly higher – two times higher – than the systems with three planets, but still too low to match the observations. During the gas disc phase, the eccentricity of the planets are damped because of their interaction with the disc, resulting in low eccentricities. After the dissipation of the disc, the planets can interact with each other leading to an increase in their final eccentricity. However, the distances between the objects are too large to allow efficient scattering events. Table 3 shows the average number of surviving planets at the end of our simulations. We see that for systems with initially four planets the occurrence of scattering events between planets is slightly higher. Hence, it is reasonable to obtain higher mean eccentricity for the four-planet systems.

3.3 Five and seven planet system

For the five planet system we studied two configurations for the initial position of the planets whereas for the seven planet system we studied only systems in compact configuration. For both systems we also tested two configurations for the initial mass. Figure 5 shows two individual systems of five (top) and seven planets (bottom). Each planet is represented by different colour. The right panel shows the evolution of the eccentricity over time and the left panel shows the evolution of the semi-major axis. We integrate each system for 40 Kyr. From Fig. 5 we see that both systems are initially unstable leading to scattering between the outermost planets. We also note that the individual systems presented in Fig. 5 are not the exception but rather a common outcome of these simulations. All the systems of five and seven planets and different simulation setups we studied exhibit this instability early on.

We can explain these initial instabilities to some level, by examining the stability of planetary systems via the Hill criterion. Gladman (1993) studied the dynamics of two small planets orbiting the Sun. He found that in order to have stable systems, the difference of the initial semi-major axis between two planets measured in mutual Hill radius R_H = + m₂)/3 M_⊙]^1/3[(a₁ + a₂)/2], must be greater than $2\sqrt{3}$. In Table 4, we show the distance between two planets in mutual Hill radius at t = 0 yr. These distances were calculated by assuming that the mass of the planets is m = M_iso + 10%. We see that the distance between the outermost planets in both systems is much lower than the critical value of $2\sqrt{3}$, explaining why the systems become unstable. This instability is clearly not prevented by the damping forces of the disc.

Fig. 5 Evolution of two different individual systems for the 5 (top) and 7 (bottom) planet system. The plots shows the evolution of the semi-major axis (left) and the eccentricity (right) over time.

4 Discussion

The presence of rings and gaps in the dust emissions of pro-toplanetary discs is often attributed to young planets and more decisive evidence of the presence of planets has been additionally mustered via gas kinematics (Pinilla et al. 2012; Flock et al. 2015; Zhang et al. 2015; Takahashi & Inutsuka 2016; Gonzalez et al. 2017; Long et al. 2018). Still, it is not clear at the present day whether these substructures should be attributed to planets as a general rule, nor whether each dip in the dust emission is to be associated with one single planet (Dong et al. 2017; Bae et al. 2017). Another complementary but related question concerns the fate of such growing planets in the context of planet formation theory and the currently observed exoplanet sample. In particular, we can ask whether these putative young forming planets represent the progenitors of the currently observed giants in the sample: if yes, we can draw information on the formation history of the currently observed planets and better inform planet formation theory; if not, we should either expect a yet-undetected population of planets, or, if future observations rule out such an unseen population, put into question our formation scenarios of planets forming efficiently in gaps. To this end, in this section we systematically compare the outcome of our simulations of growing and migrating multi-planetary systems with the observed giant exoplanets.

4.1 Stable systems

We focus here on the outcomes concerning the systems with initially three and four planets. Compared to the observed eccentricity distribution our simulations clearly show too low eccentricities of the surviving planets. A reason for this outcome could be that the planets do not interact much with each other, avoiding scattering events resulting in low eccentricities. The average number of the surviving planets after 100 Myr of integration time is 2.83 for the three planet system and 3.50 for the four planet system.

The further evolution beyond 100 Myr would most likely not affect the eccentricity and semi-major axis distribution of our systems, because the distances in units of mutual Hill radii after 100 Myr of evolution are very large indicating that scattering events should be rare.

Raymond et al. (2009) studied the eccentricity evolution of a three-planet system located initially in an unstable configuration with a separation of 4–5 R_H. They found that scattering among planets in this unstable configuration result to an eccentricity distribution that agrees with observations. Table 5 shows the distances of our simulated planets in mutual Hill radius for the different α viscosities after 100 Myr. These values correspond to both configurations for the initial position, compact and wide configuration, of the planets and for both configurations regarding the initial mass for the three planet system and for m = M_iso + 10% for the four planet system. In our simulations, the surviving planets have a larger separation than in Raymond et al. (2009), indicating why our simulations do not produce efficient scattering events and thus reflect a population of low eccentricity planets.

A way of increasing the excitation of planetary systems is the presence of leftover planetesimals. At later stages of the disc lifetime where the gas surface density drops, the dust-to-gas ratio goes up, which is a favourable condition for planetesimal formation (Carrera et al. 2015); moreover planetesimal-disc-driven instabilities are believed to have sculpted the orbital distribution of the giant planets in the outer solar system, increasing the eccentricities by a factor of 10 (Tsiganis et al. 2005; Nesvorny 2018). Another possibility is the presence of planets in the inner disc. Bitsch et al. (2020) and Bitsch & Izidoro (2023) studied the eccentricity distribution of protoplanetary embryos located in the inner disc by running N-body simulations. They found that the simulations that match best the observed eccentricity distribution of giants are those where the planets interact with each other. However, this does not necessarily rule out the possibility that planets can also exist in the outer disc. Injecting the planets according to the gaps’ positions in the currently available observations, we limited our study only to planets located in the outer disc, due to observational biases. Indeed, the inner part of the disc is not resolved with the current sensitivities and so we do not observe many gaps and rings in that region. One exception is the gap located at 1 AU in the TW Hya disc (Andrews et al. 2016).

4.2 Unstable systems

For all the different initial conditions we tested, we found that systems with five and seven planets become unstable immediately, resulting in scattering events and increased eccentricities of the surviving planets, see Fig. 5.

Lega et al. (2013) studied the planet-planet and the planet-disc interaction in a three planet system with hydrodynamical simulations. They found that when the planets reach Jupiter mass the systems become unstable. This instability disturbs the gas surface density resulting in complex patterns, where gaps and rings are not visible. However, the dust emissions of protoplanetary discs observed in the DSHARP project do not show complex gas density distribution but rather more well defined gaps and rings. This indicates that not all gaps and rings in observed protoplanetary discs are caused by planets.

In particular, we stress that the initial configuration for our seven planet system is very similar to the position of the seven gaps in the AS 209 system. This strongly indicates that for such discs it is not possible for all the observed gaps and rings to be caused by planets above the pebble isolation mass located at the positions of the gaps. On the other hand, planets below the pebble isolation mass cannot be responsible for the observed gap structures. Therefore, additional processes that can cause gaps and rings without necessarily assuming one forming planet for each gap are needed to understand disc substructure observations. For example, Bae et al. (2017) showed that a single 30 Earth-mass planet in a low (α = 5 × 10⁻⁵) viscosity disc can reproduce multi-ring features. Additionally, Dong et al. (2018) also found that the multiple gaps in the HL Tau, TW Hya and HD 163296 can be produced by a single planet in a low viscosity disc (α ≤ 10⁻⁴). Both studies align with our findings that not each gap is explained by the presence of a planet. Moreover, there exist various purely hydrodynamical, magneto-hydrodynamical and dust-evolution mechanisms that alone can produce substructures such as rings and gaps in protoplanetary discs (Rice et al. 2004; Marcus et al. 2015; Flock et al. 2015, 2017; Suriano et al. 2018; Bae et al. 2019; Li & Youdin 2021; Kuznetsova et al. 2022).

4.3 Disc model

In this study, we tested the effects of three different values of α viscosity within our simulations (e.g the viscosity influences the disc structure as well as the planetary migration and accretion rates), but we limited our study to a constant value of α with respect to time and orbital distance. If the α viscosity was radially dependent, with α increasing with distance from the star, then planets located further out in the disc would have higher migration rates. Such a radial dependence in α would lead to planets being located in a closer configuration and could possibly result to more planet-planet interactions. However, it seems that the viscosity in the majority of the disc might actually be quite low (Lesur et al. 2022; Delage et al. 2022).

We also note that the outcomes of our simulations are not affected by the surface density profile. In fact, we tested all the different simulation setups for another surface density profile where Σ₀ = 283 g cm⁻². Even with the lower value of Σ₀ the final distribution of the planets in the mass-semi major axis plot as well as the eccentricity distribution were very similar to the distribution we presented in Sect. 3. Furthermore, the systems with five and seven planets became unstable immediately in this case as well, as we would expect.

4.4 Implications for observations

Another outcome of our simulations, regardless of the level of the α-viscosity value, is the presence of far-out giant planets farther out than 5 AU which is the current RV limit. Future observations will be able to either confirm or rule out the existence of such planets. In the latter scenario, it would suggest that either our formation models are incomplete, or that the putative planets located at the gaps observed in young discs are not the progenitors of the current exoplanet sample.

5 Summary

In recent years, observations of disc substructures have drawn considerable attention due to their connection to disc evolution, planet formation and their mutual influence: how the disc environment influences the growth, migration and orbital properties of the forming planets and how these forming planets affect the disc structure. In particular, it has been speculated that gaps and rings are caused by young forming planets that have reached the pebble-isolation mass. Once they reach the pebble-isolation mass, these young forming planets perturb the disc and open a gap which generates a pressure bump outside of their orbits preventing the further inward flux of large pebbles.

Many works have investigated the fate of these putative planets as they undergo gas-driven migration and accretion according to our current understanding of planet formation, and have analysed the synthetic planet population in comparison with the observed exoplanet sample. Past works have initially only considered the case of single-planets (i.e. no planet-planet dynamical interactions within the system, although multiple gaps are commonly observed, and multiple planet systems are not uncommon; Lodato et al. 2019; Ndugu et al. 2019), or multi-planet evolution for the specific case of the HD163296 system (Müller-Horn et al. 2022; a 1.5 solar mass star).

In this work, we performed a more comprehensive N-body analysis where we simulated systems with three, four, five, and seven planets to study their growth, migration and also the gravitational interaction with the surrounding disc as well as the interactions between them. We located the planets at the positions where the gaps have been observed (Huang et al. 2018) and consider that our planets have just reached the pebble isolation mass. We investigated the effects of the viscosity parameter assuming three different values of 10⁻⁴, 3.16 × 10⁻⁴, and 10⁻³. In order to study the evolution of the planetary systems after the dissipation of the disc we continued to evolve the systems for 100 Myr in total, while modelling the evolution of the disc for a few million years.

We found that our simulations show an eccentricity distribution that is extremely low compared to the one observed for RV giants. We attribute this result to the low interaction between the planets which prevents scattering events and leads to low eccentricities. To be more precise, it is difficult to reproduce the observed eccentricity distribution by assuming three or four planets located at the positions of the gaps.

Higher values of viscosity cause the planets to migrate more and accrete at higher rates resulting in a more massive and close-in population of giants. Although higher values of viscosity lead to higher final eccentricities, even the highest value of α = 10⁻³ did not reproduce eccentricities large enough to match the observations. In addition, we did not find any strong dependence on the surface density profile.

Finally, our results show that multi-ring discs cannot be explained by a one-planet-per-gap assumption. We found that systems consist of five and seven planets become instantly unstable. This indicates that the observed gaps and rings cannot be caused only by the presence of planets above pebble isolation mass. Additional processes (such as the ones mentioned in Sect. 4.2), must be taken into account in order to understand the observed disc substructures. Thus, except in those cases where additional independent observations such as gas kinematics corroborate the presence of young accreting planets, one cannot attribute all disc substructures to planets in a bijective manner.

Future observations will be able to show if a yet-undetected long-period populations of planets actually exists, or if the putative planets located at the gaps in the observed young discs are not the progenitors of the current exoplanet sample.

Acknowledgements

B.B. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. We thank the referee Ruobing Dong for his report that helped to clarify and improve our manuscript.

Appendix A Distribution of mass and semi-major axis

We present here the cumulative distribution of mass and semimajor axis for systems with three (Fig. A.1) and four planets (Fig. A.2). We limited the dataset to include planets with final semi-major axis between [1–5] AU and final mass between [0.5–13] M_Jup, similar to the data set we presented in Sect. 3 for the cumulative distribution of the eccentricity. The left panels show the cumulative distribution of the mass while the right panels show the cumulative distribution of the semi-major axis. We show the distribution for each value of the viscosity separately as well as for the combined dataset. For both systems we see that for higher values of viscosity the planets experience gas accretion and migration more efficiently. For the combined dataset (blue line) the final mean mass and mean semi-major axis is (3.16 ± 1.98)M_Jup and (2.70 ± 1.18)AU for the three-planet system and (3.33 ± 2.31)M_Jup and (2.56 ± 1.25)AU for the four-planet system. The mean mass of the observed giant planets is (3.35 ± 2.85)M_Jup whereas the mean semi-major axis is (2.40 ± 1.07)AU. We see that both systems, with three and four planets, have similar mass and semi-major axis distribution as the observations, when taking into account the combined set. In general though, we must note that our systems cannot match the observations since our simulated systems did not reproduce the observed eccentricity distribution (see Sect. 3 and Sect. 4).

Fig. A.1 Cumulative distribution of mass (left panel) and semi-major axis (right panel) for the three-planet system. The plots include planets for which the semi-major axis is within the range of 1 ≤ a[AU] ≤ 5 and mass between 0.5 ≤ m[MJup] ≤ 13.

Fig. A.2 Similar to Fig. A.1, but for systems with initially four planets.

References

All Tables

All Figures

	Fig. 1 Power law fit between the position of the gaps observed by ALMA 1.25 mm continuum images (DSHARP; Huang et al. 2018) and the distance between them.
In the text

	Fig. 2 All the different setups we used for our simulations. These setups were tested for three α values (see Table 1).
In the text

Fig. 3 Left panel: comparison between the observed giant exoplanet population (grey dots) and all the simulations for the three-planet systems. The initial values of these planets are denoted as orange for the case where the planets located in a wide configuration while brown dots denote planets located in a compact configuration, where the small numbers indicate the position of the first, second, and third planet. Purple, green, and red data points symbolise the simulations for the different α, 10−3, 3.16 × 10−4, and 10−4, respectively, whereas the coloured numbers shows the mean value and the standard deviation of the eccentricity for each case for all the simulated three- planet systems. The horizontal lines refers to the perihelion and aphelion of the planet a(1 ± e). The vertical black dashed line represents the current RV detection limit at 5 AU. Right panel: cumulative frequency of the eccentricity for planets up to 5 AU and mass between 0.5 ≤ mass(MJup) ≤ 13. We also note that these plots include both cases for the initial mass (fixed and random value).

In the text

	Fig. 4 Similar to Fig. 3, but for systems with initially four planets.
In the text

	Fig. 5 Evolution of two different individual systems for the 5 (top) and 7 (bottom) planet system. The plots shows the evolution of the semi-major axis (left) and the eccentricity (right) over time.
In the text

	Fig. A.1 Cumulative distribution of mass (left panel) and semi-major axis (right panel) for the three-planet system. The plots include planets for which the semi-major axis is within the range of 1 ≤ a[AU] ≤ 5 and mass between 0.5 ≤ m[MJup] ≤ 13.
In the text

	Fig. A.2 Similar to Fig. A.1, but for systems with initially four planets.
In the text