The New Generation Planetary Population Synthesis (NGPPS). I. Bern global model of planet formation and evolution, model tests, and emerging planetary systems

Aims. Comparing theoretical models with observations allows one to make key step forward towards an understanding of planetary systems. It however requires a model able to (i) predict all the necessary observable quantities (not only masses and orbits, but also radii, luminosities, magnitudes, or evaporation rates) and (ii) address the large range in relevant planetary masses (from Mars mass to super-Jupiters) and distances (from stellar-grazing to wide orbits). Methods. We have developed a combined global end-to-end planetary formation and evolution model, the Generation III Bern model, based on the core accretion paradigm. This model solves as directly as possible the underlying differential equations for the structure and evolution of the gas disc, the dynamical state of the planetesimals, the internal structure of the planets yielding their planetesimal and gas accretion rates, disc-driven orbital migration, and the gravitational interaction of concurrently forming planets via a full N-body calculation. Importantly, the model also follows the long-term evolution of the planets on Gigayear timescales after formation including the effects of cooling and contraction, atmospheric escape, bloating, and stellar tides. Results. To test the model, we compared it with classical scenarios of Solar System formation. For the terrestrial planets, we find that we obtain a giant impact phase provided enough embryos (~100) are initially emplaced in the disc. For the giant planets, we find that Jupiter-mass planets must accrete their core shortly before the dispersal of the gas disc to prevent strong inward migration that would bring them to the inner edge of the disc. Conclusions. The model can form planetary systems with a wide range of properties. We find that systems with only terrestrial planets are often well-ordered while giant-planet bearing systems show no such similarity.


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(Boss , 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 deuteriumburning limit (Kratter et al. 2010). On the other extreme, for very close-in giant planets, core accretion (Perri & Cameron 1974;Mizuno 1980)  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;. 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(Baruteau et al. , 2016. The final mass and location of 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;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 implies that 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 in-stance, 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  or after the dispersal by dynamical instabilities (Inamdar & Schlichting 2016;Izidoro et al. 2017Izidoro et al. , 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).
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. (in rev., 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.

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. (2004Alibert et al. ( , 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 longterm evolution of the formed planet (Generation Ib; Mordasini et al. 2012c,b) while the other obtained the ability to form multiplanetary 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 same time, 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.

General description
We base our study 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 Section 6 regarding the impact of this specific choice). Afterwards, in the evolu-tionary phase, we follow the thermodynamical evolution of each planet individually to 10 Gyr.

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 (Section 3). These serve as sources for the accretion of the protoplanets (Section 4). The lifetime of 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 (Section 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).

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 selfconsistently, 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.

Envelope structure
The calculation of the internal structure of all planets (Section 4) during their entire formation and evolution is a crucial aspect of the Bern Model, as visible from its central position in Figure  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 A&A proofs: manuscript no. model Physical mechanisms and base assumptions included in all model generations -Formation paradigm: core accretion -Protoplanetary disc model: solution of 1D evolution equation for gas surface density in an axissymetric constant α-disc with photoevaporation -Solid accretion: rate equation (Safronov-type) from planetesimals of a single size; planetesimals are represented by a solid surface density with dynamical state -Gas accretion and planet interior structure: from solving 1D radially symmetric hydrostatic planet interior structure equations -Orbital migration: gas disc-driven, types I and II Evolution of physical mechanisms considered in various model generations Generation I (Alibert et al. 2005a): base model 1. 1 embryo per disc (no N-body), 0.6 M ⊕ 2. Formation only (to t disc ) 3. Runaway planetesimals accretion, 100 km 4. Attached phase only 5. Vertical disc structure, no stellar irradiation, no stellar evolution 6. Isothermal type I, equilibrium type II, thermal only transition criterion 7. Equilibrium gas flux in disc 8. Stellar irradiation of the disc (Fouchet et al. 2012) 9. Masses, orbital distances, bulk composition 10. Mordasini et al. (2009a,b); Alibert et al. (2011);Mordasini et al. (2012a) Generation Ib (Mordasini et al. 2012c,b): inclusion of longterm evolution 1. 1 embryo per disc (no N-body), 0.6 M ⊕ 2. Formation (to t disc ) and (thermodynamic) evolution (to 10 Gyr) 3. Runaway planetesimals accretion, 100 km 4. Attached, detached and evolutionary, with core structure 5. Vertical disc structure, with stellar irradiation, no stellar evolution 6. According to Dittkrist et al. (2014): Non-isothermal type I, non-equilibrium type II, thermal and viscous transition criterion 7. Non-equilibrium gas flux in disc 8. D-burning, atmospheric escape 9. Radii, luminosities, envelope evaporation rates 10. Mordasini et al. (2014Mordasini et al. ( , 2017; Jin & Mordasini (2018) Generation II : inclusion of N-body interaction 1. Several embryos per disc (EMPS N-body integrator), 0.01 M ⊕ 2. Formation only (to t disc ) 3. Oligarchic planetesimals accretion, 300 m  4. Attached phase only 5. Vertical disc structure, no stellar irradiation, no stellar evolution 6. According to Dittkrist et al. (2014): Non-isothermal type I, non-equilibrium type II, thermal and viscous transition criterion 7. Non-equilibrium gas flux in disc 8. Composition tracking (Thiabaud et al. 2015) 9. Multiplicity, eccentricities, MMR 10. Pfyffer et al. (2015); Alibert & Benz (2017) Generation III (this work): long-term evolution and N-body evolution 1. Many embryos per disc (Mercury N-body integrator), 0.01 M ⊕ 2. Formation (to 20 Myr) and (thermodynamic) evolution (to 10 Gyr) 3. Oligarchic planetesimals accretion, 300 m 4. Attached, detached, evolutionary (with D-burning, escape, bloating, Roche lobe overflow, core structure) 5. Vertically integrated, with stellar irradiation and stellar evolution 6. Non-isothermal type I, non-equilibrium type II, thermal and viscous transition criterion 7. Bondi-limited gas accretion 8. None 9. Combines output of Ib and II 10. This work (NGPPS series) (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 with calculation 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.
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, we now use self-consistently the iron mass fraction as given by the disc compositional model (according to Thiabaud et al. 2014, Section 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 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.
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 Figure 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 submodules visible in Figure 2.
Article number, page 5 of 45 A&A proofs: manuscript no. model

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 Nbody 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).

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), where Σ g = ∞ −∞ ρdz is the surface density of gas, andΣ g,photo andΣ g,planet are the sink terms related to photo-evaporation (Section 3.2.2) and accretion by the planets respectively. The viscosity is parametrised, following Shakura & Sunyaev (1973), with This equation is solved on a grid spaced regularly in log with 3400 points that extends from the inner location of the disc r in (an initial condition) to r max = 1000 au. At these two locations, the surface density is fixed to zero.

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 with T mid the disc mid-plane temperature, T s 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 , T mid )Σ where ρ mid = Σ/( √ 2πH) is the central density, H = c s /Ω the disc's vertical scale height, c s = k B T mid /(µm H ) the isothermal sound speed, µ = 2.24 the mean molecular weight of the gas, and m H 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-togas ratio chosen for the solids disc) and Freedman et al. (2014) (which gives molecular 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), orin 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) dustpebble 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 bẏ 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 Ω = GM /r 3 , G being the gravitational constant. The selfgravity 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 H/∂ ln r = 9/7 (Chiang & Goldreich 1997). The T irr term accounts for the direct irradiation through the disc midplane. It is computed as which is the black-body equilibrium temperature accounting for the optical depth through the midplane of the disc τ mid = ρ mid κ(ρ mid , T mid )dr. This contribution is 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 Sec. 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 T cd = 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).

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 T I ≈ 10 3 K. The corresponding sound speed is then 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 We assume that mass is removed uniformly outside of β I r g,I with β I = 0.14 (similar to Alexander & Pascucci 2012), so that the rate is given bẏ withṀ wind a parameter that provides the total mass loss rate if the disc would extend to r max = 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 r max , 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 T II ≈ 10 4 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 r 14 = β II r g,II /10 14 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 n 0 (r 14 ) = k Hol Φ 1/2 41 r −3/2 14 , where we set k Hol = 5.7 × 10 4 following the hydrodynamical simulations of Hollenbach et al. (1994) and Φ 41 = 0.1 √ M /M , which is the ionising photon luminosity in the units of 10 41 s −1 . The distance-dependent base density can then be calculated as n 0 (r) = n 0 (r 14 ) r r g,II − 5 2 . (11) We further follow Clarke et al. (2001) to getΣ g,photo,int = 2c s,II n 0 m H outside of β II r g,II .
The final photoevaporation rate is given by the sum of the effects of host star + nearby massive stars witḣ Σ g,photo =Σ g,photo,ext +Σ g,photo,int .
3.2.3. Initial gas surface density profile and example We initialise the gas surface density profile with (Veras & Armitage 2004) where r 0 = 5.2 au is the reference distance, β g = 0.9 the powerlaw index (Andrews et al. 2010), r cut,g the characteristic radius for the exponential decay and r in the inner edge of the disc. The conversion between the total mass and the normalisation surface density Σ g,0 at r 0 is obtained with It should be noted that this formula neglects the lack of gas within r in , 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.Σ g,planet = 0), 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 × 10 6 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 by the protoplanets, which means we have the radial flow rate varying with distance (bottom right panel of Fig. 3, where the radially constant 10 −1 10 0 10 1 10 2 10 3 Distance [AU]   of a protoplanetary disc. The lines represent each one snapshot the state, and are spaced by about 2 × 10 5 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. 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).

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 (e plan ) and inclination (i plan ) as dynamical state.

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 e plan and inclination i plan read aṡ 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 only Table 1. Initial conditions and parameters for the example system. The upper part contains the gas disc properties, the middle part the planetesimals disc properties, and the bottom part show planetary embryos properties.

Quantity
Value Stellar mass M 1 M Reference surface density Σ g,0 at 5.2 au 145 g cm −2 Initial gas disc mass M g 3.90 × 10 −2 M Inner edge of the gas disc r in 0.091 au (10 d) Characteristic radius of the gas disc r cut,g 66.5 au Disc viscosity parameter α 2 × 10 −3 External photoevaporation rateṀ wind 6.42 × 10 −7 M /yr Power law index of the gas disc β g 0.9 Dust-to-gas ratio 3.4 % Planetesimal disc mass 348 M ⊕ Power law index of the solids disc β s 1.5 Characteristic radius of the solids disc r cut,s r cut,g /2 Planetesimal radius 300 m Planetesimal density (rocky) 3.2 g cm −3 Planetesimal density (icy) 1 g cm −3 Embryo mass M emb,0 1 × 10 −2 M ⊕ Opacity reduction factor f opa 3 × 10 −3 evaluated 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 Re mol = v rel R plan /ν mol , where ν mol = λc s /3 is the molecular viscosity, λ = (n H 2 σ H 2 ) −1 the gas molecules' mean free path, n H 2 the number density assuming all of the gaseous molecules having hydrogen mass, σ H 2 their collisional cross-section, R plan the planetesimals' radius, v rel = v K η 2 + 5/8e 2 plan + 1/2i 2 plan (17) their relative velocity, 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 v K = Ωr the Keplerian velocity. When R plan < λ, the gas drag is assumed to be in the Epstein regime. Otherwise, if Re mol > 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 is the gas drag time scale and C D = 1. In the Stokes regimes the drag expressions arė c s ρ mid ρ plan R plan (25) (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 aṡ where the sum is over all the protoplanets present in the system, is the modulation due to separation so that the perturbation is effectively restricted to the planet's feeding zone, the planet's Hills radius, and b = 5 is the half-width of the feeding zone (see Sect. 4.3.3). The terms P VS and Q VS are given by , Here,ẽ plan = re plan /R H andĩ plan = ri plan /R H are respectively the reduced planetesimals' eccentricity and inclination, Λ = ı plan (ẽ 2 plan +ĩ 2 plan )/12, β = i plan /e plan , while for I PVS and I QVS we use the approximations obtained by Chambers (2006): The stirring by the other planetesimals is given by, following Ohtsuki et al. (2002), and M plan = 4/3πR 3 plan ρ plan , 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 e plan = 2.31 M 4/15 plan r 1/5 ρ 2/15 plan Σ 1/5 and We also compared our prescription for the dynamical state with gamma-stirring from, for instance,  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 selfstirring 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.
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 ), 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 . 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 . 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 subject of 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, with the power-law exponent is set to β s = 1.5, as in the MMSN, and r cut,s = r cut,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 e 2 plan > 0.95. 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  for instance. Finally, the planetesimals disc remains after the dispersal of the gas disc; the only difference is that the damping terms for eccentricityė 2 plan drag and inclinationi 2 plan drag vanish.

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 f s (r) in Eq. (40) 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 contribution of molecules in the solid phase are accounted for the resulting solid surface density. Thus, the value of f s 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 H 2 O-ice line in all discs, which induces the largest surface density jump because H 2 O makes up ∼60 % of all ices in mass (Marboeuf et al. 2014a).

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 Nbody 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 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 β = i plan /e plan ≈ 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.

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 10 −1 10 0 10 1 10 2 10 3 Location [AU]  , 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 × 10 5 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. 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.

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 to 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 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 f opa = 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. where is the Bondi radius, R H is the Hill's radius (Eq. 29), k 1 = 1 and k 2 = 1/4. The pressure and temperature are derived from the local properties of the disc with P(R tot ) = P neb (a planet ) and (47) and 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(R tot ) is described in Section 4.2. In the case of the mass, what is known is the core mass, that is M(R core ) = M core , while M(R tot ) = M tot is the quantity that is being searched for. We thus use an iterative method by guessing M tot , which is then used to integrate the internal structure equations until the boundary condition at the inner boundary is fulfilled, that is M(R core ) = M core . Once M tot is found, the envelope mass can be retrieved by M env = M tot − M core , and the gas accretion rate by taking the difference of the envelope mass between two successive steps of the envelope structure calculatioṅ M env = (M env (t) − M env (t − ∆t))/∆t.

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, 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 R gc around the planet. For simplicity, we compute R gc according to Eq. (45). Depending on the value of R gc with respect to H, the local disc's scale height, two different regimes occur. In the case where R gc < 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 σ cross = πR 2 gc is given bẏ with ρ ≈ Σ/H the approximate density of the gas and v rel = max (ΩR tot , c s ) the relative velocity between the gas and the planet.
On the other hand, in the case R gc > 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 havė To distinguish between the two regimes, we use the lower rate of the two, that iṡ Finally, to ensure that no more gas than available in the feeding zone M feed is accreted during one time step, we further con-strainṀ env,max < M feed /∆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 f feed = 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 for the calculation of the maximum rate and the removal of the accreted gas, witḣ 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Σ g,planet 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  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).
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.

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 P(R tot ) = P neb (a planet ) + P edd + P ram + P rad (57) with P neb (a planet ) being the pressure at the midplane of the gas disc, P edd = (2g)/(3κ) the Eddington expression for the photospheric pressure due to the material residing above the τ = 2/3 surface, P rad = (2σ SB T 4 (R tot ))/(3c) the radiation pressure, c being the speed of light in vaccum, and being the ram pressure due to the accretion shock and the freefall velocity at the surface of the planet.

Evolutionary phase
For the evolutionary phase (after the dispersal of the gas disc), the outer boundary conditions are set to where T 4 int = L tot /(4πσ SB R 2 tot ) is the intrinsic temperature, T eq = T * R /(2 * a planet ), 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).

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. 46), 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.

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 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 −3/5GM 2 core /R core . 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. (61) represents as in polytropic models the mass distribution and additionally the thermal energy content. It is retrieved from Eq. (61). The internal luminosity resulting from the accretion, cooling, and contraction L int can then be obtained as 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  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. ξ tot . To circumvent this problem, we estimate the luminosity with The correction factor C corrects for neglecting theξ tot term. The value of C can be calculated a posteriori by determining the actual total energy of the new planet, with The value of C is then used for the next time step. Marleau et al. (2017Marleau et al. ( , 2019b conducted 1D radiationhydrodynamic 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 . 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 L int (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 L radio includes contributions from the three most important long-lived radionucleides 40 K, 238 U and 232 Th (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 Q 0 ≈ 5 × 10 −7 erg g −1 s −1 of mantle material (all elements besides iron).

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): F = L /(4πa 2 planet ) 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 F exp (−τ mid ) (before the dispersal of the gas disc) is greater than 2 × 10 8 erg cm −2 s −1 (Demory & Seager 2011).

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).

Total luminosity
The final luminosity is then given by 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.

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 iṡ withΣ s the mean surface density of planetesimals in the planet's feeding zone and p coll the collision probability 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.

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 = re plan /R H and inclinationĩ plan = ri plan /R H (r is the heliocentric distance): the high-velocity regime forẽ plan ,ĩ plan 2, midvelocity for 2 ẽ plan ,ĩ plan 0.2 and low-velocity for 0.2 ẽ plan ,ĩ plan . According to Inaba et al. (2001), each regime has a different expression for the collision probability, where R cap is the planetesimal capture radius of the planet, β = i plan /e plan and the I F and I G functions can be approximated as, following Chambers (2006): The final collision is then given by In the initial stage, the capture radius R cap is the physical radius of the core R core . 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 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 10 4 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 × 10 5 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.

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, v esc = 2GM /a planet , will likely be ejected from the system. Thus, we have that the fraction of accreted-to-ejected planetesimals is (Ida & Lin 2004a) Region 1 Region 2

Before separation
After separation Central star 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. a mid2,3 and a mid3,4 are the edges of the new feeding zone.

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 For a planet on a circular orbit, conservation of the Jacobi energy implies that b = 12 + 4/3(ẽ 2 plan +ĩ 2 plan ) (e.g. Hayashi et al. 1977). So, in a quiescent disc withẽ plan ,ĩ plan 1, b = 2 √ 3 ≈ 3.5. 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 r peri − R feed to r apo + R feed , with r peri and r apo 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 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 midpoint between the two planets and found that the prescription does not significantly affect the outcomes of the simulations.

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 for the 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. (41) and (42), but with an equation of state that takes the form of a modified polytrope from Seager et al. (2007), which reads ρ(P) = ρ 0 + cP n .
We include three different materials: iron, silicates (perovskite, MgSiO 3 ) 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, be small (Grasset et al. 2009). For gas giant planets, where envelopes can reach masses of thousands of Earth masses, this can cause a significant compression of 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).

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 energylimited 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 timedependent 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 bẏ where F EUV is the EUV flux, R base the radius of the photoionisation base, calculated as in Murray On the other hand, energy-limited evaporation is not suitable when the EUV flux is high (> 10 4 erg 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 escapė M env,rr ∼ 4πρ s c s R 2 s (81) at the sonic point R s , which is calculated the same way as R acc .
Here c s = √ k B T/(m H /2) is the isothermal sound speed of ionised gas with T = 10 4 K. The density can be related to the one at the ionisation base, where τ = 1, with The photoionisation base is located where there is equilibrium between photoionisations and recombination: with n 0,base the density of neutrals at the base, hν 0 = 20 eV, σ ν 0 = 6 × 10 −18 cm 2 (hν 0 /13.6 eV) −3 , α rec = 2.7 × 10 −13 , and ρ base = n +,base m H . 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 in an 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.

Initial conditions
The simulation begin with a predetermined number of embryos whose initial mass is M emb,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 r in 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.

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.

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.

Type I migration
For Type I migration, 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) to (53) of Paardekooper et al. (2011) and (15) of Coleman & Nelson (2014), is given by 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) to (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. (84) account for the shape of the orbit. F L provides the reduction of the Lindblad torque for eccentric or inclined orbits following Cresswell & Nelson (2008), with and Here,ê = e/h = e/(H/r) andî = i/h = i/(H/r) are the planet's orbital eccentricity and inclination scaled by the disc's aspect ratio h = H/r. F e and F i provide the reduction of the corotation torques due to eccentricity and inclination (Bitsch & Kley 2010). We use as suggested by Fendyke & Nelson (2014) for the reduction due to eccentricity and for the reduction due to inclination (Coleman & Nelson 2014). Eccentricity and inclination damping time scales follow Cresswell & Nelson (2008), with τ e = t wave 0.78 1 − 0.14ê 2 + 0.06ê 3 + 0.18êî 2 and where is the characteristic time of evolution of density waves (Tanaka & Ward 2004).

Type II migration
The criterion to detect gap opening and switch migration to Type II is from Crida et al. (2006), 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, (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 v planet is given by 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 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).
Article number, page 21 of 45 A&A proofs: manuscript no. model During type II migration, the eccentricity and inclination damping time scales are set to τ e = τ i = 1 10 |τ a | = 1 10 a planet |v planet | .
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).

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 −2 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 . 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 ≈ 8 M ⊕ 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 ≈ 2 M ⊕ could reach the inner edge of the disc.

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, where x denotes the position coordinates, p the momentum coordinates, and is the Hamiltonian of the system, with ∆x i j = |x i − x j |. Here, the index i = 0 refers to the central star and M 0 = M while the subsequent are the planet with M i = M planet,i so that N is the number of planets in the system. However, while H 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 10 −1 10 0 10 1 10 2 Distance [AU] into three components, so that H = H K + H S + H I , and Here, H K represents the unperturbed Keplerian orbits of the planets about the central star, H S the kinetic energy of the star and H I 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 H K H S , H I , unless two planets come close together.
The evolution of such a system by splitting is done using a second-order method, where the notation H ... (τ) is used to represent the evolution under the given Hamiltonian for a step τ. For H I , 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, H I is extended to include additional forces representing the effect of the gas disc, see Section 5.2.1. The evolution under H S results in a shift τ/(2M ) p i while the evolution under H K is a Keplerian motion around the central star for a period τ.
As we noted, the assumption that H I is small compared to H K 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 H K so that the interaction Hamiltonian remains small. This implies that H K 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-by bodies 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 O(N 2 ), 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.

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, v K = Ωr is the Kelperian velocity.

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. (45), 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 R cap ; 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.

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 Section 4.2; so the additional energy is calculated using where µ = M tot,1 M tot,2 /(M tot,1 + M tot,2 ) is the reduced mass, and the indexes 1 and 2 refer to the quantities of the larger and smaller body respectively. v imp is the relative velocity at time of contact. Here E acc,core = G M tot,1 M core,2 R core,1 + R core,2 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 colliding at 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 viȧ where t impact is the time of the impact, τ impact = 10 4 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.

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 that are 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 -Mello et al. 2008;Jackson et al. 2009;Benítez-Llambay et al. 2011), where Q = 10 6 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).

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 ) 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.

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 r 0 = 1 au, but truncated at 2 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. (in rev., 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.01 M ⊕ ) 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 2001Chambers , 2013Emsenhuber 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.

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 embryos by 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.

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 yet more interactions. This feedback only ends when nearly all planetesimals have been accreted onto the embryos.  (Lissauer 1987). Middle panels: mass versus time; sudden increases in mass are due to embryo-embryo collisions. Bottom panels: semi-major axis versus time.
Thus, closer packed embryos lead to enhanced stirring of their eccentricities, which has two consequences: the increase of the feeding zone size 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 are capable 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 counterintuitive at first: the larger the number of embryos, the less planets remain. We observe this for instance in the bottom panel of Fig. 12.

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 10 5 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 1 au 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.

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 Fig. 14 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 than 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 moment the 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 selfstirring 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, 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. Raymond et al. (2006), which prevents the accretion into a lower number of higher mass bodies (Kokubo et al. 2006).

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 (Section 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 lowmass, 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 < 1 au 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.

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 gas accretion, 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 r cut,g of the gas disc is set as M g 2 × 10 −3 M = r cut,g 10 au (see Paper II for the motivation). We provide the remaining two parameters, the initial masses of the gas and planetesimals discs in the following.

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 to 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 planetdominated 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 . 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 oneplanet-per-disc approximation studied here than what was found by some other models using the in situ (and one-embryo-perdisc) approximation. For the latter, the favoured scenario is that a core between 10 and 20 M ⊕ forms early (less than 10 5 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 is included; 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 ) 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 Jupitermass 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).  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.

