# Language Model Metrics and Procrustes Analysis
for Improved Vector Transformation of NLP Embeddings

Thomas Conley
University of Colorado
Colorado Springs
1420 Austin Bluffs Pkwy
Colorado Springs, CO, USA
<EMAIL_ADDRESS>
Jugal Kalita
University of Colorado
Colorado Springs
1420 Austin Bluffs Pkwy
Colorado Springs, CO, USA
<EMAIL_ADDRESS>
###### Abstract

Artificial Neural networks are mathematical models at their core. This truism
presents some fundamental difficulty when networks are tasked with Natural
Language Processing. A key problem lies in measuring the similarity or
distance among vectors in NLP embedding space, since the mathematical concept
of distance does not always agree with the linguistic concept. We suggest that
the best way to measure linguistic distance among vectors is by employing the
Language Model (LM) that created them. We introduce Language Model Distance
(LMD) for measuring accuracy of vector transformations based on the
Distributional Hypothesis (_LMD_Accuracy_). We show the efficacy of this
metric by applying it to a simple neural network learning the Procrustes
algorithm for bilingual word mapping.

##  1 Introduction

The Distributional Hypothesis Firth (1961) inspired the development of
embeddings that capture the meaning of language based on how words co-occur
with each other Mikolov et al. (2013a). Natural Language Processing relies
heavily on these high dimensional vectors to represent words, phrases,
sentences or documents, in a form that can be processed by deep neural
networks which were originally designed for tasks related to computer vision.
Input embeddings are transformed by network layers into output vectors which
represent solutions to many NLP tasks Ruder et al. (2019).

In order to learn these transformations, a network must be able to calculate
the difference between predicted vectors and actual word vectors. This
distance calculation is a crucial part of measuring loss, and performing back-
propagation. These core functions of neural networks have primarily relied on
mathematical processes without regard to linguistic principles. We demonstrate
that NLP embedding transformation is better measured using linguistic
similarity functions rooted in knowledge of languages rather than concepts
such as Euclidean or angular distance, which assumes vectors to be “physical”
objects.

###  1.1 Procrustes Analysis

Matrix transformation of vector spaces has been accomplished using Generalized
Procrustes Analysis (GPA), ever since a computationally viable solution was
devised Gower (1975). In particular, GPA has been used to great effect in geo-
spatial shape manipulation Duta (2015); Crosilla et al. (2019) and qualitative
data analysis Maurício et al. (2016).

Shapes are represented by a series of landmark points in 2 or 3 dimensions.
And in survey research, qualitative opinion data is represented by a Likert
scale Likert (1932), occupying low dimensional space. In both cases, the
vector spaces must be realigned and resized for meaningful comparison.
Although these fields seem to differ, they use data structures that share
characteristics with Natural Language Processing.

The orthogonal Procrustes algorithm produces an optimal transformation matrix
$R$ for mapping one vector space to another and appears to be useful in
converting vector spaces for NLP tasks such as bilingual word mapping
Kementchedjhieva et al. (2018).

###  1.2 Procrustes Analysis for NLP tasks

Can a neural network learn to do Procrustes transformation? The answer, yes,
should be non-controversial, since every neural network performs tensor
transformation of input to output. However, tasks which require nuanced
understanding of the meaning of words, such bilingual word mapping, are
particularly difficult. Although there is some success when massive amounts of
text are available for training, the problem is more acute when resources for
learning are scarce, as in machine translation of under resourced languages.

The difficulty with vector transformations in NLP is based on the nature of
the data. NLP transformations by neural networks use distance measurements
designed to work in $L_{p}\;space$. This implies numerical data. We show that
such calculations of distance and accuracy are not as effective as
measurements based on language models.

###  1.3 Image data and language data

We consider image data as raw data with physical dimensionality where, each
dimension in a vector can be considered similar in measurement and meaning. As
such, this data occupies $L_{p}\;space$; where vectors can be added together
or multiplied by scalars without loss of their inherent meaning. For example,
a vector representing a pixel is measured the same way, and has the same
meaning, regardless of where it is in the image.

Thus, distance measurement among image vectors can use $L_{p}\;norm$ or
trigonometric calculations such as cosine distance. One specific kind of
euclidean distance measurement is called _Procrustes Distance_ and is the
basis of Procrustes Analysis Crosilla et al. (2019).

In NLP, distance measurement is less meaningful when it is based on Euclidean
axioms rather than linguistic principles. Distance is the basis of error
calculations and back-propagation, and so, the ability to calculate the
derivative of these functions is essential for classic stochastic gradient
descent (SGD) which is employed by neural networks today. Although there has
been some research in non-differentiable losses Engilberge et al. (2019) the
mathematical requirements for these functions are not always suitable for NLP.

As opposed to raw data, feature data consists of vectors in which each
dimension may have disparate meaning and measurement. Feature data does not
exist in $L_{p}\;space$, and therefore measures of distance that rely on
$L_{p}\;norm$ or trigonometric calculations may not be meaningful. We consider
NLP embeddings to be feature data, although they share some characteristics
with raw data.

##  2 Language Models and Data

As in raw data, NLP vectors dimensions typically share values that are treated
similarly and are thus undifferentiated in a sense. This seems to contradict
the assertion that each NLP embedding dimension has a specific unique meaning
like feature data. Instead, the meaning of a dimension is more like
probability, representing how often a word is used with a particular meaning,
rather than the actual meaning of the word.

Vectors with dimensions that differ in meaning, as in NLP embeddings, cannot
be used with typical spatial measurements such as $L_{p}\;norm$ and
$cosine\;distance$. We contend that NLP vector distance can best be measured
by the language models which represent the vectors. Therefore, we seek to
replace mathematic calculations with predictions from language models. We
simply rely on the language model itself to provide a distance measurement for
our custom metric.

In this research, we use the Word2Vec model Mikolov et al. (2013b) to produce
a custom bilingual word mapping dataset. This dataset, combined with the
GenSim model of keyed vectors Řehůřek and Sojka (2010), provides a
distributional distance measurement based on word movers distance Kusner et
al. (2015).

Our neural network is a simple Multilayer Perceptron (MLP) which accepts
Spanish word vectors as input and predicts English word vectors. This simple
model was chosen because it is analogous to any layer found in innumerable,
more complex, neural networks. Showing improved efficacy in this model should
demonstrate improvement in any NLP task.

Figure 1: Illustration of the Distributional Hypothesis and Language Model
Distance. The accuracy of predicted vectors $\hat{p_{i}}$ and $\hat{p_{j}}$,
is based on membership in the set of $k=2$ or $k=3$ neighbors.

##  3 Language Model Distance

An exact measurement of equality is not possible for high-dimensional NLP
embeddings. Embeddings of several hundred dimensions, and one-hot encoded
vectors on the order of tens of thousands of dimensions, are particularly
difficult to measure.

$LMD(\hat{p},t,m,k)=\begin{cases}True,&\text{if}\;t\in m.set(\hat{p},k)\\\
False,&\text{otherwise}\end{cases}$ (1)

Instead, we suggest that the true measure of NLP vector distance is best
provided by the model which defines the vectors. We present a family of
metrics, Language Model Distance (LMD), which calculates distance and equality
among NLP vectors by using the language model itself. $LMD$ is defined as in
Equation 1 where the distance between predicted vector $\hat{p}$ and known
truth vector $t$, is provided by model $m$, given neighbor threshold $k$.

The distance measure is binary because it is based on set inclusion, and not
physical or Euclidean distance. Thus, $LMD$ can be used as a measure of
accuracy, and records a true positive when $t$ is within the neighborhood of
the predicted vector ($t\in m.set(\hat{p},k)$).

###  3.1 Measuring Accuracy with Language Model Distance

Figure 1 illustrates the distributional hypothesis by showing a simple
clustering along 2 non-numeric dimensions. The circles represent neighborhoods
$m.set(\hat{p},k=2)$ and $m.set(\hat{p},k=3)$. Note that the predicted vectors
($\hat{p}$) have no words directly associated with them, because no exact
match is possible for floating point numeric vectors.

Thus we say that $LMD\\_Accuracy(k)$ measures a positive result when truth
vector ($t$) is within the $k$ sized neighborhood of the predicted vector
($t\in m.set(\hat{p},k)$). For example, $LMD\\_Accuracy(3)$ measures the
percentage of times that the true word answer was among the top 3 closest
predicted words.

Distributional distance functions can be used in neural network metrics, loss,
or activation functions, or used directly in similarity computation. However,
inserting external language models into neural networks can be difficult as
these networks are firmly rooted in mathematics which is not compatible with
linguistic processes.

We solve these difficulties by defining a simple class shown in Figure 2. By
including the language model as a static member of the class, methods of the
class may be used as network internal functions with access to external
language models.

Figure 2: Implementation of Distributional Accuracy based on Language Model
Distance. A static language model (line 2) allows linguistic functionality to
be included in purely mathematical models.

##  4 Learning Orthogonal Procrustes Analysis

The Orthogonal Procrustes Algorithm is a process for finding the optimal
mapping of one set of vectors to another. Typically, the vectors represent
points in 2 or 3 dimensional space, for image processing, or they represent
qualitative data measured in few dimensions Maurício et al. (2016). After
resizing and repositioning of vectors, an optimal rotation matrix $R$ is
produced by a method similar to singular value decomposition.

This classic approach to vector transformation has been explored as a solution
for some NLP tasks Sen et al. (2019); Kim et al. (2019). Therefore we ask: Can
a neural network be trained to perform the same optimal transformation for NLP
embeddings which occupy a much higher dimensional space?

