Open Access
Issue
A&A
Volume 656, December 2021
Article Number A69
Number of page(s) 44
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202038553
Published online 06 December 2021

© A. Emsenhuber et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Since the discovery of the first exoplanet detected around a main sequence star (Mayor & Queloz 1995), the number of known exoplanets has greatly increased. These planets span a wide range of masses and sizes, and they were detected using various techniques, such as radial velocity, transit, direct imaging, and microlensing. Despite all these observational constraints, the exact formation pathways are not yet certain. To highlight this, we first discuss possible formation mechanisms for different planet kinds.

Giant planets have been found orbiting their host star over a wide range of periods. Some are of the order of days or tens of days, which is well within the orbit of Mercury (Mayor et al. 2011; Fabrycky et al. 2014); others were detected at large separations using the direct imaging technique (Marois et al. 2008; Lagrange et al. 2010; Rameau et al. 2013; Macintosh et al. 2015; Chauvin et al. 2017; Keppler et al. 2018).

Most giant planets are thought to form via the core accretion mechanism as gravitational instability (Boss 1997, 2003) is found to work only at large separation (several tens of astronomical units, Rafikov 2005; Schib et al. 2021), though the clumps could migrate after formation (Nayakshin 2010), and for bodies above about 5 M (Schlaufman 2018) or even the deuterium-burning limit (Kratter et al. 2010). On the other extreme, for very close-in giant planets, core accretion (Perri & Cameron 1974; Mizuno 1980) requires for in situ formation a very strong pile-up of solids. While this has been proposed (Boley et al. 2016; Batygin et al. 2016; Bailey & Batygin 2018), the possibility remains heavily debated. A scenario where these planets formed further out and were subsequently moved to their final location (Lin et al. 1996) is usually considered more likely.

In the standard view, giant planets form from embryos located beyond the ice line, where solids are abundant owing to the volatiles being present in the solid form. This allows the embryo to form rapidly enough before the dispersal of the gas disc, which occurs in a time frame of several million years (Haisch et al. 2001; Fedele et al. 2010; Richert et al. 2018). Embryos initially accrete solids and a small quantity of gas. The further growth results in the accretion of gas, which is governed by the ability of the planet to radiate away the accretion energy (Pollack et al. 1996; Lee & Chiang 2015). The cooling process becomes more efficient as the mass increase, so that when the planet reaches a mass of several times that of the Earth, the amount of solids and gas are equal (the critical mass, Stevenson 1982). Once the accretion rate becomes greater than what the disc is able to supply, the envelope can no longer remain in equilibrium with the surrounding nebula and it contracts.

This process is further complicated by the implications of planetary migration (Baruteau et al. 2014, 2016). The final mass and locationof the planet depends thus on the interplay between growth and migration, not to mention the interactions with the other planets forming in the same system.

Observations show that the giant planets are divided into two sub-groups depending on the host-star metallicity (Dawson & Murray-Clay 2013; Buchhave et al. 2018). Hot-Jupiters around metal-poor stars exhibit lower stellar obliquity and eccentricity than the ones around metal-rich stars. The usual concept of inward migration due to interaction with the gas disc (Goldreich & Tremaine 1979; Ward 1997; Tanaka et al. 2002) cannot account for the obliquity of these planets, which more likely were brought there by few-body interactions combined with tidal circularisation (Dawson & Johnson 2018).

For the distant giant planets, core accretion is still favoured (Wagner et al. 2019). A possible formation pathway for some of these distant planets is accretion in the inner region of the disc followed by close encounters and scattering (Marleau et al. 2019a). This pathway is supported by evidence that it is able to reproduce the distribution of eccentricities of giant planets (Chatterjee et al. 2008; Jurić & Tremaine 2008; Ford & Rasio 2008; Raymond et al. 2010; Sotiriadis et al. 2017), and that most giant planet-harbouring system are multiple (Knutson et al. 2014; Bryan et al. 2016; Wagner et al. 2019).

Exoplanets include planets unknown in the solar system, those between the Earth and Neptune (Mayor et al. 2011; Youdin 2011; Howard et al. 2012). The density of these planets vary more than one order of magnitude (Hatzes & Rauer 2015; Otegi et al. 2020).

Sub-Neptunes exhibit a low bulk density, indicating the presence of a gaseous envelope (Weiss & Marcy 2014; Rogers 2015). This impliesthat they mostly formed in a time scale comparable with the lifetime of the protoplanetary disc. However, whether they formed early (in the same way as the core of giant planets) or towards the end of the disc (Lee & Chiang 2016) is not yet settled.

Super-Earths on the other hand are compatible with being gas-free. They are not constrained by the lifetime of the protoplanetary disc and can form over longer periods of time (Lambrechts et al. 2019; Ogihara et al. 2018). These could also have had a gaseous envelope in the past that was removed by, for instance, atmospheric escape (e.g. Jin et al. 2014) or giant impacts (Schlichting & Mukhopadhyay 2018).

Multi-planetary systems provide additional information. Many super-Earth systems have similar mass and spacing (Millholland et al. 2017; Weiss et al. 2018), though this is debated (Zhu 2020; Weiss & Petigura 2020). However, most of the planet pairs are out of mean-motion resonances (Fabrycky et al. 2014). The low number of planets in mean-motion resonances (MMR) may be surprising, as gas-driven migration is efficient at capturing the planets in MMRs. But it is possible for the resonances to be broken during the retreat of the gas disc (Liu et al. 2017) or after the dispersal by dynamical instabilities (Inamdar & Schlichting 2016; Izidoro et al. 2017, 2021). The mutual inclinations remain relatively low (Lissauer et al. 2011; Fang & Margot 2012) and they exhibit low-to-moderate eccentricities (Xie et al. 2016; Mills et al. 2019).

Many models have been developed to capture the above-mentioned effects. We may cite Pollack et al. (1996), Ida & Lin (2004a,b), Alibert et al. (2004, 2005a), Miguel et al. (2011), Coleman & Nelson (2014), Cridland et al. (2016, 2017), Liu et al. (2017), Ormel et al. (2017), Ronco et al. (2017), Ndugu et al. (2018), Chambers (2018), Alessi & Pudritz (2018), Bitsch et al. (2019), Johansen et al. (2019), Booth & Ilee (2019), or Alessi et al. (2020a,b), to mention only a few. To capture all of the above effects, models of planetary formation must include many physical processes occurring during the formation of the systems, which lead to feedbacks and non-linearities. Then, comparing the outcomes of global end-to-end models with observations is a key step for the understanding of the origins and evolution of planets systems.

As the model has many parameters, a large number of planets with different properties are required to constrain their possible values. The model must then be able to predict all the necessary observable quantities for the different observational techniques, not only masses and distances, but also radii (for transits), luminosities, magnitudes (for direct imaging), and evaporation rates. To leverage the enormous amount of statistical observational data on exoplanets, the models should also be able to make quantitative predictions which can be compared statistically with the actual planetary population. For this, planetary population synthesis (Ida & Lin 2004a; Mordasini et al. 2009a) is a frequently used approach.

In this work, we introduce a strongly improved and extended version the Bern global end-to-end model of planetary formation and evolution for multi-planetary systems. The model combines the work of Alibert et al. (2013, hereafter A13) (inclusion of N-body interactions) and the internal structure calculations and long-term thermodynamical evolution model of Mordasini et al. (2012c,b). Here, we track the planets with full N-body interactions, in contrast to Ida & Lin (2010) for instance, who introduced a semi-analytical approach, an improvement over previous works, such as Ida & Lin (2004a), to follow planet-planet interactions.The model follows the formation of many embryos, as it usually obtained from the end stage of the runaway growth of solids (Kokubo & Ida 1998), so that both terrestrial as well as giant planetary systems can be obtained.

The structure of this work is as follows: in the first part, we describe our global model. In Sect. 2, we introduce the new version of our model with a general overview of its conception, along with its relationship to previous work. Detailed description are provided in Sect. 3 for the stellar and nebular components, in Sect. 4 for the planets, and in Sect. 5 for the migration and dynamical evolution. In a second part we perform some tests for different kind of planets and show possible resulting systems. In Sect. 6, we aim at reproducing formation of terrestrial planets with our improved model to determine whether it is applicable to kind of planets. The formation of giant planets and the implications for Jupiter are discussed in Sect. 7. Finally, in Sect. 8, we apply the presented model to specific systems to assess the interaction between the different mechanisms occurring during the formation and evolution of planetary systems.

This work is the first of a series of several. In a companion paper, Emsenhuber et al. (2021, hereafter referred to as Paper II), we will use this model to compute synthetic populations of planetary systems and perform statistical analysis. In subsequent articles, we will perform more detailed comparison with observations, and analyse various parameters that we have in the present model.

2 The Bern model

2.1 History

The model presented in this work, the Generation III Bern model, combines the formation and evolution stage of planetary system. It is based on many contributions in the field that aim to study different aspects of the physics of planetary formation and evolution. We thus start by a short history of the series of model, and its different branches that we couple together in this work. A graphical sketch of the different generations of the Bern model is provided in Fig. 1.

The original model was introduced in Alibert et al. (2004, 2005a) for individual planets, then used in Mordasini et al. (2009a,b) for entire planetary populations. We refer to it as Generation I, which computed the formation on a single planet until the gas disc disperses. The model subsequently diverged into two different branches: one with the aim to follow the long-term evolution of the formed planet (Generation Ib; Mordasini et al. 2012c,b) while the other obtained the ability to form multi-planetary systems with an improved description of the planetesimals disc (Generation II; Alibert et al. 2013; Fortier et al. 2013). In this work we bring these two variants of the model back together so that we can follow the formation and the long-term evolution of multi-planetary systems. At the sametime, we extend the model with new elements, which are shown in italic on Fig. 1.

Previous versions of the model have been extensively described in referenced papers (see also Benz et al. 2014; Mordasini et al. 2015; Mordasini 2018 for reviews and the interactions between the different mechanisms involved in planet population syntheses). We nevertheless describe this new version in the remainder of this section.

2.2 General description

We base ourstudy on the Bern model of planetary formation and evolution. This global model self-consistently computes the evolution of the gas and planetesimals discs, the accretion of gas and solids by the protoplanets, their internal and atmospheric structure, as well as interactions between the protoplanets and between the gas disc and the protoplanets. We provide a diagram of the main components of the overall model as well as the most important exchanged quantities in Fig. 2.

In our coupled formation and evolution model, we first model the planets’ main formation phase for a fixed time interval (set to 20 Myr, see the related discussion in Sect. 6 regarding the impact of this specific choice). Afterwards, in the evolutionary phase, we follow the thermodynamical evolution of each planet individually to 10 Gyr.

2.2.1 Formation phase

During the formation stage (0–20 Myr), the model follows the evolution of a gaseous protoplanetary disc and the dynamical state of planetesimals (Sect. 3). These serve as sources for the accretion of the protoplanets (Sect. 4). The lifetimeof the gas disc is shorter than the simulated formation stage, so that solids accretion in a gas-free environment can also take place. The gas disc also leads to planetary migration, and interactions (scattering, collisions) between the concurrently growing proptoplanets are tracked via a N-body integrator (Sect. 5).

The formation of planets is based on the core accretion paradigm (Mizuno 1980; Pollack et al. 1996): first, a solid core is formed, and once it becomes massive enough, it starts to bind a significant H/He envelope. Core growth results from the accretion of planetesimals. Gas accretion is initially governed by the ability of the planet to radiate away the energy released by the accretion of both solids and gas. Once the gas accretion rate of the envelope exceeds the limit from the disc, the envelope can no longer maintain equilibrium with the disc; hence it subsequently contracts and passes into the detached phase (Bodenheimer et al. 2000).

2.2.2 Evolution phase

The long-term evolution of the planets (20 Myr–10 Gyr) is calculated by solving, like already in the formation phase, the standard spherically symmetric internal structure equations, but with different boundary conditions, and taking into account different physical effects like atmospheric escape, or radius inflation. In this phase, the planets evolve individually; N-body interactions and the accretion of planetesimals are no more considered. The orbits and masses of the planets may however still evolve because of effects like tides and atmospheric escape.

As described in Mordasini et al. (2012c), the coupling between the formation and evolution phases is made self-consistently, that is both the compositional information as well as the gravothermal heat content given by the formation model are given to the evolution model as initial conditions.

Regarding the temporal evolution, we now also take the thermal energy content of the planet’s core into account for a planet’s luminosity, as described in Linder et al. (2019). This is important for core-dominated low-mass planets (e.g. Lopez & Fortney 2014). As in previous calculations, the other gravothermal energy sources are the cooling and contraction of the H/He envelope, the contraction of the core, and radiogenic heating due to the presence of long-lived radionuclides in the core.

2.2.3 Envelope structure

The calculation of the internal structure of all planets (Sect. 4) during their entire formation and evolution is a crucial aspect of the Bern Model, as visible from its central position in Fig. 2. It not only yields the planetary gas accretion rate in the attached phase but is also key for the accretion of planetesimals via the drag enhanced capture radius. It also yields the radii and luminosities that on one hand enter multiple other sub-modules, and on the other hand are key observable quantities. The internal structure model assumes that planets have an onion-like spherically symmetric structure with an iron core, a silicate mantle, and depending of a planet’s accretion history, a water ice layer and a gaseous envelope made of pure H/He. In contrast to earlier syntheses predicting planetary radii (Mordasini et al. 2012c, 2014), we now use self-consistently the iron mass fraction as given by the disc compositional model (according to Thiabaud et al. 2014, Sect. 3.3.3), instead of assuming a fixed 2:1 silicate:iron mass ratio.

Physical effects that are included in the model besides the usual cooling and contraction are XUV-driven atmospheric escape (Jin et al. 2014), D-burning (Mollière & Mordasini 2012), Roche-lobe overflow, and bloating of the close-in planets (Sarkis et al. 2021).

Compared to some other 1D internal structures models in the literature (e.g. Vazan et al. 2013; Venturini et al. 2016; Valletta & Helled 2020), the model is simplified first by assuming that the gaseous envelope consists of pure H/He, while accreted solids sink to the core. In this sense, the model is similar to the original Pollack et al. (1996) model. We neglect thus the consequences of heavy element enrichment and compositional gradients in the envelope. These effects will be added in future work. One should note that also other modern models make use of the simplification of pure H/He envelopes (D’Angelo et al. 2021). Including enrichment would generally speed up gas giant formation (Venturini et al. 2016). Second, the effects of hydrodynamic flows affecting the (upper) envelope structure and cooling behaviour are currently also neglected (Ali-Dib et al. 2020; Moldenhauer et al. 2021).

On the other hand, our internal structure model allows to model the entire ‘life’ of planets from t = 0 to 10 Gyr, modelling and coupling self-consistently all phases (attached, detached, evolutionary), for both the gaseous envelope and the solid core. Importantly, the model is capable of calculating the internal structure and temporal evolution of planets ranging in mass from 10−2 M to the lithium-burning limit (about 63 Jovian masses, Burrows et al. 2001). It includes besides the standard aspects (accretion, cooling, contraction) also atmospheric escape, bloating, Roche-lobe overflow, and deuterium burning. In particular, this makes it possible to model planets that reside very close to their host star. This quite unique general applicability to very different planet types reflects the needs arising from a population synthesis calculation.

As shown in Fig. 2, atmospheric escape is only included in the evolution phase starting at 20 Myr. In reality, it would start immediately once the gas disc has dissipated and the planets start to ‘see’ the stellar XUV irradiation. This could lead to a certain under-estimation of atmospheric escape. The consequences should, however, be small, since escape continues to be important for at least the first 100 Myr when stars are in the saturated phase of high XUV emission, and not only for the first 20 Myr. The effect that atmospheric escape can destabilise resonant chains for sufficiently high mass loss (Matsumoto & Ogihara 2020) is thus not included. On the other hand, we include during the entire formation phase (also after gas disc dissipation) the accretion of planetesimals, which also changes planet masses.

In the following sections, we describe in detail all the sub-modules visible in Fig. 2.

thumbnail Fig. 1

Overview of the physical mechanisms included in the Bern model. At the top, the processes and base assumptions made in all model generations are given. The four boxes below show the four model generations with the main paper introducing each generation. The further points are: (1) number of initial embryos per disc, N-body integrator type, initial embryo mass, (2) phases simulated, (3) planetesimal accretion mode and planetesimal size, (4) phases withcalculation of the planets’ internal structure, (5) disc model characteristics, (6) orbital migration: type I, type II, transition criterion from type I to II (here thermal refers to a criterion only with the ratio between the Hill radius and the scale height of the disc, while ‘thermal and viscous’ refers to the full criterion of Crida et al. 2006), (7) disc-limited gas accretion rate, (8) later additions and improvements, (9) additional output relative to older generation, (10) population synthesis publications using this generation. In the bottom right panel, text in italic indicates new elements.

thumbnail Fig. 2

Sub-modules and most important exchanged quantities of the Generation III Bern model. The colours denote the stages at which processes are considered. Blue indicates processes active in the formation stage, but only before the dispersal of the gas disc. Green processes are considered during the entire formation stage, even after the dispersal of the gas disc. Processes in red are only considered during the evolution stage. The processes in black are included in all stages.

3 Star and protoplanetary disc

3.1 Stellar model

Instead of assuming a fixed 1 L stellar luminosity for a 1 M star as in previous model generations, stellar evolution is now considered by incorporating the stellar evolution tracks from Baraffe et al. (2015). These provide the radius R, luminosity L and temperature T for a given stellar mass M at any moment. Stellar temperature and radius are used for the outer boundary conditions of the gas disc; stellar radius is also used in the N-body integrator to detect collisions with the star and to calculate the stellar tidal migration. Finally, the stellar irradiation enters into the calculation of the outer (atmospheric) temperature (at τ = 2∕3) of the planets’ interior structure as described in Mordasini et al. (2012c) and radius bloating (Sect. 4.2.2).

3.2 Gas disc

The protoplanetary gas disc is modelled with a 1D radial axisymetric structure. The evolution is given by solving the viscous diffusion equation as function of the time t and orbital distance r (Lüst 1952; Lynden-Bell & Pringle 1974), (1)

where is the surface density of gas, and and are the sink terms related to photo-evaporation (Sect. 3.2.2) and accretion by the planets respectively. The viscosity is parametrised, following Shakura & Sunyaev (1973), with (2)

This equation is solved on a grid spaced regularly in log with 3400 points that extends from the inner location of the disc rin (an initial condition) to rmax = 1000 au. At these two locations, the surface density is fixed to zero.

3.2.1 Vertical structure

The disc’s vertical structure is computed at each step of the evolution following the approach of Nakamoto & Nakagawa (1994). This change is necessary to accommodate the new stellar model with variable quantities. With this approach, the link between the outer and midplane temperatures is given by (3)

with Tmid the disc mid-plane temperature, Ts the temperature due to irradiation (see below), σSB the Stefan-Boltzmann constant, τR and τP are the Rosseland and Planck mean optical depths respectively, and Ė is the viscous dissipation rate. This formula yields the mid-plane temperature both in the optically-thick (the term with τR) and optically-thin (the term with τP) regimes. The Rosseland optical depth is given by τR = κdisc(ρmid, Tmid)Σ where is the central density, H = cs∕Ω the disc’s vertical scale height, the isothermal sound speed, μ = 2.24 the mean molecular weight of the gas, and mH the mass of an hydrogen atom. The opacity κdisc is given by the maximum of the opacities computed according to Bell & Lin (1994) (which accounts for micrometre size with a fixed interstellar dust-to-gas ratio of 1%, independently of the dust-to-gas ratio chosen for the solids disc) and Freedman et al. (2014) (which givesmolecular opacities for a grain-free gas). For the Planck optical depth, we follow further Nakamoto & Nakagawa (1994) and set τP = 2.4τR.

It is clear that this treatment of the opacities is simplified: in reality, the evolution of the dust via coagulation, fragmentation, and drift influences via the resulting grain opacity the thermal and density structure of the disc. This structure in turn feeds back onto the dust evolution, meaning that the processes must be treated together in a self-consistently coupled way (Gorti et al. 2015; Savvidou et al. 2020).

Such a more realistic coupled model affects for example the disc lifetime, the local dust-to-gas ratio (Gorti et al. 2015), or – in the context of planet formation – the locations of the outward migration zones (Savvidou et al. 2020, see Sect. 5.1.3). They also show that the ratio of Planck and Rosseland opacity is in reality not simply a constant as currently assumed. In Voelkel et al. (2020) we have recently coupled the Birnstiel et al. (2012) dust-pebble evolution model to the Bern Model. Based on this, future version of the Bern Model will include also a more physically realistic grain opacity and therefore disc structure model. This will in particular also include the dependency of the disc opacity on the stellar metallicity, which is currently not taken into account.

In equilibrium, the radiative flux is identical as the viscous dissipation rate, which is given by (4)

with Ω being the Keplerian angular frequency at distance r from the star. The second equality holds only if purely the mass of the central star is accounted for in the Keplerian frequency, that is , being the gravitational constant. The self-gravity of the disc has been neglected.

The disc’s outer temperature due to irradiation is given by (5)

following Hueso & Guillot (2005), but also accounting for the direct irradiation through the disc’s mid plane. The first term inside the bracket is the irradiation of the star onto a flat disc. The second term in the square brackets accounts for the flaring of the disc at large separation. In our case, we do not compute this factor explicitly and instead adopt ln Hlnr = 9∕7 (Chiang & Goldreich 1997).

The Tirr term accounts for the direct irradiation through the disc midplane. It is computed as (6)

which is the black-body equilibrium temperature accounting for the optical depth through the midplane of the disc τmid = ∫ ρmidκ(ρmid, Tmid)dr. This contributionis usually important only at the very end of the disc lifetime while it clears; otherwise, the optical depth confines the contribution to the very innermost region. However, taking this contribution in account is necessary to provide a smooth transition of the temperature at the surface of the planets (see Sect. 4.1) from the time when they are embedded in the nebula to the time when they are exposed to the direct stellar irradiation.

The last term accounts for the heating by the surrounding environment (molecular cloud), which we set constant to Tcd = 10 K. We thus neglect possible variations of this background temperature depending on the stellar cluster environment in which a star and its planetary system are born (Krumholz 2006; Ndugu et al. 2018). On the other hand, different cluster environments and thus different levels of the interstellar FUV field (Fatuzzo & Adams 2008) are taken into account by varying in the population syntheses (see Paper II) the magnitude of the external photoevaporation rate wind. External photoevaporation is likely the most important environment-related factor for discs (Winter et al. 2020).

3.2.2 Disc photoevaporation

Photoevaporation in the protoplanetary discs is the principal means of controlling their lifetimes. For the prescription, we follow Mordasini et al. (2012b). In this scheme, we include contributions from both internal (due the host start itself) and external (due to nearby massive stars in the birthplace of the system) sources.

For the external photo-evaporation, we use the far-ultraviolet (FUV) description of Matsuyama et al. (2003). FUV radiation (6–13.6 eV) creates a neutral layer of dissociated hydrogen whose temperature is TI ≈ 103 K. The corresponding sound speed is then (7)

where the mean molecular weight μI = 1.35 for the dissociated gas. It corresponds to the gravitational radius (where the sound speed equals the escape velocity) of (8)

We assume that mass is removed uniformly outside of βIrg,I with βI = 0.14 (similar to Alexander & Pascucci 2012), so that the rate is given by (9)

with wind a parameter that provides the total mass loss rate if the disc would extend to rmax = 1000 au. In practice however, the actual mass loss rate due to external photoevaporation is clearly smaller than that parameter, as the disc does not extend up to rmax, but to a dynamically obtained radius which results from the interplay of viscous spreading (increasing the outer radius) and external photoevaporation (decreasing the outer radius).

