# Efficient planet formation by pebble accretion in ALMA rings

Haochang Jiang (蒋昊昌)1, Chris W. Ormel1
1Department of Astronomy, Tsinghua University, Haidian DS 100084, Beijing,
China E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In the past decade, ALMA observations have revealed that a large fraction of
protoplanetary discs contains rings in the dust continuum. These rings are the
locations where pebbles accumulate, which is beneficial for planetesimal
formation and subsequent planet assembly. We investigate the viability of
planet formation inside ALMA rings in which pebbles are trapped by either a
Gaussian-shape pressure bump or by the strong dust backreaction. Planetesimals
form at the midplane of the ring via streaming instability. By conducting
N-body simulations, we study the growth of these planetesimals by collisional
mergers and pebble accretion. Thanks to the high concentration of pebbles in
the ring, the growth of planetesimals by pebble accretion becomes efficient as
soon as they are born. We find that type-I migration plays a decisive role in
the evolution of rings and planets. For discs where planets can migrate inward
from the ring, a steady state is reached where the ring spawns
${\sim}20M_{\oplus}$ planetary cores as long as rings are fed with materials
from the outer disc. The ring acts as a long-lived planet factory and it can
explain the “fine-tuned” optical depths of the observed dust rings in the
DSHARP large program. In contrast, in the absence of a planet removal
mechanism (migration), a single massive planet will form and destroy the ring.
A wide and massive planetesimals belt will be left at the location of the
planet-forming ring. Planet formation in rings may explain the mature
planetary systems observed inside debris discs.

###### keywords:

protoplanetary discs – circumstellar matter – planets and satellites:
formation – submillimetre: planetary systems

††pubyear: 2022††pagerange: Efficient planet formation by pebble accretion in
ALMA rings–B

##  1 Introduction

Protoplanetary discs are the cradle of planets. In the past two decades, more
than 5000 exoplanets have been found, demonstrating a variety in
configurations of exoplanet systems. At the same time, thanks to the progress
in its high sensitivity and spatial resolution, the Atacama Large
Millimeter/submillimeter Array (ALMA) reveals unprecedented images of planet-
forming discs in both continuum and line-emission (e.g., Andrews et al., 2018;
Long et al., 2018; Cieza et al., 2019; Öberg et al., 2021). The ALMA continuum
observation of protoplanetary discs allows the measurement of the mass budget
of (sub)millimeter size particles, a.k.a pebbles, which are the building
blocks for planets.

One striking feature of protoplanetary discs is the ubiquitous rings and gaps,
which challenges the classical framework where discs are considered as
‘smooth’ (e.g., Weidenschilling, 1977a; Hayashi, 1981). A large fraction of
the protoplanetary discs observed with ALMA display near-perfect annular rings
in the continuum, in which large quantities of dust appear to be concentrated
(e.g., van der Marel et al., 2019). The formation of these rings is an open
question. The most popular hypothesis is that rings are pebbles trapped by a
local pressure maxima (e.g., Dullemond et al., 2018). Due to aerodynamical
drag, the radial velocity of pebbles in protoplanetary discs follows the gas
pressure gradient. A local pressure bump is, therefore, a site for dust
trapping. A number of mechanisms can form such pressure bumps in
protoplanetary discs, including dead zone boundaries (e.g., Flock et al.,
2015), zonal flows (e.g., Bai & Stone, 2014), photoevaporation (e.g., Owen &
Kollmeier, 2019), snowlines (e.g., Kretke & Lin, 2007; Okuzumi et al., 2016),
and the most widely proposed mechanism — planets (e.g., Dong et al., 2015;
Zhang et al., 2018).

The planet-disc interaction provides a plausible way of explaining the disc
structure. In particular, the detection of localized kinematic deviations from
Keplerian motion in 12CO maps in a few discs makes the planet explanation
compelling (e.g., HD 163296, Teague et al., 2018; Pinte et al., 2020). Yet,
the confirmation of direct evidence of planets inside discs is challenging.
One promising way is by the H line emission, which originates from an
accretion front at the planetary surface. Dozens of discs have been targeted
by direct imaging facilities, but PDS 70 remains the only system where
accreting exoplanets were unambiguously identified (e.g., Keppler et al.,
2018; Haffert et al., 2019). In addition, the MAPS Large program released tens
of line-emission maps of five protoplanetary discs, allowing the assessment of
the correlation between dust substructures and gas substructures (Law et al.,
2021). Yet, statistical tests on the correlation between dust continuum and
gas substructures report no significant correlation (Jiang et al., 2022).
Therefore, it is debatable whether planets, or more generally, pressure bumps,
are the primary cause for rings in discs, which questions the frequently made
link between giant planets and the rings and gaps.

In addition, forming a planet at these great distances by core accretion
requires time. Even the acclaimed pebble accretion for which the accretion
rate is amplified by gas drag (Ormel & Klahr, 2010; Lambrechts & Johansen,
2012), it takes the disc lifetime for a planetary core to grow at wide orbits
(${>}10$ au) in smooth disc (Ormel, 2017). On the other hand, from the
observational aspect, an increasing number of rings are found in young class
0/I disc (e.g., Sheehan & Eisner, 2017, 2018; de Valon et al., 2020; Segura-
Cox et al., 2020; Sheehan et al., 2020). The age of these discs are less than
$0.5$ Myr, further exacerbating the timescale problem.

Alternatively, giant planets on distant orbits could also form by
gravitational instability (GI) in the disc (e.g., Boss, 1997). Yet, a study on
the mass function of direct-imaged planets suggests GI might not be the
primary source of giant planet formation (Wagner et al., 2019). Moreover, for
the potential planets that are assumed to be associated with the ALMA rings
and gaps, the inferred mass range covers a large parameter space and strongly
depends on the unknown level of disc turbulence. In particular, some potential
"planets" are inferred to be only several earth mass (Wang et al., 2021b),
while the planet formed by GI is usually a gas giant. Recently, Deng et al.
(2021) argues that this problem can be solved if the disc is magnetized, and
the planet could be formed via the magnetically controlled disc fragmentation.
However, the model relies on the strong magnetic field which is largely
unknown in the disc. These possible intermediate-mass planets are overall hard
to explain by gravitational instability.

While the origin of the rings is not well understood, due to the local
enhancement of the dust-to-gas ratio, the rings are believed to be perfect
sites of planetesimal formation by streaming instability and self-gravity
collapse of pebbles (e.g., Carrera et al., 2015; Yang et al., 2017; Li &
Youdin, 2021; Carrera et al., 2021; Klahr & Schreiber, 2021; Carrera et al.,
2022; Xu & Bai, 2022). Tens of Earth mass in solids are accumulated within
these rings, which may allow planetesimal formation up to the same magnitude
(e.g., Stammler et al., 2019; Jiang & Ormel, 2021). Whether or not these
planetesimals can grow into planets is an intriguing question, similar to the
age-old chicken or egg dilemma—are rings in discs caused by planets that were
already formed, or are rings intrinsic properties of discs where planets form?

There are two interpretations of how pebble rings will promote the formation
and growth of planets. One group of studies focuses on the evolution of a
swarm of planetesimals formed from a pebble ring in smooth discs. By assuming
planetesimals formed from streaming instability, Liu et al. (2019) study the
growth of a planetesimal belt formed from a ring at the ice line. A super-
Earth planet can quickly form by a hybrid approach of planetesimal accretion
and pebble accretion. Jang et al. (2022) extend the study by placing the ring
at a number of different radii. However, because of the lower pebble surface
density, higher aspect ratio, and longer orbital timescale at the outer disc,
the formed planetesimals can hardly grow at wide orbital distances
(${\gtrsim}30$ au), which is the location of the ALMA rings. Yet, in both
these works, the aerodynamics of pebbles is treated in a smooth pressure
gradient where pebbles always efficiently drift inward. The enhancement of the
pebble density at the ring location is not taken into account.

Alternatively, another group of studies focuses on the pebble aerodynamics
inside a pressure bump. As a pilot study on this question, Morbidelli (2020)
analytically solves the growth of a single Mars-mass core inside the B77 ring
in Elias 24 by assuming the ring is supported by a permanent pressure bump.
They find the embryos grow very slowly at 75 au (which can only grow to
$1M_{\oplus}$ after $\sim 3$ Myr). On the other hand, the growth of the embryo
turns very efficient when putting the same ring at 5 au, where the embryos can
grow to $10M_{\oplus}$ within $0.5$ Myr. Even faster growth is reported in
Guilera et al. (2020), where a $13M_{\oplus}$ planet can form in only $0.02$
Myr at a bump located at the water snowline (3 au). Recently, Chambers (2021)
and Lee et al. (2022) revisit this question, and both of them find the growth
of embryos can become much more efficient when taking larger St and/or a more
massive ring. Yet, in these works, only one isolated embryo is considered, and
the output is sensitive to the initial mass of this embryo. However,
realistically, a large number of planetesimals will form as demonstrated in
Liu et al. (2019) and Jang et al. (2022). In addition, the mutual interaction
of planetesimals may influence the evolution significantly (e.g. Levison et
al., 2012).

In this paper, we study the planetesimal formation and the following-up growth
of the planetesimal via pebble accretion inside a pebble ring consistently.
Two different assumptions on the ring nature are considered– the clumpy ring
and pressure bump supporting ring. The clumpy ring model assumes that the ring
is a manifestation of dense, clumpy mid-plane pebbles in a smooth gas disc
(Jiang & Ormel, 2021, hereafter JO21). In this model, due to the strong back-
reaction, the clumpy medium itself hardly experiences radial drift. The
accumulation of pebbles results in positive feedback to slow down the dust
drift, and eventually leads to a traffic jam shown as a pebble ring.

In the clumpy ring model, because of the high dust-to-gas ratios, pebbles get
into the clumpy state within the ring and their drift is significantly
decelerated. A similar condition is met in the ring induced by a pressure
bump. Therefore, the pebble accretion threshold will be much lower inside the
ring. Thus, in the interval where clump pebbles and planetesimals co-exist,
(small) planetesimals could efficiently feed from pebbles. On the other hand,
mutual scattering of planetesimals would oppose the growth (e.g., Schoonenberg
et al., 2019). In this work, we investigate the viability of planet formation
in ALMA rings by conducting N-body simulations that account for planetesimal
dynamics, pebble accretion, and planet migration. We also contrast the
standard pressure bump setup with the clumpy ring interpretation proposed by
JO21.

The plan of the paper is as follows. We describe the theoretical background
and numerical setup in Section 2. We present and explain the numerical results
in Section 3, which are further discussed in Section 4. We summarize the main
results and conclusions in Section 5.

##  2 Model

Figure 1:  Sketch illustrating the differences between the clumpy ring model
(left) and the ring induced by a pressure bump (right). The clumpy ring is
where pebbles experience strong pile-ups because the dust backreaction slows
down the radial drift of pebbles. The pebble scale height is determined by the
Kelvin-Helmholtz instability at the ring location in the CR model, while the
pebbles are still well coupled with gas in PB-induced rings. Pebbles can still
leak from the clumpy ring to the inner disc ($\dot{M}_{\rm leak}$). In
contrast, we assume that the pressure bump will fully filter the pebbles from
the inner disc.

We test two different types of planetesimal-forming rings in this work where
we group them by distinguishing the nature of the pebble concentration inside
the ring:

  1. 1.

Ring induced by a permanent pressure bump (PB)

  2. 2.

Ring sustained by aerodynamical effects only, in the context of the clumpy
ring (CR) model.

Figure 1 shows a cartoon that summarizes the differences between CR and PB
models. For the PB model, Morbidelli (2020) analytically solves for the planet
growth by pebble accretion for the B77 ring in Elias 24. The setup of the dust
ring is similar in our work with the one used by them. However, they directly
initialized a Mars mass embryo inside their framework and did not account for
the preceding planetesimals coagulation process. Instead, we take into account
the formation of planetesimal self-consistently and numerically solve the
problem in N-body simulations, which produces a more realistic outcome. The
planetesimal formation approach is based on the assumption that the dust
concentration may induce streaming instability, which can fulfill the criteria
for forming planetesimals. The details of planetesimal formation will be
discussed in Sect. 2.2.

The other scenario we studied is that of the clumpy ring model. In the CR
model, it is assumed that ALMA rings are the locations where pebbles have
collapsed as their number is too great for turbulence to sustain (Youdin &
Shu, 2002). Still, as they are not contained by pressure, pebbles "leak" at
the ring front edge, requiring a continuous supply of pebbles from the
exterior disc. An intriguing feature of the clumpy ring model in JO21 is the
"outward migration" of the ring, or more accurately, of the shape of the ring.
This is due to the fact that external mass flux coming from the outer disc and
the new traffic jam always happens at the tail of the queue and the head of
the queue will be leaked the first. Different from JO21 where the rings are
formed by solving the transport equations, to directly compare the PB and CR
models, we consider a scenario where the ring features no migration.

The most important difference between the CR and the PB models is the presence
of a local pressure maxima. The long-term survival of the ring in the CR model
is mainly due to the strong dust back-reaction on the radial drift and
diffusion of pebbles. Therefore, the CR model can also operate when the gas
pressure profile is smooth, and thus the morphology of the clumpy ring is
independent of the radial gas profile and usually has some asymmetric features
(i.e. steeper at the leading front edge and smoother at the following edge,
see discussion in Section 3.4 of JO21). For simplicity, We assume that the
rings are static and long-lived. Such an assumption applies when the mass flux
feeding the ring and the mass flux losing from the ring reach a perfect
balance. This scenario is also frequently found in JO21 when the radial
transport of clumps is moderate rather than totally halted, and therefore the
outward migration of the shape of the ring will be canceled with the advection
of clumpy pebbles.

###  2.1 Disc structure

In our study, we focus on the rings observed by ALMA (e.g., Andrews et al.,
2018; Long et al., 2018). We take the B74 ring in AS 209 as an example, whose
ring location is $r_{0}=74.2\,\rm au$, ring width is $w_{\rm ring}=3.38\,\rm
au$ and the disc mid-plane temperature is $T_{0}=16\,\rm K$ at the ring
location (Dullemond et al., 2018). For simplicity, we choose the stellar mass
$M_{\star}=1M_{\odot}$. The surface density and temperature of the gas in the
smooth disc is given by

$\Sigma_{\rm g}=\Sigma_{\rm g,0}\times\left(\frac{r}{r_{0}}\right)^{-1}$ (1)
$T=T_{0}\times\left(\frac{r}{r_{0}}\right)^{-0.5}$ (2)

where $\Sigma_{\rm g,0}=3.6\,\mathrm{g}\,\mathrm{cm}^{-2}$. We assume that the
disc is isothermal in the vertical direction. The isothermal sound speed is
$c_{s}=\sqrt{k_{B}T/\mu m_{\rm H}}$, with $m_{\rm H}$ the proton mass and
$\mu=2.3$ the mean molecular weight in atomic units. The disc is flared, and
the aspect ratio is

$h_{\rm
g}\equiv\frac{H_{\mathrm{g}}}{r}\equiv\frac{c_{s}}{\Omega_{K}r}=h_{0}\times\left(\frac{r}{r_{0}}\right)^{0.25}$
(3)

where $\Omega_{K}(r)=\sqrt{GM_{\star}/r^{3}}$ is the local Keplerian frequency
and $H_{\rm g}$ is the gas scale height. With substitutions, the aspect ratio
is $h_{0}=0.07$ at the characteristic location $r_{0}$. These parameters we
choose are consistent with the classical $\alpha$-accretion disc model
(Shakura & Sunyaev, 1973) where the gas mass accretion rate
$\dot{M}_{g}={3\pi\alpha c_{s}H_{\rm g}\Sigma_{\rm g}}\equiv
10^{-9}\times(\alpha/10^{-3})M_{\odot}\rm yr^{-1}$, where $\alpha$ is the
coefficient of gas viscosity.

####  2.1.1 Pressure bump induced ring (PB)

For the case where the pebble ring is formed by a prominent pressure bump, we
directly follow the calculation in Dullemond et al. (2018) and Morbidelli
(2020) to set up the pebble ring. The key assumptions is that the radial
pressure profile follow a nominal Gaussian shape around the ring:

$p(r)=p_{0}\exp\left(-\frac{(r-r_{0})^{2}}{2w_{\rm pb}^{2}}\right)$ (4)

where $r_{0}$ and $w_{\rm pb}$ are respectively the central location and width
of the pressure bump. Since the radial pressure support, the azimuthal gas
velocity $v_{\phi}$ is slightly different from the Keplerian velocity
$v_{K}=\Omega_{K}r$. The dimensionless measure of the radial pressure gradient
is (Nakagawa et al., 1986)

$\eta\equiv=\frac{v_{K}-v_{\phi}}{v_{K}}=-\frac{1}{2}\frac{c_{s}^{2}}{v_{K}^{2}}\frac{\partial\log{P}}{\partial\log{r}}.$
(5)

For the case with a pressure bump, by substituting Eq. (4), we obtain

$\eta_{\rm pb}=\frac{1}{2}\frac{r(r-r_{0})}{w_{\rm pb}^{2}}h_{\rm g}^{2}.$ (6)

Due to gas drag, the radial velocity of pebbles follows (Weidenschilling,
1977b)

