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

# Planetary Population Synthesis and the Emergence of Four Classes of
Planetary System Architectures

Alexandre Emsenhuber<EMAIL_ADDRESS>Christoph Mordasini 22 Remo
Burn 44

###### Abstract

Planetary population synthesis is a helpful tool to understand the physics of
planetary system formation. It builds on a global model, meaning that the
model has to include a multitude of physical processes. The outcome can be
statistically compared with exoplanet observations. Here, we review the
population synthesis method and then use one population computed using the
Generation III Bern model to explore how different planetary system
architectures emerge and which conditions lead to their formation. The
emerging systems can be classified into four main architectures: Class I of
near-in situ compositionally ordered terrestrial and ice planets, Class II of
migrated sub-Neptunes, Class III of mixed low-mass and giant planets, broadly
similar to the Solar System, and Class IV of dynamically active giants without
inner low-mass planets. These four classes exhibit distinct typical formation
pathways and are characterised by certain mass scales. We find that Class I
forms from the local accretion of planetesimals followed by a giant impact
phase, and the final planet masses correspond to what is expected from such a
scenario, the ‘Goldreich mass’. Class II, the migrated sub-Neptune systems
form when planets reach the ‘equality mass’ where accretion and migration
timescales are comparable before the dispersal of the gas disc, but not large
enough to allow for rapid gas accretion. Giant planets form when the ‘equality
mass’ allows for gas accretion to proceed while the planet are migrating,
i.e., when the critical core mass is reached. The main discriminant of the
four classes is the initial mass of solids in the disc, with contributions
from the lifetime and mass of the gas disc. The distinction between mixed
Class III systems and Class IV dynamically-active giants is in part due to the
stochastic nature of dynamical interactions, such as scatterings between giant
planets, rather than the initial conditions only. The breakdown of system into
classes allows to better interpret the outcome of a complex model and
understand which physical processes are dominant. Comparison with observations
reveals differences to the actual population, pointing at limitation of
theoretical understanding. For example, the overrepresentation of synthetic
super Earths and sub-Neptunes in Class I systems causes these planets to be
found at lower metallicities than in observations.

##  1 Introduction

The advent of exoplanet discovery projects, such as the Kepler space telescope
(Borucki et al.,, 2010), the HARPS survey (Mayor et al.,, 2011), or the
California Planet Search (Howard et al., 2010a, ; Rosenthal et al.,, 2021),
has led to the discovery of a large number of exoplanets. As of writing, more
than $5000$ exoplanets have been confirmed.111According to
http://exoplanet.eu, the $5000$ exoplanets mark was passed on about 27 April
2022.

These exoplanets are diverse. Many of them have no counterpart in the Solar
System. Indeed, the first discovered exoplanet around a main sequence star, 51
Peg b (Mayor and Queloz,, 1995), is a Jupiter-like planet orbiting its host
star with a period of only $4.2\text{\,}\mathrm{d}$. Later, a broad class of
planets whose mass or size is between those of the Earth and Neptune has been
unveiled (e.g. Fulton et al.,, 2017). The large number of discovered planetary
systems permits to have statistical constraints, both on the planet properties
and their correlations with other factors (such as stellar properties or
stellar cluster environment).

But while we have thousands of evolved planets, there are only scant
observations of planets that are still forming: the PDS 70 system (Keppler et
al.,, 2018; Müller et al., 2018a, ; Haffert et al.,, 2019) and AB Aur b
(Currie et al.,, 2022), for instance. This limited number of observations
means that we have poor direct observational constraints on the physical
processes that occur during planet formation.

Parallel to the observational campaigns, there has been a wealth of works
aiming at describing the physical processes that are thought to occur during
the formation and later evolution of planetary systems (Benz et al.,, 2014;
Baruteau et al.,, 2016; Raymond and Morbidelli,, 2022). However, the different
processes cannot be tested individually because of the lack of direct data we
mentioned. Nevertheless, we can leverage the large number of the end products
of planetary formation (that is, the observed exoplanets) to provide indirect
evidence on how these processes work. This is the basic idea of planetary
population synthesis: to bridge the relative lack of observations by using the
predictions of theoretical works as tests of their effectiveness at
reproducing the statistical characteristics of the discovered exoplanets. This
requires—among other things—a global (meaning that it includes the many
interlinked physical processes occurring at each stage) and simplified (for
computational requirements) formation and evolution model, which has to
predict the necessary quantities for comparison with observations.

An assumption of planetary population synthesis is that the diversity of
exoplanets stems from the diversity of protoplanetary discs (the initial
conditions for the planetary formation process and for the numerical
calculations) and the inherent chaotic nature of planetary formation
(especially for N-body interactions), rather than different models. The global
model (for a pioneering work, see Ida and Lin, 2004a, ) must thus be able to
reproduce the diversity of planetary systems, from terrestrial planets to
giants. The interlinking of many processes included in the global model is
source of complexity of the results. As such it may become advantageous to
classify the final systems according to a certain scheme. In this work, we
explore avenues opened by such a classification.

One focus of this work is to link the birth environment to the final planetary
systems. In population synthesis, the diversity of protoplanetary discs is
represented by the initial conditions. The link between environment and
architecture is achieved by backtracking the initial conditions that lead to
the different categories.

A global model is not only used to make predictions about the final systems;
it can also be used to study many theoretical aspects of planet formation.
This includes how the various processes interact with each other, which one is
dominant in the different stages, or for mechanisms leading to the different
categories of systems. This will be the second focus of this work. Here, we
use formation pathways and analytical mass scales and match them with the
final planet masses obtained in the numerical simulations for the different
classes to determine which process is responsible for shaping the final
planetary systems.

The final aim of planetary population synthesis is to adapt the model so that
the correspondence between synthetic and observed systems is improved. To
achieve this, it is necessary to understand where the model must be corrected.
We apply this to the established correlations with stellar metallicity and
discuss pathways for how the model can modified to improve fidelity.

This work is divided as follows. We start in Sect. 2 by discussing in more
details the concepts of planetary population synthesis followed by the history
and different models in Sect. 3. The remainder of this work will focus on the
link between the properties of the protoplanetary discs and the diversity of
planetary systems and the emergence of four classes of planetary system
architectures. The methodology is presented Sect. 4 and the corresponding
results in Sect. 5.

##  2 Principles of Planetary Population Synthesis

Figure 1: Principle and flow chart of planetary population synthesis method
and its links to the overall theory of planet formation and to the
construction of astronomical instrumentation. Extended from Mordasini et al.,
(2015).

The arching goal of planetary population synthesis is to improve our
understanding of planet formation by finding the theoretical model that
provides the best match of the observed exoplanets from the properties of
protoplanetary discs. This also enables to make predictions about the planets
that have not been observed (yet). Such an endeavour encompasses many aspects,
from theory to observations.

Figure 1 show how these elements are linked. The starting point are the
detailed models of individual processes like disc structure and evolution,
planetary solid and gas accretion, orbital migration, planet-planet
interaction and so on. They inform the global end-to-end models which combine
the essence of many such detailed models into one big model. This gives end-
to-end models the capability to directly predict observable properties of
planetary systems based on the initial conditions of the planet formation
process, which are the protoplanetary disc properties.

From observations it is known that not all protoplanetary discs have the same
properties, in terms of essential disc properties like mass, size, or lifetime
(e.g. Manara et al.,, 2022; Pascucci et al.,, 2022; Miotello et al.,, 2022).
Statistical disc observations allow to derive probability distributions for
the initial and boundary conditions of planet formation.

Next, one can draw from these distribution and each time calculate the
formation of one synthetic planetary system. By doing this many times, one
obtains a synthetic planetary population. To compare quantitatively, one more
step is needed: to apply the same observational detection bias that affects
the actual population. Applying biases of different observational techniques
on the same synthetic population makes it possible to compare different
aspects of the theory, to use the full wealth of observations, and to avoid
over-fitting a specific aspect.

Then, we can compare statistically and quantitatively the actual observed
population with the biased synthetic population. Many different quantities can
be compared, like frequencies of planet types, distributions of quantities
like the mass, radius, orbital distance or the eccentricity, but also
architectures of the systems, and correlations with stellar properties like
stellar mass, age, and metallicity.

It is clear that given the complexity of planet formation, there will not be a
perfect match between the synthetic result and the multitude of concurrently
available observational constraints. These differences reveal theoretical sub-
models we need to improve, where we need to change parameters or even abandon
an old approach altogether, replacing it with new theoretical approaches. By
this loop, the understanding of planet formation is improved.

In some cases, a model is found which is satisfactory at least in some aspect.
This is of interest on the instrumental side. We can now go back to the full
underlying population and study what happens if we can improve the sensitivity
and capabilities of the instrumentation. These insights can be used to
identify and build the next generation of astronomical instrumentation that
yields the most compelling constraints for the understanding of planet
formation. This then goes into the observed population, and via this second
loop, planet formation theory is again improved.

There are numerous requirements to perform planetary populations synthesis.
From the observational side, one needs to have extensive, statistical surveys
of both protoplanetary discs and exoplanets. These survey further need to have
well-characterised biases, so that models and observations can be reliably
compared in a quantitative way. On the theoretical side, the model used for
population synthesis must satisfy multiple conditions: 1) being able to model
end-to-end formation and evolution of planetary systems, that is, from points
in time of the observations of protoplanetary discs and exoplanets to ages of
many gigayears, 2) predict the necessary observable quantities for comparison
with exoplanet observations, and 3) having sufficient performance so that many
systems can be modelled.

###  2.1 Coupling between physical processes

Conditions 1) and 2) imply that such a model must include many processes and
that their coupling is also important. Major processes occurring in planetary
formation (disc evolution, planetary accretion, orbital migration) all occur
on similar timescales and affect each other. This fact must be reflected in
models, where the processes themselves cannot be treated in isolation. This is
the basis of global models, that the physical processes must be treated
concurrently and not in isolation.

###  2.2 Low-dimensional approach

The addition of condition 3) means that the model cannot treat directly first
principles, i.e, solve the conservation laws in three dimensions. With the
current computational capabilities, detailed, three-dimensional hydrodynamical
simulations (for the protoplanetary disc, for instance) have too large
requirements to be performed for the entire evolution a protoplanetary disc,
let alone for an entire population of planetary systems with thousands of
systems. Instead, we have to rely on more simple, lower-dimensional models.

Thus, while not using them directly, a model for planetary population
synthesis relies on more sophisticated models to provide the parametrisation
required for the low-dimensional approach. In turn, the simplified global
models are able to test the specific models in a broader context encompassing
the entire picture of planetary system formation rather than what is
explorable with an isolated model, with, for instance, only orbital migration.

Another avenue for simplification of more detailed simulations is to use
heuristics derived from those. This can be used for instance to mimic N-body
interactions (Ida et al.,, 2013). In the same direction, surrogate models
obtained via machine learning could be used to reduce the computational
requirements of more complex models (e.g. as Emsenhuber et al.,, 2020, for
collisions), though this has not been yet widely used in the context of
planetary population synthesis.

###  2.3 Protoplanetary disc demographics: setting the stage for planetary
formation

The model provides only the time evolution of planetary systems. We thus need
to provide the initial conditions, which should be chosen such that they lead
to the reproduction of as many disc observations as possible. The youngest
observable targets are class 0 or class I discs, where the best measurements
of their masses and radii come from dust emission (Tychoniec et al.,, 2018;
Williams et al.,, 2019; Tobin et al.,, 2020). In contrast, measurements of
class II discs (Ansdell et al.,, 2016; Pascucci et al.,, 2016) indicate that
they are not suited as initial conditions as the dust content at this stage is
already smaller than the estimated core masses of observed exoplanets (Manara
et al.,, 2018). Dust mass measurements are however subjects to large
uncertainties, which stem from at least in part assumed temperatures and
opacities (e.g. Bergez-Casalou et al.,, 2022) that do not necessarily reflect
the reality. The obtained values should therefore be taken with some caution.

Disc dust-to-gas ratio are usually taken to be the same as stellar
metallicities (Gáspár et al.,, 2016), for which distributions have been
determined from stellar spectroscopy (Santos et al.,, 2005; Petigura et al.,,
2017). The distribution of disc lifetimes is estimated from measurements of
the IR-excess of young stars, but the interpretation is complex (coevality,
stellar ages; Haisch et al.,, 2001; Mamajek,, 2009; Fedele et al.,, 2010;
Ribas et al.,, 2014; Richert et al.,, 2018; Michel et al.,, 2021).

This approach assumes that the characteristics of the protoplanetary discs
observed today are the same as those that formed the observed exoplanets. This
a generally-accepted principle, but it should be kept in mind that observed
exoplanets usually formed billions of years ago which is not negligible in
relation to the age of the galaxy (Adibekyan et al.,, 2014). Furthermore,
current disc statistics mainly stem from small, nearby star forming regions in
which environmental effects differ from the birthplaces of most extrasolar
planets (e.g. van Terwisga et al.,, 2020, 2022, for more massive regions).

###  2.4 Reproduction of exoplanets characteristics

Determining which is the best model is achieved by statistically comparing to
specific systems, such as the Solar System of TRAPPIST-1, and the exoplanet
population as a whole. In the latter case, there are three main classes of
characteristics derived from observations that synthetic populations can be
compared to

First, there are planet-level statistics: their number, their distribution in
observed parameters (distance, mass, radius, etc.). Hot-Jupiters (planets with
a period $\lesssim$10\text{\,}\mathrm{d}$$) are found around $0.5$ to
$1\text{\,}\mathrm{\char 37\relax}$ of stars (Howard et al., 2010b, ; Mayor et
al.,, 2011), compared to $10$ to $20\text{\,}\mathrm{\char 37\relax}$ for all
giants within 5 to $10\text{\,}\mathrm{au}$ (Cumming et al.,, 2008; Johnson et
al.,, 2010). The frequency of giant planets increases with distance (Dong and
Zhu,, 2013; Santerne et al.,, 2016) followed by a peak in location near the
iceline (Fernandes et al.,, 2019; Fulton et al.,, 2021) and a decrease beyond
$3$–$10\text{\,}\mathrm{au}$ (Bryan et al.,, 2016; Nielsen et al.,, 2019). The
planet mass function derived from radial velocity observations exhibits two
peaks: one for sub-Neptunes at about $10\text{\,}\mathrm{M_{\oplus}}$ and one
for giants at $1$ to $2\text{\,}\mathrm{M_{\textrm{\jupiter}}}$ separated by a
relative ‘desert’ of intermediate-mass planets (Mayor et al.,, 2011), but the
significance of this ‘desert’ is debated (Bennett et al.,, 2021). Concerning
planet radii, there are two peaks at $1.3\text{\,}\mathrm{R_{\oplus}}$ (super
Earths) and $2.4\text{\,}\mathrm{R_{\oplus}}$ (sub-Neptunes) separated by a
gap (Fulton et al.,, 2017). Giants have a similar radius of $\sim
1\text{\,}\mathrm{R_{\textrm{\jupiter}}}$ leading to another local maximum
(Mordasini et al., 2012b, ).

