© B. Bitsch et al. 2020

1 Introduction

Since the discovery of the first close-in giant planet around a main sequence star (Mayor & Queloz 1995), many more exoplanets have been detected. These detected planets come in the following different categories: (i) terrestrial planets similar in mass to those as in our inner Solar System, (ii) close-in super-Earths and mini-Neptunes with masses of several Earth masses, and (iii) short period Jovian-mass planets (hot Jupiters, r < 0.1 au), or (iv) more distant giant planets (warm (r < 1.0 AU) and cold Jupiters (r > 1.0 au).

Based on the observational evidence, these different types of planets also show different occurrence rates. For instance, close-in super-Earths seem to be most common, with an occurrence rate of 30–50%, in inner systems (Mayor et al. 2011; Fressin et al. 2013; Mulders et al. 2018), while cold Jupiters seem to exist around 10–20% of solar-type stars (Johnson et al. 2010). Recent observations have argued for a correlation between inner super-Earths and cold gas giants, where up to 90% of cold gas giants should have inner super-Earth systems (Bryan et al. 2019; Zhu & Wu 2018). However, this trend is under debate (Barbato et al. 2018). If such a correlation exists, these inner super-Earths should have formed interior to the orbit of the gas giants (Izidoro et al. 2015).

The observed giant planets seem to show significant eccentricity (e.g. Dawson & Murray-Clay 2013; Buchhave et al. 2018), although the planet eccentricity distribution might suffer from observational biases. Recent observations (Anglada-Escudé et al. 2010; Wittenmyer et al. 2013, 2019a; Kürster et al. 2015; Boisvert et al. 2018; Hara et al. 2019) have claimed that two giant planets in circular orbits could mimic the radial velocity (RV) signature of a single eccentric giant planet if not enough RV measurements have been used to determine the true orbital configuration of the system. In addition, simulations have shown that the eccentricities of multi-planet systems can oscillate in time due to secular perturbations, with quite large variations at any given moment in time (e.g. Raymond et al. 2009a; Trifonov et al. 2014; Sotiriadis et al. 2017; Luque et al. 2019).

The observed high eccentricities of these giant planets can be explained by the following two mechanisms: (i) via gravitational planet-disc interactions (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013) or (ii) via planet–planet scattering events after gas disc dispersal (Jurić & Tremaine 2008; Raymond et al. 2009a; Sotiriadis et al. 2017). The first mechanism is related to the fact that very massive planets (above five Jupiter masses) open very deep gaps that deplete the Lindblad resonances responsible for eccentricity damping, resulting in an increase in the orbital eccentricity, which can even happen for single giant planets. The second mechanism is based on gravitational interactions between multiple planets.

Recent studies have found that eccentric giant planets are more frequently observed around metal rich stars (Dawson & Murray-Clay 2013; Buchhave et al. 2018). In the core accretion paradigm for planet formation (Pollack et al. 1996), this can be explained by the larger amount of building blocks available, which enhances the formation of pebbles, planetesimals (Bai & Stone 2010; Yang et al. 2017) and thus planetary cores via the accretion of planetesimals (Pollack et al. 1996) or pebbles (Johansen & Lacerda 2010; Ormel & Klahr 2010; Lambrechts & Johansen 2012). As more planets are formed, mutual gravitational interactions among multiple giant planets are more likely to lead to dynamical instabilities and scattering events that pump planetary eccentricities to observed levels.

Planet–planet scattering simulations dedicated to explain the observed giant planet eccentricity distribution have typically invoked fully formed giant planets from the beginning. Traditionally, these planets are artificially put close to each other to favour a prompt onset of the dynamical instability once the simulation starts (Jurić & Tremaine 2008; Raymond et al. 2009a) or have only type-II migration of already fully formed giant planet (Sotiriadis et al. 2017). These simulations are quite successful in reproducing the giant planet eccentricity distribution. In addition, all these simulations required that at least three planets are needed to trigger scattering events, if no other sources to disrupt the system are present.

In this work, we aim to reproduce the observed giant planet eccentricity distribution by including the initial step of modellingthe growth and migration from planetary embryos (of around Moon mass) to full planetary systems. We follow the planetary evolution model of Izidoro et al. (2019) and Bitsch et al. (2019). This model includes growth by pebble (Lambrechts & Johansen 2014) and gas accretion (Piso & Youdin 2014; Machida et al. 2010), migration in the type-I (Paardekooper et al. 2011)and type-II (Kanagawa et al. 2018) regime as well as an evolving disc model (Bitsch et al. 2015a). We also vary the initial number of planetary embryos that can grow and migrate to study their influence on the final planetary systems. In addition, we investigate the influence of the damping rates in the type-II regime as free parameter, motivated by the fact that damping rates have been derived in hydrodynamical simulations only for single planets, making it unclear what the exact damping rates are for multiple gap opening planets. Finally, we test two different scenarios for the growth of the giant planets. In the first one, giants planets are only allowed to grow up to one Jupiter-mass. In the second set, gas accretion onto gas giants is only limited by the gas disc flow and it only shuts down at the time of gas disc dispersal, allowing the growth of planets of up to several Jupiter masses.

The code used in this work was already introduced in our previous works, where we studied the formation of super-Earth systems (Izidoro et al. 2019) and gas giants (Bitsch et al. 2019). However, in these previous works we did not touch upon the eccentricity distribution of the giant planets and the super-Earth – cold Jupiter relation.

Our work is structured as follows. In Sect. 2, we summarise the numerical methods used in this work. In Sect. 3, we show the evolution of individual systems and in Sect. 4 we investigate the link to the eccentricity distribution of Jupiter mass planets. In Sect. 5, we investigate how faster growth and thus more massive planets influence the outcomes of our simulations. In Sect. 6, we comment on the implications of our simulations in respect to the super-Earth – cold Jupiter relation found in observations. We then discuss further implications of our results and their shortcomings in Sect. 7, while we summarise in Sect. 8.

2 Methods

In this section, we summarise the used methods and state the parameters of the model that are used in this work. The code is described in more detail in our previous work (Izidoro et al. 2019; Bitsch et al. 2019). In addition, we also discuss the initial conditions of our simulations. Our planet formation model features pebble accretion (Lambrechts & Johansen 2014), gas accretion (Piso & Youdin 2014; Machida et al. 2010), type-I and type-II migration (Paardekooper et al. 2011; Kanagawa et al. 2018) as well as a simple accretion disc model with decaying Ṁ (Bitsch et al. 2015a). We follow the planetary disc model outlined in Bitsch et al. (2015a) and used in the N-body planet formation simulations of Izidoro et al. (2019) and Bitsch et al. (2019). The disc lifetime is in total 5 Myr, and the initial disc age is 2 Myr as in Bitsch et al. (2019), meaning that the planets have 3 Myr in a gas disc environment to grow and migrate via interactions with the gas disc. In the last 100kyr we employ an exponential decay of the gas disc surface density to mimic the disc dispersal. We set the inner edge of the gas disc to be at 0.1 AU, in line with recent hydrodynamical models (Flock et al. 2019). Collisions between planetary embryos are treated as perfect mergers between the bodies.

We want to stress here that our simulations are designed to form giant planets, by choosing disc properties and pebble fluxes that allow the efficient growth of giant planets, motivated by our earlier work (Bitsch et al. 2019). We can thus not make any statements about giant planet occurrence rates.

2.1 Pebble and gas accretion

Pebble accretion is modelled following the approach in Johansen et al. (2015). This pebble accretion approach takes into account the changes of the pebble accretion rates for planetary embryos on eccentric and inclined orbits (Johansen et al. 2015). This prescription reduces the pebble accretion rates for planetary embryos on eccentric and inclined orbits. The accretion rate of the planetary core Ṁ_core is directly proportional to the pebble surface density Σ_peb.

Pebbles in the protoplanetary disc settle towards the midplane depending on their size, parameterised in this work by the dimensionless Stokes number τ_f, and depending on the level of turbulence in the protoplanetary disc described through the disc’s viscosity α_disc. Using the α prescription, Youdin & Lithwick (2007) derived the pebble scale height H_peb using the gas scale height H_g via (1)

Here α_disc corresponds to the α-viscosity value inside the gas disc. Typically the pebbles in our simulations have Stokes numbers of 0.05–0.1, which is calculated by equating the drift time-scale with the growth time-scale (Birnstiel et al. 2012; Lambrechts & Johansen 2012; Bitsch et al. 2015b). That yields a value of H_peb∕H_g ~ 0.1, in agreement with observations of protoplanetary discs (Pinte et al. 2016). In order to be comparable to our previous works (Bitsch et al. 2019) we use an α_disc value of 0.0054 for the pebble scale height¹, which is different compared to the value used for the migration (see below). A lower α value used for the pebble scale height will result in a denser midplane layer of pebbles, with faster pebble accretion rates onto the planets as a consequence (Lambrechts & Johansen 2014; Bitsch et al. 2019).

The pebble surface density Σ_peb(r_P) at the planets location can be calculated from the pebble flux Ṁ_peb via (2)

where r_P denotes the semi-major axis of the planet, v_K the Keplerian velocity, and Σ_g(r_P) stands for the gas surface density at the planets locations. The pebble flux Ṁ_peb is calculated self consistently through an equilibrium between dust growth and drift (Birnstiel et al. 2012; Lambrechts & Johansen 2014; Bitsch et al. 2018a), where these simulations predict a decrease in the pebble flux in time. Here S_peb describes the scaling factor to the pebble flux Ṁ_peb to test the influence of different pebble fluxes, where we use S_peb = 2.5 and S_peb = 5.0, which is similar to our previous works (Bitsch et al. 2019; Izidoro et al. 2019) and as explained in (Bitsch et al. 2018a). Using S_peb = 2.5 a total of 175 Earth masses of pebbles drift through the disc in the 3 Myr of integration, which is doubled for S_peb = 5.0. The pebble sticking efficiency can be taken as ɛ_P = 0.5 under the assumption of near-perfect sticking (Lambrechts & Johansen 2014).

The Stokes number of the pebbles can be related to the pebble surface density Σ_peb and gas surface density Σ_g through the following relation (3)

Here η represents a measurement of the sub-Keplerianity of the gas velocity.

At higher temperatures, water sublimates and we fix the radius of the pebbles to 1 mm for T > 170 K, corresponding to typical chondrule sizes (Morbidelli et al. 2015; Ida et al. 2016). This is consistent with the assumptions made in Morbidelli et al. (2015) to explain the dichotomy between the terrestrial planets and the gas giants in the solar system. Additionally, we reduce the pebble flux Ṁ_peb to half its nominal value to account for water loss. In our disc model, the water ice line is located at ≈1 AU at the beginning of our simulations, but moves even further inwards in time as the disc evolves (Bitsch et al. 2015a).

As a planet grows, it starts to push away material from its orbit, generating a partial gap in the protoplanetary disc, where the planet generates an inversion in the radial pressure gradient of the disc exterior its orbit, halting the inward drift of pebbles (Paardekooper & Mellema 2006; Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018b; Ataiee et al. 2018). This planetary mass is referred to as the pebble isolation mass. The pebble isolation mass in itself is a function of the local properties of the protoplanetary discs, namely the disc’s viscosity ν, aspect ratio H∕r and radial pressure gradient ∂lnP∕∂lnr as well as of the Stokes number of the particles, which can diffuse through the partial gap in the disc generated by the planet (Bitsch et al. 2018b). We follow here the fit of Bitsch et al. (2018b), who gave the pebble isolation mass including diffusion of small pebbles as (4)

with λ ≈ 0.00476∕f_fit, , and (5)

where α₃ = 0.001.

As all planets accrete pebbles at the same time, the innermost planets will feel a reduced pebble flux, compared to the case of single planets in the disc. The pebbles accreted by the outer planets are subtracted from the pebble flux that arrives at the inner planets. Once a planet reaches pebble isolation mass, we set the pebble flux to zero for all the inner planets, stopping their growth by pebble accretion. We calculate the pebble accretion rate directly from the orbital position of the planetary embryo, so a planetary embryo on an eccentric orbit that is briefly exterior to an embryo that has reached pebble isolation mass could still accrete pebbles. However, the pebble accretion rates in our model decrease strongly if the planetary eccentricity increases, preventing pebble accretion of very eccentric planetary embryos (Johansen et al. 2015).

After the planet has reached pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018b; Ataiee et al. 2018), a gaseous envelope can contract (Lambrechts et al. 2014). Once the envelope becomes as massive as the planetary core, the planet can undergo runaway gas accretion (Pollack et al. 1996). We follow here the approaches of our previous works (Bitsch et al. 2015b, 2019). Even though gas accretion rates are heavily debated in the literature (Ayliffe & Bate 2009; Machida et al. 2010; D’Angelo & Bodenheimer 2013; Schulik et al. 2019) and span several orders of magnitude in range, we keep these rates as in our previous works to allow a better comparison.

