11institutetext: School of Engineering, Physics, and Mathematics, University
of Dundee, DD1 4HN, UK
11email<EMAIL_ADDRESS>22institutetext: Earth-Life Science
Institute, Tokyo Institute of Technology, Meguro-ku, Tokyo, 152-8550, Japan

# N-body simulations of planet formation via pebble accretion II. How various
giant planets form

Soko Matsumura 11 Ramon Brasser 22 Shigeru Ida 22

###### Abstract

Aims. The connection between initial disc conditions and final orbital and
physical properties of planets is not well-understood. In this paper, we
numerically study the formation of planetary systems via pebble accretion and
investigate the effects of disc properties such as masses, dissipation
timescales, and metallicities on planet formation outcomes.

Methods. We improved the N-body code SyMBA that was modified for our paper I
by taking account of new planet-disc interaction models and type II migration.
We adopted the ‘two-$\alpha$’ disc model to mimic the effects of both the
standard disc turbulence and the mass accretion driven by the magnetic disc
wind.

Results. We successfully reproduced the overall distribution trends of semi-
major axes, eccentricities, and planetary masses of extrasolar giant planets.
There are two types of giant planet formation trends, depending on whether or
not the disc’s dissipation timescales are comparable to the planet formation
timescales. When planet formation happens fast enough, giant planets are fully
grown (Jupiter mass or higher) and are distributed widely across the disc. On
the other hand, when planet formation is limited by the disc’s dissipation,
discs generally form low-mass cold Jupiters (CJs). Our simulations also
naturally explain why hot Jupiters (HJs) tend to be alone and how the observed
eccentricity-metallicity trends arise. The low-metallicity discs tend to form
nearly circular and coplanar HJs in situ, because planet formation is slower
than high-metallicity discs, and thus protoplanetary cores migrate
significantly before gas accretion. The high-metallicity discs, on the other
hand, generate HJs in situ or via tidal circularisation of eccentric orbits.
Both pathways usually involve dynamical instabilities, and thus HJs tend to
have broader eccentricity and inclination distributions. When giant planets
with very wide orbits (’super-cold Jupiters’) are formed via pebble accretion
followed by scattering, we predict that they belong to metal-rich stars, have
eccentric orbits, and tend to have ($\sim 80\%$) companions interior to their
orbits.

Conclusions.

##  1 Introduction

Accretion of both pebbles and planetesimals is likely to contribute to planet
formation (e.g. Ida et al., 2016; Johansen & Lambrechts, 2017). The pebbles
are cm- to m-sized particles, which are very sensitive to gas drag (Adachi et
al., 1976; Weidenschilling, 1977). Since their migration timescale is short
compared to the growth timescale (Birnstiel et al., 2012; Lambrechts &
Johansen, 2014), they are expected to rapidly migrate towards the central star
without much growth. These migrating dust particles could be concentrated to
form $100-1000\,$km-sized planetesimals directly (which correspond to $\sim
10^{-6}-10^{-4}\,M_{\oplus}$) (e.g. Johansen et al., 2015; Simon et al.,
2016), either via the trapping of dust grains in pressure bumps of the
protoplanetary disc or via streaming instability (Youdin & Goodman, 2005;
Johansen et al., 2007, 2009). These planetesimals are initially likely to grow
via planetesimal-planetesimal collisions (Ida et al., 2016; Johansen &
Lambrechts, 2017), until the protoplanetary cores become massive enough to
have a large capturing radius of pebbles comparable to the Hill radius (Ormel
& Klahr, 2010; Kretke & Levison, 2014). The stage of a protoplanetary core
growth via pebbles is called ‘pebble accretion’. The pebble accretion slows
down when the capturing radius becomes larger than the pebble-disc scale
height, and thus the accretion becomes two-dimensional. At such a stage, both
pebble and planetesimal accretion may contribute to mass growth (e.g. Ida et
al., 2016).

Ormel & Klahr (2010) first developed the analytical model of particle-
protoplanet interactions in a gas disc, and pointed out the efficiency of
pebble accretion in the settling regime where gas drag effects are
significant. The pebble accretion was further studied by considering the
growth of a single core (Bitsch et al., 2015) and by using sophisticated
N-body simulations (Levison et al., 2015; Chambers, 2016). All of these
studies confirmed the efficiency of pebble accretion and showed that a variety
of planetary systems could be formed within a typical gas disc lifetime.
However, the model by Bitsch et al. (2015) was not designed to assess the
planet-planet interaction effects, which are important to determine the final
orbital architecture of planetary systems. The models by Levison et al. (2015)
and Chambers (2016) focused on the detailed physics of growth and destruction
of particles ranging from pebbles to planets, and did not take account of
planet migration effects, which are also important to determine the orbital
distribution of planets.

Matsumura et al. (2017, Paper I) implemented the analytical pebble accretion
model by Ida et al. (2016) into an N-body code SyMBA (Duncan et al., 1998),
and performed global N-body simulations of pebble accretion over a range of
disc parameters. Although the work has confirmed that pebble accretion leads
to a variety of planetary systems, it failed to reproduce the overall
distribution trends of the physical and orbital parameters of extrasolar giant
planets. In particular, few giant planets were formed, because most planetary
cores were lost to the disc’s inner edge due to migration, in a similar manner
to planetesimal accretion simulations (e.g. Coleman & Nelson, 2016b). The
problem was that a range of planetary masses leading to outward migration
(Paardekooper et al., 2011) is very narrow, so growing planets eventually
start migrating inward (also see Brasser et al., 2017). Therefore, although
pebble accretion may be more efficient compared to planetesimal accretion
(e.g. Lambrechts & Johansen, 2012; Kretke & Levison, 2014), the giant planet
formation timescale still appeared to be too long compared to the migration
timescale, and thus the migration problem remained.

There has been a significant development regarding the formation of ‘cold’
giant planets (i.e. Jupiter-like planets beyond $\sim 1\,$au) in the past few
years. The formation of such planets has been particularly difficult largely
because of the efficient migration of protoplanetary cores discussed above
(e.g. Coleman & Nelson, 2016b; Matsumura et al., 2017). Bitsch et al. (2015)
mitigated the issue and successfully formed cold Jupiters (CJs) by placing
cores far from the central star (beyond $\sim 15\,$au) in an evolved disc (up
to $\sim 2\,$Myr) with a relatively short disc lifetime of $3\,$Myr. Coleman &
Nelson (2016a) resolved the migration issue by considering the temporal planet
‘traps’ due to variations in the effective viscous stresses. More recently,
Ida et al. (2018) proposed an elegant solution by combining two recent key
developments in the field — the new type II migration formula (Kanagawa et
al., 2018) and the wind-driven accretion disc with the nearly laminar interior
(e.g. Bai & Stone, 2013; Bai, 2017), where both contribute to slowing down
type II migration and make the formation of CJs possible.

One of these key developments comes from the realisation that the migration of
a gap-opening planet is not tied to the disc evolution. This is because,
differently from the assumption of the classical type II migration model (Lin
& Papaloizou, 1986, 1993), the disc gas proved to cross the gap opened by a
planet. The new type II migration formula by Kanagawa et al. (2018) reflects
the results of recent hydrodynamic simulations (Duffell & MacFadyen, 2013;
Fung et al., 2014) as well as their own, and shows that the migration is
slower for a planet with the deeper gap (i.e. for a smaller viscosity disc, a
lower disc temperature, or a larger mass planet, also see Equation 8).
Compared to the classical type II formula, the new formula typically gives
$\sim 1$ order of magnitude longer migration timescales (Ida et al., 2018).

The other key development is the accretion mechanism of a protoplanetary disc.
In the classical picture, the disc’s viscosity drives the angular momentum
transfer in the protoplanetary disc, and the chief explanation for the
viscosity is the magnetorotational instability turbulence (e.g. Balbus &
Hawley, 1998). Recent non-ideal magnetohydrodynamic simulations, however, have
shown that the protoplanetary discs are likely to be largely laminar due to
the combined effects of ambipolar diffusion, Hall effects and Ohmic
dissipation, and the fact that the disc’s angular momentum is mostly removed
by the magnetic disc wind (e.g. Bai & Stone, 2013; Bai, 2017). Following this
idea, Ida et al. (2018) adopted two different alphas: one representing the
disc accretion driven by the mass loss due to the magnetic disc wind
$\alpha_{{\rm acc}}$, and the other representing the local disc turbulence
$\alpha_{{\rm turb}}$. Since the gap opening is controlled by the disc
turbulence and since the turbulence is expected to be very weak $\alpha_{{\rm
turb}}\lesssim 10^{-4}$ (e.g. Bai, 2017), the estimated type II migration is
slower than it is using the $\alpha$ required to explain the observed stellar
mass accretion rate ($\alpha\sim 10^{-2}$ in Hartmann et al. (1998), see
Section 2.2 as well). A similar study was conducted by Wimarsson et al. (2020)
using N-body simulations with pebble accretion. Besides the new type II
mechanism described above, they show that the transition zone between
viscously and radiatively heated disc regions could work as a temporal planet
trap and promote more efficient planetary growth via collisions.

Recently, a trilogy of papers on N-body pebble accretion simulations was
published, with two of them exploring the formation of super-Earths (SEs) or
lower mass planets (Lambrechts et al., 2019; Izidoro et al., 2019), and the
other focusing on the formation of giant planets (Bitsch et al., 2019). For
SEs and lower mass planets, they proposed that Earth-like planets were formed
in the low pebble flux disc with little migration, while SE-like planets were
formed in the higher pebble flux disc. In the latter case, their growing cores
migrated to the disc’s inner edge in a similar manner to Matsumura et al.
(2017), but these cores were more efficiently trapped in mean-motion
resonances (MMRs) because they chose the disc dissipation timescale comparable
to the migration timescale so that the migration would be halted completely
near the inner edge of the disc. As a result, their protoplanets became
dynamically unstable as the gas disc dissipated and the protoplanets collided
with one another and formed non-resonant multiple SEs (cf. Ogihara et al.,
2010). However, the disc-planet interactions near the disc edge is not well
understood, and the planet trapping efficiency depends on the sharpness of the
disc edge (Ogihara et al., 2010) and the details of the torque balance
(Brasser et al., 2018). For giant planets, on the other hand, Bitsch et al.
(2019) compared the classical $\alpha$ disc model with $\alpha=5.3\times
10^{-3}$ with a ‘two-$\alpha$’ disc model similar to that of Ida et al. (2018)
with $\alpha_{{\rm acc}}=5.3\times 10^{-3}$ and $\alpha_{{\rm turb}}=5.3\times
10^{-4}-10^{-4}$, and successfully produced CJs. However, they failed to form
highly eccentric giant planets (see discussion in Section 4.7).

The goal of this work, like that of our previous work (Matsumura et al.,
2017), is to investigate how initial disc conditions are related to physical
and orbital properties of planets by comparing planetary systems formed via
numerical simulations with observed extrasolar planetary systems. This kind of
study of connecting planetary systems with protoplanetary discs is becoming
increasingly important. ALMA observations now regularly find gaps, spirals,
and asymmetries in protoplanetary discs, which may be attributed to planet
formation processes or forming planets (e.g. Andrews et al., 2018; Huang et
al., 2018). There are also potential planetary candidates discovered in
protoplanetary discs including PDS 70 (Keppler et al., 2018), HD 100546
(Brittain et al., 2019; Casassus & Pérez, 2019), and TW Hya (Tsukagoshi et
al., 2019). The future observations would allow us to further constrain planet
formation processes occurring in protoplanetary discs.

As described in Section 2, we started with disc conditions motivated by
observations, follow planet formation and orbital evolution by using an N-body
code SyMBA (Duncan et al., 1998), and compare distributions of orbital and
physical properties of simulated planets with those of observed ones. From
this point of view, giant planets are particularly useful, because they are
observed over a wide range of orbital radii, from less than $0.01\,$au to over
$100\,$au. Super-Earths and lower-mass planets, on the other hand, are more
abundant but are limited to the inner disc region ($\lesssim 1\,$au).
Therefore, in this paper, we pay particular attention to disc conditions that
lead to different types of giant planets, ranging from close-in hot Jupiters
to very distant super-cold Jupiters, though numerical simulations will
generate all kinds of planets.

In Section 2, we introduce the numerical methods we used by highlighting the
changes made since Matsumura et al. (2017). We improve the N-body code SyMBA
(Duncan et al., 1998) by adopting a more self-consistent disc model and by
taking account of recent developments of the field. In Section 3, we show that
our numerical simulations reproduce overall trends of distributions of semi-
major axis $a$, eccentricity $e$, and mass $M_{p}$ of extrasolar giant
planets. We describe the effects of different parameters on planet formation
and also show how different types of giant planets (such as HJs and CJs) are
formed. We discuss the results further in Section 4, and summarise our
findings in Section 5.

##  2 Methods

We ran all of the simulations with the N-body integrator SyMBA (Duncan et al.,
1998), which was modified to include models of planet formation, disc-planet
interactions, and disc evolution. In this section, we highlight the updates
made since Matsumura et al. (2017).

In Section 2.1, we introduce the planet-disc interaction model adopted from
Ida et al. (2020),and we show how planet migration as well as eccentricity and
inclination damping of protoplanetary orbits are modelled in our code. We also
updated the type II migration prescription by following Kanagawa et al.
(2018). In Section 2.2, we introduce planet formation and disc evolution
models. We adopted the $\alpha$ disc model (Shakura & Sunyaev, 1973), but we
modified it to take account of a recent development in the field. Finally, in
Section 2.3, we present the initial conditions of our simulations, and show a
simple example of planet formation based on the prescriptions introduced in
this section.

###  2.1 Orbital evolution of protoplanets

The orbital evolution of protoplanets is largely determined by their
interactions with the gas disc. However, there are a few different definitions
of the equation of motion used in the planetary literature. The difference is
not only confusing but could lead to very different dynamical outcomes (e.g.
Ida et al., 2020; Matsumura et al., 2017). Recently, Ida et al. (2020)
compared a few such equations in the literature and highlighted which
equations are appropriate under which conditions. Based on this work, we will
use a new, more appropriate equation of motion.

In Matsumura et al. (2017), we adopted the equation proposed by Coleman &
Nelson (2014):

$\frac{d{\bf v}}{dt}=-\frac{{\bf v}}{\tau_{m}}-\frac{v_{r}}{\tau_{e}}{\bf
e}_{r}-\frac{0.5\left(v_{\theta}-v_{K,a}\right)}{\tau_{e}}{\bf
e}_{\theta}-\frac{\left(2.176v_{z}+0.871z\Omega_{K,a}\right)}{\tau_{i}}{\bf
e}_{z}\,,$ (1)

but we replaced the ‘migration timescale’ $\tau_{m}=-\dot{L}/L$ in the first
term with the semi-major axis evolution timescale $\tau_{a}=-\dot{a}/a$. The
motivation behind replacing $\tau_{m}$ with $\tau_{a}$ in Matsumura et al.
(2017) was to avoid an artificial outward migration at the high eccentricity,
supersonic regime. Here, $v_{K,a}$ and $\Omega_{K,a}$ are the Keplerian
orbital speed and the corresponding angular speed evaluated at the semi-major
axis $a$, while ${\bf v}$ is the planetary velocity evaluated at the
instantaneous orbital radius $r$, with $v_{r}$, $v_{\theta}$, and $v_{z}$
being its polar coordinate components defined with the unit vectors ${\bf
e}_{r}$, ${\bf e}_{\theta}$, and ${\bf e}_{z}$, respectively. Ida et al.
(2020) showed that the azimuthal component of eccentricity damping is
implicitly included in the first term of Equation 1, $-\frac{{\bf
v}}{\tau_{m}}$, and its explicit addition as the second term is not necessary.

In this work, we adopted the equation of motion proposed by Ida et al. (2020),
since the equation describes planet-disc interactions well both in subsonic
and supersonic regimes compared to the other equations proposed so far. They
introduced an intuitive model of planet-disc interactions based on dynamical
friction by combining the work of Muto et al. (2011) in the supersonic regime
with that of Tanaka & Ward (2004) in the subsonic regime. The equation reads
as follows:

$\displaystyle\frac{d{\bf v}}{dt}$ $\displaystyle=$
$\displaystyle-\frac{v_{K}}{2\tau_{a}}{\bf
e}_{\theta}-\frac{v_{r}}{\tau_{e}}{\bf
e}_{r}-\frac{v_{\theta}-v_{K}}{\tau_{e}}{\bf
e}_{\theta}-\frac{v_{z}}{\tau_{i}}{\bf e}_{z}\,.$ (2)