Second, there are correlations with stellar properties. One of the first
discovered is the correlation between the frequency of giant planets with
stellar metallicity (Gonzalez,, 1997; Santos et al.,, 2004; Fischer and
Valenti,, 2005). Low-mass planets seems to be more common around low-mass
stars (Mulders et al.,, 2015) while it is the opposite for giants (Bonfils et
al.,, 2013; Ghezzi et al.,, 2018).

Finally, we have system-level statistics. Here, we have that roughly
$50\text{\,}\mathrm{\char 37\relax}$ of the giant-planet-hosting systems have
multiple planets (Bryan et al.,, 2016). Further multiple planets in the same
system tend to have similar masses (Millholland et al.,, 2017), radii, and
orbital spacing (Weiss et al.,, 2018; Mishra et al.,, 2021; Weiss et al.,,
2022). This creates the interesting situation of similarity within the overall
diversity.

##  3 History and the different models

The pioneering series on extrasolar planetary population synthesis is the work
of Ida and Lin. The first work, Ida and Lin, 2004a , presented their
population synthesis workflow, which contains the key ingredients that we have
been discussing. They predicted the presence of a ‘planetary desert’ because
of fast growth at the onset of gas runaway (masses between $30$ and
$100\text{\,}\mathrm{M_{\oplus}}$). The second work of the series, Ida and
Lin, 2004b , reproduced the ‘metallicity effect’ (Gonzalez,, 1997; Santos et
al.,, 2004). The third work, Ida and Lin, (2005), addressed planet formation
around stars of different masses. Type I gas driven migration became part of
the model in the fourth work, Ida and Lin, 2008a  and the authors found that
migration rates have to be reduced by at least a factor 10 to make planet
occurrence rate similar to observations. Alternatively, a ‘dead zone’ near the
iceline could reduce migration and enhance accretion so that Type I migration
does not have to be artificially reduced (Ida and Lin, 2008b, ). The last two
works of the series focused on multi-planetary systems. Finally, they added a
semi-analytical approach to study planet-planet interactions without having to
perform full N-body calculation, which was applied to super-Earths (Ida and
Lin,, 2010) and giant planets (Ida et al.,, 2013). There has been some later
work, for instance Ida et al., (2018), who postulated that a reduced Type II
gas-driven migration is better able to match the location very massive planets
(super-Jupiters) near $1\text{\,}\mathrm{au}$ while less massive planets
(Saturn to Jupiter mass) being found across a wider range of distances. Also,
Miguel et al., (2020) apply an updated code to the case of late M and brown
dwarfs.

Another series of population synthesis works were performed by Mordasini and
collaborators. The first series, based on the formation model of Alibert et
al., (2004, 2005) considered only a single embryo per disc and found that
migration had to be reduced by factor between 10 and 100 to avoid planets
falling into the star (Mordasini et al., 2009a, ). Outcomes were then compared
to exoplanets discovered by radial velocity at the time (Mordasini et al.,
2009b, ). This work was later extended to different stellar masses (Alibert et
al.,, 2011). The last work of that series, (Mordasini et al., 2012a, ),
focused on the disc properties required to form giant planets. The later work
took to directions, with one branch focusing on predicting planet radii and
luminosities for comparison with transit and direct imaging surveys (Mordasini
et al., 2012c, ; Mordasini et al., 2012b, ). The other branch focused on the
formation of multi-planetary systems (Alibert et al.,, 2013), better
representation of solids accretion (Fortier et al.,, 2013), and the planets’
chemical composition (Thiabaud et al.,, 2014, 2015). These two branches were
later merged to compute multi-planet populations capable of predicting the
different properties for various observational techniques for both Solar-mass
stars (Emsenhuber et al., 2021a, ; Emsenhuber et al., 2021b, ) and M-dwarfs
(Burn et al.,, 2021). Related works include the correlation between close-in
super Earths and distant giant planets (Schlecker et al., 2021a, ), the
reproduction of the formation model model using machine learning techniques
(Schlecker et al., 2021b, ), and the reproduction of peas-in-a-pod structure
and comparison with Kepler (Mishra et al.,, 2021).

A series of population synthesis works by Alessi and Pudritz focuses on the
effect of migration traps, disc properties, and chemistry. The first work
(Alessi and Pudritz,, 2018) focuses on the conditions needed to obtain a
separation between hot and warm Jupiters in terms of dead zone and envelope
metallicity. Then, the authors study the effects of the discs initial
properties (Alessi et al., 2020b, ) and the predictions on core and envelope
composition (Alessi et al., 2020a, ). Their latest work focuses on the
inclusion of magnetised disc winds as transport mechanism (Alessi and
Pudritz,, 2022).

There are also studies based on the core accretion paradigm assuming that the
solids are accreted in the form of pebbles, but comparatively less than with
planetesimals. One is Chambers, (2018), who finds that low-viscosity and other
parameters that provide for short accretion timescales are best able to match
exoplanets observations. Another work by Liu et al., (2020) focus on low-mass
stars with a single-embryo model going from M-stars down to brown dwarfs. With
multiple interacting bodies, Lambrechts et al., (2019) focus on the assembly
of terrestrial or super-earth systems while Bitsch et al., (2019) benchmark
their models against the giant planet population. There is also a series of
work focusing on the stellar environment: Ndugu et al., (2018) investigated
how the background heating by nearby stars, while Ndugu et al., (2022)
investigated how stellar encounters affect the formation, for instance by disc
truncation. A comparison of pebble and planetesimal based models in a
population synthesis work was done in Brügger et al., (2020).

Two recent Protostars and Planets VII review chapters address planetary
population synthesis. Weiss et al., (2022) investigates compact multi-
planetary systems while Drazkowska et al., (2022) focuses on the comparison
between planetesimals and pebble accretion.

The base assumption that all extrasolar planets have formed by a bottom-up
formation pathway (i.e. by some variant of the core accretion paradigm) might
not hold (Morales et al.,, 2019; Currie et al.,, 2022). Therefore, several
population synthesis models were also built based on the gravitational
instability paradigm (Forgan and Rice,, 2013; Forgan et al.,, 2018;
Nayakshin,, 2015; Müller et al., 2018b, ). These models simulate the evolution
of the gas disc, use fragmentation criteria, and follow the evolution of the
emerging gas clumps via tidal downsizing, migration, and contraction.

##  4 Methods

In this work, we use the Generation III Bern model of planetary formation and
evolution, as it was presented in earlier works (Emsenhuber et al., 2021a, ).
Thus, we only provide a summary of the model here. The interested reader is
referred to relevant works for an in-depth description of the model and the
population synthesis methods (Emsenhuber et al., 2021a, ; Emsenhuber et al.,
2021b, ; Burn et al.,, 2021).

The Bern model (Alibert et al.,, 2004, 2005, 2013; Mordasini et al., 2009a, ;
Mordasini et al., 2012c, ) couples the formation and evolution stages. The
formation stage follows the evolution of a viscously-accreting protoplanetary
disc, the dynamical state of planetesimals in the disc, concurrent accretion
of solids and gas from the disc by the protoplanets, gas-driven planetary
migration, and dynamical (N-body) interactions between multiple protoplanets
in a each system. The evolution stage tracks the thermodynamic evolution
(cooling and contraction), atmospheric escape, and tidal migration of each
planet individually. Figure 2 shows the most important sub-modules of the
model and the key quantities that they are exchanging (the model is entirely
described in Emsenhuber et al., 2021a, ). The individual sub-modules can be of
considerable complexity, but numerically use—with the exception of the N-body
integrator—a low-dimensional approach (for example the planetary interiors are
assumed to be 1D spherically symmetric). This approximation must be made to
keep the ability to simulate the planetary systems from their origins to the
present day, which is only possible with low-dimensional models because of
computational time constraints.

Figure 2: Sub-models in the Generation III Bern Model of planet formation and
evolution. It is an example of a comprehensive end-to-end model as often used
in planetary population synthesis calculations. Arrows indicate the most
important exchanged quantities. Colors give the phase when a sub-module is
active. For computational time constraints, the ‘Formation’ phase is assumed
to last for $100\text{\,}\mathrm{Myr}$ only (extended from the
$20\text{\,}\mathrm{Myr}$ of Emsenhuber et al., 2021b, ), because simulating
with an N-body integrator thousands of planetary system over billion year-
timescales is currently prohibitively expensive. In the subsequent ‘Evolution’
phase, the planets evolve individually to a $5\text{\,}\mathrm{Gyr}$ age.
Updated from Emsenhuber et al., 2021a .

###  4.1 Formation stage

The gas component of the protoplanetary disc is modelled with a 1D radial,
axisymmetric representation that is viscously accreting (Lynden-Bell and
Pringle,, 1974). Viscosity is computed following the standard $\alpha$
parametrisation (Shakura and Sunyaev,, 1973) with a fixed value of
$\alpha=$2\text{\times}{10}^{-3}$$, while the vertical structure is computed
in a vertically-integrated approach (Nakamoto and Nakagawa,, 1994). We include
additional sink terms to the evolution to account for internal and external
photoevaporation, following Clarke et al., (2001) and Matsuyama et al.,
(2003), respectively, and the accretion by the protoplanets. The solids in the
disc are assumed to be in the form of planetesimals of a fixed radius
$r=$300\text{\,}\mathrm{m}$$. They are treated as a 1D radial, axisymmetric
surface density (similar to the gas) with a dynamical state (eccentricity and
inclination). The dynamical state is evolved under the stirring by the other
planetesimals, the protoplanets, and the damping due to the presence of a gas
disc (Ohtsuki et al.,, 2002; Chambers,, 2006; Fortier et al.,, 2013).

Planet formation follows the core accretion paradigm (Perri and Cameron,,
1974; Mizuno et al.,, 1978; Mizuno,, 1980). Initially, a given number of
embryos are randomly placed in the disc with a uniform spacing in the
logarithm of the distance, so that they are separated by a given amount of
Hill radii, as found in N-body simulations of embryo formation (Kokubo and
Ida,, 2002). They will accrete planetesimals, which is taken to be in the
oligarchic regime (Ida and Makino,, 1993; Inaba et al.,, 2001; Thommes et
al.,, 2003). The gas envelope structure is initially in contact with the
surrounding disc and accretion occurs thanks to radiating away the
gravitational binding energy released during this process, i.e., on a Kelvin-
Helmholtz timescale (Pollack et al.,, 1996; Lee and Chiang,, 2015). The gas
accretion rate is thus determined by solving the 1D spherically symmetric
internal structure equation of the envelope (Bodenheimer and Pollack,, 1986).
The presence of a gaseous envelope also enhances the capture radius of
planetesimals, which we compute following the procedure of Inaba and Ikoma,
(2003). As a protoplanet core grows, its Kelvin-Helmholtz timescale decreases
(Ikoma et al.,, 2000). The gas accretion rate thus increases until it becomes
larger than what can be supplied by the disc. When this occurs, the envelope
can no longer maintain the contact with the disc and contracts (Bodenheimer et
al.,, 2000). In this stage, the structure equations are used to compute the
planet radius instead, while the gas accretion is set by the supply from the
disc (Mordasini et al., 2012c, ; Mordasini et al., 2012b, ).

Dynamical interactions between the different protoplanets are modelled by
means of the mercury N-body package (Chambers,, 1999). Also, type I gas-driven
migration is computed according to Coleman and Nelson, (2014) while type II
and the transition between the two regimes follow Dittkrist et al., (2014).
Migration is applied as additional forces in the N-body.

###  4.2 Evolution stage

After a fixed time of $100\text{\,}\mathrm{Myr}$ (compared to
$20\text{\,}\mathrm{Myr}$ in Emsenhuber et al., 2021b,  to increase the
duration of the N-body interactions), the model transitions to the evolution
stage, where the planet are individually evolved to $5\text{\,}\mathrm{Gyr}$.
In addition to the normal thermodynamic evolution (cooling and contraction) of
the interior which is found by solving 1D internal structure equations, this
stage also tracks XUV-driven atmospheric escape (Jin et al.,, 2014; Jin and
Mordasini,, 2018) and inward migration due to stellar tides (Ferraz-Mello et
al.,, 2008; Benítez-Llambay et al.,, 2011).

###  4.3 Analysis: mass scales

Because of the numerous interaction and feedback (see the many arrows in Fig.
2), global models can produce numerical results of considerable complexity
which can be difficult to interpret. To understand what leads to the final
shape of the final systems, it is therefore helpful to analytically
investigate a number of mass scales of different processes that take place
during the formation stage. The purpose is to compare these analytical mass
scales to the numerical formation tracks to understand which process is linked
to the different features.

The first two mass scales are related to the accretion of planetesimals only.
Temporally, the first process that occurs is this accretion of planetesimals.
Absent radial motion caused by the gas disc (orbital migration) or
interactions with other bodies, the embryos can only accrete from nearby
planetesimals, which will halt once the local reservoir is depleted. This
leads to the planetesimal isolation mass (Lissauer,, 1987),

$M_{\mathrm{iso}}=\frac{\left(2\pi
ba^{2}\Sigma_{\mathrm{p}}\right)^{\frac{3}{2}}}{\sqrt{3M_{\star}}}\sim$3.6\text{\times}{10}^{-2}\text{\,}\mathrm{M_{\oplus}}$\left(\frac{a}{$1\text{\,}\mathrm{au}$}\right)^{3}\left(\frac{\Sigma_{\mathrm{p}}}{$10\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$}\right)^{\frac{3}{2}}\left(\frac{M_{\star}}{$1\text{\,}\mathrm{M_{\odot}}$}\right)^{-\frac{1}{2}},$
(1)