After the planet has reached pebble isolation mass, it can contract its envelope, where the calculate the envelope contraction via (6)

Here, f is a factor to change the accretion rate in order to match numerical and analytical results, which is normally set to f = 0.2 (Piso & Youdin 2014). For the opacity in the envelope we use the fixed value of κ_env = 0.05 cm² g⁻¹, which is very similar to the values used in the study by Movshovitz & Podolak (2008). Lower and higher envelope opacities would result in higher and lower envelope contraction rates. As soon as M_core = M_env, the envelope contraction ends and rapid gas accretion can start.

For the rapid gas accretion, we follow Machida et al. (2010), who calculated the gas accretion rates via 3D hydrodynamical simulations in shearing boxes. The accretion rates are divided into two branches (7)

and (8)

where r_H denotes the planetary hill radius. These two branches divide low mass and high mass planets, where the effective accretion rate is the minimum of these two rates. In our approach, the gap opening does not affect the gas accretion rate, because the limiting factor of the gas accretion rate is what the disc can provide to the planet, which we set to 80% of the Ṁ value of the disc, because gas accretion is not 100% efficient (Lubow & D’Angelo 2006).

2.2 Planetary migration

Planetary migration in the type-I migration regime is modelled using the equations from Paardekooper et al. (2011). This prescription includes the Lindblad torques as well as the barotropic and entropy related corotation torque. Recently Jiménez & Masset (2017) introduced a new torque formula, which includes the same effects, but should be more accurate, because it is based on 3D hydrodynamical simulations in contrast to Paardekooper et al. (2011), which was based on 2D simulations. However, the changes compared to the Paardekooper et al. (2011) torque formula seem quite small in the pebble accretion scenario (Baumann & Bitsch 2020).

The accretion of material onto the planet can change the gas dynamics around it, leading to a thermal torque (Lega et al. 2014; Benítez-Llambay et al. 2015), which can, if the accretion rates are large, lead to outward migration (Benítez-Llambay et al. 2015). In addition, this effect could also increase the planetary eccentricity (Chrenko et al. 2017). However, Baumann & Bitsch (2020) showed that these effects only become very important if the accretion rates onto the planet are very large. In fact, we do not reach the accretion rates needed for the thermal torque to become positive and thus ignore its effects in this work.

Planets that become very massive and start to open deep gaps in the protoplanetary disc, change their migration regime to type-II. Kanagawa et al. (2018) relate the type-II migration time-scale to the type-I migration time-scale (which we calculate as explained above) in the following way (9)

where Σ_up corresponds to the unperturbed gas surface density and Σ_min to the minimal gas surface density at the bottom of the gap generated by the planet. The ratio Σ_up ∕Σ_min can be expressed through (Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2015) (10)

where (11)

We use here for the migration α_mig = 10⁻⁴, as this viscosity value was found in Bitsch et al. (2019) to allow giant planets to stay exterior to 1 AU. This choice of low viscosity in the midplane of the protoplanetary disc is motivated by the fact that recent disc observations point to low levels of turbulence (Flaherty et al. 2018; Dullemond et al. 2018). In addition, new evolution models of protoplanetary discs, where most of the angular momentum is transported away through disc winds instead of a large bulk viscosity, indicate a low midplane turbulence (Suzuki et al. 2016; Bai 2016; Chambers 2019). Finally, it should be noticed that detailed studies of planet migration in3D discs driven by disc winds are still lacking, even though 2D approaches have been made (Kimmig et al. 2020) indicating the importance of disc winds for planet migration that should be included in future models. Our migration prescription results in slow inward migration for massive planets due to the low viscosity, preventing large scale migration.

This type-II migration prescription is only valid in the case of low viscosities. At higher viscosities (e.g. α_mig> 0.001), the entropy driven corotation torque could operate (Paardekooper et al. 2011), which could allow outward migration in certain regions of the protoplanetary discs (Bitsch et al. 2015a). Applying in the case of high viscosity the migration prescription of Eq. (9) would lead to an unphysical outward migration of giant planets. At low viscosities, as we use in our simulations, the entropy driven corotation torque saturates (Baruteau & Masset 2008; Paardekooper et al. 2011), so that planets in our simulations always migrate inwards.

2.3 Damping of eccentricity and inclination

Damping of the orbital eccentricity and inclination during the gas disc phase tends to increase the pebble accretion efficiency and also to avoid orbital crossing and instabilities during the gas disc phase. For small planets undergoing type-I migration, we use the type-I damping rates of Cresswell & Nelson (2008). These damping rates are then applied to each low mass planets. Small mass planets only perturb the gas disc slightly, so that the damping formulae are still valid even if multiple small mass planets are present (Pierens et al. 2013). The exact implementation of the damping formulae for small planets in our code is given in Izidoro et al. (2019).

As soon as the planets start to open a gap in the disc, they start to migrate in type-II migration, where also the damping rates onto the planet are different. Changes of eccentricity and inclination of giant gap opening planets have been studied in the past (Papaloizou et al. 2001; Kley & Dirksen 2006), where Bitsch et al. (2013) derived damping rates for eccentricity and inclination based on 3D isothermal hydrodynamical simulation. However, all these works have only considered single giant planets in discs, so the effects of damping induced by the gas disc if multiple giant planets are present have not been studied, which is why we agnostically vary the damping rates for giant planets in our simulations.

The classical K-damping prescription (Lee & Peale 2002) relate the damping rates to the type-II migration rates through (12)

In the following, we use this damping prescription (Eq. (12)) rather than the damping prescription of Bitsch et al. (2013), because the K-damping prescription makes it easier to vary the damping rates. We use here K = 5, 50, 500, and 5000. A large K value implies a fast damping rate, meaning that any orbital eccentricity and inclinations that planets eventually get during the gas disc phase will be quickly damped to low values. We then explore how different damping values influence the formation of giant planet systems. As soon as the planet reaches a gap depth of 10% of the unperturbed gas surface density, we apply type-II damping (Eq. (12)) without any smoothing function from type-I damping towards the type-II damping rate.

We use in the following for |ȧ∕a| (in Eq. (12)) not the migration rates from Kanagawa et al. (2018), but instead use the classical type-II migration rate as soon as the gas surface density inside the opening gap reduces to 10% of the unperturbed gas surface density. Here we calculate ν = α_discH²Ω, with α_disc = 0.0054 corresponding to the α value used to compute the thermal structure of the protoplanetary disc (Bitsch et al. 2015a). In addition, we include the effect of the mass inertia, which slows down type-II migration once the planet becomes more massive than the disc (e.g. Baruteau et al. 2014). This implementation of the damping rates ensures that the damping rates are easier to control in our simulations, because they just depend on the disc’s viscosity. The dependency on the planetary mass due to the mass inertia only plays a small role in our simulations with the large pebble flux.

In Fig. 1, we show how the damping rates influence the evolution of eccentricity of a single Jupiter mass planet without growth as a function of time. The disc structure and 3 Myr gas disc lifetime is the same as in all our simulations (see above).

Fig. 1 Change of the eccentricity of a one Jupiter mass planet as function of time for the different damping parameters. The K damping values are derived using the classical type-II migration timescales. The flattening of the damping at the end of the disc lifetime is due to the exponential decay of the gas disc surface in the last 100 kyr, which reduces the gas surface density and thus the damping effects.

2.4 Initial conditions and simulation parameters

Initially the planetary embryos embedded in the disc have ~0.01 Earth masses, which corresponds to the pebble transition mass, where the accretion in the Hill regime becomes dominant (Lambrechts & Johansen 2012). The initial eccentricities and inclinations of the planetary embryo are randomly selected from uniform distributions and are 0.001–0.01 and 0.01–0.5°, respectively.This is also identical to our previous simulations (Bitsch et al. 2019).

The embryos are distributed radially between ≈3 and 17 AU, with equal radial spacing as in our previous work (Bitsch et al. 2019). In addition embryos starting interior to the water ice line are normally outgrown by their counterparts exterior to the water ice line (Izidoro et al. 2019). We test three different configuration with initially 15, 30 or 60 planetary embryos. As the total radial distance over which we spread the embryos is constant, the initial distances between the embryos varies for the different configurations.

After the gas disc phase of 3 Myr, we evolve the system until 100 Myr to study its long-term dynamical stability after gas dispersal. This effect is very important, as instabilities in the inner systems can occur several 10 Myr after gas disc dissipation, shaping the structure of the inner systems (Izidoro et al. 2017, 2019; Lambrechts et al. 2019).

As mentioned before, we test four different K-damping scenarios. For each K value we run 50 simulations, where we slightly vary the initial conditions regarding the initial planetary embryos mass, as well as the initial eccentricity and inclination and the orbital elements of the planets. In the simulations with S_peb = 2.5 we limit the growth of the planet by gas accretion to 1 Jupiter mass to investigate a scenario for Jupiter mass planets. Most of the giant planets formed in these simulations reach 1 Jupiter mass only at the end of the disc’s lifetime. We relax this condition for simulations with S_peb = 5.0, where planets can grow to larger masses by accretion of pebbles and gas. For simulations with S_peb = 5.0 we constrain ourselves to simulations with 30 initial planetary embryos. In addition we do not present the results of simulations with 60 planetary embryos, S_peb = 2.5 and K = 5. In total, we present here the results of 750 N-body simulations.

3 Individual systems

In this section, we discuss the results of two selected simulations that span different K damping values. We show an additional three planetary systems in Appendix A. All the simulations presented in this section use S_peb = 2.5, where the growth by gas accretion is limited to 1 Jupiter mass, but planets can grow more massive through mutual collisions.

In Fig. 2, we show the evolution of semi-major axis, planetary mass, eccentricity and inclination of a set of 30 planetary embryos using K = 50 as a function of time. The grey lines represent small bodies, while the black line represents larger bodies that are scattered away after the gas disc phase and the coloured lines mark the surviving planets.

Initially the planetary embryos start growing by accreting pebbles, which can be seen by the smooth increase in planetary mass at the beginning of the simulation during the gas disc phase. As soon as the planets reach the pebble isolation mass, they stop accreting pebbles and undergo a phase of envelope contraction before runaway gas accretion can start. Three planetary embryos reach pebble isolation mass quite fast, after only ≈500 kyr, but they remain very small until the end of the disc’s life time. This is caused by the fact that they finish their assembly in the inner regions of the disc, where the pebble isolation mass is small (Bitsch et al. 2015b, 2018b) and in agreement with the typical masses of super-Earths (Wu 2019; Bitsch 2019), preventing them from accreting a gaseous envelope. The envelope contraction phase is a strong function of planetary mass and the opacity in the planetary envelope (Ikoma et al. 2000; Movshovitz & Podolak 2008; Piso & Youdin 2014; Lambrechts & Lega 2017).