We note that $v_{K}$ is the Keplerian orbital speed evaluated at the
instantaneous orbital radius $r$. In this model, the evolution timescales for
the semi-major axis, eccentricity, and inclination are expressed as follows:

$\displaystyle\tau_{a}$ $\displaystyle=-\frac{t_{{\rm
wave}}}{2h_{g}^{2}}\left[\frac{\Gamma_{L}}{\Gamma_{0}}\left(1-\frac{1}{\pi}\frac{\Gamma_{L}}{\Gamma_{0}}\sqrt{\hat{e}^{2}+\hat{i}^{2}}\right)^{-1}\right.$
$\displaystyle\qquad\qquad\left.+\frac{\Gamma_{C}}{\Gamma_{0}}\exp\left(-\frac{\sqrt{e^{2}+i^{2}}}{e_{f}}\right)\right]^{-1},$
(3) $\displaystyle\tau_{e}$ $\displaystyle=\frac{t_{{\rm
wave}}}{0.780}\,\left(1+\frac{1}{15}\left(\hat{e}^{2}+\hat{i}^{2}\right)^{3/2}\right),$
(4) $\displaystyle\tau_{i}$ $\displaystyle=\frac{t_{{\rm
wave}}}{0.544}\,\left(1+\frac{1}{21.5}\left(\hat{e}^{2}+\hat{i}^{2}\right)^{3/2}\right)\
,$ (5)

where $h_{g}=H_{g}/a$ is the gas disc’s aspect ratio, and eccentricities and
inclinations scaled with $h_{g}$ are defined as $\hat{e}=e/h_{g}$ and
$\hat{i}=i/h_{g}$, respectively. The exponential factor in $\tau_{a}$ is
introduced with $e_{f}=0.01+h_{g}/2$ following Fendyke & Nelson (2014), so
that the contribution from the corotation torque disappears in the supersonic
regime. Also, $\Gamma/\Gamma_{0}$ represents a normalised torque, and the
subscripts $L$ and $C$ respectively correspond to Lindblad and corotation
torques from Paardekooper et al. (2011), with the normalisation torque being

$\Gamma_{0}=\left(\frac{M_{p}}{M_{*}}\right)^{2}\Sigma_{g}a^{4}\,h_{g}^{-2}\Omega_{K,a}^{2}.$
(6)

Finally, $t_{{\rm wave}}$ is the characteristic time of the orbital evolution
by Tanaka & Ward (2004):

$t_{{\rm
wave}}=\left(\frac{M_{*}}{M_{p}}\right)\left(\frac{M_{*}}{\Sigma_{g}a^{2}}\right)h_{g}^{4}\Omega_{K,a}^{-1}\
,$ (7)

where $M_{p}$ and $M_{*}$ are planetary and stellar masses, respectively, and
$\Sigma_{g}$ is the gas disc’s surface mass density. We note that $t_{{\rm
wave}}$ is defined by using the semi-major axis $a$, not the instantaneous
orbital radius $r$, so the timescales defined in Equations 3-5 are orbit-
averaged timescales.

The above equations only hold for fully embedded type I migrators and need to
be modified for gap-opening type II migrators. For planet migration, Kanagawa
et al. (2018) have that type II migration is merely type I migration with the
reduced surface mass density in the gap, and showed that such migration is
well described by the following equation:

$\tau_{a}^{\prime}\simeq\frac{\Sigma_{g}}{\Sigma_{{\rm
gap}}}\tau_{a}\simeq\left(1+0.04K\right)\tau_{a}\,,$ (8)

where
$K=\left(\frac{M_{p}}{M_{*}}\right)^{2}\left(\frac{H_{g}}{a}\right)^{-5}\alpha_{{\rm
turb}}^{-1}$. Since there is no consensus on eccentricity and inclination
damping in the type II regime, we adjusted eccentricity and inclination
damping timescales in a similar manner ($\tau_{e}^{\prime}=(1+0.04K)\tau_{e}$
and $\tau_{i}^{\prime}=(1+0.04K)\tau_{i}$) and implemented the equation of
motion 2 with timescales replaced by $\tau_{a}^{\prime}$, $\tau_{e}^{\prime}$,
and $\tau_{i}^{\prime}$ in our code. In this prescription, the eccentricity
and inclination damping timescales are always short compared to the migration
timescale. This is an important condition for the ‘eccentricity trap’ (Ogihara
et al., 2010), where planets can be resonantly trapped near the disc’s inner
edge. The condition also avoids spurious outward migration at the disc’s edge.

Figure 1 shows these timescales as a function of a protoplanetary mass for a
circular and coplanar case. The vertical line indicates a critical mass
$M_{{\rm crit}}$ above which planet migration switches from type I to type II
regimes, and it is defined as the mass having the minimum timescale of Eq. 8:

$\displaystyle M_{{\rm crit}}$ $\displaystyle=$
$\displaystyle\left[\frac{1}{0.04}\alpha_{{\rm
turb}}\left(\frac{H_{g}}{a}\right)^{5}\right]^{1/2}M_{*}$ (9)
$\displaystyle\simeq$ $\displaystyle 9.3\left(\frac{\alpha_{{\rm
turb}}}{10^{-4}}\right)^{1/2}\left(\frac{H/a}{0.05}\right)^{5/2}\left(\frac{M_{*}}{M_{\odot}}\right)\;M_{\oplus}.$
(10)

Since this critical mass gives $K=\frac{1}{0.04}$ and thus $\frac{\Sigma_{{\rm
gap}}}{\Sigma_{g}}\simeq\frac{1}{1+0.04K}=\frac{1}{2}$, the transition occurs
when the gap depth is 50% of the nominal surface mass density (also see
Johansen et al., 2019). The transition mass is sensitive to the disc’s aspect
ratio, and we have $M_{{\rm crit}}=2.6\,M_{\oplus}$ for $\frac{H}{a}=0.03$,
$\alpha_{{\rm turb}}=10^{-4}$, and $M_{*}=M_{\odot}$.

Figure 1: Comparison of evolution timescales of semi-major axis $\tau_{a}$
(blue), eccentricity $\tau_{e}$ (orange), and inclination $\tau_{i}$ (green)
for a planet at 3.16 au with the circular and coplanar orbit in Disc 6\. Solid
and dashed blue lines are from Equations 8 and 3, respectively. The vertical
black dashed line indicates where migration timescale takes the minimum value
(Eq. 10).

###  2.2 Planet formation and disc evolution

The mass growth rate of a protoplanet $\dot{M}_{p}$ is written as follows:

$\dot{M}_{p}=\dot{M}_{{\rm gas}}+\dot{M}_{{\rm core}},$ (11)

where $\dot{M}_{{\rm gas}}$ and $\dot{M}_{{\rm core}}$ represent gas and core
accretion rates of protoplanets, respectively. In this study, we assume that a
protoplanetary core grows via pebble accretion and mutual collisions among
themselves, and we do not take the effect of planetesimal accretion into
account. This choice is justified when the pebble accretion is 3D (i.e. when
the impact radius of a core is small compared to the scale height of the
pebble disc $b<h_{p}$ ), because the timescale of 3D pebble accretion is much
shorter than that of the planetesimal accretion of the oligarchic growth stage
(Ida et al., 2016). The pebble and planetesimal accretion timescales may
become comparable to each other once the impact parameter of the core becomes
larger than the pebble disc’s scale height $b\gg h_{p}$ (Ida et al., 2016).

The core growth rate via pebble accretion is written as

$\dot{M}_{{\rm core}}=\epsilon\dot{M}_{F},$ (12)

by using the accretion efficiency $\epsilon$ and the pebble mass flux
$\dot{M}_{F}$. In the following, we describe the model of the stellar mass
accretion rate $\dot{M}_{*}$ and how $\dot{M}_{F}$ is related to it. We also
discuss the choice of pebble accretion efficiencies, and the pebble isolation
masses. Finally, we present the model of the gas accretion rate onto a
protoplanetary core $\dot{M}_{{\rm gas}}$.

####  2.2.1 Stellar mass flux $\dot{M}_{*}$

The evolution of a gas disc is represented by the stellar mass flux
$\dot{M}_{*}$, which is also related to the evolution of a dust disc
$\dot{M}_{F}$ as it is discussed in the following sub-section. In this work,
we used the same disc model as Matsumura et al. (2017), which is adopted from
Ida et al. (2016). The disc’s midplane temperature $T$ and the gas surface
mass density $\Sigma_{g}$ are modelled as the power-law functions of the
orbital radius $r$ as $T\propto r^{-q}$ and $\Sigma_{g}\propto r^{-p}$. Since
the inner disc is heated by viscous dissipation while the outer disc is heated
by the irradiation from the central star (e.g. Hueso & Guillot, 2005; Oka et
al., 2011), both $T$ and $\Sigma_{g}$ are dual power-law functions in our disc
model.

The disc evolution is described by the diffusion equation for the disc’s
surface mass density $\Sigma_{g}$:

$\frac{\partial\Sigma_{g}}{\partial t}=-\frac{1}{r}\frac{\partial}{\partial
r}\left[3r^{1/2}\frac{\partial}{\partial r}\left(\Sigma_{g}\nu
r^{1/2}\right)\right]\ ,$ (13)

where $\nu=\alpha_{{\rm acc}}c_{s}H$ is the disc’s ‘viscosity’ representing
the accretion. Assuming that the accretion is steady and that the effective
viscosity is written as a power-law function in radius as $\nu\propto r^{p}$
(Lynden-Bell & Pringle, 1974; Hartmann et al., 1998), the mass accretion rate
is related to the surface mass density as follows:

$\dot{M}_{*}\simeq 3\pi\Sigma_{g}\nu=\frac{M_{D,0}}{2(2-p)t_{{\rm
diff}}}\tilde{t}^{-(5/2-p)/(2-p)}\ .$ (14)

Here, $M_{D,0}$ is the initial disc mass, $\tilde{t}\equiv t/t_{{\rm diff}}+1$
with the diffusion timescale of

$t_{{\rm diff}}\equiv\frac{1}{3(2-p)^{2}}\frac{r_{D}^{2}}{\nu_{D}}\ ,$ (15)

where $r_{D}$ is the characteristic disc size and $\nu_{D}$ is the
corresponding viscosity.

As can be seen from the equation, there are several parameters we could use to
specify the disc evolution, such as initial values of $\dot{M}_{*}$, $M_{D}$,
and $\Sigma_{g,D;}$ as well as $r_{D}$, $\alpha_{{\rm acc}}$, and $t_{{\rm
diff}}$. However, they are dependent on one another, and we need to choose
three parameters out of these to specify a disc model. For example, Ida & Lin
(2008) chose $M_{D,0}$, $t_{{\rm diff}}$, and $\alpha_{{\rm acc}}$ as control
parameters since they are relevant to planetesimal accretion as well as gap-
opening criteria. For pebble accretion, however, the disc radius $r_{D}$ is
important because the pebble mass flux decreases sharply once the pebble
formation front reaches the outer disc radius (e.g. Sato et al., 2016). In
this work, we considered a few different sets of $M_{D,0}$ and $t_{{\rm
diff}}$, along with $r_{D}=100\,{\rm au}$. The choice of the characteristic
disc size is rather arbitrary, but it is motivated by a typical size of
observed protoplanetary discs (e.g. Andrews et al., 2010; Andrews & Williams,
2007; Vicente & Alves, 2005). The $\alpha_{{\rm acc}}$ is calculated from
these values as follows:

$\alpha_{{\rm acc}}=\frac{h_{D}^{-2}}{6\pi\left(2-p\right)^{2}}\frac{t_{{\rm
orb,D}}}{t_{{\rm diff}}},$ (16)

where $h_{D}=H_{D}/r_{D}$ is the disc’s aspect ratio at $r_{D}$, and $t_{{\rm
orb,D}}$ is the corresponding orbital period. In Matsumura et al. (2017), we
defined the stellar mass accretion rate and the alpha parameter independently
of each other and calculated the surface mass density based on these
parameters and the disc temperature models in Ida et al. (2016). In this work,
we instead calculated the stellar mass accretion rate for given disc masses
and diffusion timescales (and thus $\alpha_{{\rm acc}}$), and we used these
with the temperature model to determine the surface mass density profile. In
this sense, our current disc model is more self-consistent compared to the
previous one.

As discussed in Section 1, recent studies showed that disc accretion is driven
by the magnetic disc wind rather than by the MRI turbulence (e.g. Bai & Stone,
2013; Bai, 2017).The follow-up studies proposed to replace the classical
$\alpha$ disc model with the one incorporating the disc wind (Bai, 2016;
Suzuki et al., 2016), and the model has been adopted by N-body studies such as
Ogihara et al. (2017). Some other works, on the other hand, opt to use
two-$\alpha$ disc models by representing the mass accretion driven mainly by
the disc wind, and the local disc evolution driven by turbulence with two
different alphas (e.g. Ida et al., 2018; Johansen et al., 2019; Bitsch et al.,
2019). In this work, we also adopted this two-$\alpha$ disc model for
simplicity, and we represent the wind-driven disc accretion with $\alpha_{{\rm
acc}}$ and the disc turbulence with $\alpha_{{\rm turb}}$ . Compared to the
sophisticated wind-driven accretion models by Bai (2016) and Suzuki et al.
(2016), this corresponds to approximating the combined effects of the wind-
driven and turbulence-driven accretions with a single effective $\alpha_{{\rm
acc}}$, and ignoring the effect of the wind mass loss.

Finally, to choose our disc models, we constrained the combination of
$\left(M_{D,0},\,t_{{\rm diff}}\right)$ by using the observed stellar mass
accretion rates. Figure 2 shows eight disc models using Equation 14, along
with the observed stellar mass accretion rates from Sicilia-Aguilar et al.
(2010) (data courtesy of Sicilia-Aguilar)111We note that we assume the solar
mass for all of our simulations, while observed data include GKM stars.. Here,
we did not attempt to find the best fit to all the data, as this was done in
Hartmann et al. (1998), but instead we tried to find a set of reasonable
$\left(M_{D,0},\,t_{{\rm diff}}\right)$ that makes $\dot{M}_{*}$ go through at
least some part of the distribution of observed accretion rates. This is
because, as Hartmann et al. (1998) pointed out, the best fit to all data does
not necessarily represent a typical evolution of the stellar mass accretion
rates.

As can be seen in the figure, our chosen disc models tend to have long
lifetimes. Towards the oldest ages (about a few to several tens of Myr), the
mass accretion rate data are likely dominated by long-lived, potentially less
common protoplanetary discs. However, a recent study shows that $\sim 30\,\%$
of stars may have disc lifetimes longer than 10 Myr (Pfalzner et al., 2014),
making it possible that these long-lived discs may not be so uncommon.
Therefore, in this work, we consider a relatively wide range of disc
lifetimes. The combination of $\left(M_{D,0},\,t_{{\rm diff}}\right)$ and the
corresponding $\alpha_{{\rm acc}}$ for eight disc models are summarised in
Table 1. To mimic the photoevaporation effect, we reduced the mass accretion
rate exponentially once it became lower than the critical value of
$10^{-9}\,M_{\odot}/{\rm yr}$. The choice of this critical value is arbitrary,
but we chose it to be close to the minimum observed mass accretion rate seen
in the figure. We note that, as seen in the table, $\alpha_{{\rm acc}}$
defined with our choices of parameters ($t_{{\rm diff}}=0.1-10\,$Myr and
$r_{D}=100\,{\rm AU}$) is higher than our default turbulent viscosity alpha
value of $\alpha_{{\rm turb}}=10^{-4}$, which is consistent with our
assumption that the mass accretion is dominated by the disc’s wind effects. In
Section 3.2.1, we discuss the effects of other values of $\alpha_{{\rm
turb}}$.

Figure 2: Left: Stellar mass accretion rates $\dot{M}_{*}$ for disc models 1-8
from Table 1 (Eq 14). The accretion is exponentially reduced after the
critical mass accretion rate of $10^{-9}\,M_{\odot}/{\rm yr}$ is reached
(dashed lines). The black circles with error bars are observed stellar
accretion rates from (Sicilia-Aguilar et al., 2010) (data courtesy of Sicilia-
Aguilar). Right: Corresponding pebble mass fluxes $\dot{M}_{F}$ (Eq 25).

####  2.2.2 Pebble mass flux $\dot{M}_{F}$

To determine the pebble accretion rate by a protoplanet, we need to estimate
the pebble mass flux. Lambrechts & Johansen (2014) wrote the pebble mass flux
as follows:

