11institutetext:  Lunar and Planetary Laboratory, University of Arizona, 1629
E. University Blvd., Tucson, AZ 85721, USA
11email<EMAIL_ADDRESS>22institutetext: Physikalisches Institut,
Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland
33institutetext: Universitäts-Sternwarte München, Ludwig-Maximilians-
Universität München, Scheinerstraße 1, 81679 München, Germany 44institutetext:
Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

# The New Generation Planetary Population Synthesis (NGPPS)

I. Bern global model of planet formation and evolution, model tests,
and emerging planetary systems

Alexandre Emsenhuber The New Generation Planetary Population Synthesis
(NGPPS)The New Generation Planetary Population Synthesis (NGPPS)The New
Generation Planetary Population Synthesis (NGPPS)The New Generation Planetary
Population Synthesis (NGPPS)The New Generation Planetary Population Synthesis
(NGPPS)The New Generation Planetary Population Synthesis (NGPPS) Christoph
Mordasini The New Generation Planetary Population Synthesis (NGPPS)The New
Generation Planetary Population Synthesis (NGPPS) Remo Burn The New Generation
Planetary Population Synthesis (NGPPS)The New Generation Planetary Population
Synthesis (NGPPS)The New Generation Planetary Population Synthesis (NGPPS)The
New Generation Planetary Population Synthesis (NGPPS) Yann Alibert The New
Generation Planetary Population Synthesis (NGPPS)The New Generation Planetary
Population Synthesis (NGPPS) Willy Benz The New Generation Planetary
Population Synthesis (NGPPS)The New Generation Planetary Population Synthesis
(NGPPS) Erik Asphaug The New Generation Planetary Population Synthesis
(NGPPS)The New Generation Planetary Population Synthesis (NGPPS)

(Received 1 June 2020 / Accepted DD MMM YYYY)

###### Abstract

Context. The explosion of observational data on exoplanets gives many
constraints on theoretical models of planet formation and evolution.
Observational data probe very large areas of the parameter space and many
different planet properties.

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 of protoplanet-protoplanet collisions provided enough embryos
($\sim$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. Regarding the emergence of entire
planetary systems, many aspects can be understood with the comparison of the
timescales of growth and migration, the capture into resonances, and the
consequences of large scale dynamical instabilities caused by the
gravitational interactions of protoplanets, including the situation when a
second core starts runaway gas accretion.

Conclusions. The Generation III Bern model provides one of the most
comprehensive global end-to-end models of planetary system formation and
evolution developed so far, linking self-consistently a multitude of crucial
physical processes. The model can form planetary systems with a wide range of
properties. We find that systems with only terrestrial planets are often well-
ordered (in period, mass, and radius) while giant-planet bearing systems show
no such similarity. In a series of papers, the model will be used to perform
extensive planetary population syntheses, putting the current theoretical
understanding of planet formation and evolution to the observational test.

##  1 Introduction

Since the discovery of the first exoplanet detected around a main sequence
star (Mayor & Queloz 1995), the number of known exoplanets has exploded. These
planets span a wide range of masses and sizes, and 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 will first discuss
possible formation mechanisms for different planet kinds.

###  1.1 Giant planets

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 – well within the orbit
of Mercury (Mayor et al. 2011; Fabrycky et al. 2014); others were detected at
large separations using the direct imaging technique (Marois et al. 2008;
Lagrange et al. 2010; Rameau et al. 2013; Macintosh et al. 2015; Chauvin et
al. 2017; Keppler et al. 2018).

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

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

This process is further complicated by the implications of planetary migration
(Baruteau et al. 2014, 2016). The final mass and 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; Ford & Rasio 2008; Raymond et
al. 2010; Sotiriadis et al. 2017), and that most giant planet-harbouring
system are multiple (Knutson et al. 2014; Bryan et al. 2016; Wagner et al.
2019).

###  1.2 Low- and intermediate-mass planets

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 life time of the protoplanetary
disc. However, whether they formed early (in the same way as the core of giant
planets) or toward 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 life time of the protoplanetary disc and can form over
longer periods of time (Lambrechts et al. 2019; Ogihara et al. 2018). These
could also have had a gaseous envelope in the past that was removed by, for
instance, atmospheric escape (e.g., Jin et al. 2014) or giant impacts
(Schlichting & Mukhopadhyay 2018).

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

###  1.3 Constraining planetary formation

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

As the model have 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.

###  1.4 This work

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 e.g., Ida & Lin (2010) who introduced a semi-analytical approach,
an improvement over previous works, e.g. 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 asses 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.

##  2 The Bern model