$v_{\mathrm{dr}}=-\frac{2\mathrm{St}}{1+\mathrm{St}^{2}}\eta v_{K}$ (7)

where St is the pebble’s Stokes number. The pebble will concentrate towards
the pressure bump center. Yet, the pebble also experiences radial diffusion
with corresponding diffusivity (Youdin & Lithwick, 2007)

$D_{\mathrm{d,r}}=\frac{\delta_{\rm r}H_{\rm g}^{2}\Omega_{K}}{1+{\rm
St}^{2}}$ (8)

where we assume $\delta_{\rm r}=\alpha$ as a usual choice (e.g., Pinilla et
al., 2021). By balancing radial drift from gas drag with dust diffusion, the
surface density profile of the pebble ring reads

$\Sigma_{\rm r}(r)=\Sigma_{\rm r,0}\exp\left(-\frac{(r-r_{0})^{2}}{2w_{\rm
ring}^{2}}\right)$ (9)

where $\Sigma_{\rm r,0}$ is the peak surface density of the ring. The ring
width $w_{\rm ring}$ scales with the width of the pressure bump

$w_{\rm ring}=w_{\rm pb}\sqrt{\frac{\delta_{\rm r}}{\delta_{\rm r}+\rm St}}$
(10)

and the pebble mass of the ring can be calculated

$M_{\rm ring}=(2\pi)^{1.5}r_{0}w_{\rm ring}\Sigma_{\rm r,0}$ (11)

As argued by Morbidelli (2020), since the separation of the dust rings
observed in ALMA is typically much larger than the ring width, these formulas
should hold at least at the range of several $w_{\rm wring}$ away from the
ring center. In the background disc away from the ring, for simplicity, we do
not solve for the pebble density distribution, but use a constant pebble flux,
see Sect. 2.3 for detailed treatment. However, pebble accretion from the
pressure-supported background disc occurs at low efficiency at distant orbits.
Pebble accretion outside the ring region is significantly smaller in most of
our simulations compared with pebbles accreted from the ring.

Finally, for the PB model, we use the $\alpha$-prescription with dimensionless
diffusivity parameter $\delta_{\rm z}$ for vertical transport of the pebbles.
Assuming a vertical Gaussian distribution for the volume densities of both
pebble and gas, the scale height of the pebble can be obtained by balancing
the dust settling with vertical diffusion (Dubrulle et al., 1995; Birnstiel et
al., 2010)

$H_{\rm
p}=\sqrt{\frac{\delta_{\mathrm{z}}}{\delta_{\mathrm{z}}+\mathrm{St}}}H_{\mathrm{g}}$
(12)

For most of our simulations, we set the values of the gas viscosity $\alpha$
and pebbles’ diffusivity $\delta_{\rm r}$, $\delta_{\rm z}$ to be $10^{-3}$,
which is suggested by observation of dust distribution geometry in several
discs and widely used in theoretical studies (e.g., Pinte et al., 2016;
Dullemond et al., 2018; Villenave et al., 2020). Yet, values of these
parameters remain uncertain (see Pinilla et al., 2021; Miotello et al., 2022,
and references therein). For example, depending on the ring width and CO
kinematics of HD 163296 and AS 209, Rosotti et al. (2020) suggest $\delta_{\rm
r}=10^{-2}$ by assuming the balance between gas diffusion and dust decoupling
in the radial direction, while Villenave et al. (2022) finds $\delta_{\rm
z}=10^{-5}$ based on the thin dust layer of the edge-on disc in Oph 163131.
Moreover, it is also suggested that $\delta_{\rm r}$ and $\delta_{\rm z}$ can
be different in simulation with non-ideal MHD effects (e.g., Bai & Stone,
2014; Xu et al., 2017; Baehr & Zhu, 2021). Without loss of generality, we also
test the impact of lower disc turbulence ($\delta_{\rm z}=10^{-4}$) on our
model in parameter studies.

####  2.1.2 The clumpy ring (CR)

As explained at the start of this section, we assume a stationary ring for the
CR model in order to facilitate a direct comparison with the PB model. In
addition, we approximate the CR profile as Gaussian, Eq. (9), with the same
width $w_{\rm ring}=3.38\,$au.

In the clumpy ring model, these clumpy pebbles are concentrated in a thin
midplane of the disc. Following Sekiya (1998) and Youdin & Shu (2002), JO21
assumes that clump pebbles reside in the midplane with the scale height set by
Kelvin–Helmholtz instability

$H_{\rm c}=\sqrt{Ri_{\rm c}}\eta r\Psi(\psi)$ (13)

Hereafter referred to as the clump layer. For the critical Richardson number
$Ri_{\rm c}$, Chiang (2008) points out that when the vertically-integrated
dust-to-gas ratio is between 1 and 5 times solar abundance, the Richardson
number $Ri_{\rm c}\sim 0.1$. Thus we take $Ri_{\rm c}=0.1$ in this paper. The
dimensionless self-gravitational term
$\Psi(\psi)\equiv\sqrt{1+2\psi}-\psi\ln{[(1+\psi+\sqrt{1+2\psi})/\psi]}$ is
order of unity, where

$\psi\equiv\frac{4\pi
G\rho_{\mathrm{g}}}{\Omega_{K}^{2}}=0.16\times\left(\frac{\dot{M}_{g}}{10^{-9}M_{\odot}\,\rm
yr^{-1}}\right)\left(\frac{\alpha}{10^{-3}}\right)^{-1}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-1.5}\left(\frac{r}{r_{0}}\right)^{0.75}$
(14)

is the dimensionless self-gravity parameter of gas.

###  2.2 Planetesimal formation

We focus on the hypothesis that planetesimal are formed via streaming
instability in this work. Planetesimals formation is found to be natural
output of streaming instability in hydrodynamical simulations (Carrera et al.,
2015; Carrera et al., 2021; Yang et al., 2017; Li et al., 2019). A widely used
criterion for planetesimal formation in disc is the volume dust-to-gas ratio
exceeds unity (e.g., Johansen et al., 2006). Very recently, with improved
numerical resolution, Li & Youdin (2021) found this threshold can be slightly
lower for pebbles with $\rm St>0.01$, which make the planetesimal formation
more efficient. However, this does not qualitatively affect the picture of
planetesimal formation via SI.

We integrate the total mass located inside the planetesimal forming region
(where $\rho_{\mathrm{d}}>\rho_{\mathrm{g}}$ at the disc midplane), and the
planetesimal formation rate is (e.g., Drążkowska et al., 2016; Schoonenberg &
Ormel, 2017)

$\dot{M}_{\rm
plt}=\epsilon\frac{M(\rho_{\mathrm{d}}>\rho_{\mathrm{g}})}{t_{\rm sett}}$ (15)

where $\epsilon$ is the planetesimal formation efficiency and $t_{\rm
sett}=1/{\rm St}\Omega_{K}$ is the settling timescale.

In our N-body simulation, we assume that all planetesimals are formed with the
same initial mass $m_{\rm 0}$. For every $1000\,\rm yr$, we calculate the
planetesimal formation rate and the total number of planetesimal formed in the
past $1000\,\rm yr$. We add these planetesimals randomly inside the forming
region. We ensure that each newly formed planetesimal is separated from other
planetesimals by at least 5 Hill radii. We set the initial eccentricity
$e_{0}=10^{-5}$ and the initial inclination to be $i_{0}=e_{0}/2$. We tested
smaller values of $e_{0}$ and $i_{0}$ and find that since the orbits of
planetesimals are readily excited, the results are insensitive to the initial
values.

According to SI simulations (e.g., Schäfer et al., 2017), the mass of the
planetesimal formed by SI shows a top-heavy initial mass function with a
cutoff at the massive end. However, as demonstrated by Liu et al. (2019), the
most massive planetesimal would dominate the following mass growth and
dynamical evolution of the whole population. We choose the planetesimal masses
accordingly. By reviewing the literature on SI studies (Johansen et al., 2015;
Simon et al., 2016; Schäfer et al., 2017; Abod et al., 2019; Li et al., 2019),
Liu et al. (2020, see their Table 1) give a fit for the initial mass of the
planetesimal

$\displaystyle\frac{m_{\rm 0}}{M_{\oplus}}=$ $\displaystyle 5\times
10^{-3}\left(\frac{f_{\rm plt}}{40}\right)\left(\frac{C}{5\times
10^{-5}}\right)\left(\frac{\gamma}{\pi^{-1}}\right)^{a+1}$ (16)
$\displaystyle\left(\frac{h_{\rm
g}}{0.07}\right)^{3+b}\left(\frac{M_{\star}}{M_{\odot}}\right)$
$\displaystyle=$ $\displaystyle 2\times 10^{-3}\left(\frac{f_{\rm
plt}}{40}\right)\left(\frac{C}{5\times
10^{-5}}\right)\left(\frac{h_{0}}{0.07}\right)^{3+b}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-1.5a-0.5}$
$\displaystyle\left(\frac{\dot{M}_{g}}{10^{-9}M_{\odot}\,\rm
yr^{-1}}\right)^{a+1}\left(\frac{\alpha}{10^{-3}}\right)^{-a-1}\left(\frac{r}{r_{0}}\right)^{0.75a+0.75b+1.5}$

where the values of numerical factors ($a=0.5$, $b=0$, $C=5\times 10^{-5}$)
are chosen to be the same as in Liu et al. (2020, see their Section 2.4 for
fruitful discussion on choice of these parameters). We also control the
initial planetesimal mass through $f_{\rm plt}$, which is the ratio between
the maximum mass and the characteristic mass. The typical value falls between
tens to thousands in SI simulation (and $f_{\mathrm{plt}}=400$ in Liu et al.,
2020). We test initial masses ranging from Ceres to Mars, which corresponds to
$f_{\rm plt}=2\,\rm to\,2000$. For the fiducial value, we choose $f_{\rm
plt}=40$, which is consistent with a Pluto mass body. Surprisingly, we find
that the initial mass of planetesimals has little influence on the output in
our CR models, which we will discuss in Sect. 3.5.1.

###  2.3 Pebble accretion

Pebble accretion as a mechanism driving planet formation has seen active
development in the last decade (Johansen & Lambrechts, 2017; Ormel, 2017; Liu
& Ji, 2020; Drazkowska et al., 2022). However, in a smooth disc without
substructures, in order to trigger efficient pebble accretion at the outer
region of the disc ($>50$ au), a ${\sim}1000$ km size embryo is required
(Visser & Ormel, 2016). Thanks to the high concentration of pebbles inside the
ring, the threshold of PA is significantly lower.

We calculate the pebble accretion rate of each planetesimal depending on the
surface density of pebble and semimajor axis of

$\dot{M_{\rm PA}}={\cal R}\Sigma_{\rm peb}r^{2}\Omega_{K}$ (17)

where $\Sigma_{\rm peb}$ is the local pebble surface density. The accretion
efficiencies in the 2D and 3D limits read (converted from Liu & Ormel, 2018;
Ormel & Liu, 2018)

$\displaystyle{\cal R}$ $\displaystyle=({\cal R}_{\rm 2D}^{-2}+{\cal R}_{\rm
3D}^{-2})^{-0.5}$ (18a) $\displaystyle{\cal R}_{\rm 2D}$
$\displaystyle=4.0f_{\rm set}\sqrt{\frac{q{\rm St}\Delta v}{v_{K}}}$ (18b)
$\displaystyle{\cal R}_{\rm 3D}$ $\displaystyle=4.9f_{\rm set}^{2}\frac{q{\rm
St}}{h_{\rm peb}}$ (18c)

in which $q=m/M_{\star}$ is the mass ratio of planetesimal and $\Delta v$ (Eq.
(50)) is the approaching velocity of pebble relative to the planetesimal. The
settling factor $f_{\rm set}$ is always around unity. Detailed derivation and
formulas are in Appendix A. Since planetesimals with high inclination
$i>h_{\rm d}$ only interact with pebbles over a fraction $~{}i/h_{\rm d}$ of
their orbits. The pebble scale height coefficient should be modified as

$h_{\rm peb}=\sqrt{h_{\rm d}^{2}+\frac{\pi
i^{2}}{2}\left(1-\exp{\left[-\frac{i}{2h_{\rm d}}\right]}\right)}$ (19)

where $i$ is the inclination and