$\dot{M}_{F}=2\pi r_{{\rm pf}}\Sigma_{p}\frac{dr_{{\rm pf}}}{dt},$ (17)

which describes the pebble mass swept up by the pebble formation front per
unit time. Here, $r_{{\rm pf}}$ is the radius of the pebble formation front,
and $\Sigma_{p}$ is the pebble surface mass density. Since the dust particles
are likely to convert to pebbles over a long period of time, the dust surface
mass density $\Sigma_{d}$ should be different from that of pebbles
$\Sigma_{p}$. At the pebble formation front, however, we assume
$\Sigma_{d}=\Sigma_{p}$ so that all dust particles are converted to pebbles
there.

Another effect we need to take into account is that the pebble mass flux
decreases drastically once the pebble formation front reaches the outer disc
radius (e.g. Sato et al., 2016). Ida et al. (2019) fitted the numerical
simulations by Sato et al. (2016) and proposed the following relation:

$\Sigma_{pg}=\Sigma_{pg,0}\left(1+\frac{t}{t_{{\rm pf}}}\right)^{-\gamma},$
(18)

where $\Sigma_{pg}=\Sigma_{p}/\Sigma_{g}$, the subscript 0 indicates the
initial value, $\gamma=1+\gamma_{a}\left(\frac{300\,{\rm au}}{r_{D}}\right)$
with $\gamma_{a}\sim 0.15$, $r_{D}$ is the disc radius, and $t_{{\rm pf}}$ is
the time it takes for the pebble front to reach $r_{D}$.

In the Epstein regime (i.e. when a dust particle is small enough compared to
the mean free path: $R\lesssim\frac{9}{4}\lambda_{{\rm mfp}}$), the growth
timescale of the particle can be approximated as follows by assuming that the
relative velocities between particles are dominated by turbulence (Ida et al.,
2019):

$\displaystyle t_{{\rm grow}}$ $\displaystyle\simeq$
$\displaystyle\frac{4}{\sqrt{3\pi}}\Sigma_{pg,0}^{-1}\Omega^{-1}$ (19)
$\displaystyle=$
$\displaystyle\frac{200}{\sqrt{3\pi^{3}}}\left(\frac{\Sigma_{pg,0}}{0.01}\right)^{-1}\left(\frac{M_{*}}{M_{\odot}}\right)^{-1/2}\left(\frac{r}{{\rm
au}}\right)^{3/2}\>{\rm yr}.$

The growth timescale of pebbles from ${\rm\mu m}$ to cm sizes can be estimated
as $t_{{\rm grow,peb}}=\ln\frac{R_{{\rm peb}}}{R_{0}}\cdot t_{{\rm grow}}$
(Lambrechts & Johansen, 2014; Ida et al., 2016), where $R_{0}$ and $R_{{\rm
peb}}$ are the initial size of dust particles and the ‘final’ size of pebbles
where they start to migrate, respectively. We adopted the nominal values of
$R_{{\rm peb}}\sim 10\,{\rm cm}$ and $R_{0}\sim 1\,{\rm\mu m}$ in the above
equation.

From these, we can define the timescale for the pebble formation front to
reach the outer disc radius as follows:

$t_{{\rm pf}}\simeq\ln\frac{R_{{\rm
peb}}}{R_{0}}\cdot\frac{200}{\sqrt{3\pi^{3}}}\left(\frac{\Sigma_{pg,0}}{0.01}\right)^{-1}\left(\frac{M_{*}}{M_{\odot}}\right)^{-1/2}\left(\frac{r_{D}}{{\rm
au}}\right)^{3/2}\>{\rm yr}.$ (20)

Thus, the pebble front location is in turn given by

$\displaystyle\left(\frac{r_{{\rm pf}}}{{\rm au}}\right)$
$\displaystyle\simeq$ $\displaystyle\left(\frac{1}{\ln\frac{R_{{\rm
peb}}}{R_{0}}}\frac{\sqrt{3\pi^{3}}}{200}\right)^{2/3}\left(\frac{\Sigma_{pg,0}}{0.01}\right)^{2/3}\left(\frac{M_{*}}{M_{\odot}}\right)^{1/3}\left(\frac{t}{{\rm
yr}}\right)^{2/3},$ (21)

and we also have

$\displaystyle\frac{dr_{{\rm pf}}}{dt}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\ln\frac{R_{{\rm
peb}}}{R_{0}}}\frac{\sqrt{3\pi^{3}}}{200}\right)^{2/3}\left(\frac{\Sigma_{pg,0}}{0.01}\right)^{2/3}\left(\frac{M_{*}}{M_{\odot}}\right)^{1/3}$
(22) $\displaystyle\cdot\frac{2}{3}\left(\frac{t}{{\rm
yr}}\right)^{-1/3}\,{\rm{\rm au\,}yr^{-1}}.$

By putting these together, the pebble mass flux is written as follows:

$\displaystyle\dot{M}_{F}$ $\displaystyle=$ $\displaystyle 2\pi r_{{\rm
pf}}\Sigma_{g}\frac{dr_{{\rm
pf}}}{dt}\cdot\Sigma_{pg,0}\left(1+\frac{t}{t_{{\rm pff}}}\right)^{-\gamma}.$
(23)

Finally, we calculated the pebble mass flux in our disc model. Since the
pebble formation front is in the outer disc, we used the gas surface mass
density in the irradiation region from Ida et al. (2016):

$\displaystyle\Sigma_{g,{\rm irr}}$ $\displaystyle=$ $\displaystyle
2514\left(\frac{T_{2}}{150\,{\rm K}}\right)^{-1}L_{*0}^{-2/7}M_{*0}^{9/14}$
$\displaystyle\cdot\alpha_{3}^{-1}\dot{M}_{*8}\left(\frac{r}{{\rm
au}}\right)^{(q_{2}-3/2)}\;{\rm g/cm^{2}},$

where $T_{2}$ is a characteristic disc temperature in the irradiation region
and $T_{2}=150\,$K in our default model, $q_{2}=3/7$ is the power in $T\propto
r^{-q}$, the $L_{*0}=L_{*}/L_{\odot}$ and $M_{*0}=M_{*}/M_{\odot}$ are the
stellar luminosity and mass scaled to the solar values,
$\alpha_{3}=\alpha_{{\rm acc}}/10^{-3}$ is the disc accretion $\alpha_{{\rm
acc}}$ scaled to a typical value, and
$\dot{M}_{*8}=\dot{M}_{*}/\left(10^{-8}\,{\rm M_{\odot}/yr}\right)$.

Substituting Equations 18, 20-22, and 23 into Equation 23, we obtain the
following:

$\displaystyle\dot{M}_{F}$ $\displaystyle=$ $\displaystyle 2.066\times
10^{-2}\left(\frac{\ln\frac{R_{{\rm peb}}}{R_{0}}}{\ln
10^{4}}\right)^{-1}\left(\frac{T_{2}}{150\,{\rm K}}\right)^{-1}$ (25)
$\displaystyle
L_{*0}^{-2/7}M_{*0}^{8/7}\alpha_{3}^{-1}\dot{M}_{*8}\left(\frac{\Sigma_{pg,0}}{0.01}\right)^{2}\left(\frac{r_{D}}{{\rm
au}}\right)^{q_{2}-1}$ $\displaystyle\left(\frac{t}{t_{{\rm
pff}}}\right)^{\frac{2}{3}\left(q_{2}-1\right)}\left(1+\frac{t}{t_{{\rm
pff}}}\right)^{-\gamma}\,{\rm M_{\oplus\,}yr^{-1}}\,.$

The evolution of this equation is shown in the right panel of Figure 2.

In our simulations, we considered five different stellar metallicities:
${\rm[Fe/H]}=\left(-0.5,\,-0.3,\,0.0,\,0.3,\,0.5\right).$ Thus, we covered a
range of metallicities of planet-hosting stars that are related to the dust-
to-gas surface mass densities $\Sigma_{dg}=\Sigma_{d}/\Sigma_{g}$ as follows:

$\frac{\Sigma_{dg}}{\Sigma_{dg,\odot}}\simeq 10^{{\rm[Fe/H]}},$ (26)

where $\Sigma_{dg,\odot}=0.01$ is the value for the Solar System. Thus, in our
model, the pebble mass flux is proportional to the square of the dust-to-gas
surface mass density ratio. The dependence is not exactly the same, but
similar to that of Lambrechts & Johansen (2014), where
$\dot{M}_{F}\propto\Sigma_{dg,0}^{5/3}$.

The last column of Table 1 shows the total masses of pebbles accreting towards
the star during the simulations for the solar metallicity case
${\rm[Fe/H]}=0.0$. The total masses vary over $19.4-597\,M_{E}$ for
${\rm[Fe/H]}=0.0$ in our disc models, and they vary by a factor of 3 above and
below these values for the entire range of metallicities we considered: $\sim
6.5-1791\,M_{E}$. In comparison, Bitsch et al. (2019) obtained total pebble
masses of $70-700\,M_{E}$. Tychoniec et al. (2018) estimated that dust disc
masses of Class 0 objects vary over $10-6000\,M_{E,}$ where a typical value is
$248\,M_{E}$ and less than 10% of systems have masses above $1000\,M_{E}$.
Therefore, the range of total pebble masses we considered are consistent with
observations.

####  2.2.3 Pebble accretion efficiency $\epsilon$

The protoplanetary cores are exposed to the flux of pebbles, but they only
accrete a fraction of the incoming flux. As stated before, we write the pebble
accretion rate as $\dot{M}_{{\rm core}}=\epsilon\dot{M}_{F}$. Matsumura et al.
(2017) adopted the accretion efficiency defined by Ida et al. (2016), which is
valid in the settling regime:

$\epsilon_{{\rm IGM16}}={\rm min}\left(\frac{C_{\epsilon}\zeta^{-1}\chi
b^{2}}{4\sqrt{2\pi}\tau_{s}h_{p}}\left(1+\frac{3b}{2\chi\eta}\right),\,1\right)\
.$ (27)

The $1$ in parentheses ensures that the accretion efficiency does not become
larger than unity. In this equation, $\zeta$ and $\chi$ are functions of the
Stokes number $\tau_{s}$ and defined as

$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1+4\tau_{s}^{2}}}{1+\tau_{s}^{2}},$ (28)
$\displaystyle\zeta$ $\displaystyle=$ $\displaystyle\frac{1}{1+\tau_{s}^{2}}\
.$ (29)

The pebble aspect ratio is defined as $h_{p}=H_{p}/r,$ where $H_{p}$ is the
pebble scale height (Ida et al., 2016):

$H_{p}\simeq\left(1+\frac{\tau_{s}}{\alpha_{{\rm turb}}}\right)^{-1/2}H_{g}.$
(30)

Similarly, $b=B/r$ and $B$ is an impact parameter for a pebble passing by the
protoplanet with mass $M_{p}$:

$B\simeq\min\left(\sqrt{\frac{3\tau_{s}^{1/3}R_{H}}{\chi\eta
r}},\,1\right)\cdot 2\kappa_{{\rm OK12}}\tau_{s}^{1/3}R_{H}\,,$

where the left and the right terms in parentheses correspond to the Bondi and
Hill regimes, respectively. The
$R_{H}=\left(\frac{M_{p}}{3M_{*}}\right)^{1/3}r$ is the Hill radius of the
protoplanet, and the $\kappa_{{\rm OK12}}$ is a reduction factor of the
accretion proposed by Ormel & Kobayashi (2012) for $\tau_{s}\gg 1$ as:

$\kappa_{{\rm
OK12}}=\exp\left(-\left(\frac{\tau_{s}}{\min\left(2,\,\tau_{s}^{*}\right)}\right)^{0.65}\right),$

with $\tau_{s}^{*}=4\left(M_{p}/M_{*}\right)/\eta^{3}$. Also, $C_{\epsilon}$
is defined as

$C_{\epsilon}=\min\left(\sqrt{\frac{8}{\pi}}\frac{h_{p}}{b},\,1\right)\ ,$
(31)

where the left and right terms represent 2D and 3D accretion, respectively.
For the simulations, we scaled $C_{\epsilon}$ by a factor $C_{\epsilon,i}$ to
take account of the orbital inclination effect. By assuming that the pebble
volume density scales as $\exp\left(-z^{2}/\left(2\,H_{p}^{2}\right)\right)$,
the surface mass density between $z\pm B$ scales with:

$\displaystyle\int_{z-B}^{z+B}\exp\left(-\frac{z^{\prime
2}}{2\,H_{p}^{2}}\right)dz^{\prime}=$
$\displaystyle\frac{\sqrt{2\pi}H_{p}}{2}\left({\rm
erf\left(\frac{z+B}{\sqrt{2}H_{p}}\right)-{\rm
erf\left(\frac{z-B}{\sqrt{2}H_{p}}\right)}}\right).$ (32)

By normalising this factor with the value at $z=0$, we obtain the following:

$C_{\epsilon,i}=\frac{1}{2}\frac{\left({\rm
erf\left(\frac{z+B}{\sqrt{2}H_{p}}\right)-{\rm
erf\left(\frac{z-B}{\sqrt{2}H_{p}}\right)}}\right)}{{\rm
erf}\left(\frac{B}{\sqrt{2}H_{p}}\right)}.$ (33)

More recently, Ormel & Liu (2018) studied the 3D pebble accretion efficiency
by considering the effects of eccentricity, inclination, and disc turbulence:

$\epsilon_{{\rm OL18}}=\frac{A_{3}}{\eta h_{{\rm
p,eff}}}\left(\frac{M_{p}}{M_{*}}\right)f_{{\rm set}}^{2}\ ,$ (34)

where $A_{3}=0.39$ is a fitting constant to their pebble accretion
simulations. The effective pebble aspect ratio is defined as

$h_{{\rm
p,eff}}\sim\sqrt{h_{p}^{2}+\frac{\pi\,i_{p}^{2}}{2}\left(1-\exp\left[-\frac{i_{p}}{2h_{p}}\right]\right)}\
,$ (35)

where $i_{p}$ is the inclination of a planetary orbit. The settling fraction
is determined by integrating the accretion probability over the velocity
distribution and is written as follows for turbulence operating only in the
vertical direction:

$f_{{\rm set}}=\exp\left[-a_{{\rm set}}\left(\frac{\Delta
v_{y}^{2}}{v_{*}^{2}}+\frac{\Delta v_{z}^{2}}{v_{*}^{2}+a_{{\rm
turb}}\sigma_{P,z}^{2}}\right)\right]\frac{v_{*}}{\sqrt{v_{*}^{2}+a_{{\rm
turb}}\sigma_{P,z}^{2}}}\ ,$ (36)

where $a_{{\rm set}}=0.5$ and $a_{{\rm turb}}=0.33$ are other fit constants,
$\Delta v_{y}$ and $\Delta v_{z}$ are azimuthal and vertical approach
velocities, respectively, and $\sigma_{P,z}$ is the vertical component of
pebble root mean square velocity.

The left panel of Figure 3 compares the accretion efficiencies for a $\sim
0.1M_{\oplus}$ planet from Ida et al. (2016) (crosses) and from Ormel & Liu
(2018) (circles) as a function of $\tau_{s}$. The estimated accretion
efficiencies are similar for Ida et al. (2016) and Ormel & Liu (2018) over a
wide range of $\tau_{s}$, though the values from Ida et al. (2016) are
generally slightly higher than those from Ormel & Liu (2018). In this work, we
tested both types of the accretion efficiencies, and confirmed that the
general trends are similar in both cases except that the mass growth is less
efficient in the cases of Ormel & Liu (2018). In this paper, we focus on the
simulations with the efficiency by Ida et al. (2016) $\epsilon_{{\rm IGM16}}$,
but we discuss the effects of that by Ormel & Liu (2018) $\epsilon_{{\rm
OL18}}$ briefly in Section 3.2.2.

Figure 3: Left: Pebble accretion efficiency $\epsilon$ and the Stokes number
for different values of the vertical turbulence strength $\alpha_{{\rm
turb}}$. The figure uses the planet-to-star mass ratio of $3\times 10^{-7}$
(i.e. $\sim 0.1\,M_{E}$ for a Sun-like star), the gas disc aspect ratio of
$h_{g}=0.03$, and the disc radial pressure gradient of $\eta=1.0\times
10^{-3}$ as in Figure 4 of Ormel & Liu (2018). Circles show $\epsilon$ from
Ormel & Liu (2018) and crosses show corresponding $\epsilon$ from Ida et al.
(2016). Right: Pebble isolation masses estimated by Ataiee et al. (2018)
(orange) and Bitsch et al. (2018) (green), compared with our default case from
Ida et al. (2016, blue). Also plotted are the critical migration transition
masses (Equation 10, red).