Our task is to train a simple MLP to learn the optimal mapping $R$, between
two disparate vector spaces representing a bilingual dataset. We measure the
success of this task using $LMD$ as the basis for accuracy as in Figure 2 and
Equation 1.

We create two separate language models from a parallel corpus of European
Parliament translations, the so called EuroParl dataset Koehn (2005). We use
the Word2Vec model in continuous bag-of-words (CBOW) mode Mikolov et al.
(2013a) to build two separate distributions. By using a bilingual corpus, and
training language models separately, we ensure that the models share a common
domain, but the vector spaces remain separate. For training, we then map word
vectors from one distribution to the other, using a set of 1000 most common
words pairs, obtained from from a language learning
website111http://www.englishnspanish.com.

###  4.1 Results

Our results show that Orthogonal Procrustes Analysis can be learned for
multilingual mapping of word vectors. Furthermore, Figure 3 demonstrates that
LMD is effective as a basis for measuring the accuracy of this task.

Figure 3: Results of Learning Orthogonal Procrustes Analysis showing a better
measure of exact matches with _LMD_Accuracy_ than with _cosine similarity_.

Figure 3 indicates that _LMD_Accuracy_ is better at measuring similarity in
NLP embeddings than _cosine similarity_. In this plot, _LMD_Accuracy(1)_
indicates that the model exactly predicted the correct word in the output
language. When _LMD_Accuracy(1)_ is near 100% the value of _cosine similarity_
should be near $1$ which would indicate an exact match. The fact that _cosine
similarity_ cannot measure this exact match shows a weakness in this purely
mathematical measurement compared with our language model-based measurement.

##  5 Learning General Procrustes Analysis

To further test, we try to learn General Procrustes Analysis; a much harder
task because it requires the network to generalize.

We have just shown that a simple neural network can learn to transform
vectors. This is non-controversial since all neural networks perform this task
at every layer. However, not all networks are able to generalize. Using the
same network configuration as before, we now evaluate embeddings that we have
not seen in training, as is common. This is equivalent to learning the
_Generalized Procrustes Algorithm_.

Figure 4: Results of Learning General Procrustes Analysis showing a comparable
measure of exact matches between _LMD_Accuracy_ and _cosine similarity_ , when
generalization is required

.

Results in Figure 4 show that _LMD_Accuracy_ is more like cosine distance when
generalization is required. Note that we use _LMD_Accuracy_ only for metrics.
This model uses $cosine\;similarity$ for error calculation and back-
propagation. We conclude that such $L_{p}\;norm$ measurements can only drive
generalization as far as they are able to measure accuracy.

The local variation in _LMD_Accuracy_ , evident in Figure 4, may be
significant as it may make determining the derivative of the function
difficult. The derivative of _LMD_Accuracy_ must be worked out before it can
be incorporated into a loss function and be used in back-propagation. The
overall shape of the curve, despite irregularities is encouraging as the slope
may be computed using ordinary least squares in a calculation of rolling
regression, or by other numerical methods.

##  6 Conclusion

We suggest that language model metrics described here may be incorporated
directly into activation and loss functions, and may be used as an error
measurement for back-propagation. We suggest this basic enhancement would
improve the Generalized Procrustes Algorithm and other NLP processing in
general. This is left for future work.

## References

  * Crosilla et al. (2019) Fabio Crosilla, Alberto Beinat, Andrea Fusiello, Eleonora Maset, and Domenico Visintini. 2019.  Orthogonal procrustes analysis.  In _Advanced Procrustes Analysis Models in Photogrammetric Computer Vision_ , pages 7–28. Springer.
  * Duta (2015) Nicolae Duta. 2015.  _Procrustes Shape Distance_ , pages 1278–1279. Springer US, Boston, MA.
  * Engilberge et al. (2019) Martin Engilberge, Louis Chevallier, Patrick Pérez, and Matthieu Cord. 2019\.  Sodeep: a sorting deep net to learn ranking loss surrogates.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 10792–10801.
  * Firth (1961) John Rupert Firth. 1961.  _Papers in Linguistics 1934-1951: Repr_.  Oxford University Press.
  * Gower (1975) John C Gower. 1975.  Generalized procrustes analysis.  _Psychometrika_ , 40(1):33–51.
  * Kementchedjhieva et al. (2018) Yova Radoslavova Kementchedjhieva, Sebastian Ruder, Ryan Cotterell, and Anders Søgaard. 2018.  Generalizing procrustes analysis for better bilingual dictionary induction.  In _22nd Conference on Computational Natural Language Learning (CoNLL 2018)_ , pages 211–220. Association for Computational Linguistics.
  * Kim et al. (2019) Yunsu Kim, Petre Petrov, Pavel Petrushkov, Shahram Khadivi, and Hermann Ney. 2019\.  Pivot-based transfer learning for neural machine translation between non-english languages.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 865–875.
  * Koehn (2005) Philipp Koehn. 2005.  Europarl: A parallel corpus for statistical machine translation.  In _MT summit_ , volume 5, pages 79–86. Citeseer.
  * Kusner et al. (2015) Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger. 2015.  From word embeddings to document distances.  volume 37 of _Proceedings of Machine Learning Research_ , pages 957–966, Lille, France. PMLR.
  * Likert (1932) Rensis Likert. 1932.  A technique for the measurement of attitudes.  _Archives of psychology_.
  * Maurício et al. (2016) Angélica Maurício, A.B. Palazzo, Valeria Caselato, and Helena Bolini. 2016.  Generalized procrustes analysis and external preference map used to consumer drivers of diet gluten free product.  _Food and Nutrition Sciences_ , 07:711–723.
  * Mikolov et al. (2013a) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013a.  Distributed representations of words and phrases and their compositionality.  In _Advances in neural information processing systems_ , pages 3111–3119.
  * Mikolov et al. (2013b) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013b.  Distributed representations of words and phrases and their compositionality.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 26_ , pages 3111–3119. Curran Associates, Inc.
  * Řehůřek and Sojka (2010) Radim Řehůřek and Petr Sojka. 2010.  Software Framework for Topic Modelling with Large Corpora.  In _Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks_ , pages 45–50, Valletta, Malta. ELRA.  http://is.muni.cz/publication/884893/en.
  * Ruder et al. (2019) Sebastian Ruder, Ivan Vulić, and Anders Søgaard. 2019.  A survey of cross-lingual word embedding models.  _Journal of Artificial Intelligence Research_ , 65:569–631.
  * Sen et al. (2019) Sukanta Sen, Kamal Kumar Gupta, Asif Ekbal, and Pushpak Bhattacharyya. 2019.  Multilingual unsupervised nmt using shared encoder and language-specific decoders.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3083–3089.


11institutetext: Department of Astrophysics/IMAPP, Radboud University, P.O.
Box 9010, NL-6500 GL Nijmegen, The Netherlands 22institutetext: Lund
Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43, 22100 Lund, Sweden 33institutetext: Institute for Astronomy, School of
Physics, Zhejiang University, Hangzhou 310027, China 44institutetext: Center
for Star and Planet Formation, GLOBE Institute, University of Copenhagen,
Øster Voldgade 5-7, 1350 Copenhagen, Denmark
44email: hyerin.jang@astro.ru.nl,bbliu<EMAIL_ADDRESS>
# Growth after the streaming instability:
The radial distance dependence of the planetary growth

Hyerin Jang 1122 Beibei Liu 3⋆3⋆ Anders Johansen 4422

Streaming instability is hypothesized to be triggered at particular
protoplanetary disk locations where the volume density of the solid particles
is enriched comparable to that of the gas. A ring of planetesimals thus forms
when this condition is fulfilled locally. These planetesimals collide with
each other and accrete inward drifting pebbles from the outer disk to further
increase masses. We investigate the growth of the planetesimals that form in a
ring-belt at various disk radii. Their initial mass distributions are
calculated based on the formula summarized from the streaming instability
simulations. We simulate the subsequent dynamical evolution of the
planetesimals with a protoplanetary disk model based either on the minimum
mass solar nebula (MMSN) or on the Toomre stability criterion. For the MMSN
model, both pebble accretion and planetesimal accretion are efficient at a
close-in orbit of $0.3$ AU, resulting in the emergence of several super-Earth
mass planets after $1$ Myr. For comparison, only the most massive
planetesimals undergo substantial mass growth when they are born at $r{=}3$
AU, while the planetesimals at $r{=}30$ AU experience little or no growth. On
the other hand, in the denser Toomre disk, the most massive forming planets
can reach Earth mass at $t{=}1$ Myr and reach a mass between that of Neptune
and that of Saturn within $3$ Myr at $30$ AU and $100$ AU. Both the pebble and
planetesimal accretion rate decrease with disk radial distance. Nevertheless,
planetesimal accretion is less pronounced than pebble accretion at more
distant disk regions. Taken together, the planets acquire higher masses when
the disk has a higher gas density, a higher pebble flux, and/or a lower Stokes
number of pebbles.

###### Key Words.:

numerical – planets and satellites: formation

##  1 Introduction

####

The growth of millimeter- to centimeter-sized bodies (named pebbles hearafter)
into kilometer- or larger sized planetesimals faces a challenge in the
classical planet formation theory (see Johansen et al. (2014), Liu & Ji (2020)
and Drazkowska et al. (2022) for reviews). Their growth is mainly limited by
bouncing, fragmentation (Güttler et al., 2010; Zsom et al., 2010), and fast
radial drift (Weidenschilling, 1977). The streaming instability mechanism
provides a promising solution to overcome the above barriers, when the pebbles
self-concentrate together and collapse into planetesimals due to the
collective effect of gravity (Youdin & Goodman, 2005; Johansen & Youdin,
2007).