For the internal photoevaporation, we follow Clarke et al. (2001), which in turn is based on ‘weak stellar wind’ case of (Hollenbach et al. 1994). Here, extreme-ultraviolet (EUV; > 13.6 eV) creates a layer of ionised hydrogen whose temperature is TII ≈ 104 K and with a mean molecular weight μII = 0.68. The sound speed and gravitational radius are computed in analogy with Eqs. (7) and (8). The scaling radius r14 = βIIrg,II∕1014 cm follows Clarke et al. (2001) while we select again βII = 0.14 following Alexander & Pascucci (2012). With this, we can estimate the base density with (10)

where we set kHol = 5.7 × 104 following the hydrodynamical simulations of Hollenbach et al. (1994) and , which is the ionising photon luminosity in the units of 1041 s−1. The distance-dependent base density can then be calculated as (11)

We further follow Clarke et al. (2001) to get outside of βIIrg,II.

The final photoevaporation rate is given by the sum of the effects of host star + nearby massive stars with (12)

3.2.3 Initial gas surface density profile and example

We initialise the gas surface density profile with (Veras & Armitage 2004) (13)

where r0 = 5.2 au is the reference distance, βg =0.9 the power-law index (Andrews et al. 2010), rcut,g the characteristic radius for the exponential decay and rin the inner edge of the disc.

The conversion between the total mass and the normalisation surface density Σg,0 at r0 is obtained with (14)

It should be noted that this formula neglects the lack of gas within rin, but since the total mass is dominated by the outer disc as we have a shallow power-law, there has in practice very limited effect.

An example of evolution of such as disc, without any planets (i.e. ), is provided in Fig. 3. The initial conditions and parameter are provided in Table 1 (note that the table also contains planetesimals disc properties that are not used here). The lifetime of that disc is nearly 5.3 × 106 yr.

The temporal evolution shows overall a decrease in the surface density. A hole forms inside roughly 2 au by about 4.7 Myr. The change in the temperature profile initially between 1.5 and 3 au and that moves inwards is due to a maximum in the opacity (top right panel). This different behaviour is reflected in the surface density as the temperature affects the sound speed, hence the viscosity.

For the midplane temperature, the direct irradiation term is only important in the innermost region (within about 0.2 au) until a few hundred thousand years before the dispersal of the gas disc. The last profile before dispersal shows an increase of temperature within 0.2 au due to this contribution; otherwise the midplane temperature remains below the equilibrium temperature, apart from the inner region (≪3 au) at early times.

We compared the results obtained here with prescriptions from other works, as for instance Bitsch et al. (2015a). We find that in general, for a given stellar accretion rate (which is the parametrisation of the Bitsch et al. 2015a prescription) we obtain lower surface density profiles by 30% coupled with larger temperature by 20–40%. The two models cannot be very well compared directly due to the different underlying assumptions, like constant radial flow rate in Bitsch et al. (2015a). Our models accounts for the full evolutionary equation for the surface density including photoevaporation and gas accretion bythe protoplanets, which means we have the radial flow rate varying with distance (bottom right panel of Fig. 3, where the radially constant inflow in the inner disc, and the viscous spreading (outflow) in the outer disc can be seen). There are other model assumptions that result in the differences between the surface density and temperature in the two models: (1) the stellar luminosity, which in our case it starts with roughly 3 L as predicted by the Baraffe et al. (2015) tracks whereas Bitsch et al. (2015a) begins with 1.5 L following Baraffe et al. (1998), (2) the opacity which affects the relation between midplane and disc photospheric temperature, and (3) the different approach of including stellar irradiation (vertically integrated assuming an equilibrium for the flaring angle versus an explicit 1D vertical structure with radiative transfer).

thumbnail Fig. 3

Time evolution of the surface density (top left), opacity κ (top right), midplane temperature (bottom left), and radial flow rate (bottom right) of a protoplanetary disc. The lines represent each one snapshot the state, and are spaced by about 2 × 105 yr. The blue line in both panels shows the initial profile, which has not yet been evolved at all, and is therefore not in equilibrium. The green line in the temperature profile shows the profile at disc’s dispersal, which is given by the equilibrium temperature with the host star’s luminosity.

3.3 Planetesimal disc

Planetesimals are represented by a fluid-like description, that is they are modelled not as individual particles but on a grid as a surface density (Σs) with eccentricity (eplan) and inclination (iplan) as dynamical state.

3.3.1 Dynamical state

For the time evolution of the dynamical state, we use the approach of Fortier et al. (2013) and explicitly solve the differential equations describing the change of eccentricity and inclination. In this framework, these are stirred by both the protoplanets, and to a lesser extent the other planetesimals, and damped by drag from the gas disc. The equations for the root mean square (RMS) of the planetesimals’ eccentricity eplan and inclination iplan read as

The contributions from the aerodynamical drag, stirring by the protoplanets and the planetesimals are denoted by ‘drag’, ‘VS,M’ and ‘VS,plan’ respectively. The dynamical state is followed during the entire formation stage. The drag term is onlyevaluated while the gas disc is still present. After the dissipation of the gas disc, the term is set to 0.

The form of the drag term depends on the regime: Epstein, Stokes (laminar) or quadratic (turbulent). The distinction between those regimes is made using the criterion proposed by Rafikov (2004) using the molecular Reynolds number Remol = vrelRplanνmol, where νmol = λcs∕3 is the molecular viscosity, the gas molecules’ mean free path, the number density assuming all of the gaseous molecules having hydrogen mass, their collisional cross-section, Rplan the planetesimals’ radius, (17)

their relative velocity, (18)

the deviation between the gas and Keplerian velocities due the support of the gas by the radial pressure gradient, ρmid the midplane gas density, and vK = Ωr the Keplerian velocity. When Rplan < λ, the gas drag is assumed to be in the Epstein regime. Otherwise, if Remol > 20, the gas drag is taken to be in the quadratic (or turbulent) regime and in the Stokes regime if not.

The expressions for the drag in the quadratic regimes are (Adachi et al. 1976; Chambers 2006),

where (21)

is the gas drag time scale and CD = 1.

In the Stokes regimes the drag expressions are

while in the Epstein regime they read as

(Adachi et al. 1976; Rafikov 2004; Fortier et al. 2013). We also want to point out that we do not model the formation of gap in the gas disc by giant planets. This means that drag in the vicinity of such planets might be overestimated, resulting in lower eccentricities and inclination. As consequence, the accretion rate of planetesimals would be overestimated in this stage, which affects the heavy element contents of the planets.

As in Fortier et al. (2013), the stirring by the protoplanets follows the approach of Guilera et al. (2010), where the stirring of Ohtsuki et al. (2002) is modulated with the separation from the protoplanets. The contribution reads as

where the sum is over all the protoplanets present in the system, (28)

is the modulation due to separation so that the perturbation is effectively restricted to the planet’s feeding zone, (29)

the planet’s Hills radius, and b = 5 is the half-width of the feeding zone (see Sect. 4.3.3). The terms PVS and QVS are given by (Fortier et al. 2013),

Here, plan = replanRH and are respectively the reduced planetesimals’ eccentricity and inclination, , β = iplaneplan, while for IPVS and IQVS we use the approximations obtained by Chambers (2006):

The stirring by the other planetesimals is given by, following Ohtsuki et al. (2002),

with (36)

and , the mass of a planetesimal.

To set the initial dynamical state, we assume that the disc is initially in a cold state, that is only the self-stirring of the planetesimals contributes to their eccentricities and inclinations. In other words, this assumes that the embryos appear instantly at the beginning of the simulation. The equilibrium values can be derived by equating the contributions of self-stirring and damping (Thommes et al. 2003; Chambers 2006), which results in (37)

and (38)

We also compared our prescription for the dynamical state with gamma-stirring from, for instance, Ida et al. (2008) and Okuzumi & Ormel (2013). Although this is not straightforward due to the differences in the sources, we find that, generally, the eccentricities resulting from γ-stirring are larger than the self-stirring from the planetesimals, but lower than the stirring by the forming protoplanets. Thus, accounting for the stirring of planetesimals by turbulent diffusion in the disc would increase their eccentricities at locations far away from growing protoplanets. Close to the growing protoplanets however, where the planetesimals’ eccentricities are important for the solids accretion rate, neglecting this effect does not significantly affect planetesimals’ eccentricities.

Table 1

Initial conditions and parameters for the example system.

3.3.2 Size, initial surface density profile, and evolution

To roughly take into account the observational (e.g. Ansdell et al. 2018) and theoretical (e.g. Birnstiel & Andrews 2014) finding that solids have a more concentrated distribution than the gas, the initial surface density profile of planetesimals now follows a steeper slope than the one of the gas disc (Lenz et al. 2019; Voelkel et al. 2020). This leads to a higher concentration of solids in the inner part of the disc.

As already in the Generation II Model (Alibert et al. 2013), we assume a constant planetesimals radius of 300 m throughout the disc, which is a strong assumption and simplification. There is an ongoing discussion about the characteristic primordial planetesimal size in the literature. Observations of extrasolar debris belts (Krivov & Wyatt 2021), the presence of hypervolatile ices in comets that can only be preserved in impacts involving small bodies (Golabek & Jutzi 2021), direct size determinations by stellar occultations (Arimatsu et al. 2019) and some theoretical studies (Fraser 2009; Schlichting et al. 2013) suggest small (~1 km) characteristic planetesimals sizes. On the other hand, the absence of small craters on Pluto (Singer et al. 2019), the size distribution in the asteroid belt(Morbidelli et al. 2009), and the theoretical predictions of planetesimal formation models (e.g. Klahr & Schreiber 2020) rather point at ~100 km planetesimals. The first two points can, however, also be explained with other effects (Zheng et al. 2017; Wei et al. 2018, although the former work makes no determination about the initial size frequency distribution of planetesimals).

In the more specific context of the simulations presented here, this choice was made for the following reasons: (1) small planetesimals undergo sufficient eccentricity and inclination damping by the disc gas to sustain a planetesimal accretion rate in the oligarchic growth regime that is high enough to build giant planet cores during typical disc lifetimes (Fortier et al. 2013). We note that the Generation I and Ib Bern Models assumed in contrast runaway planetesimal accretion as Pollack et al. (1996). In the runaway regime, the eccentricities and inclinations of the planetesimals are assumed to remain low even without damping by the disc gas. Therefore, fast core growth occurred in these models also with 100 km planetesimals, which was the assumed size in these early model generations. (2) Their drift time scales are longer than typical lifetimes of gas discs (Burn et al. 2019) and (3) this size was shown to be able to reproduce several of the known exoplanet properties across a wide range of masses (Fortier et al. 2013). In any case, the constant planetesimal size is an important limitation of the model. Including in the Bern Model an explicit model for the evolution of the solid building blocks across the entire size range (dust-pebble-planetesimals) is thus subjectof ongoing research. A first important step was recently made in Voelkel et al. (2020) where we have coupled the dust-and-pebble model of Birnstiel et al. (2012) and the planetesimal formation model of Lenz et al. (2019) to our global model. These effects are, however, not yet included in the Generation III Model presented here.

To set the initial surface density profile of planetesimals, we thus use a slightly different description than for the gas, that is, (39)

with the power-law exponent is set to βs = 1.5, as in the MMSN, and rcut,s = rcut,g∕2 is the exponential cutoff radius of the solids, set half the value of the gas disc following Ansdell et al. (2018). This formula also enables us to model relatively sharp outer edges of the solids disc (Birnstiel & Andrews 2014).

The reference surface density value Σs, 0 is adjusted so that the bulk solids-to-gas ratio remains to the prescribed value (e.g. 1%).

The surface density of planetesimals is reduced by accretion onto and ejection by the protoplanets to ensure mass conservation (see Sect. 4.3), or removed entirely if . Our model only includes ejection (Sect. 4.3.2) and not scattering by the forming planets. Thus, we do not have redistribution of material to other regions of the disc by planets, as obtained by Raymond & Izidoro (2017) for instance. Finally, the planetesimals disc remains after the dispersal of the gas disc; the only difference is that the damping terms for eccentricity and inclination vanish.

3.3.3 Compositional model

The Bern model includes the simple condensation model of Thiabaud et al. (2014) and Marboeuf et al. (2014a). The initial abundance of volatile and refractory species is identical to the one given in Marboeuf et al. (2014b). Volatile species are composed of H, O, C, and S atoms whose abundance reflect solar composition (Lodders 2003). The relative abundances of the molecules are set according to interstellar medium. Then at each location in the disc at t = 0, we check whether each molecule is the solid or gas phase assuming local thermodynamical equilibrium. This yields the fraction of heavy elements that is locally condensed and thus contributes to the solid surface density (the ice line locations), and the chemical composition. This composition is tracked into the protoplanets when a propotoplanet accretes planetesimals, and in giant impacts between protoplanets. This yields in particular the final iron to silicate ratio and the volatile mass fraction of all the planets.

The factor fs(r) in Eq. (39) for the initial planetesimal surface density accounts for the mass fraction of all elements that are in the solid phase at a given location. To compute its value, we use the aforementioned condensation model. Only the contributionof molecules in the solid phase are accounted for the resulting solid surface density. Thus, the value of fs in the inner locations is the mass fraction of condensed to total solids and this value increases by small jumps each time an ice line is crossed until it becomes unity at large separation.

For the density of planetesimals, we assume that in the region where only refractory materials contributes to the solid phase ρplan = 3.2 g cm−3 while when volatiles are in the solid phase we take ρplan = 1 g cm−3. This transition corresponds to the H2O-ice line in all discs, which induces the largest surface density jump because H2 O makes up ~60% of all ices in mass (Marboeuf et al. 2014a).

3.3.4 Example

An example of the dynamical state of planetesimals is provided in Fig. 4. The initial conditions and parameters are provided in Table 1. This is the same initial disc as shown in Fig. 3, except than ten embryos were added to the disc, at the locations shown by the dashed vertical lines. In addition, both migration and N-body interactions were artificially disabled so that the embryos remain at the same location throughout the simulation.

The different jumps in the initial surface density profile are due to the crossing of the different ice lines; the most consequential one at about 3 au is due to the water-ice line. The surface density of planetesimals is equalised inside the feeding zone of each planet. It should be noted that besides this effect, we do not include planetesimals redistribution, as was found by, for instance, Levison et al. (2010). In total, the planets accreted 61 M of planetesimals (47 M of which by the giant planet) while 89 M were ejected(according to the prescription detailed in Sect. 4.3.2; virtually all of them by the giant planet). The feeding zones are all nearly depleted by the planets due to accretion, except for the giant planets where 65% of planetesimals were ejected and 35% accreted.

The stirring by the protoplanets heats the planetesimals in the surrounding region. This effect is heavily dependent on the protoplanet’s masses; the most massive one is the second outermost one (close to 10 au), which reaches a mass of about 5.4 M at the end of the formation stage. That planet has a core mass of 47 M, which corresponds (for a pure H/He envelope) to a metallicity slightly lower than that of the star (2.8% versus 3.0%). This is below the relationship found by Thorngren et al. (2016) for the planet’s mass. This, however, is not unexpected for the idealised setup used here: first, planets that form in the in-situ case tend to have lower core masses than planets that migrated (e.g. Alibert et al. 2005b). Second, with N-body interactions switched off here, giant impacts otherwise increasing the heavy element content are not possible. In the more realistic example in Sect. 8.1.2, where these effects are included, giant impacts strongly increase the solid content of the giant planets by a factor 2–3 relative to value at the moment gas runaway begins. The impacts are themselves triggered by the fast mass growth, destabilising neighbouring lower mass protoplanets.

As noted by Fortier et al. (2013), the usual assumption that β =iplaneplan ≈ 1∕2 does not hold. We find that the stirring of eccentricities takes place over larger separation to the protoplanets than for the inclinations. This can be seen for instance in the region affected by the most massive planet.

Further,the effect of the planet is not only limited to the surrounding area because of the following effect: the massive planet is able to significantly reduce the inward gas flow such that the region inside its orbit becomes gas-poor. This greatly reduces the damping of the planetesimals dynamical state to a such point that their eccentricity becomes close to unity.

4 Planet properties

4.1 Envelope structure

In the Bern model, the internal structure of the planets (and thus their gas accretion rate, radius, luminosity, and interior structure) are found at all stages (attached, detached, evolution) by directly solving the 1D structure equations. In contrast, many other global models use in contrast approximations and fits to find for example the gas accretion rate (see Alibert & Venturini 2019). While the 1D hydrostatic picture is also not the final word for low-mass planets because of multidimensional hydrodynamic effects (e.g. Ormel et al. 2015; Lambrechts & Lega 2017; Cimerman et al. 2017; Moldenhauer et al. 2021), the fits (except the deep neural networks) often fail grossly to reproduce the result of 1D structure equations that they should in principle recover (Alibert & Venturini 2019). Many fits also neglect the influence of the luminosity on the gas accretion rate (e.g. Ida & Lin 2004a; Bitsch et al. 2015b). In reality, there is an important interplay between solid accretion which is dominant for the luminosity at early stages, and gas accretion. This leads to important feedbacks that can only be captured when solving the internal structure equations (Dittkrist et al. 2014).

Also, from the point of view of guiding and interpreting astronomical observations, it is crucial to solve the internal structure equations, as this gives self-consistently at each moment in time the planet’s radius and luminosity and associated magnitudes. These are the observable quantities for transit and direct imaging surveys. By predicting them self-consistently, the output of the Bern model can be compared in population syntheses not just with methods measuring quantities depending on the planets’ mass (like RV, astrometry or microlensing), but also transit and direct imaging surveys.

The downside is that solving the internal structure for bodies ranging in mass from 10−2 M to beyond the deuterium limit requires an internal structure model that is very versatile and numerically stable in all stages of planetary formation and evolution. Solving the internal structure also comes with significant computational cost.

thumbnail Fig. 4

Time evolution of the RMS of planetesimals’ eccentricity (top left), inclination (top right), and surface density (bottom left) of a circumstellar disc that also contains 10 embryos. The lines represent temporal snapshots of the three quantities, and are spaced by about 2 × 105 yr. The blue line in both panels denote the initial profile. The dashed vertical lines represent the location of the embryos, which is fixed in this case. N-body interactions were also disabled. The lifetime of the gas disc is shorter than the case presented in Fig. 3 due to the accretion by the protoplanets.

4.1.1 Attached phase