####  2.2.4 Pebble isolation mass

The pebble isolation mass (PIM) is the protoplanetary mass at which the embryo
becomes large enough to perturb the gas disc, modify the pressure gradient so
that pebbles locally feel a tailwind rather than a headwind, and thus to halt
the inward radial drift of pebbles (Morbidelli & Nesvorny, 2012). This
critical mass marks the end of the pebble accretion stage and thus is one of
the key parameters to determine the outcome of planet formation. There are
different formulations for the PIM in the literature, which we briefly discuss
below.

Lambrechts et al. (2014) estimated the PIM as $M_{{\rm iso}}\sim
20\left(\frac{a}{5\,{\rm au}}\right)^{3/4}M_{\oplus}$ from hydrodynamic
simulations, which we adopted for Matsumura et al. (2017). However, this study
did not take account of the effects of the disc’s turbulent viscosity. More
recently, Ataiee et al. (2018) and Bitsch et al. (2018) extended this study
and investigated the dependence of the PIM on the disc aspect ratio, the
pressure gradient of the disc as well as the viscosity. Ataiee et al. (2018)
used 2D gas and dust hydrodynamic simulations, while Bitsch et al. (2018) used
3D gas-only hydrodynamic simulations with 2D integrations of particle
trajectories. As discussed in Ataiee et al. (2018) and also shown in Figure 3,
the two studies have similar PIM trends, but details are different. Moreover,
the differences become larger for the lower values of the viscosity
$\alpha_{{\rm turb}}$.

Bitsch et al. (2018) performed 3D hydro simulations over $\alpha_{{\rm
turb}}=\left[2\times 10^{-4},\>6\times 10^{-3}\right]$ to determine the PIM as
follows:

$\displaystyle M_{{\rm PIM,B18}}$
$\displaystyle\simeq\left(25+10.5\frac{\alpha_{{\rm
turb}}}{\tau_{s}}\right)f_{{\rm fit}}M_{\oplus},$ (37) $\displaystyle f_{{\rm
fit}}$
$\displaystyle=\left(\frac{h_{g}}{0.05}\right)^{3}\left(0.34\left(\frac{\log
10^{-3}}{\log\alpha_{{\rm turb}}}\right)^{4}+0.66\right)$
$\displaystyle\quad\left(1-\frac{\frac{\partial\ln P}{\partial\ln
r}+2.5}{6}\right).$ (38)

A potential problem with this fit formula is that the saturation might not
quite have been achieved for the lower end of $\alpha_{{\rm turb}}$, as
indicated in Figure 3 of Johansen et al. (2019).

On the other hand, Ataiee et al. (2018) performed 2D hydro simulations over
$\alpha_{{\rm turb}}=\left[5\times 10^{-4},\>1\times 10^{-2}\right]$ and
proposed the following PIM:

$\displaystyle M_{{\rm PIM,A18}}$ $\displaystyle\simeq
h_{g}^{3}\sqrt{37.3\alpha_{{\rm turb}}+0.01}$
$\displaystyle\quad\cdot\left(1+0.2\left(\frac{\sqrt{\alpha_{{\rm
turb}}}}{h_{g}}\sqrt{\frac{1}{\tau_{s}^{2}}+4}\right)^{0.7}\right)M_{*.}$ (39)

As seen in their Figure 8, the expression agrees well with the 2D hydro
simulations in the low-viscosity limit ($\alpha_{{\rm turb}}\sim 10^{-3}$),
while the difference is about a factor of 3 for $\alpha_{{\rm turb}}\sim
10^{-2}$.

The right panel of Figure 3 compares the pebble isolation mass model by Ataiee
et al. (2018) with that of Bitsch et al. (2018). The figure also shows the
transition mass from type I to type II migration as shown in the previous sub-
section. Since the PIM measures the first appearance of the zero pressure
gradient outside the planetary gap, the required depth is about $15\%$
(Lambrechts et al., 2014) and is shallower than that for the migration
transition of $\sim 50\%$ (Johansen et al., 2019). In other words, the planets
are expected to first obtain the PIM and then reach the migration transition
mass via gas accretion (or planetesimal accretion or core-core collisions).
This implies that, unless the gas accretion is fast enough, the protoplanets
are likely to migrate toward the central star on type I migration timescale
(Johansen et al., 2019). However, since the PIM decreases towards the central
star, it is likely for such type I migrators to eventually start gas
accretion, and for migration to slow down as a result.

The PIM by Bitsch et al. (2018) is lower than the migration transition mass
within $\sim 10\,$au for $\alpha_{{\rm turb}}=1\times 10^{-3}$ (see dashed
green and dashed red lines, respectively), but it is higher for $\alpha_{{\rm
turb}}=1\times 10^{-4}$ (see solid green and solid red lines, respectively).
The latter is counter-intuitive, since it indicates that a deeper gap is
required for migration transition than for reaching the PIM. On the other
hand, for the PIM model by Ataiee et al. (2018), the PIM is generally lower
than the migration transition mass. Thus, for this work, we adopted the PIM by
Ataiee et al. (2018)222As seen in the figure, there is a possibility that the
required gap depth becomes deeper than $50\%$ for the PIM by Ataiee et al.
(2018) as well. However, we confirm that within our simulation parameters, the
effects are not important..

####  2.2.5 Gas accretion onto a protoplanet

The rapid gas accretion onto a protoplanet starts when a protoplanetary core
becomes massive enough so that the (quasi-)hydrostatic gas pressure can no
longer support the envelope against the gravity. Following Ikoma et al.
(2000), we assume this critical core mass for gas accretion is

$M_{{\rm core,crit}}\simeq 10\left[\left(\frac{\dot{M}_{{\rm
core}}}{10^{-6}{\rm M_{E}/yr}}\right)\left(\frac{\kappa}{1\,{\rm
cm^{2}/g}}\right)\right]^{s}M_{E},$ (40)

where $\dot{M}_{{\rm core}}$ is the core accretion rate, $s=0.2-0.3$, and
$\kappa$ is the grain opacity. For our simulations, we adopted $s=0.25$ and
$\kappa=1\,{\rm cm^{2}/g}$, respectively. The equation assumes that the
envelope is heated by the planetesimal accretion, and it is not immediately
clear whether pebble accretion would provide a comparable heating. Recently,
Ogihara & Hori (2020) explored this topic by using the 1D hydrostatic model
with pebble accretion rates over $10^{-12}-10^{-2}\,M_{E}/{\rm yr}$. They
obtained the fitting function: $M_{{\rm
crit}}=13\left(\frac{\dot{M}_{F}}{10^{-6}{\rm
M_{E}/yr}}\right)^{0.23}\,M_{E}$, which is consistent with the equation
adopted here. We only became aware of their work while this manuscript was
under revision, and thus we adopted Equation 40 instead of theirs. The model
implies that the gas accretion starts once the core accretion stops either
because the PIM is reached or because the pebble flux runs out.

The gas accretion onto a protoplanet is limited both by how quickly a
protoplanet can accrete gas and by how quickly a disc can supply gas to the
protoplanet (e.g. Ida et al., 2018). The former condition is given by the
Kelvin-Helmholtz (KH) timescale $\tau_{{\rm KH}}$ as follows:

$\frac{dM_{{\rm KH}}}{dt}\simeq\frac{M_{p}}{\tau_{{\rm KH}}},$ (41)

where

$\tau_{{\rm
KH}}=10^{b}\left(\frac{M_{p}}{M_{\oplus}}\right)^{-c}\left(\frac{\kappa}{1\,{\rm
cm^{2}/g}}\right)\,{\rm yrs}.$ (42)

Parameters $b$ and $c$ are not well constrained. Ikoma et al. (2000) proposed
$\left(b,\,c\right)=\left(8,\,2.5\right)$ based on their numerical
simulations, while Ida & Lin (2004) suggested that the parameters could take
values up to $\left(b,\,c\right)=\left(10,\,3.5\right)$ depending on different
opacity tables, and they adopted $\left(b,\,c\right)=\left(9,\,3.0\right)$ for
their work. For this work, we adopted $\left(b,\,c\right)=\left(8,\,3\right),$
though we also tested other combinations. The KH timescale becomes shorter for
smaller $b$ and larger $c$, and shorter timescales lead to more massive final
giant planet masses. However, beyond a certain point, the fast KH timescale
does not help the growth any more because the gas accretion is limited by the
efficiency of gas supply by the disc.

Regarding the latter condition, the disc nominally supplies gas at a rate
$\dot{M}_{*}$ (see Equation 14). For an actively accreting protoplanet,
however, we need to take account of the reduced surface mass density effect as
proposed by Tanigawa & Tanaka (2016):

$\dot{M}_{{\rm TT16}}=\frac{\Sigma_{{\rm
gap}}}{\Sigma_{g}}\dot{M}_{*}\simeq\frac{1}{\left(1+0.04K\right)}\dot{M}_{*}.$
(43)

Putting these together, we adopted the following gas accretion rate as in Ida
et al. (2018):

$\dot{M}_{{\rm gas}}\simeq{\rm min}\left[\dot{M}_{{\rm
KH}},\;\dot{M}_{*},\;\dot{M}_{{\rm TT16}}\right].$ (44)

###  2.3 Initial conditions

| $M_{d,0}\,(M_{\odot})$ | $t_{{\rm diff}}\,({\rm Myr})$ | $\alpha_{{\rm acc}}$ | $M_{{\rm peb,tot}}\,(M_{\oplus})$
---|---|---|---|---
Disc 1 | $0.2$ | $0.1$ | $\sim 7.4\times 10^{-2}$ | $194$
Disc 2 | $0.06$ | $0.1$ | $\sim 7.4\times 10^{-2}$ | $58.2$
Disc 3 | $0.02$ | $0.1$ | $\sim 7.4\times 10^{-2}$ | $19.4$
Disc 4 | $0.2$ | $1$ | $\sim 7.4\times 10^{-3}$ | $446$
Disc 5 | $0.06$ | $1$ | $\sim 7.4\times 10^{-3}$ | $134$
Disc 6 | $0.02$ | $1$ | $\sim 7.4\times 10^{-3}$ | $44.5$
Disc 7 | $0.2$ | $10$ | $\sim 7.4\times 10^{-4}$ | $597$
Disc 8 | $0.06$ | $10$ | $\sim 7.4\times 10^{-4}$ | $179$
Table 1: Disc models 1-8. Columns 2, 3, and 4 show the initial disc mass, the
disc diffusion timescale, and the corresponding viscosity $\alpha_{{\rm
acc}}$, respectively. The last column is the total pebble mass throughout the
simulations for ${\rm[Fe/H]}=0.0$. For all the simulations, the turbulent
viscosity $\alpha_{{\rm turb}}=10^{-4}$ is used, which ensures $\alpha_{{\rm
turb}}\ll\alpha_{{\rm acc}}$ (however, see Section 3.2.1 for the cases using
other values).

In this work, we considered eight different disc models around a Sun-like star
as shown in Table 1, which lead to the stellar mass accretion rates and the
corresponding pebble accretion rates as in Figure 2. For each of these discs,
we assumed five different stellar metallicities of
${\rm[Fe/H]}=-0.5,\,-0.3,\,0.0,\,0.3,\,{\rm and\,}0.5$. As stated earlier, we
assume that the global angular momentum transfer is carried out mainly by the
disc wind with $\alpha_{{\rm acc}}$ shown in the table, while the local disc-
planet interactions are controlled by the disc’s turbulent viscosity constant
$\alpha_{{\rm turb}}=10^{-4}$ (see Section 3.2.1 for other values).

Figure 4 shows the outcome of __in situ__ planet formation (i.e. no migration)
in Disc 5 with the solar metallicity, starting with a protoplanetary core mass
of $10^{-4}M_{\oplus}$. This initial mass is comparable to a Ceres mass, which
is a characteristic mass for a planetesimal formed via streaming instability
(Johansen et al., 2015; Simon et al., 2016). Differently from simulations
presented in Section 3, the effect of the snowline is ignored here for
simplicity. The top left panel compares the Kelvin-Helmholtz gas accretion
timescale for the PIM core (orange) with the core formation timescale (i.e.
the total time required for a core to reach the PIM, blue). The total in situ
formation timescale of a gas giant planet will be comparable to the sum of
these two timescales. For example, although the core accretion time is the
shortest around 1 au, the PIM core there is too small to accrete gas
efficiently. Thus, in our model, it is likely difficult to form a gas giant in
situ at around 1 au, and the gas giant planet formation timescale becomes
shortest around several au.

The right panel shows the total mass of a protoplanet reached after 100 Myr
across the disc (orange), along with the total mass at different times of the
growth calculations (blue curves). The solid green curve shows the PIM, while
the dotted green curve shows the final core masses in the regions where the
planetary masses never reach the PIMs. The figure indicates that protoplanets
achieve the PIMs in most parts of the disc within a few au by $\sim
10^{5}\,$yr and within $\sim 10\,$au by $\sim 10^{6}\,$yr. In this disc model,
planetary growth is inefficient beyond $\sim 10\,$au, because the accretion
cross-section and thus the pebble accretion efficiency sharply decreases for a
low-mass protoplanet in the outer disc. The corresponding pebble accretion
efficiencies are shown in the bottom left panel. Although this figure is made
based on the accretion efficiency by Ida et al. (2016), the trend is similar
with the efficiency by Ormel & Liu (2018).

We compiled similar figures for other discs and found that gas giant planets
do not form in situ in our disc models beyond $\sim 20\,$au, if the initial
core mass is comparable to that of Ceres. If typical planetesimals formed are
about the size of Ceres (e.g. Johansen et al., 2015; Simon et al., 2016)
throughout the disc, planetesimals have to grow via planetesimal-planetesimal
collisions or via some other manner to reach $\sim 10^{-2}\,M_{\oplus}$ to
form gas giants via pebble accretion beyond $\sim 20\,$au. The trend shown
here is largely consistent with Johansen & Lambrechts (2017), where growing a
core from $10^{-5}\,M_{\oplus}$ to $10^{-1}\,M_{\oplus}$ via pebble accretion
takes a few Myr at $30\,$au.

However, typical planetesimal sizes may be different depending on the disc
environment. We can roughly estimate the planetesimal mass formed via
streaming instability by considering the gravitational instability in the dust
layer:

$\displaystyle M_{{\rm plsml,\,init}}$ $\displaystyle\sim$
$\displaystyle\rho_{{\rm Roche}}\cdot\frac{4\pi}{3}H_{p}^{3}\sim
3M_{*}\left(\frac{H_{p}}{H_{g}}\right)^{3}\left(\frac{H_{g}}{a}\right)^{3}$
$\displaystyle\simeq$ $\displaystyle
13.8\,M_{\oplus}\left(1+\frac{\tau_{s}}{\alpha_{{\rm
turb}}}\right)^{-3/2}L_{*0}^{3/7}M_{*0}^{-5/7}\left(\frac{a}{{\rm
au}}\right)^{6/7}.$

To get the final expression, we used the Roche density $\rho_{{\rm
Roche}}\sim\frac{9}{4\pi}\frac{M_{*}}{a^{3}}$, Equation 30 for $H_{p}/H_{g}$,
and the disc aspect ratio of the irradiation region (Ida et al., 2016):

$\left(\frac{H_{g}}{a}\right)_{{\rm irr}}=0.024\left(\frac{T_{2}}{150\,{\rm
K}}\right)^{1/2}L_{*0}^{1/7}M_{*0}^{-4/7}\left(\frac{a}{{\rm
au}}\right)^{2/7}.$ (46)

With Equation 2.3, the initial planetesimal mass becomes $10^{-4}\,M_{\oplus}$
at $\sim 2.8\,$au and increases to $10^{-2}\,M_{\oplus}$ at $\sim 83\,$au in
Disc 5 (i.e. the same disc as in Figure 4), which makes the giant planet
formation possible out to $\sim 40\,$au. For longer-lived, more massive, or
more metal-rich discs, giant planets can form out to or beyond $100\,$au. In
this paper, we choose conservative initial conditions and do not place the
planetary cores beyond 20 au.