$h_{\rm d}\equiv\frac{H_{\rm d}}{r}=\left\\{\begin{array}[]{lr}\frac{H_{\rm
p}}{r}&\rm PB\\\ \frac{H_{\rm c}}{r}&\rm CR\end{array}\right.$ (20)

For the PB model, the planetesimal typically has $i\ll h_{\rm d}$. While for
the CR model, the punishment due to the inclination on the pebble accretion
rate could be significant because of the thin pebble layer.

In our simulations, we calculate the pebble accretion rate for each
planetesimal. Thus the total mass flux of $N$ planetesimals accreting from the
ring is

$\dot{M}_{\rm PA,ring}=-\sum_{i}^{N}\dot{M}_{\rm PA,i}$ (21)

For both PB and CR, we solve for the ring mass semi-analytically. We set an
external mass flux $\dot{M}_{\rm ext}$, which is the total mass flux of
pebbles drifting from the outer disc to feed the ring location. For PB runs,
we assume that the amplitude of the pressure bump is strong enough to halt all
of the pebbles–a perfect dust filtration (e.g., Zhu et al., 2012). However,
since there is no pressure support to stop the radial drift, in the CR discs,
the rings always leak pebbles downstream. Following JO21, the leaking mass
flux is

$\displaystyle\dot{M}_{\mathrm{leak}}$
$\displaystyle\simeq\mathrm{St}\Sigma_{\mathrm{g,0}}h_{0}^{3}M_{\star}^{0.5}r^{0.25}$
(22)
$\displaystyle=50\left(\frac{\mathrm{St}}{10^{-2}}\right)\left(\frac{\Sigma_{\mathrm{g,0}}}{3.6\,\mathrm{g}\,\mathrm{cm}^{-2}}\right)\left(\frac{h_{0}}{0.07}\right)^{3}$
$\displaystyle\times\left(\frac{M_{\star}}{M_{\odot}}\right)^{0.5}\left(\frac{r}{r_{0}}\right)^{0.25}\,M_{\oplus}\,\mathrm{Myr}^{-1}$

This means that a moderate amount of pebbles can leak from the ring toward the
host star, which might further promote the growth of the embryos or the ring
in the inner disc. For completeness, we account for the pebble accretion from
the smooth pebble background disc away from the ring. We introduce a
background factor

$f_{\rm bkg}=\frac{1}{2}(1+\tanh[\frac{2(|r_{0}-r|-2w_{\rm ring})}{w_{\rm
ring}}])$ (23)

which is zero within one $w_{\rm ring}$ width of the ring center, and will
increase to unity rapidly when $>3w_{\rm ring}$ away from the central location
of the ring. Figure 2 show the radial distribution of the $f_{\rm bkg}$
factor. The pebble flux upstream is $\dot{M}_{\rm ext}$ and the pebble flux
downstream is $\dot{M}_{\rm leak}$. For PB rings, planetesimal may only be
able to feed from the upstream flux. With the the given pebble flux, we simply
use the pebble accretion efficiency $\zeta$ to calculate the accretion rate
(see Eq. (42) and Eq. (43), Ormel & Liu, 2018). The result is then

$\dot{M}_{\rm PA,disc}\equiv\left\\{\begin{array}[]{lr}\dot{M}_{\rm
PA,up}=-f_{\rm bkg}\zeta\dot{M}_{\rm ext}&\rm r>r_{0}\\\ \dot{M}_{\rm
PA,down}=-f_{\rm bkg}\zeta\dot{M}_{\rm leak}&r<r_{0},\,\rm CR\\\ \dot{M}_{\rm
PA,down}=0&r<r_{0},\,\rm PB\end{array}\right.$ (24)

In most of the simulation, $\dot{M}_{\rm PA,disc}$ is generally very low
compared with $\dot{M}_{\rm PA,ring}$ in Eq. (21). Because of mass
conservation, the mass flux that eventually feeds the ring is then

$\dot{M}_{\rm ring}=\dot{M}_{\rm ext}-\dot{M}_{\rm
leak}-\sum_{i}^{N}\dot{M}_{\rm PA,up,i}.$ (25)

Meanwhile, the ring also loses mass by planetesimal formation Eq. (15) and
pebble accretion Eq. (21). We integrate the ring mass by calculating the mass
flux change in the N-body simulation. We then change the ring mass Eq. (11)
and the pebble surface density at the ring peak Eq. (9) accordingly.

Figure 2:  Radial distribution of the $f_{\rm bkg}$ factor (orange line). Red
orange line indicate the radial profile of the ring’s surface density . The
separation of two neighbored blue dash lines is $w_{\rm ring}$.

###  2.4 Planet migration

Massive bodies feel the gravitational torques from the disc gas (so called
type I migration, e.g. Goldreich & Tremaine, 1979; Kley & Nelson, 2012;
Baruteau et al., 2014). We implement the following expressions in the N-body
code as accelerations experienced by the planets due to type I migration
torque (e.g., Papaloizou & Larwood, 2000)

$\displaystyle\bm{a}_{\rm m}$ $\displaystyle=-\frac{\bm{v}}{t_{\rm mg}}$ (26)
$\displaystyle\bm{a}_{\rm e}$
$\displaystyle=-2\frac{(\bm{v}\cdot\bm{r})\bm{r}}{r^{2}t_{\rm e}}$
$\displaystyle\bm{a}_{\rm i}$ $\displaystyle=-\frac{\bm{v}_{\rm z}}{t_{\rm
i}}$

where $\bm{v}=(v_{x},v_{y},v_{z})$ and $\bm{r}=(x,y,z)$ are the velocity and
position of the planet. The characteristic migration timescale $t_{\rm mg}$
for a planet on a circular orbit is (e.g., Cresswell & Nelson, 2008)

$\displaystyle t_{\rm mg}$ $\displaystyle=\frac{1}{2}\frac{1}{f_{\rm
mg}}\frac{1}{q}\frac{M_{\star}}{\Sigma_{\rm g}a_{\rm p}^{2}}h_{\rm
g}^{2}\frac{1}{\Omega_{K}}$ (27) $\displaystyle=40\times\left(\frac{f_{\rm
mg}}{1}\right)^{-1}\left(\frac{m_{\rm
p}}{1\,M_{\oplus}}\right)^{-1}\left(\frac{\Sigma_{\rm g,0}}{3.6\,\rm
g\,cm^{-2}}\right)^{-1}\left(\frac{h_{0}}{0.07}\right)^{2}\left(\frac{a_{\rm
p}}{r_{0}}\right)\,\rm Myr$

where $a_{\rm p}$ is the semimajor axis of the planet and $f_{\rm mg}$ the
type I migration prefactor. We take $f_{\rm mg}=1$ in most of our simulations
but also test different values for parameter studies. Similarly,

$\displaystyle t_{\rm e}$ $\displaystyle=t_{\rm i}=\frac{1}{f_{\rm
mg}}\frac{1}{q}\frac{M_{\star}}{\Sigma_{\rm g}a_{\rm p}^{2}}h_{\rm
g}^{4}\frac{1}{\Omega_{K}}$ (28) $\displaystyle=0.4\times\left(\frac{f_{\rm
mg}}{1}\right)^{-1}\left(\frac{m_{\rm
p}}{1\,M_{\oplus}}\right)^{-1}\left(\frac{\Sigma_{\rm g,0}}{3.6\,\rm
g\,cm^{-2}}\right)^{-1}\left(\frac{h_{0}}{0.07}\right)^{4}\left(\frac{a_{\rm
p}}{r_{0}}\right)\,\rm Myr$

are the eccentricity, and inclination damping timescales. This prescription is
based on the assumption that the disc is isothermal and the density gradient
of the gas is smooth (Cresswell & Nelson, 2008). The direction of planet
migration is always inwards. For the PB model, if the pressure bump is caused
by a strong gas density enhancement, due to the variation of the gas density
gradient, a migration trap can happen at the ring location (e.g., Masset et
al., 2006; Morbidelli, 2020). Alternatively, the pressure trapping could also
be caused by the local temperature variation as suggested by some observation
(e.g., Kim et al., 2020; Rosotti et al., 2021). Rather than opting for a
specific treatment for the migration, we here still apply the same simple
prescription to both CR and PB models. For comparison, we run the simulation
without migration $f_{\rm mg}\to 0$, which mimics the case where the migration
trap happens and the planet stays inside the ring, see Sect. 3.4.

###  2.5 Other numerical setups

####  2.5.1 Initial setup of the ring

We initialize the ring mass such that the mid-plane dust-to-gas ratio
$Z\equiv\rho_{\rm d}/\rho_{\rm g}=1$ at the ring peak

$\Sigma_{\rm r,0}=\Sigma_{\rm g,0}\frac{H_{\rm d}}{H_{\rm g}}$ (29)

The motivation of this setup is simply because that there is no planetesimal
formation when the dust-to-gas ratio is below unity, and all of the mass flux
can be fed to the ring.We estimate the timescale required for the ring to grow
to the initial mass by $t_{\rm load}=M_{\rm ring,0}/\dot{M}_{\rm ring}$, whose
value is listed in Table 3 for each run.

####  2.5.2 Planetesimal removal

The excitation of the eccentricity and inclination of planetesimals formed in
the ring is mainly due to planetesimals’ viscous stirring by themselves. As we
derived, highly excited planetesimal will have a low pebble accretion rate
because of their high relative velocity towards pebble and the low pebble
density out of the ring region. These small planetesimals contribute only a
small fraction of the total mass, and grow little in long time evolution but
dominate in the number of planetesimals, which significantly slows down the
N-body simulations.

In order to improve the computational efficiency, we remove planetesimals that
are scattered away from the ring and therefore cannot accrete efficiently. We
introduce a space factor $f_{\rm space}$, which is the area that the
particle’s orbit overlap with the ring region in $r-z$ plane

$\displaystyle f_{\rm space}$ $\displaystyle=\dfrac{\min((1+e)a_{\rm
p},r_{0}+w_{\rm ring})-\max((1-e)a_{\rm p},r_{0}-w_{\rm ring})}{2ea_{\rm p}}$
(30) $\displaystyle\times\min(1,\frac{h_{\rm d}}{ia_{\rm p}})$

For planetesimals whose orbits fall in the cube centered at the center the
ring with length $2w_{\rm ring}$ and height $h_{\rm d}$, the space factor
$f_{\rm space}=1$. We remove the particle only when its $f_{\rm space}$ falls
below 0.1 and the mass is lower than Mars. Planetesimals that are currently
inside the ring region will not be removed. In addition, we rank the
planetesimals by mass and never remove those whose mass ranks in the top
$N_{\rm min}=100$. We also test the case $N_{\rm min}=200,300,500,1000$ for
our default model and we find that our results are insensitive to the choice
of $N_{\rm min}$.

We set an inner boundary $r_{\rm in}=r_{0}-5w_{\rm ring}$ for our simulation
with type-I migration on, i.e., planets who migrate and pass through the inner
boundary will be removed. We discuss the post-ring evolution of these planets
in Sect. 3.6.

The simulations were conducted using the IAS15 integrator of the open-source
N-body code REBOUND, with collisions between bodies treated as perfect mergers
that conserve the linear momentum (Rein & Liu, 2012). The pebble accretion
rate expressions and type-I migration force were implemented via REBOUNDx
(Tamayo et al., 2020).

##  3 Results

The aim of this work is to test the viability of planet formation in
planetesimal-forming rings. Without loss of generality, we pick the default
parameters similar to the values used in JO21 and choose the ring location and
shape the same as the B77 ring in AS 209. We summarize these default
parameters in Table 1 in bold. Different setups between the CR and PB runs are
also listed. In Table 2, we give all runs and the corresponding key parameters
studied in each run.

Table 1: Summary of parameters studied and comparison between the Clumpy Ring (CR) and the Pressure-Bump-induced ring (PB).  name | Clumpy Ring (CR) | Pressure Bump (PB)
---|---|---
Same
$r_{0}$ |  30, 74.2, 154 au
$w_{d}$ | 3.38 au
$m_{\rm 0}$ |  $M_{\rm Ceres}$, $M_{\rm Pluto}$, $M_{\rm Moon}$, $M_{\rm Mars}$
$\dot{M}_{\rm ext}$ |  60, 80, 100, 150, 300 $M_{\oplus}\,\rm Myr^{-1}$
$\delta_{\rm z}$ |  1e-4, 1e-3
St |  1e-3, 3e-3, 1e-2, 3e-3, 1e-1
$\epsilon$ |  1e-5, 1e-4, 1e-3, 1e-2
Difference
$w_{\rm pb}$ | - | $\sqrt{\frac{St}{\alpha}}w_{d}$
$H_{\rm d}$ | $\sim\eta r$ | $\sqrt{\frac{\alpha}{\rm St}}H_{\rm g}$
$\dot{M}_{\rm leak}$ |  Eq. (22) | 0
$\eta$ |  Eq. (5) |  Eq. (6)
Table 2: Model runs input parameters. (1) Stokes number; (2) external mass flux; (3) leaking mass flux Eq. (22); (4) note on the focus of parameter studies run-ida | St | $\dot{M}_{\rm ext}$ | $\dot{M}_{\rm leak}$ | note
---|---|---|---|---
|  |  [$M_{\oplus}\rm\,Myr^{-1}$] |
| (1) | (2) | (3) | (4)
runs varying initial planetesimal mass | $m_{\rm 0}$
cr-mMars | 0.01 | 100 | 50 | $M_{\rm Mars}$
cr-mMoon | 0.01 | 100 | 50 | $M_{\rm Moon}$
cr-default | 0.01 | 100 | 50 | $M_{\rm Pluto}$
cr-m05Pluto | 0.01 | 100 | 50 | $0.5\,M_{\rm Pluto}$
cr-m03Pluto | 0.01 | 100 | 50 | $0.3\,M_{\rm Pluto}$
cr-m02Pluto | 0.01 | 100 | 50 | $0.2\,M_{\rm Pluto}$
cr-mCeres | 0.01 | 100 | 50 | $M_{\rm Ceres}$
pb-mMars | 0.01 | 100 | 0 | $M_{\rm Mars}$
pb-mMoon | 0.01 | 100 | 0 | $M_{\rm Moon}$
pb-default | 0.01 | 100 | 0 | $M_{\rm Pluto}$
cr-F60 | 0.01 | 60 | 50 |
cr-F80 | 0.01 | 80 | 50 |
cr-F150 | 0.01 | 150 | 50 |
cr-F300 | 0.01 | 300 | 50 |
cr-St001 | 0.001 | 55 | 5 |
cr-St003 | 0.003 | 65 | 15 |
cr-St030 | 0.03 | 200 | 150 |
cr-St100 | 0.1 | 550 | 500 |
pb-St030 | 0.03 | 100 | 0 |
pb-St100 | 0.1 | 100 | 0 |
cr-St001F100 | 0.001 | 100 | 5 |
cr-St003F100 | 0.003 | 100 | 15 |
cr-St006F100 | 0.006 | 100 | 30 |
cr-St015F100 | 0.015 | 100 | 75 |
cr-St018F100 | 0.018 | 100 | 90 |
runs varying planetesimal formation efficiency | $\epsilon$
cr-e2 | 0.01 | 100 | 50 | $10^{-2}$
cr-default | 0.01 | 100 | 50 | $10^{-3}$
cr-e4 | 0.01 | 100 | 50 | $10^{-4}$
cr-e5 | 0.01 | 100 | 50 | $10^{-5}$
runs varying vertical dust diffusivity | $\delta_{\rm z}$
cr-default | 0.01 | 100 | 50 | $10^{-3}$
cr-a4 | 0.01 | 100 | 50 | $10^{-4}$
pb-default | 0.01 | 100 | 0 | $10^{-3}$
pb-a4 | 0.01 | 100 | 0 | $10^{-4}$
runs varying vertical ring location | $r_{0}$
cr-r30 | 0.01 | 100 | 40 | 30 au
cr-default | 0.01 | 100 | 50 | 74.2 au
cr-r154 | 0.01 | 100 | 60 | 154 au
cr-05mg | 0.01 | 100 | 50 | $f_{\rm mg}=0.5$
cr-20mg | 0.01 | 100 | 50 | $f_{\rm mg}=2$
cr-nmg | 0.01 | 100 | 50 | $f_{\rm mg}\to 0$
pb-nmg | 0.01 | 100 | 0 | $f_{\rm mg}\to 0$
pb-St030-nmg | 0.03 | 100 | 0 | $f_{\rm mg}\to 0$
pb-St100-nmg | 0.1 | 100 | 0 | $f_{\rm mg}\to 0$
Table 3: Model runs results. (1) time required to load the ring mass to $\rho_{\rm d}/\rho_{\rm g}=1$ by the external mass flux; (2) time when the first protoplanet leaves the inner simulation domain; (3) average time interval between protoplanets leaving the box; error bars indicate the maximum and minimum intervals duration; (4) time when the sixth protoplanet passes the inner simulation domain; (5) the mass of the biggest protoplanet; (6) the average mass of protoplanet passing the inner simulation domain; (7) the predicted representative planet mass Eq. (32); (8) the initial mass of the ring; (9) the measured averaged ring mass in the saturation stage; (10) the predicted saturated mass of the ring Eq. (34); (11) section reference where the model is mainly discussed. run-ida | $t_{\rm load}$ | $t_{\rm 1}$ | $\langle\Delta\,t\rangle$ | $t_{\rm 6}$ | $m_{\rm max}$ | $m_{\rm avg}$ | $m_{\rm rp}$ | $m_{\rm ring,0}$ | $m_{\rm ring,ms}$ | $m_{\rm ring,ss}$ | Section
---|---|---|---|---|---|---|---|---|---|---|---
| [Myr] |  [$M_{\oplus}$] |
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
cr-mCeres | 0.23 | 1.06 | $0.27^{+0.28}_{-0.26}$ | 2.44 | 15.1 | 11.8 | 18 | 11.6 | 34.8 | 34 | 3.5.1
cr-m02Pluto | 0.23 | 1.04 | $0.34^{+0.43}_{-0.27}$ | 2.73 | 21.7 | 12.1 | 18 | 11.6 | 33.1 | 34 | 3.5.1
cr-m03Pluto | 0.23 | 0.97 | $0.17^{+0.31}_{-0.16}$ | 1.83 | 21.5 | 10.6 | 18 | 11.6 | 31.6 | 34 | 3.5.1
cr-m05Pluto | 0.23 | 1.01 | $0.14^{+0.25}_{-0.14}$ | 1.73 | 13.5 | 8.3 | 18 | 11.6 | 30.7 | 34 | 3.5.1
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.1
cr-mMoon | 0.23 | 0.85 | $0.35^{+0.29}_{-0.26}$ | 2.61 | 19.2 | 13.4 | 18 | 11.6 | 32.6 | 34 | 3.5.1
cr-mMars | 0.23 | 0.79 | $0.28^{+0.33}_{-0.26}$ | 2.18 | 17.9 | 10.7 | 18 | 11.6 | 28.2 | 34 | 3.5.1
pb-default | 1.61 | 0.65 | $0.04^{+0.07}_{-0.04}$ | 0.87 | 19.7 | 6.9 | 25 | 160.5 | 141.1 | 54 | 3.3
pb-mMoon | 1.61 | 0.61 | $0.04^{+0.06}_{-0.03}$ | 0.80 | 9.8 | 5.3 | 25 | 160.5 | 152.4 | 54 | 3.5.1
pb-mMars | 1.61 | 0.49 | $0.20^{+0.18}_{-0.15}$ | 1.51 | 33.6 | 17.3 | 25 | 160.5 | 167.4 | 54 | 3.5.1
cr-F60 | 1.16 | 1.71 | $0.56^{+0.84}_{-0.53}$ | 4.50 | 9.0 | 5.6 | 8 | 11.6 | 13.9 | 12 | 3.5.2
cr-F80 | 0.39 | 1.24 | $0.47^{+0.59}_{-0.36}$ | 3.60 | 21.5 | 11.3 | 14 | 11.6 | 23.9 | 24 | 3.5.2
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.5.2
cr-F150 | 0.12 | 0.77 | $0.22^{+0.26}_{-0.21}$ | 1.85 | 26.1 | 17.5 | 26 | 11.6 | 46.2 | 53 | 3.5.2
cr-F300 | 0.05 | 0.58 | $0.04^{+0.10}_{-0.04}$ | 0.79 | 29.1 | 9.7 | 41 | 11.6 | 92.4 | 99 | 3.5.2
cr-St001 | 0.23 | 2.31 | $0.56^{+0.69}_{-0.43}$ | 5.12 | 11.9 | 7.9 | 18 | 11.6 | 126.1 | 156 | 3.5.2
cr-St003 | 0.23 | 1.40 | $0.38^{+0.28}_{-0.30}$ | 3.32 | 19.7 | 8.0 | 18 | 11.6 | 82.5 | 75 | 3.5.2
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.5.2
cr-St030 | 0.23 | 0.89 | $0.22^{+0.40}_{-0.19}$ | 1.99 | 28.1 | 13.2 | 18 | 11.6 | 12.6 | 16 | 3.5.2
cr-St100 | 0.23 | 0.69 | $0.12^{+0.12}_{-0.10}$ | 1.29 | 23.4 | 13.1 | 18 | 11.6 | 5.0 | 7 | 3.5.2
pb-St030 | 0.96 | 0.29 | $0.06^{+0.05}_{-0.04}$ | 0.57 | 14.9 | 8.6 | 25 | 95.6 | 39.0 | 26 | 3.3
pb-St100 | 0.53 | 0.23 | $0.13^{+0.30}_{-0.12}$ | 0.87 | 30.1 | 14.0 | 25 | 53.0 | 16.4 | 12 | 3.3
cr-St001F100 | 0.12 | 1.71 | $0.43^{+0.37}_{-0.26}$ | 3.88 | 24.2 | 10.4 | 25 | 11.6 | 216.3 | 240 | 3.5.2
cr-St003F100 | 0.14 | 1.28 | $0.25^{+0.28}_{-0.14}$ | 2.52 | 27.1 | 12.8 | 24 | 11.6 | 109.3 | 107 | 3.5.2
cr-St006F100 | 0.17 | 1.00 | $0.33^{+0.45}_{-0.26}$ | 2.65 | 32.3 | 15.0 | 22 | 11.6 | 60.2 | 59 | 3.5.2
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.5.2
cr-St015F100 | 0.46 | 1.02 | $0.40^{+0.24}_{-0.29}$ | 3.01 | 13.2 | 11.0 | 13 | 11.6 | 14.5 | 16 | 3.5.2
cr-St018F100 | 2.32 | 1.15 | $0.51^{+0.95}_{-0.51}$ | 3.69 | 10.2 | 7.0 | 8 | 11.6 | 8.8 | 8 | 3.5.2
cr-e2 | 0.23 | 0.93 | $0.39^{+0.27}_{-0.39}$ | 2.90 | 13.2 | 9.6 | 18 | 11.6 | 22.1 | 34 | 3.5.3
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.5.3
cr-e4 | 0.23 | 0.91 | $0.17^{+0.31}_{-0.16}$ | 1.74 | 19.6 | 9.5 | 18 | 11.6 | 32.0 | 34 | 3.5.3
cr-e5 | 0.23 | 0.97 | $0.16^{+0.35}_{-0.16}$ | 1.78 | 22.0 | 12.0 | 18 | 11.6 | 35.6 | 34 | 3.5.3
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.1
cr-a4 | 0.23 | 0.74 | $0.09^{+0.07}_{-0.09}$ | 1.19 | 50.7 | 23.5 | 18 | 11.6 | 23.7 | 34 | 3.6
pb-default | 1.61 | 0.65 | $0.04^{+0.07}_{-0.04}$ | 0.87 | 19.7 | 6.9 | 25 | 160.5 | 141.1 | 54 | 3.3
pb-a4 | 1.61 | 0.63 | $0.14^{+0.33}_{-0.13}$ | 1.32 | 59.3 | 28.5 | 25 | 53.0 | 59.0 | 54 | 3.3
cr-r30 | 0.17 | 0.90 | $0.16^{+0.34}_{-0.14}$ | 1.71 | 17.2 | 10.0 | 20 | 10.4 | 15.5 | 24 | 3.5.4
cr-default | 0.23 | 0.97 | $0.18^{+0.33}_{-0.17}$ | 1.87 | 22.2 | 10.6 | 18 | 11.6 | 30.5 | 34 | 3.1
cr-r154 | 0.30 | 1.23 | $0.26^{+0.34}_{-0.16}$ | 2.51 | 22.8 | 8.3 | 16 | 12.3 | 43.2 | 42 | 3.5.4
cr-05mg | 0.23 | 1.20 | $0.40^{+0.30}_{-0.27}$ | 3.21 | 26.9 | 14.0 | 25 | 11.6 | 24.7 | 27 | 3.4
cr-20mg | 0.23 | 0.78 | $0.22^{+0.16}_{-0.16}$ | 1.89 | 13.0 | 8.1 | 13 | 11.6 | 40.5 | 43 | 3.4
cr-nmg | 0.23 |  1.57a |  |  |  60b |  |  | 11.6 |  11.0c |  | 3.4
pb-nmg | 1.61 | 1.60 |  |  | 60 |  |  | 160.5 | 155.1 |  | 3.4
pb-St030-nmg | 0.96 | 0.74 |  |  | 60 |  |  | 95.6 | 51.3 |  | 3.4
pb-St100-nmg | 0.53 | 0.84 |  |  | 60 |  |  | 53.0 | 16.6 |  | 3.4
a since there is no migration, $t_{\rm 1}$ is defined as the time when the
most massive planet open a gap, see Sect. 3.4
b gap-opening mass Eq. (38) |  c the mass of the ring at $t_{\rm 1}$, the ring should eaten by the planet eventually

###  3.1 Default model

Figure 3 shows the trajectories of protoplanets that grow to a mass higher
than $1.7\,M_{\oplus}$ during the simulation. We run the simulation until
$3\,\mathrm{Myr}$, which is long enough to spawn ten protoplanets as shown in
the Fig. 3. The criterion of $1.7\,M_{\oplus}$ is chosen so that the planets
are massive enough to migrate inward sufficiently. According to Eq. (27), a
planet as massive as $1.7\,M_{\oplus}$ can migrate over a distance equal to
the ring width in $w_{\rm ring}/r_{0}\times t_{\rm mg}=1\,\mathrm{Myr}$.
Finally, nine planets have passed the inner boundary of the simulation box. We
pick the first six as the sample for further statistic study. In order to
compare the default model with other runs equally, we continue each simulation
until six planets have crossed the inner domain boundary. This number is
sufficient to obtain statistically viable results while limiting the
computational cost.

Figure 3:  Protoplanet trajectories (solid lines) of the default model. The
colour of the solid lines corresponds to the inclination of the protoplanets.
The line width is linear with the mass of the planet. The cyan background
depicts the background planetesimal densities (other planetesimals except for
these protoplanets). The two gray vertical line marks the $r_{0}\pm w_{\rm
ring}$ width of the ring. The horizontal lines separate the three phases of
the evolution, see main texts.

Different from pebble accretion in smooth discs, thanks to the highly
concentrated dust inside the ring, the efficiency of pebble accretion is
significantly enhanced. The enhancement is due to both the higher pebble
density at the disc midplane and the lower headwind around the ring peak. The
former condition is naturally met as the planetesimals can form only in the
dense midplane with a high mass concentration inside the ring. Therefore,
pebble accretion becomes efficient as soon as the planetesimal is formed! As
Fig. 4 shows, the growth timescale of the first planet is less than $0.1$ Myr
once it was born. The blue and yellow lines in Fig. 4 show the nominal pebble
accretion timescale

$\displaystyle t_{\rm PA}=\frac{m}{\dot{M}_{\rm PA}}$
$\displaystyle=\left\\{\begin{array}[]{lr}0.35\frac{m^{1/3}M_{\star}^{2/3}}{{\rm
St}^{2/3}\Sigma_{\rm peb}r^{2}\Omega_{K}}&\rm 2D\,regime\\\
0.20\frac{M_{\star}h_{\rm peb}}{{\rm St}\Sigma_{\rm peb}r^{2}\Omega_{K}}&\rm
3D\,regime\end{array}\right.$ (31)

where we take $\Sigma_{\rm peb}$ the surface density at the ring peak and
assume the relative velocity $\Delta v$ is determined by the Keplerian shear
Eq. (45) and $f_{\rm set}=1$.

Figure 4:  Pebble accretion and migration timescales as functions of
protoplanet mass. The dark blue line represents the measured growth timescale
$m/\dot{m}$ of the first planet that reaches the inner domain boundary in the
cr-tdefault run. The light brown line indicates the type-I migration timescale
to pass the ring $w_{\rm ring}/r_{0}\times t_{\rm mg}$. The light blue,
orange, and green lines indicate the analytical calculation of 2D, 3D, and
combined pebble accretion timescale $m/\dot{M}_{\rm PA}$, assuming a surface
density equal to that of at the ring peak (Eq. (17) and Eq. (18)). The ring
mass is increasing while the first planet is growing, which initially
decreases the pebble accretion growth timescale. The deep dip in the growth
timescale around $1M_{\oplus}$ records a collisional merger. The vertical
black line marks the transition from phase I to phase II.

The evolution of protoplanets can be divided into three phases depending on
the evolution of the pebble accretion rate and the ring mass (see Fig. 5). In
the first phase, the number of the planetesimal is still small. Pebble
accretion has negligible influence on the evolution of the ring. Almost all
incoming pebbles end up in the ring, thus the ring mass linearly grows with
time at a slope that is the same as $\dot{M}_{\rm
net}=\dot{M}_{\mathrm{ext}}-\dot{M}_{\mathrm{leak}}$. Meanwhile, planetesimals
continue forming. Thus, the total pebble accretion rate monotonically
increases with time because of both the elevated number of planetesimals and
the growth of each planetesimal. Yet, as shown in the middle panel of Fig. 5,
the fraction of pebble accretion contributed by the largest particle keeps
decreasing in phase I. The accretion rate of each planetesimal differs little
since the pebble accretion is not a runaway mechanism (Ormel, 2017,
$\dot{M}_{\rm PA}\propto m$ in 3D regime, see also Eq. (18)).

Around 0.2 Myr, the total pebble accretion rate becomes non-negligible
compared with the mass flux fed to the ring, which slows down the growth of
the ring. As demonstrated in the bottom panel of Fig. 5, the slope of the ring
mass flattens. At the same time, the eccentricity and inclination of
planetesimals are excited by viscous stirring. Smaller planetesimals are more
easily stirred up. Since the small scale-height of the pebble in the clumpy
ring model, the high inclination will lower the efficiency of pebble accretion
dramatically, suppressing further growth of these planetesimals by pebble
accretion. Conversely, bigger planetesimals tend to damp their inclination by
planetesimal dynamical friction and stay closer to the midplane region where
the pebble density is the highest. As shown in Fig. 4, the first and also the
biggest protoplanet transitions from the 3D accretion regime to the 2D
accretion regime, in which the pebble accretion cross section is higher than
the pebble scale height. The measured growth timescale matches the analytical
formula very well because the planet is staying close to the ring center and
keeps a low inclination. As a result, the biggest protoplanet grows the
fastest, while the small bodies are effectively left behind (Levison et al.,
2012). As shown in Fig. 5, at the top panel, the total pebble accretion rate
does not increase significantly. The green region, which represents the
accretion rate of the largest embryos, keeps rising. The fraction of the
accretion rate of the largest embryos over the total accretion rate is plotted
in the middle panel of Fig. 5. At the beginning of phase II, it is less than
0.1 and becomes dominant as time evolves.

In Fig. 3, we show the evolution trajectories of growing protoplanets. The
colour of the trajectory line indicates the inclination of the embryos. The
first planet keeps a low inclination, while the inclination and eccentricity
of other planetesimals are excited by dynamical friction. The cyan background
corresponds to the density of planetesimals, in which the protoplanet
(trajectory lines) are not included. The semimajor axes of planetesimals are
expanding. Overall, in phase II, the biggest protoplanet stays closer to the
mid-plane region and stirs up the newly formed planetesimals. The ring mass
still grows but the growth rate becomes lower, since this protoplanet accretes
from the ring significantly.

As the first planet has already grown to ${\sim}10\,M_{\oplus}$ at $t=0.6$
Myr, it feels sufficient gravitational torques from the gas disc and therefore
leaves the ring via migration. After the planet migrates away from the massive
ring, its accretion rate drops rapidly. Even though it can still accrete from
the mass flux leaked from the ring, as discussed in previous literature, the
growth of the planet in a smooth pebble-flowing disc is slow and inefficient
at such a large distance (e.g., Ormel, 2017; Lin et al., 2018). In Fig. 6, we
stack the mass budget in different components as the function of time. The
mass contribution accreted from the smooth disc is marked in red, which is
tiny compared with those accreted from the ring. Meanwhile, as the first
planet has left the ring, a smaller embryo within the ring will get the
opportunity to grow. Since the external pebble flux that supplies the ring is
fixed, the evolution becomes cyclical. In this cyclical phase III, the biggest
planet inside the ring efficiently grows by pebble accretion and then migrates
inward and leaves the ring behind. As shown in Fig. 6, regulated by the pebble
accretion, the ring mass does not keep increasing but oscillates around a
"saturated" value. The majority of the upcoming mass flux ends up in the
largest embryos by pebble accretion from the ring. The long-lived ring acts as
a planet factory where planets sequentially form and migrate towards the inner
region of the disc. In the following section, we demonstrate the physical
model to interpret the planet masses and the ring mass.

Figure 5:  Top: Total pebble accretion flux of all planetesimals (blue line).
The green region shows the contribution from the most massive body still
inside the simulation box (beyond $r_{\rm in}$). Middle: Fraction of the
contribution from the most massive body over the total pebble accretion mass
flux. Bottom: Mass of the ring (blue) and mass of the planets that have once
been the most massive bodies (colours). When the planet migrates past the
inner boundary of the simulation box, the line terminates. The vertical line
at 2.5 Myr marks a collision between two planets. The evolution can be divided
into three phases, see main texts.

###  3.2 Mass of ring and planet in saturation phase

As we discussed in the default model, the ring mass gets into steady states
after the first planet forms and migrates away from the ring, which is marked
as phase III in Fig. 5. We find this is a robust outcome for the vast majority
of our simulations with type-I migration. We quantify the mass of the planets
and the saturated ring in this section.

In most of our model, the mass of planetesimals formed directly by streaming
instability is always small compared with the ring mass and the mass accreted
to the planetesimals by pebble accretion (see Fig. 6). In other words, the
mass flux fed into the ring is the same as the mass flux that planets are
accreting from the ring. Therefore, the ring mass saturates at a value $M_{\rm
ring,ss}$ at which the total accretion rate of planets is equal to the mass
flux fed to the ring. The mass of the planet is then determined by the balance
between planet accretion and migration, i.e, when the planet grows massive
enough, it will migrate away from the ring center and the accretion rate will
rapidly drop because of the lack of feeding materials. By assuming there is
always one planet with mass $m_{\rm rp}$ inside the ring, the typical
accretion timescale of this planet is then $m_{\rm rp}/\dot{M}_{\rm net}$. And
the timescale for this "representative" planet to migrate away from the ring
is $\Delta t_{\mathrm{ring}}=2w_{\rm ring}/r_{0}\times t_{\rm mg}(m=m_{\rm
rp})$. By equating these two timescales, we obtain the representative planet
mass

$\displaystyle m_{\rm rp}$ $\displaystyle=f_{\rm rp}\sqrt{\frac{\dot{M}_{\rm
net}}{\Sigma_{\rm g}r_{0}^{2}\Omega_{K}}\frac{1}{f_{\rm mg}}\frac{w_{\rm
ring}}{r_{0}}}h_{\rm g}M_{\star}$ (32)
$\displaystyle=18\,M_{\oplus}\times\left(\frac{f_{\rm
rp}}{1.4}\right)\left(\frac{f_{\rm
mg}}{1}\right)^{-\frac{1}{2}}\left(\frac{\dot{M}_{\rm
net}}{50\,M_{\oplus}\,\rm Myr^{-1}}\right)^{\frac{1}{2}}$
$\displaystyle\times\left(\frac{\Sigma_{\rm g}}{3.6\,\rm
g\,cm^{-2}}\right)^{-\frac{1}{2}}\left(\frac{w_{\rm ring}}{3.38\,\rm
au}\right)^{\frac{1}{2}}\left(\frac{r_{0}}{74.2\,\rm
au}\right)^{-\frac{3}{4}}\left(\frac{h_{\rm
g}}{0.07}\right)\left(\frac{M_{\star}}{M_{\odot}}\right)^{\frac{3}{4}}$

where a fit constant $f_{\rm rp}=1.4$ is used.

As discussed above, when the ring mass saturates, the accretion rate of the
planet is equal to the ring mass flux. It is shown in Fig. 4 that the pebble
accretion operates in the 2D regime for a planet with the mass of $m_{\rm
rp}$. Therefore, by equating the 2D pebble accretion rate Eq. (18) with
$\dot{M}_{\rm net}$, we obtain a typical surface density

$\displaystyle\Sigma_{\rm tp}$ $\displaystyle=\frac{1}{2.9}\frac{\dot{M}_{\rm
net}}{(Stq_{\rm rp})^{\frac{2}{3}}\Omega_{K}r_{0}^{2}}=\frac{0.35(f_{\rm
mg}\Sigma_{g})^{\frac{1}{3}}\dot{M}_{\rm net}^{\frac{2}{3}}}{(f_{\rm rp}h_{\rm
g}\mathrm{St})^{\frac{2}{3}}(w_{\rm ring}GM_{\star})^{\frac{1}{3}}}$ (33)
$\displaystyle=0.13\,{\rm g\,cm^{-2}}\times\left(\frac{f_{\rm
rp}}{1.4}\right)^{-\frac{2}{3}}\left(\frac{f_{\rm
mg}}{1}\right)^{\frac{1}{3}}\left(\frac{\dot{M}_{\rm net}}{50\,M_{\oplus}\,\rm
Myr^{-1}}\right)^{\frac{2}{3}}$ $\displaystyle\times\left(\frac{\Sigma_{\rm
g}}{3.6\,\rm g\,cm^{-2}}\right)^{\frac{1}{3}}\left(\frac{w_{\rm
ring}}{3.38\,\rm
au}\right)^{-\frac{1}{3}}\left(\frac{St}{0.01}\right)^{-\frac{2}{3}}\left(\frac{h_{\rm
g}}{0.07}\right)^{-\frac{2}{3}}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-\frac{1}{3}}$

where we have substituted Eq. (32) for $q_{\rm rp}=m_{\rm rp}/M_{\star}$ and
$\Omega_{K}=\sqrt{GM_{\star}/r_{0}^{3}}$ in the first line. As shown in Eq.
(11), this corresponds to a ring mass of

$\displaystyle M_{\rm ring,ss}$ $\displaystyle=(2\pi)^{1.5}r_{0}w_{\rm
ring}f_{\rm rs}\Sigma_{\rm tp}=5.5f_{\rm rs}\frac{(f_{\rm
mg}\Sigma_{g})^{\frac{1}{3}}(w_{r}\dot{M}_{\rm
net})^{\frac{2}{3}}r_{0}}{(f_{\rm
rp}h_{g}\mathrm{St})^{\frac{2}{3}}(GM_{\star})^{\frac{1}{3}}}$ (34)
$\displaystyle=34\,M_{\oplus}\times\left(\frac{f_{\rm
rs}}{1.8}\right)\left(\frac{f_{\rm
rp}}{1.4}\right)^{-\frac{2}{3}}\left(\frac{f_{\rm
mg}}{1}\right)^{\frac{1}{3}}$ $\displaystyle\times\left(\frac{\dot{M}_{\rm
net}}{50\,M_{\oplus}\,\rm
Myr^{-1}}\right)^{\frac{2}{3}}\left(\frac{\Sigma_{\rm g}}{3.6\,\rm
g\,cm^{-2}}\right)^{\frac{1}{3}}\left(\frac{r_{0}}{74.2\,\rm au}\right)$
$\displaystyle\times\left(\frac{w_{\rm ring}}{3.38\,\rm
au}\right)^{\frac{2}{3}}\left(\frac{h_{\rm
g}}{0.07}\right)^{-\frac{2}{3}}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-\frac{1}{3}}\left(\frac{\rm
St}{0.01}\right)^{-\frac{2}{3}}$

where $f_{\rm rs}=1.8$ is a numerical prefactor to correct the influence
caused by the no-uniform surface density inside ring. The predicted value
$M_{\rm ring,ss}$ and the measured value $M_{\rm ring,ms}$ are listed in Table
3 for each run. Equation (34) approximates the measured ring mass very well.
Substituting Eq. (34) into Eq. (32), we get

$\displaystyle m_{\rm rp}$ $\displaystyle=0.3\times f_{\rm
rs}^{-\frac{3}{4}}f_{\rm rp}^{\frac{3}{2}}f_{\rm
mg}^{-\frac{3}{4}}\left(\frac{M_{\rm ring,ss}}{\Sigma_{\rm
g}r_{0}^{2}}\right)^{\frac{3}{4}}{\rm St}^{\frac{1}{2}}h_{\rm
g}^{\frac{3}{2}}M_{\star}$ (35)
$\displaystyle=18M_{\oplus}\times\left(\frac{f_{\rm
rs}}{1.8}\right)^{-\frac{3}{4}}\left(\frac{f_{\rm
rp}}{1.4}\right)^{\frac{3}{2}}\left(\frac{f_{\rm
mg}}{1}\right)^{-\frac{3}{4}}\left(\frac{M_{\rm
ring,ss}}{30M_{\oplus}}\right)^{\frac{3}{4}}$
$\displaystyle\times\left(\frac{\Sigma_{\rm
g}r_{0}^{2}}{0.002M_{\odot}}\right)^{\frac{3}{4}}\left(\frac{\rm
St}{0.01}\right)^{\frac{1}{2}}\left(\frac{h_{\rm
g}}{0.07}\right)^{\frac{3}{2}}\left(\frac{M_{\star}}{M_{\odot}}\right)$

which is independent of the unknown feeding mass flux. Therefore, by assuming
that an ALMA ring is in steady state, if we can measure the mass of the ring,
the Stokes number, the gas surface density and scale height, we can obtain the
typical planetary core formed and delivered from the ring.

###  3.3 Similarities and differences between CR and PB models

Even though we focus on the CR model in our default run, the planet factory
picture broadly holds for any ring that meets the following conditions:

  1. 1.

active planetesimal formation inside the pebble ring;

  2. 2.

pebble inside the ring have low approaching velocity towards these formed
planetesimals;

  3. 3.

the ring is supplied by an external mass reservoir and therefore is long-
lived;

  4. 4.

planets migrate inwards away from the ring.

Conditions (i) and (ii) are required to warrant efficient pebble accretion. In
other words, the criterion for planetesimal formation is consistent with
massive volume pebble densities and the planetesimals can act as seeds of
planet embryos. The second condition ensures pebble accretion falls in the
settling regime, as otherwise (for strong headwinds) pebble accretion is
severely suppressed (Visser & Ormel, 2016). The longevity of the ring leaves
enough time for planet growth; and, finally, the inward migration prevents a
too massive planet to occupy the ring, which would otherwise destroy it, and
allows a new embryo in the ring to grow. (see Sect. 3.4 for further
assessment).

Figure 6:  Evolution of the mass budget of different components in the default
model. The blue region shows the ring mass. The brown region shows the
accumulated mass of planetesimals formed by streaming instabilities. The mass
accreted to the planetesimals and protoplanets by pebble accretion from the
ring (orange) and from the extended smooth disc component (red) are stacked on
top. Dashed region marks the accumulated mass of planets that have ever been
the largest planets inside the ring. Gray horizontal line marks the ring mass
corresponding to $Z=1$. Green dashed line gives the accumulated mass fed to
the ring, $\dot{M}_{\rm net}t$. The stacked value is a bit higher than
$\dot{M}_{\rm net}t$ due to the accretion of the planet downstream. Figure 7:
Same as Fig. 6 but for run pb-default. The evolution is similar between the
clumpy ring and the pressure-bump-supported ring. A larger initial ring mass
is required to reach $Z=1$ to start planetesimal formation. Yet, the mid-plane
dust-to-gas ratio of the ring drops below unity during the evolution,
hindering the planetesimal formation.

This picture applies to a planetesimal-forming ring held by a pressure bump as
well. To illustrate this, in the pb-default run, we use the same parameters as
the cr-default but assume that a permanent pressure bump supports the ring.
The mass budget of run pb-default is shown in Fig. 7. The ring initially grows
when the planetesimals are still small when $t<0.3$ Myr. After that, most of
the incoming mass flux ends up in the largest embryos by pebble accretion. The
ring enters a steady state and the evolution is qualitatively the same as in
the CR model.

An obvious difference between the CR model and the PB model is the ring mass.
In order to reach the mid-plane dust-to-gas ratio $Z=1$, for pebbles with $\rm
St=0.01$, one magnitude higher initial ring mass is needed in the PB model
compared with the CR model. Since the initial ring mass exceeds the predicted
saturated ring mass, the planet grows so rapidly that the mid-plane dust-to-
gas ratio drops below unity again, see Fig. 7. Planetesimal formation is,
therefore, less efficient compared to the clumpy ring model (compare the
"planetesimal" budget in Fig. 6 and Fig. 7). On the other hand, protoplanets
migrate away from the ring when they reach the representative mass Eq. (32),
at which point their accretion rate also exceeds the net feeding mass flux.
The migration prevents the planet from consuming the ring, and therefore, the
ring mass remains larger than the predicted ring mass, Eq. (34), during the
evolution as marked in Table 3.

In the CR model, pebbles reside in a thin disc whose scaleheight is
independent of St and $\delta_{\rm z}$. Thus, the timescale to reach $Z=1$,
$t_{\mathrm{load}}$, chiefly depends on $\dot{M}_{\rm net}$. However, for PB-
supported rings the pebble scaleheight and therefore $t_{\mathrm{load}}$ are
sensitive to St and $\delta_{z}$ (Eq. (12)). For pb-default ($\rm St=0.01$ and
$\rm\delta_{\rm z}=10^{-3}$), this wait amounts to $1.61$ Myr corresponding to
a ring mass of $160M_{\oplus}$. Discs with lower Stokes numbers (or higher
$\delta_{z}$) require an even longer loading time, rendering the ring very
massive but still sub-critical ($Z<1$). Planetesimal formation may never
happen within the disc lifetime. In contrast, $t_{\mathrm{load}}$ is much
shorter with a higher St or lower $\delta_{z}$ for which planetesimal
formation can be triggered at lower surface densities (a less massive ring).
For run pb-St100 ($\rm St=0.1$) and pb-a4 ($\rm\delta_{\rm z}=10^{-4}$), a
$53M_{\oplus}$ ring suffices to reach the dust-to-gas ratio $Z=1$. The lower
initial ring mass means that fewer pebbles are available for planetesimals to
accrete. Yet, the higher St and thinner pebble scaleheight enhance the pebble
accretion efficiency (Eq. (18)). Hence, the formation of the first planet is
faster (shorter $t_{\rm 1}$ in Table 3) in runs pb-St030, pb-St100 and pb-a4.
Whereas at $\dot{M}_{\rm ext}=100M_{\oplus}\,\rm Myr^{-1}$, Stokes numbers
${>}0.02$ are catastrophic for clumpy rings due to the large leaking mass flux
(see Sect. 3.5.2), PB-supported rings can accommodate high St particles.

###  3.4 Influence of type-I migration

For both CR and PB models, type I migration plays a crucial role. Massive
Planets migrate inward and eventually leave the ring, which terminates their
rapid pebble accretion from the ring but promotes the longevity of the ring.
To quantitatively understand the effect of migration on the evolution, we
change the type-I migration prefactor $f_{\rm mg}$ (Eq. (27)) in the group of
runs named with mg. We conduct runs with slower migration ($f_{\rm mg}=0.5$;
cr-05nmg) and with faster migration ($f_{\rm mg}=2$; cr-20nmg).

As expected from our physical model, the typical planet masses are higher
(lower) in the slow (fast) migration run, and the saturated ring mass is lower
(higher) accordingly, see Table 3. The outcome is intuitive. The migration
speed decides the duration of the planet inside the ring, in which the
majority of the core mass grows by the efficient pebble accretion. Thus, a
slower migration speed allows more time for a single embryo to grow and
consume the ring more, and vice versa.

In addition, in runs labeled nmg (no migration), we assume that the ring
coincides with the zero-torque location. In this case, the planet formed
inside the ring does not undergo type I migration, i.e., $f_{\rm mg}\to 0$.
Such a situation can happen especially in the PB scenario, where the pressure
bump arises from a local gas density variation and the co-rotation torque and
Lindblad torque cancel each other (e.g., Masset et al., 2006). The evolution
of all nmg runs are similar. We show the evolution of the mass budget of the
cr-nmg run in Fig. 8 as an example. For the first $0.6$ Myr, the evolution is
similar to the cr-default run. This is because the planet mass in phase I and
the beginning part of phase II are still small. However, after $0.6$ Myr, the
evolution of the no migration runs is qualitatively different from the default
runs since the most massive planet remains in the ring. Rather than leaving
the ring, this planet keeps accreting and suppresses the growth of other
embryos by scattering them away. Pebble accretion of the most massive planet
becomes so dominant that it consumes the ring.

Figure 8:  Same as Fig. 5 but for the cr-mng run. The biggest planet dominates
the pebble accretion and consumes the ring.

###  3.5 Parameter study

####  3.5.1 Initial mass of planetesimals

Figure 9:  Output data of runs with different initial masses from Ceres mass
to Mars mass from Table 3. Top: Time when the first (sixth) planet past
$r_{\rm in}$, $t_{1}$ ($t_{6}$), is marked in blue (orange). Bottom: Dots mark
the average and maximum mass of the six planets in each run when they pass
$r_{\rm in}$. Green diamonds mark the measured ring masses between $t_{1}$ and
$t_{6}$. Light orange and green lines indicate the predicted representative
planet mass Eq. (32) and the saturated ring mass Eq. (34) respectively.

Since the mass of the planetesimals formed by SI is still under debate, we
test a number of initial masses $m_{\rm 0}$ in our simulations, ranging from
Ceres ($1.5\times 10^{-4}M_{\oplus}$) to Mars ($0.1M_{\oplus}$), see Table 2.
In Fig. 9, we visualize the growth time (upper panel) and the measured masses
of planets and rings (lower panel) in runs with different initial masses.
Interestingly, the simulation results almost show no dependence on the initial
mass.

In contrast to pebble accretion in smooth discs, inside rings, the pebble
density is very high and the relative velocity between pebbles and
planetesimals is small. Therefore, planetesimals efficiently accrete pebbles
once they are born, already at Ceres-masses. For our initial ring mass $M_{\rm
ring,0}=11.6M_{\oplus}$, the e-folding timescale of pebble accretion in the 3D
regime is less than $0.1\,\mathrm{Myr}$ and independent of the mass of the
planetesimal (initial point in Fig. 4). For runs with an initial mass $m_{\rm
0}<M_{\rm Pluto}$, pebble accretion starts in the 3D regime. Due to the short
e-folding growth timescale in the 3D regime, the duration of the planet in the
3D regime is short compared with the total growth timescale. In the 2D pebble
accretion regime, the growth timescale increases with planet mass (Eq. (31)).
Therefore, the majority of the growth time is spent in the final doubling of
the planet’s mass. Consequently, the time when the first planet passes the
inner boundary after the first planetesimal formed, $t_{\rm 1}$ (blue
triangles in Fig. 9), hardly depends on the initial planetesimal mass.

For run cr-mMoon and cr-mMars (the rightmost two columns in Fig. 9) the
initial mass is so massive that pebble accretion starts in the 2D regime,
during which the pebble accretion timescale is scaling with the mass of
protoplanet. Thus $t_{\rm 1}$ shows a weak decreasing trend in these two
columns. On the other hand, for smaller planetesimals, the collision time-
scale $t_{\rm coll}\sim\rho_{\bullet}R_{\rm plts}/\Sigma\Omega_{K}$ is
shorter. Therefore, planetesimal mergers are more important in runs with
smaller $m_{\rm 0}$ and speed up the growth. Finally, we also see some
stochasticity due to the N-body interaction which results in a spread in
quantities like $t_{1}$ and $t_{6}$. Altogether, the time to form 6 planets
that can pass the inner boundary $t_{\rm 6}$ (orange triangles in Fig. 9)
shows no dependence on initial mass at all.

In the lower panel of Fig. 9, we mark the measured average and maximum mass of
the six planets that passed the inner boundary of the simulation domain with
blue and orange dots, and the measured ring mass in green. These output masses
of both rings and planets show no dependency on $m_{\rm 0}$ as well,
consistent with our prediction in Eq. (32) and Eq. (34) (horizontal lines). In
the saturation phase, the evolution is determined by the most massive planet,
which always dominates the pebble accretion. And the mass of the most massive
planet is self-regulated by migration. The initial planetesimal has a
negligible influence on this self-regulation cycle.

Although we cannot conduct experiments with even smaller planetesimal masses
due to computational constraints, the independence of the result on the
initial planetesimal mass should hold also for lower masses. The reason that
pebble accretion is so effective is that the eccentric velocity $ev_{K}\ll
v_{\ast}$. As shown in Fig. 10, planetesimals are readily excited to
eccentricities similar to the Hill value of the largest bodies $e\sim(q_{\rm
1st}/3)^{1/3}$ to enter the dispersion dominated regime (for which $i\sim
e/2$) after $t=0.1\,\rm Myr$. In the dispersion dominated regime, however, the
efficiency of the viscous stirring strongly decreases with eccentricity (e.g.,
Ida & Makino, 1993). Balancing $t_{\mathrm{vs}}$ with the growth timescale
(constant for pebble accretion in the 3D regime), an equilibrium eccentricity
for the planetesimals follows, which is not much greater than the Hill
eccentricity $e_{\rm H}\equiv v_{\rm H}/v_{K}=(q_{\rm 1st}/3)^{1/3}$. We then
obtain that $ev_{K}\sim v_{H}<v_{\ast}$ and the eccentricity of the more
massive protoplanets will still be lower due to dynamical friction. Viscous
stirring, therefore, is rather ineffective. As all velocities scale the same
with planet mass, the results are independent of the initial planetesimal
mass.

Figure 10:  Eccentricity (solid) and inclination (dashed) evolution in run cr-
mCeres. The red line represent the most massive body and the black lines
indicate rms values of the planetesimals (other bodies except the most massive
one). The gray solid line marks the critical velocity threshold for pebble
accretion $v_{\ast}$ and the gray dashed line marks the Hill eccentricity
$e_{\rm H}\equiv v_{\rm H}/v_{K}$ of the largest body. Eccentricities follow
$e_{\rm H}$ in the dispersion dominated regime and remain much lower than
$v_{\ast}$.

####  3.5.2 External mass flux and Stokes number

In our simulations, we also vary the mass flux and Stokes number. Overall, the
predicted and measured masses of planets and rings (column 6, 7, and 9, 10 of
Table 3) show good correspondence, see Fig. 11 for a visualization.

Figure 11:  Output data of runs with different Stokes number and external mass
flux from Table 3. Blue markers belong to the run with the net mass flux
feeding the ring $\dot{M}_{\rm net}=50\,M_{\oplus}\rm\,Myr^{-1}$. Orange
markers are runs with fixed external mass flux $\dot{M}_{\rm
ext}=100\,M_{\oplus}\rm\,Myr^{-1}$. Maximum planet masses $m_{\rm max}$ are
marked as "+" and measured ring masses $m_{\rm ring,ms}$ are marked as "x".
Solid and dash lines indicate the predicted representative planet mass Eq.
(32) and saturated ring mass Eq. (34) respectively.

The dependence of the evolution on the external mass flux is intuitive, i.e.,
the higher the mass flux, the more available pebbles, the denser the ring, and
the faster the growth of the planet.

The Stokes number, the aerodynamic size of pebbles, plays a crucial role in
pebble accretion. In our model, there are three aspects controlled by the
Stokes number:

  1. 1.

the settling of pebbles

  2. 2.

the leaking mass flux Eq. (22)

  3. 3.

the planetesimal formation rate Eq. (15)

In our simulations, we always start with dust-to-gas ratio unity at the ring
center. For the CR model the Stokes number will not directly influence the
initial ring mass, because the pebble scaleheight is set by the
Kelvin–Helmholtz instability (Eq. (13)). Independent of the assumption on the
ring, for a given pebble surface density, the pebble accretion rates are
higher for pebbles with large St in both 3D and 2D regimes, see Eq. (18).
However, the leaking mass flux also increases with the Stokes number, because
larger pebbles have a higher radial drift than smaller pebbles.

We investigate the role of St in CR rings by running two groups of
simulations. First, in run cr-St001 to cr-St100 (blue markers in Fig. 11), we
keep the net feeding flux $\dot{M}_{\rm net}$ constant to isolate the
influence caused by the leaking mass in our simulation. Since the leaking mass
flux is linear with Stokes number in the CR model, we change the external mass
flux accordingly to ensure that the net mass flux feeding the ring
$\dot{M}_{\rm net}=\dot{M}_{\rm ext}-\dot{M}_{\rm leak}$ is always
$50M_{\oplus}\,\rm\,Myr^{-1}$. Therefore, the external mass flux is only
$55M_{\oplus}\,\rm Myr^{-1}$ for run with $\rm St=0.001$ while it increases to
$550M_{\oplus}\,\rm Myr^{-1}$ for run with $\rm St=0.1$. Because of the low
pebble accretion efficiency in the smaller Stokes case, the ring grows much
more massive (blue dashed line in Fig. 11), in line with our prediction inEq.
(34). From both the $t_{\rm 1}$ and $t_{\rm 6}$ indicators (see Table 3), it
is clear that growth of a planet is slower in runs with lower St. However, as
we predict in Eq. (32), there is no significant dependence between the planet
mass and St (horizontal blue line in Fig. 11). This is because massive planets
migrate inward and leave the rapid-growth ring region. Termination of the
growth depends on the migration behavior, which is determined by the gas
properties of the disc rather than the pebble properties.

Then, in run cr-St001F100 to cr-St018F100 (orange markers in Fig. 11), we keep
the external mass flux fixed at $100\,M_{\oplus}\,\rm Myr^{-1}$ and vary the
Stokes number from $0.001$ to $0.018$. Now, both the feeding mass flux and the
leaking mass flux change with Stokes number. When the Stokes number is low,
the leaking mass flux is low as well and a larger portion of the external mass
flux can be fed to the ring. Yet, the growth rate is still low because of the
low pebble accretion efficiency. In contrast, when St is high, more pebbles
will leak away from the ring to the inner disc. For pebbles whose Stokes
number exceeds $0.02$, the leaking mass becomes larger than
$100\,M_{\oplus}\,\rm Myr^{-1}$, which means that the clumpy ring will
"vaporize" or never form in the first place. In Fig. 11, the measured planet
and ring masses rapidly drop when the Stokes number approaches $0.02$, as
predicted by the orange lines. Therefore, a larger Stokes number is not
necessarily beneficial to planet formation. For a fixed external mass flux of
$100\,M_{\oplus}\,\rm Myr^{-1}$, the growth rate peaks for a Stokes number
around $\rm St=0.01$, which shows smallest $t_{\rm 1}$ among all $\dot{M}_{\rm
ext}=100\,M_{\oplus}\,\rm Myr^{-1}$ runs (see Table 3).

The Stokes number in our model also influences the planetesimal formation rate
since the planetesimal formation timescale scales with the settling time
scale, see Eq. (15). Thus, pebble rings with larger Stokes numbers can form
planetesimals faster. For the default value of the planetesimal formation
efficiency parameter ($\epsilon=10^{-3}$), pebble accretion and leaking always
dominate the net mass flux fed to the ring while the flux represented by
planetesimal formation is minor.

####  3.5.3 Planetesimal formation efficiency

Figure 12:  Same as Fig. 6 but for run cr-e2. The large amount of direct
planetesimal formation (brown) limits the mass of the pebble ring
significantly.

We investigate the planetesimal formation efficiency by varying the prefactor
$\epsilon$ in Eq. (15). We test a higher $\epsilon=0.01$ in run cr-e2 and
several lower $\epsilon$ runs. As we show in Fig. 6, in the default model
where $St=0.01$, the fraction of pebbles that collapse into planetesimals (by
SI rather than ending up in planets by pebble accretion) is small. Therefore,
the mass loss of the ring due to planetesimal formation is negligible compared
with pebble accretion. This no longer holds for the run with higher
$\epsilon$, in which case planetesimal formation will significantly limit the
pebble ring mass as shown in Fig. 12. As a result, the measured saturated ring
mass in run cr-e2 is only $22.1M_{\oplus}$, significantly lower than the
predicted value $34\,M_{\oplus}$ in Eq. (34). Since planetesimal mergers is
rare at these large radii, the large planetesimal budget does not help with
the growth of the planet. The growth time of the first planet $t_{\rm 1}$ is
similar between cr-e2 and cr-default. Conversely, since the saturated ring
mass is lower in phase III, growth of later forming planets requires more
time. The average interval time $\langle\Delta\,t\rangle$ on which the ring
spawns planets doubles in the higher $\epsilon$ run (see Table 3). In the same
amount of time, more smaller planetesimals are born but fewer big planets
form. In particular, we note that further increasing the planetesimal
formation efficiency is catastrophic for the survival of the clumpy ring. In
JO21, even without accounting for pebble accretion, rings cannot survive in
runs where $\epsilon\gtrsim 0.1$, since the rapid planetesimal formation
drives the ring midplane dust-to-gas density ratio below unity. See Sect. 4.1
for further discussion.

On the other hand, for lower planetesimal formation efficiency (cr-e4, cr-e5),
there is almost no difference with the default model. Ring mass removal is
dominated by pebble accretion. Therefore, the planetesimal formation has a
negligible impact on the ring’s evolution.

####  3.5.4 Location of the ring

Finally, we investigate how the results depend on the location of the ring.
Previous literature has highlighted the influence of the disc radius on the
planet formation (e.g. Visser & Ormel, 2016; Bitsch et al., 2019; Morbidelli,
2020; Chambers, 2021; Voelkel et al., 2021; Jang et al., 2022). We setup a
ring at $30$ au in run cr-r30 and at $154$ au in run cr-r154, where the
leaking mass flux change to $40\,M_{\oplus}\,\rm Myr^{-1}$ and
$60M_{\oplus}\,\rm Myr^{-1}$ respectively (Eq. (22)). Other disc parameters
are all the same as the cr-default run.

Interestingly, both the planet mass and planet growth timescale show little
dependence on the location of the ring (see Table 3), in contrast to
Morbidelli (2020) who report a significantly efficient growth when the ring is
closer to the host star. One of the reasons for the insensitivity on disc
radius is that the saturated ring mass, as derived in Eq. (34), scales with
the orbital radius. In the saturated phase, the net pebble flux fed to the
ring will eventually be accreted by the planet. Since the planet can accrete
the ring more efficiently at closer orbits, pebbles remain a shorter time in
the ring. Therefore, in the inner disc, the ring is less massive than the ring
in the outer region, if we assume the same ring width. And the pebble
accretion rate of the planet always reaches the net mass flux $\dot{M}_{\rm
net}$ at the end. On the other hand, the timescale required to migrate through
the ring width $\Delta t_{\mathrm{ring}}$ is independent of radius for our
choice of the disc and ring model parameters, leaving a similar time for the
planet to accrete pebbles from the ring. Therefore, as we derived in Eq. (32),
the planet mass always approximates $\dot{M}_{\rm net}\Delta
t_{\mathrm{ring}}$. We find the growth of planets in the clumpy ring model
shows only a weak preference for lower radii.

###  3.6 Post-ring evolution of planets

In the CR model, planets that have migrated and left the ring can still
accrete pebbles downstream by virtue of the leaking mass flux. Compared with
pebble accreted from the ring, the fraction of pebbles accreted from the
leaking mass flux is nevertheless small around the ring region (the red region
above the green dash line in Fig. 6). We briefly discuss the subsequent
evolution of planets after they have left the simulation domain.

Figure 13:  Trajectories of the first four planets larger than $5M_{\oplus}$
that passed the inner simulation domain in cr-default. Dots mark the evolution
tracks of planets when they are still inside the N-body simulation domain.
Solid lines in the same colours indicate the predicted trajectories following
Eq. (54). Purple line "B18" marks the pebble-isolation mass from Bitsch et al.
(2018). Brown line "K18" marks the gap opening mass from Kanagawa et al.
(2018). The time when planets pass the inner boundary $r_{\rm in}$ is labeled
accordingly. The contours in the background show the expected gap depth in the
gas disc carved by the planet.

We show the evolution tracks of the first four planets larger than
$5M_{\oplus}$ that passed the simulation domain in cr-default in Fig. 13. The
dots record the evolution track when these planets are still in the N-body
simulation. Solid lines indicate their corresponding predicted trajectories
based on Eq. (54) (see Appendix B for the mathematical expressions). The blue
contours represent the expected gap depth carved by the planet at
corresponding radii (Kanagawa et al., 2018)

$\frac{\Sigma_{\rm min}}{\Sigma_{\rm un,p}}=\frac{1}{1+0.04K}$ (36)

where

$K=\left(\frac{m_{\rm p}}{M_{\star}}\right)^{2}\left(\frac{h_{\rm g}}{a_{\rm
p}}\right)^{-5}\alpha^{-1}$ (37)

These planets keep growing by pebble accretion and migrating inward until they
reach the pebble isolation mass, at which the gas gap carved by the planet is
deep enough to halt the pebble flux at the outer edge and stop the pebble
accretion (e.g., Lambrechts & Johansen, 2014; Bitsch et al., 2018). In Fig.
13, we plot the pebble isolation mass as the purple line according to Eq. (11)
in Bitsch et al. (2018). Planets that reach the pebble isolation mass then
undergo runaway gas accretion and grow to gas giants (e.g. Bitsch et al.,
2019). Meanwhile, the migration speed also slows down as the surface density
of the bottom of the gap decreases (Kanagawa et al., 2018). The brown line
indicates the gap opening mass

$M_{\rm gap}=60\left(\frac{\alpha}{10^{-3}}\right)^{0.5}\left(\frac{h_{\rm
g}}{0.07}\right)^{2.5}M_{\oplus}$ (38)

(Eq. (23) in Kanagawa et al., 2018), which marks a transition from type-I
migration to so-called type-II migration.

In addition, we test a lower disc turbulence ($\delta_{\rm z}=10^{-4}$) in run
cr-a4, in which planets grow faster and larger (see Table 3). This outcome is
not because of the planet accretion inside the ring, since the pebble scale
height of the clumpy ring is independent of $\delta_{\rm z}$. Yet, the low
disc turbulence does help with the post-ring evolution. When the turbulence is
suppressed, the disc pebble scale height is smaller. Therefore, the pebble
accretion downstream is significantly enhanced. The largest planet reaches a
mass of $50.7\,M_{\oplus}$ in cr-a4 when it arrived at the inner boundary of
the N-body simulation domain. At this mass, the planet may already have
reached the condition for gap opening if we assume that the low dust vertical
turbulence is indeed caused by the low gas viscosity (although $\alpha$ and
$\delta_{\rm z}$ can be different as we mentioned below Eq. (12)). The pebble
isolation mass and gap opening mass are only 20 $M_{\oplus}$ at the ring
location with $\alpha=10^{-4}$. Therefore, this planet should have opened a
gap in the disc already and started gas accretion.

###  3.7 Planetesimal belt

(a) $t=0.5$ Myr

(b) $t=1.5$ Myr

(c) $t=2.5$ Myr

Figure 14: Mass distribution of bodies at $t=0.5$, 1.5, 2.5 Myr in cr-default
run. Top row shows the histogram in terms of the number of bodies, and the
bottom row shows the histogram in terms of accumulated mass in each bin. Blue
bins mark the bodies classified as planets in our simulation (mass larger than
$m_{\rm crit}$ and used for statistics in Table 3), all planets eventually
migrate through the inner boundary of the simulation domain. Others are
classified as planetesimals and coloured in orange. For both planet and
planetesimals, the dark part remains in the N-body simulation domain, and the
light part has already passed the inner boundary. In cr-default, all
planetesimals start with Pluto mass, which is the cutoff on the left of each
panel.

As proposed in JO21, the large amount of planetesimals formed in the pebble
ring may evolve into a debris disc with the dispersal of the gas disc. In
JO21, we neither considered any interaction between pebbles and planetesimal
nor accounted for the dynamics of planetesimals. To balance the external
pebble flux and the ring mass loss, the planetesimal formation efficiency is
set to be $\epsilon=0.01$ in the default run, which results in a total mass up
to $100M_{\oplus}$ of planetesimals over a time-span of $5$ Myr.

In Fig. 14, we show the mass distribution of the bodies in cr-default at
different times. We colour the planet population, whose $t_{\rm
mg}(r=r_{0})\leq 1\rm Myr$, in blue and planetesimals , whose $t_{\rm
mg}(r=r_{0})>1\rm Myr$, in orange as we defined in Sect. 3.1. The majority of
the ring mass ends up in planets, while most of the planetesimals remain of
low mass. The shallow region in the histogram are those bodies that have
passed the inner simulation domain boundary. All planets migrate inward and
leave the ring region eventually as shown in Fig. 3. Due to their low mass,
planetesimals will stay around the ring location.

To understand the evolution of the planetesimal belt, we show the time-
averaged "radial profile" of planetesimals’ number density during [0,1],
[1,2], [2,3] Myr in Fig. 15. As time evolves, the radial profile of the
planetesimal belt becomes asymmetrical, characterized by a steep rise and a
gradual decline. Planetesimals experience both resonant shepherding and
outward dynamical scattering by the inward-migrating planets (e.g., Batygin &
Laughlin, 2015). As time evolves, planetesimals interior to the ring can be
shepherded by migrating planets through mean resonance (e.g., Shibata & Ikoma,
2019) and leave the inner boundary of the simulation domain. In contrast,
planetesimals are more likely scattered outward by the inward-migrating
planets and stay exterior to the ring, leading to a sharper inner edge and a
flatter outer edge.

Figure 15:  The time-averaged planetesimal number density radial profile
during 0–1, 1–2, and 2–3 Myr in cr-default run. Radial profile of the pebble
ring profile is shown in gray. As time evolves, more planetesimals are
scattered outward by the planet that migrated interior to the ring. Horizontal
segments mark the $3\sigma$ width from the pebble ring center by assuming half
Gaussian distributions for two edges of the planetesimal belt respectfully.
The outer edge is shallower than the inner.

##  4 Discussion

###  4.1 Caveats

The model presented in this work employed a number of assumptions and
simplifications.

First, in our simulation, we aim to study the local evolution around the ring
region. In order to isolate the effect caused by the uncertain disc evolution
(See Miotello et al., 2022, and references therein), we keep the mass flux
constant and ignore the global evolution of the disc. In reality, the pebble
flux could decay with time. For discs where the dispersal timescale exceeds
${\sim}\rm Myr$, the steady state solution we derived is valid. Yet, if the
disc evolves on timescales equal or shorter than ${\sim}\rm Myr$, the ring
mass and the ensuing planet masses will also decrease with the decaying pebble
mass flux as predicted by Eq. (32) and Eq. (34), limiting efficiency of the
planet factory. In future studies, a more thorough model of disc evolution
including pebble drift, coagulation, and fragmentation needs to be considered,
including, for example, the pebble-predictor prescription of Drążkowska et al.
(2021).

The planetesimal formation prescription in Eq. (15) means that streaming
instability and planetesimal formation happen when the mid-plane dust-to-gas
ratio $>\\!1$, after which the clumps will exist for a settling time-scale
(Shariff & Cuzzi, 2015). We choose a planetesimal formation efficiency
$\epsilon=0.001$ as the default value in our model, which is in line with the
choice of other works (e.g. Drążkowska et al., 2016; Schoonenberg et al.,
2018; Stammler et al., 2019). Yet, a value of $10^{-3}$ seems unphysically
low. As we showed in Sect. 3.5.3, a higher planetesimal formation efficiency
implies that the ring will spawn planetesimal much faster and a planetesimal
formation efficiency exceeding $\epsilon\gtrsim 0.1$ will even consume the
ring as shown in JO21. The low $\epsilon$ we and other works adopt is partly
driven by these reason. A possible solution to this conundrum is to argue that
the low $\epsilon$ is a manifestation of clump-planetesimal interactions,
which tears apart pebble clumps before they gravitationally collapse. Then,
although initially planetesimals form on a settling timescale, the existence
of large population of planetesimal bodies would quickly suppress further
planetesimal formation, and reduce $\epsilon$. Future work to test this
hypothesis with hydrodynamical simulations is worthwhile.

In addition, recent streaming instability simulations show that the criterion
of starting SI is more complex than the simple dust-to-gas ratio $Z=1$ (e.g.,
Carrera et al., 2015; Yang et al., 2017; Li & Youdin, 2021). Specifically, for
pebbles exceeding ${\rm St}>0.01$, the typical required dust-to-gas is only
$0.5$. Conversely, for pebbles of Stokes number less than ${\rm St}<0.01$,
they find that $Z>2$ is needed for SI (see Li & Youdin, 2021, for more
discussion). In our simulation, we do not take this effect into consideration
but use the same criterion $Z=1$ for all St. Yet, since the difference in the
threshold criterion is at most a factor two, this is not expected to
qualitatively influence our results.

Our model neglects the gravitational feedback of the planet on the gas, i.e.,
gap opening for massive planets. Such a gap will both change the pebble
dynamics and slow down the migration of the planets. An upper mass above which
our model stops being valid is when the planet strongly perturbs the gas disc
(Eq. (38)), corresponding to $60\,M_{\oplus}$ at $74.3$ au in our standard
models. Such a mass is never reached in the runs that feature migration.
However, the gravitational feedback of planets on the gas disc matters when
there is no migration since the oligarch keeps growing in situ. We discuss the
implication of the following evolution in Sect. 4.3.

Finally, a potential important omission of our model is the lack of gas
accretion. Qualitatively, gas accretion is thought to be suppressed when the
planet is still accreting pebbles because gravitational energy released by the
accreted pebble balances the heat radiated away by the envelope, i.e., the
embryos are too hot to accrete gas (e.g., Rafikov, 2006; Lee & Chiang, 2015).
In addition, gas accretion can be suppressed by the the recycling of the
atmospheres (Moldenhauer et al., 2022) and by the enhanced opacity contributed
by the accreted pebbles (Ormel et al., 2021). On the other hand, the polluted
atmospheres are characterized by an increased mean molecular weight, which
would decrease the critical core mass (Hori & Ikoma, 2011; Venturini et al.,
2015; Brouwers & Ormel, 2020; Brouwers et al., 2021). In any case, when the
planet moves away from the ring region and pebble accretion subsides, gas
accretion is expected to accelerate and the planets formed in our simulations
may enter runaway gas accretion.

###  4.2 The “fine-tuned” optical depths of the ring

One striking feature of the ALMA rings in protoplanetary discs is that their
derived peak optical depths are relatively uniform. Dullemond et al. (2018)
find that the peak optical depths of rings in DSHARP are all around $0.4$, and
none of the rings are very optically thick. This result, which is seemingly
"fine-tuned", suggests there exists mechanism(s) regulating the upper limit of
the rings’ optical depths.

Several works have attempted to explain the fine-tuned optical depths. First,
the scattering opacity may play a significant role in the interpretation of
(sub)millimeter observations (e.g., Kataoka et al., 2015; Liu, 2019). Zhu et
al. (2019) pointed out that dust scattering can reduce the emission from a
highly optically thick ring and mimic it as a marginally optically thick ring.
Yet, the albedo of dust grains in discs is uncertain (Birnstiel et al., 2018).
Therefore, in order to answer the question whether or not scattering plays an
important role, we need to conduct longer-wavelength observations with
upcoming facilities like ngVLA (e.g., Ricci et al., 2018).

It has also been suggested that the fine-tuned optical depth is the
consequence of ongoing planetesimal formation (Stammler et al., 2019). As
planetesimals form when the mid-plane dust-to-gas ratio reaches unity, the
further growth of the ring will be limited, which regulates the optical depth,
as we studied in this work. Stammler et al. (2019) studied the case where the
ring is formed by a permanent pressure bump. Several hundreds of Earth masses
in planetesimals are formed inside the ring, while the further interaction
among these planetesimals and with the pebbles in the ring is neglected. We
complement the picture by taking pebble accretion, planetesimal dynamics, and
migration into account. For the scenario where the inward migration is
applicable and a steady state emerges, the saturated ring peak density then
explains the fine-tuned ring optical depth naturally. Specifically, for the
predicted ring with peak surface density $f_{\rm rs}\Sigma_{\rm tp}$ (Eq.
(33)), the optical depth reads

$\tau_{\nu}=\kappa_{\nu}f_{\rm rs}\Sigma_{\rm tp}$ (39)

where the opacity $\kappa_{\nu}$ is expressed in terms of particle size $s$,
internal density $\rho_{\rm s}$ and absorption coefficient $Q_{\nu,\rm abs}$:

$\kappa_{\nu}=\frac{\pi s^{2}}{m_{\rm s}}Q_{\nu,\rm
abs}=\frac{3}{4}\frac{Q_{\nu,\rm abs}}{s\rho_{\rm s}}.$ (40)

Since the Stokes number $St=\frac{\pi}{2}\frac{s\rho_{\rm s}}{\Sigma_{\rm g}}$
in the Epstein drag regime, we obtain

$\displaystyle\tau_{\nu}$
$\displaystyle=\frac{1}{2.9}\frac{3\pi}{8}\frac{f_{\rm rs}Q_{\nu,\rm
abs}\dot{M}_{\rm net}}{St^{\frac{5}{3}}q_{\rm
rp}^{\frac{2}{3}}\Omega_{K}r_{0}^{2}\Sigma_{\rm g}}$ (41)
$\displaystyle=0.4\left(\frac{Q_{\nu,\rm abs}}{0.05}\right)\left(\frac{f_{\rm
rs}}{1.8}\right)\left(\frac{f_{\rm
rp}}{1.4}\right)^{-\frac{2}{3}}\left(\frac{f_{\rm
mg}}{1}\right)^{\frac{1}{3}}\left(\frac{\dot{M}_{\rm net}}{50\,M_{\oplus}\,\rm
Myr^{-1}}\right)^{\frac{2}{3}}$ $\displaystyle\times\left(\frac{\Sigma_{\rm
g}}{3.6\,\rm g\,cm^{-2}}\right)^{-\frac{2}{3}}\left(\frac{w_{\rm
ring}}{3.38\,\rm
au}\right)^{-\frac{1}{3}}\left(\frac{St}{0.01}\right)^{-\frac{5}{3}}\left(\frac{h_{\rm
g}}{0.07}\right)^{-\frac{2}{3}}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-\frac{1}{3}}.$

