scheme uses mixed-centre coordinates (also called “democratic heliocentric”):
heliocentric positions and barycentric velocities. These coordinates are
chosen so that $H_{\mathrm{K}}\gg H_{\mathrm{S}},H_{\mathrm{I}}$, unless two
planets come close together.

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

$H_{\mathrm{I}}\left(\frac{\tau}{2}\right)H_{\mathrm{S}}\left(\frac{\tau}{2}\right)H_{\mathrm{K}}\left(\tau\right)H_{\mathrm{S}}\left(\frac{\tau}{2}\right)H_{\mathrm{I}}\left(\frac{\tau}{2}\right),$
(104)

where the notation $H_{\mathrm{...}}(\tau)$ is used to represent the evolution
under the given Hamiltonian for a step $\tau$. For $H_{\mathrm{I}}$, this
means that the planets receive a kick in velocity due to the interactions with
the other bodies (except the central star). In our case, $H_{\mathrm{I}}$ is
extended to include additional forces representing the effect of the gas disc,
see Section 5.2.1. The evolution under $H_{\mathrm{S}}$ results in a shift
$\tau/(2M_{\star})\sum\mathbf{p}_{i}$ while the evolution under
$H_{\mathrm{K}}$ is a Keplerian motion around the central star for a period
$\tau$.

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

The symplectic integration scheme has a huge advantage in terms of
computational requirements compared to a standard Bulirsch-Stoer method, as
the interaction between the planets, the one part that is
$\mathcal{O}(N^{2})$, is only computed once per step.

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

####  5.2.1 Additional forces

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

$\displaystyle a_{\mathrm{\theta}}$ $\displaystyle=$
$\displaystyle-\frac{v_{\mathrm{\theta}}-v_{\mathrm{K}}}{2\tau_{\mathrm{e}}}+\frac{\Gamma}{M_{\mathrm{planet}}}$
(105) $\displaystyle a_{\mathrm{r}}$ $\displaystyle=$
$\displaystyle-\frac{v_{\mathrm{r}}}{\tau_{\mathrm{e}}}$ (106) $\displaystyle
a_{\mathrm{z}}$ $\displaystyle=$
$\displaystyle-\frac{2v_{\mathrm{z}}}{\tau_{\mathrm{i}}}$ (107)

with $a$ denoting the additional accelerations, $v$ the planet’s velocity
along each direction. Here, $v_{\mathrm{K}}=\Omega r$ is the Kelperian
velocity.

####  5.2.2 Collision detection

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

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

####  5.2.3 Collision treatment

Figure 11: Cumulative distribution of the impactor-to-target mass ratio
$\gamma=M_{\mathrm{tot,2}}/M_{\mathrm{tot,1}}$ for different target mass
ranges (as provided in the legend). Data come from the 100-embryos population
presented in Paper II.

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

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

$E_{\mathrm{impact}}=\max{\left(\frac{1}{2}\mu
v_{\mathrm{imp}}^{2}-E_{\mathrm{acc,core}},0\right)}$ (108)

where
$\mu=M_{\mathrm{tot,1}}M_{\mathrm{tot,2}}/(M_{\mathrm{tot,1}}+M_{\mathrm{tot,2}})$
is the reduced mass, and the indexes 1 and 2 refer to the quantities of the
larger and smaller body respectively. $v_{\mathrm{imp}}$ is the relative
velocity at time of contact. Here

$E_{\mathrm{acc,core}}=\mathcal{G}\frac{M_{\mathrm{tot,1}}M_{\mathrm{core,2}}}{R_{\mathrm{core,1}}+R_{\mathrm{core,2}}}$
(109)

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

The addition of the core mass and luminosity is performed via

$\displaystyle\dot{M}_{\mathrm{core,supp}}$ $\displaystyle=$
$\displaystyle\frac{M_{\mathrm{core,2}}}{\tau_{\mathrm{impact}}\sqrt{2\pi}}\exp{\left(-\frac{1}{2}\left(\frac{t-t_{\mathrm{impact}}}{\tau_{\mathrm{impact}}}-3\right)^{2}\right)}$
(110) $\displaystyle L_{\mathrm{core,supp}}$ $\displaystyle=$
$\displaystyle\frac{E_{\mathrm{impact}}}{\tau_{\mathrm{impact}}\sqrt{2\pi}}\exp{\left(-\frac{1}{2}\left(\frac{t-t_{\mathrm{impact}}}{\tau_{\mathrm{impact}}}-3\right)^{2}\right)}$
(111)

where $t_{\mathrm{impact}}$ is the time of the impact,
$\tau_{\mathrm{impact}}=${10}^{4}\text{\,}\mathrm{y}\mathrm{r}$$ is the time
scale of release taken as in Broeg & Benz (2012). These two terms are added to
the core accretion rate due to planetesimal accretion, and to the luminosity
(Sect. 4.2) used in the internal structure calculation, respectively.

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

###  5.3 Tidal evolution

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