with $b\sim 10$ the full width of the feeding zone in Hill radii of the planet
$r_{\mathrm{H}}=a\left(M/(3M_{\star})\right)^{1/3}$, $a$ the local position,
$\Sigma_{\mathrm{p}}$ the local surface density of planetesimals, and
$M_{\star}$ the mass of the central star. At large separation, however, the
accretion timescales become so slow that the isolation mass cannot be reached
in any reasonable time. There, growth is rather stalled at a time comparable
to the dispersal of the protoplanetary disc, which halts the damping of the
planetesimals random velocities and causes in turn the accretion rate to
strongly decrease. We can use the analytic accretion timescales in the
oligarchic regime (Ida and Makino,, 1993) to estimate a growth mass (see
Thommes et al.,, 2003, Eq. 14 for the derivation),

$\displaystyle M_{\mathrm{grow}}$ $\displaystyle\sim$
$\displaystyle\dot{M}\tau_{\mathrm{acc}}=k_{\mathrm{grow}}\frac{t_{\mathrm{disc}}^{3}M_{\star}^{1/2}\Sigma_{\mathrm{g}}^{6/5}\Sigma_{\mathrm{p}}^{3}}{\rho_{m}^{4/5}\rho_{M}a^{3/2}m^{2/5}(H/a)^{6/5}}$
(3) $\displaystyle\sim$
$\displaystyle$6.8\text{\,}\mathrm{M_{\oplus}}$\left(\frac{t_{\mathrm{disc}}}{$1\text{\times}{10}^{6}\text{\,}\mathrm{yr}$}\right)^{3}\left(\frac{\Sigma_{\mathrm{p}}}{$10\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$}\right)^{3}\left(\frac{a}{$1\text{\,}\mathrm{au}$}\right)^{-\frac{3}{2}}\left(\frac{M_{\star}}{$1\text{\,}\mathrm{M_{\odot}}$}\right)^{\frac{1}{2}}\left(\frac{\Sigma_{\mathrm{g}}}{$2400\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$}\right)^{6/5}$
$\displaystyle\left(\frac{\rho_{m}}{$1\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$}\right)^{-4/5}\left(\frac{\rho_{M}}{$3.2\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$}\right)^{-1}\left(\frac{m}{$1\text{\times}{10}^{18}\text{\,}\mathrm{g}$}\right)^{-1/10}\left(\frac{H/a}{0.05}\right)^{-6/5}$

with a constant prefactor
$k_{\mathrm{grow}}=\frac{\pi^{3}3^{21/15}}{4^{3/5}}\left(\frac{C_{D}b}{C_{e}}\right)^{6/5}G^{3/2}\sim
11.9\sqrt{G^{3}}$. Here, $G$ is the Newtonian gravitational constant and
reasonable values are for the drag coefficient $C_{D}\sim 1$ (Adachi et al.,,
1976), for the viscous stirring factor $C_{e}\sim 40$ (Ida and Makino,, 1993;
Ohtsuki et al.,, 2002), and for the width of the feeding zone $b\sim 10$.
Further, $\Sigma_{\mathrm{g}}$ is the local gas surface density, $H$ the disc
scale height often combined to the aspect ratio $H/a\sim 0.1$,
$t_{\mathrm{disc}}$ the lifetime of the gaseous disc, $\rho_{m}$ and
$\rho_{M}$ are the bulk densities of the planetesimals and the protoplanet,
and $m$ is the mass of one planetesimal. Note that while $M_{\mathrm{grow}}$
is much larger than the other scales at short distances (within
$1\text{\,}\mathrm{au}$), it decreases most rapidly with distance for
realistic disc and planetesimal surface density profiles. It should be noted
that $M_{\mathrm{grow}}$ only makes sense in the outer parts of the disc
(outside of $\sim 10$ AU) where it holds that $M_{\mathrm{grow}}\leq
M_{\mathrm{iso}}$.

The next two mass scales are related to orbital migration. As the mass of a
protoplanet increases, its oligarchic planetesimal accretion rate decreases,
whereas its type I migration rate increases. At one point migration becomes
faster than growth. In a mass-distance diagram this means that planetary
formation track bend from a vertical (upward) direction inward and become more
horizontal. This happens at what we call the ‘equality mass’, where the
migration timescale is equal to the accretion timescale. It can be estimated
by equating the planetesimals accretion timescale from (Mordasini,, 2018, Eq.
10, or equivalently Thommes et al.,, 2003) with the type I migration timescale
neglecting non-linear corotation torques (see Tanaka et al.,, 2002, Eq. 70,
but with non-isothermal torques from Paardekooper et al.,, 2010). This results
in

$M_{\mathrm{eq}}=k_{\mathrm{eq}}\frac{\Sigma_{\mathrm{p}}^{3/4}M_{\star}^{5/4}(H/a)^{6/5}}{\Upsilon^{3/4}\Sigma_{\mathrm{g}}^{9/20}\rho_{m}^{1/5}\rho_{M}^{1/4}a^{3/4}m^{1/10}}\,,$
(4)

where $k_{\mathrm{eq}}\simeq
3^{11/10}\pi^{3/4}\left(\frac{bC_{D}}{2C_{e}}\right)^{3/10}\simeq 4.2$,
$\Upsilon\simeq\gamma^{-1}(2.5+1.7\beta_{\mathrm{T}}-0.1\beta_{\Sigma}$) is an
order unity prefactor for the Lindblad torque which depends on the adiabatic
index of the gas $\gamma$, the radial slopes of the surface density set to
$-\beta_{\Sigma}$, and the one of temperature to $-\beta_{\mathrm{T}}$
(Paardekooper et al.,, 2010) and $m$ is the mass of an individual
planetesimal. Entering typical values, this becomes

$\displaystyle M_{\mathrm{eq}}$
$\displaystyle=$1.2\text{\,}\mathrm{M_{\oplus}}$$
$\displaystyle\left(\frac{\Upsilon}{2.23}\right)^{-3/4}\left(\frac{M_{\star}}{$\text{\,}\mathrm{M}_{\odot}$}\right)^{5/4}\left(\frac{\Sigma_{\mathrm{p}}}{$10\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$}\right)^{3/4}\left(\frac{m}{${10}^{18}\text{\,}\mathrm{g}$}\right)^{-0.1}$
(6)
$\displaystyle\left(\frac{H/a}{0.05}\right)^{6/5}\left(\frac{\Sigma_{\mathrm{g}}}{$2400\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$}\right)^{-0.45}\left(\frac{a}{$1\text{\,}\mathrm{au}$}\right)^{-3/4}\left(\frac{\rho_{m}}{$1\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$}\right)^{-1/5}\left(\frac{\rho_{M}}{$3.2\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$}\right)^{-1/4}.$

However, the classical type I Lindblad torques can be counteracted by the co-
rotation torque (Ward,, 1991; Masset,, 2001). The co-rotation region exerts a
torque on the planet if gradients in specific vorticity and entropy are
present. To not balance out those gradients and saturate the corotation
torque, the gradients need to be restored by the disc’s viscosity faster than
the time it takes for a parcel of gas to perform a full horseshoe orbit
(Paardekooper et al.,, 2011; Dittkrist et al.,, 2014). The width of the
corotation region is given by

$\Delta a_{\rm HS}=C_{\rm
HS}\frac{a}{\gamma^{1/4}}\sqrt{\frac{Ma}{M_{\star}H}},$ (7)

where $C_{\rm HS}\approx 1.1$ is a numerical factor determined from numerical
experiments (Paardekooper et al.,, 2010) and $\gamma$ is the heat capacity
ratio or adiabatic index. Since $\Delta a_{\rm HS}$ increases with planetary
mass, we can determine a critical mass at which the corotation torque
saturates by equating the libration and the viscous diffusion timescale across
the corotation region (Hellary and Nelson,, 2012)

$\frac{8\pi a}{3\Omega_{\mathrm{K}}\Delta a_{\rm
HS}}\stackrel{{\scriptstyle}}{{=}}\frac{(\Delta a_{\rm HS})^{2}}{\nu}.$ (8)

Solving for the mass results in what we call the saturation mass

$M_{\mathrm{sat}}=\left(\frac{8\pi\alpha}{3}\right)^{2/3}\sqrt{\frac{\gamma}{C_{\rm
HS}^{4}}}M_{\star}\left(\frac{H}{a}\right)^{7/3}.$ (9)

This mass scale the upper limit where planets can be trapped between outward
and inward migration regions (migration traps). Above this value, planets
always migrate towards the star. However, planets below this mass scale are
not necessarily trapped into convergence zones, as that depends on the
relative strength of the Lindblad and corotation torques. We note that this
mass scale depends on $H/a$, which is distant- and time-dependent. This leads
to time-dependent predictions, as the aspect ratio decreases as a function of
time because of decreasing viscous heating, and stellar evolution (e.g.
Baraffe et al.,, 2015). The assumption of viscous dissipation as the radial
transport method in the disc does also play a major role here, as the
inclusion of magnetically-driven disc winds combined a lower viscous
dissipation would severely reduce the corotation torques (though it may also
produce large outward corotation torques under certain conditions, Kimmig et
al.,, 2020).

The fifth mass scale is linked to giant impacts. This important mass scale
does not directly depend on orbital migration but may depend indirectly on it,
as migration leads to a radial redistribution of the protoplanets: the final
stage of the inner planets’ formation is dominated by giant impacts. This
phase can lead to planet masses significantly larger than the planetesimal
isolation mass. Assuming that the damping of the protoplanets’ eccentricities
by remaining planetesimals is weak, the resulting mass can be estimated
(Goldreich et al.,, 2004). To obtain this mass scale, we assume that the
random velocity are given by mutual gravitational interactions between the
planets only. Therefore, the derivation of this mass scale is similar to the
isolation mass, but the extent of the feeding zone is modified to include the
eccentric excursions from apastron to periastron of the growing planets, and
the building blocks are no longer directly the planetesimals, but larger
protoplanets with a surface density $\Sigma_{\mathrm{P}}$. The captial ‘P’
here stands in contrast to that of the planetesimals $\Sigma_{\mathrm{p}}$.
Here, we assumed for an order-of-magnitude estimate that, locally, all
planetesimals have been accreted by the protoplanets, for example bodies that
have reached the local planetesimal isolation mass. In detail, this assumption
does not hold in regard to the Solar System where there are planetesimal-like
bodies that remain until today, but we are here interested at the overall mass
budget, which is dominated by the planets. This assumption is also justified a
posteriori from the fact that all planetesimals are accreted in the
terrestrial region (Sect. 5.3). The radial distribution of these protoplanets
may differ from the initial planetesimals depending on whether orbital
migration lead to a significant redistribution. Let $m$ be the mass of the
planet which is given by

$m=2\pi a2\Delta a\Sigma_{\mathrm{P}}$ (10)

for a radial extent of the feeding zone $\Delta a=ea$ if the orbit is
eccentric enough such that the width used for the isolation mass (Eq. (1)) is
much smaller than the radial excursion of the planet, that is when $\Delta
a\gg r_{\mathrm{H}}$ or $e\gg\left(M/(3M_{\star})\right)^{1/3}$. Given
multiple interacting bodies, the planets’ random velocity is approximately
given by the planets’ escape velocity $v_{\mathrm{esc}}=\sqrt{2Gm/R}$ which
then leads to an eccentricity of $e\approx v_{\mathrm{esc}}/v_{\mathrm{K}}$.
Thus,

$m=4\pi a^{2}\Sigma_{\mathrm{P}}\sqrt{\frac{2Gm/R}{GM_{\star}/a}}\,,$ (11)

where we can assume a constant bulk density $\rho_{\mathrm{p}}$ for each
planet and solve for $m$ to get the Goldreich mass

$M_{\mathrm{Gold}}=16a^{3}\Sigma_{\mathrm{P}}^{3/2}\left(\frac{2\pi^{7}a^{3}\rho_{\mathrm{P}}}{3M_{\star}^{3}}\right)^{1/4}\,.$
(12)

If the effect of a previous orbital migration of the protoplanets can be
neglected, we can set $\Sigma_{\mathrm{P}}=\Sigma_{\mathrm{p}}$. It should be
noted that the Goldreich mass is inherently a mass scale that is not strictly
local, as it is obtained via a series of giant impacts. This leads to a radial
mass exchange, as can be seen in the simulations below. It is therefore useful
to calculate $M_{\mathrm{Gold}}$ not with local quantities only, but take an
average over the region where the giant impact phase takes place (like for
example inside of the iceline, as it is done in the next section).

Finally, two additional mass scales are important for giant planet formation:
the critical (core) mass for gas runaway accretion (Mizuno et al.,, 1978;
Mizuno,, 1980; Piso and Youdin,, 2014) and the mass where planets switch from
faster type I to slower type II migration (Crida et al.,, 2006). Regarding the
former, it is important to directly calculate it with the structure equations,
as semi-analytical expression can give grossly incorrect results (Alibert and
Venturini,, 2019).

##  5 Results

We conduct our main analysis for a synthetic planet population that is based
on the nominal population that was presented in Emsenhuber et al., 2021b
(NG76). However, in the simulation shown here (NG76longshot), the formation
stage which includes N-body interactions was extended from 20 to
$100\text{\,}\mathrm{Myr}$. As a simplification, the model starts with all the
solids are in the form of planetesimals, with a steeper radial slope than the
gas to account for planetesimal formation simulations (power-law index of
$\beta_{\mathrm{s}}=-1.5$ following Voelkel et al.,, 2020). The solid
component is thus always in the form of planetesimals and does not account for
the presence of dust in evolved discs. At the same time, $100$ lunar-mass
embryos are placed in each of the $1000$ simulated planetary systems. This
population contains a variety of planetary systems, from terrestrial planet
systems to multiple giants, as we will show in this section.

###  5.1 Four classes of system architectures

To analyse the formation patterns of different kinds of planetary systems and
determine how the diversity of protoplanetary discs and environment affect the
formation, we start by classifying the final systems into a number of classes
of system that show similar properties. We base the classification on mainly
the mass-distance diagram plus additionally the bulk composition. For the
latter, we distinguish five compositions: Jovian (H/He dominated in mass),
Neptunian (volatile (ice)-rich with H/He), water worlds (volatile-rich without
H/He), hyterran (silicate-iron core with H/He envelope) and terrestrial
(silicate-iron without H/He or ices).