Only the planets accreting pebbles efficiently in the outer disc reach planetary core masses (of a few Earth masses) that allow a fast-enough envelope contraction for them to reach runaway gas accretion during the lifetime of the protoplanetary disc. This phenomenon was already observed in our previous works for single bodies (Bitsch et al. 2015b; Ndugu et al. 2018) and in our N-body framework (Bitsch et al. 2019).

During the gas phase of the disc about five planets with Saturnian masses (and larger) are formed. However, towards the end of the gas disc lifetime, the eccentricity and inclination are not damped efficiently any more due to the low gas density and the system undergoes a dynamical instability shortly after gas disc dissipation. As a consequence, only the two most massive gas giants survived on very eccentric and inclined orbits. The inclination and eccentricity of both planets oscillate in time. In particular the eccentricity of the inner gas giant (marked in red in Fig. 2) varies between ~0.07 and ~0.6. This behaviour has important consequences for matching the eccentricity distribution of the observed giant planets (see below). This behaviour has also been observed in many N-body simulations that deal with giant planets (e.g. Raymond et al. 2008, 2009a).

We show in Fig. 3 a simulations with K = 5000. This large K value prevents the build-up of eccentricities and inclinations during the gas disc phase for the growing planetary embryos. As some planets grow more and more, the eccentricities of the small bodies (below 1 Earth masses), increase dramatically towards the end of the gas disc lifetime. This then results with the ejection of the small bodies at this stage. Additionally, a collision between two gas giants occurs at the end of the gas disc lifetime, but the damping is so efficient that the eccentricity of the giants is damped almost immediately.

In the end a planetary system with several inner super-Earths and several outer gas giants forms. This result is similar to the results of Bitsch et al. (2019), where a different damping formalism was used (Bitsch et al. 2013), but effectively it corresponded to a large K value. The resulting system confirms that inner super-Earths (r < 1.0 AU) and cold Jupiters can form within the same system, as predicted by observation (Bryan et al. 2019; Zhu & Wu 2018). All planets feature low eccentricities and the system is nearly co-planar. We discuss this further in Sect. 7.

We note that the here presented simulations are just examples of the full set of simulations. Also simulations with small K can harboursystems of multiple giant planets (e.g. Appendix B) and not all systems with small K undergo dynamical rearrangements. Reciprocally, also systems with a large K value can undergo dynamical instabilities. We thus discuss in the next section about the statistics of our simulations.

Fig. 2 Evolution of a system with a damping factor of K = 50. Semi major axis (top left), planetary mass (top right), eccentricity (bottom left) and inclination (bottom right) of 30 planetary embryos as function of time. The gas disc lifetime is 3 Myr after injection of the planetary embryos, where the end of the gas disc lifetime is marked by the vertical black line. The grey lines represent either small mass bodies, or bodies that are ejected during the lifetime of the gas disc. The black lines represent massive planets that are scattered away after the gas disc phase. The coloured lines represent the two surviving planets.

4 Statistics and comparison to observations

In this section, we discuss the period ratio of adjacent planets, their orbital separation and the eccentricity distribution. We also discuss the number of survived planets in our simulations at the end of the gas disc lifetime and after 100 Myr of integration. In this section, we focus on the simulations using S_peb = 2.5 with planets that can grow up to Jupiter mass by gas accretion. We show the results of our simulations with S_peb = 5.0 in Sect. 5.

In our comparison to the observations we define giant planets as planets with masses above 0.5 Jupiter masses. Theinfluence of our results is not affected if we were to use a cut at around Saturn mass (about 1∕3 the mass of Jupiter), as this would increase the number of giant planets by a maximum of 5% for simulations with K ≥ 500, but much less for simulations with lower K. In addition, the influence is also smaller, if more initial planetary embryos were used, because the scattering effects in simulations with more initial planetary embryos are stronger (see below), where the lower mass bodies are removed efficiently. We want to stress here that our simulations are designed to form giant planets, so it is no surprise that all our simulations form giant planets. In addition this prevents us from saying anything about giant planet occurrence rates.

Fig. 3 Evolution of a system of 30 initial planetary embryos with a damping factor of K = 5000. The plots and lines have the same meaning as in Fig. 2. In this case, no giant instability after the gas disc phase happened within 100 Myr and a system with inner super-Earths and several outer gas giants survived.

4.1 Evolution of the planetary systems

In this section, we discuss the different aspects of the dynamical evolution of the planetary systems formed in our simulations. We focus on the separation of the planets (Fig. 4), the period ratio between adjacent planets (Fig. 5) and how many giant planets, defined as planets with masses above 0.5 Jupiter masses, are in each system (Tables 1 and 2).

In Fig. 4, we show the distances between adjacent planetary pairs in the simulations at the end of the gas disc lifetime (3 Myr) and at 100 Myr for the simulations with 60 initial planetary embryos for the different eccentricity and inclination damping values. The separation between the planets is shown in mutual Hill radii R_H,m defined as (13)

Here M₁ and M₂ are the masses of the two planets and a₁ and a₂ their semi-major axes. In the case of two planets on circular orbits, a separation larger than in units of the mutual Hill radii is sufficient to ensure that planets avoid mutual close encounters for all time (Gladman 1993). The distances between planetary pairs for the cases of 15 and 30 initial planetary embryos follows the same trend as for 60 planetary embryos.

Comparing the distances of the planets for the different simulations at 3 Myr (end of the gas disc lifetime) shows that some planets are closer to each other in units of mutual Hill radii than , indicating that these planets cannot avoid mutual encounters for all time (Gladman 1993). However, there are basically no planets withseparation of less than 5 mutual Hill radii, if only planets with masses larger than 0.5 Jupiter masses are considered, at the end of the gas disc lifetime at 3 Myr. This indicates a certain stability of the planetary systems.

Raymond et al. (2009a) studied the eccentricity evolution of giant planets by putting already fully formed giant planets at separations of 4–5 R_H,m to induce scattering events. In our simulations, only a small fraction of giant planets is within this separation at the end of the gas disc phase at 3 Myr and thus only a small fraction of our systems formed by pebble and gas accretion reflect the initial conditions used in Raymond et al. (2009a), which has important consequences for the eccentricity distribution of the giant planets (see below).

After 100 Myr of evolution, the systems have undergone scattering events, which increases the separation between adjacent planet pairs (dotted lines in Fig. 4). In particular, for K = 5 the planets are very widely spaced at 100 Myr, indicating more scattering events compared to the simulations with high K (fast damping). The mutual separations at 3 and 100 Myr look very similar for all three sets of simulations with initially 15, 30, and 60 planetary embryos.

In Fig. 5, we show the period ratios of adjacent planetary pairs in our simulations with different initial number of embryos (top to bottom) as well as for the different damping values. At the end of the gas disc lifetime, about 50% of adjacent planet pairs are interior of the 3:2 period ratio, with no substantial difference for the different damping values for all simulations. Our simulations show a slight pile up of planets around the 3:2 and 2:1 period ratios at this phase, but the clear majority of planets are clearly not locked in first order mean motion resonances.

In our simulations, mostly the inner super-Earths are trapped in resonant configuration at the end of the gas disc lifetime. The inner two pairs of super-Earths shown in Fig. 3 are actually in a 3:2 resonance configuration. However, our simulations do not show very long chains of super-Earths as in previous simulations dedicated to the formation of super-Earths (Izidoro et al. 2017, 2019; Lambrechts et al. 2019), where the chains can accommodate even up to nine super-Earths. In our simulations we have a maximum of 4–5 super-Earths per system. This difference is caused bythe slower migration rates used in the here presented work compared to Izidoro et al. (2019) and Lambrechts et al. (2019). In our simulations the planets that form in the inner regions of the disc and stay at super-Earth mass (e.g. Fig. 3) do not always migrate to the inner edge of the disc at 0.1 AU and can thus not always form a resonant chain by this mechanism. In addition, the super-Earth systems can be destroyed by gravitational interactions between the planets (see below). Furthermore, the most massive planetary cores in our simulations start to accrete gas and become gas giants.

During the next 100 Myr of evolution, the planetary systems undergo dynamical instabilities that increase the distances between adjacent planet pairs (Fig. 4) and as a consequence also the period ratios between planetary pairs. In fact now only a very small fraction of planet pairs are closer than the 2:1 period ratio.

Doing the same cut as before, by only looking at planets with masses larger than 0.5 Jupiter masses, reveals that only a tiny fraction of systems feature giant planets with a period ratio smaller than 2:1. However, our simulations show a pile-up of planets just exteriorto the 2:1 ratio. Nevertheless, there is no clear preferred period ratio between the giant planets in our simulations. In addition, our simulations suggest that pairs of giant planets in resonance should be rare. As already noted above, a lower K damping value results in more violent instabilities, also pushing the period ratios between adjacent planet pairs to larger values.

In Table 1, we show the average number of planets with masses larger than 0.5 Jupiter masses at the end of the simulation at 100 Myr. The average number of planets slightly increases from K = 5 to K = 500, but then drops off slightly again for K = 5000.

In our simulations, we observe a complex interplay between the efficiency of eccentricity and inclination damping on the one hand and the number of planetary embryos on the other hand. The interplay between these two variables then influences how stable the formed planetary system is and thus how many giant planets survive. However, from within our simulations it is difficult to observe a very clear trend regarding the distances between the planets, their period ratios or the number of surviving planets.

Table 2 shows the average number of planets with masses larger than 0.5 Jupiter masses that could be detected by RV measurements with 1 m s⁻¹ and a distance up to 5.2 AU. Therefore, the trends of the total average number of giants with masses larger than 0.5 Jupiter masses is reflected as well. But, of course, the total number of planets detected by the RV measurements is much lower. In fact the simulations with slow damping (low K) predict an average number of less than 1.5 giant planets per system. This implies that currently observed planetary systems should have one or two giant planets, but should host on average around 2–3 planets according to our simulations (Table 1), where the maximal number of giant planets in our simulations is 7.

Most of the giant planets in our simulations are within 10 AU, they are unlikely to be observed by direct imaging. However, our simulations predict that if systems harbouring giant planets are observed for longer time by RV measurements, extending the orbital distance to which planets can be found, the average number of giant planets per system should increase by roughly 50%.

Fig. 4 Distances of adjacent planetary pairs of all planets in the simulations starting with 60 planetary embryos (including the small bodies) after 3 Myr (end of the gas disc phase, top) and 100 Myr (dashed, bottom). In addition, we show the distances after 3 Myr only for planets larger than 0.5 Jupiter masses (long-dashed, middle panel), to indicate how the stability of the system is influenced by the larger bodies. The black vertical line marks a separation of mutual Hill radii, at which two planets on circular orbits avoid mutual encounters for all time (Gladman 1993). We only show the separation in units of mutual Hill radii for simulations with initially 60 bodies. The trends are very similar for simulations with 15 and 30 initial embryos.

Fig. 5 Period ratios of adjacent planetary pairs of all planets in the simulations (including the small bodies) after 3 Myr (end of the gas disc phase) and 100 Myr (dashed). In addition, we show the period distribution including only planets with masses larger than 0.5 Jupiter masses after 100 Myr (long-dashed lines). We do not observe a significant pile-up of planets close to mean-motion resonances, in contrast to the super-Earth simulations of Izidoro et al. (2019), which we attribute to the slower migration speed in the here used simulations, which prevents the formation of resonant chains at the inner edge of the protoplanetary discs. Nevertheless, the systems become unstable and the instabilities result in a wide spacing of the planets after 100 Myr.

4.2 Eccentricity distribution