In order to initiate the streaming instability, the volume density of pebbles
needs to be enriched comparable to that of the gas, $\rho_{\rm
peb}{\simeq}\rho_{\rm g}$ (Youdin & Goodman, 2005). This criterion is also
frequently written in terms of the disk metallicity with a surface density
ratio, $Z{=}\Sigma_{\rm peb}/\Sigma_{\rm g}$. Literature numerical simulations
find that a super-solar metallicity (${\gtrsim}2{-}5\%$) is generally required
for the onset of the streaming instability (Carrera et al., 2015; Yang et al.,
2017; Li & Youdin, 2021a). Specifically, Yang et al. (2017) found particles of
$\tau_{s}=10^{-2}$ and $10^{-3}$ are concentrated by the streaming instability
when the disk metallicities are $1-2\%$ and $3-4\%$, respectively. Li & Youdin
(2021b) suggested the metallicity could be lower than 1% with $\tau_{s}>0.01$
in their particle clumping simulations. Higher $Z$ might be needed if
considering a size distribution of pebbles (Zhu & Yang, 2021). Overall, the
threshold metallicity is difficult to satisfy globally in a canonical solar-
like protoplanetary disk.

Various hypotheses have been proposed so that the enhancement of the pebble
surface density can be fulfilled locally at peculiar protoplanetary disk
locations. For instance, the ice-line location of different chemical species,
such as water, silicate, and carbon monoxide (Ros & Johansen, 2013;
Schoonenberg & Ormel, 2017; Drażkowska & Alibert, 2017; Hyodo et al., 2019)
can act as planetesimal nurseries, since large pebbles can grow there by
condensation. Pressure bumps are potential sites for triggering the streaming
instability where drifting pebbles get accumulated and the pressure gradient
varies.111 The streaming instability does not formally apply at the pressure
maximum point, since the pressure gradient is zero there. Another potential
location is the dead-zone inner edge. Since the ionization level and density
are different between the turbulent active region and quiescent dead-zone
region, their boundary can naturally trap drifting pebbles and induce
streaming instability (Chatterjee & Tan, 2014; Ueda et al., 2019).
Furthermore, pebbles might become clumped at the vortex generated by
hydrodynamical instabilities (Surville et al., 2016) or the spiral arms in
self-gravitating disks (Gibbons et al., 2012; Elbakyan et al., 2020). Taken
together, the triggering site can broadly span from the most inner sub-AU disk
region to the most outer disk regions extending to a few tens or hundreds of
AUs.

When the streaming instability is triggered at a single-site disk location, a
ring of planetesimals forms rapidly (Liu et al., 2019b). Their birth masses
follow a top-heavy distribution, which can be approximately described by a
power-law plus exponential decay function (Johansen et al., 2015; Schäfer et
al., 2017; Abod et al., 2019). The planetesimals have a characteristic mass of
${\sim}10^{-6}\ M_{\oplus}$ ($100$ km in size) when they form at the asteroid
belt region in the minimum mass solar nebula (MMSN) model (Johansen et al.,
2015). The mass increases with gas surface density, disk metallicity, and
radial distance (Liu et al., 2020).

After their formation, these planetesimals can further increase their masses
by colliding with other planetesimals in the birth ring, as well as by
accreting inwardly drifting pebbles from the outer disk. As expected,
planetesimal accretion (Kokubo & Ida, 1996) and pebble accretion (Ormel &
Klahr, 2010; Lambrechts & Johansen, 2012) operate together to contribute to
the subsequent mass growth. In this paper, we assume a streaming instability-
induced planetesimal formation at specific disk locations. We investigate the
subsequent mass growth of these planetesimals in a ring-belt structure by a
combination of planetesimal and pebble accretion. The growth at the
conventional water-ice-line location ($r_{\text{ice}}{=}2.7$ AU) has been
investigated by Liu et al. (2019b). We follow their numerical construction and
mainly focus on how the growth patterns vary as a function of radial distance.

We introduce the model setup in Section 2. In Section 3, we show different
dynamical evolution of the planetesimals at $0.3$ AU, $3$ AU, and $30$ AU in
the MMSN model. We also explore how the model parameters (planet migration,
pebble flux, Stokes number of pebbles) affect the final forming planets in
Section 4. We adopt the upper limit of a stable disk under the Toomre’s
criterion and investigate whether giant planets can form early at the distant
orbits of $30$ and $100$ AU in Section 5. The discussion and final conclusions
are summarized in Sections 6 and 7.

##  2 Method

####

We introduce the disk model in Section 2.1 and calculate the $r-$dependent
planetesimal mass distributions in Section 2.2. The N-body numerical setups
are described in Section 2.3.

###  2.1 Disk model

This study focuses on the planetesimal growth around the solar-mass central
star. We adopt the MMSN disk model (Hayashi, 1981). The disk temperature, gas
surface density, gas disk aspect ratio, and pressure gradient prefactor are
expressed as

$T=T_{0}\left(\dfrac{r}{1\text{AU}}\right)^{-1/2},$ (1)
$\Sigma_{\text{gas}}=\Sigma_{\text{gas},0}\left(\dfrac{r}{1\text{AU}}\right)^{-3/2},$
(2) $h_{\text{gas}}=h_{\text{gas},0}\left(\dfrac{r}{1\text{AU}}\right)^{1/4},$
(3)

and

$\eta{=}\eta_{0}\left(\dfrac{r}{1\text{AU}}\right)^{1/2},$ (4)

where $r$ is the disk radial distance, $\Sigma_{\text{gas},0}=1700$ g cm-2,
$T_{0}=280$ K, $h_{\text{gas},0}=0.033$, and $\eta_{0}=1.8\times 10^{-3}$ are
the corresponding values at $1$ AU.

The detailed modeling of the dust growth and radial drift is not the focus of
our study. Instead, we assume that the dust particles have already grown to
pebble-sized bodies with a unified dimensionless stopping time (or termed
Stokes number), and these pebbles sweep the disk from outside-in with a
constant flux. We chose the Stokes number of pebbles $\tau_{\rm s}{=}0.1$ and
the pebble flux of $\dot{\text{M}}_{\text{peb}}{=}100\ \rm M_{\oplus}/Myr$ in
the fiducial run. The above simplified treatment is a proper approximation
from the more advanced dust coagulation models (Birnstiel et al., 2012;
Lambrechts & Johansen, 2014; Drażkowska et al., 2021).

###  2.2 Mass distribution of planetesimals in a ring

####

Since the disk properties are correlated with radial distance, the forming
planetesimals are also expected to feature an $r$-dependence. We provide their
mass distribution calculations as follows:

Schäfer et al. (2017) showed that the initial mass function of planetesimals
that form by streaming instability can be fitted by a power law plus an
exponential decay function. The number fraction of planetesimals with
$M_{\text{p}}>M$ is given by

$\dfrac{N_{>}(M)}{N_{\text{tot}}}=\left(\dfrac{M}{M_{\text{min}}}\right)^{-p}\text{exp}\left[\left(\dfrac{M_{\text{min}}}{M_{\text{plt}}}\right)^{q}-\left(\dfrac{M}{M_{\text{plt}}}\right)^{q}\right],$
(5)

where $p{=}0.6$ and $q{=}0.35$. Equation (5) features two parameters: the
minimum mass and the characteristic mass of planetesimals ($M_{\text{min}}$
and $M_{\rm plt}$). We find that the shape of the distribution and especially
the planetesimal population that dominates the total mass are barely affected
by $M_{\text{min}}$, as long as $M_{\text{min}}{\ll}M_{\rm plt}$ (see Appendix
A). We thus set $M_{\rm min}{=}M_{\rm plt}/20$. Using a given characteristic
mass $M_{\rm plt}$ and total mass of planetesimals in the belt $M_{\rm tot}$,
we can obtain the planetesimal mass distribution from Equation (5). As we show
later, both $M_{\rm plt}$ and $M_{\rm tot}$ depend on the forming locations of
the planetesimals.

We adopt the characteristic planetesimal mass from Liu et al. (2020) based on
the extrapolation of literature on streaming instability simulations,

$\dfrac{M_{\text{plt}}}{M_{\oplus}}=5\times
10^{-5}\left(\dfrac{Z}{0.02}\right)^{a}\left(\dfrac{\gamma}{\pi^{-1}}\right)^{a+1}\left(\dfrac{h_{\text{gas}}}{0.05}\right)^{3+b},$
(6)

where the fitting coefficients are $a{=}0.5$ and $b{=}0$. The self-gravity
parameter $\gamma$ is given by

$\gamma=\dfrac{4\pi G\rho_{\text{gas}}}{\Omega_{\text{K}}^{2}},$ (7)

where the Keplerian orbital frequency
$\Omega_{\text{K}}{\equiv}\sqrt{GM_{\star}/r^{3}}$, $G$ is the gravitational
constant, $M_{\star}$ is the stellar mass, and $\rho_{\text{gas}}$ is the
midplane gas density

$\rho_{\text{gas}}=\dfrac{\Sigma_{\text{gas}}}{\sqrt{2\pi}H_{\text{gas}}},$
(8)

where $H_{\text{gas}}=h_{\rm gas}r$ is the gas disk scale height.