By visually inspecting each of the $1000$ synthetic planetary systems in the
population, we find that the system architectures emerging in the simulations
can be divided into four main classes. These are, ordered in increasing
typical mass of formed planets, 1) (compositionally) ordered Earths and ice
worlds systems, 2) migrated sub-Neptune systems, 3) mixed systems with low-
mass and giant planets, broadly speaking akin to the Solar System, and 4)
dynamical active giants. The preponderance of the different classes is as
follows: $57.6\text{\,}\mathrm{\char 37\relax}$ are part of Class I,
$22.1\text{\,}\mathrm{\char 37\relax}$ of the system are found to be part of
Class II, $12.3\text{\,}\mathrm{\char 37\relax}$ are part of Class III, and
finally $8.0\text{\,}\mathrm{\char 37\relax}$ are found to be part of Class
IV. As was already discussed by Emsenhuber et al., 2021b , most systems
contain only low-mass planets (Class I and II) while only about
$20.3\text{\,}\mathrm{\char 37\relax}$ of the systems are in classes that
contain giant planets (Class III and IV). However, these numbers are sensitive
to the choice of initial conditions, which themselves are subject to
uncertainties (Sect. 2.3).

We find that the extension of the formation stage to
$100\text{\,}\mathrm{Myr}$ (compared to $20\text{\,}\mathrm{Myr}$ in
Emsenhuber et al., 2021b, ) affects mostly the Class-I systems. The final
architecture of the systems in the other classes is only seldom affected by
this change. In the following subsections, we will review the formation
pathways of each of the four classes and provide example of system
architecture in each of them.

Our classification of planetary systems is based on a human classificator. As
becomes clear below, the classes found this way correspond to fundamentally
different temporal formation pathways which are governed by different physical
processes (like giant impacts, orbital migration, or gas runaway accretion).
These processes in turn manifest themselves through the different mass scales
introduced in Sect. 4.3. Our approach thus stresses the temporal emergence of
different system architectures and the (theoretical) physical mechanisms
included in the model leading to them. A distinct approach was recently
presented by Mishra et al., 2023a ; Mishra et al., 2023b  who also classified
the population used here in four distinct types. Their analysis is focused on
observable systems, by applying observational biases on the final systems.
They also employed a number of coefficients that can be calculated numerically
out of the physical properties to identify the classes instead of a human
classificator. This approach uses the current state of the systems and is more
observationally-driven, in contrast to the temporal/evolutionary aspect and
the underlying theoretical physical processes that are of import in the
approach used here. Interestingly enough, there is nevertheless a significant
degree of convergence in the outcomes of the two approaches, as discussed in
Sect. 5.3.2.

####  5.1.1 Class I: Compositionally ordered Earths and ice worlds systems

Figure 3: One example of formation tracks for a system that belongs to Class
I. The solid grey lines show the tracks of all protoplanets. The ones that
remain through the end ($5\text{\,}\mathrm{Gyr}$) are additionally marked with
colour symbols at the end of the track, according to their composition. Green
(blue) full circles show rocky (icy) bodies without H/He, while green (blue)
empty circles show rocky (icy) bodies with a H/He envelope, and red full
circles are gas-rich bodies (more than $50\text{\,}\mathrm{\char 37\relax}$
H/He by mass). Black stars and empty circles show target and impactor of giant
impacts, with dotted black lines connecting the impact partners. The
analytical isolation (violet, Eq. (1)), Goldreich (orange, Eq. (12)),
saturation (green, Eq. (9)), equality (cyan, Eq. (4)), and growth (pink, Eq.
(3)) mass scales are shown. The latter is only shown where it is lower than
isolation mass. In Class I systems, the masses in the inner system are
governed by the (mean) Goldreich’s mass $M_{\mathrm{Gold}}$, i.e., via the
final giant impact phase. The saturation mass is evaluated at two different
times. Figure 4: Final architecture of four example systems in Class I.
Coloured dots show the position in the mass-semimajor axis plane of different
planet composition types as indicated in the legend. The three numbers in the
top left give the initial conditions [M/H], initial gas disc mass
[$\mathrm{M_{\odot}}$], and initial mass of planetesimals
[$\mathrm{M_{\oplus}}$]. The systems are ordered according to the latter.
Horizontal black bars indicate the periastron to apastron distance, reflecting
orbital eccentricity. Grey full circles show the last position of protoplanets
that were lost. The ones lost by accretion onto the star are put on the left
y-axis of each panel and those ejected on the right one. The remaining grey
points correspond to protoplanets that were accreted by another more massive
protoplanet. The point then shows the mass and semi-major axis just before
this accretion event.

Class I contains systems of low-mass planets whose composition remain well-
ordered in the sense that there are rocky (silicate/iron) planets inside and
icy (volatile-rich) planets outside. One example of the formation tracks of
such systems is provided in Fig. 3. For the embryos that start inside the
(water) ice line, growth is divided in three phases. The first phase is in-
situ accretion of planetesimals, until the planets reach the planetesimal
isolation mass. As the isolation mass is much lower than the equality mass in
this area, migration is negligible. Only the few most massive planets, those
just beyond the iceline, become massive enough to undergo some migration. The
second phase is some inward migration without accretion. Migration is due to
the push by the more massive outer planets, meaning that large resonant
convoys form. The tight spacing of the embryos combined with the short
accretion timescales mean that all the planetesimals inside the iceline have
been accreted, rendering further solid growth impossible. The planets are also
too low mass for gas accretion. In addition, as the planetesimal accretion
timescale increases with distances, the inner planets have finished growth by
planetesimals by the time the outer planets begin to migrate, leading to
migration at constant mass. The third phase, which starts when the gas disc
vanishes, is growth by giant impacts, which allows the planets to reach masses
that are nearly two orders of magnitude larger than the isolation mass, as was
discussed by Emsenhuber et al., 2021a . The masses are not large enough to
either be substantially affected by migration (when the gas disc is still
present) or for system-wide destabilisation during the giant impact phase.
Thus, the distinction between rocky and icy planets remains, as material is
not reshuffled in the system; the only effect is an overall shift to closer
distance, which lead to the presence of icy planets down to about
$1\text{\,}\mathrm{au}$. This is inside by about 1 – $2\text{\,}\mathrm{au}$
of the location of the iceline (around $3\text{\,}\mathrm{au}$ in these
systems; its location is set from the initial disc profile and then kept
constant). In the inner system where many giant impacts have happened, the
final planet masses are comparable to, or somewhat larger than the Goldreich
mass. Its mean value (mean over the region inside of the iceline) is shown in
the figure. In the outer region, there is a sharp decrease of the planet
masses beyond the location where the growth timescale is larger than the disc
lifetime. It should be noted that we have simulated the formation phase which
in particular includes the accretion of planetesimals and the N-body
interaction including giant impacts for $100\text{\,}\mathrm{Myr}$. Further
growth may occur over very long timescales (Ida and Lin, 2004a, ). The
extension from $20\text{\,}\mathrm{Myr}$ in Emsenhuber et al., 2021b  allows
for continued growth of the terrestrial planets by giant impacts. Class-I
systems presented here have therefore on average less terrestrial planets (but
more massive) than those in the original population from Emsenhuber et al.,
2021b . However, only a small fraction of systems would have their
architecture affected in such a way that this would alter their
classification.

Four examples of final Class I architectures are shown in Fig. 4. Here we see
that the systems exhibit a similar behaviour. The inner planets tend to have
slightly increasing mass with distance. This follows the behaviour of the
planetesimal isolation mass, which for the assumed MMSN-like initial
planetesimal surface density profile $\propto a^{-1.5}$ increases with
distance. However, the mass increase is not as strong as the radial dependency
of the isolation mass that is shown with the purple line in Fig. 3. This
implies that the giant impact phase is able to reduce the effect of the
planetesimals isolation mass; rather the Goldreich mass is the relevant
scaling for the final masses at the end of the giant impact phase. One trend
that can be observed from the systems shown here is that there is an
anticorrelation between the number of planets and their masses. The systems in
the top two panels each have three terrestrial planets whose masses are $\sim
1\text{\,}\mathrm{M_{\oplus}}$ while the systems in the two bottom panels have
one or two terrestrial planets, whose masses are $\sim
4\text{\,}\mathrm{M_{\oplus}}$. We finally observe the peas-in-a-pod effect
(Mishra et al.,, 2021; Weiss et al.,, 2022) across the entire population,
where planets in a system have similar masses. The characteristic mass scale,
the Goldreich mass, depends on the solid building block surface density (Eq.
(12)). Therefore, the higher the initial mass of planetesimals in a disc
(initial mass in $\mathrm{M_{\oplus}}$ given in the figure), the higher the
Goldreich mass.

####  5.1.2 Class II: Migrated sub-Neptunes systems

Figure 5: Example of formation tracks and final architecture of a system that
belong to Class II. In this class, the governing mass scale in the inner
system is linked to orbital migration and is thus given by the saturation mass
$M_{\rm sat}$ where planets are released from migration traps (green lines,
evaluated at $t=0$ and $t=$0.1\text{\,}\mathrm{Myr}$$) and/or the equality
mass $M_{\rm{eq}}$ where the timescales of oligarchic growth and type I
orbital migration are equal (cyan line). Figure 6: Final architecture of four
systems of Class II. The plot is analogous to Fig. 4. In contrast to Class I,
all planets have a volatile-rich composition, and the mass is a decreasing or
constant function of orbital distance. One also notes the absence of planets
with H/He at small orbital distances. The numerous low-mass planets outside of
10 AU might not have yet settled into their final state after 100 Myr which is
the N-body integration time covered in our simulations.

Class II are systems with more massive planets, yet no giants. There are
usually 5 to 10 planets whose mass is larger than
$1\text{\,}\mathrm{M_{\oplus}}$. In contrast with Class I, most planets are
icy, even in the inner region. Another contrast is that the planet masses are
either independent of, or anticorrelated with their orbital distance, opposite
to Class I.

We show one example of the formation tracks of a system in this class in Fig.
5. It reveals that the inner system is populated by icy planets that
originally formed beyond the iceline and then migrate inward at constant mass
following similar tracks as the ‘horizontal branch’ planets already identified
in Mordasini et al., 2009a . The terrestrial planets that also originally from
in the system are captured in mean-motion resonances, migrate in resonant
convoys, and end up either pushed into the star or get destabilised, which
leads to a collision with an icy planets; they usually do not survive. In our
scheme, the collision between a terrestrial and an icy planet result in an icy
planet, as we use a threshold of $1\text{\,}\mathrm{\char 37\relax}$ of
volatiles by mass in the core to classify planets as icy and planets forming
beyond the iceline have a volatile mass fraction up to
$59.5\text{\,}\mathrm{\char 37\relax}$, which is radially variable and depends
on whether the individual volatile species are in condensed or gaseous form.
One way to distinguish planets that had a giant impact during the migration
phase is to see if they have a H/He envelope. Giant impacts with the
terrestrial planets will release a sufficient amount of energy to strip the
envelope from the icy planets (e.g. Denman et al.,, 2020). This is the case
for instance of the first and third innermost planets. The loss of the
envelope of the second innermost planet is due to a runaway bloating process
leading to Roche lobe overflow. In this process, the radius increases because
the heating due to bloating, which further increases bloating. This occurs for
close-in low-mass planets with a high envelope mass fraction (like for runaway
Roche lobe overflow due to tidal heating, Valsecchi et al.,, 2015, or
atmospheric escape, Baraffe et al.,, 2004; Thorngren et al.,, 2022). We have
verified that the envelope is also lost by XUV-driven atmospheric
photoevaporation in the absence of Roche lobe overflow.

Unlike systems in Class I, there is sufficient mass beyond the iceline to form
sub-Neptunes within the lifetime of the protoplanetary disc. The planets
forming just outside the iceline are thus able to reach the equality mass (the
cyan line in Fig. 5) early enough for the gas disc to be massive enough to
exert significant torques and thus migrate inward. We note that some planets
experience strong inward migration even below the saturation mass (green
lines); this is because the saturation mass decreases over time along with the
gas surface density. As consequence, the saturation mass at the time where
planet begin to rapidly migrate is lower than the value shown in the figure.
Further, the planets do not experience outward migration below the saturation
because they are not in a region where the corotation torque is large enough
to overcome the Lindblad torques. The planet on the horizontal branch have
masses that are too low for rapid gas accretion, meaning that they do not
reach the critical (core) mass. The formation of Class II systems relies on
migration. Thus, the physical process responsible for the radial transport of
gas is of special importance: in the case of magnetically-driven disc winds
combined with low viscosity, migration would likely be reduced (Suzuki et
al.,, 2016), and part of the Class II systems could behave more like Class I
systems, but with more massive planets. However, the corotation torque depends
on the viscosity (Masset,, 2001) and the transport mechanism (Kimmig et al.,,
2020), so this would have to be verified by hydrodynamical simulation.

We further provide the final architectures of four more systems in this class
in Fig. 6. These systems show a common pattern of planets of masses between
about $2$ and $10\text{\,}\mathrm{M_{\oplus}}$ planets from the inside to
about $2\text{\,}\mathrm{au}$. This is a relic of the initial growth beyond
the iceline. The equality mass (cyan line in Fig. 5) is a strongly decreasing
function of distance, which is mainly linked to the increase of the
planetesimal accretion timescale with distance. Planets just outside of the
iceline will thus be able to reach larger masses before being subject to
significant migration. Further, as the the growth timescale increases with
distance, the order of the planets is retained through the migration stage;
with the earlier and more massive planets undergoing stronger migration than
the more distant ones. Therefore, the more massive planets origination from
just outside of the iceline migrate in first, followed later on by less
massive planets originating from further away. The resulting effect of mass
decrease with distance is further compounded by the growth by giant impact
during the migration. As we pointed out above, the innermost planets are H/He-
free while more distant ones have some gaseous envelope, and this is the
consequence of giant impacts stripping the envelope of the inner planets. The
envelope represent only a fraction of the total mass, so the total mass change
by giant impact is usually dominated by the gain of core mass rather than the
loss of the envelope. Planet masses is not a strict discriminant between Class
I and Class II systems. For instance, the innermost planets in the system
shown on the top left panel of Fig. 6 are about
$10\text{\,}\mathrm{M_{\oplus}}$, while that in the bottom right panel of Fig.
4 are slightly more massive. The difference lies in that the planets in the
Class II system migrate further than the ones in Class I, and this is related
to the disc lifetimes of the respective systems. We shall discuss this effect
in greater details in Sect. 5.4.

####  5.1.3 Class III: mixed systems with low-mass and giant planets