Exercises similar to that of Figure 4 also suggest that the maximum giant
planet mass achievable in our simulations is a few to several thousand
$M_{\oplus}$. The upper mass limit is determined partly by the planet’s
efficiency of gas accretion and partly by the disc’s capability of providing
gas (see Section 2.2.5). It may be possible to achieve $\sim
10^{4}\,M_{\oplus}$ with an extreme choice of parameters such as $(b,c)=(7,3)$
in a highly turbulent disc with $\alpha_{{\rm turb}}=10^{-2}$ or in a massive,
long-lived disc (e.g. $0.2\,M_{\odot}$ with $t_{{\rm diff}}\sim 10\,$Myr).
However, the former tends to lose planetary cores to the star due to efficient
migration, while the latter-type of discs may be relatively rare. Taking a
more extreme set of parameters would not help the further mass growth because
the gas accretion is limited by the disc’s gas supply as seen in Equation 44.
By assuming the initial masses of protoplanets with Equation 2.3, the maximum
mass becomes $\sim 50\,M_{J}\sim 1.6\times 10^{4}\,M_{\oplus}$ for massive,
metal-rich, long-lived discs. We discuss this issue further in Section 4.3.

We ran $240$ multiple-planet-core simulations as well as about a hundred more
single-planet core simulations. In all of these runs, we assumed a Sun-like
star. For the main simulations, we have ten cores per system with the initial
mass of $10^{-2}\,M_{\oplus}\sim 100\,M_{{\rm Ceres}}$ over $0.5-15\,$au
333The number of cores is chosen arbitrarily. However, Bitsch et al. (2020)
recently studied pebble accretion with 15, 30, and 60 cores across $3-17\,$au
and found that the number of cores generally has little effect on the final
outcome of planetary distributions. The exception is the case when the
eccentricity and inclination damping are very slow so that the higher number
of cores leads to more random orbits. Since the number of our cores is lower
and they are distributed over a wider orbital range compared to their work,
our initial conditions are less prone to dynamical instabilities compared to
theirs.. The inner edge of the disc is set at $0.1\,$au for all of our
simulations. This is rather arbitrary, but is adopted partly because the
timesteps required to resolve an orbit become too small within this radius to
run a simulation for $100\,$Myr, and partly because tidal interactions with
the star also become important there, and we do not explicitly include this
effect in our code (however, see Section 3.1.2). Each core has a randomly
assigned initial eccentricity of $e=\left[0,\,0.01\right]$, the corresponding
inclination of $i\,{\rm(rad)}=0.5e$, and randomly chosen phase angles. Planets
are also considered ‘removed’ beyond $1000\,$au in our simulations. In
summary, we simulated the growth and orbital evolution of the ten cores in the
eight different disc models shown in Table 1 with five different stellar
metallicities. Six simulations are run for each combination of parameters,
which makes the total number of main simulations $240$.

Besides the effects of disc models and stellar metallicities, we have also
explored the effects of the Kelvin-Helmholtz timescales, the turbulent
viscosity $\alpha_{{\rm turb}}$, and the pebble accretion efficiency. We
discuss some of these in Section 3.2.1 for single-planet simulations. For main
simulations, these parameters are set as $(b,\,c)=(8,\,3)$, $\alpha_{{\rm
turb}}=10^{-4}$, and $\epsilon=\epsilon_{{\rm IGM16}}$, respectively.

Figure 4: Top left panel compares the total time to form the protoplanetary
core with the pebble isolation mass (orange, $\tau_{{\rm peb}}$) with the
Kelvin-Helmholtz gas accretion timescale for the PIM core (blue, $\tau_{{\rm
KH}}$). The right panel shows the planetary mass growth at different times
starting from $10^{-4}\,M_{E}$ in Disc 5\. The solid green line shows the PIM,
while the dotted green line represents the core masses at 100 Myr (i.e. at the
end of the simulations) for the cores not reaching PIMs. The orange line shows
the final total mass. Here, the effects of the snow line are ignored. The
bottom left panel shows the corresponding pebble accretion efficiencies by Ida
et al. (2016) at different times as a core grows from $10^{-4}\,M_{E}$.

##  3 Results

In Section 3.1, we compare the results of main numerical simulations of
multiple-planet cores with observed systems, by focusing on the effects of
disc parameters such as mass, dissipation timescale, and metallicity. In
Section 3.2.1, we discuss the effects of the turbulent viscosity alpha
$\alpha_{{\rm turb}}$ as well as the pebble accretion efficiency $\epsilon$
through single-planet simulations. Throughout this paper, we use the radial-
velocity detected, confirmed extrasolar planets obtained from the NASA
Exoplanet Archive444https://exoplanetarchive.ipac.caltech.edu as ‘observed
planets’. We do not include transit planets for comparison with simulations
unless we state otherwise, because the transit surveys have a large bias
towards close-in planets due to the detectability dependence of $\propto
a^{-1}$.

###  3.1 The effects of disc mass, dissipation timescale, and stellar
metallicity

First, we present the overall results of simulations in Section 3.1.1 and show
that the overall trends of distributions of giant planets are well-reproduced.
Then we will focus on formation of different types of giant planets in
Sections 3.1.2 and 3.1.3.

####  3.1.1 Overall results of multiple protoplanetary-core simulations

In this paper, we define planetary bodies as having a mass $\geq
0.1\,M_{\oplus}$, because this is about the Mars mass and also roughly the
smallest mass of observed exoplanets. We also divided planets into two groups
for simplicity: giant planets ($\geq 0.1\,M_{J}\sim 30\,M_{E}$) and low-mass
planets ($<0.1\,M_{J}$). This division is rather arbitrary, but it is
motivated by a mass above which the semi-major axis mass distribution of
observed exoplanets appears to change.

Figure 5 shows the overall results of our simulations compared to the observed
planetary systems. The left, middle, and right panels correspond to $M_{p}-e$,
$a-e$, and $a-M_{p}$ distributions, respectively, for all the simulated
planets (bottom panels), observable simulated planets (middle panels), and
observed planets (top panels). For the ‘observable’ planets, we applied the RV
detection limit of $1\,$m/s and $a\leq 10$ au. As seen in the figure, the
simulated planets reproduce the overall trends of $M_{p}-e$, $a-e$, and
$a-M_{p}$ distributions of extrasolar giant planets well. We note that the
disc inner edge of our simulations is set at 0.1 au and the code does not
include the tidal evolution effects. Therefore, the cluster of simulated
planets near 0.1 au and their slightly elevated eccentricities compared to
observed systems are not surprising.

As seen in the $a-M_{p}$ distribution, our simulations successfully generate
all kinds of giant planets including hot Jupiters (HJs, $a\lesssim 0.1\,$au),
warm Jupiters (WJs, $0.1\,{\rm au}\,\lesssim a\lesssim 1\,{\rm au}$), cold
Jupiters (CJs, $1\,{\rm au}\,\lesssim a\lesssim 20\,{\rm au}$), and super-cold
Jupiters (SCJs, $a\gtrsim 20\,$au). As already shown by Ida et al. (2018), CJs
are formed successfully since the new type II migration is slower than the
classical one and since the low disc turbulence $\alpha_{{\rm turb}}=10^{-4}$
is assumed in the two-$\alpha$ disc model. We discuss the formation of HJs and
SCJs further in Sections 3.1.2 and 11, respectively.

As in Matsumura et al. (2017), we still have trouble reproducing the orbital
properties of low-mass planets within $\sim 1\,$au, because many
protoplanetary cores are lost to the star since type I migration is too
efficient and the disc edge is not sharp enough to retain them efficiently.
This may not bee too surprising because, although we have improved our code
significantly, as seen in Section 2, none of those mechanisms help to slow
type I migration. Lambrechts et al. (2019) and Izidoro et al. (2019)
successfully produced rocky super-Earth systems over $0.1-1\,$au, where we
lost many such planets due to rapid migration. This is partly because their
disc ages are much shorter than ours and partly because their planets are
trapped in mean motion resonances more efficiently at the disc edge since they
terminate planet migration there. In our simulations, some gas discs’
lifetimes are as short as a few Myr like their simulations, but others are 10
Myr or longer. The long lifetimes are motivated by recent observations as
stated before, but it also leads to the loss of the majority of low-mass
planets as expected.

Having this bias of low-mass planets in mind, our simulations suggest that
there may be a group of eccentric low-mass planets in the outer disc ($\gtrsim
10\,$au). In our simulations, they are coming from Discs 2 and 3, which are
not too massive ($0.06\,M_{\odot}$ and $0.02\,M_{\odot}$) and have a short
dissipation timescale of $0.1\,$Myr. As discussed further in the next sub-
section, these are protoplanetary cores which did not have enough time to
become giant planets. Since we reproduce the overall trends of distributions
of giant planets well while those of low-mass planets are poorer due to the
uncertainties in type I migration, we focus on the outcomes of giant planets
for the rest of this paper.

Figure 5: Comparison of distributions of parameters for the observed
exoplanets and left, middle, and right panels are $M_{p}-e$, $a-e$, and
$a-M_{p}$ distributions, respectively. The top panels show observed data from
the RV-detected planets. The bottom panels show all the simulated planets at
the end of the simulations (100 Myr), while the middle panels show observable
planets with the RV detection limit of $1\,$m/s and $a\leq 10$ au. The black
dashed line shown on $a-M_{p}$ panels corresponds to $1\,$m/s limit. The
shaded areas in $a-e$ and $a-M_{p}$ distributions indicate that the inner disc
edge of our model is at $0.1\,$au, and thus we do not intend to reproduce
planet distributions there. Circles and crosses represent giant ($\gtrsim
0.1\,M_{J}$) and low-mass ($<0.1\,M_{J}$) planets, respectively, and the red,
orange, green, blue, and purple colours correspond to stellar metallicities of
-0.5, -0.3, 0.0, 0.3, and 0.5, respectively. Some planets are clustered around
0.1 au because the disc’s inner edge is set there.

Although we compare our simulations with observed giant planets, it is not our
intention to exactly reproduce the distributions of orbital and physical
parameters of extrasolar planets. Such a comparison does not make sense for
this study since it would force us to make assumptions, for example, that disc
models 1-8 exist with an equal probability and that the inner and outer disc
radii are the same for all the discs, which are unlikely to be true in
reality. Having said that, Figure 6 presents the comparison of distributions
of the semi-major axis, planetary mass, and eccentricity for observed (dashed
lines) and simulated (solid lines) giant planets. The left panels include all
the giant planets ($\geq 30\,M_{E}$) beyond $0.1\,$au, while the right panels
include only those with masses $30-3000\,M_{E}$ and semi-major axes
$0.1-10\,$au. Only the planets beyond $0.1\,$au are chosen because our inner
disc edge sits there, and because we do not include tidal evolution in our
simulations directly (however, see Section 4.2).

By eye, the agreement between observed and simulated distributions are good
for both planetary masses and eccentricities. In fact, for the mass
distribution in the left panel, the Kolmogorov-Smirnov (K-S) test could not
reject the null hypothesis with the statistics of $0.092$ and the P value of
$0.11$. The agreement is poor for semi-major axis distributions, which
indicates that we overproduced CJs compared to WJs. This may imply that we
need to understand type I migration better to produce cores of these giants
and/or that the low disc turbulence $\alpha_{{\rm turb}}$ is not appropriate
for all the systems. Nevertheless, the overall good agreement with observed
properties is encouraging.

Figure 6: Distributions of semi-major axis (top), planetary mass (middle), and
eccentricity (bottom). The solid and dashed lines correspond to simulated and
RV-detected planets, respectively. The left panels include all the giant
planets beyond $0.1\,$au, while the right ones are limited to orbital radii of
$0.1-10\,$au and masses up to $3000\,M_{E}$. The K-S test cannot reject the
null hypothesis for the mass distribution on the left with the statistics of
$0.092$ and the P value of $0.11$, but rejects the null hypothesis for all
other cases. Our simulations clearly produce more giant planets in the outer
region compared to observed systems. However, the agreements are not bad by
eye for both masses and eccentricities.

As seen in $M_{p}-e$ and $a-e$ distributions of Figure 5, as well as the
middle panels of Figure 6, our simulations also reproduce the wide range of
eccentricities observed for giant planets. Compared to the lack of highly
eccentric giant planets seen in Matsumura et al. (2017) and Bitsch et al.
(2019), this is a great improvement. We suspect that the success is partly due
to a range of disc parameters leading to a larger number of giant planets per
system than Matsumura et al. (2017), and partly due to low efficiencies of
eccentricity and inclination damping (see Section 4.7).

In Matsumura et al. (2017), the average number of giant planets formed per
system throughout the simulations was typically about two for no-migration
cases. In our new simulations, on the other hand, the average number of giant
planets can be much higher, as seen in Figure 7. The number of giant planets
decreases over time due to dynamical instabilities, which lead to the
population of giant planets with high eccentricities. Overall, the number of
giant planets tends to be higher when a disc has higher mass, higher
metallicity, and longer dissipation timescale. For Disc 3, which has the
lowest mass and the shortest dissipation timescale, giant planet formation is
a rare event, and it happens only for the highest metallicity case. The
average number of giant planets formed in high-mass and/or high-metallicity
discs is typically a few or more and decreases over time to $\sim 2$ in the
end. However, we note that the number of giant planets per system was also
very high in Bitsch et al. (2019) (typically about 5), and we discuss this
issue further in Section 4.7.

Figure 7: Evolution of the average number of giant planets (solid) and low-
mass planets (dotted) for different disc models. Red, orange, green, blue, and
purple correspond to stellar metallicities of -0.5, -0.3, 0.0, 0.3, and 0.5,
respectively.

####  3.1.2 Formation of hot Jupiters

There are three main pathways to form hot HJs: in situ formation (e.g. Batygin
et al., 2016), disc migration after the formation of a giant planet beyond the
snow line (e.g. Lin et al., 1996), and the tidal circularisation of a highly
eccentric orbit of a giant planet (e.g. Fabrycky & Tremaine, 2007; Nagasawa et
al., 2008; Naoz et al., 2011; Wu & Lithwick, 2011). Here, the in situ
formation includes both the cases where protoplanetary cores also form in situ
(e.g. Chiang & Laughlin, 2013; Boley et al., 2016) and the cases where cores
form further out, migrate to the inner disc, and then accrete gas there (e.g.
Coleman & Nelson, 2016b; Matsumura et al., 2017). In this section, we evaluate
these pathways by using our simulations and argue that differences in stellar
metallicities may naturally lead to different pathways.

Figure 8 shows two-panel plots of eight different disc models, where each disc
is run with five different stellar metallicities. For each disc model, the
bottom panel shows a histogram of giant planet formation timescales, while the
top panel shows how far from the initial orbital radii giant planets are
formed. As discussed in the previous section, we assume that a protoplanet
becomes a giant planet once its mass surpasses $0.1\,M_{J}$.

As seen in bottom panels of different disc models, almost all discs form giant
planets regularly except for Disc 3. For the disc models we considered, the
median formation timescales are short ($\sim 0.1-1\,$Myr) compared to an
often-quoted typical disc lifetime of $\sim 3\,$Myr. This confirms that pebble
accretion is an efficient way of forming giant planets within a disc’s
lifetime.

The top panels of different disc models in Figure 8 compare the ratio of the
‘formation’ semi-major axis to the initial semi-major axis: $a_{{\rm
form}}/a_{{\rm in}}$ with the formation timescales. Here, again, the formation
semi-major axis is where a planetary mass reaches $0.1\,M_{J}$. The red,
orange, green, blue, and purple symbols correspond to different metallicities
of ${\rm[Fe/H]}=-0.5$, $-0.3$, 0.0, 0.3, and 0.5, respectively. As expected,
giant planets form more quickly in higher metallicity environments within the
same disc model.

In discs with short dissipation timescales (Discs 1 to 3, $t_{{\rm
diff}}=0.1\,$Myr), we find that giant planets often form near their initial
locations. On the other hand, in discs with longer dissipation timescales
(Discs 4 to 9, $t_{{\rm diff}}=1-10\,$Myr), there are two types of giant
planets: some are forming near their initial locations (often cores beyond
several au), and others are forming __after__ some migration. In the case of
the latter, the formation radii tend to be about 1 to 2 orders of magnitude
smaller than the initial orbital radii. In Discs 4-9, when metallicities are
low ${\rm[Fe/H]}\lesssim-0.3$, cores tend to migrate and then accrete gas to
become either WJs or HJs. On the other hand, when metallicities are high
${\rm[Fe/H]}\gtrsim 0.0$, outer cores tend to become CJs, while inner cores
migrate to form WJs and/or HJs. However, in these cases, the inner cores are
often removed via planet-planet interactions as the gas disc dissipates, and
the final systems include either only CJs or HJs along with CJs. Thus, in our
simulations, when HJs form in situ after their cores migrated to the inner
disc, they tend to be in a low-metallicity disc, though they can also be
formed in a higher metallicity disc.