In Equation (6), the local metallicity $Z$ is estimated from the streaming
instability triggering condition where
$\rho_{\text{peb}}{\simeq}\rho_{\text{gas}}$ (Youdin & Goodman, 2005). This
can be equivalently expressed as
$\Sigma_{\text{peb}}/H_{\text{peb}}\simeq\Sigma_{\text{gas}}/H_{\text{gas}}$.
Hence, the threshold metallicity for the onset of streaming instability can be
written as $Z_{\rm
SI}{=}\Sigma_{\text{peb}}/\Sigma_{\text{gas}}{=}H_{\text{peb}}/H_{\text{gas}}$.
The disk scale height ratio between pebbles and gas is given by (Youdin &
Lithwick, 2007)

$\dfrac{H_{\text{peb}}}{H_{\text{gas}}}=\sqrt{\dfrac{\alpha_{t}}{\alpha_{t}+\tau_{s}}},$
(9)

where $\alpha_{t}$ is the gas turbulent diffusion coefficient. We assume that
the dust and gas have the same diffusivity for tightly coupled solid
particles. Besides that, the coefficient of turbulent gas diffusivity
approximates to the turbulent viscosity parameter for magnetorotational
instability (MRI) disks (Johansen & Klahr, 2005; Zhu et al., 2015; Yang et
al., 2018). In this paper, we do not distinguish these two and adopt one
unified parameter as $\alpha_{t}$. We set $\alpha_{t}{=}10^{-4}$ based on the
observational diagnostic of disk turbulence (Pinte et al., 2016; Flaherty et
al., 2017). Therefore, the threshold metallicity in the fiducial setup is
$0.0316$, higher than the nominal solid-to-gas column density ratio of $0.01$.

Combining Eqs. (6), (7) and (8), we can rewrite the characteristic mass as

$\dfrac{M_{\text{plt}}}{M_{\oplus}}=1.2\times
10^{-6}\left(\dfrac{Z}{0.02}\right)^{1/2}\left(\dfrac{\Sigma_{\rm
gas,0}}{1700}\right)^{3/2}\left(\dfrac{h_{\rm
gas,0}}{0.033}\right)^{3/2}\left(\dfrac{r}{\text{3 AU}}\right)^{9/8}.$ (10)

The total mass of planetesimals in the ring belt can be estimated as

$M_{\rm tot}=\epsilon M_{\rm tot,peb}=2\pi\epsilon r\Delta
r\Sigma_{\text{peb}}=2\pi\epsilon Zr\Delta r\Sigma_{\text{gas}},$ (11)

where $\epsilon$ is the pebble-to-planetesimal conversion efficiency. As
explored and discussed in Abod et al. (2019), $\epsilon$ is at approximately
ten percent and should vary with the radial pressure gradient, the Stokes
number of pebbles, and the disk metallicity. However, due to limited numerical
streaming instability simulation explorations on $\epsilon$, we followed the
same procedure as Liu et al. (2019b) and we set this efficiency at a fixed
value of $10\%$ for the sake of simplicity. The width of the planetesimal belt
can be estimated from the dense dust filament in streaming instability
simulations, which is expressed as $\Delta r=\eta r$ (Yang & Johansen, 2016;
Li et al., 2018).

In sum, the planetesimal population can be calculated from Eqs. (5), (10) and
(11), where the detailed numerical approach is illustrated in Appendix B.
Table 1 summarizes the generated planetesimal population with its
characteristic, minimum, and maximum masses in the MMSN model, while Figure 1
illustrates the discrete mass distributions at $r{=}0.3$ AU, $3$ AU, and $30$
AU, respectively.

We note that streaming instability simulations produce numerous small
planetesimals. From Equation (5), the total number of planetesimals are
dominated by low-mass bodies. Nevertheless, these bodies only contribute in a
very limited way to the subsequent mass growth (Liu et al., 2019b). Conducting
simulations with such a huge number of small planetesimals is numerically
impractical. For the sake of computational feasibility, we simulate the
massive bodies (dark bars in Figure 1), which corresponds to two-thirds of the
total masses in the belt.

We also performed a simulation with $80\%$ of the total masses in the belt to
investigate the effect of neglected masses on the simulations. Because of the
difference in number and total mass of planetesimals, the final mass of the
most massive planet is not exactly identical to the fiducial case. However,
the overall features of the growth pattern (discussed in Section 3) are very
similar. It validates our approach in the main paper without exhaustively
simulating all small planetesimals.

Table 1: Properties of planetesimals generated by streaming instability.

Description |  $r$ (AU) | $\tau_{s}$ | $N_{\text{tot}}$ | $N_{\text{sim}}$ |  $M_{\text{plt}}$ (M⊕) |  $M_{\text{min}}$ (M⊕) |  $M_{\text{max}}$ (M⊕) |  $M_{\text{tot}}$ (M⊕)
---|---|---|---|---|---|---|---|---
Fiducial | 0.3 | 0.1 | 9943 | 924 | $1.3\times 10^{-7}$ | $6.3\times 10^{-9}$ | $1.3\times 10^{-5}$ | $6.7\times 10^{-4}$
Fiducial | 3 | 0.1 | 7470 | 651 | $1.7\times 10^{-6}$ | $8.4\times 10^{-8}$ | $1.5\times 10^{-4}$ | $6.7\times 10^{-3}$
Fiducial | 30 | 0.1 | 5618 | 514 | $2.3\times 10^{-5}$ | $1.1\times 10^{-6}$ | $1.7\times 10^{-3}$ | $6.7\times 10^{-2}$
Low Stokes number | 0.3 | 0.03 | 13431 | 1224 | $1.7\times 10^{-7}$ | $8.5\times 10^{-9}$ | $1.9\times 10^{-5}$ | $1.2\times 10^{-3}$

222$N_{\text{tot}}$ is the total number of planetesimals in the belt,
$N_{\text{sim}}$ is the number of selected planetesimals for numerical
simulations, $M_{\text{plt}}$ is the characteristic mass, and $M_{\text{tot}}$
is the total mass of the planetesimals. The Stokes number ($\tau_{s}$) of
pebbles is $0.1$ in the fiducial run. For comparison, we set $\tau_{s}{=}0.03$
for the lower Stokes number run (Section 4).

Figure 1: Initial mass distributions of planetesimals generated by streaming
instability in a ring belt of the MMSN model at $r{=}0.3$ AU, $3$ AU, and $30$
AU, respectively. The darker bars represent the population of planetesimals
with $M{>}M_{\rm plt}$, which account for two-thirds of the whole population
in mass. We only consider the above subset of planetesimals (darker) in the
following numerical simulations. At 3 and 30 AU, the last darker bar is higher
than the previous bar because the decimal places have been neglected when
counting the number of planetesimals. If the number of planetesimals at the
mass just less than the maximum mass is higher than one but less than two, it
has been treated as one. Thus, the total mass of planetesimals in that bar
becomes lower than their last bars.

###  2.3 Numerical setup

####

The MERCURY N-body code (Chambers, 1999) is adopted, where we choose the
Burlirsch-Stoer integration, $3$ days for the initial timestep, and $10^{-12}$
for the integration accuracy. The code was modified to take into account of
the effects of pebble accretion for planets/planetesimals, the gas drag for
planetesimals, gravitational tidal torques for planets, and
eccentricities/inclinations damping effects (see details in Liu et al.
(2019b)).

Planetesimal accretion happens when the planetesimals collide with each other.
For the N-body part, we simply assume that the collisions are ideally
inelastic to conserve the angular momentum. Fragmentation is not taken into
account, since the eccentricities of planetesimals are relatively small and
the impact velocity is generally lower than the escape velocity.

The pebble accretion prescription is adopted from Liu & Ormel (2018) and Ormel
& Liu (2018), where they account for the stochastic motion of pebbles due to
disk turbulence. We also use the pebble isolation mass formula from Bitsch et
al. (2018), where

$M_{\rm iso}=25f_{\rm fit}\rm M_{\oplus},$ (12)

where

$f_{\rm
fit}=\left[\dfrac{H/r}{0.05}\right]^{3}\left[0.34\left(\dfrac{\log(0.001)}{\log(\alpha_{\rm
t})}\right)^{4}+0.66\right]\left[1-\dfrac{\dfrac{\partial\ln P}{\partial\ln
r}+2.5}{6}\right],$ (13)

$P$ is the gas pressure and $\alpha_{\rm t}$ stands for the turbulent
viscosity parameter.

Once the planet reaches the pebble isolation mass, its mass growth via pebble
accretion stops. Moreover, we simply assume that the other planets located
interior to this planet also terminate their pebble accretion. This is because
a pressure bump is produced outside of the orbit of the planet that reaches
the pebble isolation mass. The upcoming pebbles stop their radial drift
exterior to this pressure bump and, as a result, all the inner bodies cannot
further accrete pebbles. However, we note that a fraction of pebble fragments
at the abovementioned pressure bump can still pass the gap and replenish the
inner disk region (Liu et al., 2022). In addition, the pressure bump can also
be the formation site of the next generation planetesimals, since the inward
drifting pebbles are continuously accumulated (Eriksson et al., 2021, 2022).
Such effects are not considered here. Additionally, we also do not take into
account of the accretion of gas in this work.

Planetesimals experience an aerodynamic gas drag, and its deceleration is
expressed as

$\boldsymbol{a_{\text{drag}}}=-\left(\dfrac{3C_{\text{D}}\rho_{\text{gas}}}{8R_{\text{p}}\rho_{\bullet}}\right)v_{\text{rel}}\boldsymbol{v_{\text{rel}}},$
(14)