Figure 7: Example of formation tracks and final architecture for a system that
belongs to Class III. In this class and in contrast to Class II, a protoplanet
(with typical origin outside the water iceline) manages to grow massive enough
to start runaway gas accretion, i.e., it reaches the critical core mass. It
then passes into slower Type II migration and grows to a giant planet mass.
The final mass of the giant planet is in this system similar to the one of
Jupiter, but is found at about 0.7 AU. The system also contains several lower
mass planets in- and outside of the giant planet. During its gas accretion
phase, the giant accreted a 3 $M_{\oplus}$ planet. 8 protoplanets were ejected
out of the system and one fell into the star, as visible by grey lines
crossing the right and left y-axes, respectively. Figure 8: Final
architectures of four systems of Class III, analogous to Fig. 4. The diversity
is very large regarding both the orbital configuration, mass, and number of
the giant planets but also regarding the composition of the inner low-mass
planets which can be of a Neptunian, water world, hyterran or terrestrial
type.

Class III systems have giant planets, along with inner less massive planets.
They also have outer lower mass planets. Here again, we provide formation
tracks of one system in this class in Fig. 7. The initial formation stage is
similar to system in Class II, with the embryos inside the iceline accreting
in-situ up to the isolation mass, while the embryos beyond the iceline begin
to migrate before reaching the isolation mass. The contrast with systems in
Class II is that the planets beyond the iceline continue to grow by gas
accretion during their inward migration and end up undergoing runaway gas
accretion. At this point, the track of the forming giant planet bends strongly
upwards. The reason is the higher (equality) mass these proto-giants have
(close to $10\text{\,}\mathrm{M_{\oplus}}$, instead of a few Earth masses), as
discussed below. Once the planet undergo runaway gas accretion, the migration
rate decreases in part because in the inner disc, the migration rate is
limited by the ratio between the planet mass and that of the local gas disc.
As a consequence, the giant planets migrate less than the sub-Neptunes in
Class II, and leave the inner part of the system free from perturbation. The
inner system behaves with features of both Class I and II. The part resembling
to Class I is that there remain terrestrial planets that grow by giant
impacts, while the part resembling to Class II is the presence of two
Neptunian planets inside of the giant planets that underwent migration.

The initial disc of the system presented in Fig. 7 is more massive than that
in Fig. 5. As a consequence, the isolation mass is also larger while the
growth timescale $\tau_{\mathrm{grow}}$ is smaller, meaning faster accretion.
This pushes the equality mass to larger values, allowing the planets to reach
larger masses for a given amount of migration compared to the Class II system.
The larger core masses allow the embryos to bind proportionally even more gas
(as the Kelvin-Helmholtz cooling timescales decreases for larger body masses,
Ikoma et al.,, 2000), allowing them to undergo runaway gas accretion on a time
that is less than that of the inward migration.

Four final systems in this class are shown in Fig. 8. As they show, the inner
planets can be any combination of rocky or icy bodies. In systems with both
rocky and icy inner planets, the same compositional ordering as in Class I is
retained. Rocky planets have a formation pattern similar to those in Class I,
while the inner icy planets have similar formation patterns as for those in
Class II, migrating after having accreted from outside the iceline. This also
means that the giant planets are not necessarily the first embryos just beyond
the iceline. Outside of the giants, we see groups of several planets with
masses of a few $\mathrm{M_{\oplus}}$. They consist of mostly ices and some
H/He, reminiscent of the Solar System ice giants. They still grow slowly from
remaining planetesimals. Over long timescales, they might coalesce to form
more massive ice giants.

One common feature between the different systems in this category is that they
are dynamically cool. Planet eccentricities remain relatively low, even when
multiple giant planets are present. Still, there are low and mid-mass planets
in the outer regions that get on wide orbits or ejected (but not the giant
planets themselves). However, this is limited to protoplanets outside the
position of the giants.

####  5.1.4 Class IV: Dynamically-active giants

Figure 9: Formation tracks of a system that belongs to Class IV. The
concurrent formation of several (proto)giant planets (bodies reaching the
critical core mass for runaway gas accretion) in spatial proximity lead to
violent dynamical interactions among these massive planets, strongly
disturbing the lower-mass bodies in the system via collisions and ejections.
In the end, in this example only one very massive, distant and eccentric giant
planet remains.  Figure 10: Final architectures of four systems of Class IV,
analogous to Fig. 4. From the points next to the left and right y-axis, one
notes the significant number of planets that have either been ejected out of
these systems or that have collided with the host star. In the system in the
bottom left, collisions among giant planets have also occurred.

Systems in Class IV also have giants, like Class III; the difference lies in
that they do not contain inner low-mass planets. Overall, they have a small
number of remaining objects (sometimes even only one or two), with possible
low-mass bodies at large distance. Systems in this class tend to have very
massive giant planets ($\sim 10\text{\,}\mathrm{M_{\textrm{\jupiter}}}$),
although this is not always the case.

As for the other classes, we provide one example of the formation history of a
system in this class in Fig. 9. It reveals that the system formed a total of
three giant planets, two of which have subsequently been ejected, leaving a
single giant on a wide and eccentric orbit. These planets begin their
formation similar to the system in Class III discussed above, but the giants
dynamically interact and get destabilised. The instability arises if several
protoplanets in close-proximity start gas runaway accretion within a short
time interval. It should be noted that our model assumption of inserting all
embryos at one moment (at t=0) may artificially lead to such situations. If
the early phases of solid growth (from dust via pebbles to planetesimal and
embryos) are included, bodies may emerge more gradually (Voelkel et al.,,
2022). The dynamical instability caused by the giant planets affects the other
bodies in the system, starting with close-by bodies, some of which get
accreted by the innermost giant. The others are sent on very eccentric orbits,
ending up either colliding with the star or ejected.

We also provide the final state of four systems in this class in Fig. 10.
While they show a pattern of one or more giant planets, these planets have a
diversity of configurations, from compact pairs inside $1\text{\,}\mathrm{au}$
(top left panel) to a single body at $\sim 100\text{\,}\mathrm{au}$ (top right
panel, the same one as in Fig. 9). In addition, other low-mass objects at
large distances may remain. These are the remainder of the outer embryos which
have not been affected by the instabilities of the inner giant planets. For
instance, none remain in system 15 (top right panel) because the planet was
sent to an eccentric orbit with large separation while two other giants have
been ejected (as revealed by the grey dots at $100\text{\,}\mathrm{au}$);
however, in the other systems there are no giant planets either ejected or on
wide orbits, leaving the outer region free from perturbation that would
destabilise these bodies.

One of the main difference between systems of Class III and IV are the
dynamical instabilities. The presence of very massive giants does favour the
occurrence of planet-planet scatterings, and thus systems with the most
massive giants are more likely to be in Class IV.

###  5.2 Key mass scales in planetary systems

Figure 11: Gaussian kernel density estimates of the planet mass function for
planets within $5\text{\,}\mathrm{a}\mathrm{u}$. For each line, the same
bandwidth of 0.3 dex was chosen. The isolation (Eq. (1)) and Goldreich (Eq.
(12)) mass scales are plotted only for initially rocky planets (dotted lines).
Most of the planets in the sample were however initially icy. For comparison
to those, we show with the magenta curve the smaller of the growth $M_{\rm
grow}$ (Eq. (3)) or the equilibrium mass $M_{\rm eq}$ (Eq. (4)). Due to the
lack of statistics (only five planets), we omitted the rocky planets and their
mass estimates for Class IV.

In Sect. 4.3 we introduced a number of analytical key mass scales which can be
used to understand the main mechanisms of planetary formation. We then used
them to understand the formation tracks of individual systems. Here, we expand
this analysis to the population as a whole. Figure 11 shows the comparison of
the various mass scales with all the planets in the different classes of
planetary systems.

It is interesting to see that the typical mass scales can predict to an order
of magnitude or better the modelled distribution of solid-dominated masses.
For the icy planets, we expect them to be limited by either running out of
time to grow or by starting to migrate. The smaller of the growth or equality
mass is shown in Fig. 11 as the magenta line. This simple argument predicts a
bimodal distribution of ‘failed’ cores at the initial mass of
$0.01\text{\,}\mathrm{M_{\oplus}}$ and at a larger mass varying from
$1\text{\,}\mathrm{M_{\oplus}}10\text{\,}\mathrm{M_{\oplus}}$. The actual
upper peak shown as the dash-dotted lines are for Class I and II located a
factor of $\sim 4$ above this simple estimate. From evolution tracks it
becomes evident that planets in this category can be trapped at the edges of
outward migration regions or undergo giant impacts which enhances their mass.

The fact that collisions between embryos take place is factored into the
Goldreich mass, which we however only applied to the inner rocky planetary
system, where the simulation time is long compared to the growth timescale.
The rocky planets in Class I systems are also bimodally distributed. There, we
see that the lower peak is located close to the isolation masses while the
upper is well reproduced by the Goldreich mass. This can be interpreted as
systems or planets which went into a dynamically active stage leading to a
giant impact phase forming the upper peak while lower mass planets are also
predicted by the simulations remaining at their location without excitation
and the possibility to grow from solid material outside their feeding zones.

For Class II systems with significant migration patterns, the low-mass rocky
planets remain slightly below the mean isolation masses while the upper peak
is reduced. Also, the overall number of initially rocky planets existing at
the end of the simulations is reduced, a trend which continues also for the
classes with giant planets.

This reduces the statistical significance and meaningfulness of the comparison
for the rocky planets. Furthermore, for Classes III and IV, the final
planetary systems are dominated by the massive, gas-rich planets which leads
to expected deviations from the simple, solid accretion based mass scales. For
gas-dominated planets, see Adams and Batygin, (2022).

###  5.3 Links between initial conditions, system properties, and architecture
classes

Now, we aim to understand the conditions needed for the formation of the
different kinds of planetary systems. We find that among the different initial
conditions of our simulations, the initial mass in solids $M_{\rm p,ini}$ is
the one that can be used best to discriminate the different classes of
systems. However, it is not the only discriminant.

One way of illustrating the impact of the initial mass in solids is shown in
Fig. 12. It shows for all 1000 synthetic planetary system as a function of the
initial mass of solids (planetesimals) in the natal disc $M_{\rm p,ini}$ the
final mass of the planetary system $M_{\rm sys}$ (i.e., sum of all planets in
a system at $5\text{\,}\mathrm{Gyr}$).

Figure 12: Final mass of the planetary system (i.e., sum of all planets in a
system at $5\text{\,}\mathrm{Gyr}$) $M_{\rm sys}$ as a function of the initial
mass of the solids (planetesimals) in the parent disc $M_{\rm p,ini}$. The
colours show the number of planets in each system with a mass of at least
$1\text{\,}\mathrm{M_{\oplus}}$. From this colour code, we see how the number
of planets as a proxy of system architecture systematically changes with the
initial mass in solids. The vertical shaded region ‘MMSN’ shows the estimated
initial minimum solid content of the Solar Nebula. The horizontal dotted line
‘Solar System’ is the total system mass today ($\approx
450\text{\,}\mathrm{M_{\oplus}}$ = essentially the sum of the masses of the
four giant planets). Diagonal lines show conversion efficiencies (which can be
larger than unity because accreted H/He is included in the system masses, but
not in $M_{\rm p,ini}$).

For this system mass, one can clearly distinguish the systems without ($M_{\rm
sys}\lesssim$150\text{\,}\mathrm{M_{\oplus}}$$) and with one or several giant
planets ($M_{\rm sys}\gtrsim$150\text{\,}\mathrm{M_{\oplus}}$$). As expected
for a formation model based on the core accretion paradigm (e.g. Mordasini et
al., 2012a, ), giant planets only form in discs with sufficiently high $M_{\rm
p,ini}$ of about 150 to $200\text{\,}\mathrm{M_{\oplus}}$. For the systems
without giants, there is a clear linear scaling between $M_{\rm p,ini}$ and
$M_{\rm sys}$ with an efficiency factor of about $30\text{\,}\mathrm{\char
37\relax}$. For the systems with giants, there is still a correlation, but it
is less clear. This is expected, because for the final masses with giant
planets the disc gas mass is important and not the solid mass. There are also
some special cases with very small $M_{\rm sys}$ for high $M_{\rm p,ini}$.
They are discussed below.

The colours in Fig. 12 show the number of planets with a mass of at least
$1\text{\,}\mathrm{M_{\oplus}}$ in each system (including planets at all semi-
major axes). This planetary multiplicity is a proxy for system architecture
and exhibits an interesting pattern. At the lowest $M_{\rm p,ini}$ of about
$10\text{\,}\mathrm{M_{\oplus}}$, not even a single planet with a mass of at
least $1\text{\,}\mathrm{M_{\oplus}}$ forms, because the planets remain at
sub-Earth masses. As we move to higher $M_{\rm p,ini}$, but still below the
value needed for giant planet formation, the number of planets increases. At
$M_{\rm p,ini}$ between 100 and $200\text{\,}\mathrm{M_{\oplus}}$ there are
some planetary systems with a very high multiplicity of up to 25. In this
region, we see also a vertical stratification where at given $M_{\rm p,ini}$,
systems with a lower $M_{\rm sys}$ also have a lower number of planets. This
is a consequence of another disc property, the disc lifetime, which is
influenced mostly by the initial mass of the gas disc and the external
photoevaporation rate. At given $M_{\rm p,ini}$, system with a longer disc
lifetime have a lower $M_{\rm sys}$. The reason is not that these systems form
less planets (smaller mass or lower number) in the first place, but that
planets have more time to migrate. This has the consequence that more planets
fall into the star both by migration in resonant convoys during the presence
of the gas disc, or on long timescales after disc dissipation via tides. For
the latter the close-in parking location in systems with significant inward
migration makes the difference. This reduces both $M_{\rm sys}$ and the number
of planets. It also influences the classification of the architecture, as we
will see soon.

Moving finally to the $M_{\rm p,ini}$ that are high enough to allow giant
planet formation, we see that the number of planets again decreases. This is a
consequence of the fact that growing giant planets can destabilise lower mass
planets in their vicinity. These planets are then either ejected from the
system or accreted by the growing giant planet(s). This can be seen in Fig. 7
where the star symbols on the tracks of the growing giants indicate the
accretion of lower mass planets during the process of rapid gas accretion. At
$M_{\rm p,ini}$ less than about $400\text{\,}\mathrm{M_{\oplus}}$, there are
however still many systems with about 10 planets. Violent instabilities among
the giants are particularly efficient in removing lower-mass planets. Such
cases can be seen at the highest $M_{\rm p,ini}$. For these systems, the
number of planets reaches—as for the lowest $M_{\rm p,ini}$—again very low
values of just one or two (but now very massive) planets.