In Fig. 6, we show the eccentricity distribution for the planets formed in our simulations with different K and different initial number of planetary embryos for planets that could be detected by RV measurements with a sensitivity of 1 m s⁻¹ and up to 5.2 AU with masses above 0.5 Jupiter masses. In addition, we show the eccentricity distribution of the giant planets observed via RV fromthe exoplanet.eu database. We show here only planets with a minimum mass of 0.5 Jupiter masses and maximal M sin(i) of 1.25 Jupiter masses for the observations, but no upper mass limit of planets for our simulations, which barely reach masses larger than 1 Jupiter masses for this set of simulations. In addition, we limit ourselves to planets with orbital distance larger than 0.1 AU, in order to avoid a contamination by hot Jupiters, which are abundant in their observations due to their detectability, but are not very common in general (Mayor et al. 2011; Howard 2013; Dawson & Johnson 2018). In all our simulations we only produce about 1% of hot Jupiterswith semi-major axis smaller than 0.1 AU, which are also excluded from the following analysis, because we do not include effects of tides in our simulations.

Clearly, larger K values result in an eccentricity distribution that is too steep to explain the observations. For K ≥ 500, around ~70% of all giant planets with masses larger than 0.5 Jupiter masses have an eccentricity below 0.1. This clearly does not match the eccentricity distribution from the observations. In the case of 15 initial planetary embryos, this is also true for K = 50 and K = 5. However, when increasing the initial number of planetary embryos, the number of planets with eccentricities below 0.1 drops to around ~45% or below for K = 50 and K = 5, resulting in an eccentricity distribution closer to the observations for that particular eccentricity bin.

On the other hand, the simulations with K = 50 under predict the eccentricity distribution of the observations for 0.1 ≤ e ≤ 0.3, while they seem to slightly over predict the frequency of giant planets with e > 0.3 in the case of 30 or 60 initial planetary embryos. In the case of K = 5, this effect becomes even more prominent. Our simulations actually under produce planets with low eccentricities and predict more planets with larger eccentricities. This is caused by the slow damping during the gas disc phase which allows the planets to acquire already some small eccentricities, leading to instabilities at the end of the gas disc phase.

The comparisons in Fig. 6 already give a good clue that low K values are needed in our simulations to reproduce the eccentricity distribution of the RV observations. It seems that a K value between 5 and 50 seems to give the best results regarding the eccentricity distribution. We show in Table 3 the results of a Kolmogorov–Smirnov-test (KS-test) for the comparison between the eccentricity distribution in our simulations with the eccentricity distribution of the observed giant planets. The higher the value of the KS-test, the better the match of the simulations with the observations. Clearly, a K value between 5 and 50 seems to match the observations best. However, this does not contain the radial distribution of the planets and if the giant planets with large eccentricities originate from close-in or far away orbits.

In Fig. 7, we show the eccentricity of the planets as function of their semi-major axis for the same set of simulations as in Fig. 6 and for the same observations. However, we limit ourselves to K = 5 and K = 50, because these simulations match the eccentricity distribution of the observations best. In addition, we show the average eccentricity for five semi-major axis bins.

Our simulations allow a quite nice match to the observed eccentricity distribution of the giant planets. For K = 5, our simulations show a slightly too large frequency of giant planets with large eccentricity (see also Fig. 6). From Fig. 7, it is clear that especially planets at large orbital distances (r ≥3 AU) have a larger mean eccentricity compared to the observations for 30 or 60 initial planetary embryos. In the case of 15 initial planetary embryos the mean eccentricity at these large distances is too low compared to the observations, which is related to the reduced number of scattering events for these simulations.

Our simulations clearly show that either larger K values are inappropriate (see also Lee & Peale 2002) or the planetary masses in our simulations are not large enough to allow sufficient scattering to match the eccentricity distribution of the giant planets observed via RV detections. We relax the later assumption in the next section, where planets are allowed to grow faster and bigger. In addition, our simulations show that the initial number of planetary embryos does not play a significant role in the final number of giant planets in simulations using a slow damping rate (small K), but only seems to become slightly more important at fast damping rates (large K). We discuss the biases and limits of the observation of eccentricity giant planets in Sect. 7.

Fig. 6 Eccentricity distribution of the simulated planets within a distance of 5.2 AU from their host star with masses of at least 0.5 Jupiter masses as well as of RV observations of giant planets with masses between 0.5 and 1.25 Jupiter masses exterior to 0.1 AU. We show in colour the different damping rates and systems with initially 15, 30 and 60 embryos (from top to bottom).

Fig. 7 Eccentricity as a function of semi-major axis of simulated planetary systems with different damping rates anddifferent number of initial embryos (15, 30, 60 from top to bottom). We only show the planets with masses of at least 0.5 Jupiter masses (for simulations and observations) and only for the K values that match the observations of the eccentricity distribution best (Fig. 6). The horizontal lines depict the mean eccentricity within each radial bin. The black symbols show the data from exoplanet.eu with planetary masses ranging from 0.5 to 1.25 Jupiter masses, which is roughly the masses the planets in our simulations reach.

5 Fast growth

In this section, we relax the previous restrictions to our simulations and allow planets to grow faster (S_peb = 5.0) and beyond 1 Jupiter mass by pure gas accretion. However, we limit ourselves here to simulations with 30 initial planetary embryos, but vary the damping rates. We show the evolution of such a system in Fig. 8.

By comparing the simulation shown in Fig. 8 with the simulation shown in Fig. 2, a few differences become immediately visible. In the case of S_peb = 5.0, the planets accrete pebbles more efficiently and some planets reach pebble isolation mass well before 1 Myr. As the planets grow faster, they reach a larger pebble isolation mass (Bitsch et al. 2015b; Bitsch 2019), which in turn allows them to contract their envelope in a shorter time. This result in the planets starting to reach runaway gas accretion at about 1 Myr. In contrast the planets growing in a disc with lower pebble flux (S_peb = 2.5) planets start their runaway gas accretion at ≈2 Myr (Fig. 2).

The faster growth of the planets around 1 Myr of evolution results in stronger interactions between the planets, where the eccentricity is growing significantly. This effect is enhanced due to the slow damping (K = 50) by the gas disc, which is too slow to keep the eccentricities small. As a consequence, the system undergoes an instability after about 1.3 Myr, where then only three giant planets survive and the remaining planetary embryos are ejected from the system already during the gas phase.

After the scattering event, the eccentricities and inclinations of the remaining giant planets are damped by the interactions with the gas. This results in an eccentricity of around 0.05–0.1 for the inner two planets and of 0.65 for the outer giant planet. During the damping phase, the giant planets also continue to grow, where the inner planets reach around 3 Jupiter masses each and the outer planet is slightly below Jupiter mass. The outer planet grows very slowly, because of the low gas densities in the outer regions of the disc.

After the end of the gas disc phase, the eccentricity and inclinations of the inner planets oscillate, which has important consequences for the comparison with observations (see Sect. 7). These oscillations are also very common in the simulations with S_peb = 2.5, independently of the initial number of embryos.

In Fig. 9, we show the orbital separation between adjacent planets (top) and their period ratios (bottom) for all our simulations with S_peb = 5.0. A clear difference with respect to the simulations with S_peb = 2.5 is visible (Figs. 4 and 5). In the simulations with S_peb = 5.0 only a very tiny fraction of planets is closer than 5 mutual Hill radii at the end of the gas disc phase at 3 Myr, independently of the damping rate of eccentricity and inclination. As mentioned before, a distance of 4–5 mutual Hill radii was used in Raymond et al. (2009a) to start out their simulations in order to achieve scattering events in short times. In addition, there is no large difference if only the separations between planets with masses larger than 0.5 Jupiter masses is taken into account.

Due to the more rapid growth, the planets grow quickly and become very massive, allowing instabilities to happen already during the gas disc phase. As a result the orbital separations between the planets are increased. The surviving giant planets migrate in the slow type-II migration regime (Eq. (9)) preventing them to migrate close to each other. In particular, in the case of K = 5, the instabilities during the gas disc phase are very common and 90% of the planets have separations larger than 7 mutual Hill radii already towards the end of the gas disc phase. On the other hand, as already seen in the simulations with S_peb = 2.5, fast damping results in more tightly packed systems at the end of the gas disc phase, because the efficient damping reduces the eccentricities of the planets and thus allows planets to be close to each other². Here, the faster growth combined with the efficient damping (K ≥ 500) leads to very stable systems at wide separations.

As the majority of the instabilities for simulations with K ≤ 50 already happened during the gas disc phase, the mutual separation between the planets does not change significantly after 100 Myr of evolution. For the simulations with K ≥ 500, we also do not observe a large change in the separations between the planets after 100 Myr. This is probably caused by the fact that the planets are already further away than five mutual Hill radii from each other at the end of the gas disc phase, preventing efficient scattering on timescales of 100 Myr (Chambers et al. 1996).

This is also reflected in the period ratios between the planets (bottom in Fig. 9), where the number of planetary systems that host planets interior to the 2:1 period ratio is very small, especially when compared to the simulations with S_peb = 2.5 (Fig. 5). This is also caused by the fact that the planets grow faster and bigger, which reduces their migration speed due to the deeper planetary gap (Eq. (10)), and thus prevents them to come close to each other.

In agreement with the separation in mutual Hill radii, the period ratios of the systems do not change significantly when including only planets with masses larger than 0.5 Jupiter masses or after 100 Myr. In Fig. 9, we plot only up to period ratios of 5:1 for visibility reasons. The period ratios between the planets are largest for the simulations featuring K = 5. This is caused by the stronger scattering interactions between the planets, which separates the planets further and thus increases their period ratios.

In Table 1, we show the total number of giant planets with masses larger than 0.5 Jupiter masses in our simulations and the number of these giant planets that are detectable with RV measurements of 1 m s⁻¹ up to 5.2 AU is shown in Table 2. In comparison to the simulations with S_peb = 2.5, the simulations with S_peb = 5.0 and large K can harboursystems with more giant planets (see also Table 1). This can be explained through a larger pebble flux through the disc, which causes the planets to grow faster. This implies that more planets start initially to grow and are thus more resistant to the instabilities that will later follow at the late stages of the gas disc phase where the small bodies are ejected from the system and the larger ones remain in the disc due to the efficient damping.

In the case of slow damping (small K), the average number of giant planets is quite similar for both pebble fluxes. This could imply that the growth does not necessarily play the major role in determining the final structure of the planetary system, but the damping rates of eccentricity and inclination once the planets become massive and open gaps. In the case of slow damping, we thus only observe small variations in the number of planets that could be found via RV detections for the different pebble fluxes (Table 2).

In Fig. 10, we show the mass distribution of the giant planets formed in our simulations using S_peb = 5.0 and those of the RV observations with masses larger than 0.5 Jupiter masses, but limited to five Jupiter masses. We chose this upper limit, because planets with larger masses could predominantly be formed via gravitational instabilities (Schlaufman 2018).

The giant planets in our simulations have an average mass larger than inferred from the observations, indicating that the growth in our simulations might be too efficient compared to reality. In addition, the planets formed in the simulations with slow damping (low K) seem to be heavier than their counterparts formed in the simulations with fast damping (high K). This difference could originate from planet–planet collisions, which are more frequent in the simulations with slow damping. However, even thoughthe mean mass is higher in our simulations compared to the observations, it seems that our simulations under produce planets with masses larger than 2.5–3 Jupiter masses. This could be related to the limited disc lifetime in our simulations, where the planets only grow 3 Myr in a gas disc environment, however, discs in reality can live longer (Mamajek 2009). Planets growing in discs that live longer have more time to accrete gas and could thus grow bigger. In addition, disc masses are quite widely spread (Andrews et al. 2013), which would provide more material for the planets to grow. In addition, planets forming by gravitational instabilities could populate this area of parameter space (Schlaufman 2018).

In Fig. 11, we show the eccentricity distribution of the giant planets with masses larger than 0.5 Jupiter masses in our simulations with S_peb = 5.0 as well as ofRV observations taken from exoplanet.eu with masses between 0.5 and 5.0 Jupiter masses. Compared to Fig. 6 we now include planets with larger masses for the RV observations.