In the Rayleigh regime ($\lambda=c/\nu\gg s$), the absorption coefficient
$Q_{\nu,\rm abs}\propto s\propto\mathrm{St}$ (e.g. van de Hulst, 1957), which
mitigates the St-dependence of the above expression. Specifically, we find
that $Q_{\nu,\rm abs}=0.05$ for $200\,\mu\mathrm{m}$-sized particles
(consistent with our default model) at wavelength of $1.25\rm mm$, based on
the Mie opacity calculation provided by Birnstiel et al. (2018). This value
agrees with the typical DSHARP ring optical depth.

###  4.3 In-situ formation vs factory assembly

Depending on the strength of migration, our model outputs can be classified
into two groups: in-situ formation where one planet remains at the ring
location and factory assembly where planetary embryos constantly form and
migrate inward. For the in-situ formation case, the biggest planet reached the
gap opening mass (Eq. (38)) as shown in Fig. 8. Such changes in the gas
profile will influence the aerodynamics of pebbles around the ring, which we
do not model with our simulation. By then, the ring should have already
vaporized, due to dominant pebble accretion and the repulsion of the pebbles
due to the emerging gas gap.

Self-destruction of rings by the formation of a massive planet inside has also
been found by other works. Lee et al. (2022) run 2D hydrodynamic simulations
investigating the growth of a single planet embryo inside a pebble ring
supported by Gaussian profile pressure bump as same as we assume. As they
found, the pebble ring will be ingested by the growing planet within ${\sim}1$
Myr (see Cummins et al. (2022) for a similar conclusion). During the clean-up
stage, asymmetric substructure may form, e.g., the arc substructure in the
continuum of Oph–IRS 48 (van der Marel et al., 2013; Cummins et al., 2022).
After the first ring disappears, a second ring could be generated exterior to
the gap-opening planet formed at the location (see Paardekooper et al., 2022,
and references therein). A new loop of in-situ planet formation may happen in
the newly formed ring.