Figure 8: Top: Ratio of the semi-major axis of a giant planet when it became
giant (i.e. mass reaches $0.1\,M_{J}$) with respect to the initial semi-major
axis for different disc models. The red, orange, green, blue, and purple
symbols correspond to stellar metallicities of -0.5, -0.3, 0.0, 0.3, and 0.5,
respectively. Bottom: Formation timescales of giant planets. The vertical
dashed line indicates the median value.

Figure 9 shows the semi-major axis mass distribution of giant planets for each
disc model. Combined with the results of Figure 8, this figure shows two types
of outcomes of giant planet formation. In rapidly dissipating or low-mass
discs (Discs 1, 2, 3, and 6), low-mass (typically less than a few Jupiter
masses) and cold (orbiting beyond $\sim 1\,$au) giant planets are often
produced. In these discs, lower metallicities lead to lower mass giant planets
as expected, and the range of masses across different metallicities is about
an order of magnitude. In more slowly dissipating, moderate-to-high-mass discs
(Discs 4, 5, 7, and 8), on the other hand, giant planets have a much wider
distribution of orbital radii (two to three orders of magnitude), as expected
from Figure 8. Planetary masses per disc model are largely comparable to one
another and weakly depend on metallicities (see Section 4.6), but cores
migrated further tend to have lower masses.

This $a-M_{p}$ relation of mass increasing with orbital radius is likely to
arise from the reduced gas accretion rate onto the gap-opening core
$\dot{M}_{{\rm TT16}}$ in Equation 44. Assuming that the gas accretion
timescale can be written as $\tau_{{\rm gas}}\sim\frac{M_{p}}{\dot{M}_{{\rm
TT16}}}$, the planetary mass has the dependence of: $M_{p}\propto
h_{g}^{\frac{5}{3}}$. Since the disc aspect ratios in viscous and irradiation
regions in our disc model have $h_{g,{\rm vis}}\propto\left(\frac{a}{{\rm
au}}\right)^{\frac{1}{20}}$ and $h_{g,{\rm irr}}\propto\left(\frac{a}{{\rm
au}}\right)^{\frac{2}{7}}$, respectively (see Ida et al., 2016), the expected
mass dependences on semi-major axes are $M_{p,{\rm
vis}}\propto\left(\frac{a}{{\rm au}}\right)^{\frac{1}{12}}$ and $M_{p,{\rm
irr}}\propto\left(\frac{a}{{\rm au}}\right)^{\frac{10}{21}}$, respectively.
For Discs 4, 5, 7, and 8 in Figure 9, these trends are shown in dashed and
dotted lines, respectively, and the agreement is very good for these cases.

Figure 9: Semi-major axis mass distribution of giant planets at the end of all
the simulations for different disc models. The red, orange, green, blue, and
purple colours correspond to stellar metallicities of -0.5, -0.3, 0.0, 0.3,
and 0.5, respectively. The dashed and dotted lines in Discs 4, 5, 7, and 8 are
expected mass trends in viscous and irradiation regions, respectively (see
text).

Our simulations do not explicitly include tidal interaction effects between
the central star and planets, so we cannot directly measure the production
rate of a HJ via tidal circularisation. However, we selected planets that
could have experienced significant tidal evolution and applied the weak-
friction tidal interaction model (e.g. Hut, 1981). For giant planets that have
a pericentre radius of $a\left(1-e\right)\leq 0.1\,$au and a semi-major axis
$a\geq 0.2\,$au (i.e. an eccentricity of $e\geq 0.5$) at any time during the
simulation, tidal evolution is calculated. The results are plotted in Figure
10. As in Matsumura et al. (2010), we scaled a modified tidal quality factor
of

$Q^{\prime}=Q_{0}^{\prime}\frac{n_{0}}{n},$ (47)

where $Q^{\prime}\equiv 1.5Q/k_{2}$ is a modified tidal quality factor, $Q$ is
a tidal quality factor, $k_{2}$ is the Love number of degree 2, $n$ is the
mean motion, and the subscript $0$ represents their initial values. In our
calculations, we assumed the initial stellar and planetary modified tidal
quality factors to be $Q_{*0}^{\prime}=10^{6}$ and $Q_{p0}^{\prime}=10^{5}$,
respectively. Since a typical age of a planet hosting star ranges from $1$ to
$10\,$Gyr, we have run the tidal evolution simulations for $3\,$Gyr.

The left panel of Figure 10 shows the semi-major axis eccentricity
distribution of these planets before (open circle) and after (filled circle)
tidal evolution, where different colours correspond to different stellar
metallicities. As expected from the higher planet formation efficiency seen in
high metallicity environments (see Figure 7, for example), most planets that
achieved high enough eccentricities in the first place are in the high-
metallicity discs with ${\rm[Fe/H]\gtrsim 0.0}$.

Since we do not know what kinds of discs or disc evolution paths dominate in
reality, the fraction of HJ systems with respect to all the simulations cannot
be compared with the observed fraction of HJ systems. However, we can attempt
to estimate the fractions of __giant-planet__ systems with HJs. About $10\,\%$
of Sun-like stars are known to have giant planets, while roughly $0.5-1\,\%$
of such stars have HJs (Winn, 2018). This implies that, roughly speaking,
$\sim 5-10\,\%$ of giant-planet hosting stars have at least one HJ. In our
simulations, out of $240$ runs, $175$ formed at least one giant planet at some
point during $100\,$Myr simulations. Out of these $175$ systems, $10$ formed
HJs that have semi-major axes less than $0.1\,$au after $3\,$Gyr. Thus
$10/175\sim 5.71\,\%$ of a giant-planet-forming system led to the formation of
at least one HJ. If we relax this condition to the final semi-major axis
within 0.2 au, the frequency becomes $38/175\sim 21.7\,\%$. These are
reasonable values compared to the observed fraction of HJ systems among
systems with giant planets: $\sim 5-10\,\%$. We note that the occurrence rate
depends on the choice of the initial distribution of disc parameters, and thus
the high occurrence rate suggested by the simulations does not necessarily
mean that our models overestimate the occurrence rate.

The right panel of Figure 10 shows the corresponding eccentricity-inclination
distribution of these planets before (open circles) and after (filled circles)
applying tidal evolution. We note that the colours here indicate the final
semi-major axes rather than the stellar metallicities. The overall trend is
similar to what Nagasawa & Ida (2011) pointed out for observed exoplanets —
planets with low eccentricities have a wide range of inclinations, while those
with moderate eccentricities have low inclinations. Our simulations indicate
that this trend breaks down at larger orbital radii, because the orbits of
some planets cannot be circularised quickly enough.

From the figure, it is also clear that high-inclination planets initially had
high eccentricities as well. This outcome of tidal evolution reflects the fact
that the orbital circularisation is largely controlled by tidal dissipation
inside the planet, while the spin alignment is controlled by tidal dissipation
inside the star (i.e. the stellar spin aligns with the orbit normal rather
than the other way around). Since the former takes place faster than the
latter, the orbital eccentricity becomes small, while the orbital inclination
does not change dramatically.

These outcomes lead to an expectation about the effects of stellar
metallicities on HJ formation pathways — around low-metallicity stars, HJs
tend to form via the migration of cores followed by in situ gas accretion,
while around high-metallicity stars, they can form either via migration or via
the tidal circularisation of eccentric orbits. One way to distinguish these
two cases is to check the relation between the orbital inclination and the
stellar metallicity. The HJs formed in situ or via migration are more likely
to have low inclinations, while the tidally circularised planets are expected
to have a wide range of inclinations. We discuss this matter further in
Section 4.2.

Figure 10: Left: Semi-major axis eccentricity distribution of planets before
(open circle) and after (filled circle) tidal evolution. Different colours
correspond to the stellar metallicity. Highly eccentric planets tend to belong
to higher metallicity discs. Right: the corresponding eccentricity-inclination
distribution. However, the colour does not represent the metallicity, but the
final semi-major axis. Closer in planets tend to have either a low
eccentricity and a range of inclinations, or a moderate eccentricity with a
low inclination.

####  3.1.3 Formation of super-cold Jupiters

The observations show that some giant planets have very wide orbital radii,
ranging from a few tens of au to a few thousand au. In the pebble accretion
model, we expect either (1) an SCJ forms in situ, or (2) a gas giant forms
within $\sim 20\,$au and then gets scattered into the outer disc. Again, the
in situ formation here includes both the case where an SCJ directly forms in
the outer disc, and the case where a protoplanetary core scattered into the
outer disc accretes gas to become a giant planet there (Kikuchi et al., 2014).
Since we did not place cores beyond $20\,$au, SCJs did not form directly in
the outer disc in our simulations.

As seen in Figure 8, most giant planets in our simulations ‘form’ (i.e. mass
exceeds $0.1\,M_{J}$ in our definition) near or inside the initial orbital
radii. Therefore, most SCJs in our simulations are generated via scenario (2)
(i.e. formed and then scattered), and not via scenario (1). We note, however,
that our simulations preclude the scenario of Kikuchi et al. (2014). The key
to forming SCJs in Kikuchi et al. (2014) was that the orbits of scattered
cores get circularised efficiently via the accretion of high angular momentum
gas in the outer disc. Since our disc only stretches out to 100 au, even when
there was a scattered core, the core would repeatedly encounter the larger
planet that scattered the core in the first place, and an SCJ does not form
efficiently.

When an SCJ is formed via scenario (2), our simulations indicate the
following. First, a stellar metallicity tends to be high $\left[{\rm
Fe/H}\right]\gtrsim 0.0$ in a system with a SCJ. As seen in the left panel of
Figure 11, all but one giant planet beyond 20 au are from such systems. One
giant planet is at around 30 au in the system with $\left[{\rm
Fe/H}\right]=-0.5$, but its mass is relatively low and $\sim 40\,M_{\oplus}$.
This trend is not surprising, because higher metallicity discs tend to form
multiple giant planets without much migration (see Figure 8) and thus are
prone to dynamical instabilities that can scatter giant planets outward.

Second, the orbit of an SCJ tends to be eccentric and out to a few hundred au.
It is difficult to scatter a giant planet beyond that and retain it when it is
formed within $\sim 20\,$au. The tendency towards a non-zero eccentricity is a
direct outcome of the scattering event and thus is not surprising. In our
simulations, the eccentricity of SCJs ranges over $e\sim 0.2-0.9$ with the
mean value of 0.56.

Third, an SCJ tends to be accompanied by another giant planet. Our simulations
show $\sim 78.6\,\%$ of SCJs (22/28) are in two- or three-giant-planet
systems, while $\sim 21.4\,\%$ (6/28) are in single-planet systems. The right
panel of Figure 11 shows the semi-major axis ratio of a companion planet to
the furthermost SCJ as well as their mass ratio. In multiple-planet systems,
the furthermost SCJ tends to have a less massive companion interior to its
orbit ($15/22\sim 68.2\,\%$) rather than a more massive one ($7/22\sim
31.8\,\%$). Thus, the scattering event tends to place a more massive planet
outward and a lower-mass companion inward in our simulations.

The preference for the inner, lower-mass companion is opposite to what has
been observed in previous studies (e.g. Marzari & Weidenschilling, 2002;
Nagasawa et al., 2008; Ida et al., 2013), where a more massive planet was
scattered inward with a higher probability of $>80\,\%$. However, all of these
studies started with the configuration where the inner planet is the most
massive one $(2\,M_{J},\,M_{J},\,M_{J})$. Since the higher mass planets tend
to stay closer to the initial location upon scattering (Chatterjee et al.,
2008), it may not be surprising that the more massive planet stayed in the
innermost orbit in these studies. In our simulations, on the other hand, more
massive planets tend to form in the outer part of the disc (see Figure 4, for
example). Therefore, it is not surprising that there is a preference towards
the outermost planets being more massive than the inner ones because of their
inertia and the difficulty of switching orbits. We discuss their formation and
distributions further in Section 4.4.

Figure 11: Left: Semi-major axis mass distribution of giant planets at 100
Myr. The error bars represent the locations of pericentre and apocentre and
thus show how eccentric orbits are. As before, red, orange, green, blue, and
purple correspond to stellar metallicities of -0.5, -0.3, 0.0, 0.3, and 0.5,
respectively. The circles, squares, and crosses correspond to single-planet,
multiple-giant-planet, and multiple-planet but single giant planet systems,
respectively. Right: Mass and semi-major axis ratios for the furthermost SCJ
and its companion.

###  3.2 The effects of other parameters

In this sub-section, we show the effects of the turbulent viscosity
$\alpha_{{\rm turb}}$ as well as the pebble accretion efficiencies.

####  3.2.1 The effects of the turbulent viscosity alpha parameter
$\alpha_{{\rm turb}}$

For the bulk simulations, we assumed $\alpha_{{\rm turb}}=10^{-4}$. In this
sub-section, we investigate the dependence of pebble accretion simulations on
the turbulent alpha parameter by using single-planet-core simulations and by
comparing $\alpha_{{\rm turb}}=10^{-3}$ and $10^{-5}$ cases to the default
ones with $\alpha_{{\rm turb}}=10^{-4}$. We did not use all the disc models
but focused on representative cases of Discs 2, 5, and 8, for which the
initial disc mass is the same, $0.06\,M_{\odot,}$ and the diffusion timescales
are 0.1, 1.0, and 10 Myr, respectively. The initial semi-major axes of
protoplanets are chosen to be $\sim 3$, 5, 7, and $10\,$au, which are all
beyond the snowline initially. This is a rather narrow range compared to $\sim
0.5-15\,$au in the bulk simulations, but they are chosen so that we can test
the likelihood of forming giant planets. The stellar metallicity for all the
simulations in this sub-section is the solar value of ${\rm[Fe/H]}=0.0$. Other
parameters used for the default case are $(b,c)=(8,3)$ for the Kelvin-
Helmholtz timescale, and the pebble accretion efficiency of $\epsilon_{{\rm
IGM16}}$.

Figure 12 shows the outcomes of single-core simulations with $\alpha_{{\rm
turb}}=10^{-3},\,10^{-4},$ and $10^{-5}$. From the cases with $\alpha_{{\rm
turb}}=10^{-4}$, we can confirm the trends described in previous sections:
Disc 2 forms rather low-mass CJs, while Discs 5 and 8 form fully grown giant
planets across a range of orbital radii.

For the cases of $\alpha_{{\rm turb}}=10^{-3}$, most planets migrate to the
edge of the disc ($0.1\,$au here), except for Disc 2, where relatively low-
mass giant planets form over a range of orbits. The loss of the majority of
protoplanetary cores to the disc’s inner edge is similar to Matsumura et al.
(2017), where we have also adopted $\alpha_{{\rm turb}}=10^{-3}$. However,
differently from Matsumura et al. (2017), the formation of relatively low-mass
CJs is possible here largely due to the new migration prescription we adopted
(see Section 2.1). The results are consistent with what we would expect from
those in Section 3.1. For the higher viscosity $\alpha_{{\rm turb}}$, the PIM
is higher, and thus planet formation becomes too slow (compared to type I
migration) to form giant planets in Discs 5 and 8. The giant planets are
formed only in the most rapidly dissipating disc (Disc 2), which provides a
high mass flux for a short period of time, but their masses are low for giant
planets.

Figure 12: Outcome of single-planet formation simulations with the default
setting (i.e. $(b,c)=(8,3)$ for the Kelvin-Helmholtz timescale,
${\rm[Fe/H]}=0.0$, and the pebble accretion efficiency by IGM16) with
$\alpha_{{\rm turb}}=10^{-3}$ (left panels), $10^{-4}$ (middle panels), and
$10^{-5}$ (right panels). The top panels show the final semi-major axis of a
single planet for different initial semi-major axes and disc models. The
bottom cells marked as ‘Init’ correspond to the initial semi-major axes that
are common for all disc models. The bottom panels present the final planetary
masses in log scale. The bottom cells represent the initial protoplanetary
masses of $0.01M_{\oplus}$ for all the cases.

For the cases of $\alpha_{{\rm turb}}=10^{-5}$, the formation rate of giant
planets, especially CJs beyond 1 au, is even higher than the default case of
$\alpha_{{\rm turb}}=10^{-4}$ for Discs 5 and 8, while only super-Earths are
formed for Disc 2. The trends are again consistent with those seen in Section
3.1. For the lower $\alpha_{{\rm turb}}$ , both the PIM and the migration
transition mass are lower. Therefore, protoplanetary cores start gas accretion
at lower mass and also tend to switch from type I to type II migration earlier
and survive. In Disc 2, the cores are too small to become giant planets within
the disc’s lifetime, while in Discs 5 and 8, CJs can be formed over a wide
range of orbital radii.