Clearly, a fast damping of eccentricity and inclination (large K) results in a very small number of instabilities and thus in a very steep eccentricity distribution, similar to the simulations with S_peb = 2.5. This clearly indicates that also planets with larger masses cannot compensate for the efficient damping in this case. In particular, due to the larger masses, the planets migrate slower and have larger mutual distances at the end of the gas disc lifetime (Fig. 9) compared to the planets formed in the simulations with S_peb = 2.5, which prevents most scattering events after the gas disc phase in the 100 Myr of system evolution, resulting in systems with basically no significant eccentricity.

On the other hand, the planets formed in simulations with slow damping (small K) show a significant eccentricity distribution. For K = 50, our simulations slightly over predict the number of planets with e < 0.1, match quite nicely for 0.1 < e < 0.3, but under predict the number of planets at larger eccentricities in agreement with the simulations using S_peb = 2.5. However, the number of planets with e < 0.1 is larger for S_peb = 5.0. We again attribute this to the faster growth, resulting in slower migration and thus in wider separation of the giant planets at the end of the gas disc phase (Figs. 9 and 4), resulting in more stable systems.

But for K = 5, our simulations show an under prediction of planets with eccentricities smaller than 0.3, but over predict planets with e > 0.3 compared to the observations in line with the simulations using S_peb = 2.5. Even though the planets have wide separations at the end of the gas disc lifetime (Fig. 9), they feature a significant eccentricity distribution, which is caused by the very inefficient damping, resulting in instabilities during the gas disc phase.

We show the eccentricities of the planets with masses larger than 0.5 Jupiter masses formed in our simulations as function of their orbitaldistance in Fig. 12. We show also the eccentricities of giant planets of masses between 0.5 and 5.0 Jupiter masses detected via RV. In addition, we show the mean eccentricity for five orbital distance bins. As we include now more massive planets from the observations, the mean orbital eccentricity in the outermost bin increases significantly compared to Fig. 7. This change in the mean eccentricities is caused by the larger sample of planets, which now include more massive planets. Our sample now includes 287 planets, in contrast to the 85 planets in the mass range of 0.5 to 1.25 Jupiter masses used for Figs. 6 and 7. We discuss this more in Sect. 7.

The differences compared to the simulations with S_peb = 2.5 can already be inferred from our previous discussion. For the fast damping rates (large K), the eccentricities of the giant planets are lower compared to the simulations where planets grow with S_peb = 2.5 (Fig. 7). It is very clear that these fast damping rates fail to reproduce the eccentricity distribution of giant planets inferred from observations.

Our simulations suggest that a K value between 5 and 50 is needed to match the eccentricity distribution of the observations exterior to 1 AU. This also becomes apparent from a K–S test (see Table 3), where only the simulations with low K give a reasonable answer to match the observed eccentricity distribution. When mixing the simulations of K = 5 and K = 50 in a 3:1 fashion, we get a p-value of 0.39, indicating that the real K for these kind of simulations should be in between 5 and 50.

On the other hand, our simulations with S_peb = 5.0 clearly under predict the eccentricities interior of 1 AU. We attribute this to the low number of planets within 1 AU in our simulations. In fact the 50 simulations using K = 5 only show one planet within 0.5 and 1 AU. On the other hand, the simulations with S_peb = 2.5 and K = 5 show a nice match to the observations within 1 AU (Fig. 7).

In addition, the damping of eccentricity and inclination depends on the gap that giant planets open (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013). The gap opening in turn depends on the disc’s viscosity and aspect ratio (Crida et al. 2006; Kanagawa et al. 2018), which is smaller in the inner disc for flared disc profiles as we use here. This indicates that the faster growing planets will open their gap earlier, especially in the inner regions. As the type-II damping rate is slower than the type-I damping rate, this enhances the eccentricities of the growing planets, leading to more scattering events eventually depleting the number of giant planets in the inner regions for our simulations. On the other hand, the disc’s viscosity is set by the disc’s turbulence, which varies with orbital distance due to the operation of different instabilities (e.g. Pfeil & Klahr 2019), indicating that taking radially varying profiles of viscosity into account could alter the damping rates and thus the interactions of growing planets.

It seems that the faster growth rate and thus more massive planets have only a significant impact on the eccentricity distribution if the damping of eccentricity and inclination is slow (low K). However, the overall evolution of the planetary systems is very different. In the case of slow growth, the instabilities happen towards the end of the disc lifetime, when the damping forces reduced due to the decay of the disc. On the other hand, in the case of fast growth, the instabilities can already happen in the middle of the gas disc lifetime, with important implications for the evolution of the system and how these systems could be observed, as we discuss in detail in Sect. 7.

Fig. 8 Evolution of a system with a damping factor of K = 50 where the pebble flux is twice as large as in Fig. 2. The lines have the same meaning as in Fig. 2. Due to the massive planets, an instability already occurs during the gas phase of the disc, increasing the eccentricity of the planets. However, during the remaining gas phase of 1.5 Myr after the instability, the eccentricities of the planets get damped efficiently and the two inner planets have eccentricities of around 0.1.

Fig. 9 Distances (top) and period ratios (bottom) of planets formed in our simulations with Speb = 5.0 and different damping values. For this plot we cut the period ratio at 5.

Fig. 10 Mass distribution of the planets formed with masses larger than 0.5 Jupiter masses in our simulations and from RV observations, where we limit the observations to 5.0 Jupiter masses, as this is roughly the largest masses we reach in our simulations.

Fig. 11 Eccentricity distribution of the planets formed in our simulations with faster growth (Speb = 5.0) and the observations with planetary masses from 0.5 to 5 Jupiter masses.

6 Super-Earths – cold Jupiter relation

In inner regions of planetary systems (within 1 AU), super-Earths seem to be the most common type of planets, where 30–50% of all systems appear to harbour a super-Earth (Mulders et al. 2018). On the other hand, only 10–20% of all planetary systems should host a cold gas giant beyond 1 AU. In classical formation scenarios, the super-Earths might be the failed cores that did not make it to become giant planets (Cossou et al. 2014).

RV Observations by Barbato et al. (2018) investigated if systems hosting cold gas giants also host inner super-Earths. However, within their sample of 20 systems they did not find any super-Earth, resulting in the conclusion that at maximum 10% of cold Jupiter systems should harbour also inner super-Earths. However, Zhu & Wu (2018) and Bryan et al. (2019) concluded opposite, namely thatcold Jupiters should have inner super-Earth in up to 90% of the cases. The difference to Barbato et al. (2018) is explained inZhu & Wu (2018) through the detection limits by the RV survey that can only observe planets down to 15 Earth masses, while the Kepler satellite, whose data was used in Zhu & Wu (2018) and Bryan et al. (2019), can find much smaller objects.

The simulations by Izidoro et al. (2015) showed that gas giants that are formed interior to super-Earth can act as barriers to the faster migrating super-Earths. This scenario is built on the assumption that gas giants interior to super-Earths exist. In the solar systemcontext, this could explain the formation of Uranus and Neptune via collisions between the outer super-Earth cores (Izidoro et al. 2015).

Here, our model shows the opposite formation pathway, namely that inner super-Earths form before the giant planets emerge from gas accretion (e.g. Fig. 3) and interior to them. We now discuss within our simulations how many systems harbour outer giant planets and inner super-Earths.

We show in Table 4 the fraction of systems from our simulations that harbour outer gas giants with inner super-Earths. We define here planets that are dominated in mass by their core as super-Earths in our simulations. Our simulations show that for the fast damping cases a significant fraction of giant planet systems should host inner super-Earths or if the initial number of planetary embryos is small. The number of systems hosting inner super-Earths and outer gas giants clearly decreases when the initial number of planetary embryos is larger. This is related to the dynamical history of the system. As discussed in the previous section, a faster damping rate (larger K) and a lower number of initial planetary embryos results in planetary systems that are less likely to undergo dynamical instabilities, which then keep the inner systems intact.

On the other hand, the simulations that show a 30–40% occurrence of inner super-Earths and outer gas giants, feature 1.5–2.8 gas giants per system that could be detected via RV measurements (Table 2). Nevertheless, most observed systems with inner super-Earths feature only one gas giant, indicating a bit of discrepancy between the observations and our models. Of course, some systems like HD 160691 and HD 34445 actually feature two gas giants with orbits larger than 1 AU and with inner super-Earths, but these seem to be the exception.

Planetary systems formed in our simulations with slow damping (small K) only form a small fraction of systems with inner super-Earths and outer gas giants for the case of 30 or 60 initial planetary embryos, contradicting the super-Earth cold Jupiter relationship from observations. This is related to the dynamical instabilities that happen in these systems, where giant planets become eccentric and thus destroy the inner planetary systems, as has been shown also by pure N-body simulations (Mustill et al. 2015). On the other hand, the simulations with slow damping reproduce the eccentricity distribution of the giant planets very well. It seemed that fast damping of eccentricity and inclination is favourable to explain the super-Earth cold Jupiter relation in our simulations, but these same simulations clearly fail to reproduce the eccentricity distribution of the giant planets.

We propose a few ways to solve this apparent mystery. In the planetary systems formed from simulations with S_peb = 5.0 and K = 5 and K = 50 (Figs. B.1 and B.2), it is evident that systems of inner super-Earths mostly exist in systems where the outer gas giants have a very low eccentricity, implying that these systems have not undergone major scattering events. Our simulations feature also a significant fraction of warm Jupiters (r < 1 AU), which could have a significant influence on the stability of the inner systems once scattering events take place.

The study by Zhu & Wu (2018) suggested that up to 90% of the cold Jupiter population should have inner super Earths. However, our simulations show that violent scattering events are needed to explain the eccentricity distribution of the observed cold Jupiter population, but these scattering events destroy the inner planetary systems, so that giant planets on eccentric orbits should not harbour inner super-Earth systems. One exception in our simulations is run three in Fig. B.1 which features an inner super-Earths and a giant planet with an eccentricity of around ~0.6 at the end of the 100 Myr evolution.

We show the eccentricity distribution of the observed giant planets used in our study and those used by Barbato et al. (2018) and Zhu & Wu (2018) in Fig. 13. It is clear that the eccentricity distribution of the giant planets in the sample of Zhu & Wu (2018) features eccentricities that are on average lower than the eccentricities of the giant planets extracted from exoplanet.eu, when taking their selection criteria (the planets should have larger distances than 1 AU to their central star) into account (dashed black line).

In addition, two systems, HD 219828 and HD 125612 (both within the sample of Zhu & Wu 2018) feature highly eccentric giant planets with inner super Earths. The giant planets in these systems have several Jupiter masses. Massive planets like this can actuallyrequire their eccentricities through interactions with the protoplanetary disc itself (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013), potentially explaining their eccentricities without scattering events and thus leaving their inner super-Earth systems undisturbed. Excluding these two systems from their sample leads to a much steeper eccentricity distribution in the Zhu & Wu (2018) sample. In this case (blue dashed line vs. black dashed line), there seems to be a factor of ~2 difference inthe frequency of giant planets with eccentricities below 0.1 compared to the observations, which could significantly influence the survival rate of inner super-Earth systems as our simulations show.

The data used by Zhu & Wu (2018) suggests that giant planets on eccentric orbits (e > 0.1) do not harbour any super-Earth exterior to a 10-day orbit (≈0.1 AU), while giant planets on close to circular orbits can harbour super-Earths with much longer periods. Planets that are closer to the central star are deeper anchored within the stellar potential well and are thus harder to eject by exterior giant planets, making these planets prone to survive instabilities easier than their wider orbit companions. However, in our simulations the type-I migration speed is quite slow (due to the low surface density and low viscosity), so that super-Earths barely reach 0.1 AU, which could result in basically no survivors in the inner systems once the outer giant planet system becomes unstable. We plan to test this hypothesis in future work.

