# Planetary core formation via multi-species pebble accretion

G. Andama,1 N. Ndugu,2,3 S.K. Anguma,3 E. Jurua1
1Department of Physics, Mbarara University of Science and Technology, Mbarara,
Uganda
2Center for Space Research, North-West University, South Africa
3Department of Physics, Muni University, Arua, Uganda
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In the general classical picture of pebble-based core growth, planetary cores
grow by accretion of single pebble species. The growing planet may reach the
so-called pebble isolation mass, at which it induces a pressure bump that
blocks inward drifting pebbles exterior to its orbit, thereby stalling core
growth by pebble accretion. In recent hydrodynamic simulations, pebble
filtration by the pressure bump depends on several parameters including core
mass, disc structure, turbulent viscosity and pebble size. We investigated how
accretion of multiple, instead of single, pebble species affects core growth
rates, and how the dependence of pebble isolation mass on turbulent viscosity
and pebble size sets the final core masses. We performed numerical simulations
in viscous 1D disc, where maximal grain sizes were regulated by grain growth,
fragmentation and drift limits. We confirm that core growth rates and final
core masses are sensitive to three key parameters: threshold velocity at which
pebbles fragment on collision, turbulent viscosity and distribution of pebble
species, which yield a diversity of planetary cores. With accretion of
multiple pebble species, planetary cores can grow pretty fast, reaching over
30 – 40 $M_{\rm{E}}$ in mass. Potential cores of cold gas giants were able to
form from embryos initially implanted as far as 50 AU. Our results suggest
that accretion of multi-species pebbles could explain: the estimated 25 – 45
$M_{\rm{E}}$ heavy elements abundance inside Jupiter’s core; massive cores of
extrasolar planets; disc rings and gaps at wider orbits; early and rapid
formation of planetary bodies.

###### keywords:

planets and satellites:formation – planets and satellites: physical evolution
– planets and satellites: gaseous planets – protoplanetary discs –
hydrodynamics – stars: formation

††pubyear: 2021††pagerange: Planetary core formation via multi-species pebble
accretion–B

##  1 Introduction

The theory of planet formation in protoplanetary discs (hereafter, PPDs) has
developed by leaps and bounds ever since the monumental work of Safronov
(1969). In one school of thought, planets may form by gravitational collapse
of dense and dynamically cold gas disc (Kuiper, 1951; Cameron, 1978; Boss,
1997; Gammie, 2001; Rice et al., 2003; Tanga et al., 2004; Rafikov, 2005;
Durisen et al., 2007) followed by tidal downsizing (Nayakshin, 2010).
Gravitational collapse requires sufficiently massive discs and mainly favours
the formation of giant planets at the disc outskirts (Boss, 1997; Boley, 2009;
Armitage, 2010).

In another school of thought, planets may also form oligarchically by core
accretion paradigm (Wetherill, 1980; Kokubo & Ida, 1998; Thommes et al., 2003;
Coleman & Nelson, 2014). Here, micrometre-sized dust grains in the natal PPDs
first have to grow by coagulation into millimetre-centimetre (mm – cm) sized
particles. These mm – cm sized particle may concentrate in some regions of the
disc where they may gravitationally collapse into metre-kilometre sized bodies
called planetesimals (e.g., Youdin & Goodman, 2005; Johansen et al., 2007;
Raettig et al., 2015; Carrera et al., 2015). Planetesimals larger than 100 km
then form potential planetary embryos. These planetary embryos can also form
through a different mechanism, for example, when collisions between smaller
planetesimals result into a merger of over 100 km sized planetesimals (e.g.,
Kokubo & Ida, 2012). In the core accretion paradigm, the planetary embryos can
then accrete smaller planetesimals to grow into a full planet (Safronov, 1969;
Mizuno et al., 1978; Mizuno, 1980; Kokubo & Ida, 2012). However, core growth
via planetesimal accretion is typically slow unless most of the solid mass in
disc is converted into planetesimals less than 10 km in size (Tanaka & Ida,
1999; Thommes et al., 2003; Levison et al., 2010; Johansen & Bitsch, 2019).
Nevertheless, fast planetesimal-based core accretion rates were reported for
planetesimals with size less than 1 km (e.g., Mordasini et al., 2009a;
Mordasini et al., 2009b), even though there are no evidences in the solar
system for planetesimals of such smaller sizes (Bottke et al., 2005a, b;
Morbidelli et al., 2009; Singer et al., 2019).

A planetary core may also grow into a full planet by accreting aerodynamically
coupled bodies via gas drag, popularly known as pebble accretion (Johansen &
Lacerda, 2010; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Lambrechts et
al., 2014). Though the current framework of planet formation by core accretion
of planetesimals or pebbles is the most successful, it cannot satisfactorily
attribute the observed substructures at wider orbital locations to planets
(Lodato et al., 2019; Ndugu et al., 2019; Nayakshin et al., 2019).

The formation of planetary bodies by core accretion paradigm is strongly
shaped by the availability and size distribution of solid material in the PPDs
(see review by Johansen et al., 2014). In particular, pebble accretion is
constrained by the formidable radial drift barrier which leads to rapid loss
of centimetre-sized solids on short dynamical time-scales (Whipple, 1972;
Weidenschilling, 1977; Takeuchi & Lin, 2005; Alexander & Armitage, 2007;
Brauer et al., 2007; Brauer et al., 2008). Furthermore, Johansen et al. (2019)
demonstrated that mm-cm sized dust material with stokes number larger than 0.1
are expected to drain on shorter time scales than the disc life time.
Consequently such rapid loss of dust grains can impede formation of
planetesimals and the subsequent growth of planetary cores by pebble
accretion.

However, studies of PPDs from various surveys (e.g., Testi et al., 2003;
Wilner et al., 2005; Rodmann et al., 2006; Brauer et al., 2007; Pérez et al.,
2012; Trotta et al., 2013; Carrasco-González et al., 2016; Ansdell et al.,
2017) indicate that substantial amount of dust grains in the mm – cm range
survive even in the discs that are in their late stages of evolution, contrary
to theoretical predictions. The findings from the above surveys may not be
globally true. In fact, recent study by Tychoniec et al. (2020) revealed that
most discs contain small amount of pebbles and only few discs retain
substantial amount of pebbles. The retention of small dust grains in PPDs has
been linked to destructive collision of larger dust aggregates (Blum & Wurm,
2008) and subsequent coagulation-fragmentation equilibrium (Dominik &
Dullemond, 2008). Birnstiel et al. (2009) demonstrated that fragmentation of
grains could facilitate dust retention in the disc. Also zonal flows formed by
magnetorotational instability (MRI) (e.g., Johansen et al., 2009; Dzyurkevich
et al., 2010; Johansen et al., 2011; Uribe et al., 2011) can cause over-
densities and hence pressure bumps that act as dust traps and help to retain
dust grains in the outer region of the disc (Pinilla et al., 2012). Dust
trapped in the pressure bumps induced by massive planets may also undergo
fragmentation to produce finer grains (Drążkowska et al., 2019). The
retention, evolution and distribution of grain sizes play an important role in
planet formation models as they determine the outcome of the planetary bodies
(Barrière-Fouchet et al., 2005), as well as the disc structure (Dullemond &
Dominik, 2004).

The grain size distribution can be approximated using either the simple power-
law fits as in the MRN model (Mathis et al., 1977) or the complex analytic
approach of Birnstiel et al. (2011). Birnstiel et al. (2012) further developed
a two-dust population model of dust size distribution, classifying the
distribution into opacity bearing grains and larger grains. The subpopulation
of the smaller grains determines the temperature profile and hence the disc
structure (Dullemond & Dominik, 2004; Savvidou et al., 2020).

The discussed observational surveys of dust disc and dust distribution
simulations present piece of evidence for the existence of numerous dust
species in PPDs, each with unique spatio-temporal distribution of Stokes
number. The assortment of different grain sizes may point to the fact that
core accretion proceeds by accretion of multiple pebble species.

However, previous works that studied core growth via pebble accretion (e.g.,
Lambrechts et al., 2014; Lambrechts & Johansen, 2014; Guillot et al., 2014;
Morbidelli et al., 2015; Bitsch et al., 2015; Bitsch & Johansen, 2017; Ndugu
et al., 2018; Brügger et al., 2018; Johansen et al., 2019; Ndugu et al., 2019;
Ndugu et al., 2021) were based on two standard prescriptions. Firstly, the
studies used a single spatio-temporal dominant particle size or Stokes number,
which is assumed to carry most of the solid mass. Dust coagulation models
(e.g., Dullemond & Dominik, 2005) predict rapid conversion of most of the dust
into larger grains within very short time scales compared to the disc
lifetimes. The grain population can thus be modeled as consisting of small-
and large-size grain populations, where most of the mass is carried by the
large-size population (Birnstiel et al., 2012). Through coagulation-
fragmentation equilibrium, grain sizes may attain quasi-stationary size
distribution (Dullemond & Dominik, 2005). A quasi-stationary size means that
one is justified to use single size approach, especially when the resulting
grain sizes in the distribution are very similar. Nevertheless, in this work,
instead of a single pebble size or Stokes number, we have studied how grain
size distributions reconstructed from the dust evolution model of Birnstiel et
al. (2012) influences planetary core growth through accretion of different
pebble species. This is because considering only dominant pebble species may
not provide a complete picture of the final core masses and growth rates. This
could be important, especially when several other pebble species in the
distribution carry some considerable mass. Consequently, we can easily
underestimate final core masses or growth times if some of the species that
may also contain significant mass are neglected during core assembly, as may
be the case in single species approach.

Secondly, core growth stops whenever a pressure bump is induced by the growing
core, and hence the final core masses are fixed by the classical pebble
isolation mass (Lambrechts et al., 2014, hereafter L2014). However,
hydrodynamical simulations show that through turbulence, some pebble species
can still diffuse through the pressure bump exerted by sufficiently massive
cores and can even pass into gaps carved out by Jupiter-mass planets (Weber et
al., 2018). As a result of turbulent diffusion, Bitsch et al. (2018)
(hereafter, B2018) and Ataiee et al. (2018) further demonstrated that pebble
isolation mass might not be completely universal for all pebble species since
smaller pebble species may overcome weaker pressure bumps. This therefore
suggests that those pebble species that diffuse through the pressure bump may
sustain core accretion for an extended period.

The grain size distribution as in the two-population model of Birnstiel et al.
(2012) now provides an opportunity to study core accretion in the context of
multiple dust species (Guilera et al., 2020; Venturini et al., 2020;
Drążkowska et al., 2021; Savvidou & Bitsch, 2021; Schneider & Bitsch, 2021).
For instance, in both Guilera et al. (2020) and Venturini et al. (2020), the
authors used a mass weighted representative pebble size derived from dust
population of several species in their pebble accretion model, where the final
core masses are fixed by classical L2014 prescription. Using full grain size
distribution, Guilera et al. (2020) studied how giant planet cores can form by
hybrid pebble and planetesimal accretion at pressure maxima. In Venturini et
al. (2020), the authors focused on the formation of super – Earths inside the
snow line. On the other hand, Schneider & Bitsch (2021) studied the heavy
element content of giant planets, also based on full dust evolution, where
they focused on evolution of different chemical species rather than the
physical size distributions as considered in this work.

In Drążkowska et al. (2021), using full size distribution, the authors studied
the impact of grain growth and fragmentation on core growth rate while the
planet accreted throughout at a particular radial location. This allowed the
authors to exclusively study the actual impact of fragmentation on core
growth. Hence they performed a limited set of simulations and did not take
into account many aspects of planet formation such as gas accretion and
orbital migration, which are crucial for the final architecture of planetary
systems. For instance, inward migration could enable a core to reach pebble
isolation mass much faster in the disc regions closer to the star where the
isolation mass is lower. This then gives the planet a chance to accrete gas
and grow into a gas giant. Lastly, in their work the core masses were measured
based on the L2014 model without taking into account the diffusion of pebbles
across the pressure bump, which may impact the final core mass as discussed in
B2018.