where $C_{\text{D}}{=}0.5$ is the drag coefficient, $R_{\text{p}}$ is the
radius of the planetesimal, $\rho_{\bullet}{=}1.5\rm\ gcm^{-3}$ is the
internal density, and $\boldsymbol{v_{\text{rel}}}$ is the relative velocity
between the planetesimal and gas
$\boldsymbol{v_{\text{rel}}}{=}\boldsymbol{v_{\text{plt}}}-\boldsymbol{v_{\text{gas}}}$,
where $\boldsymbol{v_{\text{gas}}}$ is the gas rotation velocity
$\boldsymbol{v_{\text{K}}}(1-\eta)$, and $\boldsymbol{v_{\text{K}}}$ is the
Keplerian velocity.

When a planetesimal grows massive enough to excite density waves, these waves
exert torques onto the planet and result in its orbital decay (Kley & Nelson,
2012). In our simulations, the highest mass that planets can reach is
generally ${\lesssim}10\ M_{\oplus}$ (see Sections 3 and 4), and therefore
only the type I torque is considered in our study. The corresponding
accelerations are given by

$\boldsymbol{a_{m}}=-\boldsymbol{v}/t_{m},\
\boldsymbol{a_{e}}=-2\dfrac{(\boldsymbol{v\cdot
r})\boldsymbol{r}}{r^{2}t_{e}},\
\boldsymbol{a_{i}}=\dfrac{\boldsymbol{v_{z}}}{t_{i}},$ (15)

where $\boldsymbol{v}$ is the velocity vector of a migrating body, and
$t_{m}$, $t_{e}$, and $t_{i}$ are the type I migration, eccentricity, and
inclination damping timescales (see Equation (11-13) of Cresswell & Nelson
(2008)). The potential thermal radiative and dust feedback onto the type I
migration torques is neglected in this work(Jiménez & Masset, 2017; Benítez-
Llambay & Pessah, 2018).

In order to minimize the number of free parameters and investigate the in situ
growth, we set the semimajor axis damping to zero as $\boldsymbol{a_{m}}{=}0,$
while the eccentricity and inclination damping are included in Section 3 for
the fiducial run. Other setups of simulations incorporating type I migration,
and different values of pebble flux and Stokes number, are presented in
Section 4.

##  3 Results

####

In this section, we investigate the mass growth of planetesimals at three
different disk locations. The model parameters are adopted from the fiducial
run in Table 2. Figure 2, 3 and 4 show the mass and semi-major axis evolution
of the forming planetesimals at $r{=}0.3$ AU, $3$ AU, and $30$ AU for $1$ Myr.
The gradual increase in mass is due to pebble accretion, while the sudden
bumps are due to collisions among planetesimals (marked as red dots). These
figures show planetesimals that do not undergo significant mass growth, that
experience substantial mass growth but fail to reach the pebble isolation
mass, and that reach the pebble isolation mass within the simulated time,
respectively.

Figure 2:  (a) Mass growth and (b) semimajor axis evolution of planetesimals
in a ring located at $r{=}0.3$ AU. Each line represents one planetesimal and
red dots indicate planetesimal-planetesimal collisions. Gray lines represent
all planetesimals in the simulation. The yellow and blue lines correspond to
the bodies that finally reach the pebble isolation mass, and those that fail
to reach the pebble isolation mass but grow above $10^{-2}\ M_{\oplus}$,
respectively. The green horizontal line represents the time that large
planetesimals start rapid pebble accretion. The mass growth is rapid at close-
in orbits for both pebble and planetesimal accretion. Figure 3: (a) Mass
growth and (b) semimajor axis evolution of planetesimals in a ring located at
$r{=}3$ AU. The same colored plots and lines are used as for Figure 2. The
dominant mass growth channel is pebble accretion. Figure 4: (a) Mass growth
and (b) semimajor axis evolution of planetesimals in a ring located at
$r{=}30$ AU. The notations in the graph are the same as Figure 2. The growth
of planetesimal is insignificant at such distant orbits in the MMSN disk
model.

We see a fast growth of planetesimals when they form at $r{=}0.3$ AU in Figure
2a. Three large planetesimals grow their masses effectively by pebble
accretion, and one of them reaches the pebble isolation mass of $2.6\
M_{\oplus}$ at $0.31$ AU (yellow line). The other two planetesimals stop their
mass growth with $0.37\ M_{\oplus}$ at $0.25$ AU and $1.5\ M_{\oplus}$ at
$0.27$ AU (blue lines). We find in Figure 2b that the three massive growing
embryos scatter each other and shield smaller planetesimals exterior and
interior to their orbits after a few $10^{5}$ years.

The growth is slower when the planetesimals are born at $r{=}3$ AU than at
$0.3$ AU. In Figure 3a, after a few collisions, two initially massive
planetesimals (${\sim}10^{-4}\ M_{\oplus}$) start rapid pebble accretion and
finally reach around $0.07\ M_{\oplus}$ (blue lines). The pebble isolation
mass is $14\ M_{\oplus}$ at the corresponding radial distance. These two most
massive planets stay around their initial semimajor axes, while other smaller
planetesimals are continuously scattered by them into the exotic orbits (see
Figure 3b). Even though the planets do not reach the pebble isolation mass
within 1 Myr, they reach the pebble isolation mass of $14\ M_{\oplus}$ within
2-3 Myr.

Although the generated bodies are more massive at $r{=}30$ AU compared to
those at closer-in orbits (Table 1), we find that the mass growth is
substantially suppressed when the planetesimals are born at a more distant
disk location. The collisions are very rare in Figure 4a and the planetesimals
of $10^{-3}\ M_{\oplus}$ do not experience any mass increase. The width of the
initial belt expands slowly, less than $50\%$ over $1$ Myr (see Figure 4b). We
will provide the analytical explanations for such $r$-dependent growth
patterns in the following subsections.

###  3.1 Planetesimal accretion

####

The planetesimal collision rate can be expressed as

$P_{c}=n\sigma\delta v,$ (16)

where $n$ is the number density of planetesimals expressed as $n=\Sigma_{\rm
plt}/MH$, where $M$ and $\Sigma_{\rm plt}$ are the mass and surface density of
planetesimals, $\sigma$ is a collisional cross section area, and $\delta v$ is
relative velocity. Since the vertical height of planetesimals $H\simeq\delta
v/\Omega_{\rm K}$, we can rewrite Equation 16 as

$P_{\rm c}=\Sigma_{\rm plt}\Omega_{\rm K}\sigma/M=\frac{3\Sigma_{\rm
plt}\Omega_{\rm K}}{4\rho_{\bullet}R_{\rm plt}}\left(1+\frac{v_{\rm
esc}^{2}}{\delta v^{2}}\right),$ (17)

where $v_{\rm esc}$ and $R_{\rm plt}$ are the escape velocity and size of the
planetesimals. We note that the mass of the planetesimal is lower, and the
planetesimal surface density and the Keplerian orbital frequency are higher,
at a smaller $r$. Even though the cross section $\sigma$ slightly decreases
due to mutual excited orbits, the overall collision rate is much higher when
$r$ is smaller. If we assume equal-sized planetesimals and the geometric
accretion regime, the planetesimal collision timescale $\tau_{\rm
c}{=}1/P_{\rm c}\propto r^{3}$. This is why, in the simulations, more frequent
planetesimal-planetesimal collisions occur, and their mass growth via
planetesimal accretion is significant at close-in orbits. In total, there are
$902$, $111$ and $22$ collisions for the cases at $r{=}0.3$, $3,$ and $30$ AU,
respectively.

###  3.2 Pebble accretion

####

The efficiency of pebble accretion correlates with the mass of the
planetesimal/planet and the size of accreted pebbles (Ormel & Liu, 2018).
Pebble accretion only becomes efficient when the planetesimal has a relatively
high mass; otherwise the gas-aided settling effect is insignificant (Visser &
Ormel, 2016; Liu & Ji, 2020). The onset mass for the efficient pebble
accretion can be expressed as (Liu & Ji, 2020)

$M_{\text{onset}}=\tau_{s}\eta^{3}M_{\star}\simeq 10^{-3}\left(\frac{\tau_{\rm
s}}{0.1}\right)\left(\frac{\eta_{0}}{1.8\times
10^{-3}}\right)^{3}\left(\dfrac{r}{\text{3 AU}}\right)^{3/2}M_{\oplus}.$ (18)

As shown in Equation (18), the onset mass scales super-linearly with $r$
through the $\eta$ dependence. A more massive planetesimal is required to
initiate the rapid pebble accretion at a larger orbital distance. For
instance, $M_{\text{onset}}\simeq 3\times 10^{-5}\ M_{\oplus}$ at $r{=}0.3$
AU, almost corresponding to the highest-mass planetesimal in the initially
formed belt (Table 1). At $r{=}3$ AU, $M_{\rm onset}\simeq 10^{-3}\
M_{\oplus}$, a factor of $10$ higher than the highest-mass planetesimal in the
belt. The mass difference is even larger at $r{=}30$ AU, where $M_{\rm
onset}\simeq 3\times 10^{-2}\ M_{\oplus}$ and $M_{\rm max}{=}2\times 10^{-3}\
M_{\oplus}$. As a direct consequence, the planetesimals that form at a larger
orbital distance need to grow more massive in the first place, in order to
trigger the efficient pebble accretion.

As discussed in Section 3.1, the early mass increase is aided by planetesimal
accretion, which also features a radial distance dependence. Significant
mutual collisions among planetesimals occur faster at a smaller $r$.
Therefore, the largest planetesimal reaches $M_{\rm onset}$ at an earlier time
and consequently starts effective pebble accretion afterward. In Figure 2a,
the slope of the first growing planetesimal (yellow line) steepens at around
$10^{4}$ years, while in Figure 3a, the growing planetesimals steepen after
$6\times 10^{5}$ years. We do not see such a transition between two mass
growth regimes in Figure 4a.