Of the 20 planetary systems observed in the work by Barbato et al. (2018), only 6 systems feature gas giants with an eccentricity below 0.1, while 7 planets feature eccentricities larger than 0.2. Our formation simulations suggest that the search for inner super-Earth planets in the majority of the systems observed by Barbato et al. (2018) was doomed from the start formost of their systems, potentially explaining why they did not find any inner super-Earths or mini Neptunes.

Alternatively, the systems of inner super-Earths could form late and after the scattering events. Our simulations with S_peb = 5.0 and slow damping (low K) suggest that the scattering events happen early and can actually take place during the gas disc phase. If planetesimals survive in the inner disc³, these could form systems of inner super-Earths after the scattering events similar to the terrestrial planets in our solar system that finished their formation after the gas disc phase (Raymond et al. 2009b). However, this scenario requires that enough material is present in the inner disc to allow efficient planet formation.

Nevertheless, the results of our simulations suggest that systems with giant planets on very eccentric orbits should harbour only very close-in super-Earths, if at all. In addition, the more massive the outer giant planet is⁴, the less likely it is that the inner super-Earth survives the scattering event that gives the outer planet its eccentricity. Reversely, the less massive and less eccentric an outer giant planet is, the more likely it is that the inner super-Earth system stayed intact. In addition, the close a giant planet is to the inner super-Earth region, the less like it is for super-Earths to exist within the same system (Figs. B.1 and B.2). If the formation channel presented in this work (pebble accretion, gas accretion and planet migration of planetary embryos forming in the outer disc) is correct, it implies that searching for terrestrial planets in systems with outer gas giants will only be of success if the outer gas giants have very low eccentricities and that an analysis based on stability limits alone without taking the formation channels into account (Mustill et al. 2015; Agnew et al. 2018; Kokaia et al. 2020) might be misleading.

Our simulations with surviving inner super-Earths and outer gas giants harbour a maximum of three inner super-Earths in the case of slow damping, but can have up to 4–5 super-Earth in case of faster damping of eccentricity and inclination for giant planets. Chains of more inner super-Earths did not form within our simulations. However, in the simulations of Izidoro et al. (2019) chains with up to 9 super-Earths can form in the inner disc, indicating that systems with many super-Earths should not harbour outer giant planets, if both (super-Earths and gas giants) form in the outer disc and migrate inwards.

Fig. 12 Eccentricity–semi-major axis distribution of the observed giant planets (black triangles) with masses of 0.5 < MP < 5.0MJ and our simulations with different K-damping values (coloured circles) after 100 Myr of evolution. The horizontal lines show the mean eccentricity within each orbital distance bin. While our simulations with low K reproduce the e − a distribution exterior to 1 AU, the simulations show too small eccentricities for the planets interior to 1 AU. On the other hand, this might also be related to the limited number of giant planets within this region in our simulations. For planets in the order of Jupiter mass, the low K values match quite nicely (Fig. 7), due to the slightly larger number of planets within that semi-major axis bin in our simulations.

Fig. 13 Eccentricity distribution of the sample of exoplanet.eu with masses larger than 0.5 Jupiter masses where the planets have a minimal distance to their host star of 0.1 AU (solid black line) or 1.0 AU (dashed black line). In red we show the eccentricity distribution of the planets observed in Barbato et al. (2018), while we show eccentricity of the giant planets used in the study of Zhu & Wu (2018) in blue. In the dashed blue solid line we also show the sample of Zhu & Wu (2018), but exclude two planetary systems where the companions are more massive than 5 Jupiter masses, which could inherit their eccentricity from planet–disc interactions.

7 Discussion

In this section, we discuss the shortcomings of the eccentricity observations of giant planets as well as of our model. In addition, we relate our results to previous works and discuss about the implications of our results for protoplanetary disc and exoplanet observations.

7.1 Eccentricity distribution of observations

Determining the true eccentricity of giant planets from RV observations is quite difficult. From the orbital analysis of RV measurements one can determine the minimum planetary mass, the semi-major axis and the eccentricity. The success of this approach depends on the RV precision, temporal baseline of the data, their number, and also on the magnitude of the Doppler signal induced by the planet.

In particular, if only a few RV measurements have been taken, the orbital fit might result in a too large eccentricity of the planet. By taking more RV measurements of already known planetary systems, the orbital fit can be improved. It has been shown that the RV signal of single eccentric planets with only a few measurements could actually be caused by two giant planets in resonance (Wittenmyer et al. 2019a,b). This could apply to up to 30% of the giant planet systems detected via radial velocities, where the single planet has an eccentricity below 0.5. However, our simulations do not show a significant pile-up of giant planets close to first order mean motion resonances (Figs. 5 and 9). Nevertheless, the impact on the overall eccentricity distribution of exoplanets might be small, but would imply that slightly faster damping rates (slightly larger K values) might be needed to model the eccentricity distribution of giant planets.

In our simulations, the eccentricity of the giant planets are caused by scattering events either during or after the gas disc phase, which was also observed in hydrodynamical simulations (Lega et al. 2013). However, if more than one giant planet survives, the remaining giant planets can still interact, resulting in a variation of the planetary eccentricity and inclination in time (e.g. Fig. 2). The oscillations happen normally on times of several thousand (or much more) orbits. On the other hand, the observations of seemingly single giant exoplanets span a maximum of a few orbits, with usually no observational evidence of additional longer-period planets that may be part of the system. The orbital parameter uncertainties, together with the unknown physical and orbital properties of the hypothetical companions, make it impossible to estimate themagnitude of the mutual orbital perturbations, and thereafter, impossible to determine if the measured eccentricity of the observed giant planet is currently on its minimal, maximal, or on any other phase of the eccentricity evolution.

7.2 Model assumptions

Our model of planet formation is based on pebble accretion, gas accretion, protoplanetary disc evolution as well as on planet migration. In addition, our model assumes that the planetary embryos are already fully formed at the beginning of the simulations. All these ingredients have many assumptions flowing into it, which we briefly discuss here. We do not attempt to vary all parameters in our model, because our work is a proof of concept and the here presented simulations are designed to form giant planets.

In our simulations, we have used the disc model of Bitsch et al. (2015a), where the thermal structure is derived from just micrometre sized grains. In reality, these grains can grow through coagulation (Brauer et al. 2008) and condensation (Ros & Johansen 2013; Ros et al. 2019), while their maximal size is limited by fragmentation, leading to a size distribution of these grains. The first models of hydrodynamical simulations including a full grain size distribution to calculate the heating/cooling effects have just become available (Savvidou et al. 2020) and will be used in future simulations. We start our N-body simulations in a disc that is already two Myr old Bitsch et al. (2015a). We chose this approach to be consistent with our previous work (Bitsch et al. 2019). The disc’s aspect ratio in the inner disc reduces quite fast due to the reduction in viscous heating, resulting in a small pebble isolation mass in the inner disc in line with the masses of super Earths (Wu 2019; Bitsch 2019).

Recent disc evolution models including the effects of disc winds, show different radial profiles, where the gas surface density can generate pile-ups around 1 AU (Suzuki et al. 2016; Bai 2016; Chambers 2019). This could prevent the inward migration of planetary cores, generating a pile-up of planetary cores that then undergo runaway gas accretion in the outer regions of the disc. As these gas giants are in the outer disc, their scattering events could disturb the inner super-Earth population less compared to gas giants forming in the inner disc. A similar effect could be achieved by photoevaporation, which can carve a hole in the disc around 1 AU, preventing the inward migration of giant planets (Alexander & Pascucci 2012). These effects should be investigated in future simulations.

In the disc model of Bitsch et al. (2015a), the α viscosity parameter is set to 0.0054 to calculate the thermal structure. Here we keep the same disc structure, but use a reduced α value for the planet migration, α_mig = 10⁻⁴ as in our previous simulations. We motivate this choice by the theory that the accretion of discs is driven by winds on the surface, which transport the angular momentum. The gas falls in and is accreted efficiently onto the star, while the bulk viscosity in the disc midplane remains low (Suzuki et al. 2016). The value of α_mig is in line with observations of dust settling towards the midplane (Pinte et al. 2016; Dullemond et al. 2018).

The low viscosity used for planet migration, results in a slow inward migration of the growing planets. Larger values of viscosity, onthe other hand, would result in faster inward migration, because the growing planets transition later to the slow type-II migration resulting in more inward migration, which is probably not in line with the structure of the solar system (Bitsch et al. 2019).

In the classical K-damping prescriptions (Lee & Peale 2002), K ∝ 1∕h. This indicates that our simulations that match the eccentricity distribution of the giant planets imply in this simplistic way a different disc structure compared to the simulations that show a strong correlation between inner super-Earths and outer gas giants. This aspect deserves further investigation especially under the aspect that disc masses and thus disc structures can vary a lot from star-to-star (Andrews et al. 2013) and could even invoke differences in the super-Earth systems found by Kepler (Kutra & Wu 2020).

In our simulations, the planets grow by pebble accretion. We use here a model where the planets only accrete the dominating grain size (in the drift limit), without modelling a full-size distribution. However, the majority of the mass of a size distribution of pebbles is in the large grains, which is the sizes the planets accrete. In addition, our pebble flux model is based on an equilibrium between growth and radial drift, where the pebble flux reduces as the accretion rate of the disc reduces in time as well. Reducing the accretion rate all over the disc, as implied by Ṁ disc models, results in a decrease in the gas surface density and with it automatically in a reduction of the solid density (which is bound to the gas surface density through the dust-to-gas ratio). This reduction, however, results in pebble densities in the outer disc that are too low compared to observations (Bitsch et al. 2018a), which is why we increase the pebble flux S_peb by either 2.5 or5.0 in our model, as in our previous works (Bitsch et al. 2019; Izidoro et al. 2019). This increase in the pebble flux is artificial, but the resulting pebble densities are in agreement with the observations. In addition, we do not aim here to study planet population synthesis model, where simulation parameters have to be as self-consistent as possible, but investigate a simpler question regarding the origin of the eccentricities of giant planets.

In our simulations, we start with an initial distribution of planetary embryos between 3 and 18 AU, where the planetary embryos have initially all around 0.01 Earth masses. Planetesimals formed by the streaming instability, however, are only around 100 km in size (Johansen et al. 2015; Simon et al. 2016; Schreiber & Klahr 2018), much smaller than 0.01 Earth masses. At these sizes pebble accretion is also very inefficient (Visser & Ormel 2016; Johansen & Lambrechts 2017), so these planetesimals need to collide first to form bigger objects until pebble accretion can take over and becomes more efficient than planetesimal accretion, especially exterior to 1 AU (Johansen & Bitsch 2019; Voelkel et al. 2020). This growth phase by planetesimals depends crucially on the amount of available planetesimals, so the planetesimal surface density. How and where planetesimals in discs form is still under debate, where some ideas suggest that the water ice line is a favourable position (Armitage et al. 2016; Dra̧żkowska & Alibert 2017; Ida & Guillot 2016), while other simulations are based on the concept that planetesimals form in vortices that can exist all over the disc (Lenz et al. 2019). These different ideas lead to a different initial distribution of planetesimals and thus planetary embryos (Voelkel et al. 2020). Future models of planet formation should include this step self consistently.

In our model, the gas accretion process is modelled in two steps. First we follow a contraction phase of the planetary envelope, following the approach of Piso & Youdin (2014), which is based on Ikoma et al. (2000), where the planet then enters the runaway accretion phase once the planetary envelope becomes as massive as the planetary core. The contraction phase then basically decides the fate of the planet. Planets with slow contracting envelopes remain small, like super-Earths or mini-Neptunes, while planets with fast contracting envelopes can grow to gas giants. This contraction phase depends mostly on the planetary mass itself and the opacity in its envelope. We use here a constant opacity for the envelope following Movshovitz & Podolak (2008). We do not investigate here different envelope opacities, as the goal of our simulations is to generate giant planets and to study their dynamical interactions. However, our simulations show a very large abundance of warm Jupiters (r < 1 AU), which could interfere with the survival of super-Earths. Future simulations with less efficient gas accretion might shine light on this relation.