In this study, we performed similar numerical simulations, but quite different
from the one presented in Drążkowska et al. (2021). The major difference in
our work is that we self-consistently reconstruct a distribution of pebble
sizes, their corresponding mass and Stokes numbers after full dust evolution
at every time step during core growth. This allowed us to investigate the
contribution of each individual grain species as opposed to the mass averaged
pebble flux model used in previous studies. Here we focused on two major
problems: how concurrent accretion of several pebble species may impact core
growth; how the dependence of pebble isolation mass on turbulent viscosity and
pebble sizes may determine final core masses. All this is based on the two-
population model of Birnstiel et al. (2012) and the size reconstruction recipe
in Birnstiel et al. (2015). We did not perform hydrodynamic simulations of
pebble isolation mass because it is beyond the scope of this work but rather
used the formula from B2018 to calculate the final masses.

The rest of the paper is structured as follows: In Section 2, we describe the
underlying disc model, size distribution of particles and the core growth
model. In Section 3, we explain the main numerical experiments that we
performed. We present and discuss our results in Section 4. We then summarise
our findings in Section 5.

##  2 Theoretical background

###  2.1 The disc evolution model

In order to provide a complete picture of our core accretion of multi-pebble
species model, we employed the two-population model of Birnstiel et al.
(2012). Here, we only describe the key ingredients of the gas and dust
evolution model that we adopted in our numerical simulations and we refer the
reader to Birnstiel et al. (2009, 2010, 2011, 2012) for more detailed
description of dust size evolution.

In the simulations, the initial gas surface density $\Sigma_{\rm{g}}$ is
calculated using the self-similar solution of Lynden-Bell & Pringle (1974)
given by

$\Sigma_{\rm{g}}(r)\propto\left(\frac{r}{r_{\rm{c}}}\right)^{-\gamma}{\rm{exp}}\left[-\left(\frac{r}{r_{\rm{c}}}\right)^{2-\gamma}\right].$
(2.1)

Here, $\gamma$ is the viscosity power-law index and $r_{\rm{c}}$ the
characteristic radius at an initial time $t_{0}$.

For dust evolution in the disc, we adopt the two dust population model of
Birnstiel et al. (2012), where the dust surface density evolves according to
the advection-diffusion equation:

$\frac{\partial\Sigma_{\rm{p}}}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial
r}\left[r\left(\Sigma_{\rm{p}}\bar{u}-D_{\rm{g}}\Sigma_{\rm{g}}\frac{\partial}{\partial
r}\left(\frac{\Sigma_{\rm{p}}}{\Sigma_{\rm{g}}}\right)\right)\right]=0.$ (2.2)

Here, $\Sigma_{\rm{p}}$ is the dust and gas surface densities and $D_{\rm{g}}$
is the gas diffusivity. $\bar{u}$ is the velocity of the dust weighted by the
mass of the two dust populations and is given by

$\bar{u}=(1-f_{\rm{m}})u_{0}+f_{\rm{m}}u_{1},$ (2.3)

where $u_{0}$ and $u_{1}$ are respectively radial velocities of the two small-
and large-size grain populations, with representative sizes $a_{0}$ and
$a_{1}$ which are set by growth, drift and fragmentation limits as described
in Birnstiel et al. (2012). Here, the mass fraction $f_{\rm{m}}$ of the large-
size grain population at radial distance $r$ is calibrated as 0.97 and 0.75
for drift and fragmentation limited regimes, respectively.

The radial velocities $u_{0}$ and $u_{1}$ are calculated as the sum of radial
drift velocity and the radial velocity due to gas drag (Weidenschilling,
1977):

$u_{i}=-\frac{2\tau_{\rm{i}}}{1+\tau_{\rm{i}}^{2}}u_{\eta}+\frac{1}{1+\tau_{\rm{i}}^{2}}u_{\rm{g}},$
(2.4)

where $u_{\rm{g}}$ is the gas velocity, $u_{\rm{i}}~{}(i=0,1)$ represents
either $u_{0}$ or $u_{1}$, and $\tau_{\rm{i}}$ the corresponding Stokes
numbers, which will be discussed in Section 2.2. Here $u_{\eta}$ is the
headwind velocity (Weidenschilling, 1977; Nakagawa et al., 1986) given by

$u_{\eta}=-\frac{1}{2}\frac{1}{\rho_{\rm{g}}\Omega_{\rm{K}}}\frac{\partial~{}P}{\partial~{}r}.$
(2.5)

where $\rho_{\rm{g}}$ is the mid-plane gas density, $P$ the gas pressure and
the $\Omega_{\rm{K}}$ the Keplerian frequency. The gas velocity $u_{\rm{g}}$
is given by

$u_{\rm{g}}=c_{\rm{s}}\sqrt{1.5\alpha_{\rm{t}}},$ (2.6)

where $c_{\rm{s}}$ is the sound speed and $\alpha_{\rm{t}}$ is the turbulence
parameter.

Figure 2.1: The geometry of vertical particle distribution used in our core
growth model, where the particle species change from 1 – 3 with decreasing
size.

###  2.2 Particle size distribution

The distribution of solids in PPDs depends on their aerodynamic properties
characterised by friction time, $t_{\rm{f}}$, given by (Whipple, 1972;
Weidenschilling, 1977) as

$t_{\rm{f}}=\frac{mv}{||F_{\rm{D}}||},$ (2.7)

where $m$ is the mass of the particle, $v$ the terminal velocity and
$||F_{\rm{D}}||$ the gas drag force.

For small particles, it is usually convenient to express their degree of
coupling with the gas in terms of the dimensionless Stokes number
$\tau_{\rm{s}}$ given by

$\tau_{\rm{s}}=t_{\rm{f}}\Omega_{\rm{K}}=\frac{\rho_{\bullet}R}{\rho_{\rm{g}}h_{\rm{g}}}=\frac{\rho_{\bullet}\pi
R}{2\Sigma_{g}},$ (2.8)

where $\rho_{\bullet}$ is the material density, $R$ the particle size,
$\Sigma_{\rm{g}}$ the gas surface density, $h_{\rm{g}}$ the gas scale height.

Bodies with Stokes number $\tau_{\rm{s}}>1$ become increasingly less coupled
to the gas or even decouple completely, and therefore are not well suited for
core growth by aerodynamic drag in the settling regime. On the other hand,
bodies with Stokes number smaller than 0.001 have short friction time and
hence long settling time that enables them to stay coupled to the gas during
their gravitational encounter with the growing cores. Therefore, these
particles mostly follow the gas streamlines and may not be accreted
efficiently (Guillot et al., 2014; Johansen et al., 2019; Rosenthal & Murray-
Clay, 2020). Furthermore, as shown by Johansen et al. (2019), pebbles with
Stokes number $\sim 0.1$ drift faster and are lost to the central star on
short dynamical scales.

Particles settle towards the mid-plane depending on their size and material
density and may also radially drift on time-scales which equally depend on
their aerodynamic properties (Whipple, 1972; Adachi et al., 1976;
Weidenschilling, 1977; Dullemond & Dominik, 2004). For higher radial drift
velocities, larger particles tend to concentrate more in the inner disc
regions than the smaller particles (Barrière-Fouchet et al., 2005; Birnstiel
et al., 2010, 2011; Testi et al., 2014). This results into vertical and radial
stratification of particles (see Figure 2.1).

In a quiescent disc, small dust grains are found to settle to the mid-plane,
while the larger grains tend to oscillate about the mid-plane as their
oscillation amplitude decays (Garaud et al., 2004). Here, we assume a simple
model with N particle species, where particles sediment in a layered fashion
with different scale heights (Barrière-Fouchet et al., 2005). Particles with
small Stokes number stay in the upper layers of the gas disc, and those with
large Stokes number settle in the mid-plane as depicted in the schematic shown
in Figure 2.1. In this picture, the scale height, $h_{\rm{i}}$ (with i=1, 2,
…, N denoting different species) of each particle species is then calculated
from (Youdin & Lithwick, 2007) as

$h_{\rm{i}}=\left(\frac{\alpha_{\rm{t}}}{\tau_{\rm{i}}}\right)^{1/2}h_{\rm{g}},$
(2.9)

where $\tau_{\rm{i}}$ is the Stokes number of a particular species.

The differential settling of particles discussed above suggests that the upper
layers of the gas disc become relatively devoid of dust even for $\mu$m grains
(Dullemond & Dominik, 2004). This motivates us to consider N different pebble
species with surface densities, $\Sigma_{\rm{p,i}}$. The simple MRN model of
Mathis et al. (1977) describes grain size distribution quite well. However, we
used the two-population model of Birnstiel et al. (2012) in our work from
which we reconstruct the grain sizes and their surface densities,
$\Sigma_{\rm{p,i}}$ as presented in Birnstiel et al. (2015). This is because
we want to use a self-consistent size distribution which takes into account
the balance between grain growth, fragmentation and drift barrier as described
in Birnstiel et al. (2012).

Our grain sizes are distributed such that $R_{\rm{i+1}}=1.12R_{\rm{i}}$ as in
Birnstiel et al. (2011) where $R_{\rm{i}}$ is the size of the i-th species.
The dust surface density, $\Sigma_{\rm{p,i}}$, corresponding to each species
at a given radial distance and time is then reconstructed from gas surface
density ($\Sigma_{\rm{g}}$), dust surface density ($\Sigma_{\rm{p}}$),
fragmentation velocity ($\it{u}_{\rm{f}}$), turbulent strength
($\alpha_{\rm{t}}$), material density ($\rho_{\bullet}$) and the disc mid-
plane temperature (T) according to the recipe given in Birnstiel et al.
(2015).

Laboratory experiments have constrained the fragmentation velocity threshold
for silicate grains to be $\it{u}_{\rm{f}}$ = 1 m/s (Blum & Wurm, 2008), while
numerical simulations show that water-ice aggregates which can grow to
centimetre sizes tend to fragment at much higher velocities of
$\it{u}_{\rm{f}}\geq$ 10 m/s (Brauer et al., 2008; Wada et al., 2008; Gundlach
et al., 2011; Gundlach & Blum, 2015). But recent laboratory experiment by
Musiolik & Wurm (2019) seem to suggest that ice and silicate grains have
similar fragmentation velocities. However, the collisional outcomes also
depend on the turbulent strength, the internal density of the solid bodies and
the local temperature which can lead to wide range of particle sizes. Lastly,
the corresponding Stokes number, $\tau_{i}$, of each pebble species from the
reconstructed size distribution is then calculated using Equation (2.8)

###  2.3 Dominant size distribution

In the two-population model of Birnstiel et al. (2012), the surface densities
$\Sigma_{0}$ and $\Sigma_{1}$ of the small and large populations at a radial
distance $r$ are calculated as

$\displaystyle\Sigma_{0}$ $\displaystyle=\Sigma_{\rm{p}}(r)(1-f_{\rm{m}}(r))$
(2.10) $\displaystyle\Sigma_{1}$
$\displaystyle=\Sigma_{\rm{p}}(r)f_{\rm{m}}(r).$ (2.11)

In typical classical pebble accretion scenarios, a single representative size
is assumed to contain most of the dust mass, for example, in the large
population. In this classical picture we can calculate the mass averaged
dominant size $R_{\rm{d}}$ for the reconstructed size distribution as (Guilera
et al., 2020; Venturini et al., 2020)

$R_{\rm{d}}=\frac{\sum_{i}\epsilon_{\rm{i}}R_{\rm{i}}}{\sum_{\rm{i}}\epsilon_{\rm{i}}},$
(2.12)

where

$\epsilon_{\rm{i}}=\frac{\Sigma_{\rm{p,i}}}{\Sigma_{\rm{g}}}\left(\frac{\alpha_{\rm{t}}+\tau_{\rm{i}}}{\alpha_{\rm{t}}}\right)^{1/2}.$
(2.13)

The corresponding Stokes number for the dominant species with size
$R_{\rm{d}}$ can then be obtained by using Equation (2.8).

###  2.4 Core growth model