####  5.3.1 Solar System analogues

Figure 12 also shows the location of the minimum-mass Solar nebula (MMSN; with
the boundaries discussed in Emsenhuber et al., 2021b, : $M_{\rm p,ini}$ of 66
to $151\text{\,}\mathrm{M_{\oplus}}$) and the final Solar System. The Solar
System is not straightforward to form in our current model. This is mainly
because gas accretion is assumed to remain unaffected by gap formation.
Planets undergoing runaway gas accretion therefore will rapidly attain a
Jupiter mass or more, leading to systems that are more massive than the Solar
System. Obtaining more giant planets whose mass is compatible with that of
Jupiter or Saturn would require lower mass accretion rates (Suzuki et al.,,
2018; Nayakshin et al.,, 2019). Low disc viscosities (Suzuki et al.,, 2018,
and references therein) coupled with the inclusion of gap formation could
alleviate this problem. Further, the model does not produce many systems with
giant planets for discs whose mass are in the MMSN range. Rather, $M_{\rm
p,ini}$ that are a factor $\sim 2$ more massive are need. This is linked to
the overall conversion efficiency of the model, where only about a third of
the solids are accreted onto the planets, as it can be seen from the many
systems without giant planets. In the outer part of the disc, accretion
timescale is larger than the disc lifetime (see the pink curve in Figs. 3, 5,
and 7) so only a fraction of the mass is converted into planets. On the other
hand, the MMSN is a minimum mass indeed, so the actual $M_{\rm p,ini}$ of the
Solar System might had been higher.

In addition, we find that that planets in Solar-system analogues are too
close-in (Fig. 8). This is in part due to the difficultly of forming distant
planets with planetesimal accretion (e.g. Rafikov,, 2011). There are several
others aspects in the model which are missing for the specific case of the
Solar System. One is that the model does not include the Masset and
Snellgrove, (2001) mechanism that leads to outward migration of a pair of
giant planets trapped in a mean-motion resonance and sharing a common gap,
which has been invoked for the formation of Jupiter and Saturn (e.g. Walsh et
al.,, 2011). Further, we do not take into account that by blocking or at least
reducing the pebble flux, giant planet formation can affect the formation of
planetesimals, which later form the terrestrial planets (Morbidelli et al.,,
2015).

####  5.3.2 Architecture class

Figure 13: Histograms of the fractions of the different system class as
function of initial mass in solids (left), metallicity (centre), and initial
mass in gas (right). Blue represents Class I, orange represents Class II,
green represents Class III, and red represents Class IV. The mass in solids is
the quantity that most clearly discriminates the classes.

Figure 14: System architecture as function of the initial solids mass content
in the disc $M_{\rm p,ini}$ and the total mass (accounting for solids and gas)
in the final planets per system $M_{\rm sys}$ (left) and a histogram of
initial solids mass only (right). In the left panel, the three dashed diagonal
lines again denote the conversion efficiency from disc solids to planets, with
the value given next to them. In the right panel, the small horizontal lines
show the extent of each bin.

Instead of considering the number of individual planets in a system as a proxy
of system architecture, we can also study the four classes directly. We
illustrate the effect of different initial disc properties on architecture in
Fig. 13, with histograms of the fraction of system class as function of solid
mass, metallicity, and gas mass. We find the mass in solids is the quantity
that best discriminates the different categories, having the different classes
best confined to certain regions. The effect of the solid disc mass is
reflected on the metallicity and gas mass, since the three quantities are
related (only two of the three can be freely chosen). The middle panel shows
that Class I systems are anti-correlated with metallicity, having the largest
fraction of systems at low metallicity and steadily decreasing, while the
other three quantity show an initial increase following by a plateau. The gas
mass follows a similar pattern as the solid mass, but with more overlap
between the different classes.

In Fig. 14 we focus again on the solid mass, by again showing the total mass
in planets $M_{\rm sys}$ in each system is plotted against the initial content
of the solids disc $M_{\rm p,ini}$, but the systems are now colour-coded by
the four architecture classes. Here, we note a strong separation between the
first two (only low-mass planets) and last two (with giant planets) classes of
systems.

For the first two categories of systems, most of the planet masses is
comprised of solids. Therefore, the masses in planets is generally a given
fraction of the disc mass, between $20\text{\,}\mathrm{\char 37\relax}$ and
$35\text{\,}\mathrm{\char 37\relax}$, as mentioned. Usually, the whole inner
disc up to $1$ to $2\text{\,}\mathrm{au}$ is accreted by the protoplanets,
some fraction beyond that up to $\sim 10\text{\,}\mathrm{au}$, while the outer
part remain essentially intact, because the accretion timescales are too long.
It should be taken into account that extending the (numerical) integration
time (here $100\text{\,}\mathrm{Myr}$) in the simulations could further
increase the efficiency for some systems, until ejection of planetesimals will
start to limit further planet growth (Ida and Lin, 2004a, ). We also note that
there is a trend between Class I and II. The lower-mass discs usually form
Class I architectures, while more massive discs tend to result in Class II
systems.

As discussed, discs with solids content $M_{\rm p,ini}$ above $150$ to
$200\text{\,}\mathrm{M_{\oplus}}$ lead to the formation of giant planets. In
this case, planets accrete significant envelopes, which strongly increases the
conversion efficiency. However, we do not find a clear separation in the
parameter space between Class III and IV architectures, which is interesting.
We also searched the other initial conditions, including dust-to-gas ratio and
photoevaporation rate, but did not find anything clearer than what is
presented in Fig. 14. The right panel shows that Class IV systems have larger
disc masses than Class III systems, with median values of nearly
$440\text{\,}\mathrm{M_{\oplus}}$ and $563\text{\,}\mathrm{M_{\oplus}}$,
respectively, although there is significant overlap. This leaves us to
conclude that the distinction between the two classes is not only related to
initial disc characteristics and that the chaotic nature of the dynamical
interactions between the protoplanets also bears some responsibility. So for
similar disc initial conditions, Class III or IV systems can arise, depending
on the exact locations and timing of the moments the proto-giant planets in
one system start runaway accretion (which depends for example on the exact
initial positions of the embryos), and the associated presence or absence of
strong gravitational interactions. This (quasi) chaotic nature is different
than in a simpler model where only a single embryo per disc is present, whose
outcome is easy to predict (Schlecker et al., 2021b, ). It should be noted
that our full model used here is also deterministic, i.e., for given model
initial conditions it always predicts exactly the same outcome, but small
changes of the initial conditions may lead to very different outcomes because
of N-body interactions.

Finally, examples of dynamically-active Class IV systems leading to the loss
of most of the planets are visible in the bottom right corner of Fig. 14’s
left panel. Here close encounters lead to the loss of all giant planets and
only lower-mass planets remain. This leads to the (rare) points with high
$M_{\rm p,ini}$ but (very) low $M_{\rm sys}$. If we were to include in the
system mass also all planets that were ejected out into interstellar space or
that collided with the host star (and not only the planets that actually exist
at $5\text{\,}\mathrm{Gyr}$), these points would join the other ‘normal’
points at high $M_{\rm sys}$.

The classification by Mishra et al., 2023a ; Mishra et al., 2023b  was applied
to the same synthetic data as our approach. It is a complementary approach
since, in contrast to our human-based approach (that stresses underlying
physical processes and the emergence of planetary system architecture in
time), it can be applied to observational data without knowledge of the
formation history of the systems or the underlying formation processes. Visual
comparison between the left panel of Fig. 14 and the bottom right panel of
Fig. 2 of Mishra et al., 2023b  reveals the nevertheless existing
correspondence between our classification and theirs. We note that, broadly
speaking, the ‘similar’ class of Mishra et al., 2023b  corresponds to our
Class I and II. Also, our Class III systems correspond to ‘mixed’ systems and
a part of the ‘ordered’ ones, while Class IV corresponds to the other part of
‘ordered’ systems and the ‘anti-ordered’ ones. The approximate match of the
automated approach and the human classification for our multi-faceted data set
give a promising outlook for a future holistic categorisation for observed
systems. This should take into account planetary composition information as an
additional formation tracer, both regarding the bulk and the atmospheric
composition (Mollière et al.,, 2022).

###  5.4 Effect of the stellar environment

Figure 15: System architecture as function of the disc lifetime and external
photoevaporation rate as proxies of the stellar environment. Left panel:
architectural class as as a function of initial solids content of the disc
$M_{\rm p,ini}$ and the gas disc’s lifetime (which depends via external
photoevaporation on the stellar environment). Right: total mass in planets per
system $M_{\rm sys}$ as function of the external photoevaporation rate. In the
left panel, the dashed black line denotes roughly the separation between the
conditions leading to Class I and II systems. In the right panel, the dashed
black line shows the anticorrelation between the mass in planets and the
external photoevaporation rate in Class III and IV systems.

The surrounding stellar environment drives the ambient field, which in turn
drives the external photoevaporation of the protoplanetary disc. External
photoevaporation will truncate the outer disc and reduce their lifetimes
(Haworth et al.,, 2018). This in turn will affect planet formation, for
instance by limiting the formation time or reducing the migration (Winter et
al.,, 2022). Other effects of the environment include the heating of the disc
(Ndugu et al.,, 2018, e.g.) and stellar encounters (Stock et al.,, 2020), but
we do not include these effects in our model. Here, we focus on the effect of
the external photoevaporation rate and discuss two cases where we find a
correlation between the surrounding environment and either the class of system
or the planet masses.

The left panel of Fig. 15 shows that the transition between Class I and II
systems depends not only on the solids mass, but also the lifetime of the
protoplanetary disc as illustrated by the dashed black line. The longer the
disc lives, the less massive it needs to be for systems to be in Class II.
There is the following reason for this: Class II systems occur when the
planets get sufficiently massive to migrate substantially when the gas disc
still has significant mass. Thus, the planets have to reach a mass of about
$10\text{\,}\mathrm{M_{\oplus}}$ (i.e. to reach the saturation and/or the
equality mass) well before the dispersal of the gas disc. Now, the more
massive the solid disc, the faster the planets grow. This means that short-
lived discs needs to grow the massive icy planets fast, otherwise they will
not have the time to migrate. In other words: at fixed $M_{\rm p,ini}$ in the
range of about 100 to $300\text{\,}\mathrm{M_{\oplus}}$, short disc lifetimes
lead to Class I systems, while long ones lead to Class II, because there needs
to be sufficient time for orbital migration to occur to form a Class II
system. A prime example of this effect can be seen by comparing the system in
bottom right panel of Fig. 4 (Class I) and the one on the top left panel of
Fig. 6 (Class II). Both systems have a similar initial solid disc mass (120
versus $130\text{\,}\mathrm{M_{\oplus}}$), but different gas disc masses and
different disc lifetimes ($3.5$ versus $7.4\text{\,}\mathrm{Myr}$). This
results in terrestrial planets being able to survive in the former (although
icy planets migrated to within $1\text{\,}\mathrm{au}$) while they no longer
exist in the latter. The link to the environment comes from the
anticorrelation of disc lifetimes with external photoevaporation rate (Winter
and Haworth,, 2022). Absent this, disc lifetimes are correlated to their
initial mass. Thus, our model predicts that ordered and compositionally well-
separated planets are more likely in dense clusters.

The second dependency on the environment we find is between the mass of
$M_{\rm sys}$ for Class III and IV system, i.e., the mass of the giant planets
and the external photoevaporation rate, as shown in the right panel of Fig.
15. Here again, the reason is connected with the formation time of the
planets. The external photoevaporation rate strongly affect the disc
lifetimes, thus the time planets can accrete from the disc. The growth time of
a massive core is related to the mass of the solids disc, thus relatively
independent of the external photoevaporation rate. The cores will undergo
runaway gas accretion at roughly the same times. However, external
photoevaporation competes with planetary gas accretion. A similar correlation
in observation was also obtained (Winter et al.,, 2020).

###  5.5 Total number of giant planets versus surviving and loss channels

Figure 16: Correlation between the total and final number of giant planets in
each systems for those of Class III (mixed systems; left) and Class IV (giants
only; right). In each cell, the number of the top gives the percentage of
systems with the same end multiplicity originate from systems that formed that
number of giants (the sum of the top percentage on a row is unity, save for
rounding errors) while the bottom number give the percentage of systems who
formed the given number of giant planets ended with that multiplicity (the sum
of the bottom percentage on a column is unity, save for rounding errors). The
colour represent the underlying number of systems in each cell, on a
logarithmic scale.

As we saw, the formation of systems with (several) giant planets is often not
a smooth process but may involve strong dynamical events: the final system
architecture thus may not reflect the full picture of what happened during the
formation stage. In addition, we did not find any clear separation in the
initial conditions between systems in Classes III and IV. So, we want to check
if the formation history can help us to understand the final architecture. To
this effect, we select one simple metric: the total number of giant planets
that have been formed, no matter what they became (remained in the system till
the end; accreted by another more massive giant; ejected into outer space;
collided with the host star). We show the relations of the this total number
versus the number of giants that remain till the end in Fig. 16 for two
classes of systems: the mixed systems (Class III) and the ones with
essentially only giants (Class IV). For completeness, we note that one system
in Class II (migrated sub-Neptunes) has one giant planet at one point that was
then accreted by the central star; otherwise none of these have or had giants.

The results show that most of the mixed systems (Class III – in total 123
systems) did not loose any giant planets, as one might expect. The two main
situations are systems with two giants formed and two remaining (44 systems),
then one giant formed and one remaining (31 systems). Accounting for the
systems that formed and retained three or zero giants, there are in total 86
systems that did not loose any giant, compared to 37 that did loose at least
one. The largest loss channel for these giants is giant impacts (about
$53\text{\,}\mathrm{\char 37\relax}$) followed by mass loss from the envelope
(meaning planets that are still present, but were classified as giant at one
point during the formation and no longer are at the end, for example due to
the ejection of the envelope after a giant impact; about
$33\text{\,}\mathrm{\char 37\relax}$). Only about $14\text{\,}\mathrm{\char
37\relax}$ of the giants that are lost are due to ejections or collisions with
the central star.

For the systems with only giant planets (Class IV – in total 80 systems), the
story is quite different. Here, 23 systems did not loose any giant planets
while 57 did loose at least one. Further, the distribution of giants lost
across the different channels is different from Class III. About
$82\text{\,}\mathrm{\char 37\relax}$ of the losses are the result of ejections
or collisions with the central star.