$\frac{\partial a_{\mathrm{planet}}}{\partial
t}=-\frac{9}{2}\sqrt{\frac{\mathcal{G}}{M_{\star}}}\frac{R_{\star}^{5}M_{\mathrm{planet}}}{Q_{\star}}a_{\mathrm{planet}}^{-11/2}$
(112)

(Ferraz-Mello et al. 2008; Jackson et al. 2009; Benítez-Llambay et al. 2011),
where $Q_{\star}=10^{6}$ is the stellar dissipation parameter. It is clear
that his model for the tidal spiralling-in is strongly simplified. It will be
improved in future work along the lines of, e.g., Bolmont & Mathis (2016).

##  6 Terrestrial planet formation

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

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

###  6.1 Setup

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

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

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

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

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

###  6.2 Gravitational interactions

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

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

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

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

###  6.3 Interactions lead to more massive planets

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

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

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

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

###  6.4 Time needed for formation

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

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

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

###  6.5 With outer giant planets

Figure 14: Stacked eccentricity versus distance snapshots of 10 simulations
with each 46 embryos initially. The left column shows the runs with outer
giant planets whereas the right column has no outer giant planets. In each
column, the 10 systems are represented with a different colour for each one.
The bodies are shown by points whose sizes are proportional to their physical
ones. Black crosses show the solar system planets. Figure 15: Same as Fig. 14,
but showing mass versus distance stacked diagram of 10 simulations with each
46 embryos initially. As before, the black crosses show the solar system
planets.

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

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

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

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

###  6.6 Summary for terrestrial planets

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

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

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

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

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

##  7 Giant planets

Figure 16: Formation and evolution tracks of four giant planets with final
masses in the $1/3$ to $2\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$ range
in discs with a single embryo. The top panels present the formation tracks
with total mass $M_{\mathrm{tot}}$ versus distance (times goes toward the top)
and total mass $M_{\mathrm{tot}}$ (solid lines) and core mass $M_{\rm core}$
(dashed lines) versus time. The three panels on the bottom row show the time
dependence of the outer luminosity $L_{\mathrm{tot}}$ (bottom left), the
distance (bottom centre) and the total radius $R_{\mathrm{tot}}$ (bottom
right). For all panels except for the mass versus time (top right), the line
styles denote the phase: dashed lines for the attached phase, solid line for
the detached phase during formation and dash-dotted lines for the evolution
stage. Line widths denote the migration regime, with tick lines for Type I and
think lines for Type II. The legend in the top right panel shows the gas (in
Solar masses) and planetesimals (in Earth mass) disc masses.

Figure 17: Formation and evolution tracks of two groups of four giant planets
with final masses between $2$ and
$10\text{\,}\mathrm{M}_{\textrm{\tiny\jupiter}}$ in discs with a single
embryo. Panel and line descriptions are the same as Fig. 16.

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

$\frac{M_{\mathrm{g}}}{$2\text{\times}{10}^{-3}\text{\,}\mathrm{M}_{\odot}$}=\left(\frac{r_{\mathrm{cut,g}}}{$10\text{\,}\mathrm{au}$}\right)^{1.6}.$
(113)

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

###  7.1 Formation and evolution of Jupiter-mass planets

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

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

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

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

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

Figure 18: Formation and evolution tracks of one giant planet that ends up
being accreted by the star during the evolution stage due to tidal migration.
Top left: total mass $M_{\mathrm{tot}}$ versus distance; Top right: total mass
$M_{\mathrm{tot}}$ (black) and core mass $M_{\rm core}$ (red) versus time;
Bottom left: total luminosity $L_{\mathrm{tot}}$ (black) and the bloating
contribution $L_{\mathrm{bloat}}$ (red) versus time; Bottom centre: distance
versus time; Bottom right: total radius $R_{\mathrm{tot}}$ (black) and outer
(attached phase) or Hill (detached phase) radius (red) versus time.

###  7.2 More massive planets

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

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

###  7.3 Giant planets ending in the star by tidal migration

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

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

###  7.4 Summary for giant planets

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

##  8 Individual systems

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

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

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

####  8.1.1 Low initial solid content

Figure 19: Example of the formation of a planetary system from initially 100
lunar-mass embryos in low gas mass initial mass
$0.023\text{\,}\mathrm{M}_{\odot}$), low metallicity ([Fe/H]=-0.13) disc. The
initial mass of plantesimals is $65.1\text{\,}\mathrm{M}_{\oplus}$. Four
moments in time (in years) are shown. Lines shows the growth tracks in the
semimajor axis-mass plane. Black points show (proto)planets existing at a
given epoch. Gray open circles show the last position of protoplanets that
were accreted by another more massive body in a giant impact. The colours of
the lines are $|\tau_{\rm mig}/\tau_{\rm grow}|=|d\ln{m}/d\ln{a}|$. Figure 20:
Same system as in Fig. 19, but now showing the semi-major axes $a$ of the
planets as a function of time, colour coding in panel (a) the planets’ mass,
in (b) the planetesimal surface density in the planets’ feeding zone, and in
(c) the local gas surface density. Here, the vertical line indicates the
moment of gas disc dissipation. Panel (d) shows mass as a function of time,
colour coding the semi-major axis. Small black circles indicate giant impacts,
by showing the position or mass of the target (the more massive collision
partner) at the moment of the impact.

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

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