In the initial phase, known as the attached phase, the envelope is in equilibrium with the gas disc and the gas density smoothly transitions from the value in the protoplanetary envelope to the one in the background nebula. The planets do not yet have a well-defined outer radius. During this phase, the gas accretion rate is governed by the ability to radiate the gravitational energy liberated by the accretion of solids and gas, and the envelope’s contraction. For the forming giant planets, this phase generally lasts until the planets reach a total mass in the range of 30–100 M where envelope contraction becomes fast, depending on the conditions. There is no fixed mass boundary; the transition occurs when the gas accretion rate obtained from solving the internal structure equations (that is the envelope’s Kelvin-Helmholtz contraction) becomes larger than the disc-limited rate (Sect. 4.1.2. For low-mass planets which have very low gas accretion rates (very long Kelvin-Helmholtz timescales), the attached phase lasts (almost) until the gas disc dissipates.

Gas accretion is calculated by solving the classical 1D radially symmetric internal structure equations (Bodenheimer & Pollack 1986),

with M the mass enclosed in the radius R, P the pressure, T the temperature, ρ = ρ(P, T) the density, computed using the SCvH equation of state (Saumon et al. 1995), and ∇ad and ∇rad the adiabatic and radiative gradients respectively. The minimum of the two indexes is the Schwarzschild criterion (e.g. Kippenhahn & Weigert 1994), and is used to ensure stability against convection. The adiabatic gradient comes from the equation of state, while the radiative gradient is given by (43)

with L being the luminosity.

The opacity in the envelope κ is obtained in similar way as for the gas disc, but following Mordasini et al. (2014), the interstellar medium (ISM) grain opacity contribution in Bell & Lin (1994) is multiplied by a factor fopa = 0.003. This value was found in Mordasini et al. (2014) to fit best the detailed simulations by Movshovitz & Podolak (2008) and Movshovitz et al. (2010) of the grain dynamics in protoplanetary atmospheres (growth, settling) and the resulting dust opacities. Using one global reduction factor of the ISM opacity can of course not reproduce the full complex behaviour of the grain opacity which depends on planetary properties like the core or envelope mass as found in grain dynamics models (Movshovitz & Podolak 2008). But as shown in Mordasini et al. (2014), it still provides a useful first approximation. The value is not increased when a planetary system with higher metallicity is simulated. The reason is that a higher dust input in the outermost layer (as possibly associated with a high metallicity system) does not lead to a strong increase of the opacity. This was found numerically in Movshovitz & Podolak (2008) and explained analytically in Mordasini (2014): a higher dust input leads to a higher dust-to-gas mass ratio (which increases the opacity), but also larger grains (which decreases the opacity). These effects cancel each other out in the dominating growth regime of differential settling.

The boundary conditions for the integration are taken as follows: the outer radius is given by, following Lissauer et al. (2009), (44)

where (45)

is the Bondi radius, RH is the Hill’s radius (Eq. (29)), k1 = 1 and k2 = 1∕4. The pressure and temperature are derived from the local properties of the disc with

and (48)

being the optical depth at the surface of the planet (Mordasini et al. 2012c), using the reduced opacities for the grains. The more complex parts come from the luminosity and the mass. The calculation of the outer luminosity L(Rtot) is described in Sect. 4.2. In the case of the mass, what is known is the core mass, that is M(Rcore) = Mcore, while M(Rtot) = Mtot is the quantity that is being searched for. We thus use an iterative method by guessing Mtot, which is then used tointegrate the internal structure equations until the boundary condition at the inner boundary is fulfilled, that is M(Rcore) = Mcore. Once Mtot is found, the envelope mass can be retrieved by Menv = MtotMcore, and the gas accretion rate by taking the difference of the envelope mass between two successive steps of the envelope structure calculation env = (Menv(t) − Menv(t − Δt))∕Δt.

4.1.2 Maximum gas accretion rate

In the initial stages, the gas accretion is limited by the planet’s ability to radiate away the potential energy provided of the accretion material, that is the Kelvin-Helmholtz process. The rate at which gas can be accreted is set by the Kelvin-Helmholtz time scale, (49)

However, as the planet’s core reaches a mass of about 10 M, the value of τKH becomes so low that the planet undergo runaway gas accretion. In this phase, the amount of gas that the planet can accrete is constrained by the supply from the gas disc. Therefore, we compute the quantity env,max, which is used to limit the value of env found by solving the internal structure equations.

Our approach to compute the maximum rate is similar to Mordasini et al. (2012c) but using only the ‘local reservoir’ component. This a major difference from the previous versions of the Bern model, where gas accretion was constrained from the radial flow of the gas. Following D’Angelo & Lubow (2008) and Zhou & Lin (2007), we adopt a Bond- or Hill-like accretion in a region of size Rgc around the planet. For simplicity, we compute Rgc according to Eq. (44). Depending on the value of Rgc with respect to H, the local disc’s scale height, two different regimes occur. In the case where Rgc< H, the planet will not accrete from the full vertical extent of the disc, and so the gas flow through the gas capture cross section is given by (50)

with ρ ≈ Σ∕H the approximate density of the gas and the relative velocity between the gas and the planet.

On the other hand, in the case Rgc > H, the planet will accrete from the whole gas column and the approximation of constant gas density breaks down. In this situation, the gas flow through the planet’s capture radius is provided only by the radial extension of the gas capture area, hence we have (51)

To distinguish between the two regimes, we use the lower rate of the two, that is (52)

Finally, to ensure that no more gas than available in the feeding zone Mfeed is accreted during one time step, we further constrain env,max < Mfeed∕Δt. We consider the limiting case to be that gap formation does not reduce the planetary gas accretion rate. Such a situation arises if the eccentric instability (Papaloizou et al. 2001; Kley & Dirksen 2006) allows the planets to efficiently access disc material even after a gap has formed. For circular orbits, gap formation would in contrast strongly reduce the gas accretion rate (Lubow et al. 1999; Bryden et al. 1999), and limit planetary masses to ~ 5−10 M.

The radial extent of the feeding zone is set by

with ffeed = 0.5 so that the overall extent is a half a Hill radius larger than the radial excursion of the planet’s orbit. This radial extent provides the location over which the disc properties (Σ, H, etc.) are averaged forthe calculation of the maximum rate and the removal of the accreted gas, with (55)

The planet’s eccentricity consequently does not directly affect the maximum gas accretion rate, but only indirectly through the size of the feeding zone. The self-limitation of gas accretion by removal of local disc gas by the planet, which then needs to be replenished by the inflow from more distant disc regions (i.e. mass conservation) is fully taken into account in our scheme via the term entering the evolutionary equation of the disc gas surface density. We also take into account that for planets of any mass growing in multi-planet systems, the eccentricity can be increased via gravitational planet-planet interactions, which then affects the feeding zone width and thus indirectly the gas accretion rate. On the other hand, we currently do not take into account that the eccentric instability (i.e. the increase of a single giant’s eccentricity because of gravitational interaction with the gas disc) in reality only acts for sufficiently massive planets (Papaloizou et al. 2001; Kley & Dirksen 2006). This could lead to an overestimation of gas accretion at lower to intermediate giant planet masses. This could potentially explain why our current model of disc-limited gas accretion seems to too strongly reduce the stellar gas accretion rate (Manara et al. 2019; Bergez-Casalou et al. 2020). Gap formation would reduce this effect, but could potentially lead to another issue: observationally, the giant planet mass function seems to extend smoothly to about 30 M (Sahlmann et al. 2011; Adibekyan 2019) (though Santos et al. 2017 and Schlaufman 2018 found a change in the metallicity dependency at about 4 to 5 M and concluded that the planets above that threshold formed predominantly by gravitational instability). Reaching such high masses could be difficult given the expected reduction of gas accretion because of gap formation in the circular case.

The reduction of gas inflow into the inner disc because of an accreting giant planet can result in the clearing of the inner region of the protoplanetary disc by photoevaporation (Rosotti et al. 2013). This effect is also automatically taken into account by our model.

To compare the prescription presented here with previous work, we provide in Fig. 5 the comparison of the gas accretion rate for the second outermost planet from the case shown in Fig. 4. The previous methodology, using the radial gas flow and taking into account the geometry was described in Mordasini et al. (2012c), with a limit of 0.9 of the radial flow to allow some gas to flow through the gap (Lubow & D’Angelo 2006). The results show that using the Bondi rate, as we presented here, gives a somewhat stronger limitation of gas accretion by the forming planet, especially during the onset of the runaway gas accretion. As a result, the final planet’s mass is a bit lower when using the Bondi rate.

thumbnail Fig. 5

Comparison of two prescriptions for the maximum (i.e. disc limited) gas accretion rate, the one presented in this work (labelled ‘Bondi rate’) with that of Mordasini et al. (2012c) (labelled ‘Flow rate’). These are two different simulations (one for each prescription) whose initial conditions represent the second outermost planet in Fig. 4. Top panel: maximum value that can be supplied by the gas disc (labelled ‘Max.’) and effective accretion rates (labelled ‘Eff.’), which is given by intrinsic cooling in the initial attached phase and the maximum rate in the detached phase. Bottom panel: corresponding enveloppe mass (i.e. total gas accreted).

4.1.3 Detached phase

Once the gas accretion rate exceeds the maximum that can be provided by the disc – which includes the planet no longer being in a region where gas is present – the accretion regimes changes to the detached phase (Bodenheimer et al. 2000). In the detached phase, the solid and gas accretion rate are known (for the gas, it is given by the disc-limited rate), but not the planet’s radius. The radius is determined following the approach of Mordasini et al. (2012c,b), that is by using the same internal structure equations as in the attached, but iterating on the radius until convergence is reached.

The pressure outer boundary conditions are modified to take into account that the disc and the envelope are no longer connected, and that the gas free-falls onto the surface of the planet (56)

with Pneb(aplanet) being the pressure at the midplane of the gas disc, Pedd = (2g)∕(3κ) the Eddington expression for the photospheric pressure due to the material residing above the τ = 2∕3 surface, Prad = (2σSBT4(Rtot))∕(3c) the radiation pressure, c being the speed of light in vaccum, and (57)

being the ram pressure due to the accretion shock and the free-fall velocity at the surface of the planet.

4.1.4 Evolutionary phase

For the evolutionary phase (after the dispersal of the gas disc), the outer boundary conditions are set to

where is the intrinsic temperature, , and A = 0.343 is the albedo, which is taken be the same as Jupiter (Guillot 2005). This value was selected for simplicity, although hot-Jupiter planets may have lower values (e.g. Mallonn et al. 2019).

We thus use an Eddington grey boundary condition taking the stellar irradiation into account, as described in Mordasini et al. (2012c). During evolution, we assume a solar-composition condensate-free gas for the opacities, using the opacity tables of Freedman et al. (2014). Nebular grain opacity is neglected, at they are found to rain out quickly once gas accretion stops (Movshovitz & Podolak 2008). The identical envelope and atmospheric composition (pure H/He, solar composition opacities) in all planets means that for planets with identical bulk properties (orbital distance, core and envelope mass), the predicted radii will exhibit an artificially reduced spread. In reality, planets have different enrichment levels of heavy elements in the envelope (e.g. Fortney et al. 2013). This affects the equation of state and opacity, resulting in particular in a larger spread of the radii (e.g. Burrows et al. 2011; Müller et al. 2020).

4.1.5 Example

To illustrate the calculation of the internal structure, we provide snapshots of envelope structures in Fig. 6 and the time evolution of the radius and luminosity in Fig. 7. These are taken from the second outermost planet of the system shown in Fig. 4, which is a giant planet whose final mass is 6.4 M. Due to the different scales involved in the attached, detached, and evolutionary phases, they are shown in different panels. During the attached phase, the structure extends to the Bondi radius (Eq. (45)), which is much larger than the core radius. Therefore, the structure spans a wider range of pressure. The upper part of the envelope is radiative while the lower part is convective, with several transitions in the mid region. The red profile shows the internal structure at the beginning of the transition from the attached to detached phase (the time marked with a dashed vertical line in the insert in Fig. 7).

Note that the planet is still accreting during the initial stages of the detached phase.

4.2 Luminosity

4.2.1 Accretion and contraction

The luminosity calculation suffers from the same problem as the total mass in the attached phase, or the outer radius in the detached phase; that is that the new structure needs to be known to retrieve its energy, hence the luminosity. This means that the total energy of the new structure needs to be estimated for a luminosity to be obtained.

The model uses the approach from Mordasini et al. (2012c). The total energy is given as (60)

with u being the specific internal energy of the gas, as obtained from the equation of state. The gravitational binding energy term includes the contribution from the core. For simplicity, we assume that it has a constant density, so its contribution is taken as . It should be noted that this is not strictly self-consistent with our model to determine its density or radius, which assumes differentiation (Mordasini et al. 2012b); however, the difference remains small (Linder et al. 2019). The parameter ξ in Eq. (60) represents as in polytropic models the mass distribution and additionally the thermal energy content. It is retrieved from Eq. (60). The internal luminosity resulting from the accretion, cooling, and contraction Lint can then be obtained as (61)

with tot = core + env being the total accretion rate of the planet (solids and gas). The value tot in the attached phase and of tot in the detached phase are determined from the guess for the mass or radius during the iterations. The same is not true for . To circumvent this problem, we estimate the luminosity with (62)

The correction factor C corrects for neglecting the term. The value of C can be calculated a posteriori by determining the actual total energy of the new planet, with C = −(Etot(t) − Etot(t − Δt))∕(Lint ⋅ Δt). The value of C is then used for the next time step.

Marleau et al. (2017, 2019b) conducted 1D radiation-hydrodynamic simulations of the planetary gas accretion shock, a feature that is seen in various 3D radiation-hydrodynamic simulations of accreting protoplanets of sufficiently high mass (e.g. Szulágyi & Mordasini 2017; Schulik et al. 2020). High postshock entropies were found, suggesting that warm or hot gas accretion is more plausible than cold accretion (see also Berardo et al. 2017; Berardo & Cumming 2017). We therefore assume in our model that gas accretion in the detached phase is hot, which means that we do not subtract the accretion shock luminosity from Lint (see Mordasini et al. 2012c).

In addition to the accretion and contraction luminosity, we include the luminosity from radioactive decay, bloating for close-in planets, and, in the case of brown-dwarfs, deuterium fusion. The radiogenic luminosity Lradio includes contributions from the three most important long-lived radionucleides 40K, 238U and 232Th (Wasserburg et al. 1964). To compute the luminosity contributions, we follow the procedure of Mordasini et al. (2012b): we assume the mantle of the protoplanets has a chrondritic composition and the energy production rate are retrieved from meteoritic values of William (2007). The initial radiogenic contribution is Q0 ≈ 5 ×−7 erg g−1 s−1 of mantle material (all elements besides iron).

thumbnail Fig. 6

Snapshots of the internal structure of the second outermost planet of Fig. 4. The structures are split according tothe phases, with attached (top left), transition (the initial stage of the detached phase; top right), detached (bottom left) and evolutionary (bottom right). The red line shows the first profile of the detached phase and is shown of both panels. The green and blue profiles lie at the transition between two stages and are shown of two panels each. In each profile, thin lines show the part where energy transport is radiative and thick lines for convective.

4.2.2 Bloating of close-in planets

Massive, close-in planets exhibit anomalously large radii (Laughlin et al. 2011). To reproduce this effect, we include a bloating mechanism based on Sarkis et al. (2021). For planets which are in the detached and evolutionary phase and directly irradiated by the host star, we include an additional luminosity contribution that is based on the best empirical fit formula obtained by Thorngren & Fortney (2018): (63)

with (64)

the total stellar flux at the planet’s location, and τmid is the optical depth from the star to the planet location through the mid-plane of the disc, as in Eq. (6). We only apply bloating if the stellar flux F (in the evolutionary phase) or the stellar flux multiplied by the optical depth (before the dispersal of the gas disc) is greater than 2 × 108 erg cm−2 s−1 (Demory & Seager 2011).

thumbnail Fig. 7

Time evolution of the planet’s radius and luminosity of the second outermost planet of Fig. 4; the same as in Fig. 6. The insert on the left panel shows the contraction at the transition between attached and detached phase. The exact time where the model switches between the two phases is shown with the vertical dashed line.

4.2.3 Deuterium-burning

For the calculation of the luminosity due to deuterium fusion, we follow the procedure of Mollière & Mordasini (2012). In this framework, the energy generation rate (per unit mass and time) is given by Kippenhahn & Weigert (1994), with the assumption that nuclei are fully ionised and non-degenerate. The energy released in each process is computed according to Fowler et al. (1967). The specific deuterium burning luminosity of a planet depends on the conditions in the planet’s gaseous envelope, most notably the density, temperature, and the remaining deuterium nuclei. This implies that there is no universal mass at which deuterium burning starts, but as already found in Mollière & Mordasini (2012) (see also Bodenheimer et al. 2013), the mass where burning becomes important clusters around about 13 M. The presence of a solid core does thus not significantly alter the mass where burning starts relative to (coreless) brown dwarfs (Chabrier & Baraffe 2000). We use an initial deuterium number fraction [D/H] = 2 × 10−5, which is the primordial value.

Our model also includes the enhancing of the reaction rate by screening, that is the shielding of the positive charges by the surrounding electron. In turn, screening is affected by the electron degeneracy, as we are dealing with objects of high central densities. This procedure follows the work of Dewitt et al. (1973) and Graboske et al. (1973).

4.2.4 Total luminosity

The final luminosity is then given by (65)

We assume that at a given time, the luminosity does not change within the envelope, that is ∂L∂r = 0. This approximation is fine under most circumstances because energy transport is due to convection and the luminosity enters only in the radiative gradient. During rapid gas accretion in the detached phase, under the effect of hot accretion, the interior may become radiative (Berardo et al. 2017; Berardo & Cumming 2017) and we do not account for the decrease of the luminosity with depth. This will be addressed in future work.

4.3 Accretion of solids

The growth of the astrophysical core of the planets can occur via three channels: (1) the accretion of planetesimals (e.g. Greenzweig & Lissauer 1992; Thommes et al. 2003), (2) the accretion of pebbles (e.g. Ormel & Klahr 2010; Johansen & Lacerda 2010; Lambrechts & Johansen 2012), and (3) by the collision with other embryos (which we call giant impacts). In the Generation III model, we consider accretion by planetesimals and giant impacts; the inclusion of pebble accretion is subject of ongoing work (Voelkel et al. 2020).

For planetesimals accretion, core growth is given by the probability of collisions with planetesimals in the oligarchic regime (Ida & Makino 1993), as described in Fortier et al. (2013). This is a major difference to the first generation of the Bern model which followed Pollack et al. (1996) for the planetesimal accretion rate. According to Chambers (2006), the core growth can be computed assuming a particle-in-a-box approximation is (66)

with the mean surface density of planetesimals in the planet’s feeding zone and pcoll the collisionprobability with planetesimals. As Ida & Lin (2008), we use the same prescription to calculate the planetesimal accretion rate independently of a protoplanet’s orbital migration rate. In addition, we address the possible impact that orbital migration could have in the context of the shepherd/predator regimes proposed by Tanaka & Ida (1999): in the idealised situation studied by Tanaka & Ida (1999) (single protoplanet per disc, no local reservoir of planetesimals, no growth via collisions with other protoplanets), shepherding was found to significantly reduce the planetesimal accretion rate for protoplanets migrating sufficiently slowly. However, in the more realistic N-body simulations by Daisaka et al. (2006) where multiple protoplanets (oligarchs) form and grow concurrently as expected in the oligarchic regime (Kokubo & Ida 1998), the trapping of planetesimals by the protoplanets is only tentative and does not significantly reduce their accretion rates. We similarly find that in the more realistic situation we consider here with many embryos per disc, the existence of a local reservoir of planetesimals in a protoplanet’s initial feeding zone accessible without migration and a time sequence of a solid accretion dominated initially by planetesimals and later on collisions with other protoplanets (Sect. 8.1), shepherding should only be of limited importance. We discuss these points further in Appendix A.

thumbnail Fig. 8

Evolution of the planet’s core, planetesimals capture, and total radii as function of time (left panel) and core mass (right panel) for the second outermost planet of Fig. 4. In the right panel, we also include two time indicators at 105 and 106 yr with dashedvertical lines. The upper-right portion of the plot for the total radius in the left panel is also shown on the left panel of Fig. 7.

4.3.1 Capture probability

We distinguish three different accretion regimes depending on the random velocities: low-, mid- and high-velocity. The distinction is based on the reduced planetesimals’ eccentricity plan = replanRH and inclination (r is the heliocentric distance): the high-velocity regime for , mid-velocity for and low-velocity for . According to Inaba et al. (2001), each regime has a different expression for the collision probability,

where Rcap is the planetesimal capture radius of the planet, β = iplaneplan and the IF and IG functions can be approximated as, following Chambers (2006):

The final collision is then given by (72)

In the initial stage, the capture radius Rcap is the physical radius of the core Rcore. Once the planet has sufficiently massive H/He envelope, it will enhance the capture cross-section of planetesimals. As in Fortier et al. (2013), the capture radius is obtained following Inaba & Ikoma (2003) by solving the implicit equation (73)

The enhancement of the capture radius over the physical radius is very important for increasing the planet’s planetesimals accretion rate (Podolak et al. 1988; Venturini & Helled 2020). We highlight this in Fig. 8, which compares the planetesimals capture radius to that of the core for the same planet we highlighted in Fig. 7. The calculation of the envelope structure begins at about 104 yr, before that, the capture radius is equal to that of the core. At that moment, the core mass is 9 × 10−2 M. By the time the core reaches 1 M at 4.8 × 105 yr, the capture radius is 9 times the core radius. Therefore, for small roughly km-sized planetesimals as in our case, the enhancement of the capture radius is already important for low-mass bodies (starting about × 10−1 M), and the calculation of gaseous envelopes cannot be omitted at any stage. Besides the factor that the eccentricity and inclination damping by nebular gas drag is more efficient for smaller planetesimals which leads to a larger gravitational cross section (a larger Safronov factor), the larger envelope drag enhancing the planet capture radius further is the second effect making the accretion of small planetesimal more efficient. This reflects that the accretion of km-sized planetesimals is not a pure gravitational process.

thumbnail Fig. 9

Illustration of the procedure to separate planetesimals’ feeding zones when zones would otherwise overlap. The horizontal axis represents the separation to the central star and four planets are shown. The light colour areas below the horizontal line show the initial feeding zones while the ones above show the final zones. amid2,3 and amid3,4 are the edges of the new feeding zone.

4.3.2 Ejection of planetesimals

Planets not only accrete material; they also induce gravitational perturbations on the planetesimals that come close-by but are not accreted. These planetesimals, if they receive a sufficient velocity kick from a close approach by a planet, can be ejected from the system. To estimate this effect, we follow a procedure similar to Ida & Lin (2004a). The planetesimals that receive a velocity kick greater than the escape velocity from the primary, , will likely be ejected from the system. Thus, we have that the fraction of accreted-to-ejected planetesimals is (Ida & Lin 2004a) (74)

with the characteristic surface velocity. The rate at which planetesimals are removed from the disc is then (75)

It should be noted that our model does not include the redistribution of planetesimals by scattering (Raymond & Izidoro 2017).

4.3.3 Feeding zone

To obtain the mean surface density of planetesimals in the feeding zone, we must determine its extent. The half-width of the feeding zone(centred at the planet’s location) is usually given in terms of the Hill radius with (76)

For a planet on a circular orbit, conservation of the Jacobi energy implies that (e.g. Hayashi et al. 1977). So, in a quiescent disc with , . For numerical stability reasons, however, we assume b = 5, as in Fortier et al. (2013).

In the general case, to account for a non-circular orbit of the planet, we take the feeding zone to span from rperiRfeed to rapo+ Rfeed, with rperi and rapo being the peri- and apocentre of the planet’s orbit respectively.

When multiple planets are present in the same disc, their feeding zones may overlap. To avoid problems with two planets accreting from the same location, such as mass-conservation issues, we separate the feeding zones so that there is at most one planet accreting at any location the disc. A graphical representation of the following procedure is provided in Fig. 9. First, we compute regions in the disc from where planets accrete. In the case a region contains a single planet, then the feeding zone is the same as in the single planet case (as in Region 1 on Fig. 9). If there are multiple planets in one region (as in Region 2 on that figure), the inner edge of the innermost planet and the outer edge of the outermost planet are set to the edges of the region. For the other edges, we sort planets by distance, and for each pair, we compute the location of the limit between their feeding zones with (77)

where the subscripts indicate the inner (in) and outer (out) planets of the pair. We scale with the square root of the planet masses because the area of the feeding zone scales with the square of the distance. This scaling keeps the area of the feeding zones related to the planet masses. We tested alternative prescriptions, like using the cubic root of the mass (as in the Hill sphere) or the mid-point between the two planets and found that the prescription does not significantly affect the outcomes of the simulations.

4.3.4 Core radius

To obtain the radius of the core (and its density), we applied a methodology similar to Mordasini et al. (2012b). This model also accounts forthe composition of the core and the pressure burden exerted by the envelope.

The principle is to solve similar structure equations as for the envelope, that is Eqs. (40) and (41), but with an equation of state that takes the form of a modified polytrope from Seager et al. (2007), which reads (78)

We include three different materials: iron, silicates (perovskite, MgSiO3) and ice, whose parameters ρ0, c and n are taken from Seager et al. (2007). Because of the small thermal expansion coefficient of these materials compared to H/He, we neglect via the temperature-independent modified polytropic EOS a possible temperature dependency of the radius of the core. It should, in any case, besmall (Grasset et al. 2009).

For gas giant planets, where envelopes can reach masses of thousands of Earth masses, this can cause a significant compressionof the core (Baraffe et al. 2008). Thus, the pressure on the core’s surface is taken as boundary condition of the calculation to include this effect. Core compression can be observed in Fig. 8, where the core radius shrinks after the envelope contracts at 1.06 Myr.

The core composition is retrieved from the accreted planetesimals described in Sect. 3.3.3 and other embryos in case of giant impacts. The chemical composition is used to obtain the fraction of the different elements to compute the core radius. While in the chemistry model includes 32 (Thiabaud et al. 2014) refractory and 8 volatile (Marboeuf et al. 2014b) chemical species, the core radius calculation groups them into only three types: iron, silicates, and water ice. Thus, we map all ice species to water ice when calculating the core structure and all refractories except iron to the silicate mantle. The reason for this is that first, equations of state are only available for a limited number of species. Second, the differences between different types of, for instance, silicates is not very large (Seager et al. 2007).

4.4 Atmospheric escape

During the evolutionary phase, that is after the dissipation of the gaseous disc, planets at small distances of their host star (~ 0.1 au) receive intense XUV stellar irradiation, which will drive atmospheric escape. This effect is especially important for the low-mass planets, that can loose the whole of their gaseous envelope due to their low gravitational binding energy (e.g. Lammer et al. 2009; Lopez et al. 2012; Owen & Jackson 2012; Jin et al. 2014; Jin & Mordasini 2018).

The stripping of the whole envelope has a significant effect on the planets radius. Due to the low density of gas, the presence of an envelope will result a significant increase of the planets’ sizes even if the envelope mass is only on a percent level of the total planet mass. Bare cores are thus clearly separated from object that retain a gaseous envelope, and a gap is observed in the distribution of planetary radii (Owen & Wu 2013; Lopez & Fortney 2013; Jin et al. 2014; Chen & Rogers 2016; Fulton & Petigura 2018).

The evaporation model is based on Jin et al. (2014). It takes into account contributions from X-ray and extreme-ultraviolet (XUV) irradiation. At the early stages, the evaporation is typically X-ray driven. We describe this regime using the energy-limited rate from Jackson et al. (2012) using the flux in the 1 to 20 Å range from Ribas et al. (2005) and assuming an efficiency factor ϵ = 0.1.

At later stages, the evaporation from EUV takes over. We also use the work of Ribas et al. (2005) to obtain the time-dependent EUV stellar luminosity for a Sun-like star. EUV evaporation can be divided into two sub-regimes (Murray-Clay et al. 2009). At low EUV fluxes, the same energy-limited approximation as for the X-ray flux is used. In this case, the escape flux is given by (79)

where FEUV is the EUV flux, Rbase the radius of the photoionisation base, calculated as in Murray-Clay et al. (2009), and ϵ = 0.3 is the heating efficiency, taken as in Murray-Clay et al. (2009).

On the other hand, energy-limited evaporation is not suitable when the EUV flux is high (> 104erg cm−2 s−1), as a substantial part of the heating is lost in cooling radiation. In this regime, we adopt the radiation-recombination-limited approximation of Murray-Clay et al. (2009). The mass loss rate is given by wind due to escape (80)

at the sonic point Rs, which is calculated the same way as Racc. Here is the isothermal sound speed of ionised gas with T = 104 K. The density can be related to the one at the ionisation base, where τ = 1, with (81)

The photoionisation base is located where there is equilibrium between photoionisations and recombination: (82)

with n0,base the density of neutrals at the base, 0 = 20 eV, , αrec= 2.7 × 10−13, and ρbase = n+,basemH.

The model also includes the effect of Roche lobe overflow. When solving the internal structure equations, there are sometimes solutions found in the detached and evolutionary phase where the radius is larger than the Hill sphere. This occurs in two situations: first, for close-in low-mass planets with a high envelope mass fraction. At the moment when the nebula dissipates (and thus the ambient pressure vanishes), and when the star starts to irradiate the planets directly (resulting inan increase of the temperature, see Fig. 3), these planets bloat. Second, giant planets that get very close to the star because of tidal spiral in (see Sect. 5.3) can also overflow their Roche lobe. In this case, we remove at each time step the part of the H/He envelope that is outside of the Hill sphere.

4.5 Initial conditions

The simulation begin with a predetermined number of embryos whose initial mass is Memb,0 = 10−2 M (approximately the mass of our Moon). They are randomly placed with an uniform probability in log a, where a is the semi-major axis, between rin and 40 au. The starting location zone is slightly more extended that in the previous studies, where the upper boundary was set to 20 au. Also, two embryos cannot be placed within 10 Hill radii from each other. It should be noted that for the simulations with largest initial number of embryos, 100, this represents an average spacing of 28 Hill radii.

The presence of a number of embryos right at the beginning of the simulations is a strong assumption we made because the model does not track the formation of the embryos themselves. This shortcoming of the model will be addressed in future evolutions of the model (Voelkel et al. 2020), where the evolution of the dust, pebbles and planetesimal and embryo formation is followed.

5 Dynamical evolution: orbital migration, N-body interaction, and tides

As the planet mass increases, it will generate a stronger perturbation in the density of the gas around the planet. This perturbation will cause the nebula to no longer be axis-symmetric, and as a consequence produces a torque back on the planet, leading to planetary migration. At the same time, convergent migration can result in capture in mean-motion resonances or orbital destabilisation. Hence migration and dynamical evolution must be performed together to capture all the effects.

5.1 Planetary migration

We include two types of migration, Type I for low mass planets embedded in the gas disc and Type II for planets massive enough to open a gap in the disc.

5.1.1 Type I migration

For Type Imigration, our model follows the approach of Coleman & Nelson (2014). This includes the torques formulation from Paardekooper et al. (2011), modified to consider that orbital eccentricity and inclinations attenuate the co-rotation torques (Bitsch & Kley 2010, 2011).

The total Type I torque on a planet, following Eqs. (50)–(53) of Paardekooper et al. (2011) and (15) of Coleman & Nelson (2014), is given by (83)

with

where ΓL, Γhs,baro, Γhs,ent, Γc,lin,baro and Γc,lin,baro are the Linblad torque, barotropic and entropy part of the horseshoe drag and linear corotation torque respectively. They are given by Eqs. (3)–(7) of Paardekooper et al. (2011). The function F governs saturation, while G and K provide the cutoff at high viscosity, and are given by Eqs. (22), (30) and (31) of Paardekooper et al. (2011).

The other factors in Eq. (83) account for the shape of the orbit. FL provides thereduction of the Lindblad torque for eccentric or inclined orbits following Cresswell & Nelson (2008), with (86)

and (87)

Here, ê = eh = e∕(Hr) and are the planet’s orbital eccentricity and inclination scaled by the disc’s aspect ratio h = Hr. Fe and Fi provide the reduction of the corotation torques due to eccentricity and inclination (Bitsch & Kley 2010). We use (88)

as suggested by Fendyke & Nelson (2014) for the reduction due to eccentricity and (89)

for the reduction due to inclination (Coleman & Nelson 2014).

Eccentricity and inclination damping time scales follow Cresswell & Nelson (2008), with (90)

and (91)

where (92)

is the characteristic time of evolution of density waves (Tanaka & Ward 2004).

5.1.2 Type II migration

The criterion to detect gap opening and switch migration to Type II is from Crida et al. (2006), (93)

with ν is the viscosity from Eq. (2).

Type II orbital migration follows the non-equilibrium approach from Dittkrist et al. (2014). Here, the planet follows the radial velocity of the gas, (94)

(Pringle 1981), but is limited if the planet’s mass is much larger than the local disc mass (the fully suppressed case, see Alexander & Armitage 2009). The radial velocity of the planet vplanet is given by (95)

For the larger planet masses, when the migration rate is constrained by the disc-to-planet mass ratio, this expression result in a similar behaviour as the formula obtained by Kanagawa et al. (2018), although it does not take into account the aspect ratio of the disc h.

For our migration scheme, we convert the radial velocity into a torque according to (96)

This prescription allows in principle planets in Type II to migrate outwards if the disc is decreting (Veras & Armitage 2004). However, in practice this mechanism is limited by the restriction to planets that are already at large distances or during the final moments of the disc, and limited by the small surface density (Dittkrist et al. 2014).

During type II migration, the eccentricity and inclination damping time scales are set to (97)

This relationship was selected because hydrodynamical simulations of migrating planets in this regime have shown that eccentricity and inclination damping act on time scales that are shorter than migration (Kley et al. 2004; Kley 2019).

5.1.3 Migration map

An example of the outcome of the whole migration scheme for one disc profile is provided in Fig. 10. The disc is the same as the example shown in Fig. 3 at 1 Myr; at this time the disc mass is 1.46× 10−02 M. Its outer radius is 123 au, so we cut the figure at 200 au since there is no migration outside this distance.

Migration is most efficient for intermediate mass planets, above about 10 M up to the transition to Type II migration (shown with the dashed black line on the migration map). The outward migration at large separation for the type II migration regime is due to the outward spreading of the gas disc. We also note two convergence zones for low- to mid-mass planets. These are due to opacity transitions (Lyra et al. 2010) or structures in the gas disc (Kretke & Lin 2012) such as the increase of the surface density close to the inner edge of the disc (Masset et al. 2006). These are the locations where, for a given planet mass, outward migration happens on the inner side and inward migration on the outer side. Hence, at this moment of evolution, planets with masses less than ≈ 8M cannot reach the inner edge of the disc by migration only. However, as time goes and gas becomes scarcer, the zones of outward migration (hence the convergence zones) shift towards lower planetary masses. Thus, by the end of the gas disc, planet with masses down to ≈ 2M could reach the inner edge of the disc.

thumbnail Fig. 10

Radial surface density profile (top panel), temperature profile (middle panel), and migration map (bottom) as function of the planet mass (assuming zero eccentricity and inclination), for the same disc presented in Fig. 3 at t = 1 Myr. The value plotted in the bottom panel is the relative migration rate 1∕τa = −vplanetaplanet; blue regions indicate inward migration, red regions outward migration. For both directions, the locations in bright colours are where migration is inefficient while dark tones indicate efficient migration. The dashed black line shows the boundary between type I (below) and type II (above) migration regimes.

5.2 N-body integration

Gravitational interactions between the protoplanets are now modelled with the mercury N-body code (Chambers 1999) using the hybrid method. Unlike the direct resolution of the equation of motion (as performed in A13), this use a symplectic integration scheme (see e.g. Sanz-Serna 1992, for a review). The basic principle is to use the solution of Hamilton’s equations, (98)

where x denotes the position coordinates, p the momentum coordinates, and (99)

is the Hamiltonian of the system, with Δxij = |xixj|. Here, the index i = 0 refers to the central star and M0 = M while the subsequent are the planet with Mi = Mplanet,i so that N is the number of planets in the system.

However, while has no analytical solution for N > 1, it is possible to split the Hamiltonian into several pieces, solving the simpler problems to finally combine them back so that a solution close to that of the original system. The Hamiltonian is divided into three components, so that , and

Here, HK represents the unperturbed Keplerian orbits of the planets about the central star, HS the kinetic energy of the star and HI the interactions between the planets. The separation into three different Hamiltonians (rather than two) is required because the scheme uses mixed-centre coordinates (also called ‘democratic heliocentric’): heliocentric positions and barycentric velocities. These coordinates are chosen so that HKHS, HI, unless two planets come close together.

The evolution of such a system by splitting is done using a second-order method, (103)

where the notation H...(τ) is used to represent the evolution under the given Hamiltonian for a step τ. For HI, this means that the planets receive a kick in velocity due to the interactions with the other bodies (except the central star). In our case, HI is extended to include additional forces representing the effect of the gas disc, see Sect. 5.2.1. The evolution under HS results in a shift τ∕(2M)∑ pi while the evolution under HK is a Keplerian motion around the central star for a period τ.

As we noted, the assumption that HI is small compared to HK is no longer valid when two bodies become close together. In that situation, the idea is to bring the interaction between the two close-by bodies into HK so that the interaction Hamiltonian remains small. This implies that HK is no longer analytically integrable during that period, but only for the orbits of the involved bodies. In practice, the orbits of the two close-bybodies are integrated with a conventional Bulirsch-Stoer method (Stoer & Bulirsch 1980) for the duration of the encounter. That algorithm is described in detail in Chambers (1999).

The symplectic integration scheme has a huge advantage in terms of computational requirements compared to a standard Bulirsch-Stoer method, as the interaction between the planets, the one part that is , is only computed once per step.

We do not use the N-body when there is only one protoplanet in a system as the solution is analytical. This happens either for populations with one embryo per system or in the unlikely case that only one planet survives in a planetary system with initially multiple embryos per system.

5.2.1 Additional forces

Migration and damping are included as additional forces in the N-body. The contributions from migration and eccentricity damping apply in the orbital plane and are split into tangential (θ) and radial (r) components, while the inclination damping acts on the vertical component (z), resulting in

with a denoting the additional accelerations, v the planet’s velocity along each direction. Here, vK = Ωr is the Kelperian velocity.

5.2.2 Collision detection

Collisions are detected when two planets come closer than a predetermined distance, which is the sum of their radii. When the closet approach is found inside to be during one of the substeps of the N-body, the minimum distance is retrieved by fitting a third-degree polynomial equation whose condition are set by the relative separation and their radial velocity at the beginning and end of the substep (similar to A13).

For planets with a significant and extended envelope (like during the attached phase), the assumption that planets have a unique radius which decides whether a collisions occurs or not is no very accurate, as the outcome is determined by gas dynamics inside the merging envelopes. As we do not have the full envelope structure in the N-body, we nevertheless remaining with a unique radius approach. In the attached phase, the envelope transitions smoothly to surrounding nebula. The outer radius, as provided by Eq. (44), is unsuitable for the detection of collisions, as it corresponds to very low gas densities. Thus, the radius used to detect collisions is computed assuming that the whole planet mass has the same density as its core. This is an approximation, but reflects that the gas density in the envelope is much higher close to the (solid) core surface. In the detached phase, we use the planetesimals’ capture radius Rcap; this is normally an overestimation of the effective collision radius, larger bodies needing to penetrate deeper down in the envelope to be captured. However, in this phase the envelope scale height is small compared to the radius except for the very short time directly after detachment,so the actual error is small.

5.2.3 Collision treatment

When a collision is detected, the following procedure is applied: the cores merge, the eventual envelope of the less massive body is deemed to be ejected, and the impact energy is added as a additional contribution to the luminosity for the structure calculation of the new body.

The merger of the cores will make that a part of the impact energy will already be taken into account consistently with the luminosity calculation described in Sect. 4.2; so the additional energy is calculated using (107)

where μ = Mtot,1Mtot,2∕(Mtot,1 + Mtot,2) is the reduced mass, and the indexes 1 and 2 refer to the quantities of the larger and smaller body respectively. vimp is the relative velocity at time of contact. Here (108)

is the centre-of-mass impact energy of two bodies with the total mass of the target and the core mass of the impactor colliding at their mutual escape velocity. Also, we restrict the supplementary energy to positive values. Negative value can arise if the bodies are collidingat below the mutual escape velocity, which is possible due to the drag by the gas disc or in the case of specific configuration, such as co-orbitals. However, the impact velocity is never quite lower than the mutual escape velocity, so that the error remains small.

The addition of the core mass and luminosity is performed via

where timpact is the time of the impact, τimpact = 104 yr is the time scale of release taken as in Broeg & Benz (2012). These two terms are added to the core accretion rate due to planetesimal accretion, and to the luminosity (Sect. 4.2) used in the internal structure calculation, respectively.

This impact model was tailored for the most common collisions that we find in our simulations. We highlight this by showing cumulative distrubutions of impactor-to-target mass ratio γ for different ranges of target masses in Fig. 11. At low masses, most target/impactor pairs are of similar masses, thus this source of growth cannot be neglected. In contrast, most collisions involving giant planets are with much smaller impactors (the red curve in Fig. 11). Our models neglect the envelope of the impactor, but there are only few collisions where this could provide significant source of mass.

thumbnail Fig. 11

Cumulative distribution of the impactor-to-target mass ratio γ = Mtot,2Mtot,1 for different target mass ranges (as provided in the legend). Data come from the 100-embryos population presented in Paper II.

5.3 Tidal evolution

During the evolution phase we include the inward migration of planets due to tides they raise onto the central star. In addition to planets thatare pushed inwards due to capture in mean-motion resonances, this gives another channel to obtain planets that are within the inner boundary of the gas disc. For the tidal migration rate, we compute the rate according to (111)

(Ferraz-Mello et al. 2008; Jackson et al. 2009; Benítez-Llambay et al. 2011), where Q = 106 is the stellar dissipation parameter. It is clear that this model for the tidal spiralling-in is strongly simplified. It will be improved in future work along the lines of, for example, Bolmont & Mathis (2016).

6 Terrestrial planet formation

We begin by studying whether the new generation of the Bern model with a higher initial number of embryos, but which still includes a statistical description of planetesimals, is capable of reproducing models of terrestrial planets that use purely N-body (e.g. Chambers 2001), that is where the planetesimals are represented as individual (test) particles. This test is crucial to assess whether we can reach our goal of having a formation model which is able to simulate the growth of planets with a very large mass range from about that of Mars, to brown dwarfs. This is in contrast with earlier generations of the Bern Model, where mainly more massive planets were at the focus (or more specifically, planets for which the giant impact phase after disc dissipation is not very important).

The formation of terrestrial planets does not have the same time constraint as for gas giants. In the case of planets with a significant H/He envelope, a sufficiently massive core must be formed before the dispersal of the gas disc, but this does not apply to terrestrial planets. Indeed, in the case of the Earth, cosmochemical evidences point to a formation time between a few tens Myr (Yin et al. 2002; Kleine et al. 2002) to roughly 100 Myr (Touboul et al. 2007; Allègre et al. 2008; Kleine et al. 2009). This is longer than the expected lifetime of the solar system’s nebula of 4 Myr (Wang et al. 2017) by about an order of magnitude or more. Hence the modelling of formation of planetary systems with terrestrial planets needs to span a longer time period for dynamical effects (i.e. the ‘late stage’) than for gas-dominated planets.

6.1 Setup

For this test cases, we performed a few modifications to our main model to mimic earlier work like Chambers (2001) and Raymond et al. (2005). Orbital migration has been disabled; as for the envelope structure calculation and the evolution phase, here, all planets are treated as purely rocky. We adopt an initial surface density profile close the minimum-mass solar nebula (MMSN; Weidenschilling 1977; Hayashi 1981), with a reference surface density of Σ0,s= 7.1 g cm−2 at r0 = 1 au, but truncated at2 au, as we are primarily interested in the inner planets. This also helps to determine more precisely which fraction of the planetesimal disc has been accreted by the terrestrial planets during their formation. This gives a solids mass of 3.67 M. The initial number of embryos is selected to have a similar spacing as the two populations presented Emsenhuber et al. (2021, Paper II) with most embryos per system, which means that we have initially 23 (correspond to 50 in Paper II) and 46 (corresponding to 100) lunar-mass (0.01M) embryos. In addition to that, we perform one run with 9 embryos initially in Sect. 6.3, which corresponds to 20 embryos in Paper II.

It should be noted that the model lacks the ‘dynamical friction’ obtained in N-body simulations with a large number of small bodies (O’Brien et al. 2006; Raymond et al. 2006) because we do not include the effect of the damping of eccentricities and inclinations of the embryos by the planetesimals. However, after all material has been accreted onto the planets, the remainder of the formation process is similar to pure N-body simulations of terrestrial planet accretion, as all the mass is now contained in bodies that are directly followed by the N-body.

For some simulations, we include Jupiter and Saturn to determine the effects they have on the formation of the inner planets. In that case, Jupiter and Saturn are on their present-day orbits, but they are rotated so that their invariant plane coincides with that of the disc (as in Chambers 2001, 2013; Emsenhuber et al. 2020). We do not model the formation of these planets, because they form over a period that is much shorter than the terrestrial planets.

To obtain a better overview of the influence of the parameters we are studying, and to reduce (and better understand) the stochastic effects of N-body interactions, we perform 10 simulations for each combination of parameters (initial number of embryos and presence of the outer planets). The only differences between the 10 simulations are the initial position of the terrestrial planet embryos. For the 10 simulations, we consider the average outcomes as being representative (e.g. Fig. 12).

The simulations starts with a gas disc, which lives for roughly 4.4 Myr. Its only effect however is to damp the eccentricities and inclinations of the planetesimals. Planetesimals accretion continues after the dispersal of the gas disc. As the planets do not have envelopes, we perform only the formation stage of the calculation. However, the duration of that stage has been extended to 400 Myr to account for the much longer time needed for the solar system’s terrestrial planets to converge.

thumbnail Fig. 12

Average over the 10 simulations for each set of parameters. The top panels show the masses of solids (excepting the outer giant planets in the relevant cases) in the protoplanetary disc, that is still in planetesimals (solid lines), accreted by the embryos (dashed lines) and ejected (dotted lines). The dashed black line denotes the total mass of solids in each simulation. The bottom panel shows the number of embryos that remain. The ‘J/S’ in the legend refers to Jupiter and Saturn.

6.2 Gravitational interactions

If the embryos remain at their initial locations during the whole formation process, then they grow to their isolation mass (Lissauer 1987). In our model, we obtain this behaviour if we artificially remove the N-body interactions, unless the feeding zones of two adjacent embryos overlap at some point, in which case the masses become slightly lower. When using this mode, the runs starting with 46 embryos have accreted roughly half of the disc’s mass onto the embryos by about 4 Myr (the time at which the gas disc disperses) and accrete very slowly thereafter. For the runs starting with 23 embryos, only a quarter of the mass ends in the embryosby 4 Myr.

For the other parameter sets (all with gravitational interactions), Fig. 12 provides the averaged results over 10 simulations, for the masses of solids and the number of embryos. The story is quite different when N-body interactions are included. We see for instance that in the case with 46 embryos and no outer giant planets, nearly all the planetesimals have been accreted onto the embryos. For the case with 23 embryos initially and no outer giant planets, more than half of the planetesimals end up accreted.

There are two aspects we point out here. First, in the figure, the planetesimal mass accreted by embryos that have been later ejected is accounted as accreted. Second, our planetesimal model does not include redistribution of material by interactions with the embryos. For instance, in their less realistic setup where embryos only populate a limited orbital distance range in the disc, Levison et al. (2010) found that planetesimals can be redistributed to locations outside of the embryos’ feeding zone rather than be accreted. However, when they add the mechanism that embryos can reside in all parts of the disc (which is more realistic, Levison et al. 2010) no gap in the planetesimal disc opens, as embryos mutually scatter planetesimal into their vicinity and accrete them eventually. This leads to an efficient formation of massive planets. Thus, feeding zones overlap in these simulations, therefore the effect of planetesimal redistribution should play little role in our case as there are very few locations in the disc where planetesimals would not be accreted by the local embryo and/or scattered back into the feeding zone of other embryos.

thumbnail Fig. 13

Comparison of formation tracks for a system with a MMSN-like surface density of planetesimals, and with 9 embryos (left), 23 (centre), and 46 (right), and no outer giant planets. Each line represents one embryo. Top panels: mass versus semi-major axis; embryos start at the bottom and move upwards as they grow. The final positions of the remaining planets are shown by dots. The dashed black line denotes the isolation mass (Lissauer 1987). Middle panels: mass versus time; sudden increases in mass are due to embryo-embryo collisions. Bottom panels: semi-major axis versus time.

6.3 Interactions lead to more massive planets

To understand how the embryo-embryo interactions lead to a quasi-complete accretion of the planetesimals disc, we show the formation tracks for one particular system with a varying number of embryos in Fig. 13.

We can easily observe that the larger the number of embryos, the more and the sooner they start to move around. In the system with only 9 embryos, they basically remain where they started and grow slightly above their isolation mass. For the other two simulations, however, the local isolation mass is sufficient to trigger significant embryo-embryo interactions that will change their positions in the disc. This in turn enables them to accrete from regions that would otherwise inaccessible, which creates a positive feedback since more massive planets will result in yetmore interactions. This feedback only ends when nearly all planetesimals have been accreted onto the embryos.

Thus, closer packed embryos lead to enhanced stirring of their eccentricities, which has two consequences: the increase of the feeding zonesize because of radial excursion for eccentric orbits, and collisions between embryos. Embryos having a greater eccentricity can sample a broader region of the disc, thus grow to a larger mass before depleting the disc. Collisions with other embryos arecapable to bring material from more distant regions of the disc that would otherwise not be accessible to one embryo. At the end, we arrive at a result that is maybe counter-intuitive at first: the larger the number of embryos, the less planets remain. We observe this for instance in the bottom panel of Fig. 12.

6.4 Time needed for formation

We find a similar pattern for the timing at which interactions start in the two simulations with the higher number of embryos of Fig. 13 (23 and 46 embryos). In the early phase (a few 105 yr), no dynamical interactions occur, because the embryos need to reach a certain mass before the eccentricities can be significantly excited. Then, the first embryos to show an increased eccentricity are located at ~ 0.3 au, and then this propagates both inwards and outwards. In the inner part of the system, collisions happen rather rapidly so that the system has essentially obtained its final configuration by several Myr.

On the other hand, in the outer region we observe that embryos remain on eccentric orbits for a certain amount of time before suffering from collisions. It takes more than 10 Myr for the planets located at about 1au to reach their final mass. In the even more distant regions, it takes even longer, and we see the phase with several embryos on eccentric orbits remaining for more than 100 Myr. Such a growth wave travelling from the inside to the outside is expected, as the growth process scales with the local Keplerian frequency.

Therefore, our choice of the integration time dictates the location where and how accurately the model can follow the formation of the terrestrial planets. With our choice of an integration time limited to 20 Myr for the formation phase, the model can only track most of the giant impact stage inside of roughly 1 au for systems that have a MMSN-like surface density of solids. Even within 1 au, the giant impact stage is not entirely finished within our set time frame, as it can be see in the innermost planet by about 300 Myr in the bottom right panel of Fig. 13. Nevertheless, these events remain rare. Locations further away or systems with a lower amount of solids (as formation is slower for less massive systems, Kokubo et al. 2006; Dawson et al. 2015) will, however, not have reached a final state by end of the formation stage at 20 Myr.

6.5 With outer giant planets

As the final stage of terrestrial planet formation (the giant impact stage) takes longer than formation of the giant planets, we also want to consider the effects of their presence on terrestrial planet formation. Here we perform the same simulations again, each time with the addition of two outer giant planets that represent Jupiter and Saturn.

To provide a better comparison point between the two cases, we provide in Figs. 14 and 15 several snapshots of the simulations. One general consequence at earlier times is that there is slower growth for the embryos beyond 1.5 au. We see in the two top rows of Fig. 14 that the outermost embryos remain smaller in the runs with outer giant planets. Also, their eccentricities have already increased in the first snapshot, while this is not the case at all for the runs without giant planets. The underlying cause is stirring of planetesimal’s eccentricity and inclination by the giant planets; this heavily reduces the collision probability with the low-mass protoplanets (Inaba et al. 2001) and hence the accretion rate.

A consequence of the longer timescales of accretion in the outer part of the disc is the state at the moment of the dispersal of the gas disc. In the runs with giant planets, a larger percentage of the planetesimals remains unaccreted at the momentthe gas disperses. In addition, after that point, there is no longer gas present to counterbalance the effects of the stirring by the giant planets. This means that after a short moment, the planetesimals will reach eccentricities of the order of unity and will be ejected. This can be observed in Fig. 12, where we see that up to a quarter of the original initial mass is ejected from the planetesimals disc. The final eccentricities of the terrestrial bodies are similar in both cases (Fig. 14), as the inner region is subject to the self-stirring while in the outer region, excitation by the outer planets makes up from a weaker self-stirring as the masses are lower.

Thus, the outer giant planets will limit and delay growth of the terrestrial planets in the outer region. The number of objects is a bit higher than the one obtained by pure N-body simulations of terrestrial planet formation, but we are using a somewhat smaller initial surface density profile compared to, for example, Raymond et al. (2006), which prevents the accretion into a lower number of higher mass bodies (Kokubo et al. 2006).

thumbnail Fig. 14

Stacked eccentricity versus distance snapshots of 10 simulations with each 46 embryos initially. The left column shows the runs with outer giant planets whereas the right column has no outer giant planets. In each column, the 10 systems are represented with a different colour for each one. The bodies are shown by points whose sizes are proportional to their physical ones. Black crosses show the solar system planets.

thumbnail Fig. 15

Same as Fig. 14, but showing mass versus distance stacked diagram of 10 simulations with each 46 embryos initially. As before, the black crosses show the solar system planets.

6.6 Summary for terrestrial planets

To summarise, we have just seen that as long as the separation between the embryos is sufficiently small that dynamical interactions are triggered before the embryos reach their local isolation mass, the model is capable of reproducing the main features of the formation of terrestrial planets in good agreement with pure N-body models. This is due to embryo-embryo interactions being able to increase the eccentricities, so that the embryos can move out of their original locations, and almost entirely depletes the planetesimals.

An integration period (for the formation stage) longer than the lifetime of the protoplanetary disc is necessary to follow the giant impact phase. The time required for the bodies to obtain their final characteristics increase with distance (as shown here) and with decreasing initial amount of solids (e.g. Kokubo et al. 2006; Dawson et al. 2015). The limitation of the formation stage to 20 Myr (Sect. 2.2) permits to capture all the accretion of planetesimals (provided there are enough embryos initially) and most of the dynamical interactions of Earth-mass and larger planets forming via giant impacts up to roughly 1 au, and sub-Earths planets in the first few tenths of an au (corresponding to periods of roughly 100 days). For the population syntheses in Paper II, we estimate from tracking major changes of the planets’ orbits, that for orbital distances of ≲ 1 au around 90% of the major instabilities should have been captured when integrating the systems for 20 Myr.

The integration time needed to capture most instabilities within a given orbital distance range is a function of the architecture of the planetary systems that results from the previous growth stages. If the growth and migration during the presence of the gas disc leads for example to very closely packed systems of massive planets, instabilities will often occur shortly after disc dispersal. On the other hand, if at gas disc dispersal only low-mass, widely-space planets are present, they will first have to grow further via accretion of remaining planetesimals and embryos – which can take a very long time – to eventually (or also never) become unstable. For larger planet masses, gravitational interactions can extend further: even on distant orbits, massive planets can destabilise the system as noted by Bitsch et al. (2020) and Matsumura et al. (2021). This could explain why Izidoro et al. (2021) find that by 20 Myr only a fraction of the instabilities between the planets have happened in their setup, whereas Mulders et al. (2020) on the contrary find that increasing the integration time from 10 Myr to 100 Myr only leads to minor further evolution in their simulations. For systems lacking outer planets, Izidoro et al. (2021) also found a convergence after ~30 Myr (their Model III). This is in better agreement to the analysis for planets with a < 1au done in Mulders et al. (2020).

The purpose of the model here is to obtain planetary systems that can be compared with observations at the population level. For this, it is important to see that the region where long-term growth will be most important (distant low-mass planets) represents at the same time the parameter space currently not accessible to most detection techniques of extrasolar planets (radial velocity, transits, and direct imaging). This should minimise the impact of this limitation. We acknowledge, however, that generally speaking, not all dynamical interactions will have taken place by the end of the integration time of 20 Myr in the model. While the later evolution should not be substantial enough to strongly affect the statistical results in the inner systems, this limitations must be critically kept in mind when comparing for example to microlensing surveys (e.g. Suzuki et al. 2018) that probe more affected regions.

Nevertheless, we conclude that the new generation of syntheses can be used to describe in a much more comprehensive way planetary sub-populations ranging from sub-Earths to super-Jupiters.

7 Giant planets

The formation of giant planets is quite different. Cores must form before the dispersal of the gas disc so that they can undergo runaway gasaccretion, and since we have massive cores in a gas disc, migration is efficient. To gain an understanding of the interplay of accretion and migration, we here show some illustrative cases with a single embryo per disc. For this case, we use the model without modifications, but the N-body is not used. The following examples are taken from the single-embryo population of Paper II. Simulations parameters are the same as provided in Table 1, except for disc masses and surface density (both gas and planetesimals), inner edge, characteristic radius, and external photoevaporation rate of the gas disc. In the following simulations, the inner radius has negligible effect on the final outcome, as we do not study close-in planets, and so we do not mention it. The characteristic radius rcut,g of the gas disc is set as (112)

(see Paper II for the motivation). We provide the remaining two parameters, the initial masses of the gas and planetesimals discs in the following.

thumbnail Fig. 16

Formation and evolution tracks of four giant planets with final masses in the 1∕3 to 2 M range in discs with a single embryo. Top panels: formation tracks with total mass Mtot versus distance (time goes towards the top) and total mass Mtot (solid lines) and core mass Mcore (dashed lines) versus time. Three panels on the bottom row: time dependence of the outer luminosity Ltot (bottom left), the distance (bottom centre) and the total radius Rtot (bottom right). For all panels except for the mass versus time (top right), the line styles denote the phase: dashed lines for the attached phase, solid line for the detached phase during formation and dash-dotted lines for the evolution stage. Line widths denote the migration regime, with tick lines for Type I and think lines for Type II. The legend in the top right panel shows the gas (in Solar masses) and planetesimals (in Earth mass) disc masses.

7.1 Formation and evolution of Jupiter-mass planets

We show in Fig. 16 the formation tracks of a few synthetic giant planets whose masses are in the 100–500 M range and have a wide range of final positions. Due the inclusion of migration in the model, we observe that the final position of these planets is closer-in that the initial location of the embryo: all the embryos start beyond 10 au, with one close to 30 au, while all the planets end up inside 10 au.

During the initial stage, both accretion and migration are slow, but accretion is still faster. As the planets grow, migration becomes more efficient; we observe that most of the migration occurs while the planets are close to the transition to gas giants, with masses between 20 and 50 M. The innermost planet shows a strong inward migration at this stage, but this is due to limited accretion while migration remains at the same rate. Once the planets undergo the runaway accretion of gas and switch to type II migration, accretion is strong, and they experience limited migration. This leads again to near-vertical tracks. The two changes (from an attached to a detached envelope and from type I to type II migration) happen in the same period, not always in the same order. In one case (the inner most planet shown in red), the change of the migration regime occurs first, while in the three other cases it is the reverse. Once the migration regime changes to type II, the rate slows down (bottom centre panel of Fig. 17) but the accretion remains mostly constant. Thus, accretion dominates at the onset of this stage, but this reverses at the end. In contrast, Mordasini et al. (2009a) used the equilibrium values of the radial gas flow for both gas accretion and migration. Thus, the slope of detached planet migrating with the planet-dominated case of the Type II regime exhibited a common slope in the mass-distance diagram. It should also be noted that for planets inside roughly 1 au, it happens that the criterion limiting the gas accretion rate changes to the mass in the feeding zone which leads to a reduction of the rate at the end of the formation. It can also be noted that our model allows for the growth of embryos at large separation (up to about 30 au, unlike the work of Johansen & Bitsch (2019). The difference is mainly related to the planetesimals size. Smaller planetesimals have lower eccentricities and inclinations because of more efficient damping by the disc gas, and in addition, a larger capture probability by the planets for a given surface density because of the more strongly drag-enhanced capture radius for small planetesimals. This results in a larger accretion rate of solids, which enables planets to sufficiently grow to undergo runaway gas accretion before the dispersal of the gas disc also in the outer parts of the disc.

The formation of Jupiter-like planets with migration and planetesimals accretion follows a different pattern in the one-planet-per-disc approximation studied here than what was found by some other models using the in situ (and one-embryo-per-disc) approximation. For the latter, the favoured scenario is that a core between 10 and 20 M forms early (less than 105 yr) and undergoes runaway gas accretion only close to the dispersal of the gas disc (Pollack et al. 1996; Alibert et al. 2018). The slow accretion of planetesimals, resulting in a steady luminosity, is able to prevent runaway gas accretion during the intermediate stage. This intermediate stage is the problematic part when migration isincluded; the reason being that migration is most efficient for planets that are between 10 and 50 M (see Sect. 5.1.3 and Fig. 10). Hybrid pebble-planetesimals (Alibert et al. 2018) or pure pebble (Bitsch et al. 2019) accretion models can account for the migration during the intermediate phase, as the cores are able to form at larger separation, provided that far out, bodies emerge early and massive enough to be able to efficiently capture pebbles.

The simulations presented here show a situation with fast migration where no intermediate stage is possible, because the planets would otherwise end up at the inner edge of the gas without the opportunity to undergo runaway gas accretion. This means that the cores must form just at the time to undergo runaway gas accretion. The usual picture of the formation of Jupiter-mass planets in our model is then more similar that what was found by Alibert et al. (2005a), with an almost nonexistent intermediate phase (left panel on the second row of Fig. 17). As the accretion time scale are longer at large separation, the embryos will accrete their mass over a longer period. At the same time, the inward migration experienced by the protoplanets means that their feeding zone is not depleted as in the in situ formation scenario.

In contrast, for the multi-embryos simulations that we present in Paper II (see also Sect. 8.1), migration can be significantly altered by mean-motion resonances chains. In that case, the torque acting on one planet must be spread over all the bodies, meaning that the planet with the largest specific torque will migrate slower than it would were it not in a resonance chain. This provides a way to obtain an intermediate stage and less overall efficiency of migration, as we show in that work. This effect leaves open the possibility to have an intermediate stage for the formation of giant planets, as obtained in Alibert et al. (2018). Thus, once multiplicity is included, the simulations here become more similar to Jupiter formation models like Alibert et al. (2018).

thumbnail Fig. 17

Formation and evolution tracks of two groups of four giant planets with final masses between 2 and 10 M in discs with a single embryo. Panel and line descriptions are the same as Fig. 16.

thumbnail Fig. 18

Formation and evolution tracks of one giant planet that ends up being accreted by the star during the evolution stage due to tidal migration. Top left: total mass Mtot versus distance; top right: total mass Mtot (black) and core mass Mcore (red) versus time; bottom left: total luminosity Ltot (black) and the bloating contribution Lbloat (red) versus time; bottom centre: distance versus time; bottom right: total radius Rtot (black) and outer (attached phase) or Hill (detached phase) radius (red) versus time.

7.2 More massive planets

Figure 17 shows the formation tracks of planets that are in the 2 to 10 M range. Compared to the planets previously discussed, these ones show a greater range of initial locations (from 6 to 40 au) and overall effect of migration. The planet shown in orange is the quickest to accrete a massive core and undergo runaway gas accretion, due to both the more massive disc and the inner location. The latter is made possible due to the disc’s mass. This is also the one to migrate the least before reaching 10 M because (1) the fast formation limits the effect of migration and (2) enters a convergence zone (see Fig. 10 and the discussion in Sect. 5.1.3). As the boundary of convergence zone moves inward (Lyra et al. 2010; Dittkrist et al. 2014)and to lower planetary masses over time, the planet shown in red will encounter the convergence zone at a different location, which will not affect the planet as much.

Unlike the Jupiter-mass planets, all the ones of this group first switch to type II migration before going to the detached phase. This is seen on the top right panel of Fig. 17, where the tracks become dashed and thin during a brief section. Theslope break that was discussed in for the Jupiter-mass is stronger for the two innermost planets. Comparing the time evolution of the two, it can be noted that the migration rate remains mostly constant while in the type II regime, while the accretion rate decreases. Concerning the radius and luminosity, we observe that all the planets show a similar behaviour even with the difference in the final location.

7.3 Giant planets ending in the star by tidal migration

As an illustration how close-in planets are affected by the newly added physical processes during evolution, we finally discuss the formation and evolution of a close-in giant planet. These will raise tides onto the star, which will result in tidal migration. The consequenceis that the planet can be accreted by the star at some point during its evolution. We show such a case in Fig. 18. The formation stage looks quite similar to the previous example, with the difference that the planet ends at a close-in location, 0.04 au. The radius shrinks already before the planet goes to the detached phase, because it experiences a strong inward migration at the same time (as it can be seen in the lower right panel of Fig. 18). As the planet migrates inward, the Hill radius shrinks. Once the detached phase begins, the Hill radius continues th shrink as further inward migration continues.

As this planet is close to the star (0.04 au), the evolution stage is different from the case shown previously. The luminosity increases over time time, the envelope gradually expands and looses mass due to atmospheric escape and the planet migrates further inward due to the tides raised onto the star. The migration rate increases over time due to its strong dependence on the distance between the planet and the star (see Eq. (111)). To determine the reason for the luminosity increase, we print alongside the total value, the contribution from bloating (Eq. (63)). We see that from late in the formation stage until the end, this contributes to nearly all the planet’s luminosity. And as it goes with the stellar flux, it increases at late times due to tidal migration. The luminosity increase in turns leads to an expansion of the envelope, which increases the loss rate by atmospheric escape. But rather than this being the main cause of gas loss, we see that the bulk of the envelope is removed because it overflows the Hill sphere. This occurs suddenly at the end of the planet’s life, once the outer radius gets larger than the Hill sphere. Only a bare core remains, which get accreted by the star shortly thereafter.

7.4 Summary for giant planets

The formation and evolution of giant planets involves multiple concurrent processes. Migration being most efficient during the onset of the gas runaway accretion, this phase must occur in a relatively short time for the planets to not end up at the inner edge of the disc, in the absence of another planet to prevent migration. This also means that the cores must form late (i.e. shortly before the dispersal of the disc) to prevent a massive envelope from being accreted. Close-in planets will experience addition effects during their evolution, such as atmospheric escape and inward tidal migration that can lead to accretion by the star. In the latter case, it is possible for Hill sphere overflow to cause the loss of most of the envelope.

8 Individual systems

After discussing formation pathways of terrestrial under idealised conditions, and of single giant planets, we finally show results obtained with the full model. Using many embryos per system, the model is able to produce a very large variety of planetary systems. These range from terrestrial planets (as we saw in the previous section) to giant planets. We first provide two examples of the temporal emergence of planetary systems and then show the variety of the final architecture of 23 systems.

8.1 Two examples of the temporal evolution of the emergence of planetary systems

Figures 19 and 20 show the formation of a planetary system in a protoplanetary disc with a low initial content of solids (System 30 in NG76, see Paper II). The initial disc gas mass is 0.023 M while [Fe/H] is −0.13, corresponding to a dust to gas ratio of 0.011. This results in a low initial solid content of 65.1 M in the disc of planetesimals. The disc is seeded with 100 lunar mass embryos at t = 0, distributed uniformly in the logarithm of the semi-major axis inside of 40 au (see Paper II for more details on the initial conditions).

8.1.1 Low initial solid content

Many aspects of the emergence of the planetary systems can be understood with the comparison of the timescales of growth and migration, and the consequences of (large-scale) dynamical instabilities caused by the gravitational interactions of protoplanets. Therefore we colour code in Fig. 19 showing the temporal evolution of the system in the aM plane the tracks of the planets by the ratio |τmigτgrow| = |d lnm∕d lna|. Regarding the timescales, it is of fundamental importance that the oligarchic planetesimal accretion timescale increases with increasing planet mass (e.g. Thommes et al. 2003), whereas the orbital migration timescale in the Type I regime decreases with planet mass (e.g. Ward 1989).

At the beginning (105 yr, top left panel of Fig. 19), the quasi in-situ accretion of planetesimals present in the initial)feeding zone of the embryosis the dominating process. Migration occurs at these very low masses on a much longer timescale, leading to nearly vertical upward tracks. We note that the model does not include any artificial reduction factors of Type I migration. The specific distance dependency of the mass to which the protoplanets have grown by 105 yr is given by the following interplay of growth timescale as a function of orbital distance and the local availability of solids: from the innermost embryo at about 0.03 au to the one at about 0.6 au, the protoplanets have already grown to the local planetesimal isolation mass (Lissauer 1987). Given the planetesimal surface density scaling with r−3∕2, the isolation mass increases with orbital distance. As can be seen in panel b of Fig. 20 which shows the mean planetesimal surface density in the feeding zone of the planets, at 105 yr, the surface density is already strongly depleted in the inner parts of the disc. Between the local maximum at 0.7 au and the water iceline at 2.7 au, the mass is in contrast decreasing with distance because protoplanets further out grow slower. The next feature is a sharp increase of the protoplanets’ mass by about a factor 2 across the water iceline because of the increase of the solid surface density. One protoplanets grows in the transition zone, giving it an intermediate mass. Outside of the iceline, the masses decrease again with distance because of the longer growth timescales. For the protoplanets in the inner part that have already reached the isolation mass, the growth is temporally stalled. Because of the very low (isolation) masses of these protoplanets, orbital migration is nevertheless negligible.

At the very beginning, all protoplanets grow as if they were the only bodies in the disc, not feeling the influence of the other protoplanets. With increasing mass, the interaction (either directly via N-body interactions) or indirectly via resonant migration, become important. By × 105 yr, the first dynamical interactions have started among some of the more massive protoplanets, which is visible as a ‘jitter’ in some tracks, and two collisions, which are shown by two open grey circles.

At 1 Myr (top right panel in Fig. 19), inside of the iceline, the character of growth has changed from planetesimal-dominated, to some first growth via giant impacts (embryo-embryo collisions) for some protoplanets or stalled growth for others. As can be seen in panel a of Fig. 20 which shows the semi-major axis of the (proto)planets as a function of time colour coding the mass, about 10 further giant impacts have occurred. This has allowed the protoplanets in the inner disc to grow beyond the local isolation mass. As visible in panel b of Fig. 19, at 1 Myr, the planetesimal disc is now depleted out to about 1.3 au, and as time proceeds, the depletion moves even further out. We thus see a growth wave moving outwards (Thommes et al. 2003). All solid mass has been transferred into the embryos in this part, and their mutual interaction (giant impacts) governs the further mass growth. This implies that the accretion of planetesimals is only important at the early phases when the planets grow mostly in-situ. In the outer disc beyond the iceline, growth in contrast still proceeds mainly via planetesimal accretion, as there is a larger mass reservoir available. Between 2 and 4 au, a group of about 10 protoplantes with a mass of about 1 M has formed, meaning that the most massive planets are now found further out than before. These protoplanets originate from (just) beyond the iceline. The colours of the lines in Fig. 19 show that migration is still much slower than accretion for these planets at 1 Myr, but some slight inward migration is now occurring, causing the tracks to bend inwards. This applies also to the inner disc, where horizontal tracks are visible. They result from the depletion of the planetesimal disc, and the fact that the cores are of such a low mass that virtually no gas accretion is possible.

At 3 Myr (bottom left panel of Fig. 19), in the inner disc, the dominant effect is further growth via giant impacts. About 25 protoplanets with masses between the one of Mars and Earth are now present. In the outer disc, beyond the iceline, the aforementioned group of the about 10 most massive protoplanets has grown further now reaching a maximum mass of 3 M, and has also migrated further inward. As these planets migrate into zones that have been previously depleted by inner planets (in particular inside of the iceline), planetesimal accretion is quickly stalled. This means that planetesimal accretion for migrating planets is usually limited in low-mass multiple systems like the one present here. This means that a possible shepherding effect (Tanaka & Ida 1999) that we do not include in the model should not affect the outcome very much, except for a transition phase where τmigτacc for some planets. This phase can be seen for the outer group from the cyan line colours. As can be seen in panel a of Fig. 20, the planets capture each other in very large resonant convoys and migrate together (e.g. Cresswell & Nelson 2008; Alibert et al. 2013). In this configuration, outer more massive planets push inner smaller planets.

As visible in Fig. 20 by the small black circles, many giant impacts seem to occur in groups (i.e. at similar moments in time in fast sequence): a first group occurs at about 3 Myr, a next one at 4 Myr, and again one at the moment when the disc inside of about 2 au becomes free of gas. This is visible in panel c of Fig. 20, which colour codes the gas surface density at the planets’ position. This moment corresponds to the opening of the inner hole in the gas disc because of internal photoevaporation (cf. Fig. 3). At this moment, the damping effect of the gas vanishes, allowing orbit crossings and collisions (e.g. Ida & Lin 2010). The outer gas disc dissipates a bit later, at 5.1 Myr, shown by the vertical line in panel c of Fig. 20. After the dissipation of the disc, only 3 more giant impacts occur in this system to 20 Myr.

The temporal evolution of the eccentricities is shown in Fig. 21. The colours show the planet mass. For clarity, only planets with a mass of at least 0.1 M were included. One can clearly see the increase of the typical values of the eccentricities near the time the gas discdissipates at around 5 Myr. Before, typical values of the eccentricities are of the order of 10−3 to a few 10−2. After disc dissipation, they increase to values between about 0.02 to 0.2. Such values are expected from the increase of the velocity dispersion of the orbits until they are comparable to the escape velocity from their surfaces resulting from close encounters, once the damping by the gas is gone (Goldreich et al. 2004). One also sees that more massive bodies tend to be less eccentric, likely a consequenceof energy equipartition.

In our model, dynamical friction by residual planetesimals is neglected. This would reduce the eccentricities and inclination of the protoplanets. This implies that our model tends to overestimate the eccentricities and inclinations of lower mass planets for which dynamical friction by planetesimals would play a role.

The general sequence of solid growth that is first dominated by the near in-situ accretion of planetesimals followed by the second phase of growth via giant impacts is well visible in panel d of Fig. 20. It shows the mass of the protoplanets as a function of time. The line colours show the semi-major axis. We note how the transition between the two regimes occurs the later the more distant a planet is. At the largest orbital distances where embryos were inserted into the disc (maximum starting distance is 40 au), nearly no growth at all has occurred during the simulated period. As described in Sect. 5.2.3, numerically speaking, we add the mass of the impactor in a giant impact over a timescale of 104 yr to the target. This is the reason why the vertical steps in the curves corresponding to giant impacts (indicated with the black circles) are not strictly vertical. This is visible particularly at the early ages.

The bottom right panel of Fig. 19 shows the system at 20 Myr, which corresponds to the time where we stop the N-body integration and planetesimal accretion. Between 3 and 20 Myr, numerous giant impacts have reduced the number of planets and destroyed the mean motion resonances (see also Fig. 20). The inner system now contains 8 roughly Earth-mass planets, exhibiting a certain inter-system similarity of the mass scale (Millholland et al. 2017) with an increase towards the exterior (Weiss et al. 2018). At 0.7 au, there is a sudden increase in the typical mass, corresponding to the transition from volatile-poor planets that have formed inside of the iceline, to very volatile-rich planets originating from beyond the iceline. Compared to the original location of the iceline at 2.7 au, there was thus an inward shift in this transition by about 2 au because of orbital migration.

In the end, the planetesimal disc is depleted out to about 5 au. Outside, about 35 M remain in the form of planetesimals. This corresponds to a fraction of about 46% of the initial planetesimal mass that was converted into planets. Since we follow the accretion for only 20 Myr, this remaining planetesimal mass must be considered an upper limit for the actual mass of remaining planetesimals, as over longer timescales, the distant protoplanets would continue to accrete. However, since the accretion timescales at several 10 au in the absence of eccentricity damping (corresponding to orderly growth) become extremely long (Ida & Lin 2004a), at least some part of these planetesimals could remain to eventually form a debris disc, in analogy to the Kuiper belt beyond the orbit of Neptune in the Solar System.

thumbnail Fig. 19

Example of the formation of a planetary system from initially 100 lunar-mass embryos in low gas mass initial mass 0.023 M), low metallicity ([Fe/H] = −0.13) disc. The initial mass of planetesimals is 65.1 M. Four moments in time (in years) are shown. Lines shows the growth tracks in the semi-major axis-mass plane. Black points show (proto)planets existing at a given epoch. Grey open circles show the last position of protoplanets that were accreted by another more massive body in a giant impact. The colours of the lines are |τmigτgrow| = |d lnm∕d lna|.