The timescales of pebble accretion in 2D/3D regimes can be expressed as (Liu
et al., 2019b)

$\begin{split}\tau_{\text{PA,2D}}\simeq&6.5\times
10^{4}\left(\dfrac{M}{0.05\text{M}_{\oplus}}\right)^{1/3}\left(\dfrac{\tau_{s}}{0.1}\right)^{1/3}\left(\dfrac{\eta}{1.8\times
10^{-3}}\right)\\\
&\left(\dfrac{\dot{M}_{\text{peb}}}{100\text{M}_{\oplus}/\text{Myr}}\right)^{-1}\text{yr}\end{split}$
(19)

and

$\tau_{\text{PA,3D}}\simeq 1.5\times
10^{4}\left(\dfrac{h_{\text{peb}}}{1\times
10^{-3}}\right)\left(\dfrac{\eta}{1.8\times
10^{-3}}\right)\left(\dfrac{\dot{M}_{\text{peb}}}{100\text{M}_{\oplus}/\text{Myr}}\right)^{-1}\text{yr},$
(20)

where $h_{\rm peb}{=}H_{\rm peb}/r$ is the pebble disk aspect ratio (Equation
(9)). The above equations are valid when the planet mass is much higher than
$M_{\rm onset}$. Since the pressure gradient prefactor $\eta$, the
planetesimal mass, and the aspect ratio of pebble disk all increase with $r$,
the adopted timescales are $\tau_{\text{PA,2D}}\propto r^{1/2}$ and
$\tau_{\text{PA,3D}}\propto r^{3/4}$. Therefore, after reaching $M_{\rm
onset}$, the planetesimals also take more time to grow their masses at a
larger orbital distance through pebble accretion.

Accounting for the above factors, we find that 1) although the initial birth
mass and the total mass in the planetesimal belt are lower, the subsequent
planetesimal growth pattern is still faster when they form at closer-in
orbits, and 2) planetesimal accretion has a steeper $r$-dependence than pebble
accretion, making it even more inefficient at larger orbital distances.

Table 2: Simulation setups for the parameter study.

Description |  | type I
---
migration
| Pebble flux
---
[M⊕ Myr-1]
$\tau_{\rm s}$
1\. fiducial | No | 100 | 0.1
2\. migration | Yes | 100 | 0.1
3\. high pebble flux | No | 200 | 0.1
4\. low Stokes number | No | 100 | 0.03
333Only one model parameter is varied for different cases compared to the
fiducial setup. The semimajor axis decay due to type I migration is
implemented in the migration case, the pebble flux is increased to $200\
M_{\oplus}\ \rm Myr^{-1}$ in the high pebble flux case, and the Stokes number
is reduced to $0.03$ in the low Stokes number case.

##  4 Parameter study

####

In Section 3 we only conduct simulations with the fiducial setup. Here, we
vary different model parameters to explore their influence on the planetesimal
growth. Only one parameter is varied in each case compared to the fiducial
case. For instance, we implement the type I migration and vary the pebble flux
and the Stokes number of pebbles (see Table 2). This parameter study is only
performed at $r{=}0.3$ AU, since it is the disk location where the
planetesimals actively grow. The parameter study at 2.7 AU, similar to 3 AU,
is performed in Liu et al. (2019b). In order to account for the statistical
variations, we run three different realizations for each case by randomizing
the initial locations and phase angles of forming planetesimals. The resultant
planets that are more massive than $0.01\ M_{\oplus}$ are shown in Figure 5.

Figure 5: Final planet masses vs. time when they reach these masses around
$r=0.3$ AU of the MMSN model. The planets with masses larger than 0.1 M⊕ are
presented, including the fiducial case (red), the type I migration case
(yellow), the high pebble flux case (blue), and the low Stokes number case
(black). The filled circles are bodies that reach their pebble isolation
masses within 1 Myr, while unfilled circles are bodies that experience
effective mass growth above $0.1$ M⊕ but fail to reach the pebble isolation
mass. Time is plotted in ranges, such as 0.1-0.2 Myrs or 0.2-0.3 Myrs. When
the plots are overlapped, they are replaced horizontally on the side in the
same time range.

First, we include type I migration to see how this process affects the growth
of planetesimals. We assume that planets stop their migration at $r{=}0.05$ AU
because of the gas clearing due to the central stellar magnetic torque (Liu et
al., 2017). In the migration case, the embryos migrate out of the birth
planetesimal belt once their masses grow beyond $0.01\ M_{\oplus}$ (Figure 6).
The random velocities of the remaining planetesimals are less excited
afterward, and these low-mass bodies can grow faster. The migrating embryos
grow their masses by pebble accretion and quickly reach the inner disk edge.
Since the pebble isolation mass is lower at a smaller radial distance, the
final planets in Figure 5 tend to have lower masses in the migration case
(yellow) compared to those in the fiducial case (red). Those embryos can
collide with each other actively through giant impacts (Figure 6b).

Figure 6: (a) Mass growth and (b) semimajor axis evolution of planetesimals in
a ring located at $r{=}0.3$ AU with type I migration. The same colored plots
and lines are used as for Figure 2. The embryos migrate out of the birth ring
when their masses exceed $10^{-3}\ M_{\oplus}$ and stop pebble accretion at
the inner edge of the disk. They can further increase their masses via giant
impacts.

The pebble flux is a representative of the available building blocks for
planet growth. We increase the pebble flux to $200$ M⊕ Myr-1 in the high
pebble flux case. We find in Figure 5 that Earth-mass planets reach the pebble
isolation mass at an earlier time in this case (blue), and their masses are
similar to those formed in the fiducial case. For instance, in the fiducial
case, the planets reach the pebble isolation mass between $t{\sim}0.3{-}0.4$
Myr, while in the high pebble flux case, the planets reach the pebble
isolation mass as early as $t{\sim}0.1{-}0.2$ Myr. This can be understood from
Equation (19) and Equation (20), since the growth timescale is inversely
proportional to the pebble flux.

We chose $\tau_{\rm s}{=}0.03$ in the low Stokes number case. In Figure 5,
planets also reach the pebble isolation mass at earlier times (black) of
$t{\sim}0.2{-}0.3$ Myr. The reasons are explained as follows. First, since the
threshold metallicity for triggering streaming instability $Z_{\rm SI}$
depends on the Stokes number, a lower $\tau_{\rm s}$ leads to a higher
metallicity through Equation (9) and, therefore, more massive birth
planetesimals (see Table 1). Second, the onset mass for pebble accretion is
also lower in the low Stokes number case. In Equation (18), the onset mass
with $\tau_{s}=0.03$ at $0.3$ AU is $M_{\rm onset}\simeq 9.5\times 10^{-6}$,
which is similar to $M_{\rm max}=1.9\times 10^{-5}$ (the maximum body from the
birth population). Thus, large planetesimals exceed the onset mass threshold
earlier for rapid pebble accretion in the low Stokes number case.

##  5 Planetesimal growth in the Toomre limit

####

The Atacama Large Millimeter/submillimeter Array (ALMA) has revealed that
rings and gaps are common in young protoplanetary disks at the orbital
distances of a few tens to hundreds of AUs (Avenhaus et al., 2018; Huang et
al., 2018; Long et al., 2019; Cieza et al., 2019). One leading explanation is
that these substructures are caused by embedded planets due to planet-disk
interactions (but also see Jiang & Ormel (2021)). This explanation requires
that the corresponding planetary mass is around Neptunian to Jovian mass
(Zhang et al., 2018; Bae et al., 2018; Liu et al., 2018; Lodato et al., 2019).
Here, we also attempt to explore whether massive protoplanets can form at
young ages of $1{-}3$ Myr from our described configuration, where the
streaming instability is triggered at a single disk location and a ring of
planetesimals form. In the very early disk phase, the disk surface density can
be much higher than in the fiducial MMSN model.

We estimate the upper limit of the gas surface density for a stable disk from
the Toomre criterion, where the Toomre parameter is given by (Toomre, 1964)

$Q\equiv\dfrac{c_{s}\Omega_{\rm K}}{\pi G\Sigma_{\text{gas}}}=h_{\rm
gas}\frac{M_{\star}}{\pi r^{2}\Sigma_{\rm gas}},$ (21)

where $c_{s}{=}H_{\rm gas}\Omega_{\rm K}$ is the gaseous sound speed, and
$Q{<}1$ corresponds to an unstable disk. Thus, the upper limit of the gas
surface density for a stable disk can be written as

$\Sigma_{\text{gas, Q}}=240\left(\frac{h_{\rm
gas,0}}{0.033}\right)\left(\dfrac{r}{30\rm\ AU}\right)^{-7/4}\text{g
cm${}^{-2}$}.$ (22)

The Toomre gas surface densities are 240 g/cm2 and 30 g/cm2 at $30$ AU and
$100$ AU, which are $24$ and $18$ times higher than the corresponding values
in the MMSN disk model. The planetesimals could have the highest masses if
streaming instability is triggered at such a high threshold disk density
condition. We test planetesimal formation in this extreme regime. The
properties of the formed planetesimals are summarized in Table 3. Since the
disk density is higher at such an early stage, we also explore cases with a
different pebble flux of $100$, $200,$ and $300\ \rm M_{\oplus}\ Myr^{-1}$. We
note that the mass of the planetesimal in this case can be treated as an upper
limit, and more realistic final masses of the planets should likely lie
somewhere between the Toomre and MMSN circumstances.