At the very beginning, all protoplanets grow as if they were the only bodies
in the disc, not feeling the influence of the other protoplanets. With
increasing mass, the interaction (either directly via N-body interactions) or
indirectly via resonant migration, become important. By
${10}^{5}\text{\,}\mathrm{yr}$, the first dynamical interactions have started
among some of the more massive protoplanets, which is visible as a “jitter” in
some tracks, and two collisions, which are shown by two open gray circles.

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

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

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

Figure 21: Temporal evolution of the eccentricities of the planets of the
system emerging in the low-mass disc shown in Fig. 19. Colours indicate the
planet mass. For better visibility, only planets more massive than
$0.1\text{\,}\mathrm{M}_{\oplus}$ are shown. The curves were also temporally
smoothed such that one mainly sees the mean values instead of the rapid
variations of the eccentricities. The increase of the eccentricities at around
$5\text{\,}\mathrm{Myr}$ when the gas disc dissipates is visible.

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

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

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

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

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

####  8.1.2 High initial solid content

Figure 22: Example of the formation of a planetary system from initially 100
lunar-mass embryos in a high gas mass (initial mass
$0.066\text{\,}\mathrm{M}_{\odot}$), high metallicity ([Fe/H]=0.23) disc. The
initial mass of planetesimals is $432\text{\,}\mathrm{M}_{\oplus}$. The plot
is analogous to Fig. 19, but the y-axis now extends to much higher masses, and
the moments in time that are shown are different. At the end of the simulation
at $20\text{\,}\mathrm{Myr}$, this system contains one close-in sub-Neptunian
planet, three giant planets, and a group of outer very low-mass planets.
Figure 23: Temporal evolution of the system shown in Fig. 22. The four panels
have the same meaning as in Figure 20.

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

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

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

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

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

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

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

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

###  8.2 Overview of the diversity of system architectures

Figure 24: Mass-distance diagrams of specific systems with 100 embryos
initially (panels a to w), which are taken from the nominal population
predicted for a $1\text{\,}\mathrm{M}_{\odot}$ star (NG76). Symbols are as
follows: red points show gas-rich planets where $M_{\mathrm{env}}/M_{\rm
core}>1$. Blue symbols are planets that have accreted some volatile material
(ices) outside of the ice line(s). Green symbols are planets that have only
accreted refractory solids. Open green and blue circles have 0.1$\leq
M_{\mathrm{env}}/M_{\rm core}\leq 1$ while filled green points and blue
crosses have $M_{\mathrm{env}}/M_{\rm core}\leq 0.1$. For all these bodies,
the grey horizontal bars go from $a-e$ to $a+e$. The top left panel with black
crosses shows the solar system. Bodies lost because of collisions or ejections
are shown in light grey. Planets accreted by the central star are show in the
very left of each panel, the ejected ones on the very right and planets that
collided with another (more massive) planet are shown at their last position
on the diagram. The dotted vertical line in each system shows the location of
the ice line. The number after each panel name is the metallicity [M/H] of the
system expressed in dex, while the value on the top right is the initial mass
of the planetesimals disc. Figure 25: Radius-distance diagrams for a subset of
the systems shown in Fig. 24. Here, we focus on systems with multiple low-mass
planets; panel letters correspond to the same system in Fig. 24. In contrast
to Fig. 24, lost planets are not shown for clarity.

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

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

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

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

####  8.2.1 Mass and final number of planets

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

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

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

####  8.2.2 Similarity in the low-mass systems

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

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

####  8.2.3 Composition of the close-in planets

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

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

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

##  9 Summary and conclusions

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

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

The model now includes a multitude of physical processes (see Fig. 2). During
both the formation and evolution phase, these are:

  * •

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

  * •

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

  * •

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

  * •

the evolution of the host star from tabulated stellar evolution tracks
(Baraffe et al. 2015).

During formation, the following processes are included:

  * •

a 1D radial (axis-symmetric) constant $\alpha$-disc model for radial structure
of the protoplanetary gas disc. The effects of internal and external
photoevaporation are included.

  * •

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

  * •

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

  * •

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

  * •

a model for the planetary gas accretion rate in the disc-limited regime that
is based on Bondi- and Hill-type gas accretion in the 2D and 3D case.

  * •

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

  * •

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

During the evolutionary phase we include:

  * •

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

  * •

for close-in planets, the addition of a bloating luminosity modelled with the
empirical relation of Thorngren & Fortney (2018). \- tidal spiral-in because
of stellar tides (Benítez-Llambay et al. 2011). Roche-lobe overflow is also
included.

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

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

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

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

  * •

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

  * •

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

  * •

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

  * •

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

  * •

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

  * •

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

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

###### Acknowledgements.

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

## References

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

##  Appendix A Possible impact of migration on the accretion of planetesimals

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

###  A.1 General considerations

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

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

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

###  A.2 Applicability of Tanaka & Ida (1999)

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

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

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

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

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

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

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

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

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

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

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