However, Lee et al. (2022) reports the concern that the rapid evaporation of
the ring will shorten the observational window of the pebble ring, which may
conflict with the fact that rings are ubiquitous. Thus, they believe that the
formation of planets in ALMA rings must be either a rare or a slow process.
Accounting for migration, we show that a ring can be both long-lived and spawn
planets, which reconciles the high frequency of rings and the rapid growth of
planets.

###  4.4 Implication of the planet formation model

In this work, we investigated a scenario where planets form in distant rings.
Observational signatures of protoplanets are found in several discs. Here, we
compare our simulation outputs with several systems.

####  4.4.1 Planet-ring connection

Since we have modeled our setup on the B74 dust ring in AS 209 – a typical
DSHARP ring – it is justified to ask where its planets are? The hunt for
detecting planets in DSHARP discs has been attempted through direct imaging of
the planets (e.g., VLT/MUSE, Xie et al. (2020); VLT/NaCo, Jorquera et al.
(2021)) and resolving the potential circumplanetary disc (Andrews et al.,
2021). Yet, currently, no clear candidates have been found. Indirectly, by
studying the CO kinematic feature, it has been suggested that planets are
responsible for nine localized deviations from Keplerian rotation in the
kinematic map of eight DSHARP discs (e.g., Pinte et al., 2020). Except for one
candidate kink in HD 163296, all the other eight are located interior to a
pebble ring. It is natural to conceive, therefore, that pebble rings are
manifestations of pressure bumps produced by the planet. However, the non-
Keplerian motion associated with these pressure bumps is not detected in these
DSHARP rings (Pinte et al., 2020). Thus, the causal relationship "planet
first, then ring" may not always be appropriate. This leaves open the
possibility of "ring first, then planet", as proposed by our model.