Thus, in total, nearly half of the systems that formed giant planets did loose
at least one giant (94 out of 203). However, these are disproportionally found
in systems that are left with giant planets only (Class IV). On top of that,
the loss mechanisms are different between the classes. The loss of giant
planet by ejection or collision with the central star is usually sufficient to
destabilise the inner system, so that the resulting planetary system ends up
in Class IV. Conversely, the loss of a giant planet by collision (a giant
impact between giant planets) is in itself usually not sufficient to
destabilise the system. Planetary systems who did not loose a giant planet or
did so by collisions are more likely to be mixed systems (Class III).

Nevertheless, the ejection of a giant planet is not a complete indicator of
dynamical activity. We may give as example system 4 shown in Fig. 10, which is
classified as giants-only (Class IV). It formed two giants yet did not loose
either of them. Here, we have a pathway that the two giant planets interacted
dynamically and that alone led to collisions and ejections in the inner part
of the system.

###  5.6 Effect of stellar metallicity

####  5.6.1 Planet eccentricities

Figure 17: Eccentricity distribution of the giant planets (mass
$>100\text{\,}\mathrm{M_{\oplus}}$). Left: eccentricity versus initial
metallicity. The black line shows the mean value in each interval. Right:
cumulative distribution of giant planets eccentricities in Class III and Class
IV systems.

There are several correlations between host star metallicity and exoplanet
characteristics, as discussed in Sect. 2.4. One is a positive trend between
stellar metallicity and giant planet eccentricities (Dawson and Murray-Clay,,
2013; Buchhave et al.,, 2018). To check for this effect in our model, we
provide in the left panels of Fig. 17 the eccentricities of giant planets as
function of metallicity. We first see that overall giant planet eccentricities
are rather low in our model, with a mean value of $0.17$. Still, we get a
slight correlation between stellar metallicity and the mean eccentricity
(represented by the black line).

In addition, we see a different trend in the eccentricities as function of the
system class. The highly-eccentric planets ($e\gtrsim 0.6$) are mostly found
in Class IV systems and the mid-eccentricities ($0.6\gtrsim e\gtrsim 0.3$) are
also dominated by Class IV systems, though a higher proportion of planets in
Class III systems is also present. Only in the low eccentricity regime
($0.3\gtrsim e$) do planets in Class III systems dominate. To support this, we
provide in the right panel of Fig. 17 distinct cumulative functions for the
giant planet eccentricities as function of class. There, we see that the
eccentricity distribution is different between the two classes, with planets
in Class IV systems having overall larger eccentricities.

We also have a different distribution of stellar metallicities between Class
III and Class IV systems. We already saw in Fig. 14 that Class IV systems
tended to have larger initial disc masses and this is also reflected in the
metallicity distribution. The median metallicity of Class III systems is
$0.11$ while for Class IV systems it is $0.17$. Thus, our model can reproduce
the general trend of increasing planet eccentricities with stellar
metallicities, albeit eccentricities are overall low.

####  5.6.2 Planet and system kinds

Figure 18: Metallicity dependence of certain planet and system kinds in the
synthetic population, which can be compared to observations. Left: occurrence
of planet classes at periods between $10$ and $100\text{\,}\mathrm{d}$
(including an observability criterion similar to Petigura et al.,, 2018) as
function of stellar metallicity (transit detections). Super-Earths are between
$1.0$ and $1.7\text{\,}\mathrm{R_{\oplus}}$, mini-Neptunes between $1.7$ and
$4.0\text{\,}\mathrm{R_{\oplus}}$, sub-Saturns between $4.0$ and
$8.0\text{\,}\mathrm{R_{\oplus}}$, and Jupiters above
$8.0\text{\,}\mathrm{R_{\oplus}}$. Right: cumulative distribution of
metallicities of systems harbouring an inner Super-Earth (SE; defined as
distance between $0.02$ and $1\text{\,}\mathrm{au}$ and a mass between $2$ and
$30\text{\,}\mathrm{M_{\oplus}}$) and with (without) outer giant planet – or
cold Jupiter (CJ) – in green (blue), defined as distance between $0.23$ and
$10\text{\,}\mathrm{au}$ and mass between $30$ and
$6000\text{\,}\mathrm{M_{\oplus}}$. Only planets with a radial velocity semi-
amplitude $K>$2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ are accounted
for.

There are observational correlations between kinds of either planets or system
architectures and stellar metallicity, both in transit and radial velocity
surveys. To assess our model against such trends, we perform two comparisons.
First, we compare the metallicity effect of different kind of planets
following Petigura et al., (2018). Here we select planets whose orbital period
is between $10$ and $100\text{\,}\mathrm{d}$, apply an selection similar to
transit surveys, as described in Emsenhuber et al., 2021b , and categories the
planets according to their radii following the definitions of Petigura et
al.,, 2018: super Earths from $1.0$ to $1.7\text{\,}\mathrm{R_{\oplus}}$,
mini-Neptunes from $1.7$ to $4.0\text{\,}\mathrm{R_{\oplus}}$, sub-Saturns
from $4.0$ to $8.0\text{\,}\mathrm{R_{\oplus}}$, and Jupiters above
$8.0\text{\,}\mathrm{R_{\oplus}}$. The occurrence of small-radii planets,
‘super Earths’, show an anti-correlation with stellar metallicity. For
comparison, Petigura et al., (2018) found that the occurrence of the same
planets (panel b of their Fig. 10) has overall no significant correlation with
metallicity. However, the occurrence rate of the largest metallicity bin
(0.2–0.4) is lower than the others, as in our simulations. Our finding is also
consistent with the occurrence of low-mass planets reported by Emsenhuber et
al., 2021b , where systems having Earth-like (mass from $0.5$ to
$2\text{\,}\mathrm{M_{\oplus}}$) planets within $1\text{\,}\mathrm{au}$ have a
median metallicity of $-0.09$ (the same value for the population presented
here is $-0.11$). All other categories, from mini-Neptunes to Jupiters show a
positive correlation between occurrence rate and metallicity. For the mini-
Neptunes and sub-Saturns this is consistent with the trend found by Petigura
et al., (2018), while the authors have been unable to find any trend for
Jupiters from the Kepler mission, as there were only few such planets
observed. Our model is thus able to reproduce the general trends correlation
with metallicity found in the Kepler survey.

87% of super Earths are found in Class I systems, whose occurrence rate is
anticorrelated with metallicities (Fig. 13, centre panel). This, combined with
their total absence from Class IV systems, explains the trend of anti-
correlation between the occurrence of super Earths with stellar metallicity.
Mini-Neptune planets originate principally from Class II (49%) and Class I
(40%) systems. As the occurrence of Class II systems is positively correlated
with stellar metallicity, this explains the shift in the metallicity
correlation between super Earths and mini-Neptunes. From this discussion, we
also see that Class I systems may produce large planets (both in mass and
radii) which enter the mini-Neptune category. Having them originate more from
Class II systems would be a way to reinforce the metallicity correlation.

A second comparison can be done on the correlation between inner super Earths
and distant giant planets (or cold Jupiters). Such a correlation was observed
(Zhu and Wu,, 2018; Bryan et al.,, 2019) and the metallicity effect of these
systems was studied by Rosenthal et al., (2022). Schlecker et al., 2021a
determined there was a correlation between the occurrence of super Earths and
cold Jupiters in the same model as presented in this work (although assuming
$20\text{\,}\mathrm{Myr}$ for the formation stage), though weaker than in
observations. Here we extend the comparison for the trend in metallicity
between systems hosting only inner super Earths and both super Earths and cold
Jupiters in radial velocity survey. For this, we first restrict the synthetic
planets to those with a radial velocity semi-amplitude
$K>$2\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$, as in Schlecker et al.,
2021a . We use the same definitions as Rosenthal et al., (2022), name super
Earths have distances between $0.02$ and $1\text{\,}\mathrm{au}$ and masses
between $2$ and $30\text{\,}\mathrm{M_{\oplus}}$ while cold Jupiters have
distances between $0.23$ and $10\text{\,}\mathrm{au}$ and masses between $30$
and $6000\text{\,}\mathrm{M_{\oplus}}$. The cumulative distributions of
metallicity for system hosting only super Earths and both super Earths and
cold Jupiters are shown with the blue and green curves in the right panel of
Fig. 18, respectively. The p-value of a Kolmogorov-Smirnov (KS) test between
the two distributions to assess whether they are drawn from the same sample
yields $3.7\text{\times}{10}^{-3}$, implying that the two distributions are
distinct. Our results are compatible with those Rosenthal et al., (2022) in
that 1) systems with only super Earths have lower metallicities than those
with both super Earths and cold Jupiters, 2) a statistical test assesses that
the two groups have a distinct underlying distribution, and 3) the metallicity
distributions are broadly similar. In Rosenthal et al., (2022), the
metallicity of systems with only super Earths spans roughly $-0.4$ to $0.3$
with most of the systems between $-0.2$ and $0.2$, while the metallicity of
systems with both super Earths and cold Jupiters spans roughly $0.05$ to
$0.35$ (their Fig. 8). In our model, the values are similar although the the
metallicities of the latter type are shifted towards lower values, with values
mainly between $-0.1$ and $0.3$. We find that $77\text{\,}\mathrm{\char
37\relax}$ of the systems with both super Earths and cold Jupiters have
$[\mathrm{Fe/H}]>0$, and $41\text{\,}\mathrm{\char 37\relax}$ of them have
$[\mathrm{Fe/H}]>0.15$. Thus, while the systems containing both kinds of
planets are not as strongly biased towards large metallicities in our
synthetic population compared to what Rosenthal et al., (2022) obtained, we
nonetheless obtain a bias.

When comparing these types of systems with our classes, we find that systems
with only detectable super Earths are $77\text{\,}\mathrm{\char 37\relax}$
from Class II, $20\text{\,}\mathrm{\char 37\relax}$ from Class I and
$3\text{\,}\mathrm{\char 37\relax}$ from Class III. Systems with both super
Earths and cold Jupiters are $94\text{\,}\mathrm{\char 37\relax}$ from Class
III, with contributions from all others categories. We note that even a system
of Class I is included here; this is because that system contains a rocky
planet of $18\text{\,}\mathrm{M_{\oplus}}$ at $0.15\text{\,}\mathrm{au}$ and
an icy planet of $36\text{\,}\mathrm{M_{\oplus}}$ at
$0.39\text{\,}\mathrm{au}$. Although the planet are massive for their kind
(and hence enter the definitions of super Earth and cold Jupiters of Rosenthal
et al.,, 2022), the formation pattern of the system is similar to that of
Class I with only a small effect of migration and was thus categorised as
such. Thus, the difference in metallicities between super Earth systems and
with or without cold Jupiters is consistent with Class III systems having
overall larger metallicities than Class II systems.

We find different metallicity effects for planets defined as ‘super Earths’ in
transit or radial velocity surveys. The super Earths according to the
definition of Petigura et al., (2018) show an anticorrelation, or a negative
metallicity effect. Conversely, using the definitions of the radial velocity
survey of Rosenthal et al., (2022), the median metallicity of systems hosting
super Earths is $0.05$, thus showing a positive metallicity effect. The
difference is due to $75\text{\,}\mathrm{\char 37\relax}$ of the defined as
super Earths by Rosenthal et al., (2022) having radii between $1.7$ and
$4.0\text{\,}\mathrm{R_{\oplus}}$, the range defined as ‘mini-Neptunes’ by
Petigura et al., (2018). We thus caution that there is an inconsistency
between the naming of ‘super Earths’ in transit and radial velocity, and that
‘super Earths’ in radial velocity are similar to ‘mini-Neptunes’ in transit
surveys, which both show a positive metallicity effect.

We further find that the median metallicity of systems where a giant planet
has been accreted by the central star is $0.15$. This is an indication that
our model could reproduce the metallicity effect of the high-eccentricity
migration channel of hot-Jupiters (Dawson and Murray-Clay,, 2013; Dawson and
Johnson,, 2018). However, we cannot fully confirm this because the model does
not account for tidal circularisation.

##  6 Discussion and conclusion

In this work, we discussed how planetary population synthesis can be used to
infer how planetary formation may work although we have at the moment only few
direct observations of forming planets. The general idea is to find a
theoretical end-to-end model that can reproduce the observed planet population
from the diversity of known protoplanetary discs in a quantitative statistical
way. The underlying theoretical model has to be global, so that it includes
the essence of all (currently known) important governing physical processes
and treat their interactions. We can also use such a global model to
investigate the physical processes that lead to different classes of
architectures of synthetic planetary systems and how the disc and stellar
environmental properties affects the outcome. Hence, we classified the $1000$
synthetic systems from the nominal NGPPS population of Emsenhuber et al.,
2021b  based on classes of formation pathways seen in the simulations and
resulting planetary system architectures. This approach stresses the temporal
emergence of system architectures and the underlying physical processes
included in the theoretical model. This is in contrast to the approach of
Mishra et al., 2023a ; Mishra et al., 2023b , who classified the systems based
on the observable present day properties of the final systems.

We find four classes of formation and architecture: (compositionally) ordered
Earths and ice worlds systems (Class I), migrated sub-Neptunes systems (Class
II), mixed systems with inner low-mass and giant planets (Class III), and
systems of dynamically-active giants (Class IV). Then, three different types
of analysis of these classes are presented: first, using (analytical) mass
scales, namely the growth, equality, saturation, Goldreich, and the critical
core mass which can be used to understand which physical process is mainly
responsible for shaping the final systems architecture and the characteristic
mass scale in a system. Second we studied the disc initial conditions and the
stellar environment that lead to the formation of systems in the four classes.
Third, we assess how the model reproduces the observational trends with
stellar metallicity on planet eccentricities and the occurrence of certain
kinds of planet and architectures.

The compositionally ordered Earths and ice worlds systems, which consist of
low-mass rocky planets inside, and icy ones outside, form from the lowest-mass
discs mostly by in-situ accretion of planetesimals followed by giants impacts
among the protoplanets. The relevant mass scaling for the final planet masses
is the Goldreich mass, which applies to systems ending with a giant impact
phase. Intermediate-mass discs lead to migrated sub-Neptune systems, where
sub-Neptunian planets migrate all the way to the inner disc edge. The relevant
mass scaling is here the equality mass (where the accretion and migration
timescales are the same) and/or the saturation mass (where planets are
released from Type I migration traps because the positive corrotation torque
saturates). It is interesting to note that the disc solid mass required to
transition from Class I to Class II depends on the disc lifetime. This is
because more massive discs form sub-Neptunes more rapidly and then undergo
gas-driven migration towards the inner region.