We performed simulations for two different cases: either with a planetesimal
mass distribution based on Equation (5) or with only the single largest
planetesimal from that distribution. In the former case, both the N-body
interaction and pebble accretion are taken into account, while in the latter
case, only pebble accretion is considered. We compare these two cases at
$t{=}1$ Myr. Figure 7 shows the masses of the embryos that grow by $50\%$ from
their initial masses and are greater than 0.1 M⊕ at the end of the
simulations.

We find, in Figure 7, that the planetesimal-planetesimal interaction has
little contribution to the mass growth of the most massive body. This is in
line with our expectation that pebble accretion is more predominant for the
subsequent planetesimal growth at the more distant disk location (Section 3).
We continue to simulate the growth of the single planetesimal for another $2$
Myr. Their masses at $t{=}2$ Myr and $3$ Myr are plotted in Figure 7 with
larger unfilled circles. For each case, the empty plots are denoted next to
each other, but they all have the same semimajor axes of 30 AU or 100 AU. For
comparison, the gap opening masses based on Eq. (35) of Liu et al. (2019a),
followed by Kanagawa et al. (2015), are also depicted with horizontal dotted
lines. With the explored parameters in the Toomre’s disk, we find that the
forming planets can have highest masses up to $54\ M_{\oplus}$ at $r{=}30$ AU,
which exceeds the gap opening mass, and up to $24\ M_{\oplus}$ at $r{=}100$ AU
within $3$ Myr. In a low turbulent disk of $\alpha_{\rm t}{=}10^{-4}$, the
pebble isolation mass obtained from Bitsch et al. (2018) is higher than the
gap opening mass calculated from Kanagawa et al. (2015). For illustrative
purposes, we continue our simulations in Figure 7, even when the planets reach
the gap opening mass.

Table 3:  Properties of planetesimals generated by streaming instability at $r{=}3$ and $30$ AU in an upper limit of the stable disk under the Toomre’s criterion. Radial distance | $N_{\text{tot}}$ | $N_{\text{sim}}$ |  $M_{\text{plt}}$ [M⊕] |  $M_{\text{min}}$ [M⊕] |  $M_{\text{max}}$ [M⊕] |  $M_{\text{tot}}$ [M⊕]
---|---|---|---|---|---|---
30 AU | 2225 | 199 | $2.6\times 10^{-3}$ | $1.3\times 10^{-4}$ | $1.29\times 10^{-1}$ | $1.7$
100 AU | 3096 | 273 | $6.4\times 10^{-3}$ | $3.2\times 10^{-4}$ | $3.7\times 10^{-1}$ | $4.2$
Figure 7: Final masses and semimajor axes of formed planets from planetesimal
rings at (a) $30$ AU and (b) $100$ AU orbital distances based on Toomre’s
surface density criterion. The filled dots represent the final bodies from the
planetesimal population (Eq. (5)) that grow by $50\%$ from their initial
masses and are greater than $0.1\rm M_{\oplus}$ at $t{=}1$ Myr. The unfilled
circles denote the single largest planetesimal that only undergoes pebble
accretion, for comparison. The three unfilled circles from small to large
refer to the mass at $t{=}1$, $2,$ and $3$ Myr, respectively. Red, blue, and
yellow correspond to the pebble flux of $100$, $200,$ and $300\ \rm
M_{\oplus}\ Myr^{-1}$, respectively. In order to distinguish the mass more
clearly, we artificially shift the blue and yellow unfilled circles slightly
away from the birth locations. The black dotted horizontal lines represent the
gap-opening masses based on Eq. (35) of Liu et al. (2019a) at $30$ AU and
$100$ AU.

##  6 Discussion

####

The disk turbulence plays a crucial role in regulating the rate of
planetesimal formation via streaming instability. Many theoretical studies
investigate how the growth of the streaming instability is influenced by disk
turbulence (Umurhan et al., 2020; Gole et al., 2020; Chen & Lin, 2020; Schäfer
et al., 2020; Xu & Bai, 2022). For instance, Chen & Lin (2020) find that,
since the dust particles get stirred by turbulence, the growth rate of
streaming instability is largely reduced. This suppression becomes
insignificant when $\alpha_{\rm t}$ is low and $\tau_{\rm s}$ is high. On the
other hand, Schäfer et al. (2020) find that the streaming instability becomes
more pronounced when it coexists with the vertical shear instability. A
similar positive effect is also confirmed by Xu & Bai (2022) in their
magnetohydrodynamic simulations. In these circumstances, depending on the
nature of turbulence, the disk exhibits local, temporary pressure bumps,
promoting the concentration of dust particles into clumps. Nevertheless, there
is a lack of systematic work on how the disk turbulence affects the birth mass
of forming planetesimals. Although our study assumes a weakly turbulent disk
of $\alpha_{\rm t}$, the initial planetesimal mass function is adopted from
Schäfer et al. (2017), where turbulence is not accounted for in their
simulations.

The formation of massive giant planets at distant orbits has also been studied
in other recent work (Rosenthal & Murray-Clay, 2018; Najita et al., 2021;
Chambers, 2021). For instance, Najita et al. (2021) consider various initial
disk solid mass and planetesimal formation efficiency and find that Pluto-,
Mars-, or super-Earth-mass planets form at greater than $30{-}40$ AU in their
evolutionary models. Chambers (2021) studies multiple pressure bumps in the
inner and outer regions of the disk with different birth masses of the
planetesimals and disk turbulent strength. The author finds that a Jupiter-
mass planet can form at $60$ AU within $0.5{-}3$ Myr in a $0.03\ M_{\odot}$
disk. Rosenthal & Murray-Clay (2018) also study the formation of a gas giant
at a wide orbit via a combination of planetesimal and pebble accretions with
different levels of turbulence of the disk. The authors find that a gas giant
can form at $\lesssim$ 40AU in a strongly turbulent disk and at $\lesssim$
70AU in a weakly turbulent disk.

Similarly, we also tested whether massive planets could form in the outer disk
region of $30$ and $100$ AU. We consider the upper limit of a stable disk with
Toorme’s criterion and tested different pebble fluxes. We find that several
super-Earth planets form within $1$ Myr when $100{-}300\ M_{\oplus}$ pebbles
drift across the orbit of the planetesimal ring. The giant planets of
$20{-}60\ M_{\oplus}$ form after $3$ Myr at the corresponding disk radii when
$900\ M_{\oplus}$ pebbles sweep by.

Noticeably, Najita et al. (2021) and Chambers (2021) investigate planet
formation with pressure bumps. They infer the total mass of forming
planetesimals from dust mass in the observed rings. Our study is, however, not
limited to this circumstance. Specifically, we consider the situation when
streaming instability is triggered at a single disk location. This is
irrelevant with the presence of pressure bumps. Meanwhile, in our work, the
total mass of planetesimals is calculated according to the streaming
instability triggering condition (Youdin & Goodman, 2005).

In addition, the starting mass of the planetesimal is the key free parameter
in other studies (Chambers, 2021). We self-consistently calculate the
planetesimal mass distribution based on realistic streaming instability
simulations (Schäfer et al., 2017; Liu et al., 2020). The masses of birth
planetesimals vary with disk radial distances. We fixed the turbulent strength
at $10^{-4}$. This corresponds to the optimized configuration of giant planet
formation in Chambers (2021).

Rosenthal & Murray-Clay (2018) suggest that the formation of gas giants at
wide orbits is difficult, due to inefficient planetesimal and pebble
accretion; but pebble accretion is still powerful enough to form smaller
planets. This is in line with our results. As pebble accretion is still
effective at forming a few Earth-mass planets in Rosenthal & Murray-Clay
(2018), we also find that super-Earth planets naturally form at $30$ and $100$
AU in a quiescent and denser disk. Although Rosenthal & Murray-Clay (2018) do
not find that massive cores trigger gas accretion, we find tens of Earth-mass
planets with a high pebble flux at the upper limit of Toomre’s stable disk,
within a growth time of $3$ Myr.

##  7 Conclusion

####

In this work, we investigate the planetesimal growth by a combination of
pebble and planetesimal accretion after their formation via the streaming
instability mechanism at a single disk location. A ring belt of planetesimals
is assumed to form in such a configuration (Liu et al., 2019b). We
particularly focus on how the growth pattern depends on radial distance $r$.
The planetesimals with initial mass distributions based on the streaming
instability simulations are generated at different disk radii of $0.3$, $3$,
$30,$ and $100$ AU. The characteristic mass of the birth planetesimals
increases with orbital distance and gas surface density (Liu et al., 2020). We
simulate their subsequent mass growth and dynamical evolution. The key
conclusions are listed as follows:

  * $\bullet$

We find that both planetesimal and pebble accretion are more active at smaller
orbital distances. Although they both become slower at more distant disk
locations, the radial distance dependence of pebble accretion is nevertheless
weaker than that of planetesimal accretion (Section 3).

  * $\bullet$

In the fiducial MMSN disk model, several planets with masses of
${\sim}0.1{-}3\ M_{\oplus}$ emerge after $1$ Myr at the close-in orbit of
$0.3$ AU (Figure 2). The planets have lower final masses at larger orbital
distances, with ${\sim}0.1\ M_{\oplus}$ at $3$ AU and ${\sim}10^{-3}\
M_{\oplus}$ at $30$ AU, respectively (Figures 3 and 4). The planetesimal
growth is strongly limited at $r{=}30$ AU for the MMSN model.

  * $\bullet$

The influence of different model parameters on the growth pattern is explored
in Section 4. Compared to the fiducial model, we find that when planetary
orbital migration is considered, a greater number of Earth-mass planets form
earlier around the inner edge of the disk (Figure 6).

  * $\bullet$

The planets reach their pebble isolation masses at earlier times when the
disks have a higher pebble flux or a lower Stokes number of pebbles (Figure
5).

  * $\bullet$