thumbnail Fig. 20

Same system as in Fig. 19, but now showing the semi-major axes a of the planets as a function of time, colour coding in panel a the planets’ mass, in (b) the planetesimal surface density in the planets’ feeding zone, and in (c) the local gas surface density. Here, the vertical line indicates the moment of gas disc dissipation. Panel d: mass as a function of time, colour coding the semi-major axis. Small black circles indicate giant impacts, by showing the position or mass of the target (the more massive collision partner) at the moment of the impact.

thumbnail Fig. 21

Temporal evolution of the eccentricities of the planets of the system emerging in the low-mass disc shown in Fig. 19. Colours indicate the planet mass. For better visibility, only planets more massive than 0.1 M are shown. The curves are running averages such that one sees more clearly the mean values instead of rapid variations of the eccentricities. The thick black line is the mass of the gas disc relative to the value at 105 yr, which is in turn very similar to the initial value. The increase of the eccentricities at around 5 Myr when the gas disc dissipates is visible.

8.1.2 High initial solid content

The second system we consider is System 852 in NG76 (Paper II). The initial conditions are here a disc mass of 0.066 M and a metallicity of [Fe/H] = 0.23. This leads to an initial planetesimal mass of 432 M, 6.6 times as much as in the first example. As in the previous case, 100 lunar mass embryos are put into the disc at the beginning, uniform in log of the semi-major axis out to an orbital distance of 40 au. The evolution in the aM plane is shown in Fig. 22. The semi-major axis and mass as function of time is shown in Fig. 23.