When assessing the correlation between planet kind and metallicity, we find
that the observed trends are reproduced, although overall super Earths and
sub-Neptunes are found in lower-metallicity environment than in observations.
In our synthetic populations, a significant number of observable sub-Neptunes
originate from Class I systems. This tends to favour systems with lower
metallicities. Having these planets originate predominantly from Class II
systems would favour higher metallicities, providing a better correspondence
between the model and observations. Whether this is possible should be
investigated in the future.

The other two classes correspond to systems with giants planets that have
(Class III) or do not have (Class IV) inner lower-mass companions. Here, the
relevant formation pathway is that the equality mass must be of the order of
$10\text{\,}\mathrm{M_{\oplus}}$, which corresponds to the critical core mass,
to allow for strong gas accretion while the planets migrate. The distinction
between the two classes is whether some strong dynamical instability has taken
place between the giants, which lead to the loss of the inner companions.
While these two categories together are well separated from the others in the
space of initial conditions (Fig. 14), there is no such clear separation in
the initial disc properties (like initial mass of solids) between the systems
in Class III versus Class IV. Thus, we attribute the discerning factor mainly
to the stochasticity of the N-body interactions: for some initial placings of
the embryos, the forming giants do undergo dynamical instabilities, for others
not. There are nevertheless some trends, with Class IV systems forming
preferentially from higher-mass discs (Fig. 14) and where photoevaporation is
lower and thus disc lifetimes are longer (Fig. 15).

We further investigate whether the total number of giant planets that ever
existed in a system can explain the distinction between Classes III and IV. We
find that the loss of a giant planet due to an ejection or a collision with
the central star (as the result of a close encounters between two giant
planets) will likely result in the destabilisation of the inner system,
thereby resulting in a system with only giant planets. However, if the two
giant planets collide rather than scatter themselves, this is not sufficient
to destabilise the inner system. Thus, the stochastic nature of close
encounters plays an important role in the distinction between Class III and
IV, rather than only the initial conditions.

The discovery space of extrasolar planet detection methods expands
continuously. Several upcoming missions and instruments like GAIA, NIRPS,
PLATO, Roman, or Ariel will foster this trend even more, and allow to build
observationally a more complete statistical picture of the demographics of
extrasolar planets. The times are therefore promising for a method like
planetary population synthesis which can take full advantage of the wealth of
all these new observational constraints. It will, however, be important to
compare the same model with all different constraints and methods
simultaneously to uncover the shortcomings in our theoretical understanding,
and to avoid building models to reproduce only a certain observation.
Furthermore, it is time to conduct comparisons of models and observations in a
quantitative, and no longer a qualitative fashion only, given that the
observational data is fully on a level to allow for this. Combined with direct
observations of ongoing planet formation, this should make it possible to
markedly improve our understanding of the planet formation process in the
coming decade.

## Acknowledgements

We thank the anonymous reviewer for the helpful comments and suggestions. This
work was funded in part by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - 362051796. It has been carried out partially within the
framework of the NCCR PlanetS supported by the Swiss National Science
Foundation under grants 51NF40_182901 and 51NF40_205606. C.M. acknowledges the
funding from the Swiss National Science Foundation under grant 200021_204847
‘PlanetsInTime’. Most plots shown in this work were generated using matplotlib
(Hunter,, 2007).

## Data availability

The datasets analysed during the current study are available from the
corresponding author on reasonable request.

## References

  * Adachi et al., (1976) Adachi, I., Hayashi, C., and Nakazawa, K. (1976).  The gas drag effect on the elliptical motion of a solid body in the primordial solar nebula.  Progress of Theoretical Physics, 56:1756–1771.
  * Adams and Batygin, (2022) Adams, F. C. and Batygin, K. (2022).  Analytic Approach to the Late Stages of Giant Planet Formation.  , 934(2):111.
  * Adibekyan et al., (2014) Adibekyan, V. Z., González Hernández, J. I., Delgado Mena, E., Sousa, S. G., Santos, N. C., Israelian, G., Figueira, P., and Bertran de Lis, S. (2014).  On the origin of stars with and without planets. Tc trends and clues to Galactic evolution.  , 564:L15.
  * (4) Alessi, M., Inglis, J., and Pudritz, R. E. (2020a).  Formation of planetary populations - III. Core composition and atmospheric evaporation.  , 497(4):4814–4833.
  * Alessi and Pudritz, (2018) Alessi, M. and Pudritz, R. E. (2018).  Formation of planetary populations - I. Metallicity and envelope opacity effects.  , 478(2):2599–2617.
  * Alessi and Pudritz, (2022) Alessi, M. and Pudritz, R. E. (2022).  Combined effects of disc winds and turbulence-driven accretion on planet populations.  , 515(2):2548–2577.
  * (7) Alessi, M., Pudritz, R. E., and Cridland, A. J. (2020b).  Formation of planetary populations - II. Effects of initial disc size and radial dust drift.  , 493(1):1013–1033.
  * Alibert et al., (2013) Alibert, Y., Carron, F., Fortier, A., Pfyffer, S., Benz, W., Mordasini, C., and Swoboda, D. (2013).  Theoretical models of planetary system formation: mass vs. semi-major axis.  , 558:A109.
  * Alibert et al., (2004) Alibert, Y., Mordasini, C., and Benz, W. (2004).  Migration and giant planet formation.  , 417:L25–L28.
  * Alibert et al., (2011) Alibert, Y., Mordasini, C., and Benz, W. (2011).  Extrasolar planet population synthesis. III. Formation of planets around stars of different masses.  , 526:A63.
  * Alibert et al., (2005) Alibert, Y., Mordasini, C., Benz, W., and Winisdoerffer, C. (2005).  Models of giant planet formation with migration and disc evolution.  , 434(1):343–353.
  * Alibert and Venturini, (2019) Alibert, Y. and Venturini, J. (2019).  Using deep neural networks to compute the mass of forming planets.  , 626:A21.
  * Ansdell et al., (2016) Ansdell, M., Williams, J. P., van der Marel, N., Carpenter, J. M., Guidi, G., Hogerheijde, M., Mathews, G. S., Manara, C. F., Miotello, A., Natta, A., Oliveira, I., Tazzari, M., Testi, L., van Dishoeck, E. F., and van Terwisga, S. E. (2016).  ALMA Survey of Lupus Protoplanetary Disks. I. Dust and Gas Masses.  , 828(1):46.
  * Baraffe et al., (2015) Baraffe, I., Homeier, D., Allard, F., and Chabrier, G. (2015).  New evolutionary models for pre-main sequence and main sequence low-mass stars down to the hydrogen-burning limit.  , 577:A42.
  * Baraffe et al., (2004) Baraffe, I., Selsis, F., Chabrier, G., Barman, T. S., Allard, F., Hauschildt, P. H., and Lammer, H. (2004).  The effect of evaporation on the evolution of close-in giant planets.  , 419:L13–L16.
  * Baruteau et al., (2016) Baruteau, C., Bai, X., Mordasini, C., and Mollière, P. (2016).  Formation, Orbital and Internal Evolutions of Young Planetary Systems.  , 205(1-4):77–124.
  * Benítez-Llambay et al., (2011) Benítez-Llambay, P., Masset, F., and Beaugé, C. (2011).  The mass-period distribution of close-in exoplanets.  , 528:A2.
  * Bennett et al., (2021) Bennett, D. P., Ranc, C., and Fernandes, R. B. (2021).  No Sub-Saturn-mass Planet Desert in the CORALIE/HARPS Radial-velocity Sample.  , 162(6):243.
  * Benz et al., (2014) Benz, W., Ida, S., Alibert, Y., Lin, D., and Mordasini, C. (2014).  Planet Population Synthesis.  In Beuther, H., Klessen, R. S., Dullemond, C. P., and Henning, T., editors, Protostars and Planets VI, page 691.
  * Bergez-Casalou et al., (2022) Bergez-Casalou, C., Bitsch, B., Kurtovic, N. T., and Pinilla, P. (2022).  Constraining giant planet formation with synthetic ALMA images of the Solar System’s natal protoplanetary disk.  , 659:A6.
  * Bitsch et al., (2019) Bitsch, B., Izidoro, A., Johansen, A., Raymond, S. N., Morbidelli, A., Lambrechts, M., and Jacobson, S. A. (2019).  Formation of planetary systems by pebble accretion and migration: growth of gas giants.  , 623:A88.
  * Bodenheimer et al., (2000) Bodenheimer, P., Hubickyj, O., and Lissauer, J. J. (2000).  Models of the in Situ Formation of Detected Extrasolar Giant Planets.  , 143(1):2–14.
  * Bodenheimer and Pollack, (1986) Bodenheimer, P. and Pollack, J. B. (1986).  Calculations of the accretion and evolution of giant planets: The effects of solid cores.  , 67(3):391–408.
  * Bonfils et al., (2013) Bonfils, X., Delfosse, X., Udry, S., Forveille, T., Mayor, M., Perrier, C., Bouchy, F., Gillon, M., Lovis, C., Pepe, F., Queloz, D., Santos, N. C., Ségransan, D., and Bertaux, J. L. (2013).  The HARPS search for southern extra-solar planets. XXXI. The M-dwarf sample.  , 549:A109.
  * Borucki et al., (2010) Borucki, W. J., Koch, D., Basri, G., Batalha, N., Brown, T., Caldwell, D., Caldwell, J., Christensen-Dalsgaard, J., Cochran, W. D., DeVore, E., Dunham, E. W., Dupree, A. K., Gautier, T. N., Geary, J. C., Gilliland, R., Gould, A., Howell, S. B., Jenkins, J. M., Kondo, Y., Latham, D. W., Marcy, G. W., Meibom, S., Kjeldsen, H., Lissauer, J. J., Monet, D. G., Morrison, D., Sasselov, D., Tarter, J., Boss, A., Brownlee, D., Owen, T., Buzasi, D., Charbonneau, D., Doyle, L., Fortney, J., Ford, E. B., Holman, M. J., Seager, S., Steffen, J. H., Welsh, W. F., Rowe, J., Anderson, H., Buchhave, L., Ciardi, D., Walkowicz, L., Sherry, W., Horch, E., Isaacson, H., Everett, M. E., Fischer, D., Torres, G., Johnson, J. A., Endl, M., MacQueen, P., Bryson, S. T., Dotson, J., Haas, M., Kolodziejczak, J., Van Cleve, J., Chandrasekaran, H., Twicken, J. D., Quintana, E. V., Clarke, B. D., Allen, C., Li, J., Wu, H., Tenenbaum, P., Verner, E., Bruhweiler, F., Barnes, J., and Prsa, A. (2010).  Kepler Planet-Detection Mission: Introduction and First Results.  Science, 327(5968):977.
  * Brügger et al., (2020) Brügger, N., Burn, R., Coleman, G. A. L., Alibert, Y., and Benz, W. (2020).  Pebbles versus planetesimals. The outcomes of population synthesis models.  , 640:A21.
  * Bryan et al., (2016) Bryan, M. L., Knutson, H. A., Howard, A. W., Ngo, H., Batygin, K., Crepp, J. R., Fulton, B. J., Hinkley, S., Isaacson, H., Johnson, J. A., Marcy, G. W., and Wright, J. T. (2016).  Statistics of Long Period Gas Giant Planets in Known Planetary Systems.  , 821(2):89.
  * Bryan et al., (2019) Bryan, M. L., Knutson, H. A., Lee, E. J., Fulton, B. J., Batygin, K., Ngo, H., and Meshkat, T. (2019).  An Excess of Jupiter Analogs in Super-Earth Systems.  , 157(2):52.
  * Buchhave et al., (2018) Buchhave, L. A., Bitsch, B., Johansen, A., Latham, D. W., Bizzarro, M., Bieryla, A., and Kipping, D. M. (2018).  Jupiter Analogs Orbit Stars with an Average Metallicity Close to That of the Sun.  , 856(1):37.
  * Burn et al., (2021) Burn, R., Schlecker, M., Mordasini, C., Emsenhuber, A., Alibert, Y., Henning, T., Klahr, H., and Benz, W. (2021).  The New Generation Planetary Population Synthesis (NGPPS). IV. Planetary systems around low-mass stars.  , 656:A72.
  * Chambers, (2006) Chambers, J. (2006).  A semi-analytic model for oligarchic growth.  , 180(2):496–513.
  * Chambers, (2018) Chambers, J. (2018).  Planet Formation: An Optimized Population-synthesis Approach.  , 865(1):30.
  * Chambers, (1999) Chambers, J. E. (1999).  A hybrid symplectic integrator that permits close encounters between massive bodies.  , 304(4):793–799.
  * Clarke et al., (2001) Clarke, C. J., Gendrin, A., and Sotomayor, M. (2001).  The dispersal of circumstellar discs: the role of the ultraviolet switch.  , 328(2):485–491.
  * Coleman and Nelson, (2014) Coleman, G. A. L. and Nelson, R. P. (2014).  On the formation of planetary systems via oligarchic growth in thermally evolving viscous discs.  , 445(1):479–499.
  * Crida et al., (2006) Crida, A., Morbidelli, A., and Masset, F. (2006).  On the width and shape of gaps in protoplanetary disks.  , 181(2):587–604.
  * Cumming et al., (2008) Cumming, A., Butler, R. P., Marcy, G. W., Vogt, S. S., Wright, J. T., and Fischer, D. A. (2008).  The Keck Planet Search: Detectability and the Minimum Mass and Orbital Period Distribution of Extrasolar Planets.  , 120(867):531.
  * Currie et al., (2022) Currie, T., Lawson, K., Schneider, G., Lyra, W., Wisniewski, J., Grady, C., Guyon, O., Tamura, M., Kotani, T., Kawahara, H., Brandt, T., Uyama, T., Muto, T., Dong, R., Kudo, T., Hashimoto, J., Fukagawa, M., Wagner, K., Lozi, J., Chilcote, J., Tobin, T., Groff, T., Ward-Duong, K., Januszewski, W., Norris, B., Tuthill, P., van der Marel, N., Sitko, M., Deo, V., Vievard, S., Jovanovic, N., Martinache, F., and Skaf, N. (2022).  Images of embedded Jovian planet formation at a wide separation around AB Aurigae.  Nature Astronomy.
  * Dawson and Johnson, (2018) Dawson, R. I. and Johnson, J. A. (2018).  Origins of Hot Jupiters.  , 56:175–221.