7.3 Previous simulations

Since the first giant exoplanets on eccentric orbits have been discovered, scattering among a few of these objects was proposed as an explanation for their orbital properties (Ford et al. 2001, 2005; Jurić & Tremaine 2008). Many simulations since then have been undertaken to study how the eccentricity distribution could be originating from scattering events (Raymond et al. 2009a; Sotiriadis et al. 2017).

The simulations by Raymond et al. (2009a) start with three giant planets with an initial separation of 4–5 mutual Hill radii. These small distances between the planets lead very soon to interactions, resulting in scattering events. The resulting eccentricity distribution matched the observations of exoplanets quite nicely. In addition, Raymond et al. (2009a) include the effects of an outer planetesimal disc, which can induce scattering, similar to what is proposed to have happened in the solar system (Tsiganis et al. 2005).

Sotiriadis et al. (2017) improved on this concept a step further, by including the type-II migration phase of the giant planets during the gas disc. In their simulations, the giant planet were also already fully formed. In addition, they kept the innermost of their three giant planet fixed, which allowed the outer planets to migrate closer, resulting in recreating the initial conditions of Raymond et al. (2009a) and thus also the eccentricity distribution of the giant planets.

In our simulations, we include also the growth phase of the planet from a planetary embryo all the way to a gas giant. In addition, our simulations include also self-consistent the formation of inner super-Earths, not included in the previous simulations.

7.4 Resonances of giant planets

The growth of planets by pebble accretion happens inside-out in our simulations (e.g. Fig. 2), where planets in the inner disc reach the pebble isolation mass first. These planets can then start to contract their gaseous envelope and become gas giants. Once these planets reach masses of several 10 Earth masses, they start to open gaps in the disc and eventually migrate in the slower type-II migration regime. Planets starting further out, grow slower and are thus longer in the faster type-I migration regime. This faster type-I migration allows them to catch up to the inner, slower migrating giant planets. This behaviour is a typical outcome of planet migration simulations (Masset & Snellgrove 2001), if the outer planet is smaller than the inner planet and can lead to planets getting trapped into a resonance configuration. This has also been applied to the solar system either as starting conditions for the Nice-model (Tsiganis et al. 2005) or also as an ingredient of the Grand Tack model⁵ (Walsh et al. 2011), where the giant planets then migrate outwards.

However, this migration behaviour does not necessarily need to result in capture in resonance (Pierens et al. 2014), but can also lead to instabilities during the end of the gas phase (e.g. Fig. 3 and Lega et al. 2013). In most of our simulations, if the outer planets are in resonance configuration in the outer disc, these resonances of giant planets are broken when the system becomes dynamically unstable. The instability could also be aided by the growth of the planets in resonance, which can destabilise the system (Matsumoto & Ogihara 2020). In systems that do not undergo any instability, these resonant configurations are maintained (e.g. system 5 in Fig. B.2, where the outer three planets are in 3:2 and 2:1 resonance configuration). However, the majority of the systems in our simulations do not show resonant configuration of the outer giant planets.

7.5 Metallicity – eccentricity correlation

Recent observations have revealed that the eccentricity of cold gas giants is related to the host star metallicity (Dawson & Murray-Clay 2013; Buchhave et al. 2018). This can be explained by the fact that giant planets are more abundant around metal rich stars⁶ (Santos et al. 2004; Fischer & Valenti 2005; Johnson et al. 2010; Narang et al. 2018), which enhances the probability that these giant planets scatter after the gas disc dissipated, resulting in eccentric giant planets. On the other hand, if the host star is less metal rich, giant planet formation is hindered and maybe only one gas giant might form, which will then remain probably on a very circular orbit, due to the lack of partners to scatter with.

In this work, we have tested the outcome of our simulations with two different pebble fluxes and three different number of initial planetary embryos, which can function as a proxy of the host star metallicity. The final eccentricities of the giant planets formed in the low pebble flux simulations and for the same initial number of embryos are generally a bit lower (Figs. 6 and 7) than in the simulations with larger pebble flux (Figs. 11 and 12). In addition, the simulations with less initial embryos show a steeper eccentricity distribution, implying that not enough planets are available to lead to massive scattering events (Fig. 6). However the largest influence on the eccentricities of the giant planets originates from the different damping rates of eccentricity and inclination. In general the trend in our simulations thus seem to confirm the idea that scattering is mostly responsible to explain the different eccentricities of giant planets around stars with different metallicities.

7.6 Hot Jupiters

The first exoplanet found around a main sequence star was a Jupiter type planet orbiting its host star in just a few days (Mayor & Queloz 1995). These planets were named hot Jupiters and even though they are very easy to detect, they should only exist around 1–2% of stars (Mayor et al. 2011). In addition, the host star metallicity, which is a proxy for the amount of building blocks available to form planets, is similar to the host star metallicity of giant planets on eccentric orbits (Buchhave et al. 2018), hinting potentially at a common origin of these planets.

In our simulations, the inner edge of the disc is at 0.1 AU, corresponding roughly to a 10-day orbit. In our simulations only less than 1% of the giant planets end up closer than 0.1 AU. Considering the limitations of our model, this result is actually quite encouraging, as future models that can also include more physics (e.g. tides).

7.7 Heavy element content of giant planets

The large heavy element content of giant planets presents a puzzle for planet formation theories, where a Jupiter mass planet on average should harbour around 60 Earth masses of heavy elements (Thorngren et al. 2016). This result is based on the match between observed planetary masses and their radii with interior models. With pure core accretion this is hard to achieve, because the pebble isolation mass in the inner disc is much lower than that and of the order of 10–20 Earth masses (Bitsch et al. 2018b).

In a recent work by Ginzburg & Chiang (2020), it was proposed that collisions between the giant planets could increase the heavy element fraction of the surviving planet. In their model, the planets need to merge multiple times in order to achieve the large heavy element contents predicted by the observations. However, Ginzburg & Chiang (2020) did not model the growth, migration or scattering of these giant planets directly.

In our simulations we observe an average of 2.5 collisions per system in the case of slow damping, which matches the eccentricity distribution of the giant planets best. In about 35% of our systems do planet experience more than one merger event and only about 30% of those experience more than two merger events⁷. However, in most of the cases these events happen for super-Earths planets with small cores. In addition, not all the planets that underwent merger events survive the dynamical instabilities in the system. This implies that collisions could only account for a very tiny fraction of planets with large heavy element content and we thus deem that the scenario of Ginzburg & Chiang (2020) is probably not the main reason for the heavy element content of the giant planets.

Another mechanism that could explain the heavy element content of giant planets is the evaporation of drifting pebbles, which locally enriches the gas by their volatile contributions released into the case (Booth et al. 2017). This could account for large heavy element contents in the giant planet atmospheres, but needs to be tested within a framework of multi planet formation.

7.8 Implications for protoplanetary disc observations

Recent observations of protoplanetary discs with ALMA have revealed an amazing level of substructures in these discs (Andrews et al. 2018), where basically all of these large discs have rings and gaps in the mm emissions. There are many ideas what could cause these rings and gaps, for instance, ice lines (Zhang et al. 2015), MHD effects (Flock et al. 2015), but also massive planets that generate pressure perturbations in the protoplanetary discs where pebbles accumulate (Paardekooper & Mellema 2006; Pinilla et al. 2012; Zhang et al. 2018).

Recent observations of the gas velocity dispersions in protoplanetary discs have revealed that some discs could indeed harbour giant planets (Teague et al. 2018; Pinte et al. 2018), which cause rings and gaps. However, the real cause of these rings and gaps is still under debate.

In our simulations, the growing planets generate pressure perturbations, where dust could accumulate and form the rings we observe with ALMA. However, in the simulations with slow damping (low K) and S_peb = 5.0, the planets already interact gravitationally and scatter during the gas disc phase. This would destroy the nearly axisymmetric rings and gaps observed by ALMA. In particular multiple gas giants on eccentric orbits can result in a very chaotic gas distribution (Lega et al. 2013), which might not be in line with the observations of clean rings and gaps in protoplanetary discs. On the other hand, the structures generated by multiple giant planets are in line with the observation of the PDS70 system (Keppler et al. 2018; Bae et al. 2019).

These recent ALMA observations now give another constraint on planet formation theories. Not only do the models have to match the exoplanet observations, but they should also be in line with the observations of their birth phase, namely their natal protoplanetary discs. We think that this avenue needs to be explored in much more detail in the future to constrain planet formation models.

8 Summary

In this work we have combined an N-body framework with pebble and gas accretion as well as planet migration in an evolving protoplanetary discs following our previous works (Bitsch et al. 2019). We investigated the influence of different gas damping rates of eccentricity and inclination of giant planets as well as the influence of a different number of initial planetary embryos to study the eccentricity distribution of giant planets. In addition, our simulations formed self consistently systems with inner super-Earths and outer gas giants, which was not included in past simulations aimed to study the eccentricity distribution of giant planets (Jurić & Tremaine 2008; Raymond et al. 2009a; Sotiriadis et al. 2017). We show some of the systems formed in our simulations in Appendix B.

Our simulations show that fast damping (large K) results in planetary systems that do not undergo large scattering events during the gas disc phase and afterwards. The resulting eccentricities are too low to be close to the observations of giant planets. This result seems to be independent of how fast the planets actually grow in the disc.

Slow damping (low K) of eccentricity and inclination, on the other hand, results in scattering events already during the gas disc phase. In the end, these simulations reproduce the eccentricity distribution of the giant exoplanets, where K needs to be between 5 and 50. These results are independent of how fast the planets grow in the disc as well, implying that the damping rates are more important in determining the fate of the planetary system. Future hydrodynamical simulations with multiple planets are needed to determine conclusively how damping of eccentricity and inclination evolves for systems of multiple giant planets.

In addition, we find only a small dependency on the number of giant planets formed in our simulation on the initial number of planetary embryos. If many embryos are implanted within the disc, the gravitational interactions are already initially quite strong, scattering the majority of objects and thus preventing their growth. In addition, our simulations suggest that the giant planets on average should not be single, but should come in multiples.

The simulations that match the eccentricity distribution of the giant planets best, undergo large scattering events. As a consequence of these scattering events, the inner super-Earth systems are destroyed and only a very small fraction of the systems formed in our simulations harbour inner super-Earths and outer gas giants. Schlecker et al. (2020) studied the super-Earth - cold Jupiter relation as well, but used instead 300m planetesimals for the growth of their planetary embryos, not in agreement with the evidence in the solar system (Bottke et al. 2005; Morbidelli et al. 2009; Singer et al. 2019; Stern et al. 2019). However, they recover a similar result regarding the eccentricity of outer gas giants and surviving inner super-Earths.

On the first look this seems to contradict the observations using the Kepler data (Zhu & Wu 2018; Bryan et al. 2019), while RV observations seem to agree with our simulation results (Barbato et al. 2018). However, these studies did not take the eccentricity of the giant planets into account. Our simulations predict that systems with lower mass giant planets on nearly circular orbits are more likely to harbour systems of inner super Earths. On the other hand, systems harbouring massive giant planets on eccentric orbits should not harbour any inner super-Earths. More detailed statistics of the super-Earth–cold Jupiter relation including the eccentricity distribution of the giant planets could verify if our planet formation approach is correct and constrain future theories.

Acknowledgements

B.B. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. Andre Izidoro acknowledges NASA grant 80NSSC18K0828 to Rajdeep Dasgupta, during preparation and submission of the work. A.I. also thanks financial support from FAPESP via grants 16/12686-2 and 16/19556-7, and CNPq via grant 313998/2018-3. We thank the referee John Chambers, whose comments helped a lot to improve the manuscript.

Appendix A Evolution of individual planetary systems