####  3.2.2 The effects of the pebble accretion efficiency

Figure 13 shows the same as Figure 12, but with the pebble accretion
efficiency by Ormel & Liu (2018) $\epsilon_{{\rm OL18}}$ (see Equation 34)
instead of that by Ida et al. (2016) $\epsilon_{{\rm IGM16}}$ (see Equation
27). Despite the small differences between the two accretion efficiencies seen
in Figure 3, these two figures look dramatically different. There are a couple
of reasons here. First, a factor of a few difference in accretion efficiencies
actually leads to a large difference. For example, a factor of 3 difference
corresponds to the metallicity difference between ${\rm[Fe/H]=0.0}$ and
${\rm[Fe/H]=0.5}$. Second, the accretion efficiency by Ormel & Liu (2018) is
more sensitive to the inclination. For example, among cases shown in Figure
13, a protoplanet at $\sim 5\,$au consistently produces a much lower mass
planet compared to the neighbouring ones. This is because the randomly chosen
initial inclination for this core is $i\sim 0.26$ degree (i.e. $e\sim 0.01$)
as opposed to $i\lesssim 0.14$ degrees (i.e. $e\lesssim 0.005$) for the
others. Compared to the circular, coplanar case, the accretion efficiency
$\epsilon_{{\rm OL18}}$ decreases by more than an order of magnitude with
$e\sim 0.01$ and $i\sim 0.26$ degree, while the difference is about a factor
of a few with $e\lesssim 0.005$ and $i\lesssim 0.14$ degrees.

Besides the slow growth of a protoplanet at $\sim 5\,$au, the outcomes shown
in Figure 13 are consistent with Figure 12 and the expectations from Section
3.1. With $\alpha_{{\rm turb}}=10^{-5}$, giant planets are formed across a
range of disc radii for Discs 5 and 8. Compared to the corresponding cases in
Figure 12, although the final masses of giant planets are comparable, the
final semi-major axes are smaller in these cases. This is because it takes
longer to form planets with lower efficiencies, and planetary cores migrate
more before becoming giant planets in these cases. In Disc 2, only very low-
mass planets are formed, which is also consistent with this expectation.

With $\alpha_{{\rm turb}}=10^{-4}$, none of the cases with the solar
metallicity lead to giant planet formation. With $\alpha_{{\rm
turb}}=10^{-3}$, trends are similar to the cases in $\alpha_{{\rm
turb}}=10^{-4}$, but HJs are formed at the disc’s inner edges in Disc 5. Since
no giant planets are formed in the corresponding cases in Figure 12, these
cases may appear surprising. What happened here is that protoplanets grew
slower with the lower pebble accretion efficiencies, stayed longer beyond the
snow line, and obtained slightly more massive cores before migrating into the
inner disc region, where they become HJs by accreting gas.

Figure 13: Same as Figure 12, but with the pebble accretion efficiency of
Ormel & Liu (2018) rather than that of Ida et al. (2016).

For completion, we also ran a set of multiple-core simulations with
$\epsilon_{{\rm OL18}}$ and showed the results in Figure 14. Here, we only use
Discs 2, 5, and 8, but the other parameters are the same as those shown in
Section 3.1 with the default accretion efficiency $\epsilon_{{\rm IGM16}}$.
The overall distributions are similar to Figure 5 for all the parameters,
except that giant planets are mostly formed with high metallicities
${\rm[Fe/H]}\gtrsim 0.3$. This is consistent with the expectation from single-
core simulations shown above. The lack of very massive giant planets ($\gtrsim
2000\,M_{E}$) compared to Figure 5 is unlikely to be the real feature. As seen
in Figure 9, most massive giant planets are formed in most massive discs
($0.2\,M_{\odot}$ in our case), while all the discs used in this figure have
$0.06\,M_{\odot}$. Therefore, to have distributions similar to Figure 5 with
the pebble accretion efficiency $\epsilon_{{\rm OL18}}$ by Ormel & Liu (2018),
we would probably need to assume lower turbulence alpha $\alpha_{{\rm
turb}}<10^{-4}$ and very low inclinations for initial protoplanetary cores.

Figure 14: Same as Figure 5, but with the pebble accretion efficiency
$\epsilon_{{\rm OL18}}$ from Ormel & Liu (2018). Only Discs 2, 5, and 8 are
used here. The overall distributions are similar to Figure 5, but giant
planets are formed mostly with high metallicities ${\rm[Fe/H]}\gtrsim 0.3$
(see text for further discussion).

##  4 Discussion

In this section, we further discuss our results by focusing on several key
topics. First, we compare planet occurrence rates obtained from our
simulations with observations in Section 4.1. Then, we discuss formation of
HJs (Section 4.2) and SCJs (Section 4.4). We also discuss to what orbital
radii planets can be formed via pebble accretion and how massive planets can
become in Section 4.3, and discuss metal and gas mass fractions of planets in
Section 4.5. Finally, we compare our results with planetesimal accretion work
(Section 4.6) as well as other pebble accretion work (Section 4.7).

###  4.1 Planet occurrence rates

Figure 15: Top: Distribution of orbital periods of giant planets at the end of
the simulations (blue) and those of giant planets whose orbits are tidally
evolved for 3 Gyr (orange, see Section 3.1.2). Bottom: Corresponding
occurrence rates of giant planets as a function of orbital period. The blue
symbols only include simulated giant planets, while the orange ones include
both simulated and tidally evolved giant planets. The rates are calculated for
the number of planets per star per $\Delta\log P=0.23$ to be compared with
Figure 3 of Winn (2018). The black slopes are occurrence rates estimated from
observed planets. The solid, dashed, dotted, and dot-dashed lines correspond
to $dn/d\ln P\propto P^{0.26\pm 0.1}$ for giant planets with $>0.3\,M_{J}$
within $2000\,$days from Cumming et al. (2008), $dn/d\ln P\propto
P^{-1.45_{-0.44}^{+0.51}}$ for giant planets $5-5000\,$au from Baron et al.
(2019), an asymmetric broken power law of $dn/d\ln P\propto P^{0.53\pm 0.09}$
and $dn/d\ln P\propto P^{-1.22\pm 0.47}$ with a break period of $P_{{\rm
break}}=1717\pm 432\,$days for giant planets with $0.1-20\,M_{J}$ with periods
$\sim 1-10^{4}\,$days from Fernandes et al. (2019), and $dn/d\ln P\propto
P^{-0.453}$ for giant planets with $5-13\,M_{J}$ with periods $\sim 1.16\times
10^{4}-3.65\times 10^{5}\,$days ($\sim 10-100\,$au) from Nielsen et al.
(2019), respectively.

The top panel of Figure 15 shows the overall distribution of orbital periods
of giant planets from our work. The blue and orange histograms correspond to
the distributions of simulated planets and tidally evolved planets from
Section 3.1.2, respectively. The strong peak seen around $10\,$days comes from
the inner disc edge being set at $0.1\,$au in our simulations, and the actual
distribution around this period should be broader depending on the locations
of the inner disc edges. Otherwise, the overall distribution shows that CJs
near $10^{3}-10^{4}\,$days are most abundant, and the number decreases for
planets interior and exterior to these orbits.

The bottom panel of Figure 15 shows the corresponding occurrence rates
compared with observed slopes. Following the approach by Winn (2018), we write
the occurrence rate as: $\frac{\Delta n}{\Delta\log P}$, where $n=N_{{\rm
pl}}/N_{*}$ is the ratio of the number of giant planets to that of the host
star, while $\Delta\log P=0.23$ is chosen so that the comparison with Figure 3
of Winn (2018) is easier. The blue and orange symbols correspond to the
occurrence rates calculated for simulated planets only and simulated and
tidally evolved planets together, respectively. The total number of simulated
giant planets is $N_{{\rm pl,\,sim}}=274,$ and the corresponding number of
stars is $N_{{\rm*,\,sim}}=240$, while we just assume the number of the
tidally evolved planets and that of stars are the same: $N_{{\rm
pl,\,tide}}=N_{{\rm*,\,tide}}=76$.

The lines shown in the figure indicate the slopes of $\frac{dn}{d\ln P}\propto
P^{-\beta}$ obtained from four studies of observed systems (Cumming et al.,
2008; Baron et al., 2019; Fernandes et al., 2019; Nielsen et al., 2019).
Cumming et al. (2008) used RV-detected planets around FGK stars and suggested
$\beta=0.26\pm 0.1$ for giant planets with masses $>0.3\,M_{J}$ and periods
$<2000\,$days, while Fernandes et al. (2019) investigated both RV and Kepler
studies and found a broken power-law function with $\beta=0.53\pm 0.09$ and
$\beta=-1.22\pm 0.47$ with the break period of $1717\pm\,432$days. The other
two studies focus on outer planets observed by direct imaging methods. Baron
et al. (2019) suggested $\beta\sim-1.45_{-0.44}^{+0.51}$ for giant planets
with masses $1-20\,M_{J}$ and orbital radii of $5-5000\,$au, while Nielsen et
al. (2019) found $\beta\sim-0.453$ for planets with $5-13\,M_{J}$ and orbital
radii of $10-100\,$au. We note that the scaling factors used for these lines
are chosen to make the comparison with our data easier, and the actual
observed occurrence rates of these planets are much lower since $N_{*}$ is
larger and the detection probability is not $100\,\%$ (e.g. Winn, 2018). Here,
we only attempt to compare the slopes to see the relative occurrence rates of
different types of giant planets.

The overall trends of observed slopes agree well with those seen for simulated
planets. The occurrence rates estimated from our simulations increase with
orbital periods out to $\sim 750\,{\rm days}\,(\sim 1.6\,{\rm au})$, decrease
with periods slightly beyond $\sim 10^{4}\,{\rm days}\,(\sim 9\,{\rm au})$,
and have a broad peak between these periods. As seen in the figure, the inner
ascending trend is well explained by the slopes estimated by Cumming et al.
(2008) and Fernandes et al. (2019). The fact that the occurrence rates
increase beyond $\sim 1\,$au is also consistent with estimates both from RV
surveys (Mayor et al., 2011) and the Kepler transit survey (Santerne et al.,
2016). The outer descending trend is also broadly consistent with the
estimates by Nielsen et al. (2019), Baron et al. (2019), and Fernandes et al.
(2019), though details are different. Baron et al. (2019) proposed a much
steeper slope compared to our predictions, and Fernandes et al. (2019)
estimated the break point to be at a shorter orbital period. The overall shape
is also consistent with Bryan et al. (2016), who suggested a broad peak in the
distribution between 3 and 10 au and positive and negative $\beta$ values
interior and exterior to these regions. However, where the occurrence rates
start decreasing is still unclear. The future observations such as Gaia may
provide a population of planets that help to constrain this trend better.

The total occurrence rates are estimated to be $\sim 10\,\%$ for giant planets
with $\gtrsim 0.3\,M_{J}$ and $\lesssim 3-4.6\,$au (Cumming et al., 2008;
Mayor et al., 2011), $\sim 52.4\,\%$ for planets with $1-20\,M_{J}$ and
$5-20\,$au (Bryan et al., 2016), and $2.61_{-1.00}^{+6.97}\,\%$ for giant
planets with $1-20\,M_{J}$ and $20-5000\,$au. The trends of these occurrence
rates are broadly consistent with those seen in our simulations.

Table 2 summarises the numbers of HJs, WJs, CJs, and SCJs formed in our
simulations. We note that the fractions shown here are determined with respect
to the total number of giant planets and are different from the occurrence
rates. We find that $>50\,\%$ of all the giant planets are CJs, and thus more
than half the giant planets stay near the formation region.

| HJs | WJs | CJs | SCJs | Total
---|---|---|---|---|---
|  $\leq 0.1\,$au |  $0.1-1\,$au |  $1-20\,$au |  $>20\,$au |
$N_{{\rm pl,\,sim}}$ | $33$ | $43$ | $174$ | $24$ | $274$
$\%$ | $12.0$ | $15.7$ | $63.5$ | $8.8$ | $100$
$N_{{\rm pl,\,sim}}+N_{{\rm pl,\,tide}}$ | $50$ | $96$ | $180$ | $24$ | $350$
$\%$ | $14.3$ | $27.4$ | $51.4$ | $6.9$ | $100$
Table 2: Numbers and fractions of HJs, WJs, CJs, and SCJs as defined in the
semi-major axis range shown. $N_{{\rm sim}}$ corresponds to the number of
simulated giant planets, while $N_{{\rm tide}}$ corresponds to the number of
giant planets for which tidal evolution is further calculated for 3 Gyr (see
Section 3.1.2).

###  4.2 Formation of hot Jupiters

Any HJ formation model needs to explain that (1) HJs tend not to have nearby
companions (e.g. Wright et al., 2009; Steffen et al., 2012; Huang et al.,
2016), and (2) planets around higher metallicity stars tend to have a wider
range of eccentricities (e.g. Dawson & Murray-Clay, 2013; Shabram et al.,
2016; Buchhave et al., 2018). Our simulations naturally explain both of these
trends.

When HJs are formed via the tidal circularisation of highly eccentric orbits,
the underlying assumption is that planetary systems have undergone some sort
of dynamical evolution that increased eccentricities, such as planet-planet
scattering (e.g. Ford & Rasio, 2008; Chatterjee et al., 2008; Jurić &
Tremaine, 2008), the secular chaos (e.g. Wu & Lithwick, 2011), or the Kozai-
Lidov mechanism (e.g. Fabrycky & Tremaine, 2007; Naoz, 2016). The evolution of
such planets often removes neighbouring planets, and the resulting HJs tend to
be alone and have a range of eccentricities and inclinations. In our
simulations, these HJs often belong to high-metallicity and/or massive discs,
which typically produce a few giant planets that are prone to dynamical
instabilities. On the other hand, when HJs form in situ (i.e. via the
migration of cores followed by gas accretion in our simulations), they often
belong to low-metallicity discs, which typically produce 1-2 giant planets.
Therefore, the resulting HJs tend to be alone and have low eccentricities and
inclinations due to damping by the disc. HJs can also be formed in situ or via
migration in high-metallicity discs. Such discs usually form multiple giant
planets initially, including HJs, WJs, and CJs, but HJs and/or WJs are often
removed via dynamical instabilities later on, leaving either (i) only CJs, or
(ii) HJs and CJs (also see Mustill et al., 2015). Therefore, even in these
cases, HJs usually have no neighbouring giant planets (but could be
accompanied by CJs further out).

There may be some difficulties for the scenario in which HJs are formed in
situ, in particular in low-metallicity discs. First, gas accretion may be
stalled near the central stars, because the gas envelope becomes nearly
isothermal and the timescale of replenishing the atmosphere with the disc gas
becomes shorter than the Kelvin-Helmholtz gas accretion timescale (e.g. Ormel
et al., 2015a, b; Ali-Dib et al., 2020). However, the gas accretion may still
be possible when the disc is heated by the stellar radiation alone and
temperature and entropy in the inner disc region is low (Ali-Dib et al.,
2020), which may be the case towards the end of the disc’s lifetime. Second,
the near-infrared disc fraction study by Yasui et al. (2010) suggested that
the disc fraction of the low metallicity clusters declines within 1 Myr as
opposed to several Myr for the solar-metallicity clusters. If the trend were
to be confirmed by longer-wavelength studies, this may place a strong
constraint on planet formation timescales in the low-metallicity environments.
In our simulations, giant planets could be formed within $\sim 1\,$Myr even in
the low-metallicity discs as long as disc masses are high enough (see Figure
8), but the disc’s lifetimes are much longer than that. Since our current
studies do not take account of the feedback of the planets onto discs, the
future study needs to investigate this further.

###  4.3 Upper limits on masses and planet formation radii

It is unclear what the maximum planetary mass generated by the core accretion
scenario would be. Recently, Schlaufman (2018) pointed out that giant planets
with $\lesssim 4\,M_{J}$ tend to be around metal-rich dwarfs, while the same
trend is not observed for $\gtrsim 10\,M_{J}$. He suggested that these two
distinct populations may correspond to core accretion and gravitational
instability scenarios, with a transition mass being $\sim 4-10\,M_{J}$. The
transition mass is also comparable to planetary masses generated by the
gravitational instability ($\gtrsim 3-5\,M_{J}$) that are reported in the
semi-analytic studies (e.g. Kratter et al., 2010; Forgan & Rice, 2011) and the
numerical studies (e.g. Stamatellos & Whitworth, 2008; Hall et al., 2017).
Most of our simulations do not generate planets above $\sim 10\,M_{J}$, and
planetary masses of $\sim 20\,M_{J}$ were only found in massive, long-lived
discs (Disc 7). As discussed in Section 2.3, it will probably be difficult to
go above $\sim 50\,M_{J}$ with our model. Thus, although the highest-mass
planets or brown dwarfs may be formed by gravitational instability, there is
probably a significant overlap for masses of bodies both scenarios can
generate.