In the top left panel of Fig. 22 we see that at 105 yr, the basic picture regarding the (relative) mass of the protoplanets as a function of orbital distances is analogous to the one in the low-mass disc at the same time. In absolute terms, the planet masses are, however, about one order of magnitude larger. As can be seen in panel b of Fig. 23, the planetesimal disc is already strongly depleted out to about 1 au. Some giant impacts have also already occurred in the inner disc. This fast development away from the initial conditions is a sign that the early phase of solid growth (from dust to embryos) should be treated more explicitly (e.g. Voelkel et al. 2021).

The situation at 0.5 Myr is already quite different, as a first core has undergone runaway gas accretion, at about 0.35 Myr. By 0.5 Myr, its mass has already grown to about 350 M. In the end, it will have a mass close to 750 M and be the innermost giant planet. The starting position of this embryo was 4.5 au. The water iceline in this system is for comparison found at about 3.4 au. The formation of this first giant planet does not yet strongly affect the rest of the system, at least at this moment. In the inner system, we in particular see a similar development as in the low-mass disc: the formation of very large resonant convoys and some giant impacts. However, shortly after 0.5 Myr, a second core, located about 0.5 au outside of the first giant, also starts runaway gas accretion. The embryo of this planets started at about 5.3 au, and was for some time in a resonant configuration with the first giant-to-be. It will eventually become the most massive giant planet in the system (about 2100 M) at 1.2 au.