###  2.1 History

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 $\alpha$-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\text{\,}\mathrm{M}_{\oplus}$ 2. Formation only (to
$t_{\mathrm{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 long-term evolution 1.
1 embryo per disc (no N-body), $0.6\text{\,}\mathrm{M}_{\oplus}$ 2. Formation
(to $t_{\mathrm{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-burining, atmospheric escape 9. Radii,
luminosities, envelope evaporation rates 10. Mordasini et al. (2014, 2017);
Jin & Mordasini (2018) Generation II (Alibert et al. 2013): inclusion of
N-body interaction 1. Several embryos per disc (EMPS N-body integrator),
$0.01\text{\,}\mathrm{M}_{\oplus}$ 2. Formation only (to $t_{\mathrm{disc}}$)
3. Oligarchic planetesimals accretion, 300 m (Fortier et al. 2013) 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\text{\,}\mathrm{M}_{\oplus}$ 2. Formation (to $20\text{\,}\mathrm{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)

Figure 1: Overview of the physical mechanisms included in the Bern model. At
the top, the processes and base assumptions made in all model generations are
given. The four boxes below show the four model generations with the main
paper introducing each generation. The further points are: (1) number of
initial embryos per disc, N-body integrator type, initial embryo mass, (2)
phases simulated, (3) planetesimal accretion mode and planetesimal size, (4)
phases 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.

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

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

###  2.2 General description

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

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\text{\,}\mathrm{Myr}$,
see the related discussion in Section 6 regarding the impact of this specific
choice). Afterwards, in the evolutionary phase, we follow the thermodynamical
evolution of each planet individually to $10\text{\,}\mathrm{Gyr}$.

####  2.2.1 Formation phase

During the formation stage (0–$20\text{\,}\mathrm{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 life time 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).

####  2.2.2 Evolution phase

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

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

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

####  2.2.3 Envelope structure

The calculation of the internal structure of all planets (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 depending of a planet’s
accretion history, a water ice layer and a gaseous envelope made of pure H/He.
In contrast to earlier syntheses predicting planetary radii (Mordasini et al.
2012c, 2014), we now use self-consistently the iron mass fraction as given by
the disc compositional model (according to Thiabaud et al. 2014, 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 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\text{\,}\mathrm{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}\text{\,}\mathrm{M}_{\oplus}$ 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. This
makes it in particular 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\text{\,}\mathrm{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\text{\,}\mathrm{Myr}$
when stars are in the saturated phase of high XUV emission, and not only for
the first $20\text{\,}\mathrm{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 (in particular also after gas disc dissipation) the
accretion of planetesimals, which also changes planet masses.

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

##  3 Star and protoplanetary disc

###  3.1 Stellar model

Instead of assuming a fixed $1\text{\,}\mathrm{L}_{\odot}$ stellar luminosity
for a $1\text{\,}\mathrm{M}_{\odot}$ 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_{\star}$,
luminosity $L_{\star}$ and temperature $T_{\star}$ for a given stellar mass
$M_{\star}$ at any moment. Stellar temperature and radius are used for the
outer boundary conditions of the gas disc; stellar radius is also used in the
N-body integrator to detect collisions with the star and to calculate the
stellar tidal migration. Finally, the stellar irradiation enters into the
calculation of the outer (atmospheric) temperature (at $\tau=2/3$) of the
planets’ interior structure as described in Mordasini et al. (2012c) and
radius bloating (Sect. 4.2.2).

###  3.2 Gas disc

The protoplanetary gas disc is modelled with a 1D radial axissymetric
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),

$\frac{\partial\Sigma_{\mathrm{g}}}{\partial
t}=\frac{1}{r}\frac{\partial}{\partial
r}\left[3r^{1/2}\frac{\partial}{\partial
r}\left(r^{1/2}\nu\Sigma_{\mathrm{g}}\right)\right]-\dot{\Sigma}_{\mathrm{g,photo}}-\dot{\Sigma}_{\mathrm{g,planet}},$
(1)

where $\Sigma_{\mathrm{g}}=\int_{-\infty}^{\infty}\rho\mathrm{d}z$ is the
surface density of gas, and $\dot{\Sigma}_{\mathrm{g,photo}}$ and
$\dot{\Sigma}_{\mathrm{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

$\nu=\alpha c_{s}H$ (2)

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

####  3.2.1 Vertical structure

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

$\sigma_{\mathrm{SB}}T_{\mathrm{mid}}^{4}=\frac{1}{2}\left(\frac{3}{8}\tau_{R}+\frac{1}{2\tau_{P}}\right)\dot{E}+\sigma_{\mathrm{SB}}T_{\mathrm{s}}^{4}$
(3)

with $T_{\mathrm{mid}}$ the disc mid-plane temperature, $T_{\mathrm{s}}$ the
temperature due to irradiation (see below), $\sigma_{\mathrm{SB}}$ the Stefan-
Boltzmann constant, $\tau_{R}$ and $\tau_{P}$ are the Rosseland and Planck
mean optical depths respectively, and $\dot{E}$ is the viscous dissipation
rate. This formula yields the mid-plane temperature both in the optically-
thick (the term with $\tau_{R}$) and optically-thin (the term with $\tau_{P}$)
regimes. The Rosseland optical depth is given by
$\tau_{R}=\kappa_{\mathrm{disc}}(\rho_{\mathrm{mid}},T_{\mathrm{mid}})\Sigma$
where $\rho_{\mathrm{mid}}=\Sigma/(\sqrt{2\pi}H)$ is the central density,
$H=c_{s}/\Omega$ the disc’s vertical scale height,
$c_{s}=\sqrt{k_{\mathrm{B}}T_{\mathrm{mid}}/(\mu m_{\mathrm{H}})}$ the
isothermal sound speed, $\mu=2.24$ the mean molecular weight of the gas, and
$m_{\mathrm{H}}$ the mass of an hydrogen atom. The opacity
$\kappa_{\mathrm{disc}}$ is given by the maximum of the opacities computed
according to Bell & Lin (1994) (which accounts for micrometer size with a
fixed interstellar dust-to-gas ratio of $1\text{\,}\mathrm{\char 37\relax}$,
independently of the dust-to-gas 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 $\tau_{P}=2.4\tau_{R}$.

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

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

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

$\dot{E}=\Sigma\nu\left(r\frac{\partial\Omega}{\partial
r}\right)^{2}=\frac{9}{4}\Sigma\nu\Omega^{2}$ (4)

with $\Omega$ 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, i.e.
$\Omega=\sqrt{\mathcal{G}M_{\star}/r^{3}}$, $\mathcal{G}$ being the
gravitational constant. The self-gravity of the disc has been neglected.

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

$T_{\mathrm{s}}^{4}=T_{\star}^{4}\left[\frac{2}{3\pi}\left(\frac{R_{\star}}{r}\right)^{3}+\frac{1}{2}\left(\frac{R_{\star}}{r}\right)^{2}\frac{H}{r}\left(\frac{\partial\ln{H}}{\partial\ln{r}}-1\right)\right]+T_{\mathrm{irr}}^{4}+T_{\mathrm{cd}}^{4}$
(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
$\partial\ln{H}/\partial\ln{r}=9/7$ (Chiang & Goldreich 1997).

The $T_{\mathrm{irr}}$ term accounts for the direct irradiation through the
disc midplane. It is computed as

$T_{\mathrm{irr}}^{4}=\frac{L_{\star}}{16\pi
r^{2}\sigma_{\mathrm{SB}}}e^{-\tau_{\mathrm{mid}}},$ (6)

which is the black-body equilibrium temperature accounting for the optical
depth through the midplane of the disc
$\tau_{\mathrm{mid}}=\int\rho_{\mathrm{mid}}\kappa(\rho_{\mathrm{mid}},T_{\mathrm{mid}})\mathrm{d}r$.
This contribution is usually important only at the very end of the disc life
time 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_{\mathrm{cd}}=$10\text{\,}\mathrm{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 $\dot{M}_{\mathrm{wind}}$. External
photoevaporation is likely the most important environment-related factor for
discs (Winter et al. 2020).

####  3.2.2 Disc photoevaporation

Figure 3: Time evolution of the surface density (top left), opacity $\kappa$ (top right), midplane temperature (bottom left), and radial flow rate (bottom right) of a protoplanetary disc. The lines represent each one snapshot the state, and are spaced by about $2\text{\times}{10}^{5}\text{\,}\mathrm{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. 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_{\star}$ | $1\text{\,}\mathrm{M}_{\odot}$
Reference surface density $\Sigma_{g,0}$ at $5.2\text{\,}\mathrm{au}$ | $145\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$
Initial gas disc mass $M_{\mathrm{g}}$ | $3.90\text{\times}{10}^{-2}\text{\,}\mathrm{M}_{\odot}$
Inner edge of the gas disc $r_{\mathrm{in}}$ |  $0.091\text{\,}\mathrm{au}$ ($10\text{\,}\mathrm{d}$)
Characteristic radius of the gas disc $r_{\mathrm{cut,g}}$ | $66.5\text{\,}\mathrm{au}$
Disc viscosity parameter $\alpha$ | $2\text{\times}{10}^{-3}$
External photoevaporation rate $\dot{M}_{\mathrm{wind}}$ | $6.42\text{\times}{10}^{-7}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{yr}$
Power law index of the gas disc $\beta_{\mathrm{g}}$ | 0.9
Dust-to-gas ratio | $3.4\text{\,}\mathrm{\char 37\relax}$
Planetesimal disc mass | $348\text{\,}\mathrm{M}_{\oplus}$
Power law index of the solids disc $\beta_{\mathrm{s}}$ | 1.5
Characteristic radius of the solids disc $r_{\mathrm{cut,s}}$ | $r_{\mathrm{cut,g}}/2$
Planetesimal radius | $300\text{\,}\mathrm{m}$
Planetesimal density (rocky) | $3.2\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$
Planetesimal density (icy) | $1\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$
Embryo mass $M_{\mathrm{emb,0}}$ | $1\text{\times}{10}^{-2}\text{\,}\mathrm{M}_{\oplus}$
Opacity reduction factor $f_{\rm opa}$ | $3\text{\times}{10}^{-3}$

Photoevaporation in the protoplanetary discs is the principal means of
controlling their life times. 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 birth
place 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\text{\,}\mathrm{e}\mathrm{V}$) creates a neutral layer of dissociated
hydrogen whose temperature is
$T_{\mathrm{I}}\approx${10}^{3}\text{\,}\mathrm{K}$$. The corresponding sound
speed is then

$c^{2}_{\mathrm{s,I}}=\frac{k_{\mathrm{B}}T_{\mathrm{I}}}{\mu_{\mathrm{I}}m_{\mathrm{H}}},$
(7)

where the mean molecular weight $\mu_{\mathrm{I}}=1.35$ for the dissociated
gas. It corresponds to the gravitational radius (where the sound speed equals
the escape velocity) of

$r_{\mathrm{g,I}}=\frac{\mathcal{G}M_{\star}}{c^{2}_{\mathrm{s,I}}}.$ (8)

We assume that mass is removed uniformly outside of
$\beta_{\mathrm{I}}r_{\mathrm{g,I}}$ with $\beta_{\mathrm{I}}=0.14$ (similar
to Alexander & Pascucci 2012), so that the rate is given by

$\dot{\Sigma}_{\mathrm{g,photo,ext}}(r)=\left\\{\begin{array}[]{ll}0&\text{for}\leavevmode\nobreak\
r<\beta_{\mathrm{I}}r_{\mathrm{g,I}}\\\
\frac{\dot{M}_{\mathrm{wind}}}{\pi(r_{\mathrm{max}}^{2}-\beta^{2}_{\mathrm{I}}r^{2}_{\mathrm{g,I}})}&\text{otherwise}\end{array}\right.$
(9)

with $\dot{M}_{\mathrm{wind}}$ a parameter that provides the total mass loss
rate if the disc would extend to
$r_{\mathrm{max}}=$1000\text{\,}\mathrm{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_{\mathrm{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\text{\,}\mathrm{e}\mathrm{V}$$) creates a
layer of ionised hydrogen whose temperature is
$T_{\mathrm{II}}\approx${10}^{4}\text{\,}\mathrm{K}$$ and with a mean
molecular weight $\mu_{\mathrm{II}}=0.68$. The sound speed and gravitational
radius are computed in analogy with Eqs. (7) and (8). The scaling radius
$r_{14}=\beta_{\mathrm{II}}r_{\mathrm{g,II}}/${10}^{14}\text{\,}\mathrm{cm}$$
follows Clarke et al. (2001) while we select again $\beta_{\mathrm{II}}=0.14$
following Alexander & Pascucci (2012). With this, we can estimate the base
density with

$n_{0}(r_{14})=k_{\mathrm{Hol}}\Phi_{41}^{1/2}r_{14}^{-3/2},$ (10)

where we set $k_{\mathrm{Hol}}=$5.7\text{\times}{10}^{4}$$ following the
hydrodynamical simulations of Hollenbach et al. (1994) and
$\Phi_{41}=0.1\sqrt{M_{\star}/M_{\odot}}$, which is the ionising photon
luminosity in the units of ${10}^{41}\text{\,}{\mathrm{s}}^{-1}$. The
distance-dependent base density can then be calculated as

$n_{0}(r)=n_{0}(r_{14})\left(\frac{r}{r_{\mathrm{g,II}}}\right)^{-\frac{5}{2}}.$
(11)

We further follow Clarke et al. (2001) to get
$\dot{\Sigma}_{\mathrm{g,photo,int}}=2c_{\mathrm{s,II}}n_{0}m_{\mathrm{H}}$
outside of $\beta_{\mathrm{II}}r_{\mathrm{g,II}}$.

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

$\dot{\Sigma}_{\mathrm{g,photo}}=\dot{\Sigma}_{\mathrm{g,photo,ext}}+\dot{\Sigma}_{\mathrm{g,photo,int}}.$
(12)

####  3.2.3 Initial gas surface density profile and example

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

$\Sigma_{\mathrm{g}}(t=0)=\Sigma_{g,0}\left(\frac{r}{r_{\mathrm{0}}}\right)^{-\beta_{\mathrm{g}}}\exp{\left(-\left(\frac{r}{r_{\mathrm{cut,g}}}\right)^{2-\beta_{\mathrm{g}}}\right)}\left(1-\sqrt{\frac{r_{\mathrm{in}}}{r}}\right)$
(13)

where $r_{\mathrm{0}}=$5.2\text{\,}\mathrm{au}$$ is the reference distance,
$\beta_{\mathrm{g}}=0.9$ the power-law index (Andrews et al. 2010),
$r_{\mathrm{cut,g}}$ the characteristic radius for the exponential decay and
$r_{\mathrm{in}}$ the inner edge of the disc.

The conversion between the total mass and the normalisation surface density
$\Sigma_{g,0}$ at $r_{\mathrm{0}}$ is obtained with

$M_{\mathrm{g}}=\frac{2\pi\Sigma_{g,0}}{2-\beta_{\mathrm{g}}}\left(r_{\mathrm{0}}\right)^{\beta_{\mathrm{g}}}\left(r_{\mathrm{cut,g}}\right)^{2-\beta_{\mathrm{g}}}.$
(14)

It should be noted that this formula neglects the lack of gas within
$r_{\mathrm{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.
$\dot{\Sigma}_{\mathrm{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 life time
of that disc is nearly
$5.3\text{\times}{10}^{6}\text{\,}\mathrm{y}\mathrm{r}$.

The temporal evolution shows overall a decrease in the surface density. A hole
forms inside roughly $2\text{\,}\mathrm{au}$ by about
$4.7\text{\,}\mathrm{Myr}$. The change in the temperature profile initially
between 1.5 and $3\text{\,}\mathrm{au}$ and that moves inward 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\text{\,}\mathrm{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\text{\,}\mathrm{au}$ due to this contribution; otherwise the midplane
temperature remains below the equilibrium temperature, apart from the inner
region (¡$<3\text{\,}\mathrm{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\text{\,}\mathrm{\char 37\relax}$ coupled with larger temperature by
$20$–$40\text{\,}\mathrm{\char 37\relax}$. 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 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\text{\,}\mathrm{L}_{\odot}$ as predicted by the Baraffe et al. (2015)
tracks whereas Bitsch et al. (2015a) begins with
$1.5\text{\,}\mathrm{L}_{\odot}$ 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).

###  3.3 Planetesimal disc

Planetesimals are represented by a fluid-like description, i.e., they are
modelled not as individual particles but on a grid as a surface density
($\Sigma_{\mathrm{s}}$) with eccentricity ($e_{\mathrm{plan}}$) and
inclination ($i_{\mathrm{plan}}$) as dynamical state.

####  3.3.1 Dynamical state

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

$\displaystyle\dot{e}_{\mathrm{plan}}^{2}$ $\displaystyle=$
$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}+\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,M}}+\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,plan}}$
(15) $\displaystyle\dot{i}_{\mathrm{plan}}^{2}$ $\displaystyle=$
$\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}+\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,M}}+\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,plan}}.$
(16)

The contributions from the aerodynamical drag, stirring by the proto-planets
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 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
$\mathrm{Re}_{\mathrm{mol}}=v_{\mathrm{rel}}R_{\mathrm{plan}}/\nu_{\mathrm{mol}}$,
where $\nu_{\mathrm{mol}}=\lambda c_{\mathrm{s}}/3$ is the molecular
viscosity, $\lambda=(n_{\mathrm{H_{2}}}\sigma_{\mathrm{H_{2}}})^{-1}$ the gas
molecules’ mean free path, $n_{\mathrm{H_{2}}}$ the number density assuming
all of the gaseous molecules having hydrogen mass, $\sigma_{\mathrm{H_{2}}}$
their collisional cross-section, $R_{\mathrm{plan}}$ the planetesimals’
radius,

$v_{\mathrm{rel}}=v_{\mathrm{K}}\sqrt{\eta^{2}+5/8e_{\mathrm{plan}}^{2}+1/2i_{\mathrm{plan}}^{2}}$
(17)

their relative velocity,

$\eta=-\frac{1}{2\Omega r\rho_{\mathrm{mid}}}\frac{\partial p}{\partial r}$
(18)

the deviation between the gas and Keplerian velocities due the support of the
gas by the radial pressure gradient, $\rho_{\mathrm{mid}}$ the midplane gas
density, and $v_{\mathrm{K}}=\Omega r$ the Keplerian velocity. When
$R_{\mathrm{plan}}<\lambda$, the gas drag is assumed to be in the Epstein
regime. Otherwise, if $\mathrm{Re}_{\mathrm{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),

$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$
$\displaystyle-\frac{e_{\mathrm{plan}}^{2}}{\tau_{\mathrm{drag}}}\sqrt{\eta^{2}+\frac{5}{8}e_{\mathrm{plan}}^{2}+\frac{1}{2}i_{\mathrm{plan}}^{2}}$
(19) $\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$
$\displaystyle-\frac{i_{\mathrm{plan}}^{2}}{2\tau_{\mathrm{drag}}}\sqrt{\eta^{2}+\frac{5}{8}e_{\mathrm{plan}}^{2}+\frac{1}{2}i_{\mathrm{plan}}^{2}},$
(20)

where

$\tau_{\mathrm{drag}}=\frac{8\rho_{\mathrm{plan}}R_{\mathrm{plan}}}{3C_{\mathrm{D}}\rho_{\mathrm{mid}}v_{\mathrm{K}}}$
(21)

is the gas drag time scale and $C_{\mathrm{D}}=1$.

In the Stokes regimes the drag expressions are

$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$
$\displaystyle-\frac{3}{2}\frac{\lambda\rho_{\mathrm{mid}}e_{\mathrm{plan}}^{2}}{\rho_{\mathrm{plan}}R_{\mathrm{plan}}^{2}}$
(22) $\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$
$\displaystyle-\frac{3}{4}\frac{\lambda\rho_{\mathrm{mid}}i_{\mathrm{plan}}^{2}}{\rho_{\mathrm{plan}}R_{\mathrm{plan}}^{2}},$
(23)

while in the Epstein regime they read as

$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$ $\displaystyle-
e_{\mathrm{plan}}^{2}\frac{c_{\mathrm{s}}\rho_{\mathrm{mid}}}{\rho_{\mathrm{plan}}R_{\mathrm{plan}}}$
(24) $\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$
$\displaystyle=$
$\displaystyle-\frac{i_{\mathrm{plan}}^{2}}{2}\frac{c_{\mathrm{s}}\rho_{\mathrm{mid}}}{\rho_{\mathrm{plan}}R_{\mathrm{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 in Fortier et al. (2013), the stirring by the protoplanets follows the
approach of Guilera et al. (2010), where the stirring of Ohtsuki et al. (2002)
is modulated with the separation from the protoplanets. The contribution reads
as

$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,M}}$
$\displaystyle=$
$\displaystyle\sum_{\jmath=1}^{n}f_{\Delta,\jmath}\left(\frac{\Omega
M_{\mathrm{planet,\jmath}}}{6\pi bM_{\star}}\right)P_{\mathrm{VS}}$ (26)
$\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,M}}$
$\displaystyle=$
$\displaystyle\sum_{\jmath=1}^{n}f_{\Delta,\jmath}\left(\frac{\Omega
M_{\mathrm{planet,\jmath}}}{6\pi bM_{\star}}\right)Q_{\mathrm{VS}}$ (27)

where the sum is over all the protoplanets present in the system,

$f_{\Delta,\jmath}^{-1}=1+\left(\frac{|r-a_{\mathrm{planet,\jmath}}|}{bR_{\mathrm{H}}}\right)^{5}$
(28)

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

$R_{\mathrm{H}}=a_{\mathrm{planet,\jmath}}\sqrt[3]{\frac{M_{\mathrm{planet,\jmath}}}{3M_{\star}}}$
(29)

the planet’s Hills radius, and $b=5$ is the half-width of the feeding zone
(see Sect. 4.3.3). The terms $P_{\mathrm{VS}}$ and $Q_{\mathrm{VS}}$ are given
by (Fortier et al. 2013),

$\displaystyle P_{\mathrm{VS}}$ $\displaystyle=$
$\displaystyle\left[\frac{73\tilde{e}_{\mathrm{plan}}^{2}}{10\Lambda^{2}}\right]\ln{\left(1+10\Lambda^{2}/\tilde{e}_{\mathrm{plan}}^{2}\right)}$
(30)
$\displaystyle+\left[\frac{72I_{\mathrm{PVS}}(\beta)}{\pi\tilde{e}_{\mathrm{plan}}\tilde{\imath}_{\mathrm{plan}}}\right]\ln{\left(1+\Lambda^{2}\right)}$
$\displaystyle Q_{\mathrm{VS}}$ $\displaystyle=$
$\displaystyle\left[\frac{4\tilde{\imath}_{\mathrm{plan}}^{2}+0.2\tilde{\imath}_{\mathrm{plan}}\tilde{e}_{\mathrm{plan}}^{3}}{10\Lambda^{2}\tilde{e}_{\mathrm{plan}}}\right]\ln{\left(1+10\Lambda^{2}\tilde{e}_{\mathrm{plan}}^{2}\right)}$
(31)
$\displaystyle+\left[\frac{72I_{\mathrm{QVS}}(\beta)}{\pi\tilde{e}_{\mathrm{plan}}\tilde{\imath}_{\mathrm{plan}}}\right]\ln{\left(1+\Lambda^{2}\right)}$

Here, $\tilde{e}_{\mathrm{plan}}=re_{\mathrm{plan}}/R_{\mathrm{H}}$ and
$\tilde{\imath}_{\mathrm{plan}}=ri_{\mathrm{plan}}/R_{\mathrm{H}}$ are
respectively the reduced planetesimals’ eccentricity and inclination,
$\Lambda=\tilde{\imath}_{\mathrm{plan}}(\tilde{e}_{\mathrm{plan}}^{2}+\tilde{\imath}_{\mathrm{plan}}^{2})/12$,
$\beta=i_{\mathrm{plan}}/e_{\mathrm{plan}}$, while for $I_{\mathrm{PVS}}$ and
$I_{\mathrm{QVS}}$ we use the approximations obtained by Chambers (2006):

$\displaystyle I_{\mathrm{PVS}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\frac{\beta-0.36251}{0.061547+0.16112\beta+0.054473\beta^{2}},$
(33) $\displaystyle I_{\mathrm{QVS}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\frac{0.71946-\beta}{0.21239+0.49764\beta+0.14369\beta^{2}}.$
(34)

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

$\displaystyle\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,plan}}$
$\displaystyle=$
$\displaystyle\frac{1}{6}\sqrt{\frac{\mathcal{G}r}{M_{\star}}}\Sigma_{\mathrm{s}}h_{\mathrm{plan}}P_{\mathrm{VS}}$
(35)
$\displaystyle\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{VS,plan}}$
$\displaystyle=$
$\displaystyle\frac{1}{6}\sqrt{\frac{\mathcal{G}r}{M_{\star}}}\Sigma_{\mathrm{s}}h_{\mathrm{plan}}Q_{\mathrm{VS}}$
(36)

with

$h_{\mathrm{plan}}=\sqrt[3]{\frac{2M_{\mathrm{plan}}}{3M_{\star}}},$ (37)

and $M_{\mathrm{plan}}=4/3\pi R_{\mathrm{plan}}^{3}\rho_{\mathrm{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_{\mathrm{plan}}=2.31\frac{M_{\mathrm{plan}}^{4/15}r^{1/5}\rho_{\mathrm{plan}}^{2/15}\Sigma_{\mathrm{g}}^{1/5}}{C_{D}^{1/5}\rho^{1/5}M_{\star}^{2/5}}$
(38)

and

$i_{\mathrm{plan}}=\frac{1}{2}e_{\mathrm{plan}}.$ (39)

We also compared our prescription for the dynamical state with gamma-stirring
from, e.g., Ida et al. (2008) and Okuzumi & Ormel (2013). Although this is not
straightforward due to the differences in the sources, we find that,
generally, the eccentricities resulting from $\gamma$-stirring are larger than
the self-stirring from the planetesimals, but lower than the stirring by the
forming protoplanets. Thus, taking into account 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 (Alibert et al. 2013), we assume a
constant planetesimals radius of $300\text{\,}\mathrm{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 ($\sim
1\text{\,}\mathrm{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 $\sim 100\text{\,}\mathrm{km}$ planetesimals.
The first two points can, however, also be explained with other effects (Zheng
et al. 2017; Wei et al. 2018, although the former work makes no determination
about the initial size frequency distribution of planetesimals).

In the more specific context of the simulations presented here, this choice
was made for the following reasons: 1) small planetesimals undergo sufficient
eccentricity and inclination damping by the disc gas to sustain a planetesimal
accretion rate in the oligarchic growth regime that is high enough to build
giant planet cores during typical disc lifetimes (Fortier et al. 2013). 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 life times of gas discs (Burn et al. 2019) and 3) this size was shown
to be able to reproduce several of the known exoplanet properties across a
wide range of masses (Fortier et al. 2013). In any case, the constant
planetesimal size is an important limitation of the model. Including in the
Bern Model an explicit model for the evolution of the solid building blocks
across the entire size range (dust-pebble-planetesimals) is thus subject of
ongoing research. A first important step was recently made in Voelkel et al.
(2020) where we have coupled the dust/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,

$\Sigma_{s}(t=0)=\Sigma_{\mathrm{s,0}}f_{\mathrm{s}}(r)\left(\frac{r}{r_{0}}\right)^{-\beta_{\mathrm{s}}}\exp{\left(-\left(\frac{r}{r_{\mathrm{cut,s}}}\right)^{2}\right)}$
(40)

with the power-law exponent is set to $\beta_{\mathrm{s}}=1.5$, as in the
MMSN, and $r_{\mathrm{cut,s}}=r_{\mathrm{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 $\Sigma_{\mathrm{s,0}}$ is adjusted so
that the bulk solids-to-gas ratio remains to the prescribed value (e.g.,
$1\text{\,}\mathrm{\char 37\relax}$).

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_{\mathrm{plan}}^{2}>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 Raymond & Izidoro (2017) for instance. Finally, the planetesimals
disc remains after the dispersal of the gas disc; the only difference is that
the damping terms for eccentricity
$\left.\dot{e}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$ and inclination
$\left.\dot{i}_{\mathrm{plan}}^{2}\right|_{\mathrm{drag}}$ vanish.

####  3.3.3 Compositional model

Figure 4: Time evolution of the RMS of planetesimals’ eccentricity (top left),
inclination (top right), and surface density (bottom left) of a circumstellar
disc that also contains 10 embryos. The lines represent temporal snapshots of
the three quantities, and are spaced by about
$2\text{\times}{10}^{5}\text{\,}\mathrm{y}\mathrm{r}$. 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 life time of the gas disc is shorter than the case
presented in Fig. 3 due to the accretion by the protoplanets.

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).
Volatiles 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_{\mathrm{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 taken into account for the resulting solid surface density.
Thus, the value of $f_{\mathrm{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
$\rho_{\mathrm{plan}}=$3.2\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$
while when volatiles are in the solid phase we take
$\rho_{\mathrm{plan}}=$1\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$. This
transition corresponds to the H2O-ice line in all discs, which induces the
largest surface density jump because H2O makes up
$\sim$$60\text{\,}\mathrm{\char 37\relax}$ of all ices in mass (Marboeuf et
al. 2014a).

####  3.3.4 Example

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

The different jumps in the initial surface density profile are due to the
crossing of the different ice lines; the most consequential one at at about
$3\text{\,}\mathrm{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, e.g., Levison et al. (2010). In total, the
planets accreted $61\text{\,}\mathrm{M}_{\oplus}$ of planetesimals
($47\text{\,}\mathrm{M}_{\oplus}$ of which by the giant planet) while
$89\text{\,}\mathrm{M}_{\oplus}$ 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\text{\,}\mathrm{\char 37\relax}$ of
planetesimals were ejected and $35\text{\,}\mathrm{\char 37\relax}$ accreted.

The stirring by the proto-planets 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\text{\,}\mathrm{au}$),
which reaches a mass of about
$5.4\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$ at the end of the formation
stage. That planet has a core mass of $47\text{\,}\mathrm{M}_{\oplus}$, which
corresponds (for a pure H/He envelope) to a metallicity slightly lower than
that of the star ($2.8\text{\,}\mathrm{\char 37\relax}$ versus
$3.0\text{\,}\mathrm{\char 37\relax}$). 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
$\beta=i_{\mathrm{plan}}/e_{\mathrm{plan}}\approx 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 significant
reduce the inward gas flow such that the region inside its orbit becomes gas-
poor. This greatly reduces the damping of the planetesimals dynamical state to
a such point that their eccentricity becomes close to unity.

##  4 Planet properties

###  4.1 Envelope structure

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

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

The downside is that solving the internal structure for bodies ranging in mass
from ${10}^{-2}\text{\,}\mathrm{M}_{\oplus}$ 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.

####  4.1.1 Attached phase

In the initial phase, known as the attached phase, the envelope is in
equilibrium with the gas disc and the gas density smoothly transitions from
the value in the protoplanetary envelope to the one in the background nebula.
The planets do not yet have a well defined outer radius. During this phase,
the gas accretion rate is governed by the ability to radiate the gravitational
energy liberated by the accretion of solids and gas, and the envelope’s
contraction. For the forming giant planets, this phase generally lasts until
the planets reach a total mass in the range of $30$ to
$100\text{\,}\mathrm{M}_{\oplus}$ 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 (i.e., 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),

$\displaystyle\partial M/\partial R$ $\displaystyle=$ $\displaystyle 4\pi
R^{2}\rho$ (41) $\displaystyle\partial P/\partial R$ $\displaystyle=$
$\displaystyle-\mathcal{G}M\rho/R^{2}$ (42) $\displaystyle\partial T/\partial
R$ $\displaystyle=$ $\displaystyle\partial P/\partial
R\min{(\nabla_{\mathrm{ad}},\nabla_{\mathrm{rad}})},$ (43)

with $M$ the mass enclosed in the radius $R$, $P$ the pressure, $T$ the
temperature, $\rho=\rho(P,T)$ the density, computed using the SCvH equation of
state (Saumon et al. 1995), and $\nabla_{\mathrm{ad}}$ and
$\nabla_{\mathrm{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

$\nabla_{\mathrm{rad}}=\frac{3\kappa
L}{64\pi\sigma_{\mathrm{SB}}\mathcal{G}MT^{3}},$ (44)

with $L$ being the luminosity.

The opacity in the envelope $\kappa$ 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_{\mathrm{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. (2009),

$\frac{1}{R_{\mathrm{tot}}}=\frac{1}{k_{1}R_{\mathrm{acc}}}+\frac{1}{k_{2}R_{\mathrm{H}}},$
(45)

where

$R_{\mathrm{acc}}=\frac{\mathcal{G}M_{\mathrm{tot}}}{c_{s}^{2}}$ (46)

is the Bondi radius, $R_{\mathrm{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

$\displaystyle P(R_{\mathrm{tot}})$ $\displaystyle=$ $\displaystyle
P_{\mathrm{neb}}(a_{\mathrm{planet}})\ \ \mathrm{and}$ (47) $\displaystyle
T^{4}(R_{\mathrm{tot}})$ $\displaystyle=$ $\displaystyle
T_{\mathrm{neb}}^{4}(a_{\mathrm{planet}})+\frac{3\tau_{\mathrm{out}}L(R_{\mathrm{tot}})}{8\pi\sigma_{\mathrm{SB}}R_{\mathrm{tot}}^{2}},$
(48)

and

$\tau_{\mathrm{out}}=\max{\left(\kappa(\rho_{\mathrm{neb}},T_{\mathrm{neb}})\rho_{\mathrm{neb}}R_{\mathrm{tot}},\frac{2}{3}\right)}$
(49)

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_{\mathrm{tot}})$ is described in Section 4.2. In the case of the mass,
what is known is the core mass, i.e. $M(R_{\mathrm{core}})=M_{\rm core}$,
while $M(R_{\mathrm{tot}})=M_{\mathrm{tot}}$ is the quantity that is being
searched for. We thus use an iterative method by guessing $M_{\mathrm{tot}}$,
which is then used to integrate the internal structure equations until the
boundary condition at the inner boundary is fulfilled, i.e.
$M(R_{\mathrm{core}})=M_{\rm core}$. Once $M_{\mathrm{tot}}$ is found, the
envelope mass can be retrieved by $M_{\mathrm{env}}=M_{\mathrm{tot}}-M_{\rm
core}$, and the gas accretion rate by taking the difference of the envelope
mass between two successive steps of the envelope structure calculation
$\dot{M}_{\mathrm{env}}=(M_{\mathrm{env}}(t)-M_{\mathrm{env}}(t-\Delta
t))/\Delta t$.

####  4.1.2 Maximum gas accretion rate

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

$\tau_{\rm
KH}=\frac{\mathcal{G}M_{\mathrm{tot}}^{2}}{R_{\mathrm{tot}}L_{\mathrm{tot}}}.$
(50)

However, as the planet’s core reaches a mass of about
$10\text{\,}\mathrm{M}_{\oplus}$, the value of $\tau_{\rm 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 $\dot{M}_{\mathrm{env,max}}$, which
is used to limit the value of $\dot{M}_{\mathrm{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_{\mathrm{gc}}$ around the planet. For simplicity, we compute
$R_{\mathrm{gc}}$ according to Eq. (45). Depending on the value of
$R_{\mathrm{gc}}$ with respect to $H$, the local disc’s scale height, two
different regimes occur. In the case where $R_{\mathrm{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 $\sigma_{\mathrm{cross}}=\pi
R_{\mathrm{gc}}^{2}$ is given by

$\dot{M}_{\mathrm{env,max,3D}}=\rho\sigma_{\mathrm{cross}}v_{\mathrm{rel}}$
(51)

with $\rho\approx\Sigma/H$ the approximate density of the gas and
$v_{\mathrm{rel}}=\max{\left(\Omega R_{\mathrm{tot}},c_{\mathrm{s}}\right)}$
the relative velocity between the gas and the planet.

On the other hand, in the case $R_{\mathrm{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 have

$\dot{M}_{\mathrm{env,max,2D}}=2R_{\mathrm{gc}}\Sigma v_{\mathrm{rel}}.$ (52)

To distinguish between the two regimes, we use the lower rate of the two, i.e.

$\dot{M}_{\mathrm{env,max}}=\min{\left(\dot{M}_{\mathrm{env,max,2D}},\dot{M}_{\mathrm{env,max,3D}}\right)}$
(53)

Finally, to ensure that no more gas than available in the feeding zone
$M_{\mathrm{feed}}$ is accreted during one time step, we further constrain
$\dot{M}_{\mathrm{env,max}}<M_{\mathrm{feed}}/\Delta 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 $\sim
5-$10\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$$.

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

The radial extent of the feeding zone is set by

$\displaystyle r_{\mathrm{feed,inf}}$ $\displaystyle=$ $\displaystyle
a_{\mathrm{planet}}\left(1-e\right)\left(1-f_{\mathrm{feed}}\sqrt{\frac{M_{\mathrm{tot}}}{3M_{\star}}}\right)$
(54) $\displaystyle r_{\mathrm{feed,sup}}$ $\displaystyle=$ $\displaystyle
a_{\mathrm{planet}}\left(1+e\right)\left(1+f_{\mathrm{feed}}\sqrt{\frac{M_{\mathrm{tot}}}{3M_{\star}}}\right)$
(55)

with $f_{\mathrm{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 provide the location over which the disc properties ($\Sigma$, $H$,
etc.) are averaged for the calculation of the maximum rate and the removal of
the accreted gas, with

$\dot{\Sigma}_{\mathrm{g,planet}}=\frac{\dot{M}_{\mathrm{env}}}{\pi\left(r_{\mathrm{feed,sup}}^{2}-r_{\mathrm{feed,inf}}^{2}\right)}.$
(56)

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 $\dot{\Sigma}_{\mathrm{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). This could lead to an overestimation of gas accretion at lower to
intermediate giant planet masses. This could potentially explain why our
current model of disc-limited gas accretion seems to too strongly reduce the
stellar gas accretion rate (Manara et al. 2019; Bergez-Casalou et al. 2020).
Gap formation would reduce this effect, but could potentially lead to another
issue: observationally, the giant planet mass function seems to extend
smoothly to about $30\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$ (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\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$ 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.

####  4.1.3 Detached phase

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

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

$P(R_{\mathrm{tot}})=P_{\mathrm{neb}}(a_{\mathrm{planet}})+P_{\mathrm{edd}}+P_{\mathrm{ram}}+P_{\mathrm{rad}}$
(57)

with $P_{\mathrm{neb}}(a_{\mathrm{planet}})$ being the pressure at the
midplane of the gas disc, $P_{\mathrm{edd}}=(2g)/(3\kappa)$ the Eddington
expression for the photospheric pressure due to the material residing above
the $\tau=2/3$ surface,
$P_{\mathrm{rad}}=(2\sigma_{\mathrm{SB}}T^{4}(R_{\mathrm{tot}}))/(3c)$ the
radiation pressure, $c$ being the speed of light in vaccum, and

$\displaystyle P_{\mathrm{ram}}=\frac{\dot{M}_{\mathrm{env}}}{4\pi
R_{\mathrm{tot}}^{2}}v_{\mathrm{ff}}$ ; $\displaystyle
v_{\mathrm{ff}}^{2}=2\mathcal{G}M_{\mathrm{tot}}\left(\frac{1}{R_{\mathrm{tot}}}-\frac{1}{R_{\mathrm{H}}}\right)$
(58)

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

####  4.1.4 Evolutionary phase

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

$\displaystyle P(R_{\mathrm{tot}})$ $\displaystyle=$ $\displaystyle
P_{\mathrm{edd}}+P_{\mathrm{rad}}$ (59) $\displaystyle
T^{4}(R_{\mathrm{tot}})$ $\displaystyle=$ $\displaystyle
T^{4}_{\mathrm{int}}+(1-A)T^{4}_{\mathrm{eq}}$ (60)

where
$T_{\mathrm{int}}^{4}=L_{\mathrm{tot}}/(4\pi\sigma_{\mathrm{SB}}R_{\mathrm{tot}}^{2})$
is the intrinsic temperature,
$T_{\mathrm{eq}}=T_{\star}*\sqrt{R_{\star}/(2*a_{\mathrm{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).

####  4.1.5 Example

Figure 6: Snapshots of the internal structure of the second outermost planet
of Fig. 4. The structures are split according to the phases, with attached
(top left), transition (the initial stage of the detached phase; top right),
detached (bottom left) and evolutionary (bottom right). The red line shows the
first profile of the detached phase and is shown of both panels. The green and
blue profiles lie at the transition between two stages and are shown of two
panels each. In each profile, thin line show the part where energy transport
is radiative and thick lines for convective. Figure 7: Time evolution of the
planet’s radius and luminosity of the second outermost planet of Fig. 4; the
same as in Fig. 6. The insert on the left panel shows the contraction at the
transition between attached and detached phase. The exact time where the model
switches between the two phases is shown with the vertical dashed line.

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\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$. 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.

###  4.2 Luminosity

####  4.2.1 Accretion and contraction

The luminosity calculation suffers from the same problem as the total mass in
the attached phase, or the outer radius in the detached phase; i.e. 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

$E_{\mathrm{tot}}=-\mathcal{G}\int_{0}^{M_{\mathrm{tot}}}\frac{M}{R}\mathrm{d}M+\int_{M_{\rm
core}}^{M_{\mathrm{tot}}}u\mathrm{d}M=-\xi\frac{\mathcal{G}M_{\mathrm{tot}}^{2}}{2R_{\mathrm{tot}}}$
(61)

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/5\mathcal{G}M_{\rm
core}^{2}/R_{\mathrm{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 $\xi$ 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 form the accretion, cooling, and contraction $L_{\mathrm{int}}$ can
then be obtained as

$L_{\mathrm{int}}=\dot{E}_{\mathrm{tot}}=\frac{\xi\mathcal{G}M_{\mathrm{tot}}}{R_{\mathrm{tot}}}\dot{M}_{\mathrm{tot}}-\frac{\xi\mathcal{G}M_{\mathrm{tot}}^{2}}{2R_{\mathrm{tot}}^{2}}\dot{R}_{\mathrm{tot}}+\frac{\mathcal{G}M_{\mathrm{tot}}^{2}}{2R_{\mathrm{tot}}}\dot{\xi}_{\mathrm{tot}}.$
(62)

with $\dot{M}_{\mathrm{tot}}=\dot{M}_{\rm core}+\dot{M}_{\mathrm{env}}$ being
the total accretion rate of the planet (solids and gas). The value
$\dot{M}_{\mathrm{tot}}$ in the attached phase and of $\dot{R}_{\mathrm{tot}}$
in the detached phase are determined from the guess for the mass or radius
during the iterations. The same is not true for $\dot{\xi}_{\mathrm{tot}}$. To
circumvent this problem, we estimate the luminosity with

$L_{\mathrm{int}}\approx
C\left(\frac{\xi\mathcal{G}M_{\mathrm{tot}}}{R_{\mathrm{tot}}}\dot{M}_{\mathrm{tot}}-\frac{\xi\mathcal{G}M_{\mathrm{tot}}^{2}}{2R_{\mathrm{tot}}^{2}}\dot{R}_{\mathrm{tot}}\right).$
(63)

The correction factor $C$ corrects for neglecting the
$\dot{\xi}_{\mathrm{tot}}$ term. The value of $C$ can be calculated a
posteriori by determining the actual total energy of the new planet, with
$C=-(E_{\mathrm{tot}}(t)-E_{\mathrm{tot}}(t-\Delta
t))/(L_{\mathrm{int}}\cdot\Delta t)$. The value of $C$ is then used for the
next time step.

Marleau et al. (2017, 2019b) conducted 1D radiation-hydrodynamic simulations
of the planetary gas accretion shock, a feature that is seen in various 3D
radiation-hydrodynamic simulations of accreting protoplanets of sufficiently
high mass (e.g., Szulágyi & Mordasini 2017; Schulik et al. 2020). High
postshock entropies were found, suggesting that warm or hot gas accretion is
more plausible than cold accretion (see also Berardo et al. 2017; Berardo &
Cumming 2017). We therefore assume in our model that gas accretion in the
detached phase is hot, which means that we do not subtract the accretion shock
luminosity from $L_{\mathrm{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_{\mathrm{radio}}$ includes contributions from the three most important
long-lived radionucleides 40K, 238U and 232Th (Wasserburg et al. 1964). To
compute the luminosity contributions, we follow the procedure of Mordasini et
al. (2012b): we assume the mantle of the protoplanets has a chrondritic
composition and the energy production rate are retrieved from meteoritic
values of William (2007). The initial radiogenic contribution is
$Q_{0}\approx$5\text{\times}{10}^{-7}\text{\,}\mathrm{erg}\text{\,}{\mathrm{g}}^{-1}\text{\,}{\mathrm{s}}^{-1}$$
of mantle material (all elements besides iron).

####  4.2.2 Bloating of close-in planets

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

$L_{\mathrm{bloat}}=\epsilon F_{\star}e^{-\tau_{\mathrm{mid}}}\pi
R_{\mathrm{tot}}^{2}$ (64)

with

$\epsilon=2.37\exp{\left(-\frac{\left(\log{(F_{\star}/${10}^{9}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$)}-0.14\right)^{2}}{2\cdot
0.37^{2}}\right)},$ (65)

$F_{\star}=L_{\star}/(4\pi a_{\mathrm{planet}}^{2})$ the total stellar flux at
the planet’s location, and $\tau_{\mathrm{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_{\star}$ (in the evolutionary
phase) or the stellar flux multiplied by the optical depth
$F_{\star}\exp{\left(-\tau_{\mathrm{mid}}\right)}$ (before the dispersal of
the gas disc) is greater than
$2\text{\times}{10}^{8}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$
(Demory & Seager 2011).

####  4.2.3 Deuterium-burning

For the calculation of the luminosity due to deuterium fusion, we follow the
procedure of Mollière & Mordasini (2012). In this framework, the energy
generation rate (per unit mass and time) is given by Kippenhahn & Weigert
(1994), with the assumption that nuclei are fully ionised and non-degenerate.
The energy released in each process is computed according to Fowler et al.
(1967). The specific deuterium burning luminosity of a planet depends on the
conditions in the planet’s gaseous envelope, most notably the density,
temperature, and the remaining deuterium nuclei. This implies that there is no
universal mass at which deuterium burning starts, but as already found in
Mollière & Mordasini (2012) (see also Bodenheimer et al. 2013), the mass where
burning becomes important clusters around about
$13\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$. 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\text{\times}{10}^{-5}$, which is the primordial
value.

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

####  4.2.4 Total luminosity

The final luminosity is then given by

$L(R_{\mathrm{tot}})=L_{\mathrm{tot}}=L_{\mathrm{int}}+L_{\mathrm{radio}}+L_{\mathrm{bloat}}+L_{\mathrm{D-burn}}.$
(66)

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

###  4.3 Accretion of solids

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

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

$\dot{M}_{\rm
core}=\Omega\bar{\Sigma}_{\mathrm{s}}R_{\mathrm{H}}^{2}p_{\mathrm{coll}},$
(67)

with $\bar{\Sigma}_{\mathrm{s}}$ the mean surface density of planetesimals in
the planet’s feeding zone and $p_{\mathrm{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.

####  4.3.1 Capture probability

We distinguish three different accretion regimes depending on the random
velocities: low-, mid- and high-velocity. The distinction is based on the
reduced planetesimals’ eccentricity
$\tilde{e}_{\mathrm{plan}}=re_{\mathrm{plan}}/R_{\mathrm{H}}$ and inclination
$\tilde{\imath}_{\mathrm{plan}}=ri_{\mathrm{plan}}/R_{\mathrm{H}}$ ($r$ is the
heliocentric distance): the high-velocity regime for
$\tilde{e}_{\mathrm{plan}},\,\tilde{\imath}_{\mathrm{plan}}\gtrsim 2$, mid-
velocity for
$2\gtrsim\tilde{e}_{\mathrm{plan}},\,\tilde{\imath}_{\mathrm{plan}}\gtrsim
0.2$ and low-velocity for
$0.2\gtrsim\tilde{e}_{\mathrm{plan}},\,\tilde{\imath}_{\mathrm{plan}}$.
According to Inaba et al. (2001), each regime has a different expression for
the collision probability,

$\displaystyle p_{\mathrm{high}}$ $\displaystyle=$
$\displaystyle\frac{(R_{\mathrm{cap}}+R_{\mathrm{plan}})^{2}}{2\pi
R_{\mathrm{H}}}\left(I_{\mathrm{F}}(\beta)+\frac{6R_{\mathrm{H}}I_{\mathrm{G}}(\beta)}{(R_{\mathrm{cap}}+R_{\mathrm{plan}})\tilde{e}_{\mathrm{plan}}^{2}}\right)$
(68) $\displaystyle p_{\mathrm{mid}}$ $\displaystyle=$
$\displaystyle\frac{(R_{\mathrm{cap}}+R_{\mathrm{plan}})^{2}}{4\pi
R_{\mathrm{H}}\tilde{\imath}_{\mathrm{plan}}}\left(17.3+\frac{232R_{\mathrm{H}}}{R_{\mathrm{cap}}+R_{\mathrm{plan}}}\right)$
(69) $\displaystyle p_{\mathrm{low}}$ $\displaystyle=$ $\displaystyle
11.3\left(\frac{R_{\mathrm{cap}}+R_{\mathrm{plan}}}{R_{\mathrm{H}}}\right),$
(70)

where $R_{\mathrm{cap}}$ is the planetesimal capture radius of the planet,
$\beta=i_{\mathrm{plan}}/e_{\mathrm{plan}}$ and the $I_{\mathrm{F}}$ and
$I_{\mathrm{G}}$ functions can be approximated as, following Chambers (2006):

$\displaystyle I_{\mathrm{F}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\frac{1+0.95925\beta+0.77251\beta^{2}}{\beta\left(0.13142+0.12295\beta\right)}$
(71) $\displaystyle I_{\mathrm{G}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\frac{1+0.39960\beta}{\beta\left(0.0369+0.048333\beta+0.006874\beta^{2}\right)}.$
(72)

The final collision is then given by

$p_{\mathrm{coll}}=\min{\left(p_{\mathrm{med}},\left(p_{\mathrm{high}}^{-2}+p_{\mathrm{low}}^{-2}\right)^{-1/2}\right)}$
(73)

Figure 8: Evolution of the planet’s core, planetesimals capture, and total
radii as function of time (left panel) and core mass (right panel) for the
second outermost planet of Fig. 4. In the right panel, we also include two
time indicators at ${10}^{5}\text{\,}\mathrm{yr}$ and
${10}^{6}\text{\,}\mathrm{yr}$ with dashed vertical lines. The upper-right
portion of the plot for the total radius in the left panel is also shown on
the left panel of Fig. 7.

In the initial stage, the capture radius $R_{\mathrm{cap}}$ is the physical
radius of the core $R_{\mathrm{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

$R_{\mathrm{plan}}=\frac{3}{2}\frac{\rho(R_{\mathrm{cap}})R_{\mathrm{cap}}}{\rho_{\mathrm{plan}}}\left(\frac{v_{\mathrm{rel}}^{2}+2\mathcal{G}M(R_{\mathrm{cap}})/R_{\mathrm{cap}}}{v_{\mathrm{rel}}^{2}+2\mathcal{G}M(R_{\mathrm{cap}})/R_{\mathrm{H}}}\right).$
(74)

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}\text{\,}\mathrm{yr}$, before that, the capture
radius is equal to that of the core. At that moment, the core mass is
$9\text{\times}{10}^{-2}\text{\,}\mathrm{M}_{\oplus}$. By the time the core
reaches $1\text{\,}\mathrm{M}_{\oplus}$ at
$4.8\text{\times}{10}^{5}\text{\,}\mathrm{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}\text{\,}\mathrm{M}_{\oplus}$), and the
calculation of gaseous envelopes can not 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.

####  4.3.2 Ejection of planetesimals

Planets not only accrete material, they also induce gravitational
perturbations on the planetesimals that come close-by but are not accreted.
These planetesimals, if they receive a sufficient velocity kick from a close
approach by a planet, can be ejected from the system. To estimate this effect,
we follow a procedure similar to Ida & Lin (2004a). The planetesimals that
receive a velocity kick greater than the escape velocity from the primary,
$v_{\mathrm{esc}}=\sqrt{2\mathcal{G}M_{\star}/a_{\mathrm{planet}}}$, will
likely be ejected from the system. Thus we have that the fraction of accreted-
to-ejected planetesimals is (Ida & Lin 2004a)

$\frac{\dot{M}_{\rm
core}}{\dot{M}_{\mathrm{ejec}}}=\left(\frac{v_{\mathrm{esc}}}{v_{\mathrm{surf}}}\right)^{4}$
(75)

with $v_{\mathrm{surf}}=\sqrt{\mathcal{G}M_{\mathrm{tot}}/R_{\mathrm{cap}}}$
the characteristic surface velocity. The rate at which planetesimals are
removed from the disc is then

$\dot{M}_{\mathrm{plan}}=\dot{M}_{\rm core}+\dot{M}_{\mathrm{ejec}}.$ (76)

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

####  4.3.3 Feeding zone

Figure 9: Illustration of the procedure to separate planetesimals’ feeding
zones when zones would otherwise overlap. The horizontal axis represent 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_{\mathrm{mid2,3}}$ and $a_{\mathrm{mid3,4}}$
are the edges of the new 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

$R_{\mathrm{feed}}=bR_{\mathrm{H}}.$ (77)

For a planet on a circular orbit, conservation of the Jacobi energy implies
that
$b=\sqrt{12+4/3(\tilde{e}_{\mathrm{plan}}^{2}+\tilde{\imath}_{\mathrm{plan}}^{2})}$
(e.g. Hayashi et al. 1977). So, in a quiescent disc with
$\tilde{e}_{\mathrm{plan}},\tilde{\imath}_{\mathrm{plan}}\ll 1$,
$b=2\sqrt{3}\approx 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_{\mathrm{peri}}-R_{\mathrm{feed}}$ to
$r_{\mathrm{apo}}+R_{\mathrm{feed}}$, with $r_{\mathrm{peri}}$ and
$r_{\mathrm{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

$a_{\mathrm{mid}}=\frac{a_{\mathrm{out}}\sqrt{M_{\mathrm{in}}}+a_{\mathrm{in}}\sqrt{M_{\mathrm{out}}}}{\sqrt{M_{\mathrm{in}}}+\sqrt{M_{\mathrm{out}}}}$
(78)

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

####  4.3.4 Core radius

To obtain the radius of the core (and its density), we apply a methodology
similar to Mordasini et al. (2012b). This model also takes into account 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,
i.e. 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

$\rho(P)=\rho_{0}+cP^{n}.$ (79)

We include three different materials: iron, silicates (perovskite, MgSiO3) and
ice, whose parameters $\rho_{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\text{\,}\mathrm{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, e.g., silicates is not very large (Seager et al.
2007).

###  4.4 Atmospheric escape

During the evolutionary phase, i.e., after the dissipation of the gaseous
disc, planets at small distances of their host star
($\sim$0.1\text{\,}\mathrm{au}$$) receive intense XUV stellar irradiation,
which will drive atmospheric escape. This effect is specially 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 of the planets
radius. Due to the low density of gas, the presence of an envelope will result
a significant increase of the planets’ sizes even if the envelope mass is only
on a percent level of the total planet mass. Bare cores are thus clearly
separated from object that retain a gaseous envelope, and a gap is observed in
the distribution of planetary radii (Owen & Wu 2013; Lopez & Fortney 2013; Jin
et al. 2014; Chen & Rogers 2016; Fulton & Petigura 2018).

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

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

$\dot{M}_{\mathrm{env,e}}=\epsilon\frac{\pi
F_{\mathrm{EUV}}R^{3}_{\mathrm{base}}}{\mathcal{G}M_{\mathrm{tot}}}$ (80)

where $F_{\mathrm{EUV}}$ is the EUV flux, $R_{\mathrm{base}}$ the radius of
the photoionisation base, calculated as in Murray-Clay et al. (2009), and
$\epsilon=0.3$ is the heating efficiency, taken as in Murray-Clay et al.
(2009).

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

$\dot{M}_{\mathrm{env,rr}}\sim
4\pi\rho_{\mathrm{s}}c_{\mathrm{s}}R_{\mathrm{s}}^{2}$ (81)

at the sonic point $R_{\mathrm{s}}$, which is calculated the same way as
$R_{\mathrm{acc}}$. Here
$c_{\mathrm{s}}=\sqrt{k_{\mathrm{B}}T/(m_{\mathrm{H}}/2)}$ is the isothermal
sound speed of ionised gas with $T=${10}^{4}\text{\,}\mathrm{K}$$. The density
can be related to the one at the ionisation base, where $\tau=1$, with

$\rho_{\mathrm{s}}\sim\rho_{\mathrm{base}}\exp{\left[\frac{\mathcal{G}M_{\mathrm{tot}}}{R_{\mathrm{tot}}c^{2}_{\mathrm{s}}}\left(\frac{R_{\mathrm{tot}}}{R_{\mathrm{s}}}-1\right)\right]}.$
(82)

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

$\frac{F_{\mathrm{UV}}}{h\nu_{0}}\sigma_{\nu_{0}}n_{\mathrm{0,base}}\sim
n^{2}_{\mathrm{+,base}}\alpha_{\mathrm{rec}}$ (83)

with $n_{\mathrm{0,base}}$ the density of neutrals at the base,
$h\nu_{0}=$20\text{\,}\mathrm{e}\mathrm{V}$$,
$\sigma_{\nu_{0}}=$6\text{\times}{10}^{-18}\text{\,}{\mathrm{cm}}^{2}$(h\nu_{0}/$13.6\text{\,}\mathrm{e}\mathrm{V}$)^{-3}$,
$\alpha_{\mathrm{rec}}=$2.7\text{\times}{10}^{-13}$$, and
$\rho_{\mathrm{base}}=n_{\mathrm{+,base}}m_{\mathrm{H}}$.

The model also include 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.

###  4.5 Initial conditions

The simulation begin with a predetermined number of embryos whose initial mass
is $M_{\mathrm{emb,0}}=${10}^{-2}\text{\,}\mathrm{M}_{\oplus}$$ (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_{\mathrm{in}}$ and
$40\text{\,}\mathrm{au}$. The starting location zone is slightly more extended
that in the previous studies, where the upper boundary was set to
$20\text{\,}\mathrm{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
120 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 shortcomming of the model will be addressed in
future evolutions of the model (Voelkel et al. 2020), where the evolution of
the dust, pebbles and planetesimal and embryo formation is followed.

##  5 Dynamical evolution: orbital migration, N-body interaction, 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 produce a torque back on
the planet, leading to planetary migration. At the same time, convergent
migration can result in capture in mean-motion resonances or orbital
destabilisation. Hence migration and dynamical evolution must be performed
together to capture all the effects.

###  5.1 Planetary migration

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

####  5.1.1 Type I migration

For Type I migration, our model follows the approach of Coleman & Nelson
(2014). This includes the torques formulation from Paardekooper et al. (2011),
modified to take into account 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

$\Gamma_{1}=F_{\mathrm{L}}\Gamma_{\mathrm{L}}+F_{\mathrm{e}}F_{\mathrm{i}}\left(\Gamma_{\mathrm{c,baro}}+\Gamma_{\mathrm{c,ent}}\right),$
(84)

with

$\displaystyle\Gamma_{\mathrm{c,baro}}$ $\displaystyle=$
$\displaystyle\Gamma_{\mathrm{hs,baro}}F(p_{\nu})G(p_{\nu})+\Gamma_{\mathrm{c,lin,baro}}(1-K(p_{\nu}))$
(85) $\displaystyle\Gamma_{\mathrm{c,ent}}$ $\displaystyle=$
$\displaystyle\Gamma_{\mathrm{hs,ent}}F(p_{\nu})F(p_{\chi})\sqrt{G(p_{\nu})G(p_{\chi})}+$
(86)
$\displaystyle\Gamma_{\mathrm{c,lin,baro}}\sqrt{(1-K(p_{\nu}))(1-K(p_{\chi}))},$

where $\Gamma_{\mathrm{L}}$, $\Gamma_{\mathrm{hs,baro}}$,
$\Gamma_{\mathrm{hs,ent}}$, $\Gamma_{\mathrm{c,lin,baro}}$ and
$\Gamma_{\mathrm{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_{\mathrm{L}}$ provides the reduction of the Lindblad torque for eccentric
or inclined orbits following Cresswell & Nelson (2008), with

$F_{\mathrm{L}}^{-1}=P_{\mathrm{e}}+\left(\frac{P_{\mathrm{e}}}{|P_{\mathrm{e}}|}\right)\left(0.07\hat{\imath}+0.085\hat{\imath}^{4}-0.08\hat{e}\hat{\imath}^{2}\right)$
(87)

and

$P_{\mathrm{e}}=\frac{1+\left(\frac{\hat{e}}{2.25}\right)^{1/2}+\left(\frac{\hat{e}}{2.84}\right)^{6}}{1-\left(\frac{\hat{e}}{2.02}\right)^{4}}.$
(88)

Here, $\hat{e}=e/h=e/(H/r)$ and $\hat{\imath}=i/h=i/(H/r)$ are the planet’s
orbital eccentricity and inclination scaled by the disc’s aspect ratio
$h=H/r$. $F_{\mathrm{e}}$ and $F_{\mathrm{i}}$ provide the reduction of the
corotation torques due to eccentricity and inclination (Bitsch & Kley 2010).
We use

$F_{\mathrm{e}}=\exp{\left(-\frac{e}{h/2+0.01}\right)}$ (89)

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

$F_{\mathrm{i}}=1-\tanh{(\hat{\imath})}$ (90)

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

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

$\tau_{\mathrm{e}}=\frac{t_{\mathrm{wave}}}{0.78}\left(1-0.14\hat{e}^{2}+0.06\hat{e}^{3}+0.18\hat{e}\hat{\imath}^{2}\right)$
(91)

and

$\tau_{\mathrm{i}}=\frac{t_{\mathrm{wave}}}{0.544}\left(1-0.3\hat{\imath}^{2}+0.24\hat{\imath}^{3}+0.14\hat{e}^{2}\hat{\imath}\right),$
(92)

where

$t_{\mathrm{wave}}=\left(\frac{M_{\star}}{M_{\mathrm{planet}}}\right)\left(\frac{M_{\star}}{\Sigma
a_{\mathrm{planet}}^{2}}\right)h^{4}\Omega^{-1}$ (93)

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

####  5.1.2 Type II migration

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

$\frac{3H}{4R_{\mathrm{H}}}+\frac{50\nu
M_{\star}}{M_{\mathrm{planet}}a_{\mathrm{planet}}^{2}\Omega}\leq 1,$ (94)

with $\nu$ 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,

$v_{\mathrm{rad}}=\frac{1}{\Sigma_{\mathrm{g}}\sqrt{r}}\frac{\partial}{\partial
r}\left(\nu\Sigma_{\mathrm{g}}\sqrt{r}\right)$ (95)

(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_{\mathrm{planet}}$ is given by

$\frac{v_{\mathrm{planet}}}{v_{\mathrm{rad}}}=\min{\left(1,\frac{2\Sigma_{\mathrm{g}}a_{\mathrm{planet}}^{2}}{M_{\mathrm{planet}}}\right)}.$
(96)

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

$\Gamma_{2}=\frac{1}{2}M_{\mathrm{planet}}\Omega
a_{\mathrm{planet}}v_{\mathrm{planet}}.$ (97)

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

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

$\tau_{e}=\tau_{i}=\frac{1}{10}|\tau_{\mathrm{a}}|=\frac{1}{10}\frac{a_{\mathrm{planet}}}{|v_{\mathrm{planet}}|}.$
(98)

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

####  5.1.3 Migration map

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

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\text{\,}\mathrm{Myr}$; at this time the disc mass is
$1.46\text{\times}{10}^{-02}\text{\,}\mathrm{M}_{\odot}$. Its outer radius is
$123\text{\,}\mathrm{au}$, so we cut the figure at $200\text{\,}\mathrm{au}$
since there is no migration outside this distance.

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

###  5.2 N-body integration

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

$\dot{x_{i}}=\frac{\partial\mathcal{H}}{\partial p_{i}},\hskip
10.00002pt\dot{p_{i}}=-\frac{\partial\mathcal{H}}{\partial x_{i}}$ (99)

where $x$ denotes the position coordinates, $p$ the momentum coordinates, and

$\mathcal{H}=\sum_{i=0}^{N}\frac{p_{i}^{2}}{2M_{i}}-\mathcal{G}\sum_{i=0}^{N}M_{i}\sum_{j=i+1}^{N}\frac{M_{j}}{\Delta
x_{ij}}$ (100)

is the Hamiltonian of the system, with $\Delta x_{ij}=|x_{i}-x_{j}|$. Here,
the index $i=0$ refers to the central star and $M_{0}=M_{\star}$ while the
subsequent are the planet with $M_{i}=M_{\mathrm{planet,i}}$ so that $N$ is
the number of planets in the system.

However, while $\mathcal{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 close to that of the
original system. The Hamiltonian is divided into three components, so that
$\mathcal{H}=H_{\mathrm{K}}+H_{\mathrm{S}}+H_{\mathrm{I}}$, and

$\displaystyle H_{\mathrm{K}}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left(\frac{p_{i}^{2}}{2M_{i}}-\mathcal{G}\frac{M_{\star}M_{i}}{\Delta
x_{i0}}\right)$ (101) $\displaystyle H_{\mathrm{S}}$ $\displaystyle=$
$\displaystyle\frac{1}{2M_{\star}}\left(\sum_{i=1}^{N}\mathbf{p}_{i}\right)^{2}$
(102) $\displaystyle H_{\mathrm{I}}$ $\displaystyle=$
$\displaystyle-\mathcal{G}\sum_{i=1}^{N}\sum_{j=i+1}^{N}\frac{M_{i}M_{j}}{\Delta
x_{ii}}.$ (103)

Here, $H_{\mathrm{K}}$ represents the unperturbed Keplerian orbits of the
planets about the central star, $H_{\mathrm{S}}$ the kinetic energy of the
star and $H_{\mathrm{I}}$ the interactions between the planets. The separation
into three different Hamiltonians (rather than two) is required because the