Particularly, it has been proposed that an accreting planet is present in the
continuum ring of HD 100546. Casassus & Pérez (2019) detected a Doppler-flip
in the ${}^{12}\rm CO$ 2–1 line centroid, which Casassus et al. (2022)
interpreted in the context of a ${\sim}10\,M_{\oplus}$ accreting planet
launching an magnetocentrifugally driven outflow. This can cause a disc
eruption, and therefore could lead to a surface disturbance to the Keplerian
flow, explaining the flip and localized SO line emission. Such an embedded
planet would be in agreement with our model.

####  4.4.2 PDS 70

With two Jovian mass planets detected in the inner cavity, the transitional
disc PDS 70 is currently the only protoplanetary disc with unambiguous
evidence of harboring planets (Keppler et al., 2018; Haffert et al., 2019).
Bae et al. (2019) studied the architecture of the PDS 70 system by
hydrodynamical simulation and reproduces the orbital elements of the two
planets as well as the ALMA continuum. However, in their work, they initially
set the mass of PDS 70 b as $5\,M_{\rm Jup}$ and put it directly at 20 au
regardless of the formation history. Considering the 5.4 Myr-old age of PDS
70, the core formation of PDS 70 b, c is inefficient in a smooth disc. For
example, if the critical core mass of the PDS 70 protoplanet was
$10M_{\oplus}$ and the pebble mass flux amounted to $100\,M_{\oplus}\,\rm
Myr^{-1}$, it requires ${\sim}10$ Myr for a Pluto-mass embryo at 20 au to grow
in situ by pebble accretion in the 3D accretion regime. The growth timescale
is even longer for PDS 70 c.