In the stable disks with Toomre’s gas surface density at $30$ and $100$ AU
(Section 5), planets can approach the gap opening mass within $2{-}3$ Myr
after their formation, depending on the pebble flux. The final planet can
attain $54\ M_{\oplus}$ at $30$ AU and $24\ M_{\oplus}$ at $100$ AU within $3$
Myr by a given pebble flux of $300\ M_{\oplus}\ \rm Myr^{-1}$ (Figure 7).
Since the generated planetesimals are massive enough to trigger rapid pebble
accretion in this high density disk condition, the growth is significantly
dominated by efficient pebble accretion.

In this paper, we focus on the radial distance dependence. In future work, we
will extend the study to look at the growth of planetesimals after the
streaming instability around different mass stars.

## Acknowledgments

We thank the anonymous referee for their insightful comments. H. J.
acknowledges the finical support from A.J. for the summer student internship,
which simulates the idea of this work. B. L. is supported by the start-up
grant of the Bairen program from Zhejiang University, National Natural Science
Foundation of China (No.12173035 and 12111530175), the Fundamental Research
Funds for the Central Universities (2022-KYY-506107-0001). A.J. acknowledges
funding from the European Research Foundation (ERC Consolidator Grant
724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg
Scholar Grant 2019.0442), the Swedish Research Council (Project Grant
2018-04867), the Danish National Research Foundation (DNRF Chair Grant
DNRF159) and the Göran Gustafsson Foundation.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44
  * Bae et al. (2018) Bae, J., Pinilla, P., & Birnstiel, T. 2018, ApJ, 864, L26
  * Benítez-Llambay & Pessah (2018) Benítez-Llambay, P. & Pessah, M. E. 2018, ApJ, 855, L28
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  * Bitsch et al. (2018) Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Chambers (2021) Chambers, J. 2021, ApJ, 914, 102
  * Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
  * Chatterjee & Tan (2014) Chatterjee, S. & Tan, J. C. 2014, ApJ, 780, 53
  * Chen & Lin (2020) Chen, K. & Lin, M.-K. 2020, ApJ, 891, 132
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698
  * Cresswell & Nelson (2008) Cresswell, P. & Nelson, R. P. 2008, A&A, 482, 677
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv e-prints, arXiv:2203.09759
  * Drażkowska & Alibert (2017) Drażkowska, J. & Alibert, Y. 2017, A&A, 608, A92
  * Drażkowska et al. (2021) Drażkowska, J., Stammler, S. M., & Birnstiel, T. 2021, A&A, 647, A15
  * Elbakyan et al. (2020) Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V., & Vorobyov, E. I. 2020, A&A, 637, A5
  * Eriksson et al. (2021) Eriksson, L. E. J., Ronnet, T., & Johansen, A. 2021, A&A, 648, A112
  * Eriksson et al. (2022) Eriksson, L. E. J., Ronnet, T., Johansen, A., et al. 2022, arXiv e-prints, arXiv:2202.12723
  * Flaherty et al. (2017) Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150
  * Gibbons et al. (2012) Gibbons, P. G., Rice, W. K. M., & Mamatsashvili, G. R. 2012, MNRAS, 426, 1444
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42
  * Hyodo et al. (2019) Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90
  * Jiang & Ormel (2021) Jiang, H. & Ormel, C. W. 2021, MNRAS, 505, 1162
  * Jiménez & Masset (2017) Jiménez, M. A. & Masset, F. S. 2017, MNRAS, 471, 4917
  * Johansen et al. (2014) Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 547
  * Johansen & Klahr (2005) Johansen, A. & Klahr, H. 2005, ApJ, 634, 1353
  * Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109
  * Johansen & Youdin (2007) Johansen, A. & Youdin, A. 2007, ApJ, 662, 627
  * Kanagawa et al. (2015) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2015, ApJ, 806, L15
  * Kley & Nelson (2012) Kley, W. & Nelson, R. P. 2012, ARA&A, 50, 211
  * Kokubo & Ida (1996) Kokubo, E. & Ida, S. 1996, Icarus, 123, 180
  * Lambrechts & Johansen (2012) Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32
  * Lambrechts & Johansen (2014) Lambrechts, M. & Johansen, A. 2014, A&A, 572, A107
  * Li & Youdin (2021a) Li, R. & Youdin, A. N. 2021a, ApJ, 919, 107
  * Li & Youdin (2021b) Li, R. & Youdin, A. N. 2021b, ApJ, 919, 107
  * Li et al. (2018) Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14
  * Liu & Ji (2020) Liu, B. & Ji, J. 2020, Research in Astronomy and Astrophysics, 20, 164
  * Liu et al. (2022) Liu, B., Johansen, A., Lambrechts, M., Bizzarro, M., & Haugbølle, T. 2022, Science Advances, 8, eabm3045
  * Liu et al. (2019a) Liu, B., Lambrechts, M., Johansen, A., & Liu, F. 2019a, A&A, 632, A7
  * Liu et al. (2020) Liu, B., Lambrechts, M., Johansen, A., Pascucci, I., & Henning, T. 2020, A&A, 638, A88
  * Liu & Ormel (2018) Liu, B. & Ormel, C. W. 2018, A&A, 615, A138
  * Liu et al. (2019b) Liu, B., Ormel, C. W., & Johansen, A. 2019b, A&A, 624, A114
  * Liu et al. (2017) Liu, B., Ormel, C. W., & Lin, D. N. C. 2017, A&A, 601, A15
  * Liu et al. (2018) Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87
  * Lodato et al. (2019) Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453
  * Long et al. (2019) Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49
  * Najita et al. (2021) Najita, J. R., Kenyon, S. J., & Bromley, B. C. 2021, arXiv e-prints, arXiv:2111.06406
  * 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
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
  * Ros & Johansen (2013) Ros, K. & Johansen, A. 2013, A&A, 552, A137
  * Rosenthal & Murray-Clay (2018) Rosenthal, M. M. & Murray-Clay, R. A. 2018, ApJ, 864, 66
  * Schäfer et al. (2020) Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190
  * 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
  * Surville et al. (2016) Surville, C., Mayer, L., & Lin, D. N. C. 2016, ApJ, 831, 82
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
  * Ueda et al. (2019) Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4
  * Visser & Ormel (2016) Visser, R. G. & Ormel, C. W. 2016, A&A, 586, A66
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Xu & Bai (2022) Xu, Z. & Bai, X.-N. 2022, ApJ, 924, 3
  * Yang & Johansen (2016) Yang, C.-C. & Johansen, A. 2016, ApJS, 224, 39
  * Yang et al. (2017) Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81
  * Zhu & Yang (2021) Zhu, Z. & Yang, C.-C. 2021, MNRAS, 501, 467
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57

##  Appendix A Variation of $M_{\rm min}$ in the planetesimal mass
distribution

####

In the main paper, $M_{\rm min}$ generated from the streaming instability is
set to be $20$ times smaller than the characteristic mass $M_{\rm plt}$. We
note that this factor (i.e., how low the mass of the planetesimal can be) is
actually less well constrained. For instance, the smallest size of the
planetesimal is largely limited by the grid resolution in the streaming
instability numerical simulations (Schäfer et al. 2017). Here, we vary
different values of $M_{\rm min}$ and re-generate the planetesimal populations
by a fixed total mass of planetesimal in the belt. The results of $M_{\rm
min}{=}0.01M_{\rm plt}$ are listed in Table 4. When comparing these
populations with those presented in Table 1, we find that the highest mass of
the planetesimal, as well as the mass distribution, are barely affected by the
variation of $M_{\rm min}$. This is valid as long as $M_{\rm min}{\ll}M_{\rm
plt}$.

Table 4: Planetesimal population based on Equation (5) by adopting $M_{\rm
min}{=}0.01M_{\rm plt}$.

$r$ |  $M_{\text{min}}$ [M⊕] |  $M_{\text{max}}$ [M⊕]
---|---|---
0.3 AU | $1.3\times 10^{-9}$ | $1.2\times 10^{-5}$
3 AU | $1.7\times 10^{-8}$ | $1.4\times 10^{-4}$
30 AU | $2.3\times 10^{-7}$ | $1.7\times 10^{-3}$
444The mass of the maximum body and the shape of the mass distribution are
still the same compared to the cases ($M_{\rm min}{=}0.05M_{\rm plt}$)
explored in the main text (see Table 1).

##  Appendix B Calculation of the mass distribution of planetesimals

####

Equation (5) is a cumulative number fraction of planetesimals, so we converted
the fractions to be noncumulative $N(m)/N_{\text{tot}}$ by subtracting next
number fraction bin. Then, we were able to obtain the number of planetesimals
from the noncumulative number fractions by multiplying the total number of
planetesimals. The total number of planetesimals is estimated from the sum of
all the mass fractions for all the mass bins, which is represented by
$M_{\text{tot}}/N_{\text{tot}}$. Therefore, the total number of planetesimals
is

$N_{\text{tot}}=\dfrac{M_{\text{tot}}}{M_{\text{tot}}/N_{\text{tot}}}=\dfrac{M_{\text{tot}}}{\sum\dfrac{N(m)}{N_{\text{tot}}}\times
m}.$ (23)

With this total number, we were able to convert the number fraction
($N/N_{\text{tot}}$) of planetesimals to the number of planetesimals for each
mass bin. Then, we determined the maximum mass of the local distribution to be
the mass that contains only one planetesimal. Therefore, we have the local
mass distributions from the minimum mass to the maximum mass. Since the
simulations only take integer number of planetesimals, we ignored the numbers
in decimal places.