In our planetary growth model, planetary embryos start accreting pebbles at
the transition mass where Hill accretion becomes more efficient than the Bondi
regime (Lambrechts & Johansen, 2012). For the governing equations, we start
from the classical pebble accretion rate of a dominant pebble species $i$ in
the 2-D Hill regime given by (Morbidelli et al., 2015):

$\dot{M}_{\rm{2D}}=\left\\{\begin{array}[]{ll}2\left(\tau_{\rm{i}}/{0.1}\right)^{2/3}\Omega_{\rm{K}}r_{H}^{2}\Sigma_{\rm{p,i}}&(\tau_{\rm{i}}<0.1)\\\
\\\ 2\Omega_{\rm{K}}r_{H}^{2}\Sigma_{\rm{p,i}}&(\tau_{\rm{i}}\geq
0.1)\end{array}\right.$ (2.14)

where $r_{\rm{H}}$ is the Hill radius.

The 2-D solid accretion takes place when the effective accretion radius is
greater than the particle scale height, otherwise the core grows by 3-D
accretion. Since the protoplanets are initially small, their gravitational
reach is mostly below the particle scale height and hence the accretion rate
follows the 3-D mode. The 3-D and 2-D accretion modes are then related as in
Morbidelli et al. (2015) by

$\dot{M}_{\rm{3D}}=\left[\sqrt{\frac{\pi}{8}}\left(\frac{\tau_{\rm{i}}}{0.1}\right)^{1/3}\frac{r_{\rm{H}}}{h_{\rm{i}}}\right]\dot{M}_{\rm{2D}},$
(2.15)

where we can calculate the critical core masses at which accretion switches
from 3-D to 2-D as

$M_{\rm{3D2D}}=4.06\times
10^{5}\times\alpha_{\rm{t}}^{3/2}\tau_{\rm{i}}^{-5/2}\left(\frac{H}{r}\right)^{3}~{}M_{\rm{E}},$
(2.16)

Where $H/r$ is the disc aspect ratio. Hence the core accretion rate of the
$i{\rm{-th}}$ pebble species is given by

$\dot{M}_{\rm{core,i}}=\left\\{\begin{array}[]{ll}\dot{M}_{\rm{2D}}&\rm{for~{}}~{}\sqrt{\frac{\pi}{8}}\left(\frac{\tau_{\rm{i}}}{0.1}\right)^{1/3}r_{\rm{H}}>h_{\rm{i}}\\\
\dot{M}_{\rm{3D}}&~{}\text{otherwise}\end{array}\right.$ (2.17)

Figure 2.2: The planet mass at which pebble accretion switches from 3-D to 2-D
as function of the Stokes number, for $\alpha_{\mathrm{t}}=10^{-3}$ and
$\alpha_{\mathrm{t}}=10^{-4}$. The colour scale shows the mass
$M_{\rm{core,i}}$ accumulated by the planet from the $i\rm{-th}$ pebble
species, which is what each pebble species contributed to the final planet
mass.

In Figure 2.2, we illustrate the mass at which a growing planet switches from
3-D to 2-D accretion as a function of pebble size at 5 AU, 10 AU and 15 AU for
turbulent strengths of $\alpha_{\mathrm{t}}=10^{-3}$ and
$\alpha_{\mathrm{t}}=10^{-4}$. In the first place, the switch from 3-D to 2-D
in discs with low turbulence occurs at a lower core mass, and the converse is
true for high disc turbulences. This is because in less turbulent disc the
pebble scale heights are low, which increases the
$r_{\mathrm{H}}/h_{\mathrm{i}}$ ratio at a fixed $r_{\mathrm{H}}$ and hence
constant core mass. Secondly, the accretion of smaller pebbles transitions
from 3-D to 2-D mode at higher core masses. For the case of larger pebbles,
this transition occurs at much lower core masses. This is because, for the
same turbulence strength, larger pebbles typically settle closer to the mid-
plane and hence have low scale heights compared with the smaller pebbles.

Also in Figure 2.2 the contribution of particles of a given Stokes number to
the total mass is shown by the colour plot, computed using Equation (2.17).
Here, pebble distributions were taken from grain size reconstruction method
described in more detail in Section 3. From the illustration in Figure 2.2,
pebbles with smaller Stokes number do not significantly contribute to the
total core mass. For example, at 5 AU, pebbles with Stokes number less than
0.001 contribute roughly less than $10^{-4}~{}M_{\rm{E}}$ to a planetary core
of 40 $M_{\rm{E}}$. This is because the smaller pebbles are not easily
accreted compared with the large pebbles which carry most of the mass.

###  2.5 Pebble isolation mass

As the core grows massive enough, at some point it will begin to open a gap
and induce a pressure bump at the outer edge of the gap where pebbles get
trapped (Morbidelli & Nesvorny, 2012; Paardekooper & Mellema, 2006).
Consequently, at a critical mass usually referred to as the pebble-isolation
mass, the core stops accreting pebbles (; ; Ataiee et al., 2018).

However, recent hydrodynamical simulations suggest that pebbles with small
Stokes number can still cross to the inner disc through a gap carved out by a
Jupiter-mass planet (Weber et al., 2018, and references therein), where
pebbles with larger Stokes number are more efficiently filtered out. Since the
pressure bump scales with planet mass, where typical pebbles isolation masses
are an order of magnitude lower than a Jupiter-mass planet (Morbidelli &
Nesvorny, 2012; Lambrechts et al., 2014; Bitsch et al., 2018), the growing
cores may not efficiently filter smaller pebbles. The smaller pebbles may thus
overcome the pressure and hence be accreted.

The pebble isolation mass as originally developed by L2014 without
consideration of turbulent effects is given by

$M_{\rm{iso}}=20\left(\frac{H/r}{0.05}\right)^{3}\rm{M_{\rm{E}}}\text{ .}$
(2.18)

B2018 and Ataiee et al. (2018) built on L2014 and investigated pebble
isolation mass taking into account turbulent diffusion and obtained similar
results where turbulence can significantly change pebble isolation mass. One
could use either formulation to study core growth by pebble accretion, where
the resulting core masses should only differ by a factor of 1.5 – 2 as
discussed in Ataiee et al. (2018). However, the B2018 and L2014 prescriptions
give very similar results in the limit of weak turbulence as compared with
Ataiee et al. (2018). Hence for consistence we follow B2018 since we have used
classical model of L2014 to compare the final core masses with and without
diffusion of pebbles across the pressure bump. B2018 derived an expression for
pebble isolation mass with diffusion as

$M_{\rm{iso}}=M_{\rm{iso}}^{\dagger}+\frac{\Pi_{\rm{crit}}}{\lambda}~{}M_{\rm{E}}\text{
,}$ (2.19)

where $M_{\rm{iso}}^{\dagger}$ is the pebble isolation mass without diffusion,
$\Pi_{\rm{crit}}$ is the critical pressure gradient parameter, $\lambda$
defines the change in slope of pressure gradient inside the pressure bump
generated by the planet. $M_{\rm{iso}}^{\dagger}$, $\Pi_{\rm{crit}}$ and
$\lambda$ are given by

$M_{\rm{iso}}^{\dagger}=25f~{}M_{\rm{E}},~{}\lambda\approx\frac{0.00476}{f},~{}\Pi_{\rm{crit}}=\frac{\alpha_{\rm{t}}}{2\tau_{\rm{i}}}.$
(2.20)

Here $f$ is a fit to the isolation mass given by

$f=\left[\frac{H/r}{0.05}\right]^{3}\left[0.34\left(\frac{\rm{log}(\alpha_{3})}{\rm{log}(\alpha_{\rm{t}})}\right)^{4}+0.66\right]\left[1-\frac{\frac{\partial{\rm{ln}}P}{\partial{\rm{ln}}r}+2.5}{6}\right]\text{,}$
(2.21)

where $\alpha_{3}=0.001$ is the scaling factor and
$\partial{\rm{ln}}P/\partial{\rm{ln}}r$ is the pressure gradient.

From Equations (2.19), (2.20) and (2.21), we can calculate the pebble
isolation mass with turbulent diffusion, $M_{\rm{iso,i}}$, for each pebble
species as

$M_{\rm{iso,i}}=17.51\left(\frac{H/r}{0.05}\right)^{3}\left[0.34\left(\frac{\rm{log\alpha_{3}}}{\rm{log\alpha_{\rm{t}}}}\right)^{4}+0.66\right]\\\
\times\left(3.5-\frac{\partial{\rm{ln}}P}{\partial{\rm{ln}}r}\right)\left(0.238+\frac{\alpha_{\rm{t}}}{\tau_{\rm{i}}}\right)~{}M_{\rm{E}}.$
(2.22)

Figure 2.3: The pebble isolation mass as a function of pebble size and orbital
distance for nominal disc turbulence parameters $\alpha_{\mathrm{t}}=10^{-3}$
and $\alpha_{\mathrm{t}}=10^{-4}$.

To put Equation (2.22) in perspective, the dependence of pebble isolation mass
on grain size, turbulence viscosity and orbital distance is illustrated in
Figure 2.3. For each grain species, the isolation mass increases with orbital
distance, and with turbulence levels. Additionally, at a particular radial
distance, pebble isolation mass increases with decreasing Stokes number since
from Equation (2.22), pebble isolation mass is inversely related to the Stokes
number.

As a consequence of Equation (2.22), a growing planet may block pebbles at
different stages as illustrated in Figure 2.4, which shows what pebble species
a non-migrating planet accretes as it reaches different masses. For example,
at 5 AU a 20 $M_{\mathrm{E}}$ planet is accreting all pebble species. When the
planet reaches 60 $M_{\mathrm{E}}$, it has blocked all pebble species to the
right of the dashed vertical line and can only accrete pebbles to the left of
the vertical line. At higher orbital distances, the planet needs to grow
bigger before it can start filtering larger pebbles.

Figure 2.4: An illustration of how, using Equation (2.22), a growing and non-
migrating planet may filter pebbles with different Stokes numbers at 5 AU, 10
AU and 15 AU for a nominal disc turbulence parameter of
$\alpha_{\mathrm{t}}=10^{-3}$, fragmentation velocity of 10 m/s and dust-to-
gas ratio of 0.01, where grain sizes were obtained from the method described
in Section 3.

###  2.6 Planetary migration scheme

We implement orbital evolution during core growth and follow the Paardekooper
et al. (2011) prescription, where cores experience Lindblad and corotation
torques $\Gamma_{\rm{L}}$ and $\Gamma_{\rm{C}}$, respectively with the total
torque $\Gamma_{\rm{tot}}$ given by

$\Gamma_{\rm{tot}}=\Gamma_{\rm{L}}+\Gamma_{\rm{C}}.$ (2.23)

Here, the Lindblad torque is expressed as

$\frac{\gamma\Gamma_{\rm{L}}}{\Gamma_{0}}=-2.5-1.7\beta+0.1s,$ (2.24)

where $\gamma=1.4$ is the adiabatic index, $s$ and $\beta$ are the negatives
of radial gradients of gas surface density, $\Sigma_{\rm{pla}}$, and the
temperature, $T$, respectively calculated at the planet’s location,
$r_{\rm{pla}}$. $\Gamma_{0}$ is expressed as

$\Gamma_{0}=\left(\frac{q}{H/r}\right)^{2}\Sigma_{\rm{pla}}r_{\rm{pla}}^{4}\Omega_{\rm{pla}}^{2},$
(2.25)

where $q$ is the planet-star mass ratio and $\Omega_{\rm{pla}}$ the Keplerian
frequency and the disc aspect ratio $H/r$ is calculated at planet’s location.
The corotation torque which is induced by material corotating with the
planetary body is composed of barotropic and entropy related parts. For
detailed discussion, we refer the reader to the review by Baruteau et al.
(2014). The corotation torque is calculated using the formula (Paardekooper et
al., 2010)

$\frac{\gamma\Gamma_{\rm{C}}}{\Gamma_{0}}=1.1\left(\frac{3}{2}-s\right)+7.9\frac{\xi}{\gamma},$
(2.26)

where the first and second terms are the barotropic and entropy components of
the corotation torque and $\xi=\beta-(\gamma-1)s$ is the radial entropy
gradient.

In our model, we only incorporated type I migration but not type II because
our model is limited to solid core growth just before pebble isolation mass,
making type II migration not necessary.

##  3 Numerical simulations

The full dust size population in Birnstiel et al. (2010, 2011, 2012) features
a large number of different dust species, broadly classified into small and
large dust size population. Some of these species are either accretable or
just contribute to the disc opacity, which then determines the temperature
profile of the disc.

We incorporated the two-population dust evolution code of Birnstiel et al.
(2012)111https://github.com/birnstiel/two-pop-py in our numerical code for
dust evolution. In the simulations, particles evolve both in time and space as
the disc evolves, governed by balance between grain growth, fragmentation and
drift size limits. We then reconstructed the surface density of each particle
species for a sample of 150 species at every time step, using the size
distribution reconstruction code of Birnstiel et al.
(2015)222https://github.com/birnstiel/Birnstiel2015_scripts.

We performed our simulations in axisymmetric 1-D disc, with 200
logarithmically spaced radial grid points. For the global disc evolution, the
computational grid extends from 0.05 – 3000 AU with characteristic radius,
$r_{\rm{c}}=200$ AU, while we implant 0.01 $M_{\rm{E}}$ planetary embryos
between 1 – 50 AU. The central star has mass, $M_{*}$ = 1.0 $M_{\odot}$,
temperature, $T_{*}$ = 5778 K, radius, $R_{*}$ = 1.0 $R_{\odot}$. We assume a
disc mass, $M_{\rm{disc}}$ = 0.1 $M_{\odot}$, which gives $\sim
330~{}M_{\rm{E}}$ of dust mass for a nominal solid-to-gas ratio of 0.01. This
is within the range of dust masses measured in different star-forming regions,
especially for some Class 0 disc systems (see, e.g., Manara et al., 2019;
Tychoniec et al., 2020).

In the the simulations, we tested three different initial dust-to-gas ratios,
$f_{\rm{DG}}$ = 0.01,0.015, 0.02 with the following combinations of
fragmentation velocity and turbulence parameter:

  1. –

$u_{\rm{f}}=10~{}\rm{m/s},~{}\alpha_{t}=10^{-3}$;

  2. –

$u_{\rm{f}}=10~{}\rm{m/s},~{}\alpha_{t}=10^{-4}$;

  3. –

$u_{\rm{f}}=1~{}\rm{m/s},~{}\alpha_{t}=10^{-3}$;

  4. –

$u_{\rm{f}}=1~{}\rm{m/s},~{}\alpha_{t}=10^{-4}$.

Here, the turbulence parameter $\alpha_{\mathrm{t}}$ regulates grain size in
the fragmentation regime, pebble scale height and migration, which all change
at the same time when $\alpha_{\mathrm{t}}$ changes. The disc temperature is
assumed to be constant in time and varies only with radial distance as defined
in the two-population code.

We performed two sets of simulations. In the first part, we adopted the
classical core accretion of dominant pebble size where the isolation mass is
governed by the classical Equation (2.18) of L2014. In our dominant species
model, we used the full size distribution to calculate the mass averaged
Stokes number and the surface density corresponding to the dominant size, as
described in Section 2.3.

In the second part, we feed the full size distribution into the core accretion
routine, where pebbles of a given Stokes number are accreted independently. In
this scenario, pebble isolation mass of the individual species is set by
Equation (2.22) of B2018, and we calculate the accretion rates according to
the recipe described in Appendix B. Furthermore, we terminate the accretion of
full size distribution and measure the planet mass when the core has reached
the isolation mass of pebble species with Stokes number $\approx$ 0.001. This
is because, as we saw in Figure 2.2, the contribution of pebbles with Stokes
number less than 0.001 is small and stopping their accretion after the
isolation mass of pebbles with Stokes number $\approx$ 0.001 is reached will
not significantly affect our results.

Migration of planetary cores also take place during core growth. We terminate
core growth when the planet has migrated to 0.05 AU or when it has reached the
pebble isolation mass. If the core fails to reach pebble isolation mass,
growth continues until the end of time evolution of the disc.

##  4 Results and discussion

Figure 4.1: Top panel: The final core masses as a function of the starting
position resulting from accretion of dominant pebble species. The dominant
pebble size is calculated from the reconstructed pebble size distribution as
derived from the full dust evolution model prescribed in Birnstiel et al.
(2012). Core growth is measured within a period of 2 Myr of gaseous disc life
span for different combinations of fragmentation velocity and turbulence
levels. The final planetary cores become increasingly more massive if the
dust-to-gas ratio $f_{\rm{DG}}$ increases from 0.01 to 0.02. Bottom panel: The
corresponding total growth times of the planetary cores for the same values of
parameters on the left panel. The planetary cores with growth times below 2
Myrs reached their pebble isolation mass, and the growth time to reach
isolation mass reduces with increasing dust-to-gas ratios.

We present and discuss results for the two pebble accretion models. In the
first part, we start by discussing the classical scenario of core growth with
dominant size distribution, where the isolation mass is determined by Equation
(2.18). Here, the dominant sizes were obtained by computing the mass averaged
dust size and Stokes number from the full grain size distribution
reconstructed from the two-population model (Birnstiel et al., 2012). In the
second part we present results of our model featuring concurrent accretion of
multi-pebble species, where pebble isolation mass is controlled by Equation
(2.22). We then point out the differences between the two models and discuss
what it means for formation of planetary systems by core accretion paradigm in
relation to current observations. Although the transition mass is sensitive to
initial position of the embryo, throughout our work, we assumed embryos with
transition mass $0.01M_{\rm{E}}$ for all initial orbital positions considered
here.

###  4.1 Core growth with dominant pebble species

Figure 4.1 represents core growth timescales and final masses that were
obtained from core growth by accretion of dominant size. Here, we show results
for different sets of fragmentation velocity and turbulence parameter. The
plots in top and bottom rows respectively represent the final core masses and
total growth times as a function of starting positions of the embryos. Here,
we first determined at what initial orbital positions our planetary embryos
can grow and reach pebble isolation mass. From the plots, the planetary cores
that took 2 Myrs had not reached their pebble isolation mass. We note here
that migration of planetary cores also took place during the entire process of
core growth (their final orbital distances and growth times will be discussed
in more detail later on in Section 4.3).

From Figure 4.1, our simulations with $u_{\rm{f}}=1~{}\rm{m/s}$ and
$\alpha_{\rm{t}}=10^{-3}$ produced no significant growth of planetary embryos
over the 2 Myr of disc evolution. This is because low fragmentation velocities
in turbulent discs keep overall grain sizes small and more coupled to the gas,
which makes their accretion difficult. Accretion becomes even more difficult
when planetary cores are very small, where small dust grains simply drift past
the embryo (Guillot et al., 2014). This is because the small embryos do not
have strong enough gravitational force to pull off small grains that are
strongly attached to the gas.

For the simulations performed with fragmentation velocity of $10~{}\rm{m/s}$,
and turbulence strengths of $10^{-3}$ and $10^{-3}$, planetary cores reached
pebble isolation mass only when growth of the embryos started at a radial
distance within 20 AU for nominal dust-to-gas ratio of 0.01. The corresponding
growth timescales for the cores implanted within 20 AU to reach their
isolation mass ranges between 0.1 – 1.5 Myrs as shown in the bottom panels of
Figure 4.1.

With increasing dust-to-gas ratio, core growth is boosted and planetary
embryos that are introduced as far as 35 AU can now reach the pebble isolation
mass. At each initial radial distance, the growth timescales are also greatly
improved when the dust-to-gas ratio increases. The improved growth rates at
higher dust-to-gas ratios can be related to the higher dust surface densities
that allow cores to grow faster and bigger. Furthermore, radial drift motions
of dust also depend on dust-to-gas ratio, which when increased, can bring
about fast inward transport of solid material from the outskirts of the disc.
This then increases pebble flux and accretion efficiency at the planet’s
location, and hence the embryos may effectively grow bigger with reduced
growth time.

From Figure 4.1 where $\alpha_{\rm{t}}=10^{-4}$ and $u_{\rm{f}}=1~{}\rm{m/s}$,
we registered much better final core masses in comparison with the other
models. Here, planetary cores can grow to 20 – 30 $M_{\rm{E}}$ when the cores
start growing at orbital distances beyond 20 AU. This is contrary to what has
so far been reported in most pebble accretion models based on the dominant
size approach, where growing such massive cores when planetary embryos start
at wider orbits is difficult. Furthermore, many of these studies have
predicted that pebble accretion is either too efficient or inefficient,
depending on the prevailing physical conditions in the disc. However, in our
simulations, we may attribute such massive core sizes to the fact that low
fragmentation velocities keep grain sizes small. These small size grains
migrate to the inner disc regions slowly compared with larger grains and may
last longer in many parts of the disc, thereby promoting core growth
especially large orbital distances.

In our simulations, the growth timescales to reach pebble isolation mass
typically span from $\sim$ 0.2 to 1.8 Myr as shown in the bottom panel of
Figure 4.1. This depends on the initial location of embryo and the dust-to-gas
ratio, where in the inner disc regions the cores take shorter time to grow to
their isolation mass. This is because the pebble isolation mass as found by
L2014 and B2018 is a cubic function of the disc aspect ratio, which increases
with orbital distances. Thus, cores that grow at far orbital distances require
more time and material to reach isolation mass (see Bitsch et al., 2015, for a
more detailed discussion).

###  4.2 Core growth with full grain size distribution

Figure 4.2: The final core masses (top panel) and the total growth time
(bottom panel) as a function of the starting position, resulting from
concurrent accretion of multiple pebble species considering full pebble size
distribution. The plots have the same meaning as in Figure 4.1. Here the same
single species model from Figure 4.1 is overplotted for comparison with
concurrent accretion of full size distribution model. The planets in these
plots migrate and their growth tracks are shown in Figure 4.5. Figure 4.3:
Grain size distribution reconstructed after 0.5 Myr (top row), 1 Myr (middle
row) and 2 Myr (bottom row) of disc evolution for nominal $f_{\rm{DG}}=0.01$.
We derived the grain size distributions using $\alpha_{\rm{t}}=10^{-3}$ with
$u_{\rm{f}}$ = 10 m/s (solid line), $\alpha_{\rm{t}}=10^{-4}$ with
$u_{\rm{f}}$ = 10 m/s (dashed line) and $\alpha_{\rm{t}}=10^{-4}$ with
$u_{\rm{f}}$ = 1 m/s (dotted line). The Stokes numbers corresponding to the
grain sizes are indicated on the upper axis. The grain sizes were
reconstructed from the two-population model of Birnstiel et al. (2012) where
dust evolution is governed by growth, fragmentation and drift limits. Here,
the large population carries 0.75 and 0.97 of the solid mass in the
fragmentation and drift limits, respectively (for details, see Appendix A).
The vertical lines show pebble species that would be blocked at a planetary
mass labelled on each line. The line labels are pebble isolation masses with
diffusion corresponding to $\alpha_{\mathrm{t}}=10^{-4}$. Here, all pebble
species on the right of the vertical line would be blocked. Figure 4.4: The
growth time of the cores in Figure 4.2 that grow by accretion of multi-pebble
species paradigm for nominal $f_{\rm{DG}}=0.01$, with the starting positions
of planetary embryos shown on each plot. Some planetary cores grow very fast
initially, and then the growth quickly stalls as shown by the horizontal
lines. Along the horizontal line, the planet is accreting pebble species with
the smallest Stokes numbers which have insignificant contribution to the core
mass.

In Figure 4.2, we present the planetary core masses resulting from accretion
of full dust size distribution, where the size of different dust species were
reconstructed using the grain size reconstruction recipe of Birnstiel et al.
(2015). In the simulations, we used the pebble isolation mass with diffusion
to separate pebbles species for which the isolation mass has been attained.
This ensures that when the planet has reached isolation mass of a particular
pebble species, that species is not accreted again during the rest of core
growth period. In Figure 4.2, we overplotted the results of the dominant
species model presented in Figure 4.1 for purpose of comparison with the
multi-species model.

From Figure 4.2, in the multi-pebble species approach, the core masses
increase substantially in comparison with our dominant species approach. This
is because the cores accrete not just a single but a variety of pebble
species, where growth is sustained by accretion of smaller pebbles even after
the larger pebbles have been isolated. We note here that the sustained growth
is a consequence of each pebble species being isolated at different core
masses due to the dependence of pebble isolation on pebbles Stokes number (see
Equation (2.22)). Consequently, the core takes more time to grow than in the
case of single species model. In the latter case, accretion stops as soon as
the planet just reaches its classical pebble isolation mass, and hence the
planets have shorter growth times. We shall discuss this in detail in Section
4.4.

Keeping the fragmentation velocity at $10~{}\rm{m/s}$ and reducing the
turbulence strength to $\alpha_{\rm{t}}=10^{-4}$, results into final core
masses that are much smaller than those obtained from the model with
turbulence level $\alpha_{\rm{t}}=10^{-3}$. Here, the growth patterns in both
the single and multi-pebble species models are very similar. This is because
the low turbulence level has three main effects: first, it allows larger grain
size distribution which settle more efficiently toward the mid-plane;
secondly, the larger grains drift much faster than the smaller grains produced
in a more turbulent disc. The former effect would result in better accretion
efficiency, but then this could be counteracted by fast drifting grains which
could result in significant loss of the larger grain. Thus, all in all, for
$10~{}\rm{m/s}$ and $\alpha_{\rm{t}}=10^{-4}$, the cores may accumulate much
smaller material and take a lot more time to reach pebble isolation mass in
such an environment. Thirdly, low turbulence levels reduce turbulent diffusion
of grains across the pressure bump. Hence in this case, the pebble isolation
mass with diffusion approaches the classical case without diffusion (Bitsch et
al., 2018; Ataiee et al., 2018). Here, the pebble isolation mass is smaller
than the case where turbulence strength is greater, hence resulting into
smaller core masses. This is further illustrated in Figure 4.3, where for
$\alpha_{t}=10^{-4}$, pebbles with Stokes number greater 0.01 can be isolated
by very similar planet masses (shown by red labels on the vertical lines).
This also explains why the final core masses in both the dominant and multi-
species models in the middle panel of Figure 4.2 are very similar.

With $u_{\mathrm{f}}=1~{}\rm{m/s}$ and $\alpha_{t}=10^{-4}$, we obtained
supermassive cores in the same way as the case of single species model as
shown on the right panel of Figure 4.2. To further understand this trend, we
refer to the grain size distributions shown in Figure 4.3. Firstly, as
demonstrated in Figure 4.3, a fragmentation velocity of 1 m/s keeps pebbles at
rather much smaller sizes, and so drift at lower speeds than larger pebbles
produced when the fragmentation limit is set to 10 m/s. Consequently, the
small and slow drifting pebbles live much longer (possibly everywhere in the
disc) than the larger pebbles which may drain much more quickly. This creates
opportunity for core growth at wider orbits. For example, from Figure 4.3, the
grain size distribution remains more stable after 1 Myr of disc evolution for
our disc model with $\alpha_{\rm{t}}=10^{-4}$ and $u_{\rm{f}}=1~{}\rm{m/s}$
compared with the other two models. Secondly, the low turbulence ensures that
the small pebbles are not stirred too far from the mid-plane.

From the above discussion, it implies that in the disc model with
$\alpha_{\rm{t}}=10^{-4}$ and $u_{\rm{f}}=1~{}\rm{m/s}$, pebbles are kept
within the feeding zone of the core, possibly most of the time. The cores can
then continue to accrete the small pebbles for an extended period of time.
This suggests that at low disc turbulence and low fragmentation velocities,
core growth from pebble accretion might be possible even in the last stages of
the disc lifetime.

We note that the classical dominant species approach relies on the premise
that the grain population containing the biggest solid mass budget constitutes
grains of very similar aerodynamic size. This can be true if turbulence,
fragmentation and grain growth conspire rightly to keep grains at very similar
sizes. Consequently, in such an environment where grain size distribution can
be similar, our multi-pebble species approach should yield very similar core
masses as the classical dominant species model, possibly at faster growth
rates.

In addition, our multi-species model requires two conditions: (i) the grains
need to have varied size distribution with dissimilar aerodynamic properties
and (ii) turbulence is required to operate in such a way as to enforce pebble
isolation mass dependence on pebble Stokes number. The first condition leads
to a departure from the classical dominant species model because the different
grain species are now subjected to different gas drag, which underpins pebble
accretion. The second condition is required for sustaining core growth from
smaller pebbles after the larger ones have been isolated. This holds only if
pebble isolation mass is not supposedly universal for all pebble species at a
given radial location.

In our simulations, core growth rates and the final core masses are sensitive
to dust-to-gas ratio, fragmentation velocity and turbulence strength with the
incorporation of the full grain size distributions. This gives a diverse
outcome of final core masses as shown in Figure 4.2.

It should be noted that the masses in Figure 4.3 do not directly correspond to
the masses in Figure 4.2. This is because in Figure 4.2, we show the planet
mass as a function of the starting position. On the other hand, the orbital
distances in Figure 4.3 are sample distances at which we constructed grain
size distributions, and at which we showed sample isolation masses to
illustrate at what mass pebble species may be blocked. Therefore, since we
incorporated orbital migration, we need the final orbital distances from
Figure 4.5 if we wish to relate the core masses in Figure 4.2 to the isolation
masses in Figure 4.3. Furthermore, according to our multi-species accretion
model, accretion stops once pebbles with $\tau_{\rm{i}}=0.001$ are blocked as
earlier explained in Section 3. This means that the final core mass is
determined by the pebbles with Stokes number closest to 0.001. Because this
Stokes number may be slightly greater than 0.001, the final core mass may not
necessarily correspond to the isolation mass for the Stokes number 0.001
indicated in Figure 4.3, but closer to it.

###  4.3 Growth tracks of planetary cores

Figure 4.5: The growth tracks of the cores in Figure 4.2 for nominal
$f_{\rm{DG}}=0.01$. The solid lines show the evolution of the orbital position
and the filled circles indicate the final masses and positions of the cores.
The horizontal migration tracks show that the planet is accreting pebble
species with smallest Stokes numbers and does not significantly increase in
mass, as it continues to migrate.

Figure 4.4 represents the time evolution of the planetary cores that grow via
the multi-species accretion paradigm and initially implanted at orbital
positions shown in Figure 4.2. The corresponding orbital evolution is shown in
Figure 4.5 for a nominal dust-to-gas ratio of 0.01.

In Figures 4.4 and 4.5, the points where the growth tracks flatten show the
points of growth saturation, and consequently there is no significant increase
in mass. This occurs when the majority of pebbles species have been isolated,
but the core continues to slowly accrete the remaining pebble species with the
smallest Stokes number. It is important to remember that in our simulations,
the final mass of the planet is determined when it reaches the isolation mass
of pebbles with Stokes number $\approx$ 0.001. This then sets the final core
growth time of the planets. In our results the constribution of pebbles with
small Stokes number close to 0.001 is negligibly small as reported in several
studies (e.g., Guillot et al., 2014; Johansen et al., 2019). This is shown by
the flattening of mass-time curves for the growing cores (see Figure 4.4). The
growth tracks which end at 2 Myr of disc evolution pertain to cores that have
not yet reached their pebble isolation masses.

In Figure 4.4, planetary cores accrete more efficiently in the moderately
turbulent disc where we set $u_{\rm{f}}=10~{}\rm{m/s}$ and
$\alpha_{\rm{t}}=10^{-3}$ compared with the models for which the turbulence
parameter is set to $\alpha_{\rm{t}}=10^{-4}$. The apparent reasons for the
differences in the growth times for these models are associated with the
particle size distribution regulated by growth, fragmentation and turbulence
as previously discussed in Subsection 4.2.

With $u_{\rm{f}}=10~{}\rm{m/s}$ and $\alpha_{\rm{t}}=10^{-4}$, the grain size
distribution mostly constitutes larger grains that drift faster compared with
the other two models. This means smaller material availability and hence
slower core growth rates as most of the larger grains are lost on short
timescales via radial drift. This results into a relatively slow growth and
hence longer growth times as compared with the other two models as shown in
Figure 4.4, especially for initial orbital positions between 10 – 20 AU.

The size distribution produced by $u_{\rm{f}}=1~{}\rm{m/s}$ and
$\alpha_{\rm{t}}=10^{-4}$ are typically in millimetre range. These smaller
sized dust grains are accreted less efficiently, resulting into extended
growth times before the cores reach pebble isolation mass. However, concurrent
accretion of these smaller multi-pebble species allows the planet to reach
pebble isolation mass well before the 2 Myr of disc evolution, even for wider
initial orbital locations.

Since the final core masses are determined by the pebble isolation mass for
the smaller pebbles, one would expect longer growth times because of the slow
accretion rates of these small pebbles. However, initially the accretion of
larger pebbles is more efficient than the smaller pebbles whose accretion rate
increases as the core grows bigger. Also the concurrent accretion of different
species helps the core to rapidly increase in mass in a much shorter time.
Thus the smaller pebbles can now be accreted much more efficiently, thereby
shortening the time for the planet to reach the isolation of the smallest
species.

In Figure 4.5, we show the core growth trajectories for the same sets of the
parameters as in Figures 4.4. Here, planetary cores that were able to reach
the pebble isolation mass migrated significantly inward. However, the fast
growth rates of some cores propel them to reach pebble isolation mass after
migrating over relatively short orbital distances. For example, for the models
with low turbulence level of $\alpha_{\rm{t}}=10^{-4}$, planetary cores that
start at 5 AU migrated relatively short distances (by $\approx$ 2.5 AU) by the
time they reached their isolation mass. Thus there may be a possibility that
planetary core can form almost in-situ in some parts of the disc if core
growth proceeds by concurrent accretion of multiple pebble species.

Analogous to the core growth saturation points discussed above, the points
where the curves in Figure 4.5 flatten mark the orbital distance at which the
cores reached their isolation mass. We remind the reader here that we did not
include type II migration which might significantly change the final orbits of
the planets.

In our simulations, we were interested in core growth and thus did not model
gas accretion as well as transition from type I to type II migration. Type II
migration could be triggered during gas accretion phase after the cores have
reached their full pebble isolation mass. We envision that with the inclusion
of type II migration, orbital decay is enhanced and could be significant to
the extent that planets are lost to the host star, unless slower migration
mechanisms like the dynamical torque (Paardekooper, 2014) and slower type II
migrations (as in Crida & Bitsch, 2017; Crida et al., 2017; Robert et al.,
2018; Bergez-Casalou et al., 2020; Ndugu et al., 2021) are considered.
Therefore, even though our simple multiple-species accretion model shortens
the growth time of cores, the orbital dilemma of gas giant planets migration
needs to be explored further in the accretion of multi-pebble species
scenario. We shall investigate this in our upcoming work. However, a promising
solution is the heating torque (Benítez-Llambay et al., 2015) which could
prevent loss of planets to the central star by stopping inward migration, and
this allows planet formation closer to the star.

###  4.4 The role of pebble isolation mass

Figure 4.6: The dependence of core masses on turbulence strength and
fragmentation velocity. B2018_diff and B2018_ndiff denote the pebble isolation
mass models of B2018 with and without diffusion, respectively.

The new formulation of pebble isolation mass in Equation (2.22), functionally
depends on several parameters, including the turbulence strength and pebble
Stokes number; it relates with pebble Stokes number in such way that
$M_{\rm{iso,~{}i}}\propto(1/\tau_{\rm{i}})$. As a consequence of Equation
(2.22), we assume in our model that different pebble species are sequentially
isolated by the growing planet when it starts to exert a pressure bump. This
means that in our model, pebbles with large Stokes numbers are trapped
exterior the planet’s orbit first, and those with small Stokes numbers are
blocked last. We remind the reader that we did not model pressure bump due to
the growing planet. This would require modelling a rather complex gap profile
that is anchored self-consitently to viscous disc evolution as in Dullemond et
al. (2018), which is out of scope of this study and will be considered in
future works.

As we illustrated in Figures 2.2 and 4.3, pebbles may likely be blocked
preferentially according to their sizes as a consequence of the pebble
isolation prescription of B2018. For example, from both Figures 2.2 and 4.3,
all pebble species to the right of the vertical line are blocked at the
planetary mass indicated on the line.

From the above notion, the planet then reaches different isolation masses
corresponding to each pebble species. At each pebble isolation stage, the
planet continues to accrete smaller pebbles that are able to diffuse through
the pressure bump that it generates, resulting in sustained core growth.
However, when most of the pebbles have been blocked, the pebble accretion rate
drops as the planet continues to slowly accrete the remaining small pebble
species. For example, from Figure 4.5, for the parameter set of
$u_{\rm{f}}=1~{}\rm{m/s}$ and $\alpha_{\rm{t}}=10^{-4}$, a planetary seed
planted at 50 AU reaches its final core mass of $\sim 40~{}M_{\rm{E}}$ at
around 20 AU, at which pebble isolation mass ranges from 33 – 47 $M_{\rm{E}}$
for pebbles with Stokes number in the range 0.001 – 0.1 (see Figure 4.3).
Before that, as seen in Figure 4.5, the accretion rate drops at around 30 AU
without any significant increase in mass. This is because from 30 AU down to
20 AU the planet is now slowly accreting small pebble species whose isolation
mass is still high, while majority of the other pebble species have been
blocked. Finally, the planet reaches an isolation mass that roughly
corresponds to isolation mass of $\sim 40~{}M_{\rm{E}}$ for approximately 0.02
– 0.1 cm pebble species at 20 AU as shown in Figure 4.3.

As shown in Figure 4.3, pebble isolation mass spans a wide range of values. We
observe here that relatively large pebble species in the distribution need
very similar pebble isolation mass in order to be blocked. On the other hand,
the millimetre-sized pebbles have a broader range of pebble isolation mass
needed to block each species. This suggests that a planet may reach pebble
isolation mass only once if it is accreting predominantly large-size pebble
species. On the contrary, the planet may need to reach different isolation
mass corresponding to each pebble species if the distribution is dominated by
small-size pebble species.

The distributions in Figure 4.3 further suggest that centimetre-sized pebbles
produce relatively smaller cores since they are isolated at relatively low
planetary masses compared with the smaller millimetre-sized pebbles. Thus,
this suggests that building larger cores requires the presence of small-size
pebbles, typically in the millimetre range, but more growth time will be
needed to reach pebble isolation mass (see Figure 4.4). This would require
either moderately turbulent discs with high fragmentation velocities or
quiescent discs with low fragmentation velocities, both of which produce
pebbles in millimetre-centimetre ranges.

Furthermore, building massive cores requires that we implant planetary cores
at far orbital distances as can be seen in Figure 4.5. This enables the planet
to continue accreting before it reaches pebble isolation mass, which scales
not only with pebble size but also orbital distance as in Figures 4.3.
However, there appears a problem of fast planet migration that mostly drives
the planet to the inner discs, unless type I migration rates are slowed via
mechanisms like dynamical torques as in Ndugu et al. (2021). As demonstrated
here, fast core growth via multi-pebble accretion paradigm could also provide
another pathway for the planet to reach high pebble isolation mass before it
migrates rapidly to inner disc regions where the isolation mass is smaller.
Such fast core growths can allow the planet to beat the rapid type I migration
as suggested in Johansen et al. (2019).

The combination of millimetre-centimetre size distribution thus has two main
advantages: the centimetre-size species initially offer high accretion rates
that enables the planet to rapidly become massive enough to accrete the
smaller species more efficiently; the small species survive much longer in the
disc and hence steady supply of material. This can then result in massive
cores if multi-pebble accretion is taken into account together with diffusion
as a determinant for pebble isolation mass.

The influence of turbulent diffusion in the context of multi-pebble accretion
model is demonstrated in Figure 4.6. From the top panel of Figure 4.6, there
is a significant difference in final core masses resulting from models with
and without diffusion, especially in disc regions inside 20 AU. However, at
wider initial growth locations, the planet has the same mass regardless of
pebble isolation mass model used as shown in the top panel of Figure 4.6. In
this case, the grain sizes may have similar isolation mass as previous
discussed.

In the case of $u_{\mathrm{f}}=10~{}\mathrm{m/s}$ and
$\alpha_{\mathrm{t}}=10^{-4}$, as shown in the middle panel of Figure 4.6, the
difference between the pebble isolation mass schemes with and without
diffusion is marginal. This is because, as already mentioned before, grain
sizes are generally bigger at high fragmentation velocities, and with low
turbulence, such large pebbles may not easily diffuse across pressure bump
generated by the planet. This means the dependence of the final core mass on
isolation mass for larger pebbles species may turn out to be inconsequential.
However, for the same turbulence strength of $\alpha_{\mathrm{t}}=10^{-4}$ and
a lower fragmentation velocity of for $u_{\mathrm{f}}=1~{}\mathrm{m/s}$, the
effect of diffusion of pebble isolation mass is still manifested, as shown in
the bottom panel of Figure 4.6. This is because low fragmentation velocities
keep grain sizes small, which may be able to diffuse much more easily in less
turbulent disc compared with the larger grains.

In consideration of sequential isolation of pebbles as discussed above, as
more and more pebbles are trapped outside the planet’s orbit, the dust-to-gas
ratio and pebble flux would ultimately reduce in the inner disc. This would
affect growth of planetary cores interior to the planet’s orbit. In the outer
orbit of the planet, the pile up of pebbles enhances dust-to-gas ratio at the
pressure bump. This could trigger streaming instability and subsequent
formation of planetesimals, thus reducing the amount of drifting pebbles,
including those that would possibly diffuse through the pressure bump. This
can then lead to reduced multi-pebble accretion rates and hence low mass
cores.

If the idea that smaller pebbles may overcome the pressure bump holds, then
their conversion into planetesimals at the pressure bump together with the
other species would affect the accretion rate of the smaller pebbles. But the
planetesimals formed may be potential targets for accretion by the planet.
However, in our simulations we ignored the possible conversion of pebbles into
planetesimals at the pressure bump, which may play a key role in the mass
budget of the inward drifting pebbles. Here we instead focused on how the
diffusion of smaller pebbles affects core growth rates.

In the light of the novel flow isolation mass (Rosenthal & Murray-Clay, 2020),
we envisage that our results could substantially change. In particular, the
flow isolation mass would limit the accretion of smaller, tightly coupled
pebbles as they begin to interact with the planet’s atmosphere, and may simply
flow past the planet without being accreted. Consequently, the final core
masses would then be restricted to the isolation mass of the larger pebbles,
which could be lower than the values reported here. Future comprehensive
planet population synthesis studies should consider detailed comparison
between the existing pebble isolation criterion and the flow mass paradigm in
explaining the existing features of the observed gas giant planets.

###  4.5 Giant planet core formation

Our results suggest that it is very challenging to form cores of giant planets
when pebbles are much more coupled to the gas, for instance with
$\tau_{\rm{i}}\lesssim 0.001$, in agreement with Johansen et al. (2019). Such
pebbles tend to follow gas streamlines and interact gravitationally poorly
with the protoplanets, making their capture difficult (see Guillot et al.,
2014, for pebble capturing efficiencies by protoplanets). Furthermore, B2018
have demonstrated that reaching pebble isolation mass for tightly coupled
pebbles requires cores as massive as 50 $M_{\rm{E}}$. Reaching such a high
core mass may prove difficult in the framework of dominant species, whose
performance we have not rigorously tested here.

The Jupiter Near-polar Orbiter (JUNO) mission provided precise measurements of
Jupiter’s gravitational field, which has enabled better estimates of core
masses. For instance, Jupiter’s core mass was estimated to be 7 – 25
$M_{\rm{E}}$ (Wahl et al., 2017) while Debras & Chabrier (2019) estimated 25 –
30 $M_{\rm{E}}$ and 30 – 45 $M_{\rm{E}}$ for non-compact and compact Jupiter
cores, respectively. Other studies have inferred similar core masses for solar
system gas giants: Jupiter (10 –40 $M_{\rm{E}}$ (Guillot, 1999) and 37
$M_{\rm{E}}$ (Thorngren et al., 2016)); Saturn (20 –30 $M_{\rm{E}}$ (Guillot,
1999) and 27 $M_{\rm{E}}$ (Thorngren et al., 2016)). Many models have
considered giant impacts to account for Jupiter’s highly enriched core with
heavy elements (e.g., Liu et al., 2019; Ginzburg & Chiang, 2020). We have
tested that it is quite challenging to obtain such high core masses using the
classical pebble accretion model adopted in this work, unless we at least
unphysically increase pebble concentration and consider discs with longer life
times. Moreover, the cores of giant planets that migrated to $<10$ AU could
have started their growth beyond between 15 – 30 AU (Johansen et al., 2019).
In this orbital domain, we found it difficult to grow cores to masses above 25
– 30 $M_{\rm{E}}$ without evoking accretion of multiple pebble species. We
therefore think that through our multi-pebble accretion scheme, massive cores
of gas planets that match the results the JUNO mission can form.

##  5 Conclusion

In this work, we have extended the pebble accretion paradigm, where core
growth proceeds by concurrent accretion of multiple pebble species. We took
into account the dependence of pebble isolation mass on turbulence parameter
and pebble size. In our model, the final planetary core masses were set by the
planet mass needed to block the smallest sized pebble species in the
distribution of grain sizes under consideration. We self-consistently
reconstructed grain distribution throughout the core growth process using the
reconstruction model of Birnstiel et al. (2015). For the dust evolution, we
employed the two-population model of Birnstiel et al. (2012).

Our simulations yielded a diversity of planetary cores with a wide range of
core masses. The final outcome of core masses is primarily dictated by a
combination of dust-to-gas ratio, fragmentation velocities and turbulence
strength that underpin the size distributions. Under a suitable set of disc
conditions, cores of gas giants can form at orbits as far out as 50 AU.

Our work has several other important consequences. Firstly, recent studies
(Wahl et al., 2017; Debras & Chabrier, 2019) modelled Jupiter’s core to
contain heavy elements totalling to 25 – 45 $M_{\rm{E}}$, adding a further
constraint on core accretion models. It is a difficult to invoke pebble
accretion to explain such massive core build up. Moreover, within the
limitations of our numerical simulations, we could not easily build massive
cores well above $\sim 30~{}M_{\rm{E}}$ using the classical scenario of
dominant species. If we are to invoke pebble accretion to explain assembly of
such massive cores, then our model of concurrent accretion of multi-pebble
species would provide a possible mechanism to explain this enigma.

Secondly, in the study of metallicity correlation of extrasolar planets,
Guillot et al. (2006) found that their sample planets contained $\sim$ 20 –
100 $M_{\rm{E}}$ in heavy elements. The work of Thorngren et al. (2016) also
suggests a strong correlation between the planetary mass and amount of heavy
elements accreted by the planet. If these results are confirmed, then the
accumulation of such large amount of heavy material is in conflict with
classical core accretion model of planet formation. However, some recent
solutions have been proposed, for example, the merger of giant planets via
collisions (Liu et al., 2019; Ginzburg & Chiang, 2020) or the accretion of the
evaporated material may increase the heavy element content of the giant
planets (Schneider & Bitsch, 2021).

Another important aspect of classical core accretion that we did not account
for in this study is the accretion of planetesimals. Although gas giant planet
cores form with difficulties via planetesimal accretion, the capture of
planetesimals by the planet can also enrich the planet with heavy elements
(Shibata & Ikoma, 2019; Shibata et al., 2020; Venturini & Helled, 2020).
Nevertheless, as we have demonstrated here, it is possible that these planets
consumed several different dust species if their formation proceeded through
pebble accretion paradigm, and can form such massive gas giant planet cores
even at the disc’s outskirts.

Thirdly, a diversity of gaps and rings over a wide orbital regime has been
discovered in many discs through the ALMA program (e.g., Huang et al., 2018a,
b, c; Long et al., 2018) and the DSHARP campaign (e.g., Andrews et al., 2018).
Planet-disc interactions have been evoked as one possible way to explain the
occurrence of these structures, though this is subject to discussion (Ndugu et
al., 2019). Recently, in a sample of 16 ALMA discs studied by van der Marel et
al. (2019), the authors found some substructures at orbital distances far
beyond 30 AU, and found evidence for presence of massive planets at those
orbital distances. An important question as to whether or not these
substructures are indeed caused by growing planets has been put forth (e.g.,
Lodato et al., 2019; Ndugu et al., 2019). If it is true that planets may be
responsible for opening gaps at such wider orbital distances, then the
existing core accretion models are missing important ingredients for
explaining the formation of planets at these remote locations. It is, however,
possible to form gas giant planets at these wide orbits if: (i) concurrent
accretion of multiple pebble species is accounted for, (ii) it is assumed that
the giant planets at such wide orbits form via gravitational instability
paradigm (Boss, 1997; Boley, 2009; Armitage, 2010) and (iii) assuming the
observed rings are caused by other phenomena such as MRI other than an on-
going planet formation. The rings could in turn be hot spots for planet
formation at such wide orbits (Morbidelli, 2020). Thus the precise explanation
for the origin of the observed gaps/rings at wider orbits is a complicated
one, more like a Chicken-Egg problem.

Although our model has shed some light on the formation of massive giant
planets at wider orbits of discs, it did not capture important planet
formation aspects like competition of multiple cores for the available
building block and orbital manipulation by N-body interactions between cores.
We therefore recommend future models that study multi-pebble species accretion
to consider core growth competition and the N-body paradigm for detailed and
substantial explanation of the formation of gas giant planets at wider orbits.

## Acknowledgments

We thank the anonymous referee for the insightful reports that helped us to
improve this paper. G.A. further thanks Michiel Lambrechts for constructive
and inspirational comments. We thank Swedish International Development
Cooperation Agency (SIDA) for financial support through International Science
Programme (ISP) Uppsala University Sweden to the East Africa Astronomical
Research Network (EAARN).

## Data availability

For the purpose of reproducibility, the code used to obtain results in this
paper will be provided upon reasonable request.

## References

  * Adachi et al. (1976) Adachi I., Hayashi C., Nakazawa K., 1976, Progress of Theoretical Physics, 56, 1756
  * Alexander & Armitage (2007) Alexander R. D., Armitage P. J., 2007, MNRAS, 375, 500
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Ansdell et al. (2017) Ansdell M., Williams J. P., Manara C. F., Miotello A., Facchini S., van der Marel N., Testi L., van Dishoeck E. F., 2017, AJ, 153, 240
  * Armitage (2010) Armitage P., 2010, Astrophysics of Planet Formation. Cambridge University Press
  * Ataiee et al. (2018) Ataiee S., Baruteau C., Alibert Y., Benz W., 2018, A&A, 615, A110
  * Barrière-Fouchet et al. (2005) Barrière-Fouchet L., Gonzalez J. F., Murray J. R., Humble R. J., Maddison S. T., 2005, A&A, 443, 185
  * Baruteau et al. (2014) Baruteau C., et al., 2014, Protostars and Planets VI, pp 667–689
  * Benítez-Llambay et al. (2015) Benítez-Llambay P., Masset F., Koenigsberger G., Szulágyi J., 2015, Nature, 520, 63
  * Bergez-Casalou et al. (2020) Bergez-Casalou C., Bitsch B., Pierens A., Crida A., Raymond S. N., 2020, A&A, 643, A133
  * Birnstiel et al. (2009) Birnstiel T., Dullemond C. P., Brauer F., 2009, A&A, 503, L5
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Birnstiel et al. (2011) Birnstiel T., Ormel C. W., Dullemond C. P., 2011, A&A, 525, A11
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Birnstiel et al. (2015) Birnstiel T., Andrews S. M., Pinilla P., Kama M., 2015, The Astrophysical Journal, 813, L14
  * Bitsch & Johansen (2017) Bitsch B., Johansen A., 2017. Springer International Publishing, Cham, pp 339–366, doi:10.1007/978-3-319-60609-5_12
  * Bitsch et al. (2015) Bitsch B., Lambrechts M., Johansen A., 2015, A&A, 582, A112
  * Bitsch et al. (2018) Bitsch B., Morbidelli A., Johansen A., Lega E., Lambrechts M., Crida A., 2018, A&A, 612, A30
  * Blum & Wurm (2008) Blum J., Wurm G., 2008, ARA&A, 46, 21
  * Boley (2009) Boley A. C., 2009, ApJ, 695, L53
  * Boss (1997) Boss A. P., 1997, Science, 276, 1836
  * Bottke et al. (2005a) Bottke W. F., Durda D. D., Nesvorný D., Jedicke R., Morbidelli A., Vokrouhlický D., Levison H., 2005a, Icarus, 175, 111
  * Bottke et al. (2005b) Bottke W. F., Durda D. D., Nesvorný D., Jedicke R., Morbidelli A., Vokrouhlický D., Levison H. F., 2005b, Icarus, 179, 63
  * Brauer et al. (2007) Brauer F., Dullemond C. P., Johansen A., Henning T., Klahr H., Natta A., 2007, A&A, 469, 1169
  * Brauer et al. (2008) Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859
  * Brügger et al. (2018) Brügger N., Alibert Y., Ataiee S., Benz W., 2018, A&A, 619, A174
  * Cameron (1978) Cameron A. G. W., 1978, Moon and Planets, 18, 5
  * Carrasco-González et al. (2016) Carrasco-González C., et al., 2016, ApJ, 821, L16
  * Carrera et al. (2015) Carrera D., Johansen A., Davies M. B., 2015, A&A, 579, A43
  * Coleman & Nelson (2014) Coleman G. A. L., Nelson R. P., 2014, Monthly Notices of the Royal Astronomical Society, 445, 479
  * Crida & Bitsch (2017) Crida A., Bitsch B., 2017, Icarus, 285, 145
  * Crida et al. (2017) Crida A., Bitsch B., Ndugu N., Morbidelli A., 2017, in European Planetary Science Congress. pp EPSC2017–44
  * Debras & Chabrier (2019) Debras F., Chabrier G., 2019, ApJ, 872, 100
  * Dominik & Dullemond (2008) Dominik C., Dullemond C. P., 2008, A&A, 491, 663
  * Drążkowska et al. (2019) Drążkowska J., Li S., Birnstiel T., Stammler S. M., Li H., 2019, ApJ, 885, 91
  * Drążkowska et al. (2021) Drążkowska J., Stammler S. M., Birnstiel T., 2021, arXiv e-prints, p. arXiv:2101.01728
  * Dullemond & Dominik (2004) Dullemond C. P., Dominik C., 2004, A&A, 421, 1075
  * Dullemond & Dominik (2005) Dullemond C. P., Dominik C., 2005, A&A, 434, 971
  * Dullemond et al. (2018) Dullemond C. P., et al., 2018, ApJ, 869, L46
  * Durisen et al. (2007) Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T., Rice W. K. M., 2007, Protostars and Planets V, pp 607–622
  * Dzyurkevich et al. (2010) Dzyurkevich N., Flock M., Turner N. J., Klahr H., Henning T., 2010, A&A, 515, A70
  * Gammie (2001) Gammie C. F., 2001, The Astrophysical Journal, 553, 174
  * Garaud et al. (2004) Garaud P., Barrière-Fouchet L., Lin D. N. C., 2004, ApJ, 603, 292
  * Ginzburg & Chiang (2020) Ginzburg S., Chiang E., 2020, MNRAS, 498, 680
  * Guilera et al. (2020) Guilera O. M., Sándor Z., Ronco M. P., Venturini J., Miller Bertolami M. M., 2020, A&A, 642, A140
  * Guillot (1999) Guillot T., 1999, Planetary and Space Science, 47, 1183
  * Guillot et al. (2006) Guillot T., Santos N. C., Pont F., Iro N., Melo C., Ribas I., 2006, A&A, 453, L21
  * Guillot et al. (2014) Guillot T., Ida S., Ormel C. W., 2014, A&A, 572, A72
  * Gundlach & Blum (2015) Gundlach B., Blum J., 2015, ApJ, 798, 34
  * Gundlach et al. (2011) Gundlach B., Kilias S., Beitz E., Blum J., 2011, Icarus, 214, 717
  * Huang et al. (2018a) Huang J., et al., 2018a, The Astrophysical Journal, 852, 122
  * Huang et al. (2018b) Huang J., et al., 2018b, ApJ, 869, L42
  * Huang et al. (2018c) Huang J., et al., 2018c, ApJ, 869, L43
  * Johansen & Bitsch (2019) Johansen A., Bitsch B., 2019, A&A, 631, A70
  * Johansen & Lacerda (2010) Johansen A., Lacerda P., 2010, MNRAS, 404, 475
  * Johansen et al. (2007) Johansen A., Oishi J. S., Mac Low M.-M., Klahr H., Henning T., Youdin A., 2007, Nature, 448, 1022
  * Johansen et al. (2009) Johansen A., Youdin A., Klahr H., 2009, ApJ, 697, 1269
  * Johansen et al. (2011) Johansen A., Klahr H., Henning T., 2011, A&A, 529, A62
  * Johansen et al. (2014) Johansen A., Blum J., Tanaka H., Ormel C., Bizzarro M., Rickman H., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 547 (arXiv:1402.1344), doi:10.2458/azu_uapress_9780816531240-ch024
  * Johansen et al. (2019) Johansen A., Ida S., Brasser R., 2019, A&A, 622, A202
  * Kokubo & Ida (1998) Kokubo E., Ida S., 1998, Icarus, 131, 171
  * Kokubo & Ida (2012) Kokubo E., Ida S., 2012, Progress of Theoretical and Experimental Physics, 2012, 01A308
  * Kuiper (1951) Kuiper G. P., 1951, pnas, 37, 1
  * Lambrechts & Johansen (2012) Lambrechts M., Johansen A., 2012, A&A, 544, A32
  * Lambrechts & Johansen (2014) Lambrechts M., Johansen A., 2014, A&A, 572, A107
  * Lambrechts et al. (2014) Lambrechts M., Johansen A., Morbidelli A., 2014, A&A, 572, A35
  * Levison et al. (2010) Levison H. F., Thommes E., Duncan M. J., 2010, AJ, 139, 1297
  * Liu et al. (2019) Liu S.-F., Hori Y., Müller S., Zheng X., Helled R., Lin D., Isella A., 2019, Nature, 572, 355
  * Lodato et al. (2019) Lodato G., et al., 2019, MNRAS, 486, 453
  * Long et al. (2018) Long F., et al., 2018, ApJ, 869, 17
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Manara et al. (2019) Manara C. F., Mordasini C., Testi L., Williams J. P., Miotello A., Lodato G., Emsenhuber A., 2019, A&A, 631, L2
  * Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  * Mizuno (1980) Mizuno H., 1980, Progress of Theoretical Physics, 64, 544
  * Mizuno et al. (1978) Mizuno H., Nakazawa K., Hayashi C., 1978, Progress of Theoretical Physics, 60, 699
  * Morbidelli (2020) Morbidelli A., 2020, A&A, 638, A1
  * Morbidelli & Nesvorny (2012) Morbidelli A., Nesvorny D., 2012, A&A, 546, A18
  * Morbidelli et al. (2009) Morbidelli A., Bottke W. F., Nesvorný D., Levison H. F., 2009, Icarus, 204, 558
  * Morbidelli et al. (2015) Morbidelli A., Lambrechts M., Jacobson S., Bitsch B., 2015, Icarus, 258, 418
  * Mordasini et al. (2009a) Mordasini C., Alibert Y., Benz W., 2009a, A&A, 501, 1139
  * Mordasini et al. (2009b) Mordasini C., Alibert Y., Benz W., Naef D., 2009b, A&A, 501, 1161
  * Musiolik & Wurm (2019) Musiolik G., Wurm G., 2019, ApJ, 873, 58
  * Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
  * Nayakshin (2010) Nayakshin S., 2010, MNRAS, 408, L36
  * Nayakshin et al. (2019) Nayakshin S., Dipierro G., Szulágyi J., 2019, MNRAS, 488, L12
  * Ndugu et al. (2018) Ndugu N., Jurua E., Bitsch B., 2018, MNRAS, 474, 886
  * Ndugu et al. (2019) Ndugu N., Bitsch B., Jurua E., 2019, MNRAS, 488, 3625
  * Ndugu et al. (2021) Ndugu N., Bitsch B., Morbidelli A., Crida A., Jurua E., 2021, MNRAS, 501, 2017
  * Ormel & Klahr (2010) Ormel C. W., Klahr H. H., 2010, A&A, 520, A43
  * Paardekooper (2014) Paardekooper S.-J., 2014, MNRAS, 444, 2031
  * Paardekooper & Mellema (2006) Paardekooper S.-J., Mellema G., 2006, A&A, 459, L17
  * Paardekooper et al. (2010) Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  * Paardekooper et al. (2011) Paardekooper S.-J., Baruteau C., Meru F., 2011, MNRAS, 416, L65
  * Pérez et al. (2012) Pérez L. M., et al., 2012, ApJ, 760, L17
  * Pinilla et al. (2012) Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L., Testi L., Natta A., 2012, A&A, 538, A114
  * Raettig et al. (2015) Raettig N., Klahr H., Lyra W., 2015, ApJ, 804, 35
  * Rafikov (2005) Rafikov R. R., 2005, ApJ, 621, L69
  * Rice et al. (2003) Rice W. K. M., Armitage P. J., Bonnell I. A., Bate M. R., Jeffers S. V., Vine S. G., 2003, MNRAS, 346, L36
  * Robert et al. (2018) Robert C. M. T., Crida A., Lega E., Méheut H., Morbidelli A., 2018, A&A, 617, A98
  * Rodmann et al. (2006) Rodmann J., Henning T., Chandler C. J., Mundy L. G., Wilner D. J., 2006, A&A, 446, 211
  * Rosenthal & Murray-Clay (2020) Rosenthal M. M., Murray-Clay R. A., 2020, ApJ, 898, 108
  * Safronov (1969) Safronov V., 1969, Evolution of Protoplanetary Cloud and Formation of the Earth and the Planets
  * Savvidou & Bitsch (2021) Savvidou S., Bitsch B., 2021, arXiv e-prints, p. arXiv:2103.11995
  * Savvidou et al. (2020) Savvidou S., Bitsch B., Lambrechts M., 2020, A&A, 640, A63
  * Schneider & Bitsch (2021) Schneider A. D., Bitsch B., 2021, arXiv e-prints, p. arXiv:2105.13267
  * Shibata & Ikoma (2019) Shibata S., Ikoma M., 2019, MNRAS, 487, 4510
  * Shibata et al. (2020) Shibata S., Helled R., Ikoma M., 2020, A&A, 633, A33
  * Singer et al. (2019) Singer K. N., et al., 2019, Science, 363, 955
  * Takeuchi & Lin (2005) Takeuchi T., Lin D. N. C., 2005, ApJ, 623, 482
  * Tanaka & Ida (1999) Tanaka H., Ida S., 1999, Icarus, 139, 350
  * Tanga et al. (2004) Tanga P., Weidenschilling S. J., Michel P., Richardson D. C., 2004, A&A, 427, 1105
  * Testi et al. (2003) Testi L., Natta A., Shepherd D. S., Wilner D. J., 2003, A&A, 403, 323
  * Testi et al. (2014) Testi L., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 339 (arXiv:1402.1354), doi:10.2458/azu_uapress_9780816531240-ch015
  * Thommes et al. (2003) Thommes E. W., Duncan M. J., Levison H. F., 2003, Icarus, 161, 431
  * Thorngren et al. (2016) Thorngren D. P., Fortney J. J., Murray-Clay R. A., Lopez E. D., 2016, ApJ, 831, 64
  * Trotta et al. (2013) Trotta F., Testi L., Natta A., Isella A., Ricci L., 2013, A&A, 558, A64
  * Tychoniec et al. (2020) Tychoniec Ł., et al., 2020, A&A, 640, A19
  * Uribe et al. (2011) Uribe A. L., Klahr H., Flock M., Henning T., 2011, ApJ, 736, 85
  * Venturini & Helled (2020) Venturini J., Helled R., 2020, A&A, 634, A31
  * Venturini et al. (2020) Venturini J., Guilera O. M., Ronco M. P., Mordasini C., 2020, A&A, 644, A174
  * Wada et al. (2008) Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2008, ApJ, 677, 1296
  * Wahl et al. (2017) Wahl S. M., et al., 2017, Geophys. Res. Lett., 44, 4649
  * Weber et al. (2018) Weber P., Benítez-Llambay P., Gressel O., Krapp L., Pessah M. E., 2018, ApJ, 854, 153
  * Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
  * Wetherill (1980) Wetherill G. W., 1980, ARA&A, 18, 77
  * Whipple (1972) Whipple F. L., 1972, in Elvius A., ed., From Plasma to Planet. p. 211
  * Wilner et al. (2005) Wilner D. J., D’Alessio P., Calvet N., Claussen M. J., Hartmann L., 2005, ApJ, 626, L109
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588
  * van der Marel et al. (2019) van der Marel N., Dong R., di Francesco J., Williams J. P., Tobin J., 2019, The Astrophysical Journal, 872, 112

##  Appendix A The dust distributions and mass fractions

Figure A.1: The grain surface density and Stokes number as a function of grain
size reconstructed from the two-population model at 5 AU, 10 AU, 20 AU and 50
AU at time 0.1 Myr, 0.5 Myr, 1 Myr and 2 Myr. On the plots, we show the
fractions of small and large grains that have Stokes numbers less and greater
than 0.001, respectively. The shaded regions show large grains with Stokes
number 0.001, that were mostly accreted in our growth model. Here, 2pop is the
fraction of the large dust population from the two-population simulations that
contains most of the solid mass. The dust evolution was performed with
$u_{\mathrm{f}}=10\mathrm{m/s}$, $\alpha_{\mathrm{t}}=10^{-3}$ and
$f_{\mathrm{DG}}=0.01$ for 2 Myrs. Figure A.2: The radial mass fractions of
small (with Stokes number < 0.001) and large ( with Stokes number $\geq
0.001$) grains. $\Sigma_{0}$ and $\Sigma_{1}$ are the small and large dust
populations from the two-dust population model. Each row shows simulations
with different combinations of fragmentation velocity and turbulence strength,
and initial dust-to-gas ratio of $f_{\mathrm{DG}}=0.01$.

In Figure A.1, we show the variation of pebble surface density and the Stokes
number as a function of grain for the dust size distribution reconstructed at
radial distances 5 AU, 10 AU, 20 AU and 50 AU, for
$\alpha_{\mathrm{t}}=10^{-3}$, u${}_{\mathrm{f}}=10~{}\mathrm{m/s}$ and
f${}_{\mathrm{DG}}=0.01$. On each plot, we show the fractions of pebbles that
have Stokes number greater and less than 0.001. Grains with Stokes number
$\geq 0.001$ are shown by the shaded region. We also show the fractions of the
large and small populations from the two-population model. In the two-
population scheme, grain sizes are calculated and fixed by growth,
fragmentation and drift limits. Here, in the inner and outer disc regions, the
fragmentation and drift limits dominate, respectively with 0.75 and 0.97 of
the mass constituting the large population. As shown in the plots, the
reconstructed grains distributions with Stokes number $\geq 0.001$ dominate in
the inner disc regions, $\leq$ 30 AU, which carry over 80% of the mass,
similar to the fractions of the two dust population model.

Furthermore, as shown in Figure A.1, very few grains have Stokes number
greater 0.1, and most of the mass in reconstructed grain distribution is
carried by the population of grains with Stokes number $\geq 0.001$. In
addition, the accretion of pebbles with Stokes number < 0.001 is generally
inefficient since these pebbles are usually well coupled to the gas (Guillot
et al., 2014), and we observed the same trend in our simulations. Thus in our
simulations, the embryos mostly accrete pebbles with Stokes number $\geq
0.001$.

In Figure A.2, we show dust fractions as a function of radial distance for the
different combinations of fragmentation velocities and turbulence strengths,
and for both the simulated two-population and the reconstructed grain size
distributions. In the simulations with fragmentation velocity of 10 m/s, most
of the mass is carried by grains with Stokes number $\geq 0.001$ within 50 AU
as shown by the panels in the first and second rows of Figure A.2. On the
other hand, for the case of fragmentation velocity of 1 m/s, the fraction of
pebbles with Stokes number $\geq 0.001$ begins to approach that of the large
population simulated using the two dust model at orbital distances greater
than 5 AU. Here, within 5 AU, pebbles with Stokes number less than 0.001
dominate. This may suggest a stronger fragmentation inside 5 AU that produces
smaller grains, which then have small Stokes numbers compared with other parts
of the disc.

##  Appendix B The recipe for concurrent accretion of multi-pebble species

In this Section we give a brief description of the recipe that we used in our
numerical simulations involving concurrent accretion of multi-pebble species
in reference to the full-size distribution as in Birnstiel et al. (2012). It
is a simple extension of the implementation of single species model that can
readily be adopted for other models of particle distribution. At each time
snapshot and orbital distance referenced by j and k respectively, we implement
core accretion of multiple pebble species as follows:

  1. Step 1.

Initialise planet mass, $M_{\rm{j=0,k=0}}$, orbital distance, $r_{\rm{k=0}}$.

  2. Step 2.

Calculate the disc parameters at rk and time tj:

  3. Step 3.

Obtain the logarithmic distribution of N pebble sizes
$R_{\rm{i+1}}=1.12R_{\rm{i}}$ as in Birnstiel et al. (2011).

  4. Step 4.

From the logarithmic particle size distribution, calculate the pebble surface
density, $\Sigma_{\rm{peb,i}}$, for each of the N species according the
reconstruction scheme in Birnstiel et al. (2015) and calculate the Stokes
number, $\tau_{\rm{i}}$, for each species in Epstein regime as

$\tau_{\rm{i}}=\frac{\rho_{\bullet}\pi R_{\rm{i}}}{2\Sigma_{g}}.$ (B.1)

  5. Step 5.

Calculate pebble isolation mass for each species according to Equation 2.22.

  6. Step 6.

If the planet mass is larger than the pebble isolation mass of pebble species
with the smallest Stokes number in the size distribution, stop time
integration. This is because as the planet grows, pebbles with larger Stokes
numbers are isolated first, while the ones with smaller Stokes numbers are
isolated last, depending on the mass of the planet (see Equation (2.22) and
more explanation in B2018). Hence the core mass is determined by the isolation
mass of pebbles with smallest Stokes numbers.

  7. Step 7.

If planet mass is larger than the isolation mass of the i-th species, stop
accretion of the i-th species. Otherwise accrete the i-th species by
calculating the core accretion rate $\dot{M}_{\rm{core,i}}$ using Equations
(2.14), (2.15) and (2.17).

  8. Step 8.

Calculate the total core growth rate by consolidating contributions from each
species:

$\dot{M}_{\rm{j,k}}=\sum_{\rm{i}}\dot{M}_{\rm{core,i}}.$ (B.2)

  9. Step 9.

Calculate the new planet mass from

$M_{\rm{j+1,k}}=M_{\rm{j,k}}+\dot{M}_{\rm{j,k}}\delta t,$ (B.3)

where $\delta t$ is the timestep.

  10. Step 10.

Calculate the migration rate and the new orbital location of the growing core.

  11. Step 11.

Increase time by $\delta t$ and repeat from Step 2 until the condition in Step
6 is fulfilled.