In addition, the orbits of PDS 70 b, c are consistent with a 2:1 mean motion
resonance configuration (e.g., Wang et al., 2021a). The appearance of this
first-order mean motion resonance hints that they may form in sequence in the
more outer region and experienced convergent migration (Terquem & Papaloizou,
2007), which will require an even longer pebble accretion timescale.
Therefore, mechanisms to accelerate the growth of the protoplanet are
required.

Our model of planet assembly in rings has the potential to explain the
formation history of PDS 70 b and c, i.e., the two planetary cores are formed
from an outer first-generation ring, which may still appear as the continuum
peak at ${\sim}74$ au. Since planet b formed earlier than planet c, it also
started gas accretion earlier and transform into the slower type-II migration
before planet c. Together with the fact that PDS 70 c is more massive than PDS
70 b and thus experienced stronger type-I migration torque, it is reasonable
to assume that they will undergo convergent migration, and then form a common
gap (e.g., Dong & Fung, 2017).

####  4.4.3 HR 8799 and other planet systems in debris disc

Apart from planet-forming discs, another line of evidence for planet formation
in rings comes from the planets detected in debris discs. By statistically
comparing the dust mass budget in protoplanetary discs and debris disc, Michel
et al. (2021) propose that debris disc are the leftovers of planetesimal-
forming rings. This idea is supported by recent theoretical works (e.g., Jiang
& Ormel, 2021; Najita et al., 2022).

Among debris discs, a number of sources are also confirmed with mature planet
system inside the debris disc belt, e.g., the four giant planet system HR
8799. With years of follow-up by direct imaging, the orbital configuration of
the four giant planets in HR 8799 has been shown to be consistent with a
1:2:4:8 mean motion resonance chain (Wang et al., 2018), which, again, hints
at the convergent migration of planets. In addition, a clear dust ring is
resolved in ALMA (sub)millimeter continuum, whose radial intensity profile
peaks around 150 au (e.g., Booth et al., 2016; Faramaz et al., 2021). By
modeling the radiation transfer in ALMA band 6 and band 7, and reproducing
both the Hershel far-infrared map and ALMA continuum, Geiler et al. (2019)
infer a total mass of ${\sim}100\,M_{\oplus}$ of planetesimals up to
${\sim}100\rm\,km$ by collisional simulations, which qualitatively matches the
mass budget in the remaining planetesimals of our model (Sect. 3.7).

In the context of our model, the four giant planet cores form from a
planetesimal-forming pebble ring around 150 au. The cores migrate inward one
by one, evolve into gas giants, and form the Laplace resonance chain. After
the gas disc has dissipated, the leftover planetesimals belt from the birth
ring stays in situ and now constitutes a debris disc belt. By a collisional
cascade, the planetesimals belt will generate a moderate amount of dust which
might be observed in far-IR or (sub)millimeter wavelength (e.g., Wyatt & Dent,
2002). Detailed modeling of the formation history of PDS 70 and HR 8799 with
proper treatment of gas accretion and type-II migration will be studied in a
future work.

In addition, observationally, resolved debris disc belts are typically much
wider than rings in protoplanetary discs. If rings in protoplanetary discs and
debris disc are related, it leaves open the question of how to reconcile their
different widths. Miller et al. (2021) proposed that this can be explained by
a planetesimal-forming ring migrating inward together with the planet. In
JO21, we also showed that the clumpy ring can migrate outward since the pebble
supply is upstream. In this work, with fixed ring location , we show that the
dynamical interaction between planet and planetesimal can also widen the
debris disc ring width (see Fig. 15). Debris disc formed in this way should
have a sharper inner edge and flatter outer edge (see Sect. 3.7). Such an
asymmetry profile is found in several debris disc systems in the ALMA
continuum, e.g., HR 8799 (Faramaz et al., 2021) and q1 Eri (Lovell et al.,
2021).

##  5 Conclusions

Rings are the locations where pebbles accumulate to elevate the dust-to-gas
ratio, to create conditions conducive to planetesimal formation and further
planet assembly. We considered two designs for the rings. One is the clumpy
ring (CR) supported by robust aerodynamic backreaction, and the other is the
pressure bump (BP) ring which features a Gaussian profile pressure bump. Both
designs can reproduce the observations by ALMA in the (sub)millimeter
continuum. We have followed the growth of planetesimals formed in these
massive pebble ring by N-body simulations.Type-I migration plays a decisive
role in the evolution of the ring. Our main findings are as follows:

  1. 1.

Pebble rings are ideal locations for planet formation at large orbital
distances as the restraining conditions that normally limit pebble accretion
no longer apply. The high concentration of pebbles and the low headwind
velocity inside the ring render pebble accretion very efficient.

  2. 2.

Consequently, our results hardly depend on the initial planetesimal mass. We
test that planetesimals whose mass is down to Ceres immediately enjoy
effective pebble accretion, while the conclusion can be extended to lower
initial mass generally.

  3. 3.

At a distance of 74.3 au and an external mass flux of $100\,M_{\oplus}\,\rm
Myr^{-1}$, a ${\sim}10\,M_{\oplus}$ planet grows from planetesimals formed
inside the ring within 0.5 Myr.

  4. 4.

Planet formation in rings exhibit similar characteristics, independent of the
mechanism supporting the ring (aerodynamically-supported clumpy ring model or
pressure bump-induced). For the CR model the mass required to trigger
planetesimal formation is relatively modest due to its small scaleheight,
while for PB rings this depends on the pebble scaleheight
($\delta_{z}/\mathrm{St}$) and could be high.

  5. 5.

For discs in which inward type-I migration operates across rings, planets
migrate away from the pebble ring after reaching a typical mass of
${\sim}20\,M_{\oplus}$ (Eq. (32)). If external conditions will remain similar,
the process will repeat and the ring will spawn a new planet.

  6. 6.

The ring’s mass is a balance between the consumption by pebble accretion and
the external supply. The mass of the ring is typically ${\sim}30\,M_{\oplus}$
(Eq. (34)) at $74.3\,\rm au$, consistent with the mass of annular rings
observed in ALMA.

  7. 7.

Planet-forming rings at larger orbital radii are more massive. Yet, the mass
of the protoplanet only weakly depends on the locations of the ring. Equation
(35) describes the relation between the observed ring mass and the typical
planetary core mass that can form inside the ring.

  8. 8.

On the other hand, for disc in which type-I migration is inefficient or rings
where planets are trapped, a single big planet will emerge and consume the
ring by pebble accretion. In that situation, the window for it to appear as a
ring in ALMA imagery could be shorter than (${\lesssim}1\,\mathrm{Myr}$).

  9. 9.