The growth of this second giant planet has important system-wide consequences, as can be seen in the panel at 1 Myr. It not only destabilises several Neptunian planets in the vicinity of the forming giants, but it also sends a protoplanet of about 3 M from about 0.9 au into the inner system (close to 0.1 au). The orbit of this planet is eccentric, and triggers numerous giant impacts among the protoplanets in the inner system (see panel d of Fig. 23). These orbit crossings and impacts are facilitated because the runaway gas accretion by the two forming giants strongly reduces the gas surface density in the inner disc temporarily, reducing eccentric damping (panel c of Fig. 23). In the end, only the intruder from the exterior remains, the mass of which has increased to about 13 M by accreting the local protoplanets. The formation of the second giant also scatters an initially low-mass protoplanet (0.7 M) from about 2 au onto a very eccentric orbit with a semi-major axis of about 15 au. This protoplanet grows then out there (potentially in a monarchical growth mode, Weidenschilling 2005), reaching a mass of about 3 M by 1 Myr.

By about 1.4 Myr, its mass has increased to 8 M, and a phase of rapid inward migration sets in. It then runs from outside into a group of 7 protoplanets at about 2 to 4 au that are captured in MMRs with the giant planet that had formed second (see panel a of Fig. 23. A series of giant impacts occur, and at 1.8 Myr, the protoplanet coming from the outside starts runaway gas accretion. It eventually becomes the third giant planet in the system with a mass of about 630 M at 2.5 au. Interestingly enough, this implies that giant planets in a system need not to be strictly coeval, which could be of importance for example for direct imaging observations. Here, the outermost giant is nearly 1.5 Myr younger, and starts runaway accretion only when the inner two planets have already reached nearly their final mass. Actually, the fact that this third outer planet forms strongly reduces the gas accretion rate of the middle giant, by reducing the gas surface density in the inner system (see panel c of Fig. 23). So, more precisely speaking, the formation of this third giant actually sets the final mass of the giant planet inside of it. A comparable, transient depletion of the inner gas disc is already also seen when the inner two giants form, as mentioned. It should be noted that the degree of depletion of the inner disc because of gas accreting giant planets might be overestimated in our model (Manara et al. 2019; Nayakshin et al. 2019; Bergez-Casalou et al. 2020). Then, this indirect interactions via the disc would be reduced. The lifetime of the gas disc is in this example about 3.4 Myr. This is less than the lifetime of the low mass system studied in the previous section, despite the higher starting mass. The difference is mainly a consequence of the higher external photoevaporation by nearly a factor 5 (it is an independent initial condition, see Paper II). The gas accretion of the giant planets also contributes to the dispersal by them containing in the end about 0.01 M of gas (out of the initial disc mass of 0.066 M).

The temporal evolution of this system shows how the growth of multiple giant planets strongly affects the overall system architecture. This also has important consequences for the giant planets themselves (see panel d of Fig. 23): while they accrete their gas envelopes, they get hit by several lower-mass protoplanets that they destabilise. This increases the core mass of the three giants from about 24, 14 and 10 M at the onset of gas runaway accretion to clearly higher finales values of 64, 26, and 21 M, respectively.Such giants impacts thus strongly influence the final heavy element content (Thorngren & Fortney 2018), and could potentially lead to the existence of diluted cores as found in Jupiter (Liu et al. 2019).

At the end of the simulation at 20 Myr, the system contains four planets more massive than one Earth mass. During the emergence of the system, eight protoplanets have collided with the host star and four were ejected. About 244 M of planetesimals remain out of the starting value of 432 M, corresponding to a difference of 188 M. However, the planets actually existing at the end contain only 123 M, meaning that about 65 M of planetesimals were ‘lost’ because they were either directly ejected or contained in planets that were themselves ejected or fell into the star. This correspond to a solid conversion efficiency of planetesimals into planets of about 28%. Over gigayear timescales, atmospheric escape reduces the mass of the close-in planet at 0.08 au from 13.2 to 11.6 M, but it does retain a remaining H/He envelope. Under the effect of bloating, the planet therefore has a relatively large radius of 5.3 R at 5 Gyr. It is an example of an inner sub-Neptunian planet in a system with outer giant planets (see Paper III).

Finally, it is worth mentioning that systems with three giant planets are statistically a very rare outcome in the population synthesis (Paper II): there are only five such systems among the 1000 synthesised in the nominal population NG76. Systems with one or two giants are in comparison much more common (each about 100 systems). In the system at hand, orbital stability is provided by the giant planets residing in the 3:1 mean motion resonance for both pairs of planets. This allows them to remain stable (Alves et al. 2016) despite their relative proximity to each other, corresponding for both pairs to about 6 to 7 mutual Hill radii, and their significant eccentricities (about 0.08, 0.18, and 0.40 for the inner, middle, and outer planet). We have further tested the stability of this system by extending the orbital integration (including all bodies in the system) from 20 to 100 Myr. At least on this timescale, the system remained stable without secular growth of the eccentricities.

thumbnail Fig. 22

Example of the formation of a planetary system from initially 100 lunar-mass embryos in a high gas mass (initial mass 0.066 M), high metallicity ([Fe/H] = 0.23) disc. The initial mass of planetesimals is 432 M. The plot is analogous to Fig. 19, but the y-axis now extends to much higher masses, and the moments in time that are shown are different. At the end of the simulation at 20 Myr, this system contains one close-in sub-Neptunian planet, three giant planets, and a group of outer very low-mass planets.

thumbnail Fig. 23

Temporal evolution of the system shown in Fig. 22. The four panels have the same meaning as in Fig. 20.

8.2 Overview of the diversity of system architectures

Figures 24 and 25 show the mass-distance and radius-distance of 23 synthetic systems. The solar system is shown in the top-left panel for comparison. All these systems are again taken from the nominal synthetic population NG76 for 1 M star that will be presented in Paper II. However, here we study these as individual systems without taking into account the likelihood of such systems in populations. Hereafter we give an overview of some major correlations that we find. For quantitative results, we refer to the next papers of the series.

It is important to point out a potential complication concerning the formation of the outer two giant planets in the system shown in the panel q of Fig. 24. These planets accreted their cores and started undergoing runaway gas accretion in the inner region of the system. They were subsequently moved to the outer region of the disc where they continue to accrete gas. However, the reason for their final distant locations are not planet-planet scatterings. The presence of a inner massive giant planet (in this case, the one at 0.2 au with nearly 30 M, which corresponds to 3% of the mass of the central star, for an initial disc mass of 9% of the stellar mass) results here in the outer planets obtaining large eccentricities. This, in turn, causes the prescription for the modulation of the torque (Eqs. (86) and (87)) to reverse its sign. Via the additional forces added to the N-body integrator (Sect. 5.2.1), this if found to lead to outward migration in the present case.

Generally speaking, a positive torque means that the angular momentum of a planet has to grow. For an eccentric planet, this can occur via two ways (Cresswell et al. 2007): by eccentricity reduction (circularisation) or outward migration (increase of the semi-major axis). The different approaches how to translate the positive torque found in hydrodynamical simulations into the additional N-body forces have been inconsistent with one another in the literature in the past in this regard (Ida et al. 2020). A re-assessment was recently made in Ida et al. (2020), but is not yet included in the simulations shown here.

The problem we encounter in the special setup here (the presence of an inner very massive giant planet) is likely that the eccentricity state towards which eccentricity damping is acting is becoming ill-defined. The setups used to derive the eccentricity and inclination damping expressions and their translation into additional N-body forces (e.g. Papaloizou & Larwood 2000; Cresswell & Nelson 2008; Bitsch et al. 2013; Ida et al. 2020) assume that the disc orbits on a nearly circular orbit centred on the star. However, in the case here, the planet and outer disc will tend to orbit the barycenter of the star-inner giant pair, which means that the eccentricity can likely not be stabilised near zero. Figuring out the consequences for the orbital evolution of the different involved planets would likely require dedicated hydrodynamical simulations. This shows a limitation of our N-body approach with additional forces instead of direct hydrodynamical simulations. This implies that the model results for distant giant planets with an inner massive planet must be taken with caution.

thumbnail Fig. 24

Mass-distance diagrams of specific systems with 100 embryos initially (panels a–w), which are taken from the nominal population predicted for a 1 M star (NG76). Symbols are as follows: red points show gas-rich planets where MenvMcore > 1. Blue symbols are planets that have accreted some volatile material (ices) outside of the ice line(s). Green symbols are planets that have only accreted refractory solids. Open green and blue circles have 0.1envMcore ≤ 1 while filled green points and blue crosses have MenvMcore ≤ 0.1. For all these bodies, the grey horizontal bars go from ae to a + e. The top leftpanel with black crosses shows the solar system. Bodies lost because of collisions or ejections are shown in light grey. Planets accreted by the central star are show in the very left of each panel, the ejected ones on the very right and planets that collided with another (more massive) planet are shown at their last position on the diagram. The dotted vertical line in each system shows the location of the ice line. The number after each panel name is the metallicity [M/H] of the system expressed in dex, while the value on the top right is the initial mass of the planetesimals disc.

thumbnail Fig. 25

Radius-distance diagrams for a subset of the systems shown in Fig. 24. Here, we focus on systems with multiple low-mass planets; panel letters correspond to the same system in Fig. 24. In contrast to Fig. 24, lost planets are not shown for clarity.

8.2.1 Mass and final number of planets

The number of planets that remain past the formation stage is anti-correlated to the mass of the formed planets. Systems forming giant planets loose more embryos than the ones forming low-mass planets only. We obtained some systems where only one giant planets remains, for instance in panels a and c (including a single one in the latter case), where all the other embryos were removed during the formation stage. When this occurs, at least one of the final planets remains on a wide orbit, as it needs to clear the outer embryos. If this is not the case, then we observe that some embryos with low masses remain in the outer region (e.g. panels e and i).

Systems that still form giant planets, but of lower masses, are able to retain more bodies. We have a few examples that have an architecture in the fashion the solar system, with terrestrial planets inside of giants, such as in panels e, i, m and p. However, those are not comparable to the solar system for several reasons. First, the giant planets are quite more massive than in the solar system; it is not uncommon to find masses of the order of 5–10 M. Likewise, the terrestrial planets are many Earth masses. Further, the location of the giants is much closer in that Jupiter, with distances that are around 1 au. These findings indicate that (1) the gas accretion rate in the disc-limited regime could be too high and (2) the simple Type II migration model we employ in this work leads to too much inward migration.

Finally, systems that form low-mass planets only remain with the largest number of bodies. This is seen for instance in panels l, n, r, t, v and w, where many ice-free bodies (shown in green circles) are present at the end.

8.2.2 Similarity in the low-mass systems

Systems where only terrestrial planets are present have planets with similar properties. It can be seen in panels d, g, h, l, n, r, s, t, v, and w. This is result consistent with observational results about masses and spacing (Millholland et al. 2017). To provide a comparison point with the similarity of planet radii (Weiss et al. 2018), we provide a radius-distance diagram in Fig. 25. For the rocky planets, both masses and radii show the same similarity. The transition from rocky to icy planets affects the radii only slightly. More important is the presence of (remaining) H/He envelopes that were not removed by photoevaporation.

We observe a general slight increase of mass with distance, at least in the inner region. This is most likely linked to the surface density profile of solids. The isolation mass Misor(1.5(2−βs)) (Lissauer 1987), and so since we have βs = 1.5, the value increases with distance. This increase stops at locations usually slightly outside of 1 au, which could be due to our limited integration time, as we discussed in the previous section.

8.2.3 Composition of the close-in planets

We find that close-in terrestrial planets are likely to be rocky, which is in agreement with inferences from observations (Jin & Mordasini 2018). This is especially the case for systems where no planets grow to more than a few Earth masses. We observe in all systems that icy planets are found inside the location of the ice line (the dotted vertical line). Nevertheless, the innermost planets only accrete from the inner region of the disc where the planetesimals are rocky. This indicates that these planets neither migrate from beyond the ice line to their current position, nor get moved to other locations by mean of dynamical instability. It should be noted that in our simulations, planetesimals compositionis set from the initial temperature and pressure profile of the gas disc (Sect. 3.3.3).

Nevertheless, there are systems without giant planet that consist of only ice-bearing bodies; these are shown in panels d, g and u. These systems form planets that are more massive than the previous ones, with most of them having at least one planet above 10 M. The Type I migration timescale decreases with increasing mass; therefore these more massive planets can migrate from outside of the water ice line to their current position, increasing the compositional diversity of the systems (Raymond et al. 2018).

Systems with giant planets exhibit different behaviours. Some have only ice-bearing planets (panels b, f and q) while others have also terrestrial planets. In the latter case, the giant planets do not necessarily separate rocky bodies from icy ones. Panels e and o show systems where rocky and icy planets are separated by giants, while in panels i, m and p icy planets are present both inside and outside of the gas giants. This points at a high diversity of the composition of planets in systems containing both giant and low-mass planets. Correlations between the occurrences of giant planets and others in planetary systems will be investigated in more details in Paper II. Schlecker et al. (2021a, hereafter Paper III) will look thoroughly at correlations between close-in Super-Earth planets and long-period giants.

9 Summary and conclusions

In this work, we presented the Generation III version of the Bern global model of planetary formation and evolution. In this generation, the following two main aspects were improved. First is the ability to simulate planets with a mass range from Mars to deuterium-burning planets. Older generations of the Bern model could not address terrestrial planets, as they we lacking the giant-impact stage. To reach this goal, we improved the N-body integrator so that per disc, hundreds of concurrently forming embryos can now be included. This is crucial for the formation of low-mass planets in general and the Solar System. We also added several new physical processes to take into account the consequences of stellar proximity, allowing us to simulate with the new model planets that cover the widest range of orbital separations, from star-grazing to distant and even rogue planets. Second, the ability to predict self-consistently for multi-planet systems as many directly observable quantities as possible: not only masses and orbital elements as in the past, but also other key observables like luminosities, magnitudes, transit radii, or evaporation rates. To achieve this, we coupled our planet formation model (to 20 Myr) to our planet evolution model (20 Myr–10 Gyr). Thanks to this, we can now self-consistently and statistically compare the same population to all important observational techniques, as will be done in the series of NGPPS papers. This is crucial, as different methods probe distinct planetary sub-populations. This combined comparison puts extremely compelling and powerful constraints on any theoretical model.

The formation and evolution model follows the envelope structure of the giant planets during they entire lifetime. This allows for example to study the luminosities at any time (Mordasini et al. 2017), and enables the comparison with directly-imaged exoplanets (e.g. Vigan et al. 2017).

The model now includes a multitude of physical processes (see Fig. 2). The following are included during both the formation and evolution phase:

  • A solution of 1D radially symmetric internal structure equations (Bodenheimer & Pollack 1986) is used to calculate the internal structure of the H/He envelope and thus the gas accretion rate (during the attached phase), radius and luminosity, which includes Deuterium burning (Mollière & Mordasini 2012) and bloating of close-in planets.

  • The solution of the 1D internal core structure is used to obtain the radius of the solid core with a modified polytropic EOS (Seager et al. 2007).

  • An atmospheric model yields the outer boundary conditions during the attached, detached, and evolutionary phase. For the detached phase, we assume hot gas accretion. For the evolutionary phase, we use a simple grey atmosphere.

  • The host star properties are retrieved from tabulated stellar evolution tracks (Baraffe et al. 2015).

During formation, the following processes are included:

  • The radial structure of the protoplanetary gas disc is computed with a 1D radial (axis-symmetric) constant α-disc model. The effects of internal and external photoevaporation are included.

  • The vertical structure of the disc is modelled by building on radiative equilibrium (Nakamoto & Nakagawa 1994), including viscous heating and stellar irradiation (Fouchet et al. 2012). Irradiation now also includes the direct irradiation in the disc midlplane important when the disc becomes optically thin.

  • Planetesimals are presented by a 1D radial (axis-symmetric) disc, with a surface density and a dynamical state (eccentricity, inclination). The temporal evolution of e and i are explicitly followed, including the dynamic excitation by protoplanets and planetesimals, and damping from gas drag (Fortier et al. 2013). The composition of the planetesimal and the position of ice lines is found from an equilibrium condensation model (Thiabaud et al. 2015).

  • The equation for the planetesimal accretion rate of the protoplanet is computed assuming the oligarchic regime (Chambers 2006). The enhancement of the planetesimal capture radius because of the planetary H/He envelope is included (Inaba & Ikoma 2003).

  • A prescription based on Bondi- and Hill-type gas accretion in the 2D and 3D cases limits the planetary gas accretion rate in the disc-limited regime.

  • Gas-driven Type I and Type II orbital migration are computed including the effects of non-isothermality and of the planet’s eccentricity and inclination (Paardekooper et al. 2011; Coleman & Nelson 2014; Dittkrist et al. 2014).

  • Full N-body interaction between all the embryos forming concurrently in one disc are tracked using the mercury integrator (Chambers 1999). Orbital migration and the damping of eccentricity and inclination are input in the integrator via additional forces. In case of a collision, the impact energy is added as an additionally luminosity term (Broeg & Benz 2012) to the internal structure model. This can lead to the loss of the H/He envelope.

During the evolutionary phase we include:

  • XUV-driven atmospheric photoevaporation in the energy and radiation-recombination-limited approximation (Jin et al. 2014),

  • for close-in planets, the addition of a bloating luminosity modelled with the empirical relation of Thorngren & Fortney (2018), and

  • tidal spiral-in because of stellar tides (Benítez-Llambay et al. 2011), along with Roche-lobe overflow.

We show in Sect. 6 where we study the formation of terrestrial planets that provided there are initially enough embryo in each disc, mutual gravitational interactions will stir their eccentricities. Due to the radial excursions, embryos will have access to more material until all the planetesimals are accreted. Afterwards, a phase of giant impacts sets in. Thus, despite the use a fluid-like description for the planetesimals, the model is able to reproduce the giant impact phase of terrestrial planet formation. Due to the limitation of the integration time (20 Myr), this is only completely modelled within a distance of roughly 1 au. Giant planets, in contrast, are not affected by the integration time limitation as they must anyway form before the dispersal of the gas disc. The model is then able to track the formation of all planets in the inner part of planetary systems.

After the description of the model, we study how the many different sub-models included in the Bern Generation III Model interact in the full end-to-end model by simulating the formation of two planetary systems. To understand the results, it is helpful to compare the timescales of growth and migration, revealing which process is dominant. It is also helpful to study the planetesimal surface density, revealing the solid accretion mode (planetesimal accretion versus growth by giant impacts). Other key processes occurring during the emergence of the planetary systems include the capture of many protoplanets into large resonant convoys, and the consequences of dynamical instabilities caused by the gravitational interactions between the protoplanets. This includes the destabilisation of other protoplanets at the moment a giant planet (especially a second one in the system) starts runaway gas accretion as well as series of giant impacts at the moment the gas disc dissipates.

We also give a short overview of the diversity of planetary systems that were obtained using the model. We find that systems containing giant planets can have a great diversity of configurations, while for systems forming only low-mass (Earth-like) planets exhibit arranged planets with similar masses.

This work is the first of a series. Here we present the outline of the series:

  • Paper II will introduce the methods to calculate population syntheses. Several populations for Solar-mass stars with different numbers of initial embryos per system are computed. The effects of this parameter at the population level will be investigated.

  • Paper III will look for correlations between of the occurrence of inner low-mass and outer giant planets.

  • Paper IV (Burn et al. 2021) will extend the population synthesis to lower-mass stars (down to late M-dwarfs) and analyse the effects of the stellar mass.

  • Paper V(Schlecker et al. 2021b) will study the mapping of disc initial conditions to planet properties with machine learning.

  • Paper VI (Mishra et al. 2021) will look for the diversity between planets in each system compared to diversity of the overall population (Weiss et al. 2018).

  • There are then three papers on the quantitative comparison with various observational techniques: radial velocity with HARPS and CARMENES, and transits with Kepler.

The discrepancies uncovered in these comprehensive and multi-aspect comparisons with observations will be helpful to improve the understanding of planet formation and evolution.

Acknowledgements

The authors thank David Swoboda, Natacha Brügger, Martin Schlecker, and the members of the Theoretical Astrophysics and Planetary Science (TAPS) group at the University of Bern for fruitful discussions. We also thank the anonymous reviewer, whose remarks and suggestions helped improve the quality of this manuscript. A.E. and E.A. acknowledge the support from The University of Arizona. A.E. and C.M. acknowledge the support from the Swiss National Science Foundation under grant BSSGI0_155816 ‘PlanetsInTime’. Parts of this work have been carried out within the frame of the National Center for Competence in Research PlanetS supported by the SNSF. The plots shown in this work were generated using matplotlib (Hunter 2007).

Appendix A Possible impact of migration on the accretion of planetesimals

A.1 General considerations

Ward (1986) and Ward (1989) were the first to suggested that the orbital migration of a protoplanet could strongly accelerate its planetesimal accretion if the planet is able to catch most planetesimals through which it is sweeping during its migration. In the terminology of Ward, the protoplanet would then act as a predator.

On the other hand, Tanaka & Ida (1999) found that in the case of a slow migration timescale τmig, the protoplanet rather carves a gap in the planetesimals disc around its orbit and shepherds the planetesimals, which stalls the accretion of the protoplanet. They derived a criterion in terms of a critical migration timescale τmig,crit only below which the protoplanets can act as predator, τmigτmig,crit. Otherwise, the protoplanets acts as shepherd.

Alibert et al. (2005a) studied this effect with the Generation I Bern Model and found that the migration rates were generally high enough for the protoplanets to be predators, and thus ignored the shepherding effect. Since the Generation 3 Model differs in numerous aspects (like the oligarchic growth mode, the disc model, the planetesimal size, the migration model, or the multiplicity of forming planets) from the model of 2005, we here re-assess this question. We do this based on existing simulations and published criteria. It is clear that for a more definitive assessment, direct N-body simulations of thousands of concurrently growing and migrating oligarchs accreting planetesimals and gas in the setup that we consider here would be needed. In these simulations, the planetesimals would have to be included as individual fully interacting bodies in the N-body. This differs from our current approach of representing planetesimals statistically as a surface density with a dynamic state (eccentricity and inclination, see Sect. 3.3).

A.2 Applicability of Tanaka & Ida (1999)

Coming to the general applicability of the work of Tanaka & Ida (1999), one should note that they studied a highly special setup. In what follows, we discuss several important effects they neglected. When considering them, it becomes much less clear whether the shepherding effect can at all represent a major impediment to planetesimal accretion. We mention these points also in view of the specific results discussed in Sect. A.3 where we try to identify planets for which shepherding could potentaially be relevant.

First, the predator-shepherd mechanism of Tanaka & Ida (1999) was found for a single core growing alone in a disc of planetesimals. This seems a highly unlikely setup given that many protoplanets (oligarchs) emerge concurrently from the pre-dating runaway planetesimal accretion phase and that subsequently grow further in lockstep (e.g. Kokubo & Ida 1998). As Tanaka & Ida (1999) state actually themselves, in the case of such multiple protoplanets, ‘the protoplanets push planetesimals into the feeding zones of others and they can grow.’

Indeed, the more recent and realistic work of Daisaka et al. (2006) who run N-body simulation including type-I migration and tidal damping of eccentricity and inclination, starting from 7000-14’000 equal-mass self-gravitating planetesimals whose size is roughly 1000 km showed a different picture than the single protoplanet simulations of Tanaka & Ida (1999): namely that in the multiple protoplanet situation, the trapping of planetesimals by cores is only tentative and does not significantly reduce their accretion rates. This was already noted by Ida & Lin (2008), motivating them to not modify the core accretion rate for migrating protoplanets.

Second, according to the Tanaka & Ida (1999) timescale criterion one finds (at least formally) that whenever the migration timescale is longer than the critical value, no planetesimal accretion occurs. This includes in particular the case that a planet does not migrate at all. However, as shown by many works (Kokubo & Ida 1998; Levison et al. 2010; Alibert et al. 2013), multiple protoplanets forming concurrently can grow in-situ by planetesimals that are already in their feeding zones. When a protoplanet grows, its feeding zone which is proportional to the Hill sphere and thus the mass of the protoplanet, also expands, bringing new planetesimals in reach of the protoplanet in a process that is unaffected by shepherding. The protoplanet’s mass growth can occur via the accretion of solids, but also via the accretion of gas. In fact, the interplay of gas accretion leading to an extension of the solid feeding zone which in turn increases the core mass, and then again the gas accretion rate is the underlying process for Phase II in Pollack et al. (1996). Because of the absence of a gap in planetesimal disc around the planet if many protoplanets grow concurrently thanks to scattering (Levison et al. 2010; Alibert et al. 2013), each protoplanet has a local reservoir of planetesimals from which it grows. Thus, when judging whether shepherding could affect growth, we do not need to consider planets that have not (yet) migrated outside of the planetesimal feeding zone around their starting location, since they still accrete their local reservoir.

Third, the temporal sequence of how solid accretion proceeds in forming multi-planet systems itself reduces the impact of shepherding. Here it is important to consider that solid growth of protoplanets occurs via two channels: the accretion of planetesimals, but also via protoplanet-protoplanet collisions (giant impacts). As can be seen in the growth tracks in Figures 19 and 22, at the lower masses, when the planetesimal accretion timescale is less than the migration timescale, the protoplanets grow nearly in-situ by accreting planetesimals from their local feeding zone. As explained before, in this case shepherding is not important. In the subsequent phase, when the planets have grown to a mass where they start migrating, they also start accreting other protoplanets via giant impacts, which is also unaffected by shepherding. They also often migrate into parts of the disc which were previously depleted in planetesimals by another protoplanet that has formed further in. Here, shepherding would not have an effect. This insight also further justifies why using the in-situ prescription for the planetesimal accretion rate is appropriate. This also implies that when identifying planets that could be affected by shepherding, we do not need to consider planets that accrete in any case only very slowly planetesimals even when neglecting shepherding.

Finally, we also do not need to consider protoplanets where the gas disc has already dissipated. In this case, no gas disc-driven orbital migration occurs and thus no shepherding.

A.3 Assessment of the population-wide importance of shepherding

Except for the first point (the effect of protoplanet multiplicity) which questions the existence of shepherding in a fundamental way, the considerations from the previous section can - in an approximate way - be cast into a set of conditions where shepherding could be important. This allows us to identify these protoplanets which might at least in principle be affected by shepherding.

thumbnail Fig. A.1

Mass-distance diagram of the nominal synthetic population NG76 of solar-like starts which stars with initially 100 moon-mass embryos per disc (see Paper II). The epochs of 0.1, 0.5, and 1 Myr are shown. Coloured points show protoplanets that can no longer accrete planetesimals of the initial local reservoir, which have a planetesimal accretion timescale of less than 3 Myr, and which are still embedded in the parent gaseous disc. When , the planetesimal accretion of these planets could in principle be affected by shepherding if they would be the only protoplanets growing in the disc.

Figure A.1 shows the mass-distance diagram of the nominal synthetic population NG76 from Paper II. Three moments in time are shown where planetesimal accretion is in general important. We colour code the absolute value of the ratio of the normalised migration timescale of a planet to the normalised critical migration timescale which both are calculated as in Tanaka & Ida (1999). When this ratio is larger than approximately unity, shepherding would occur for a single protoplanet migrating alone through a disc of planetesimals. Only protoplanets which could in principle be affected by shepherding by fulfilling the following criteria are colour coded: First, the distance a planet has migrated away from its starting location is larger than five times the size of its Hill sphere. This means that it can no longer accrete from its initial local reservoir of planetesimals. Second, the planetesimal accretion timescale is less than three million years (the typical disc lifetime), meaning that planetesimal accretion (as opposed to growth via giant impacts) is still relevant. Third, the gas disc has not yet dissipated. Other protoplanets should in any case not be significantly affected by shepherding and are shown in grey.

The plot first shows that the large majority of protoplanets are grey, meaning that shepherding should not be important for them in any case. Then, more specifically, at 0.1 Myr, there is a group of Mars- to Earth-mass protoplanets inside of the ice line where . At 0.5 Myr, there is a radial interval from about 1 to 4 au where is longer than , however in most cases by less than one order of magnitude. These are regions where usually groups of tens of protoplanets form together (see Sect. 8.1), so that it is not clear if shepherding would occur at all. Planets where the ratio is clearly larger, and thus where the effect could in principle be particularly strong, are rare. At 1 Myr, a similar pattern is seen, but the potentially affected region is reduced.

It is clear that this simple a posteriori analysis cannot be seen as a final result - for this, simulations where planetesimal are included directly in the N-body would be necessary. Nevertheless, together with the finding of Daisaka et al. (2006) that shepherding is by principle not important when several protoplanets form concurrently, they indicate that shepherding can only affect a relatively limited part of all growing protoplanets.

References

  1. Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 [Google Scholar]
  2. Adibekyan, V. 2019, Geosciences, 9, 105 [Google Scholar]
  3. Alessi, M., & Pudritz, R. E. 2018, MNRAS, 478, 2599 [Google Scholar]
  4. Alessi, M., Pudritz, R. E., & Cridland, A. J. 2020a, MNRAS, 493, 1013 [Google Scholar]
  5. Alessi, M., Inglis, J., & Pudritz, R. E. 2020b, MNRAS, 497, 4814 [Google Scholar]
  6. Alexander, R. D., & Armitage, P. J. 2009, ApJ, 704, 989 [Google Scholar]
  7. Alexander, R. D., & Pascucci, I. 2012, MNRAS, 422, L82 [Google Scholar]
  8. Alibert, Y., & Benz, W. 2017, A&A, 598, L5 [Google Scholar]
  9. Alibert, Y., & Venturini, J. 2019, A&A, 626, A21 [Google Scholar]
  10. Alibert, Y., Mordasini, C., & Benz, W. 2004, A&A, 417, L25 [Google Scholar]
  11. Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005a, A&A, 434, 343 [Google Scholar]
  12. Alibert, Y., Mousis, O., Mordasini, C., & Benz, W. 2005b, ApJ, 626, L57 [Google Scholar]
  13. Alibert, Y., Mordasini, C., & Benz, W. 2011, A&A, 526, A63 [Google Scholar]
  14. Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [Google Scholar]
  15. Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nat. Astron., 2, 873 [Google Scholar]
  16. Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440 [Google Scholar]
  17. Allègre, C. J., Manhès, G., & Göpel, C. 2008, Earth Planet. Sci. Lett., 267, 386 [Google Scholar]
  18. Alves, A. J., Michtchenko, T. A., & Tadeu dos Santos, M. 2016, Celest. Mech. Dyn. Astron., 124, 311 [Google Scholar]
  19. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241 [Google Scholar]
  20. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [Google Scholar]
  21. Arimatsu, K., Tsumura, K., Usui, F., et al. 2019, Nat. Astron., 3, 301 [Google Scholar]
  22. Bailey, E., & Batygin, K. 2018, ApJ, 866, L2 [Google Scholar]
  23. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403 [Google Scholar]
  24. Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315 [Google Scholar]
  25. Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 [Google Scholar]
  26. Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 667 [Google Scholar]
  27. Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77 [Google Scholar]
  28. Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114 [Google Scholar]
  29. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [Google Scholar]
  30. Benítez-Llambay, P., Masset, F., & Beaugé, C. 2011, A&A, 528, A2 [Google Scholar]
  31. Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 691 [Google Scholar]
  32. Berardo, D., & Cumming, A. 2017, ApJ, 846, L17 [Google Scholar]
  33. Berardo, D., Cumming, A., & Marleau, G.-D. 2017, ApJ, 834, 149 [Google Scholar]
  34. Bergez-Casalou, C., Bitsch, B., Pierens, A., Crida, A., & Raymond, S. N. 2020, A&A, 643, A133 [Google Scholar]
  35. Birnstiel, T., & Andrews, S. M. 2014, ApJ, 780, 153 [Google Scholar]
  36. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [Google Scholar]
  37. Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [Google Scholar]
  38. Bitsch, B., & Kley, W. 2011, A&A, 530, A41 [Google Scholar]
  39. Bitsch, B., Crida, A., Libert, A. S., & Lega, E. 2013, A&A, 555, A124 [Google Scholar]
  40. Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015a, A&A, 575, A28 [Google Scholar]
  41. Bitsch, B., Lambrechts, M., & Johansen, A. 2015b, A&A, 582, A112 [Google Scholar]
  42. Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, A&A, 623, A88 [Google Scholar]
  43. Bitsch, B., Trifonov, T., & Izidoro, A. 2020, A&A, 643, A66 [EDP Sciences] [Google Scholar]
  44. Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 [Google Scholar]
  45. Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2 [Google Scholar]
  46. Bodenheimer, P., D’Angelo, G., Lissauer, J. J., Fortney, J. J., & Saumon, D. 2013, ApJ, 770, 120 [Google Scholar]
  47. Boley, A. C., Granados Contreras, A. P., & Gladman, B. 2016, ApJ, 817, L17 [Google Scholar]
  48. Bolmont, E., & Mathis, S. 2016, Celest. Mech. Dyn. Astron., 126, 275 [Google Scholar]
  49. Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998 [Google Scholar]
  50. Boss, A. P. 1997, Science, 276, 1836 [Google Scholar]
  51. Boss, A. P. 2003, ApJ, 599, 577 [Google Scholar]
  52. Broeg, C. H., & Benz, W. 2012, A&A, 538, A90 [Google Scholar]
  53. Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89 [Google Scholar]
  54. Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344 [Google Scholar]
  55. Buchhave, L. A., Bitsch, B., Johansen, A., et al. 2018, ApJ, 856, 37 [Google Scholar]
  56. Burn, R., Marboeuf, U., Alibert, Y., & Benz, W. 2019, A&A, 629, A64 [Google Scholar]
  57. Burn, R., Schlecker, M., Mordasini, C., et al. 2021, A&A, 656, A72 (Paper IV) [Google Scholar]
  58. Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Rev. Mod. Phys., 73, 719 [Google Scholar]
  59. Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47 [Google Scholar]
  60. Chabrier, G., & Baraffe, I. 2000, ARA&A, 38, 337 [Google Scholar]
  61. Chambers, J. E. 1999, MNRAS, 304, 793 [Google Scholar]
  62. Chambers, J. E. 2001, Icarus, 152, 205 [Google Scholar]
  63. Chambers, J. 2006, Icarus, 180, 496 [Google Scholar]
  64. Chambers, J. E. 2013, Icarus, 224, 43 [Google Scholar]
  65. Chambers, J. 2018, ApJ, 865, 30 [Google Scholar]
  66. Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580 [Google Scholar]
  67. Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, L9 [Google Scholar]
  68. Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180 [Google Scholar]
  69. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
  70. Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 [Google Scholar]
  71. Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 [Google Scholar]
  72. Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479 [Google Scholar]
  73. Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677 [Google Scholar]
  74. Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329 [Google Scholar]
  75. Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [Google Scholar]
  76. Cridland, A. J., Pudritz, R. E., & Alessi, M. 2016, MNRAS, 461, 3274 [Google Scholar]
  77. Cridland, A. J., Pudritz, R. E., Birnstiel, T., Cleeves, L. I., & Bergin, E. A. 2017, MNRAS, 469, 3910 [Google Scholar]
  78. Daisaka, J. K., Tanaka, H., & Ida, S. 2006, Icarus, 185, 492 [Google Scholar]
  79. D’Angelo, G., & Lubow, S. H. 2008, ApJ, 685, 560 [Google Scholar]
  80. D’Angelo, G., Weidenschilling, S. J., Lissauer, J. J., & Bodenheimer, P. 2021, Icarus, 355, 114087 [Google Scholar]
  81. Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175 [Google Scholar]
  82. Dawson, R. I., & Murray-Clay, R. A. 2013, ApJ, 767, L24 [Google Scholar]
  83. Dawson, R. I., Chiang, E., & Lee, E. J. 2015, MNRAS, 453, 1471 [Google Scholar]
  84. Demory, B.-O., & Seager, S. 2011, ApJS, 197, 12 [Google Scholar]
  85. Dewitt, H. E., Graboske, H. C., & Cooper, M. S. 1973, ApJ, 181, 439 [Google Scholar]
  86. Dittkrist, K. M., Mordasini, C., Klahr, H., Alibert, Y., & Henning, T. 2014, A&A, 567, A121 [Google Scholar]
  87. Emsenhuber, A., Cambioni, S., Asphaug, E., et al. 2020, ApJ, 891, 6 [Google Scholar]
  88. Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A70 (Paper II) [Google Scholar]
  89. Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [Google Scholar]
  90. Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92 [Google Scholar]
  91. Fatuzzo, M., & Adams, F. C. 2008, ApJ, 675, 1361 [Google Scholar]
  92. Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72 [Google Scholar]
  93. Fendyke, S. M., & Nelson, R. P. 2014, MNRAS, 437, 96 [Google Scholar]
  94. Ferraz-Mello, S., Rodríguez, A., & Hussmann, H. 2008, Celest. Mech. Dyn. Astron., 101, 171 [Google Scholar]
  95. Ford, E. B., & Rasio, F. A. 2008, ApJ, 686, 621 [Google Scholar]
  96. Fortier, A., Alibert, Y., Carron, F., Benz, W., & Dittkrist, K. M. 2013, A&A, 549, A44 [Google Scholar]
  97. Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80 [Google Scholar]
  98. Fouchet, L., Alibert, Y., Mordasini, C., & Benz, W. 2012, A&A, 540, A107 [Google Scholar]
  99. Fowler, W. A., Caughlan, G. R., & Zimmerman, B. A. 1967, ARA&A, 5, 525 [Google Scholar]
  100. Fraser, W. C. 2009, ApJ, 706, 119 [Google Scholar]
  101. Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25 [Google Scholar]
  102. Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [Google Scholar]
  103. Golabek, G. J., & Jutzi, M. 2021, Icarus, 363, 114437 [Google Scholar]
  104. Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [Google Scholar]
  105. Goldreich, P., Lithwick, Y., & Sari, R. 2004, ApJ, 614, 497 [Google Scholar]
  106. Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29 [Google Scholar]
  107. Graboske, H. C., Dewitt, H. E., Grossman, A. S., & Cooper, M. S. 1973, ApJ, 181, 457 [Google Scholar]
  108. Grasset, O., Schneider, J., & Sotin, C. 2009, ApJ, 693, 722 [Google Scholar]
  109. Greenzweig, Y., & Lissauer, J. J. 1992, Icarus, 100, 440 [Google Scholar]
  110. Guilera, O. M., Brunini, A., & Benvenuto, O. G. 2010, A&A, 521, A50 [Google Scholar]
  111. Guillot, T. 2005, AREPS, 33, 493 [Google Scholar]
  112. Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [Google Scholar]
  113. Hatzes, A. P., & Rauer, H. 2015, ApJ, 810, L25 [Google Scholar]
  114. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [Google Scholar]
  115. Hayashi, C., Nakazawa, K., & Adachi, I. 1977, PASJ, 29, 163 [NASA ADS] [Google Scholar]
  116. Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654 [Google Scholar]
  117. Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [Google Scholar]
  118. Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [Google Scholar]
  119. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  120. Ida, S., & Lin, D. N. C. 2004a, ApJ, 604, 388 [Google Scholar]
  121. Ida, S., & Lin, D. N. C. 2004b, ApJ, 616, 567 [Google Scholar]
  122. Ida, S., & Lin, D. N. C. 2008, ApJ, 673, 487 [Google Scholar]
  123. Ida, S., & Lin, D. N. C. 2010, ApJ, 719, 810 [Google Scholar]
  124. Ida, S., & Makino, J. 1993, Icarus, 106, 210 [Google Scholar]
  125. Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292 [Google Scholar]
  126. Ida, S., Muto, T., Matsumura, S., & Brasser, R. 2020, MNRAS, 494, 5666 [Google Scholar]
  127. Inaba, S., & Ikoma, M. 2003, A&A, 410, 711 [Google Scholar]
  128. Inaba, S., Tanaka, H., Nakazawa, K., Wetherill, G. W., & Kokubo, E. 2001, Icarus, 149, 235 [Google Scholar]
  129. Inamdar, N. K., & Schlichting, H. E. 2016, ApJ, 817, L13 [Google Scholar]
  130. Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [Google Scholar]
  131. Izidoro, A., Bitsch, B., Raymond, S. N., et al. 2021, A&A, 650, A152 [EDP Sciences] [Google Scholar]
  132. Jackson, B., Barnes, R., & Greenberg, R. 2009, ApJ, 698, 1357 [Google Scholar]
  133. Jackson, A. P., Davis, T. A., & Wheatley, P. J. 2012, MNRAS, 422, 2024 [Google Scholar]
  134. Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [Google Scholar]
  135. Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [Google Scholar]
  136. Johansen, A., & Bitsch, B. 2019, A&A, 631, A70 [Google Scholar]
  137. Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475 [NASA ADS] [Google Scholar]
  138. Johansen, A., Ida, S., & Brasser, R. 2019, A&A, 622, A202 [Google Scholar]
  139. Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603 [Google Scholar]
  140. Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [Google Scholar]
  141. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [Google Scholar]
  142. Kippenhahn, R., & Weiagert, A. 1994, Stellar Structure and Evolution (Berlin: Springer), 192 [Google Scholar]
  143. Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54 [Google Scholar]
  144. Kleine, T., Münker, C., Mezger, K., & Palme, H. 2002, Nature, 418, 952 [Google Scholar]
  145. Kleine, T., Touboul, M., Bourdon, B., et al. 2009, Geochim. Cosmochim. Acta, 73, 5150 [Google Scholar]
  146. Kley, W. 2019, Planet Formation and Disk-Planet Interactions (Berlin, Heidelberg: Springer Berlin Heidelberg), 151 [Google Scholar]
  147. Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [Google Scholar]
  148. Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 [Google Scholar]
  149. Knutson, H. A., Fulton, B. J., Montet, B. T., et al. 2014, ApJ, 785, 126 [Google Scholar]
  150. Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [Google Scholar]
  151. Kokubo, E., Kominami, J., & Ida, S. 2006, ApJ, 642, 1131 [Google Scholar]
  152. Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375 [Google Scholar]
  153. Kretke, K. A., & Lin, D. N. C. 2012, ApJ, 755, 74 [Google Scholar]
  154. Krivov, A. V., & Wyatt, M. C. 2021, MNRAS, 500, 718 [Google Scholar]
  155. Krumholz, M. R. 2006, ApJ, 641, L45 [Google Scholar]
  156. Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [Google Scholar]
  157. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
  158. Lambrechts, M., & Lega, E. 2017, A&A, 606, A146 [Google Scholar]
  159. Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83 [Google Scholar]
  160. Lammer, H., Odert, P., Leitzinger, M., et al. 2009, A&A, 506, 399 [Google Scholar]
  161. Laughlin, G., Crismani, M., & Adams, F. C. 2011, ApJ, 729, L7 [Google Scholar]
  162. Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41 [Google Scholar]
  163. Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90 [Google Scholar]
  164. Lenz, C. T., Klahr, H., & Birnstiel, T. 2019, ApJ, 874, 36 [Google Scholar]
  165. Levison, H. F., Thommes, E., & Duncan, M. J. 2010, AJ, 139, 1297 [Google Scholar]
  166. Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606 [Google Scholar]
  167. Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85 [Google Scholar]
  168. Lissauer, J. J. 1987, Icarus, 69, 249 [Google Scholar]
  169. Lissauer, J. J., Hubickyj, O., D’Angelo, G., & Bodenheimer, P. 2009, Icarus, 199, 338 [Google Scholar]
  170. Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8 [Google Scholar]
  171. Liu, B., Ormel, C. W., & Lin, D. N. C. 2017, A&A, 601, A15 [Google Scholar]
  172. Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355 [Google Scholar]
  173. Lodders, K. 2003, ApJ, 591, 1220 [Google Scholar]
  174. Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [Google Scholar]
  175. Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [Google Scholar]
  176. Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 [Google Scholar]
  177. Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526 [Google Scholar]
  178. Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 [Google Scholar]
  179. Lüst, R. 1952, Zeitschrift Naturforschung Teil A, 7, 87 [Google Scholar]
  180. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
  181. Lyra, W., Paardekooper, S.-J., & Mac Low, M.-M. 2010, ApJ, 715, L68 [Google Scholar]
  182. Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64 [Google Scholar]
  183. Mallonn, M., Köhler, J., Alexoudi, X., et al. 2019, A&A, 624, A62 [Google Scholar]
  184. Manara, C. F., Mordasini, C., Testi, L., et al. 2019, A&A, 631, L2 [Google Scholar]
  185. Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014a, A&A, 570, A36 [Google Scholar]
  186. Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014b, A&A, 570, A35 [Google Scholar]
  187. Marleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C. 2017, ApJ, 836, 221 [Google Scholar]
  188. Marleau, G.-D., Coleman, G. A. L., Leleu, A., & Mordasini, C. 2019a, A&A, 624, A20 [Google Scholar]
  189. Marleau, G.-D., Mordasini, C., & Kuiper, R. 2019b, ApJ, 881, 144 [Google Scholar]
  190. Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348 [Google Scholar]
  191. Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478 [Google Scholar]
  192. Matsumoto, Y., & Ogihara, M. 2020, ApJ, 893, 43 [Google Scholar]
  193. Matsumura, S., Brasser, R., & Ida, S. 2021, A&A, 650, A116 [Google Scholar]
  194. Matsuyama, I., Johnstone, D., & Hartmann, L. 2003, ApJ, 582, 893 [Google Scholar]
  195. Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
  196. Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv e-prints [arXiv:1109.2497] [Google Scholar]
  197. Miguel, Y., Guilera, O. M., & Brunini, A. 2011, MNRAS, 417, 314 [Google Scholar]
  198. Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33 [Google Scholar]
  199. Mills, S. M., Howard, A. W., Petigura, E. A., et al. 2019, AJ, 157, 198 [Google Scholar]
  200. Mishra, L., Alibert, Y., Leleu, A., et al. 2021, A&A, 656, A74 (Paper VI) [Google Scholar]
  201. Mizuno, H. 1980, Prog. Theor. Phys., 64, 544 [Google Scholar]
  202. Moldenhauer, T. W., Kuiper, R., Kley, W., & Ormel, C. W. 2021, A&A, 646, L11 [Google Scholar]
  203. Mollière, P., & Mordasini, C. 2012, A&A, 547, A105 [Google Scholar]
  204. Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558 [Google Scholar]
  205. Mordasini, C. 2014, A&A, 572, A118 [Google Scholar]
  206. Mordasini, C. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte (Berlin: Springer Living Reference Work), 143 [Google Scholar]
  207. Mordasini, C., Alibert, Y., & Benz, W. 2009a, A&A, 501, 1139 [Google Scholar]
  208. Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009b, A&A, 501, 1161 [Google Scholar]
  209. Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012a, A&A, 541, A97 [Google Scholar]
  210. Mordasini, C., Alibert, Y., Georgy, C., et al. 2012b, A&A, 547, A112 [Google Scholar]
  211. Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012c, A&A, 547, A111 [Google Scholar]
  212. Mordasini, C., Klahr, H., Alibert, Y., Miller, N., & Henning, T. 2014, A&A, 566, A141 [Google Scholar]
  213. Mordasini, C., Mollière, P., Dittkrist, K. M., Jin, S., & Alibert, Y. 2015, Int. J. Astrobiol., 14, 201 [Google Scholar]
  214. Mordasini, C., Marleau, G. D., & Mollière, P. 2017, A&A, 608, A72 [Google Scholar]
  215. Movshovitz, N., & Podolak, M. 2008, Icarus, 194, 368 [Google Scholar]
  216. Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [Google Scholar]
  217. Mulders, G. D., O’Brien, D. P., Ciesla, F. J., Apai, D., & Pascucci, I. 2020, ApJ, 897, 72 [Google Scholar]
  218. Müller, S., Ben-Yami, M., & Helled, R. 2020, ApJ, 903, 147 [Google Scholar]
  219. Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [Google Scholar]
  220. Nakamoto, T.,& Nakagawa, Y. 1994, ApJ, 421, 640 [Google Scholar]
  221. Nayakshin, S. 2010, MNRAS, 408, L36 [Google Scholar]
  222. Nayakshin, S., Dipierro, G., & Szulágyi, J. 2019, MNRAS, 488, L12 [Google Scholar]
  223. Ndugu, N., Bitsch, B., & Jurua, E. 2018, MNRAS, 474, 886 [Google Scholar]
  224. O’Brien, D. P., Morbidelli, A., & Levison, H. F. 2006, Icarus, 184, 39 [Google Scholar]
  225. Ogihara, M., Kokubo, E., Suzuki, T. K., & Morbidelli, A. 2018, A&A, 615, A63 [Google Scholar]
  226. Ohtsuki, K., Stewart, G. R., & Ida, S. 2002, Icarus, 155, 436 [Google Scholar]
  227. Okuzumi, S., & Ormel, C. W. 2013, ApJ, 771, 43 [Google Scholar]
  228. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [Google Scholar]
  229. Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 [Google Scholar]
  230. Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [Google Scholar]
  231. Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43 [Google Scholar]
  232. Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [Google Scholar]
  233. Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
  234. Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293 [Google Scholar]
  235. Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823 [Google Scholar]
  236. Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [Google Scholar]
  237. Perri, F., & Cameron, A. G. W. 1974, Icarus, 22, 416 [Google Scholar]
  238. Pfyffer, S., Alibert, Y., Benz, W., & Swoboda, D. 2015, A&A, 579, A37 [Google Scholar]
  239. Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 [Google Scholar]
  240. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [Google Scholar]
  241. Pringle, J. E. 1981, ARA&A, 19, 137 [Google Scholar]
  242. Rafikov, R. R. 2004, AJ, 128, 1348 [Google Scholar]
  243. Rafikov, R. R. 2005, ApJ, 621, L69 [Google Scholar]
  244. Rameau, J., Chauvin, G., Lagrange, A. M., et al. 2013, ApJ, 772, L15 [Google Scholar]
  245. Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134 [Google Scholar]
  246. Raymond, S. N., Quinn, T., & Lunine, J. I. 2005, ApJ, 632, 670 [Google Scholar]
  247. Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265 [Google Scholar]
  248. Raymond, S. N., Armitage, P. J., & Gorelick, N. 2010, ApJ, 711, 772 [Google Scholar]
  249. Raymond, S. N., Boulet, T., Izidoro, A., Esteves, L., & Bitsch, B. 2018, MNRAS, 479, L81 [Google Scholar]
  250. Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680 [Google Scholar]
  251. Richert, A. J. W., Getman, K. V., Feigelson, E. D., et al. 2018, MNRAS, 477, 5191 [Google Scholar]
  252. Rogers, L. A. 2015, ApJ, 801, 41 [Google Scholar]
  253. Ronco, M. P., Guilera, O. M., & de Elía, G. C. 2017, MNRAS, 471, 2753 [Google Scholar]
  254. Rosotti, G. P., Ercolano, B., Owen, J. E., & Armitage, P. J. 2013, MNRAS, 430, 1392 [Google Scholar]
  255. Sahlmann, J., Ségransan, D., Queloz, D., et al. 2011, A&A, 525, A95 [Google Scholar]
  256. Santos, N. C., Adibekyan, V., Figueira, P., et al. 2017, A&A, 603, A30 [Google Scholar]
  257. Sanz-Serna, J. M. 1992, Acta Numerica, 1, 243 [Google Scholar]
  258. Sarkis, P., Mordasini, C., Henning, T., Marleau, G. D., & Mollière, P. 2021, A&A, 645, A79 [Google Scholar]
  259. Saumon, D., Chabrier, G., & van Horn H. M. 1995, ApJS, 99, 713 [Google Scholar]
  260. Savvidou, S., Bitsch, B., & Lambrechts, M. 2020, A&A, 640, A63 [Google Scholar]
  261. Schib, O., Mordasini, C., Wenger, N., Marleau, G. D., & Helled, R. 2021, A&A, 645, A43 [Google Scholar]
  262. Schlaufman, K. C. 2018, ApJ, 853, 37 [Google Scholar]
  263. Schlecker, M., Mordasini, C., Emsenhuber, A., et al. 2021a, A&A, 656, A71 (Paper III) [Google Scholar]
  264. Schlecker, M., Pham, D., Burn, R., et al. 2021b, A&A, 656, A73 (Paper V) [Google Scholar]
  265. Schlichting, H. E., & Mukhopadhyay, S. 2018, Space Sci. Rev., 214, 34 [Google Scholar]
  266. Schlichting, H. E., Fuentes, C. I., & Trilling, D. E. 2013, AJ, 146, 36 [Google Scholar]
  267. Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187 [Google Scholar]
  268. Seager, S., Kuchner, M., Hier-Majumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279 [Google Scholar]
  269. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
  270. Singer, K. N., McKinnon, W. B., Gladman, B., et al. 2019, Science, 363, 955 [Google Scholar]
  271. Sotiriadis, S., Libert, A.-S., Bitsch, B., & Crida, A. 2017, A&A, 598, A70 [Google Scholar]
  272. Stevenson, D. J. 1982, Planet. Space Sci., 30, 755 [Google Scholar]
  273. Stoer, J., & Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag) [Google Scholar]
  274. Suzuki, D., Bennett, D. P., Ida, S., et al. 2018, ApJ, 869, L34 [Google Scholar]
  275. Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64 [Google Scholar]
  276. Tanaka, H., & Ida, S. 1999, Icarus, 139, 350 [Google Scholar]
  277. Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [Google Scholar]
  278. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [Google Scholar]
  279. Thiabaud, A., Marboeuf, U., Alibert, Y., et al. 2014, A&A, 562, A27 [Google Scholar]
  280. Thiabaud, A., Marboeuf, U., Alibert, Y., Leya, I., & Mezger, K. 2015, A&A, 574, A138 [Google Scholar]
  281. Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 161, 431 [Google Scholar]
  282. Thorngren, D. P., & Fortney, J. J. 2018, AJ, 155, 214 [Google Scholar]
  283. Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64 [Google Scholar]
  284. Touboul, M., Kleine, T., Bourdon, B., Palme, H., & Wieler, R. 2007, Nature, 450, 1206 [Google Scholar]
  285. Valletta, C., & Helled, R. 2020, ApJ, 900, 133 [Google Scholar]
  286. Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283 [Google Scholar]
  287. Venturini, J., & Helled, R. 2020, A&A, 634, A31 [Google Scholar]
  288. Venturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90 [Google Scholar]
  289. Veras, D., & Armitage, P. J. 2004, MNRAS, 347, 613 [Google Scholar]
  290. Vigan, A., Bonavita, M., Biller, B., et al. 2017, A&A, 603, A3 [Google Scholar]
  291. Voelkel, O., Klahr, H., Mordasini, C., Emsenhuber, A., & Lenz, C. 2020, A&A, 642, A75 [Google Scholar]
  292. Voelkel, O., Deienno, R., Kretke, K., & Klahr, H. 2021, A&A, 645, A132 [Google Scholar]
  293. Wagner, K., Apai, D., & Kratter, K. M. 2019, ApJ, 877, 46 [Google Scholar]
  294. Wang, H., Weiss, B. P., Bai, X.-N., et al. 2017, Science, 355, 623 [Google Scholar]
  295. Ward, W. R. 1986, Icarus, 67, 164 [Google Scholar]
  296. Ward, W. R. 1989, ApJ, 345, L99 [Google Scholar]
  297. Ward, W. R. 1997, Icarus, 126, 261 [Google Scholar]
  298. Wasserburg, G. J., MacDonald, G. J. F., Hoyle, F., & Fowler, W. A. 1964, Science, 143, 465 [Google Scholar]
  299. Wei, Q., Hu, Y., Liu, Y., et al. 2018, ApJ, 856, L14 [Google Scholar]
  300. Weidenschilling, S. J. 1977, Ap&SS, 51, 153 [Google Scholar]
  301. Weidenschilling, S. J. 2005, Space Sci. Rev., 116, 53 [Google Scholar]
  302. Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6 [Google Scholar]
  303. Weiss, L. M., & Petigura, E. A. 2020, ApJ, 893, L1 [Google Scholar]
  304. Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [Google Scholar]
  305. William, L. 2007, Fundamentals of Geophysics (Cambridge: Cambridge University Press) [Google Scholar]
  306. Winter, A. J., Kruijssen, J. M. D., Chevance, M., Keller, B. W., & Longmore, S. N. 2020, MNRAS, 491, 903 [Google Scholar]
  307. Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proc. Natl. Acad. Sci. U.S.A., 113, 11431 [Google Scholar]
  308. Yin, Q., Jacobsen, S. B., Yamashita, K., et al. 2002, Nature, 418, 949 [Google Scholar]
  309. Youdin, A. N. 2011, ApJ, 742, 38 [Google Scholar]
  310. Zheng, X., Lin, D. N. C., & Kouwenhoven, M. B. N. 2017, ApJ, 836, 207 [Google Scholar]
  311. Zhou, J.-L., & Lin, D. N. C. 2007, ApJ, 666, 447 [Google Scholar]
  312. Zhu, W. 2020, AJ, 159, 188 [Google Scholar]