How far away a planet can be formed via core accretion is not well understood
either. In this work, we did not place planetary cores beyond $\sim 20\,$au
because pebble accretion becomes inefficient for a Ceres-size planetesimal to
grow in the outer disc in our model (see Section 2.3). However, as we also
show in Section 2.3, pebble accretion may still lead to giant planet formation
out to and beyond $100\,$au if much larger planetesimals form there. If there
are massive enough planetesimals available, pebble accretion (along with
planetesimal accretion) may allow the growth of giant planets far beyond
$20\,$au. Recent disc observations including ALMA regularly find gaps and
rings over a wide range of disc radii, which may be attributed to growing
planets (e.g. Zhang et al., 2016; Huang et al., 2018; Long et al., 2018).
Lodato et al. (2019) estimated planetary masses from 48 observed gaps located
around $\sim 10-130\,$au and found planetary masses ranging from
$10^{-3}\,M_{J}\sim 0.3\,M_{\oplus}$ to over $10\,M_{J}$. If these gaps are
indeed signposts of planets, the fact that a range of planetary masses is
observed in the outer disc may indicate that core accretion is active out to
these radii because gravitational instability is unlikely to form lower end
mass planets. There is, however, a discrepancy between the frequency of
observed SCJs (¡4.1 for FGK stars, Bowler, 2016) and the fact that the
substructures such as rings and gaps are ubiquitous among observed discs (e.g.
Andrews et al., 2018; Huang et al., 2018). Thus, if all of these rings and
gaps correspond to planets, such planets would have to be removed from the
outer disc region efficiently. This could be achieved if forming planets are
redistributed due to migration (Lodato et al., 2019) or removed due to
episodic accretion to the star (Brittain et al., 2020). The future studies
should investigate the connection between forming planets and evolving
protoplanetary discs further.

###  4.4 Formation of super cold Jupiters

In Section 3.1.3, we discuss two different pathways of forming SCJs via core
accretion: (1) SCJs form in situ, either directly in the outer disc or via
scattering of protoplanetary cores to the outer disc followed by gas accretion
there (e.g. Kikuchi et al., 2014), or (2) giant planets first form within
$\sim 20\,$au and are then scattered into the outer disc. These two different
scenarios are likely to lead to different distributions of physical and
orbital parameters.

First, the in situ formation scenario of SCJs may lead to the $a-M_{p}$
correlation where the planetary mass increases with the semi-major axis,
because the disc aspect ratio increases with $a,$ and thus the critical gap-
opening mass for giant planets also becomes higher (Ida et al., 2013). On the
other hand, the formation-then-scattering scenario might lead to an opposite
trend because a lower-mass planet is more easily scattered further. We checked
the Spearman’s rank correlation for simulated giant planets beyond $1\,$au
(see the left panel of Figure 11) and found the correlation coefficient of
$r_{s}\sim-0.18$ with the p-value of $p=0.011$ ($r_{s}\sim-0.19$ with the
p-value of $p=0.0015$ for all). Since the critical value for $p<0.05$ is $\mid
r_{s}\mid\sim 0.2$ for the size of our data, there is little or very weak
negative correlation between the semi-major axis and the mass for giant
planets. Due to the observational bias for high-mass planets, it is difficult
to assess whether there is any trend for observed planets.

Second, the eccentricity distributions may also be different between these two
scenarios. The formation-then-scattering scenario is likely to lead to an
isolated eccentric giant planet, while the in situ formation scenario is more
likely to lead to single- or multiple-giant planets on nearly circular orbits.
In our simulations, giant planets scattered beyond $20\,$au have
eccentricities ranging from $e\sim 0.2-0.8$ with the mean value of $0.56$ (see
Figure 5). Kikuchi et al. (2014), on the other hand, showed that the
eccentricity of the scattered core is reduced to $e<0.2$ before reaching the
Saturn mass via gas accretion in the outer disc. The eccentricities of
multiple giant planets could increase via dynamical instabilities after
formation, as long as they are not too far from each other (e.g. Chatterjee et
al., 2008).

Recently, Bowler et al. (2020) studied the eccentricity distributions of 18
brown dwarfs ($15-75\,M_{J}$) and nine giant planets ($2-15\,M_{J}$) with
orbital radii $5-100\,$au, and they found that the estimated eccentricity
distribution is very broad and ranges over $e=0-1$ for brown dwarfs and is
much narrower with a peak at $\bar{e}=0.13$ for giant planets, which may
support the in situ formation scenario of SCJs more than the formation-then-
scattering one. However, the peak eccentricity becomes $\bar{e}=0.23$ for five
giant planets excluding four HR 8799 planets on nearly circular orbits. The
future observations will further constrain the eccentricity distribution and
provide us with an important clue regarding the formation pathways of SCJs.

###  4.5 Mass fractions and metal contents

The left panel of Figure 16 shows fractions of a planetary mass in the core
(blue) and in the envelope (orange) for all the simulated planets. The
envelope should usually also contain some metals, but here we added all the
accreted dust particles to the core mass. Also plotted are the fitted metal-
mass fractions of extrasolar planets from Thorngren et al. (2016) (dashed
line) and the estimated metal fractions in Jupiter, Uranus, and Neptune (stars
with error bars, Fortney & Nettelmann, 2010; Helled et al., 2020). As we can
see, mass fractions of planetary cores of our simulations are lower than the
estimate by Thorngren et al. (2016) by a factor of a few or more, though the
distribution is consistent with the estimated metal-mass fraction of Jupiter.

This underestimation of the metal-mass fractions in our simulations may result
from our assumption of gas accretion. For our simulations, the core mass and
the envelope mass become comparable to each other: $M_{{\rm core}}\sim M_{{\rm
env}}$ when $M_{p}\sim 4\,M_{E}$. As seen in the figure, this is much lower
than the crossover mass expected from the observed planets ($\sim 15\,M_{E}$).
The crossover mass represents a group of planets that are on the verge of
becoming giant planets but failed to do so because the gas disc dissipated
(i.e. $t_{{\rm diff}}\sim\tau_{KH}$):

$\frac{M_{p}}{M_{E}}\sim 10^{\frac{b-6}{c}}\left(\frac{t_{{\rm
diff}}}{10^{6}\,{\rm yr}}\right)^{-\frac{1}{c}}\sim 4.64\left(\frac{t_{{\rm
diff}}}{10^{6}\,{\rm yr}}\right)^{-\frac{1}{3}}.$ (48)

In Matsumura et al. (2017), the crossover mass was $M_{{\rm core}}\sim M_{{\rm
env}}\sim 10\,M_{\oplus}$ as seen in Figure 11, where we adopted
$(b,\,c)=(9,\,3)$ rather than $(8,\,3)$. Thus, it is likely that we need to
adopt $(b,\,c)=(9,\,3)$ to form planets with more appropriate rock-gas mass
ratios. The reason why we adopted $(b,\,c)=(8,\,3)$ instead was to speed up
gas accretion and to avoid the loss of planets due to type I migration. As
discussed in Section 2.1, the pebble isolation mass is reached before the mass
of migration transition from type I to type II. Therefore, the gas accretion
timescale in this phase needs to be short compared to the migration timescale
for a growing planet to switch from type I to type II regime. This further
confirms the need for a proper mechanism of type I migration.

The right panel of Figure 16 compares the mass-metallicity relation for
simulated (circles) and observed (crosses) giant planets. The apparent lack of
simulated planets with $\sim 0.1-1\,M_{J}$ for the solar metallicity case is
most likely due to the choice of disc models rather than the real feature,
because planets of this mass range are undergoing rapid gas accretion and have
difficulties in forming unless the formation timescale is comparable to the
disc dissipation timescale. There is an indication that the maximum masses of
giant planets may increase for higher metallicity environments or even peak
near ${\rm[Fe/H]}\sim 0$, which is also observed by planetesimal accretion
work by Mordasini et al. (2012). Otherwise, as stated in Section 3.1.2, masses
of giant planets do not strongly depend on stellar metallicities. Thorngren et
al. (2016) studied RV and transit planets and found a similar overall trend in
the mass-metallicity relation (see their Figure 9).

Figure 16: Left: Mass fractions of the rocky core (blue) $\frac{M_{{\rm
core}}}{M_{p}}$ and the gaseous envelope (orange) $\frac{M_{{\rm
env}}}{M_{p}}$ for all the planets formed in our simulations. The core and the
envelope have comparable masses $M_{{\rm core}}\sim M_{{\rm env}}$ when a
planetary mass is $\sim 4\,M_{\oplus}$. The purple dashed line shows the metal
fractions of observed giant planets estimated by Thorngren et al. (2016).The
red stars with error bars indicate the metal-mass fractions estimated for
Jupiter, Uranus, and Neptune (Fortney & Nettelmann, 2010). Right: Planetary
mass and metallicity for simulated (circles) and observed (crosses) giant
planets. The colours indicate the final and observed semi-major axes,
respectively. Masses are generally independent of the stellar metallicity.

###  4.6 Planetesimal accretion and pebble accretion

Planet formation via planetesimal accretion has been studied and compared with
observed systems by many groups (e.g. Ida & Lin, 2004; Ida et al., 2013;
Mordasini et al., 2009, 2012; Thommes et al., 2008; Coleman & Nelson, 2016a).
Here, we compare our pebble accretion work with that of Mordasini et al.
(2012) and Ida et al. (2013). Although a direct comparison is difficult, we
discuss a few similarities and differences in planet formation outcomes.

Mordasini et al. (2012) used the single-planet formation population synthesis
model and studied the effects of the disc’s metallicity, mass, and lifetime;
while Ida et al. (2013) used a similar population synthesis model but took the
effects of multiple-planet formation and dynamical interactions into account.
The disc properties adopted in these studies are different, but they do
overlap with ours. We considered the disc radii over $0.1-100\,$au with disc
masses $0.02-0.2\,M_{\odot}$ as well as stellar metallicities
${\rm[Fe/H]}=[-0.5,\,0.5]$. Mordasini et al. (2012) considered the disc outer
radii of $30\,$au with disc masses $0.004-0.09\,M_{\odot}$ and metallicities
${\rm[Fe/H]}=[-0.5,\,0.5]$. Applying their surface mass density out to
$100\,$au, the disc masses become $0.008-0.16\,M_{\odot}$. Ida et al. (2013),
on the other hand, considered the disc outer radii of $30\,$au with disc
masses $1.6\times 10^{-4}-0.16\,M_{\odot}$ and metallicities
${\rm[Fe/H]}=[-0.2,\,0.2]$. Again, applying their surface mass density to the
disc out to $100\,$au, the disc masses correspond to $0.0053-0.53\,M_{\odot}$.

First, we highlight the difference in dependences of planet formation rates on
disc properties such as metallicities, masses, and dissipation timescales. In
terms of these disc parameters, the pebble mass accretion rate has the
following dependency (see Equations 25, 14, and 16):

$\dot{M}_{{\rm
core}}=\epsilon\dot{M}_{F}\propto\left(10^{{\rm[Fe/H]}}\right)^{2}\alpha_{{\rm
acc}}^{-1}\dot{M}_{*}\propto\left(10^{{\rm[Fe/H]}}\right)^{2}M_{D},$ (49)

where $M_{D}$ is the disc mass and the function of both the initial disc mass
$M_{D,0}$ and the dissipation timescale $t_{{\rm diff}}$. On the other hand, a
similar equation for planetesimal accretion is proportional to the dust
surface mass density $\Sigma_{d}=10^{{\rm[Fe/H]}}\Sigma_{g}$ and has the
following dependency (Mordasini et al., 2012):

$\dot{M}_{{\rm core,plsml}}\propto\Sigma_{d}\propto 10^{{\rm[Fe/H]}}M_{D}.$
(50)

Ida et al. (2013) took account of the gas drag effect on planetesimals (Kokubo
& Ida, 2002) and thus had a slightly different dependency of:

$\dot{M}_{{\rm core,plsml}}\propto 10^{{\rm[Fe/H]}}M_{D}^{7/5}.$ (51)

Therefore, in both pebble and planetesimal accretion models, increasing either
metallicity or disc mass makes planet formation more efficient. This
complementary effect of the metallicity and the disc mass was also pointed out
by Mordasini et al. (2012).

In pebble accretion, the metallicity effect is more significant than the disc
mass effect, as seen in Equation 49. More specifically, changing the
metallicity effect by an order of magnitude (e.g. from ${\rm[Fe/H]=-0.5}$ to
$0.5$ in $10^{{\rm[Fe/H]}}$) has a much larger effect on the accretion rate
than changing the disc mass by an order of magnitude (e.g. from
$0.02\,M_{\odot}$ to $0.2\,M_{\odot}$). Indeed, Figure 8 shows that the giant
planet formation timescale is different by $\sim 2$ orders of magnitude for
the same disc mass with ${\rm[Fe/H]}=-0.5$ and $0.5$, while the timescale is
different by $\sim 1$ order of magnitude for the same metallicity with
$M_{D,0}=0.02\,M_{\odot}$ and $0.2\,M_{\odot}$.

The final mass, on the other hand, is less affected by the stellar metallicity
both in planetesimal accretion (Mordasini et al., 2012) and in pebble
accretion, though the reasons may be slightly different. For pebble accretion,
the metallicity does not greatly affect final masses, partly because the PIM
is independent of the stellar metallicity (see Equation 39) and partly because
the final planetary mass is determined by the disc-limited gas accretion and
thus is controlled by the available gas disc mass. For planetesimal accretion,
however, the isolation mass depends on the dust density and thus masses of
low-mass planets or cores of giant planets are likely to be affected by
metallicities. Similarly to us, both Mordasini et al. (2012) and Ida et al.
(2013) assumed that gas accretion is limited by the disc supply towards the
end:

$\dot{M}_{{\rm gas}}\propto\dot{M}_{*}\propto M_{D},$

though we took the reduced surface mass density effect in the gap into
consideration (Equation 43), as shown in Equation 44. Thus, the final mass of
giant planets in particular is likely determined by the balance between the
disc mass and the disc dissipation timescale, rather than by the metallicity.

In our work, this can be seen in Figure 9, where planetary masses are higher
for the higher disc mass and for the longer diffusion timescale. We can also
see that planetary masses are higher for the same initial disc mass if the
dissipation time is longer (i.e. when the disc experiences the lower mass
accretion rate for a longer period of time, rather than the higher mass
accretion rate for a shorter period of time). On the other hand, planetary
masses are similar for the same disc model with different metallicities (see
Section 3.1.2 for the discussion on the increasing mass trend in $a-M_{p}$).

One major difference between pebble and planetesimal accretion is that the
formation timescale of giant planets in pebble accretion is short compared to
the disc’s lifetime over a range of disc radii (e.g. Lambrechts & Johansen,
2012). This makes it easier for giant planets to form even in less-optimal
environments such as lower-mass, lower-metallicity, and/or more rapidly
dissipating discs. Furthermore, as is also discussed in Section 4.3, pebble
accretion may allow giant planet formation far beyond the radii in which we
placed the cores for our simulations.

Finally, our Figure 5 can be compared with Figure 8 of Ida et al. (2013).
Although adopted models and initial conditions are very different between
these two studies, both reproduce the observed $M_{p}-e$, $a-e$, and $a-M_{p}$
distributions well and show similar trends in simulated distributions. Perhaps
the main difference is that, compared to Ida et al. (2013), our work has
produced more eccentric planets. Ida et al. (2013) reported that only $\sim
10\,\%$ of their planets have moderate eccentricities of $e\gtrsim 0.2$. In
our work, $\sim 26.2\,\%$ of all simulated planets (bottom panels) and $\sim
40.3\,\%$ of observable planets (middle panels) have $e\geq 0.2,$ as opposed
to $\sim 42\,\%$ of observed planets (top panels). Thus, our work agrees much
better with observations, at least for eccentricities (also see Figure 6).

There may be a few factors leading to this difference. The dynamical
instabilities generally require a few giant planets that are close enough to
each other (e.g. Chatterjee et al., 2008). Ida et al. (2013) adopted
planetesimal accretion, which is slower than pebble accretion, and therefore