We show here additional evolution of systems with K = 500 and 30 initial planetary embryos as well as simulations with 15 or 60 initial planetary embryos and K = 50. All the here presented simulations feature S_peb = 2.5.

Fig. A.1 Evolution of a system with a damping factor of K = 500. The plots and lines have the same meaning as in Fig. 2. In the end three gas giants survive with a significant eccentricity, originating from a scattering event about 15 Myr after gas disc dispersal.

We show the evolution of a planetary system using K = 500 in Fig. A.1. The evolution of the planetary system is very similar as in the K = 50 case during the early stages of the gas disc lifetime. However, due to the larger K value, the damping is still efficient at the end of the gas disc lifetime, resulting in a multi planetary system with inner super-Earths and outer gas giants emerging from the gas disc. The eccentricities of these planets is of the order of a few percent.

But at about 20 Myr, the planetary system undergoes an instability event, including a collision between two planets, resulting in a planet of about 2 Jupiter masses. The smaller planets of the system have been ejected and only three gas giants on eccentric orbits (with an average eccentricity of 0.2–0.3) survive.

In Fig. A.2, we show a simulation using K = 50, but in contrast to Fig. 2 with only 15 initial planetary embryos. The overall evolution during the gas disc phase is again very similar to the simulation with 30 embryos, except that even the small planetary embryos start to grow to at least ≈0.1 Earth masses. This is caused by the initially larger separation between the bodies, which prevents initially planet-planet interactions and thus keeps the eccentricities and inclinations low during this phase, allowing all planetary embryos to accrete some pebbles. However, as a few dominating bodies emerge, the eccentricities and inclinations start to increase and the small bodies stop growing.

Similar to the simulation with 30 bodies (Fig. 2), the damping of the gas disc reduces towards the end of the disc lifetime and the planetary systems undergoes an instability event, which removes the smaller bodies from the disc, where only two gas giants on highly eccentric orbits remain. However, these two gas giants scatter again after 20 Myr and only one gas giant on a highly eccentric orbit survives.

Finally,we show in Fig. A.3 a simulation with initially 60 planetaryembryos and again with K = 50. The initiallylarge number of planetary embryos increases the interactions between them compared to the previous simulations. This results in an increase of eccentricity and inclination of the small bodies as soon as the first bodies start to reach masses above Earthmass. As a consequence, as before, only a few dominating bodies emerge. As in the previous simulation, towards the end of the disc lifetime the system becomes dynamically unstable and only two gas giants on highly eccentric orbits survive.

Even though the simulations start with a different number of initial planetary embryos, the resulting planetary systems are very similar for simulations with the same damping values. The reason for that is that only a few dominating bodies emerge in each simulation, similar to the N-body simulations with pebble accretion from Levison et al. (2015), who, however, did not include planet migration. If many bodies are initially present, many of them acquire an eccentricity of a few percent that prevents the planets to accrete pebbles efficiently and thus stopping their growth (Johansen et al. 2015). As a result only a few bodies start to grow efficiently. On the other hand, if only initially a small number of embryos is present, then obviously only a few bodies can grow. This is not only evident in the here presented simulations (Figs. A.2, 2 and A.3), but is true for the full sets of our simulations.

Fig. A.2 Evolution of a system with a damping factor of K = 50 and 15 initial embryos. The plots and lines have the same meaning as in Fig. 2. After the scattering event shortly about 15 Myr after the gas disc phase, only one gas giant ona highly eccentric orbit survives.

Fig. A.3 Evolution of a system with a damping factor of K = 50 and 60 initial embryos. The plots and lines have the same meaning as in Fig. 2. After the scattering event shortly after the end of the gas disc phase, two gas giants on eccentric orbits survive.

However, if K is large and thus damping of eccentricity and inclination is efficient, the initially larger number of planetary embryos can lead to a larger number of surviving planets, because the planets remain on nearly circular orbits during the gas disc phase (Figs. A.1 and 3), increasing the stability of the systems also after the gas disc phase. In this case, accretion is very efficient and systems with many planets can emerge.

Appendix B Structures of planetary systems

We show here the structure of planetary systems formed in discs with K = 5 and K = 50 in the case of 30 initial embryos and S_peb = 5.0 at 100 Myr in Figs. B.1 and B.2. The size of the circle is proportional to the total planetary mass (green) by the 3rd root and to the mass of the planetary core (black) also by the 3rd root. The black arrows indicate the aphelion and perihelion positions of the planet calculated through r_P ± e × r_P.

Clearly, a larger K value, whichimplies more efficient damping of eccentricity and inclination during the gas disc phase, results in planetary systems with more giant planets. In addition, it is also clear that the eccentricities of the giant planets are larger for K = 5 compared to K = 50.

Fig. B.1 Final configurations after 100 Myr of integration of all our simulations with Speb = 5.0 and K = 5. The size of the circle is proportional to the total planetary mass (green) by the 3rd root and to the mass of the planetary core (black) also by the 3rd root. The black arrows indicate the aphelion and perihelion positions of the planet calculated through rP ± e × rP.

Fig. B.2 Same as Fig. B.1, except that K = 50 was used.

References

All Tables

All Figures

Fig. 1 Change of the eccentricity of a one Jupiter mass planet as function of time for the different damping parameters. The K damping values are derived using the classical type-II migration timescales. The flattening of the damping at the end of the disc lifetime is due to the exponential decay of the gas disc surface in the last 100 kyr, which reduces the gas surface density and thus the damping effects.

In the text

Fig. 2 Evolution of a system with a damping factor of K = 50. Semi major axis (top left), planetary mass (top right), eccentricity (bottom left) and inclination (bottom right) of 30 planetary embryos as function of time. The gas disc lifetime is 3 Myr after injection of the planetary embryos, where the end of the gas disc lifetime is marked by the vertical black line. The grey lines represent either small mass bodies, or bodies that are ejected during the lifetime of the gas disc. The black lines represent massive planets that are scattered away after the gas disc phase. The coloured lines represent the two surviving planets.

In the text

	Fig. 3 Evolution of a system of 30 initial planetary embryos with a damping factor of K = 5000. The plots and lines have the same meaning as in Fig. 2. In this case, no giant instability after the gas disc phase happened within 100 Myr and a system with inner super-Earths and several outer gas giants survived.
In the text

Fig. 4 Distances of adjacent planetary pairs of all planets in the simulations starting with 60 planetary embryos (including the small bodies) after 3 Myr (end of the gas disc phase, top) and 100 Myr (dashed, bottom). In addition, we show the distances after 3 Myr only for planets larger than 0.5 Jupiter masses (long-dashed, middle panel), to indicate how the stability of the system is influenced by the larger bodies. The black vertical line marks a separation of mutual Hill radii, at which two planets on circular orbits avoid mutual encounters for all time (Gladman 1993). We only show the separation in units of mutual Hill radii for simulations with initially 60 bodies. The trends are very similar for simulations with 15 and 30 initial embryos.

In the text

Fig. 5 Period ratios of adjacent planetary pairs of all planets in the simulations (including the small bodies) after 3 Myr (end of the gas disc phase) and 100 Myr (dashed). In addition, we show the period distribution including only planets with masses larger than 0.5 Jupiter masses after 100 Myr (long-dashed lines). We do not observe a significant pile-up of planets close to mean-motion resonances, in contrast to the super-Earth simulations of Izidoro et al. (2019), which we attribute to the slower migration speed in the here used simulations, which prevents the formation of resonant chains at the inner edge of the protoplanetary discs. Nevertheless, the systems become unstable and the instabilities result in a wide spacing of the planets after 100 Myr.

In the text

	Fig. 6 Eccentricity distribution of the simulated planets within a distance of 5.2 AU from their host star with masses of at least 0.5 Jupiter masses as well as of RV observations of giant planets with masses between 0.5 and 1.25 Jupiter masses exterior to 0.1 AU. We show in colour the different damping rates and systems with initially 15, 30 and 60 embryos (from top to bottom).
In the text

Fig. 7 Eccentricity as a function of semi-major axis of simulated planetary systems with different damping rates anddifferent number of initial embryos (15, 30, 60 from top to bottom). We only show the planets with masses of at least 0.5 Jupiter masses (for simulations and observations) and only for the K values that match the observations of the eccentricity distribution best (Fig. 6). The horizontal lines depict the mean eccentricity within each radial bin. The black symbols show the data from exoplanet.eu with planetary masses ranging from 0.5 to 1.25 Jupiter masses, which is roughly the masses the planets in our simulations reach.

In the text

Fig. 8 Evolution of a system with a damping factor of K = 50 where the pebble flux is twice as large as in Fig. 2. The lines have the same meaning as in Fig. 2. Due to the massive planets, an instability already occurs during the gas phase of the disc, increasing the eccentricity of the planets. However, during the remaining gas phase of 1.5 Myr after the instability, the eccentricities of the planets get damped efficiently and the two inner planets have eccentricities of around 0.1.

In the text

	Fig. 9 Distances (top) and period ratios (bottom) of planets formed in our simulations with Speb = 5.0 and different damping values. For this plot we cut the period ratio at 5.
In the text

	Fig. 10 Mass distribution of the planets formed with masses larger than 0.5 Jupiter masses in our simulations and from RV observations, where we limit the observations to 5.0 Jupiter masses, as this is roughly the largest masses we reach in our simulations.
In the text

	Fig. 11 Eccentricity distribution of the planets formed in our simulations with faster growth (Speb = 5.0) and the observations with planetary masses from 0.5 to 5 Jupiter masses.
In the text

Fig. 12 Eccentricity–semi-major axis distribution of the observed giant planets (black triangles) with masses of 0.5 < MP < 5.0MJ and our simulations with different K-damping values (coloured circles) after 100 Myr of evolution. The horizontal lines show the mean eccentricity within each orbital distance bin. While our simulations with low K reproduce the e − a distribution exterior to 1 AU, the simulations show too small eccentricities for the planets interior to 1 AU. On the other hand, this might also be related to the limited number of giant planets within this region in our simulations. For planets in the order of Jupiter mass, the low K values match quite nicely (Fig. 7), due to the slightly larger number of planets within that semi-major axis bin in our simulations.

In the text

Fig. 13 Eccentricity distribution of the sample of exoplanet.eu with masses larger than 0.5 Jupiter masses where the planets have a minimal distance to their host star of 0.1 AU (solid black line) or 1.0 AU (dashed black line). In red we show the eccentricity distribution of the planets observed in Barbato et al. (2018), while we show eccentricity of the giant planets used in the study of Zhu & Wu (2018) in blue. In the dashed blue solid line we also show the sample of Zhu & Wu (2018), but exclude two planetary systems where the companions are more massive than 5 Jupiter masses, which could inherit their eccentricity from planet–disc interactions.

In the text

	Fig. A.1 Evolution of a system with a damping factor of K = 500. The plots and lines have the same meaning as in Fig. 2. In the end three gas giants survive with a significant eccentricity, originating from a scattering event about 15 Myr after gas disc dispersal.
In the text

	Fig. A.2 Evolution of a system with a damping factor of K = 50 and 15 initial embryos. The plots and lines have the same meaning as in Fig. 2. After the scattering event shortly about 15 Myr after the gas disc phase, only one gas giant ona highly eccentric orbit survives.
In the text

	Fig. A.3 Evolution of a system with a damping factor of K = 50 and 60 initial embryos. The plots and lines have the same meaning as in Fig. 2. After the scattering event shortly after the end of the gas disc phase, two gas giants on eccentric orbits survive.
In the text

	Fig. B.1 Final configurations after 100 Myr of integration of all our simulations with Speb = 5.0 and K = 5. The size of the circle is proportional to the total planetary mass (green) by the 3rd root and to the mass of the planetary core (black) also by the 3rd root. The black arrows indicate the aphelion and perihelion positions of the planet calculated through rP ± e × rP.
In the text

	Fig. B.2 Same as Fig. B.1, except that K = 50 was used.
In the text