All Tables

Table 1

Initial conditions and parameters for the example system.

All Figures

thumbnail Fig. 1

Overview of the physical mechanisms included in the Bern model. At the top, the processes and base assumptions made in all model generations are given. The four boxes below show the four model generations with the main paper introducing each generation. The further points are: (1) number of initial embryos per disc, N-body integrator type, initial embryo mass, (2) phases simulated, (3) planetesimal accretion mode and planetesimal size, (4) phases withcalculation of the planets’ internal structure, (5) disc model characteristics, (6) orbital migration: type I, type II, transition criterion from type I to II (here thermal refers to a criterion only with the ratio between the Hill radius and the scale height of the disc, while ‘thermal and viscous’ refers to the full criterion of Crida et al. 2006), (7) disc-limited gas accretion rate, (8) later additions and improvements, (9) additional output relative to older generation, (10) population synthesis publications using this generation. In the bottom right panel, text in italic indicates new elements.

In the text
thumbnail Fig. 2

Sub-modules and most important exchanged quantities of the Generation III Bern model. The colours denote the stages at which processes are considered. Blue indicates processes active in the formation stage, but only before the dispersal of the gas disc. Green processes are considered during the entire formation stage, even after the dispersal of the gas disc. Processes in red are only considered during the evolution stage. The processes in black are included in all stages.