After the planet migrates away from the ring, its growth by pebble accretion
drops rapidly. In the CR design, planets can still grow from the downstream
pebble mass flux and it is conceivable that these planets will trigger runaway
gas accretion after exiting the ring region.

  10. 10.

The remnant planetesimal belts are candidates for the cold debris discs.
Scattered or shepherded by the formed planet, the radial profile of the
planetesimal belt shows an asymmetry with a steep rise closer at its interior
side and a more gradual decline towards the exterior regions – a feature that
may be observable.

  11. 11.

We obtain that the (sub))millimeter optical depth of the planetesimal-forming
ring is around ${\sim}0.4$ with minimal dependence on disc and pebble
properties (Eq. (41)). This conclusion is in line with observations that discs
are marginally optically thick (cf. Dullemond et al., 2018; Stammler et al.,
2019).

The roles of ring and planet recall the classical “chicken or egg” dilemma.
According to our model, the massive rings in the protoplanetary disc perform
the role of egg-laying chicken. Our model provides an alternative view on the
connection between planets and rings than the mainstream model, in which
planets induce the formation of the rings.

This leaves the open question of how the first-generation rings form. Due to
the great variety in discs properties, it is likely that rings observed in
ALMA originate from different mechanisms, e.g., sintering (Okuzumi et al.,
2016), secular gravitational instability (Tominaga et al., 2019), coagulation
front of pebbles (Ohashi et al., 2021), anisotropic infall (Kuznetsova et al.,
2022). Such a primordial ring can be the birthplace of the first-generation
planet, which can later sculpt the disc to possibly generate new rings.

The model where planets form in a pebble ring and then migrate inward has the
potential to explain certain systems involving multiple giant planets in young
protoplanetary discs and debris discs. Application of this model to the
formation history of realistic planetary architecture will be conducted in a
follow-up work.

## Acknowledgements

H.J. would like to thank Beibei Liu and Douglas N. C. Lin for useful
discussions. This work has used Matplotlib (Hunter, 2007), Numpy (Harris et
al., 2020), Scipy (Virtanen et al., 2020) software packages.

## Data Availability

The data underlying this article will be shared on reasonable requests to the
corresponding author.

## References

  * Abod et al. (2019) Abod C. P., Simon J. B., Li R., Armitage P. J., Youdin A. N., Kretke K. A., 2019, ApJ, 883, 192
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Andrews et al. (2021) Andrews S. M., et al., 2021, ApJ, 916, 51
  * Bae et al. (2019) Bae J., et al., 2019, ApJ, 884, L41
  * Baehr & Zhu (2021) Baehr H., Zhu Z., 2021, ApJ, 909, 136
  * Bai & Stone (2014) Bai X.-N., Stone J. M., 2014, ApJ, 796, 31
  * Baruteau et al. (2014) Baruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 667 (arXiv:1312.4293), doi:10.2458/azu_uapress_9780816531240-ch029
  * Batygin & Laughlin (2015) Batygin K., Laughlin G., 2015, Proceedings of the National Academy of Science, 112, 4214
  * Birnstiel et al. (2010) Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79
  * Birnstiel et al. (2018) Birnstiel T., et al., 2018, ApJ, 869, L45
  * Bitsch et al. (2018) Bitsch B., Morbidelli A., Johansen A., Lega E., Lambrechts M., Crida A., 2018, A&A, 612, A30
  * Bitsch et al. (2019) Bitsch B., Izidoro A., Johansen A., Raymond S. N., Morbidelli A., Lambrechts M., Jacobson S. A., 2019, A&A, 623, A88
  * Booth et al. (2016) Booth M., et al., 2016, MNRAS, 460, L10
  * Boss (1997) Boss A. P., 1997, Science, 276, 1836
  * Brouwers & Ormel (2020) Brouwers M. G., Ormel C. W., 2020, A&A, 634, A15
  * Brouwers et al. (2021) Brouwers M. G., Ormel C. W., Bonsor A., Vazan A., 2021, A&A, 653, A103
  * Carrera et al. (2015) Carrera D., Johansen A., Davies M. B., 2015, A&A, 579, A43
  * Carrera et al. (2021) Carrera D., Simon J. B., Li R., Kretke K. A., Klahr H., 2021, AJ, 161, 96
  * Carrera et al. (2022) Carrera D., Thomas A. J., Simon J. B., Small M. A., Kretke K. A., Klahr H., 2022, ApJ, 927, 52
  * Casassus & Pérez (2019) Casassus S., Pérez S., 2019, ApJ, 883, L41
  * Casassus et al. (2022) Casassus S., Cárcamo M., Hales A., Weber P., Dent B., 2022, ApJ, 933, L4
  * Chambers (2021) Chambers J., 2021, ApJ, 914, 102
  * Chiang (2008) Chiang E., 2008, ApJ, 675, 1549
  * Cieza et al. (2019) Cieza L. A., et al., 2019, MNRAS, 482, 698
  * Cresswell & Nelson (2008) Cresswell P., Nelson R. P., 2008, A&A, 482, 677
  * Cummins et al. (2022) Cummins D. P., Owen J. E., Booth R. A., 2022, arXiv e-prints, p. arXiv:2206.13866
  * Deng et al. (2021) Deng H., Mayer L., Helled R., 2021, Nature Astronomy, 5, 440
  * Dong & Fung (2017) Dong R., Fung J., 2017, ApJ, 835, 146
  * Dong et al. (2015) Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5
  * Drazkowska et al. (2022) Drazkowska J., et al., 2022, arXiv e-prints, p. arXiv:2203.09759
  * Drążkowska et al. (2016) Drążkowska J., Alibert Y., Moore B., 2016, A&A, 594, A105
  * Drążkowska et al. (2021) Drążkowska J., Stammler S. M., Birnstiel T., 2021, A&A, 647, A15
  * Dubrulle et al. (1995) Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237
  * Dullemond et al. (2018) Dullemond C. P., et al., 2018, ApJ, 869, L46
  * Faramaz et al. (2021) Faramaz V., et al., 2021, AJ, 161, 271
  * Flock et al. (2015) Flock M., Ruge J. P., Dzyurkevich N., Henning T., Klahr H., Wolf S., 2015, A&A, 574, A68
  * Geiler et al. (2019) Geiler F., Krivov A. V., Booth M., Löhne T., 2019, MNRAS, 483, 332
  * Goldreich & Tremaine (1979) Goldreich P., Tremaine S., 1979, ApJ, 233, 857
  * Guilera et al. (2020) Guilera O. M., Sándor Z., Ronco M. P., Venturini J., Miller Bertolami M. M., 2020, A&A, 642, A140
  * Haffert et al. (2019) Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., Girard J. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hayashi (1981) Hayashi C., 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Hori & Ikoma (2011) Hori Y., Ikoma M., 2011, MNRAS, 416, 1419
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Ida & Makino (1993) Ida S., Makino J., 1993, Icarus, 106, 210
  * Jang et al. (2022) Jang H., Liu B., Johansen A., 2022, arXiv e-prints, p. arXiv:2205.10103
  * Jiang & Ormel (2021) Jiang H., Ormel C. W., 2021, MNRAS, 505, 1162
  * Jiang et al. (2022) Jiang H., Zhu W., Ormel C. W., 2022, ApJ, 924, L31
  * Johansen & Lambrechts (2017) Johansen A., Lambrechts M., 2017, Annual Review of Earth and Planetary Sciences, 45, 359
  * Johansen et al. (2006) Johansen A., Henning T., Klahr H., 2006, ApJ, 643, 1219
  * Johansen et al. (2015) Johansen A., Mac Low M.-M., Lacerda P., Bizzarro M., 2015, Science Advances, 1, 1500109
  * Jorquera et al. (2021) Jorquera S., et al., 2021, AJ, 161, 146
  * Kanagawa et al. (2018) Kanagawa K. D., Tanaka H., Szuszkiewicz E., 2018, ApJ, 861, 140
  * Kataoka et al. (2015) Kataoka A., et al., 2015, ApJ, 809, 78
  * Keppler et al. (2018) Keppler M., et al., 2018, A&A, 617, A44
  * Kim et al. (2020) Kim S., et al., 2020, ApJ, 888, 72
  * Klahr & Schreiber (2021) Klahr H., Schreiber A., 2021, ApJ, 911, 9
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Kretke & Lin (2007) Kretke K. A., Lin D. N. C., 2007, ApJ, 664, L55
  * Kuznetsova et al. (2022) Kuznetsova A., Bae J., Hartmann L., Low M.-M. M., 2022, ApJ, 928, 92
  * Lambrechts & Johansen (2012) Lambrechts M., Johansen A., 2012, A&A, 544, A32
  * Lambrechts & Johansen (2014) Lambrechts M., Johansen A., 2014, A&A, 572, A107
  * Law et al. (2021) Law C. J., et al., 2021, ApJS, 257, 3
  * Lee & Chiang (2015) Lee E. J., Chiang E., 2015, ApJ, 811, 41
  * Lee et al. (2022) Lee E. J., Fuentes J. R., Hopkins P. F., 2022, arXiv e-prints, p. arXiv:2206.01219
  * Levison et al. (2012) Levison H. F., Duncan M. J., Thommes E., 2012, AJ, 144, 119
  * Li & Youdin (2021) Li R., Youdin A. N., 2021, ApJ, 919, 107
  * Li et al. (2019) Li R., Youdin A. N., Simon J. B., 2019, ApJ, 885, 69
  * Lin et al. (2018) Lin J. W., Lee E. J., Chiang E., 2018, MNRAS, 480, 4338
  * Liu (2019) Liu H. B., 2019, ApJ, 877, L22
  * Liu & Ji (2020) Liu B., Ji J., 2020, Research in Astronomy and Astrophysics, 20, 164
  * Liu & Ormel (2018) Liu B., Ormel C. W., 2018, A&A, 615, A138
  * Liu et al. (2019) Liu B., Ormel C. W., Johansen A., 2019, A&A, 624, A114
  * Liu et al. (2020) Liu B., Lambrechts M., Johansen A., Pascucci I., Henning T., 2020, A&A, 638, A88
  * Long et al. (2018) Long F., et al., 2018, ApJ, 869, 17
  * Lovell et al. (2021) Lovell J. B., et al., 2021, MNRAS, 506, 1978
  * Masset et al. (2006) Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, ApJ, 642, 478
  * Michel et al. (2021) Michel A., van der Marel N., Matthews B. C., 2021, ApJ, 921, 72
  * Miller et al. (2021) Miller E., Marino S., Stammler S. M., Pinilla P., Lenz C., Birnstiel T., Henning T., 2021, MNRAS, 508, 5638
  * Miotello et al. (2022) Miotello A., Kamp I., Birnstiel T., Cleeves L. I., Kataoka A., 2022, arXiv e-prints, p. arXiv:2203.09818
  * Moldenhauer et al. (2022) Moldenhauer T. W., Kuiper R., Kley W., Ormel C. W., 2022, A&A, 661, A142
  * Morbidelli (2020) Morbidelli A., 2020, A&A, 638, A1
  * Najita et al. (2022) Najita J. R., Kenyon S. J., Bromley B. C., 2022, ApJ, 925, 45
  * Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
  * Öberg et al. (2021) Öberg K. I., et al., 2021, ApJS, 257, 1
  * Ohashi et al. (2021) Ohashi S., et al., 2021, ApJ, 907, 80
  * Okuzumi et al. (2016) Okuzumi S., Momose M., Sirono S.-i., Kobayashi H., Tanaka H., 2016, ApJ, 821, 82
  * Ormel (2017) Ormel C. W., 2017, The Emerging Paradigm of Pebble Accretion. p. 197, doi:10.1007/978-3-319-60609-5_7
  * Ormel & Klahr (2010) Ormel C. W., Klahr H. H., 2010, A&A, 520, A43
  * Ormel & Liu (2018) Ormel C. W., Liu B., 2018, A&A, 615, A178
  * Ormel et al. (2021) Ormel C. W., Vazan A., Brouwers M. G., 2021, A&A, 647, A175
  * Owen & Kollmeier (2019) Owen J. E., Kollmeier J. A., 2019, MNRAS, 487, 3702
  * Paardekooper et al. (2022) Paardekooper S.-J., Dong R., Duffell P., Fung J., Masset F. S., Ogilvie G., Tanaka H., 2022, arXiv e-prints, p. arXiv:2203.09595
  * Papaloizou & Larwood (2000) Papaloizou J. C. B., Larwood J. D., 2000, MNRAS, 315, 823
  * Pinilla et al. (2021) Pinilla P., Lenz C. T., Stammler S. M., 2021, A&A, 645, A70
  * Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Rafikov (2006) Rafikov R. R., 2006, ApJ, 648, 666
  * Rein & Liu (2012) Rein H., Liu S. F., 2012, A&A, 537, A128
  * Ricci et al. (2018) Ricci L., Liu S.-F., Isella A., Li H., 2018, ApJ, 853, 110
  * Rosotti et al. (2020) Rosotti G. P., Teague R., Dullemond C., Booth R. A., Clarke C. J., 2020, MNRAS, 495, 173
  * Rosotti et al. (2021) Rosotti G. P., Ilee J. D., Facchini S., Tazzari M., Booth R. A., Clarke C., Kama M., 2021, MNRAS, 501, 3427
  * Schäfer et al. (2017) Schäfer U., Yang C.-C., Johansen A., 2017, A&A, 597, A69
  * Schoonenberg & Ormel (2017) Schoonenberg D., Ormel C. W., 2017, A&A, 602, A21
  * Schoonenberg et al. (2018) Schoonenberg D., Ormel C. W., Krijt S., 2018, A&A, 620, A134
  * Schoonenberg et al. (2019) Schoonenberg D., Liu B., Ormel C. W., Dorn C., 2019, A&A, 627, A149
  * Segura-Cox et al. (2020) Segura-Cox D. M., et al., 2020, Nature, 586, 228
  * Sekiya (1998) Sekiya M., 1998, Icarus, 133, 298
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Shariff & Cuzzi (2015) Shariff K., Cuzzi J. N., 2015, ApJ, 805, 42
  * Sheehan & Eisner (2017) Sheehan P. D., Eisner J. A., 2017, ApJ, 840, L12
  * Sheehan & Eisner (2018) Sheehan P. D., Eisner J. A., 2018, ApJ, 857, 18
  * Sheehan et al. (2020) Sheehan P. D., Tobin J. J., Federman S., Megeath S. T., Looney L. W., 2020, ApJ, 902, 141
  * Shibata & Ikoma (2019) Shibata S., Ikoma M., 2019, MNRAS, 487, 4510
  * Simon et al. (2016) Simon J. B., Armitage P. J., Li R., Youdin A. N., 2016, ApJ, 822, 55
  * Stammler et al. (2019) Stammler S. M., Drążkowska J., Birnstiel T., Klahr H., Dullemond C. P., Andrews S. M., 2019, ApJ, 884, L5
  * Tamayo et al. (2020) Tamayo D., Rein H., Shi P., Hernandez D. M., 2020, MNRAS, 491, 2885
  * Tanaka et al. (2002) Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
  * Teague et al. (2018) Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ, 860, L12
  * Terquem & Papaloizou (2007) Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110
  * Tominaga et al. (2019) Tominaga R. T., Takahashi S. Z., Inutsuka S.-i., 2019, ApJ, 881, 53
  * Venturini et al. (2015) Venturini J., Alibert Y., Benz W., Ikoma M., 2015, A&A, 576, A114
  * Villenave et al. (2020) Villenave M., et al., 2020, A&A, 642, A164
  * Villenave et al. (2022) Villenave M., et al., 2022, arXiv e-prints, p. arXiv:2204.00640
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Visser & Ormel (2016) Visser R. G., Ormel C. W., 2016, A&A, 586, A66
  * Voelkel et al. (2021) Voelkel O., Deienno R., Kretke K., Klahr H., 2021, A&A, 645, A132
  * Wagner et al. (2019) Wagner K., Apai D., Kratter K. M., 2019, ApJ, 877, 46
  * Wang et al. (2018) Wang J. J., et al., 2018, AJ, 156, 192
  * Wang et al. (2021a) Wang J. J., et al., 2021a, AJ, 161, 148
  * Wang et al. (2021b) Wang S., Kanagawa K. D., Suto Y., 2021b, ApJ, 923, 165
  * Weidenschilling (1977a) Weidenschilling S. J., 1977a, Ap&SS, 51, 153
  * Weidenschilling (1977b) Weidenschilling S. J., 1977b, MNRAS, 180, 57
  * Wyatt & Dent (2002) Wyatt M. C., Dent W. R. F., 2002, MNRAS, 334, 589
  * Xie et al. (2020) Xie C., Haffert S. Y., de Boer J., Kenworthy M. A., Brinchmann J., Girard J., Snellen I. A. G., Keller C. U., 2020, A&A, 644, A149
  * Xu & Bai (2022) Xu Z., Bai X.-N., 2022, arXiv e-prints, p. arXiv:2207.08858
  * Xu et al. (2017) Xu Z., Bai X.-N., Murray-Clay R. A., 2017, ApJ, 847, 52
  * Yang et al. (2017) Yang C. C., Johansen A., Carrera D., 2017, A&A, 606, A80
  * Youdin & Lithwick (2007) Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588
  * Youdin & Shu (2002) Youdin A. N., Shu F. H., 2002, ApJ, 580, 494