More massive planets
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. The slope 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.

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 consequence is 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. (112)). To determine the reason for the luminosity increase, we print alongside the total value, the contribution from bloating (Eq. (64)). 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.

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.

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. 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).

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 a − M plane the tracks of the planets by the ratio |τ mig /τ grow | = |d ln m/d ln a|. Regarding the timescales, it is of fundamental importance that the oli-garchic 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 (10 5 yr, top left panel of Fig. 19), the quasi in-situ accretion of planetesimals present in the initial)feeding zone of the embryos is 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 10 5 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 10 5 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 A&A proofs: manuscript no. model 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 10 5 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- 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) shows 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. 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 to 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 in-Article number, page 33 of 45  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 10 5 years, 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.
ner 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 Figure 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 disc dissipates 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 consequence of 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 10 4 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 Nbody 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 volatilerich 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.

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 a − M 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 10 5 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 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 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 systemwide 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 ex-A&A proofs: manuscript no. model ample 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-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.

Overview of the diversity of system architectures
Figures Fig. 24 and 25 show the mass-distance and radiusdistance 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. 87 and 88) 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 reassessment 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 eccen-tricity 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.

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 lowmass 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 to 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.

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 M iso ∝ r (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.

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 composition is 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. (in pressa, hereafter Paper III) will look thoroughly at correlations between close-in Super-Earth planets and long-period giants.

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 giantimpact 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 to 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 , and enables the comparison with directlyimaged 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 . 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 (Earthlike) 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. in press) 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. in pressb) will study the mapping of disc initial conditions to planet properties with machine learning. -Paper VI (Mishra et al. subm.) 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 multiaspect comparisons with observations will be helpful to improve the understanding of planet formation and evolution. Mass-distance diagrams of specific systems with 100 embryos initially (panels a to 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 M env /M core > 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.1≤ M env /M core ≤ 1 while filled green points and blue crosses have M env /M core ≤ 0.1. For all these bodies, the grey horizontal bars go from a − e to a + e. The top left panel 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. 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τ mig τ mig,c , the planetesimal accretion of these planets could in principle be affected by shepherding if they would be the only protoplanets growing in the disc.
toplanets which might at least in principle be affected by shepherding. 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τ mig to the normalised critical migration timescaleτ mig,c 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τ mig τ mig,c . At 0.5 Myr, there is a radial interval from about 1 to 4 au whereτ mig is longer thanτ mig,c , 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.