In the text
thumbnail Fig. 3

Time evolution of the surface density (top left), opacity κ (top right), midplane temperature (bottom left), and radial flow rate (bottom right) of a protoplanetary disc. The lines represent each one snapshot the state, and are spaced by about 2 × 105 yr. The blue line in both panels shows the initial profile, which has not yet been evolved at all, and is therefore not in equilibrium. The green line in the temperature profile shows the profile at disc’s dispersal, which is given by the equilibrium temperature with the host star’s luminosity.

In the text
thumbnail Fig. 4

Time evolution of the RMS of planetesimals’ eccentricity (top left), inclination (top right), and surface density (bottom left) of a circumstellar disc that also contains 10 embryos. The lines represent temporal snapshots of the three quantities, and are spaced by about 2 × 105 yr. The blue line in both panels denote the initial profile. The dashed vertical lines represent the location of the embryos, which is fixed in this case. N-body interactions were also disabled. The lifetime of the gas disc is shorter than the case presented in Fig. 3 due to the accretion by the protoplanets.

In the text
thumbnail Fig. 5

Comparison of two prescriptions for the maximum (i.e. disc limited) gas accretion rate, the one presented in this work (labelled ‘Bondi rate’) with that of Mordasini et al. (2012c) (labelled ‘Flow rate’). These are two different simulations (one for each prescription) whose initial conditions represent the second outermost planet in Fig. 4. Top panel: maximum value that can be supplied by the gas disc (labelled ‘Max.’) and effective accretion rates (labelled ‘Eff.’), which is given by intrinsic cooling in the initial attached phase and the maximum rate in the detached phase. Bottom panel: corresponding enveloppe mass (i.e. total gas accreted).

In the text
thumbnail Fig. 6

Snapshots of the internal structure of the second outermost planet of Fig. 4. The structures are split according tothe phases, with attached (top left), transition (the initial stage of the detached phase; top right), detached (bottom left) and evolutionary (bottom right). The red line shows the first profile of the detached phase and is shown of both panels. The green and blue profiles lie at the transition between two stages and are shown of two panels each. In each profile, thin lines show the part where energy transport is radiative and thick lines for convective.

In the text
thumbnail Fig. 7

Time evolution of the planet’s radius and luminosity of the second outermost planet of Fig. 4; the same as in Fig. 6. The insert on the left panel shows the contraction at the transition between attached and detached phase. The exact time where the model switches between the two phases is shown with the vertical dashed line.

In the text
thumbnail Fig. 8

Evolution of the planet’s core, planetesimals capture, and total radii as function of time (left panel) and core mass (right panel) for the second outermost planet of Fig. 4. In the right panel, we also include two time indicators at 105 and 106 yr with dashedvertical lines. The upper-right portion of the plot for the total radius in the left panel is also shown on the left panel of Fig. 7.

In the text
thumbnail Fig. 9

Illustration of the procedure to separate planetesimals’ feeding zones when zones would otherwise overlap. The horizontal axis represents the separation to the central star and four planets are shown. The light colour areas below the horizontal line show the initial feeding zones while the ones above show the final zones. amid2,3 and amid3,4 are the edges of the new feeding zone.

In the text
thumbnail Fig. 10

Radial surface density profile (top panel), temperature profile (middle panel), and migration map (bottom) as function of the planet mass (assuming zero eccentricity and inclination), for the same disc presented in Fig. 3 at t = 1 Myr. The value plotted in the bottom panel is the relative migration rate 1∕τa = −vplanetaplanet; blue regions indicate inward migration, red regions outward migration. For both directions, the locations in bright colours are where migration is inefficient while dark tones indicate efficient migration. The dashed black line shows the boundary between type I (below) and type II (above) migration regimes.

In the text
thumbnail Fig. 11

Cumulative distribution of the impactor-to-target mass ratio γ = Mtot,2Mtot,1 for different target mass ranges (as provided in the legend). Data come from the 100-embryos population presented in Paper II.

In the text
thumbnail Fig. 12

Average over the 10 simulations for each set of parameters. The top panels show the masses of solids (excepting the outer giant planets in the relevant cases) in the protoplanetary disc, that is still in planetesimals (solid lines), accreted by the embryos (dashed lines) and ejected (dotted lines). The dashed black line denotes the total mass of solids in each simulation. The bottom panel shows the number of embryos that remain. The ‘J/S’ in the legend refers to Jupiter and Saturn.

In the text
thumbnail Fig. 13

Comparison of formation tracks for a system with a MMSN-like surface density of planetesimals, and with 9 embryos (left), 23 (centre), and 46 (right), and no outer giant planets. Each line represents one embryo. Top panels: mass versus semi-major axis; embryos start at the bottom and move upwards as they grow. The final positions of the remaining planets are shown by dots. The dashed black line denotes the isolation mass (Lissauer 1987). Middle panels: mass versus time; sudden increases in mass are due to embryo-embryo collisions. Bottom panels: semi-major axis versus time.

In the text
thumbnail Fig. 14

Stacked eccentricity versus distance snapshots of 10 simulations with each 46 embryos initially. The left column shows the runs with outer giant planets whereas the right column has no outer giant planets. In each column, the 10 systems are represented with a different colour for each one. The bodies are shown by points whose sizes are proportional to their physical ones. Black crosses show the solar system planets.

In the text
thumbnail Fig. 15

Same as Fig. 14, but showing mass versus distance stacked diagram of 10 simulations with each 46 embryos initially. As before, the black crosses show the solar system planets.

In the text
thumbnail Fig. 16

Formation and evolution tracks of four giant planets with final masses in the 1∕3 to 2 M range in discs with a single embryo. Top panels: formation tracks with total mass Mtot versus distance (time goes towards the top) and total mass Mtot (solid lines) and core mass Mcore (dashed lines) versus time. Three panels on the bottom row: time dependence of the outer luminosity Ltot (bottom left), the distance (bottom centre) and the total radius Rtot (bottom right). For all panels except for the mass versus time (top right), the line styles denote the phase: dashed lines for the attached phase, solid line for the detached phase during formation and dash-dotted lines for the evolution stage. Line widths denote the migration regime, with tick lines for Type I and think lines for Type II. The legend in the top right panel shows the gas (in Solar masses) and planetesimals (in Earth mass) disc masses.

In the text
thumbnail Fig. 17

Formation and evolution tracks of two groups of four giant planets with final masses between 2 and 10 M in discs with a single embryo. Panel and line descriptions are the same as Fig. 16.

In the text
thumbnail Fig. 18

Formation and evolution tracks of one giant planet that ends up being accreted by the star during the evolution stage due to tidal migration. Top left: total mass Mtot versus distance; top right: total mass Mtot (black) and core mass Mcore (red) versus time; bottom left: total luminosity Ltot (black) and the bloating contribution Lbloat (red) versus time; bottom centre: distance versus time; bottom right: total radius Rtot (black) and outer (attached phase) or Hill (detached phase) radius (red) versus time.

In the text
thumbnail Fig. 19

Example of the formation of a planetary system from initially 100 lunar-mass embryos in low gas mass initial mass 0.023 M), low metallicity ([Fe/H] = −0.13) disc. The initial mass of planetesimals is 65.1 M. Four moments in time (in years) are shown. Lines shows the growth tracks in the semi-major axis-mass plane. Black points show (proto)planets existing at a given epoch. Grey open circles show the last position of protoplanets that were accreted by another more massive body in a giant impact. The colours of the lines are |τmigτgrow| = |d lnm∕d lna|.

In the text
thumbnail Fig. 20

Same system as in Fig. 19, but now showing the semi-major axes a of the planets as a function of time, colour coding in panel a the planets’ mass, in (b) the planetesimal surface density in the planets’ feeding zone, and in (c) the local gas surface density. Here, the vertical line indicates the moment of gas disc dissipation. Panel d: mass as a function of time, colour coding the semi-major axis. Small black circles indicate giant impacts, by showing the position or mass of the target (the more massive collision partner) at the moment of the impact.

In the text
thumbnail Fig. 21

Temporal evolution of the eccentricities of the planets of the system emerging in the low-mass disc shown in Fig. 19. Colours indicate the planet mass. For better visibility, only planets more massive than 0.1 M are shown. The curves are running averages such that one sees more clearly the mean values instead of rapid variations of the eccentricities. The thick black line is the mass of the gas disc relative to the value at 105 yr, which is in turn very similar to the initial value. The increase of the eccentricities at around 5 Myr when the gas disc dissipates is visible.

In the text
thumbnail Fig. 22

Example of the formation of a planetary system from initially 100 lunar-mass embryos in a high gas mass (initial mass 0.066 M), high metallicity ([Fe/H] = 0.23) disc. The initial mass of planetesimals is 432 M. The plot is analogous to Fig. 19, but the y-axis now extends to much higher masses, and the moments in time that are shown are different. At the end of the simulation at 20 Myr, this system contains one close-in sub-Neptunian planet, three giant planets, and a group of outer very low-mass planets.

In the text
thumbnail Fig. 23

Temporal evolution of the system shown in Fig. 22. The four panels have the same meaning as in Fig. 20.

In the text
thumbnail Fig. 24

Mass-distance diagrams of specific systems with 100 embryos initially (panels a–w), which are taken from the nominal population predicted for a 1 M star (NG76). Symbols are as follows: red points show gas-rich planets where MenvMcore > 1. Blue symbols are planets that have accreted some volatile material (ices) outside of the ice line(s). Green symbols are planets that have only accreted refractory solids. Open green and blue circles have 0.1envMcore ≤ 1 while filled green points and blue crosses have MenvMcore ≤ 0.1. For all these bodies, the grey horizontal bars go from ae to a + e. The top leftpanel with black crosses shows the solar system. Bodies lost because of collisions or ejections are shown in light grey. Planets accreted by the central star are show in the very left of each panel, the ejected ones on the very right and planets that collided with another (more massive) planet are shown at their last position on the diagram. The dotted vertical line in each system shows the location of the ice line. The number after each panel name is the metallicity [M/H] of the system expressed in dex, while the value on the top right is the initial mass of the planetesimals disc.

In the text
thumbnail Fig. 25

Radius-distance diagrams for a subset of the systems shown in Fig. 24. Here, we focus on systems with multiple low-mass planets; panel letters correspond to the same system in Fig. 24. In contrast to Fig. 24, lost planets are not shown for clarity.

In the text
thumbnail Fig. A.1

Mass-distance diagram of the nominal synthetic population NG76 of solar-like starts which stars with initially 100 moon-mass embryos per disc (see Paper II). The epochs of 0.1, 0.5, and 1 Myr are shown. Coloured points show protoplanets that can no longer accrete planetesimals of the initial local reservoir, which have a planetesimal accretion timescale of less than 3 Myr, and which are still embedded in the parent gaseous disc. When , the planetesimal accretion of these planets could in principle be affected by shepherding if they would be the only protoplanets growing in the disc.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.