11institutetext: Earth-Life Science Institute, Tokyo Institute of Technology,
Meguro-ku, Tokyo 152-8550, Japan
11email<EMAIL_ADDRESS>22institutetext: Laboratoire J.-L. Lagrange, Université
Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, F-06304 Nice, France
33institutetext: ISAS/JAXA, Sagamihara, Kanagawa, Japan
44institutetext: Department of Eartha and Planetary Sciences, Tokyo Institute
of Technology, Meguro-ku, Tokyo 152-8551, Japan
55institutetext: Steward Observatory/The Lunar and Planetary Laboratory,
University of Arizona, Tucson, Arizona 85721, USA.

# Planetesimal formation around the snow line: I. Monte Carlo simulations of
silicate dust pile-up in a turbulent disk

Shigeru Ida 11 Tristan Guillot 22 Ryuki Hyodo 33 Satoshi Okuzumi 44 Andrew N.
Youdin 55

(DRAFT: )

###### Abstract

Context.  The formation of rocky planetesimals is a long-standing problem in
planet formation theory. One of the possibilities is that it results from
gravitational instability as a result of pile-up of small silicate dust
particles released from sublimating icy pebbles that pass the snow line.

Aims.  We want to understand and quantify the role of the water snow line for
the formation of rock-rich and ice-rich planetesimals. In this paper, we focus
on the formation of rock-rich planetesimals. A companion paper examines the
combined formation of both rock-rich and ice-rich planetesimals.

Methods.  We develop a new Monte Carlo code to calculate the radial evolution
of silicate particles in a turbulent accretion disk, accounting for the back-
reaction (i.e., inertia) of the particles on their radial drift velocity and
diffusion. Results depend in particular on the particle injection width
(determined from the radial sublimation width of icy pebbles), the pebble
scale height and the pebble mass flux through the disk. The scale height
evolution of the silicate particles, which is the most important factor for
the runaway pile-up, is automatically calculated in this Lagrange method.

Results.  From the numerical results, we derive semi-analytical relations for
the scale height of the silicate dust particles and the particles-to-gas
density ratio at the midplane, as functions of a pebble-to-gas mass flux ratio
and the $\alpha$ parameters for disk gas accretion and vertical/radial
diffusion. We find that the runaway pile-up of the silicate particles
(formation of rocky planetesimals) occurs if the pebble-to-gas mass flux ratio
is $\ \lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\ [(\alpha_{{\rm
D}z}/\alpha_{\rm acc})/3\times 10^{-2}]^{1/2}$ where $\alpha_{{\rm D}z}$ and
$\alpha_{\rm acc}$ are the $\alpha$ parameters for vertical turbulent
diffusion and disk gas accretion.

###### Key Words.:

Planets and satellites: formation, Planet-disk interactions, Accretion,
accretion disks

##  1 Introduction

In order for planets to form, micron-size grains must grow to planetesimal
size (i.e., kilometers or more) before they are lost by gas drag into the
central star. One possibility is through streaming instabilities (SI), a
mechanism to concentrate small solid particles to high enough densities to
trigger gravitational collapse (e.g., Youdin & Goodman 2005). A midplane
particle-to-gas ratio of order unity is required to trigger strong clumping,
but this condition is not easy to achieve, even when considering the most
favorable particles, i.e. icy pebbles (e.g., Krijt et al. 2016). It is even
more difficult to achieve for smaller silicate dust particles, especially when
other sources of turbulence are present (Carrera et al. 2015; Yang et al.
2017; Gole et al. 2020).

The role of ice lines, and particularly of the water snow line, has long been
recognized as potentially important: The sublimation of pebbles across the
snow line leads to an outward diffusion of vapor and its recondensation, thus
increasing the pebble surface density (e.g. Stevenson & Lunine 1988; Ciesla &
Cuzzi 2006; Schoonenberg & Ormel 2017; Drażkowska & Alibert 2017).
Recondensation of outwardly diffused water vapor onto pebbles also induces
local pebble growth (Ros & Johansen 2013; Schoonenberg & Ormel 2017;
Drażkowska & Alibert 2017). A local pressure bump due to the ionization rate
change across the snow line caused by the sublimation of icy pebbles could
halt the pebble drift to result in their pile-up, under favorable disk
conditions (Kretke & Lin 2007; Brauer et al. 2008; Ida & Lin 2008). These
mechanisms potentially lead to the formation of ice-rich planetesimals111Our
subsequent papers, Hyodo et al. (2020, submitted: referred to as Paper II) and
Hyodo et al. (2020, submitted), investigate pile-ups of both icy pebble and
silicate dust particles.. But in the present work, our aim is to show that the
sublimation of icy pebbles can also result in the formation of silicate-rich
(rocky) planetesimals inside the snow line.

Many small silicate dust particles that are released by the sublimation of
individual icy pebbles are coupled with the disk gas (Section 2.4). Because
incoming icy pebbles drift fast by gas drag and pebbles are supplied with a
relatively high flux, the silicate dust particles can also pile up inside the
snow line, potentially leading to the conditions for a formation of silicate-
rich planetesimals by direct gravitational instability (Saito & Sirono 2011;
Ida & Guillot 2016). Determining the reality of the processes and their
relative contributions is of course crucial to understand the formation of
planetesimals, the composition of planets and ultimately account for the
global ice-to-rock ratio in the solar system (see Kunitomo et al. 2018).

However, the studies of this process (Ida & Guillot 2016; Schoonenberg & Ormel
2017; Drażkowska & Alibert 2017; Hyodo et al. 2019) lead to different,
seemingly contradictory results. In order to examine the problem in a new
light, we develop two tools. A first one is a 2D ($r$-$z$) Monte Carlo code to
simulate the pile-up of small silicate particles injected from sublimating icy
pebbles in a turbulent protoplanetary disk, and using the Lagrangian
advection-diffusion method by Ciesla (2010, 2011) with our new addition of the
back-reaction to the radial velocity and diffusion of the silicate particles.
A second one is a 1D diffusion-advection grid simulation based on the work of
Schoonenberg & Ormel (2017) as updated by Hyodo et al. (2019) and including
input results obtained from our Monte Carlo code. These tools enable us to
study the fate of silicate dust particles and icy pebbles for a wide range of
conditions, including the possibilities that turbulent mixing may differ in
the radial and vertical directions (Zhu et al. 2015; Yang et al. 2018).

We first focus on the 2D Monte Carlo code and its results. In a companion
paper (Paper II), we will apply the results to the diffusion-advection
simulation, allowing us to include the formation of both silicate-rich and
water-rich planetesimals. We will treat the case of complex protoplanetary
disk models in a third paper.

The present article is organized as follows: In Section 2, we summarize the
results of Ida & Guillot (2016), Schoonenberg & Ormel (2017), and Hyodo et al.
(2019), focusing on the silicate dust particle pile-up near the snow line. In
Section 3, we describe our Monte Carlo simulation code. In Section 4, after
testing our simulations, we present our results and derive semi-analytical
relations for the silicate dust scale height and the pile-up. We conclude in
Section 5.

##  2 Previous studies and analytical solution

###  2.1 Setting

The outskirts of protoplanetary disks are a large reservoir of condensed
particules, initially in the form of submicron-sized particles. Their
progressive growth and drift due to gas drag (Adachi et al. 1976;
Weidenschilling 1977) can lead to a wave of “pebbles”, i.e., centimeter- to
meter-sized particles that drift rapidly inward from the outer to the inner
regions of the disk (e.g. Garaud 2007; Lambrechts et al. 2014). The mass flux
of incoming pebbles can increase significantly above the standard value of
$\sim 1/100$ the gas mass flux that would be obtained for micron-size dust and
a solar composition gas. For example, Appelgren et al. (2020) obtain that, a
decrease in the drift rate of pebbles from the outer to the inner region can
lead to a $\sim 4$ times increase in the dust-to-gas ratio. Similarly, Mousis
et al. (2019) obtain an enrichment factor of the inner disk in volatiles that
can reach $\sim 20$ to $30$.

Figure 1:  Histogram of the time-dependent values of $F_{\rm
p/g}\equiv\dot{M}_{\rm p}/\dot{M}_{\rm g}$, the ratio of the pebble to gas
mass flux at the inner boundary of three protoplanetary disks (models L, M and
S, respectively) from Elbakyan et al. (2020). For this setting, irregularities
in the radial and azimuthal structures of the disks yield a highly time-
dependent $F_{\rm p/g}$ value.

The values of $F_{\rm p/g}$, the ratio of the mass flux of pebbles
$\dot{M}_{\rm p}$ to that of the gas $\dot{M}_{\rm g}$, depend on the disk
models and underlying assumptions on particle growth. Using 2D simulations of
the formation of a protoplanetary disk from the collapse of a molecular cloud
core, Elbakyan et al. (2020) obtained a highly time-dependent $F_{\rm p/g}$
that is shaped by the perturbations of the disk, as exemplified in Fig. 1. The
high values, much beyond the canonical value ($\sim 1/100$) are therefore
possible during the time-variation in each disk, either as a progressive
increase during the later evolution phases of the disk (Garaud 2007; Mousis et
al. 2019; Appelgren et al. 2020), or as short pulses (Elbakyan et al. 2020),
presumably due to the sudden release of pebbles.

Given the possible presence of such a pebble flux, Saito & Sirono (2011) had
pointed out that small silicate dust particles released by incoming icy
pebbles would pile up, potentially leading to the conditions for a direct
gravitational collapse of the dust subdisk. However, they had assumed a
stationary disk without inward disk gas accretion nor diffusion. Our goal is
to understand, depending on the value of $F_{\rm p/g}$, whether indeed
planetesimals may form from such a pile-up.

###  2.2 A comparison of previous studies

Several studies have thus far considered the problem in the framework of an
evolving accretion disk. Using an analytical approach, Ida & Guillot (2016)
found that even though the disk gas accretion tends to smooth out the radial
distribution of silicate dust particles, a runaway pile-up may be achieved
because of the back-reaction to inward drift of the particles, that is, the
slowing-down of the radial drift of the silicate particles due to the
increasing inertia of the piled-up particles. They derived a simple analytical
criterion linking the critical pebble-to-gas mass flux ratio for runaway pile-
up to a ratio between the dust-to-gas scale height and the fraction of dust
contained by icy pebbles (Section 2.3 and 4.1.2). Because Ida & Guillot (2016)
assumed that the scale height of silicate particles is similar to that of
incoming icy pebbles (i.e., much smaller than the disk gas scale height), and
neglected radial and vertical turbulent diffusion, their critical pebble
accretion flux corresponds to the most favorable limit to pile-up silicate
dust particles and form rock-rich planetesimals by gravitational instability.

Schoonenberg & Ormel (2017) performed detailed 1D diffusion-advection
simulations of gas and icy pebbles with an Eulerian grid code to study pile-
ups of the icy pebbles in an accretion disk, taking account of the back-
reaction to their drift velocity, sublimation process, and both radial and
vertical turbulent diffusion diffusion. Near the snow line, re-condensation of
water vapor diffused from the region inside the snow line enhances the ice
surface density. For icy pebbles with Stokes number of $10^{-2}-1$, SI can
occur for the solid-to-gas surface density ratio $Z_{\Sigma}\ \lower
3.0pt\hbox{${\buildrel>\over{\sim}}$}\ 0.02-0.03$ (Johansen et al. 2009;
Carrera et al. 2015; Yang et al. 2017). They found that the ice surface
density is enhanced by a factor of several by the recycling process in the
case of relatively vigorous turbulence ($\alpha\sim 10^{-3}-10^{-2}$) and a
relatively high pebble-to-gas mass accretion rate ($F_{\rm p/g}\sim 0.8$).

Schoonenberg & Ormel (2017) also calculated the radial distribution of
silicate particles released from drifting pebbles inside the snow line.
However, because they were focusing on the pebbles’ SI, they did not consider
the back-reaction to the radial drift velocity of the silicate particles. They
also assumed immediate turbulent stirring such that the scale height of
silicate particles is the same as the disk gas scale height. As a result, they
did not find a runaway pile-up of the silicate particles (for details, see
Section 2.3 and 4.1.2).

Hyodo et al. (2019) updated Schoonenberg & Ormel (2017)’s 1D diffusion-
advection grid simulations, adding the back-reactions to radial drift and
diffusion of silicate particles. They also considered the evolution of the
scale height of silicate dust particles. When the silicate particles are
released from icy pebbles passing the snow line, their scale height is similar
to that of icy pebbles. It is gradually increased by the vertical turbulent
diffusion. Hyodo et al. (2019) analytically modeled the dust scale height
evolution and incorporated the model into the grid simulation. They treated
the diffusion $\alpha$-parameter and the effective $\alpha$-parameter for the
disk accretion independently. The turbulent mixing and the angular momentum
transfer can be different mechanisms. For example, the angular momentum
transfer of the disk can be dominated by a process other than radial turbulent
diffusion such as disk winds (e.g., Suzuki et al. 2016; Bai et al. 2016).
Therefore, it is reasonable to distinguish the diffusion and the accretion
$\alpha$-parameters. They identified the parameter regions of the runaway
pile-up of the silicate particles on the plane of the pebble mass flux and the
diffusion $\alpha$-parameter.

As Ida & Guillot (2016) suggested, the silicate particle pile-up condition
depends directly on their scale height. Hyodo et al. (2019) estimated the
scale height evolution only by accounting for vertical diffusion, and not
taking radial diffusion into account.

###  2.3 Analytical derivation

As we described in Section 1, we assume that icy pebbles consists of a large
number of small silicate particles covered by icy mantles. When the pebbles
drift inward and pass the snow line, the icy mantles are sublimated and
consequently the small silicate particles are released. We set that the Stokes
number of icy pebbles is $\tau_{\rm s,p}=0.1$, which is relevant to the
Epstein regime (Okuzumi et al. 2012; Ida et al. 2016), and that of silicate
particles is $\tau_{\rm s,d}\ll 1$, where the Stokes number is defined by the
stopping time due to gas drag $t_{\rm stop}$ and Keplerian frequency
$\Omega_{\rm K}$ as $\tau_{\rm s}=t_{\rm stop}\Omega_{\rm K}$.

We consider a steady accretion disk. We show that our formulation and the
results are scaled by the pebble-to-gas mass flux ratio and we do not need to
specify the magnitude of the disk gas accretion rate. Actually, the pebble
mass flux through the disk ($\dot{M}_{\rm p}$) is proportional to disk gas
accretion rate ($\dot{M}_{\rm g}$) in the case of the steady accretion before
the pebble formation front reaches the disk outer edge (Ida et al. 2016).
Kanagawa et al. (2017) pointed out that the vertically averaged disk gas flow
can be outward, when $\dot{M}_{\rm p}$ is so large that the advective angular
momentum carried by pebbles dominates over the viscous angular momentum
transfer. In that case, the steady accretion disk model breaks down. However,
as we show in Appendix A, such a situation is not realized in the parameter
range we cover in this paper.

Table 1: List of notations notation | definition | eq. #
---|---|---
$f_{\rm d/p}$ | silicate dust mass fraction in a pebble |
| beyond the snow line [$=0.5$ (nominal)] |
$F_{\rm p/g}$ | pebble-to-gas mass flux in the disk | (4)
| $(=\dot{M}_{\rm p}/\dot{M}_{\rm g})$ |
$Z_{\Sigma}$ | vertically averaged dust-to-gas mass ratio | (3)
| $(=\Sigma_{\rm d}/\Sigma_{\rm g})$ |
$Z$ | local dust-to-gas mass ratio | (5)
| ($=\rho_{\rm d}/\rho_{\rm g}=Z_{\Sigma}/h_{\rm d/g}$) |
$\Lambda$ | parameter for the solid particle inertia | (10)
| $[=\rho_{\rm g}/(\rho_{\rm d}+\rho_{\rm g})=1/(1+Z)]$ |
$C_{\eta}$ | power index of pressure gradient | (8)
|  $[=-(1/2)(\partial\ln P/\partial\ln r);=11/8$(nominal)] |
$\alpha_{\rm acc}$ |  $\alpha$-parameter for disk gas accretion | (9)
$\alpha_{{\rm D}r}$ |  $\alpha$-parameter for radial mixing |
$\alpha_{{\rm D}z}$ |  $\alpha$-parameter for vertical stirring |
$(\rho_{\rm d}/\rho_{\rm g})_{0}$ | maximum value of $\rho_{\rm d}/\rho_{\rm g}$ in all the grids |
| (usually, at the grid of $(x,z)\sim(0,0)$) |
$h_{\rm d/g}$ | dust-to-gas scale height ratio $(=H_{\rm d}/H_{\rm g})$ | (12)
$h_{\rm d/g,0}$ |  $h_{\rm d/g}$ at $x$ with $(\rho_{\rm d}/\rho_{\rm g})_{0}$ |
$h_{\rm p/g}$ | pebble-to-gas scale height ratio $(=H_{\rm p}/H_{\rm g})$ | (19)
$\tau_{\rm s,p}$ | Stokes number of pebbles [$=0.1$ (nominal)] |
$\tau_{\rm s,d}$ | Stokes number of dust [$=10^{-5}$ (nominal)] |

First we derive the vertically averaged solid-to-gas mass ratio (metallicity),
inside the snow line. In this paper, the subscripts, “g”,“p”, and “d”
represent the disk H-He gas, icy pebbles, and silicate dust particles,
respectively. In the steady state, the surface density of the silicate dust
particles $\Sigma_{\rm d}$ and the disk gas $\Sigma_{\rm g}$ are given by

$\displaystyle\Sigma_{\rm d}$ $\displaystyle=f_{\rm d/p}\dot{M}_{\rm p}/2\pi
r\varv_{r},$ (1) $\displaystyle\Sigma_{\rm g}$ $\displaystyle=\dot{M}_{\rm
g}/2\pi ru_{r},$ (2)

where $\varv_{r}$ and $u_{r}$ are the drift velocity of the silicate particles
and the gas accretion velocity, and $f_{\rm d/p}$ is the silicate dust mass
fraction in the icy pebbles. We adopt $f_{\rm d/p}=0.5$ as a nominal value.
The vertically averaged metallicity is

$Z_{\Sigma}\equiv\frac{\Sigma_{\rm d}}{\Sigma_{\rm g}}=f_{\rm d/p}F_{\rm
p/g}\frac{u_{r}}{\varv_{r}},$ (3)

where

$F_{\rm p/g}\equiv\frac{\dot{M}_{\rm p}}{\dot{M}_{\rm g}}.$ (4)

We also use the local solid-to-gas ratio defined by

$Z\equiv\frac{\rho_{\rm d}}{\rho_{\rm g}},$ (5)

where $\rho_{\rm g}$ and $\rho_{\rm d}$ are the local densities of gas and
silicate dust particles.

The radial drift velocity of the particles and disk gas is given by Ida &
Guillot (2016) and Schoonenberg & Ormel (2017),

$\displaystyle\varv_{r}$ $\displaystyle=-\Lambda^{2}\frac{2\tau_{\rm
s,d}}{1+\Lambda^{2}\tau_{\rm s,d}^{2}}\eta\varv_{\rm
K}-\Lambda\frac{1}{1+\Lambda^{2}\tau_{\rm s,d}^{2}}u_{\nu},$ (6)

where $\varv_{\rm K}$ is Keplerian velocity, $\tau_{\rm s,d}$ is Stokes number
of the particles, $\eta$ is the degree of deviation of the gas rotation
angular velocity ($\Omega$) from Keplerian one ($\Omega_{\rm K}$), given by

$\displaystyle\eta$ $\displaystyle\equiv\frac{\Omega_{\rm
K}-\Omega}{\Omega_{\rm K}}=C_{\eta}\left(\frac{H_{\rm g}}{r}\right)^{2},$ (7)
$\displaystyle C_{\eta}$ $\displaystyle\equiv\frac{1}{2}\frac{\partial\ln
P}{\partial\ln r},$ (8)

where $C_{\eta}=1.3-1.4$, depending on the disk structure (Ida et al. 2016).
We use $C_{\eta}=11/8$ for $\Sigma_{\rm g}\propto r^{-1}$ and $T\propto
r^{-1/2}$ in this paper. In the above equation, $u_{\nu}$ is an unperturbed
disk gas accretion (advection) velocity given by

$\displaystyle u_{\nu}\simeq\frac{3\nu_{\rm acc}}{2r}\simeq\frac{3\alpha_{\rm
acc}H_{\rm g}^{2}\Omega}{2r}\simeq\frac{3\alpha_{\rm
acc}}{2}\left(\frac{H_{\rm g}}{r}\right)^{2}\varv_{\rm K},$ (9)

where $\nu_{\rm acc}$ is the effective viscosity for the angular momentum
transfer of the disk gas, $H_{\rm g}$ is the disk gas scale height, and
$u_{\nu}$ is defined to be positive for inward flow. The effect of back-
reaction from the piled-up particles is represented by $\Lambda$ that is
defined by

$\Lambda\equiv\frac{\rho_{\rm g}}{\rho_{\rm g}+\rho_{\rm d}}=\frac{1}{1+Z}.$
(10)

Because $\rho_{\rm d}\simeq\Sigma_{\rm d}/\sqrt{2\pi}H_{\rm d}$, where $H_{\rm
d}$ is the scale hight of the silicate dust particles, and $\rho_{\rm
g}\simeq\Sigma_{\rm g}/\sqrt{2\pi}H_{\rm g}$, Eq. (10) for $z<H_{\rm d}$ is
approximated for the particle motion by

$\Lambda^{-1}=1+\frac{\rho_{\rm d}}{\rho_{\rm g}}\simeq
1+Z_{\Sigma}\frac{H_{\rm g}}{H_{\rm d}}\simeq 1+h_{\rm d/g}^{-1}Z_{\Sigma}.$
(11)

We define the scale height ratio,

$h_{\rm d/g}\equiv H_{\rm d}/H_{\rm g}.$ (12)

In the case of silicate particle pile-up due to sublimation of icy pebbles,
the pile-up is radially localized just inside of the snow line. If we apply
Eq. (47), the gas accretion velocity ($u_{r}$) is locally slowed down from
$u_{\nu}$ by the back-reaction to the gas motion. If the steady accretion with
constant $\alpha_{\rm acc}$ is assumed, the gas surface density given by Eq.
(2) is locally increased. However, in reality, the assumption for vertically
constant $\alpha_{\rm acc}$ may be broken down in that case, or some local
instability may smooth out the local gas concentration. Therefore, we assume
that $u_{r}\simeq u_{\nu}$ also in the silicate dust pile-up region. Even if
we take account of the back-reaction to the gas motion, the condition for the
runaway pile-up of the silicate dust particles does not change significantly,
as shown in Appendix A.

Because $\tau_{\rm s,d}\ll 1$, silicate particle motions are strongly coupled
with disk gas, but the back-reaction to the drift velocity of the silicate
particles must be be taken into account. Equation (6) with $\tau_{\rm s,d}\ll
1$ implies $\varv_{r}\simeq\Lambda u_{\nu}$. Substituting $u_{r}\simeq
u_{\nu}$, $\varv_{r}\simeq\Lambda u_{\nu}$, and Eq. (11) into Eq. (3), we
obtain

$\displaystyle Z_{\Sigma}\simeq\frac{f_{\rm d/p}}{\Lambda}\,F_{\rm p/g}\simeq
f_{\rm d/p}(1+h_{\rm d/g}^{-1}Z_{\Sigma})F_{\rm p/g},$ (13)

which is solved as (Ida & Guillot 2016)

$\displaystyle Z_{\Sigma}=\frac{f_{\rm d/p}F_{\rm p/g}}{1-h_{\rm
d/g}^{-1}f_{\rm d/p}F_{\rm p/g}}.$ (14)

The midplane $Z=\rho_{\rm d}/\rho_{\rm g}$ is given by

$\displaystyle Z=\frac{Z_{\Sigma}}{h_{\rm d/g}}=\frac{f_{\rm d/p}F_{\rm
p/g}}{h_{\rm d/g}-f_{\rm d/p}F_{\rm p/g}}.$ (15)

If we consider the back-reaction to gas motion, the divergence condition is
only slightly modified to (see Appendix A):

$\displaystyle Z=\frac{f_{\rm d/p}F_{\rm p/g}}{h_{\rm d/g}-(1-h_{\rm
p/g})f_{\rm d/p}F_{\rm p/g}}.$ (16)

The metallicity $Z$ diverges, which would lead to rocky planetesimal
formation, for

$\displaystyle F_{\rm p/g}$ $\displaystyle>\frac{h_{\rm d/g}}{f_{\rm
d/p}},\hskip 14.22636pt[\mbox{w/o the back rection to gas}]$ (17)
$\displaystyle F_{\rm p/g}$ $\displaystyle>\frac{h_{\rm d/g}}{(1-h_{\rm
p/g})f_{\rm d/p}}.\hskip 14.22636pt[\mbox{w/ the back rection to gas}]$ (18)

The former condition was derived by Ida & Guillot (2016). We mostly refer to
the former condition.

Schoonenberg & Ormel (2017) neglected the back-reaction to the silicate dust
drift velocity and assumed $\Lambda=1$ in Eq. (13). In this case, Eq. (13) is
reduced to $Z_{\Sigma}\simeq f_{\rm d/p}F_{\rm p/g}$, which shows no
divergence. This demonstrates that the back-reaction to the silicate particle
radial velocity plays an essential role in the occurrence of the runaway pile-
up.

In the runaway pile-up condition, $F_{\rm p/g}>f_{\rm d/p}^{-1}\,h_{\rm d/g}$,
the estimation of $H_{\rm d}$ is essentially important. The scale height of
pebbles ($H_{\rm p}$) is given with their Stokes number $\tau_{\rm s,p}$ and
the vertical mixing parameter $\alpha_{{\rm D}z}$ as (Dubrulle et al. 1995;
Youdin & Lithwick 2007)

$\displaystyle h_{\rm p/g}=\frac{H_{\rm p}}{H_{\rm g}}=\left(1+\frac{\tau_{\rm
s,p}}{\alpha_{{\rm D}z}}\right)^{-1/2}.$ (19)

We will discuss the increase of $H_{\rm p}$ by Kelvin-Helmholtz instability
for the vertical shear due to pile-up of particles in Paper II.

The scale height of silicate dust particles ($H_{\rm d}$) must be the same as
$H_{\rm p}$ just after the release from the icy pebbles. Because $\tau_{\rm
s,d}\ll\tau_{\rm s,p}$, $H_{\rm d}$ is increased by the vertical turbulent
stirring afterward. Ida & Guillot (2016) assumed that $H_{\rm d}\simeq H_{\rm
p}$ ($h_{\rm d/g}=h_{\rm p/g}$), until the runaway pile-up develops. For
$h_{\rm d/g}\sim 0.1$ ($\tau_{\rm s,peb}/\alpha_{{\rm D}z}\sim 10^{2}$) and
$f_{\rm d/p}\sim 0.5$, the optimistic estimate by Eq. (15) predicts that the
runaway pile-up occurs for $F_{\rm p/g}>0.2$, which would be available for
relatively early phase of disk evolution.

On the other hand, Schoonenberg & Ormel (2017) assumed immediate vertical
stirring by turbulence, that is, $H_{\rm d}\simeq H_{\rm g}$ ($h_{\rm
d/g}\simeq 1$). Because they neglected the back-reaction to $\varv_{r}$, they
ever found the runway pile-up in their simulations. Even if we apply Eq. (15),
which was derived with the back-reaction, for the case of $h_{\rm d/g}\simeq
1$ and $f_{\rm d/p}\sim 0.5$, the runaway pile-up condition is $F_{\rm
p/g}>2$. The high value of $F_{\rm p/g}$ may not be completely ruled out, but
not easy to be established (Ida & Guillot 2016, and see Fig. 1 in this paper).

Hyodo et al. (2019) treated the effects of the vertical and radial diffusion
more carefully. They considered the vertical stirring characterized by
$\alpha_{{\rm D}z}$ to derive a model for $h_{\rm d/g}$ as a function of
$x=r-r_{\rm snow}$. They also introduced the back-reaction to the diffusion
coefficients $D_{r}$ and $D_{z}$ by the pile-up of the particles such that

$\displaystyle D_{r}$ $\displaystyle=\alpha_{{\rm D}r}H_{\rm
g}^{2}\Omega\Lambda^{K},$ (20) $\displaystyle D_{z}$
$\displaystyle=\alpha_{{\rm D}z}H_{\rm g}^{2}\Omega\Lambda^{K}.$ (21)

A simple argument that $D_{z},D_{r}\propto H_{\rm g}^{2}\propto
c_{s}^{2}\propto P/(\rho_{\rm g}+\rho_{\rm d})\propto\Lambda$ (Laibe & Price
2014; Lin & Youdin 2017) suggests $K=1$. Another argument is based on
Kolmogorov theory (e.g. Cho et al. 2003). In the eddy cascade, the energy
transfer rate, $(\rho_{\rm g}+\rho_{\rm d})\,\varv_{\ell}^{2}/t_{\rm eddy}$ is
independent of $\ell$, where $\ell,\varv_{\ell}$, and $t_{\rm eddy}$ are the
eddy size, velocity, and turnover time $\sim\ell/\varv_{\ell}$. Because
$D_{z},D_{r}\sim\ell_{0}\varv_{\ell 0}$ where $\ell_{0}$ and $\varv_{\ell 0}$
are the eddy size and velocity controlling the turbulence,
$D_{z},D_{r}\propto\ell_{0}\,[\ell_{0}/(\rho_{\rm g}+\rho_{\rm
d})]^{1/3}\propto\Lambda^{1/3}$. This argument thus implies $K=1/3$. This
should be a lower limit because an increase in the solid-to-gas is expected to
suppress turbulence itself. Given these uncertainties, Hyodo et al. (2019)
performed runs with $K=$ 0, 1, and 2. They found that the runaway pile-up is
caused by inclusion of the back-reaction to $\varv_{r}$, and that the
condition for the runaway pile-up depends on the relations among $\alpha_{\rm
acc}$, $\alpha_{{\rm D}r}$ and $\alpha_{{\rm D}z}$, while Ida & Guillot (2016)
considered the limits of $\alpha_{{\rm D}r},\alpha_{{\rm D}z}\rightarrow 0$
and their condition, Eq. (15), does not include $\alpha_{\rm acc}$. The
condition derived by Hyodo et al. (2019) is more severe than the prediction by
Ida & Guillot (2016). The runaway pile-up occurs for $F_{\rm p/g}>0.2-0.3$ and
$\alpha_{{\rm D}r}=\alpha_{{\rm D}z}\simeq 10^{-3}$ in the case of $K=1$ and 2
with $\alpha_{\rm acc}=10^{-2}$, while it does not occur with $K=0$ as long as
in the parameter space of $\alpha_{{\rm D}r}=\alpha_{{\rm
D}z}=10^{-3}-10^{-2}$ and $F_{\rm p/g}<0.6$. In this paper, we will show this
back-reaction effect on the diffusion coefficients more clearly.

###  2.4 Stokes number of silicate particles

The models discussed thus far and the model that we present are based on the
assumption that many silicate dust particles released by the individual icy
pebbles remain small compared to pebbles (i.e., with a Stokes number that is
much smaller than unity). In this case, the radial drift is much faster for
icy pebbles than for silicate particles. Because the gas accretion velocity is
given by $u_{\nu}\sim(3/2)\alpha_{\rm acc}(H_{\rm g}/r)^{2}\varv_{\rm K}$ (Eq.
(9)) and the drift speed relative to the gas is given by
$\varv_{r}\simeq\tau_{\rm s,d}(H_{\rm g}/r)^{2}\varv_{\rm K}$ (Eq. (6)), the
total drift velocity is dominated by $u_{\nu}$ if $\tau_{\rm s,d}\ \lower
3.0pt\hbox{${\buildrel<\over{\sim}}$}\ \alpha_{\rm acc}$. Therefore, the
assumption here is expressed as $\tau_{\rm s,d}\ \lower
3.0pt\hbox{${\buildrel<\over{\sim}}$}\ \alpha_{\rm acc}<\tau_{\rm s,p}$.

In Section 2.3, we assumed that $\tau_{\rm s,p}\sim 0.1$, $\alpha_{\rm
acc}\sim 10^{-3}-10^{-2}$, and $\tau_{\rm s,d}\ \lower
3.0pt\hbox{${\buildrel<\over{\sim}}$}\ \alpha_{\rm acc}$. Saito & Sirono
(2011) assumed that many $\mu$m-sized silicate particles are embedded in
pebbles. Morbidelli et al. (2015) considered chondrule-sized ($\sim$ mm)
particles. Because the Stokes number is proportional to the particle size in
the Epstein regime and its square in Stokes regime, it is reasonable to assume
that $\tau_{\rm s,d}<\alpha_{\rm acc}$. This model corresponds to the “many-
seeds model” of Schoonenberg & Ormel (2017). In the Monte Carlo simulations
here, we adopt $\tau_{\rm s,d}=10^{-5}$.

In the silicate particle pile-up region, the particles may quickly grow up to
sizes determined by the threshold collision velocity for
fragmentation/rebound. The collision velocity is set by the velocity
dispersion induced by the turbulence $\sim(3(\alpha_{{\rm
D}r}^{2}+\alpha_{{\rm D}z}^{2})^{1/2}\tau_{\rm s,d})^{1/2}c_{s}$ or by the
drift velocity difference between the particles $\sim 2\tau_{\rm
s,d}\,\eta\,\varv_{\rm K}$ (e.g., Sato et al. 2016). The collision velocity
between silicate particles is therefore $(\tau_{\rm s,d}/\tau_{\rm
s,p})^{-(0.5-1)}$ times of that between icy pebbles. If the threshold
collision velocity is 10 times lower for the silicate particles than for the
icy particles (Blum & Wurm 2000; Zsom et al. 2010, 2011; Wada et al. 2011;
Weidling et al. 2012; Wada et al. 2013), the silicate particles can grow only
up to sizes such that $\tau_{\rm s,d}\sim(10^{-1}$–$10^{-2})\tau_{\rm s,p}\sim
10^{-3}$–$10^{-2}$. Therefore, for $\alpha_{\rm acc}=10^{-3}$–$10^{-2}$, the
drift velocity of the silicate particles is comparable to $u_{\nu}$ and the
pile-up would not be significantly affected by the value of $\tau_{\rm s,d}$.

We note that this conventional view has been challenged recently (Kimura et
al. 2015; Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019).
A higher fragmentation threshold for silicates would open the possibility for
dust to coagulate and grow. We hence performed runs with larger values of
$\tau_{\rm s,d}$ to find that the results are not affected for $\tau_{\rm
s,d}\ \lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\ \alpha_{\rm acc}$ and that
the runaway pile-up region is removed from lower $F_{\rm p/g}$ regions as
$\tau_{\rm s,d}$ becomes larger.

Our model hence remains valid if the fragmentation velocity for rock-rich dust
remains smaller than for ice. We also note that an exploration of the
consequence of a different change of fragmentation velocity across ice lines
is explored by Vericel & Gonzalez (2019) (see also Gonzalez et al. 2017). In
this case, planetesimal formation may still be possible, but through a slow-
down of the drift rate due to the gas back-reaction and to a progressive
growth of the particles.

##  3 The Monte Carlo approach

In order to model precisely the effects of drift due to gas drag, gas
advection, radial and vertical diffusion, we develop a Monte Carlo simulation
of silicate dust particles in a turbulent accretion disk. The dust particles
released from sublimating icy pebbles are injected near the snow line. The
back-reactions to $\varv_{r}$ and the diffusion of the particles are included
using the super-particle approximation. We set radial and vertical
coordinates, $(x,z)$, where $x\equiv r-r_{\rm snow}$ is the radial distance
from the snow line at $r_{\rm snow}$ and $z$ is the distance from the
midplane.

We consider a radially local region near the snow line and neglect the
$r$-dependence of the Keplerian frequency ($\Omega$) and the disk gas scale
height ($H_{\rm g}$). As discussed in Section 2.4, we adopt the conventional
view to set $\tau_{\rm s,p}=0.1$ and $\tau_{\rm s,d}=10^{-5}$ in our
simulations.

###  3.1 Injection of particles

To model the release from drifting icy pebbles due to sublimation, at each
time step $\delta t$, we randomly inject a new silicate super-particle with
mass $m$ with a radially uniform distribution near the snow line in the range
of $-0.5\Delta x_{\rm subl}<x<0.5\Delta x_{\rm subl}$ and a Gaussian
distribution of $z$ of the root mean square $\Delta z_{\rm subl}$ as

$\displaystyle x$ $\displaystyle=0.5\;{\cal R}\;\Delta x_{\rm subl},$ (22)
$\displaystyle z$ $\displaystyle=\sqrt{2}\;{\rm erf}^{-1}(|{\cal
R}|)\;\frac{{\cal R}}{|{\cal R}|}\Delta z_{\rm subl},$ (23)

where ${\cal R}$ is a random number in the range of $[-1,1]$, $\Delta x_{\rm
subl}$ is the characteristic sublimation radial length, and $\Delta z_{\rm
subl}$ corresponds to the scale height of the incoming icy pebbles, $H_{\rm
p}$ given by Eq. (19). We adopt $\Delta x_{\rm subl}=0.1H_{\rm g}$ as a
nominal case. In Section 4.4.2, we also take a more complicated function that
fits the numerical result of the grid code simulation obtained by Paper II.

Here we distinguish the effective viscosity parameters for advection, radial
mixing, and vertical mixing, denoted by $\alpha_{\rm acc}$, $\alpha_{{\rm
D}r}$, and $\alpha_{{\rm D}z}$, respectively. According to different values of
$\alpha_{{\rm D}z}$, we use the consistent value of the scale height of the
injected silicate particles given by Eq. (19). Although Hasegawa et al. (2017)
suggested $\alpha_{{\rm D}r},\alpha_{{\rm D}z}\sim 0.1\,\alpha_{\rm acc}$, the
relations among $\alpha_{{\rm D}r},\alpha_{{\rm D}z}$ and $\alpha_{\rm acc}$
are not clear (e.g., Armitage et al. 2013). Therefore, we survey broad
parameter ranges of $\alpha_{{\rm D}r},\alpha_{{\rm D}z}$ and $\alpha_{\rm
acc}$; We mostly show the results with $\alpha_{\rm acc}\geq\alpha_{{\rm
D}r}=\alpha_{{\rm D}z}$. In some cases, the results with $\alpha_{{\rm
D}r}\neq\alpha_{{\rm D}z}$ ($\alpha_{\rm acc}\geq\alpha_{{\rm
D}r},\alpha_{{\rm D}z}$) are also shown.

###  3.2 Advection and diffusion

At each time step $\delta t$, we change the $r$–$z$ locations of the particles
that were injected before, following Ciesla (2010, 2011) by

$\displaystyle\delta x$ $\displaystyle=v_{{\rm adv},r}\delta t+{\cal
R}\,(6D_{r}\delta t)^{1/2},$ (24) $\displaystyle\delta z$
$\displaystyle=v_{{\rm adv},z}\delta t+{\cal R}\,(6D_{z}\delta t)^{1/2},$ (25)

where the 1st and 2nd terms in the right hand side represent advection (drift)
and diffusion, respectively. The root mean square of ${\cal R}$ is
$1/\sqrt{3}$. In our simulations, we adopt $\delta t=\Omega^{-1}$. Fromang &
Papaloizou (2006) suggested a turbulent correlation time in a protoplaneatry
disk is $\sim 0.15T_{\rm K}\sim\Omega^{-1}$. Ciesla (2011) adopted planar
($x$–$y$) random walks, in addition to vertical one ($z$), to take account of
the global curvature effect. In our case, because we only consider a local
radial range, we adopt a simpler $r$–$z$ random walks.

We perform runs with $K=0$ and $K=1$ for the diffusion coefficients (Eqs. (20)
and (21)), with $K=1$ being our nominal case. We also performed runs with
$K=1/3$ (the lower limit) and found that the overall features of the results
do not depend on the value of $K$ (only the runaway pile-up timescale is
different), as long as back-reaction to diffusion is considered (i.e., $K>0$).
Hereafter, we show the results with $K=1$ as representative of the cases
including back-reaction to diffusion.

The radial drift (advection) velocity $v_{{\rm adv},r}$ is the same as
$\varv_{r}$ in Eq. (6). The vertical advection velocity is (Ciesla 2010)

$v_{{\rm adv},z}=-(\alpha_{{\rm D}z}+\tau_{\rm s,d})\Omega\;z.$ (26)

The 1st term ($\propto\alpha_{{\rm D}z}$) comes from the effect that diffusion
acts to smooth the concentration ($\rho_{\rm d}/\rho_{\rm g}$), but not
$\rho_{\rm d}$ itself.

###  3.3 Inclusion of back-reaction effects

The advection and diffusion scheme was originally developed by Ciesla (2011)
and Ciesla (2011). Here, we also include the effects of back-reaction due to
the silicate particle pile-up as well as a scheme to model the injection of
silicate particles from sublimating icy pebbles.

We utilize the “super-particle” method, where one super-particle represents
the mass of a large number of particles while it suffers the same specific
drag force as the individual particles. We inject one super-particle at every
timestep $\delta t=\Omega^{-1}$. For given $\dot{M}_{\rm p}$, we can determine
the individual mass of the super-particles as

$m=f_{\rm d/p}\dot{M}_{\rm p}\,\delta t=f_{\rm d/p}\dot{M}_{\rm
p}\,\Omega^{-1},$ (27)

where $f_{\rm d/p}$ is a silicate fraction in migrating pebbles ($f_{\rm
d/p}=0.5$ in the nominal case). To locally average the silicate dust mass
density, we use linear grids, $\Delta x$ and $\Delta z$. In the nominal cases,
we adopt $\Delta x=0.1\,H_{\rm g}$. We use $\Delta z\sim\Delta z_{\rm subl}$
near the midplane and larger $\Delta z$ for upper regions to keep
statistically enough number of particles in the individual grids.

The silicate particle mass density at a grid of $[x,x+\Delta x]$ and
$[|z|,|z|+\Delta z]$ is

$\rho_{\rm d}=\frac{m\,\Delta N_{x,z}}{2\pi r\Delta x\times 2\Delta
z}=\frac{H_{\rm g}}{r}\frac{\Delta N_{x,z}\,f_{\rm d/p}\dot{M}_{\rm
p}}{4\pi\,(\Delta x/H_{\rm g})\,(\Delta z/H_{\rm g})\,H_{\rm g}^{3}\,\Omega},$
(28)

where $\Delta N_{x,z}$ is a total number of the particles in the grid. We
simply assume the vertical isothermal hydrodynamical equilibrium. Then the gas
density is given by

$\displaystyle\rho_{\rm g}$ $\displaystyle=\frac{\Sigma_{\rm
g}}{\sqrt{2\pi}H_{\rm g}}\exp\left(-\frac{z^{2}}{2H_{\rm g}^{2}}\right)$
$\displaystyle=\frac{\dot{M}_{\rm g}}{\sqrt{2\pi}\,3\pi\alpha_{\rm acc}H_{\rm
g}^{3}\,\Omega}\exp\left(-\frac{z^{2}}{2H_{\rm g}^{2}}\right).$ (29)

The silicate dust to gas ratio at $(x,z)$ is

$\displaystyle Z$ $\displaystyle=\frac{\rho_{\rm d}}{\rho_{\rm
g}}=\frac{3\sqrt{2\pi}}{4}\frac{H_{\rm g}}{r}\,\alpha_{\rm acc}\frac{f_{\rm
d/p}\dot{M}_{\rm p}}{\dot{M}_{\rm g}}\frac{\Delta N_{x,z}}{(\Delta x/H_{\rm
g})(\Delta z/H_{\rm g})}\exp\left(\frac{z^{2}}{2H_{\rm g}^{2}}\right)$
$\displaystyle=3.75\times 10^{-2}\left(\frac{H_{\rm
g}/r}{0.04}\right)\left(\frac{\alpha_{\rm
acc}}{10^{-2}}\right)\left(\frac{f_{\rm d/p}}{0.5}\right)F_{\rm p/g}$
$\displaystyle\hskip 14.22636pt\times\left(\frac{\Delta x}{0.1H_{\rm
g}}\right)^{-1}\left(\frac{\Delta z}{0.1H_{\rm g}}\right)^{-1}\Delta
N_{x,z}\exp\left(\frac{z^{2}}{2H_{\rm g}^{2}}\right).$ (30)

We count $\Delta N_{x,z}$ in the simulation bins of $[x,x+\Delta x]$ and
$[|z|,|z|+\Delta z]$ to calculate $Z$ according to Eq. (30). With calculated
$Z$, we update $\Lambda=1/(1+Z)$ in Eqs. (24) and (25) through Eqs. (6), (20),
and (21) for the next step.

In the simulations, we scale all the length by $H_{\rm g}$, including the
range of particle injection, $\Delta x_{\rm subl}$ and $\Delta z_{\rm subl}$,
and the grid sizes, $\Delta x$ and $\Delta z$. The time is scaled by
$\Omega^{-1}$, that is, each super-particle is injected at every unit scaled
time. Therefore, our results can be applied for any location of $r_{\rm
snow}$, as long as $H_{\rm g}/r$ is the same.

In the simulations here, we adopt $f_{\rm d/p}=0.5$ and $H_{\rm g}/r=0.04$ at
$r=r_{\rm snow}$, as nominal values. We do not need to specify the values of
$H_{\rm g}$ and $\dot{M}_{\rm g}$ as well as the value of $r_{\rm snow}$,
because the only scaled values, $h_{\rm d/g}(=H_{\rm d}/H_{\rm g}),h_{\rm
p/g}(=H_{\rm p}/H_{\rm g}),H_{\rm g}/r$, and $F_{\rm p/g}(=\dot{M}_{\rm
p}/\dot{M}_{\rm g})$, are used in the simulations. For the viscous
$\alpha$-parameters, we will show that only the ratios, $\alpha_{{\rm
D}z}/\alpha_{\rm acc}$ and $\alpha_{{\rm D}r}/\alpha_{\rm acc}$, are
important.

##  4 Results

###  4.1 Code check

We test our code by comparing the simulated particle scale heights with the
existing analytical argument, as Ciesla (2010) already did, and by reproducing
the simple analytical result by Ida & Guillot (2016) and the numerical result
without the back-reaction to silicate particles by Schoonenberg & Ormel
(2017).

We will also do detailed comparison with the results obtained by an updated
code of Hyodo et al. (2019) in Paper II. While Hyodo et al. (2019) included
both non-zero $\alpha_{{\rm D}r}$ and $\alpha_{{\rm D}z}$, they assumed the
evolution of $H_{\rm d}$ analytically estimated from only the vertical
stirring with $\alpha_{{\rm D}z}$. In our simulation, the evolution of $H_{\rm
d}$ is self-consistently calculated and we found that the radial mixing with
$\alpha_{{\rm D}r}$ is also important for $H_{\rm d}$. On the other hand, the
sublimation width $\Delta x_{\rm subl}$ is self-consistently calculated in
Hyodo et al. (2019), while it needs to be assumed in the simulation here.

####  4.1.1 Confirmation of the theoretically predicted particle scale height

Figure 2:  Snapshots of the particle distribution (red dots) on the $x$-$z$
plane at time $t=10^{4}\Omega^{-1}$ and the root mean square of $z$ (blue
solid curve) as a function of the radial direction $x$ (centered on the snow
line location and in units of the gas pressure scale height $H_{\rm g}$). We
plot $\mid z\mid$, while $z$ takes values of either positive or negative. We
set $\alpha_{\rm acc}=10^{-2}$, $\alpha_{{\rm D}r}=0$ and $\tau_{\rm
s}=10^{-3}$. The vertical mixing parameter is (a) $\alpha_{{\rm D}z}=10^{-5}$,
(b) $\alpha_{{\rm D}z}=10^{-3}$, and (c) $\alpha_{{\rm D}z}=10^{-1}$. The
analytical estimates of the equilibrium particle scale height for each
$\alpha_{{\rm D}z}$ are shown by the dashed green lines. The back-reactions
are not included to $\varv_{r}$ ($\Lambda=1$) and $\alpha_{{\rm D}z}$ ($K=0$),
We inject the particles in the range of $x=[-0.05,0.05]H_{\rm g}$ and
$z=[-0.1,0.1]H_{\rm g}$.

Figure 2 shows the distribution of particles and their scale height (the local
root mean square of $z$ of the particles) in steady state at
$t=10^{4}\Omega^{-1}$ (after $10^{4}$ particles are injected), obtained by our
Monte Carlo simulations for (a) $\alpha_{{\rm D}z}=10^{-5}$, (b) $10^{-3}$,
and (c) $10^{-1}$. For this runs, we neglected the back-reactions to
$\varv_{r}$ ($\Lambda=1$) and $\alpha_{{\rm D}z}$ ($K=0$). To highlight the
effect of vertical stirring, we set $\alpha_{{\rm D}r}=0$. The other parameter
are the particles’ Stokes number $\tau_{\rm s}=10^{-3}$ and $\alpha_{\rm
acc}=10^{-2}$. The the equilibrium scale height predicted from Eq. (19) with
$\tau_{\rm s,p}$ replaced by $\tau_{\rm s}=10^{-3}$ is (a) $0.1H_{\rm g}$, (b)
$0.71H_{\rm g}$ and (c) $1.0H_{\rm g}$, respectively, which are shown in the
green dashed lines in the plots. As the particles drift inward, the local root
mean square of $z$ (represented by the blue solid curves) asymptotically
approaches the theoretical values in each panel. Thus, we confirm that the
analytically derived Eq. (19) is reproduced by our Monte Carlo simulations.

####  4.1.2 Reproduction of Schoonenberg & Ormel (2017)’s and Ida & Guillot
(2016)’s results

Ida & Guillot (2016) adopted the following setup:

  * •

The turbulent diffusion is neglected for simplicity.

  * •

The silicate particle pile-up is fast enough that $H_{\rm d}$ is not increased
from $\Delta z_{\rm subl}(=H_{\rm p})$ by the vertical stirring.

To mimic the Ida & Guillot (2016)’s setup, we set the following in the Monte
Carlo simulations:

  * •

We adopt $\alpha_{{\rm D}r}=0$, $K=0$, and $\Delta z_{\rm subl}=0.03H_{\rm
g}$.

  * •

To keep $H_{\rm d}=\Delta z_{\rm subl}$ as an equilibrium, we adopt
$\alpha_{{\rm D}z}=10^{-5}$ and $\tau_{\rm s,d}=10^{-2}$.

  * •

The dust particle drift speed is set to be $\varv_{r}=-\Lambda u_{\nu}$, to
which the artificially enlarged $\tau_{\rm s,d}$ is not reflected (Eq. (6)).
For comparison, $\varv_{r}=-u_{\nu}$ (without the back-reaction to
$\varv_{r}$) is also examined. For gas, $u_{r}=-u_{\nu}$ is always assumed.

In the later results, we always set $\tau_{\rm s,p}=10^{-1}$ and $\tau_{\rm
s,d}=10^{-5}$ and consistently calculate $H_{\rm p}(=\Delta z_{\rm subl})$ by
Eq. (19) and the evolution of $H_{\rm d}$ from $\Delta z_{\rm subl}$ by the
vertical and radial diffusion, with given $\alpha_{{\rm D}z}$ and
$\alpha_{{\rm D}r}$. The above artificial setting is used only in the
particular runs in Figure 3.

Figure 3 shows the time evolution of the maximum value of $\rho_{\rm
d}/\rho_{\rm g}$ (the value on the midplane near the snow line) for different
values of $F_{\rm p/g}$, in the case of $\alpha_{\rm acc}=10^{-2}$, $\Delta
x_{\rm subl}=0.1H_{\rm g}$, and $\Delta z_{\rm subl}=0.03H_{\rm g}$. In panel
(a), the back-reaction to $\varv_{r}$ is not included ($\Lambda=1$) in the
Monte Carlo calculations, that is, $\varv_{r}=u_{r}=-u_{\nu}$, which
corresponds to the result in Schoonenberg & Ormel (2017) except that they
assumed $h_{\rm d/g}=1$. The maximum $\rho_{\rm d}/\rho_{\rm g}$ quickly
reaches the equilibrium values that are given by Eq. (3) with
$\varv_{r}=u_{r}=-u_{\nu}$ and $Z_{\Sigma}=Z\,h_{\rm d/g}$ as

$\frac{\rho_{\rm d}}{\rho_{\rm g}}\simeq f_{\rm d/p}h_{\rm d/g}^{-1}F_{\rm
p/g}\simeq 16.7\,F_{\rm p/g}.$ (31)

The open circles in panel (a) represent this analytical solution for each of
$F_{\rm p/g}$, which completely agree with the numerical results here.

In Figure 3b, the back-reaction to $\varv_{r}$ is included ($\Lambda<$),
corresponding to Ida & Guillot (2016)’s setup. They predicted that runaway
pile-up occurs for $F_{\rm p/g}>h_{\rm d/g}f_{\rm d/p}^{-1}$, which is $F_{\rm
p/g}>0.06$ for the parameters adopted here. This panel shows a clear runaway
pile-up for $F_{\rm p/g}\ \lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\ 0.1$,
which is consistent with the prediction.

Figure 3:  The time evolution of the maximum $\rho_{\rm d}/\rho_{\rm g}$. The
back-reaction to $\varv_{r}$ is included in the results in panel (b), but not
in panel (b). The light blue, magenta, blue, green, and red curves are the
results with $F_{\rm p/g}=0.025,0.05,0.1,0.2,0.4$ and 0.8, respectively. The
other parameters are fixed as $\alpha_{\rm acc}=10^{-2}$ and $\Delta x_{\rm
subl}=0.1H_{\rm g}$ and $\Delta z_{\rm subl}=0.03H_{\rm g}$. To mimic the
settings of Ida & Guillot (2016), we artificailly set $\alpha_{{\rm D}r}=0$,
$\alpha_{{\rm D}z}=10^{-5}$ and $\tau_{\rm s,d}=10^{-2}$ (however, the
relatively large $\tau_{\rm s,d}$ for the silicate dust is not reflected to
the dust drift speed). The latter two parameters are to maintain $H_{\rm
d}=\Delta z_{\rm subl}$. The circles in panel (a) represent the analytical
solution given by Eq. (31). That in panel (b) is the analytical solution for
non-divergent case ($F_{\rm p/g}=0.025$) given by Eq. (15).

Comparison between Figure 3a and b clearly shows that the back-reaction to
$\varv_{r}$ of the silicate particles plays an essential role in the
occurrence of the runaway pile-up. The increase in the pile-up slows down the
drift velocity and accordingly increases the pile-up itself. Because it is
certain that the back-reaction to $\varv_{r}$, we include it in the results in
the rest of the paper.

If the supply rate of silicate particles ($f_{\rm d/p}\dot{M}_{\rm p}$)
exceeds a threshold value, the mode of the pile-up becomes runaway. The
threshold value must be regulated by the local value of $\rho_{\rm
d}/\rho_{\rm g}$, because the back-reaction becomes effective once $\rho_{\rm
d}/\rho_{\rm g}$ exceeds unity. While $\dot{M}_{\rm p}$, more exactly $F_{\rm
p/g}$, determines $Z_{\Sigma}=\Sigma_{\rm d}/\Sigma_{\rm g}$, the local value
of $Z=\rho_{\rm d}/\rho_{\rm g}$ is regulated by $H_{\rm d}$ and $Z_{\Sigma}$
($Z=Z_{\Sigma}/h_{\rm d/g}$) near the snow line. The silicate dust particle
scale height $H_{\rm d}$ is regulated by the pebble scale height $H_{\rm p}$
($=\Delta z_{\rm subl}$ in our simulation) and turbulent mixing parmeters,
$\alpha_{{\rm D}z}$ and $\alpha_{{\rm D}r}$, as shown in Section 4.3.

###  4.2 Typical results

Figure 4:  Snapshots of silicate particles at $t=10^{4}\Omega^{-1}$ (red
dots), and the root mean square of $z$ (blue solid curves) at each grid of
$x$. The dashed blue curves are given by Eq. (35) with $\Delta x_{\rm subl}$
replaced by $\max(\Delta x_{\rm subl},|x|)$ (see the discussion at the end of
Section 4.3.2). The magenta curve represents $\log_{10}(\rho_{\rm d}/\rho_{\rm
g})$ near the midplane. We set $\alpha_{\rm acc}=10^{-2}$, $\tau_{\rm
s,d}=10^{-5}$, and $f_{\rm d/p}=0.5$, and $\dot{M}_{\rm
g}=10^{-8}M_{\odot}/{\rm yr}$. The nine panels adopt different values of
$F_{\rm p/g}$ and $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}$, as indicated in each
panel.  Figure 5:  The same as Fig. 4 except for $K=1$.

We first show typical results in our simulations for different $F_{\rm p/g}$
and $\alpha_{{\rm D}z}\,(=\alpha_{{\rm D}r})$. As we will show below, the
radial diffusion also contributes to $H_{\rm d}$. The back-reaction to
$\varv_{r}$ is included. In these results, $\Delta x_{\rm subl}=0.1H_{\rm g}$
is adopted, as well as $\alpha_{\rm acc}=10^{-2}$.

Figure 4 show the snapshots of silicate particles at $10^{4}\Omega^{-1}$ for
the cases without the back-reaction to diffusion coefficients ($K=0$). In the
right panels with $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}=10^{-2}$ ($\alpha_{{\rm
D}z}/\tau_{\rm s,d}=10^{3}$), the vertical mixing is fast enough to realize
the upper limit of $H_{\rm d}$ ($H_{\rm d}\simeq H_{\rm g}$) for all of
$F_{\rm p/g}=0.1,0.3$ and 1. Even for high mass flux of pebbles ($F_{\rm
p/g}=1.0$), no pile-up of the silicate particles is found for $\alpha_{{\rm
D}z}=\alpha_{{\rm D}r}=10^{-2}$. In the case of $\alpha_{{\rm
D}z}=\alpha_{{\rm D}r}=10^{-6}$, on the other hand, the mixing is so weak that
$H_{\rm d}$ is not increased to the equilibrium value [$\simeq(\alpha_{{\rm
D}r}/\tau_{\rm s,d})^{1/2}H_{\rm g}\sim 0.3H_{\rm g}$] in the range of
$x>-2\,H_{\rm g}$. Nevertheless, $H_{\rm d}$ near the injection point ($x\sim
0$) is much higher than $\Delta z_{\rm subl}=H_{\rm p}=0.003H_{\rm g}$. This
means that $H_{\rm d}$ is regulated by $\alpha_{{\rm D}r}$ as well as by
$\alpha_{{\rm D}z}$, which will be discussed in Section 4.3.

Figure 5 show the results with the back-reaction to diffusion coefficients
($K=1$). In general, the runaway pile-up is more pronounced for smaller
$\alpha_{{\rm D}z}=\alpha_{{\rm D}r}$ and larger $F_{\rm p/g}$. The magenta
curves represent $\log_{10}(\rho_{\rm d}/\rho_{\rm g})$ near the midplane. The
midplane $\rho_{\rm d}/\rho_{\rm g}$ usually takes the maximum value, at the
injection (sublimation) region, $x\sim 0$, which is denoted by $(\rho_{\rm
d}/\rho_{\rm g})_{0}$. In the panels with $(\rho_{\rm d}/\rho_{\rm g})_{0}\
\lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\ 1$ (the panels with
$\alpha_{{\rm D}z}=10^{-2}$ and those with $\alpha_{{\rm D}z}=10^{-4}$ and
$F_{\rm p/g}=0.3$ and 0.1), the results are similar to those in Figure 4. If
$\rho_{\rm d}/\rho_{\rm g}\ \lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\ {\rm
several}$ in Figure 4, however, the pile-up is much more pronounced in this
case, because the radial and vertical diffusions become much weaker as the
pile-up proceeds.

Figure 6: The time evolution of $(\rho_{\rm d}/\rho_{\rm g})_{0}$ and the dust
scale height at the maximum $\rho_{\rm d}/\rho_{\rm g}$ locations scaled by
the initial one ($=H_{\rm p}$). The back-reaction to diffusion is included
with $K=1$ in panel (b), but not ($K=0$) in panel (a), while the back-reaction
to $\varv_{r}$ is included in both panel (a) and (b). The light blue, magenta,
blue, green, and red curves are the results with $F_{\rm
p/g}=0.05,0.1,0.2,0.4$ and 0.8, respectively. The other parameters are fixed
as $\alpha_{\rm acc}=10^{-2}$, $\Delta x_{\rm subl}=0.1$, and $\Delta z_{\rm
subl}=0.03$. The diffusion parameters are $\alpha_{{\rm D}z}=\alpha_{{\rm
D}r}=10^{-5}$ in panel (a) and $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}=10^{-4}$
in panel (b).

Figure 6 show the time evolution of $(\rho_{\rm d}/\rho_{\rm g})_{0}$ and
$H_{\rm d}/H_{\rm p}$ at the maximum $\rho_{\rm d}/\rho_{\rm g}$ locations
($h_{\rm d/g,0}$), for different $F_{\rm p/g}$. In panel (a), the back-
reaction to diffusion is not included ($K=0$). As Figures 4 and 5 suggest, the
transient $\alpha_{{\rm D}z}$ $(=\alpha_{{\rm D}r})$ values between the
runaway pile-up and non-pile-up cases are $\alpha_{{\rm D}z}\simeq 10^{-5}$
for $K=0$ and $\alpha_{{\rm D}z}\simeq 10^{-4}$ for $K=1$, so that we plot the
results with $\alpha_{{\rm D}z}=10^{-5}$ in panel (a) and $\alpha_{{\rm
D}z}=10^{-4}$ in panel (b). These figures show that transition to the runaway
pile-up is much clearer in the case of $K=1$.

This difference comes from the $H_{\rm d}$ evolution. The comparison with the
result of Ida & Guillot (2016) in Section 2 suggests that $H_{\rm d}$ is
directly related to a threshold value of $F_{\rm p/g}$ for the runaway pile-
up: A lower $H_{\rm d}$ leads to a lower threshold value of $F_{\rm p/g}$.
However, Figure 6 show that $H_{\rm d}$ increases as $\rho_{\rm d}/\rho_{\rm
g}$ increases in the case of $K=0$, while it decreases for $K=1$. Therefore,
the pile-up suffers a negative feedback for $K=0$, but it suffers a positive
feedback, resulting in much clearer transition to the runaway pile-up for
$K=1$.

In both cases, the back-reaction to $\varv_{r}$ is included. As the pile-up
proceeds, $\varv_{r}$ decreases, which enhances the pile-up. However, we also
include the vertical stirring. As $\varv_{r}$ decreases, the silicate
particles are stirred up for longer time before they leave the sublimation
region, resulting in higher $H_{\rm d}$ near the sublimation region. In panel
(a), $H_{\rm d}$ near the sublimation region quickly asymptotes to an
equilibrium value for $F_{\rm p/g}\leq 0.2$ (lower $(\rho_{\rm d}/\rho_{\rm
g})_{0}$ cases); it increases with time for $F_{\rm p/g}=0.8$ (the runaway
pile-up cases). In panel (b), in the runaway pile-up cases ($F_{\rm p/g}=0.4$
and 0.8), $H_{\rm d}$ initially increases with time. However, after $\rho_{\rm
d}/\rho_{\rm g}$ exceeds $\sim 1$, the diffusion coefficients becomes so small
that the settling of the particles overwhelms the vertical stirring and
$H_{\rm d}$ decreases. Thus, the different responses to the pile-up between
$K=0$ and $K=1$ result in different $H_{\rm d}$ and $\rho_{\rm d}/\rho_{\rm
g}$ evolution in the pile-up cases. As we discussed in Section 3.2, $K=1$
(with the diffusion back-reaction) is more realistic.

###  4.3 Silicate particle scale height

####  4.3.1 Monte Carlo simulation results

Figure 7: The maximum $\rho_{\rm d}/\rho_{\rm g}$ and $H_{\rm d}/H_{\rm g}$ at
the maximum $\rho_{\rm d}/\rho_{\rm g}$ locations at $t=15000\Omega^{-1}$ as a
function of $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}$ for different $F_{\rm p/g}$.
The back-reaction to diffusion is included with $K=1$ in panel (b), but not
($K=0$) in panel (a). The other parameters are fixed as $\alpha_{\rm
acc}=10^{-2}$, $\Delta x_{\rm subl}=0.1H_{\rm g}$, and $\Delta z_{\rm
subl}H_{\rm p}$ given by Eq. (19). In the upper and bottom panels, the red
dashed curves are the analytical formulas given by Eqs. (67) and (35) with Eq.
(36). In the lower panels, the green dashed lines represents $\Delta z_{\rm
subl}/H_{\rm g}$.

Figure 7 show $(\rho_{\rm d}/\rho_{\rm g})_{0}$ and $h_{\rm d/g,0}$ at
$t=15000\Omega^{-1}$ as a function of $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}$
for different $F_{\rm p/g}$. This figure more clearly shows that the
asymptotic equilibrium values of $h_{\rm d/g,0}$ are independent of $F_{\rm
p/g}$ in the non runaway pile-up cases which is suggested by Figure 6. The
envelope curves in the plots for (a) $K=0$ and (b) $K=1$ match the analytical
formula derived in the next subsubsection (Eqs. (67) and (35) with Eq. (36)).
The equilibrium values are generally larger than the scale height of the
injected silicate particles $\Delta z_{\rm subl}$, due to the effects of
vertical and radial diffusion. In panels (a) and (b), the values of $h_{\rm
d/g,0}$ that deviate from the envelope curves correspond to the runaway pile-
up; $h_{\rm d/g,0}$ deviates to higher values for $K=0$ and lower values for
$K=1$. If we calculate on longer timescales, the deviations increase.

####  4.3.2 Analytical formula for the silicate scale height

The red dashed curves in Figure 7 are analytical formulas for the non pile-up
cases, which is derived as follows. The timescale to drift by $\Delta x_{\rm
subl}$ is

$\displaystyle t_{\rm drift}$ $\displaystyle\simeq\frac{\Delta x_{\rm
subl}}{\varv_{r}}$ $\displaystyle\simeq\frac{\Delta x_{\rm
subl}}{(3/2)\Lambda\alpha_{\rm acc}(H_{\rm g}/r)^{2}\varv_{\rm
K}}\simeq\frac{\Delta x_{\rm subl}/H_{\rm g}}{(3/2)\Lambda\alpha_{\rm
acc}(H_{\rm g}/r)\Omega}.$ (32)

As the pile-up proceeds ($\Lambda=1/(1+Z)$ decreases), the particle drift
becomes slower and accordingly $t_{\rm drift}$ becomes longer. The vertical
diffusion length during $t_{\rm drift}$ is

$\displaystyle z_{\rm diff}\simeq\sqrt{D_{z}t_{\rm drift}}$
$\displaystyle\simeq\left(\frac{2}{3}\alpha_{D,z}\Lambda^{K}H_{\rm
g}^{2}\Omega\times\frac{\Delta x_{\rm subl}/H_{\rm g}}{\Lambda\alpha_{\rm
acc}(H_{\rm g}/r)\Omega}\right)^{1/2}$
$\displaystyle\simeq\left(\frac{2}{3}\Lambda^{K-1}\frac{\Delta x_{\rm
subl}/H_{\rm g}}{H_{\rm g}/r}\frac{\alpha_{D,z}}{\alpha_{\rm
acc}}\right)^{1/2}H_{\rm g}.$ (33)

This expression suggests that $z_{\rm diff}$ becomes higher as the pile-up
proceeds ($\Lambda$ decreases) for the case of $K=0$. Hereafter we derive the
equilibrium value of $H_{\rm d}/H_{\rm g}$ for the non pile-up cases with
$\Lambda\simeq 1$ (the envelope curves in Figure 7). If $Z\gg 1$ for the
equilibrium value of $H_{\rm d}/H_{\rm g}$, the runaway pile-up occurs, so
that $Z$ keeps increasing and the equilibrium value of $H_{\rm d}/H_{\rm g}$
no more exists.

If only vertical diffusion is considered, the silicate particle scale height
scaled by $H_{\rm g}$ is predicted as

$\displaystyle h_{\rm d/g}=\frac{H_{\rm d}}{H_{\rm g}}$
$\displaystyle\simeq\frac{\max\left(\Delta z_{\rm subl},z_{\rm
diff}\right)}{H_{\rm g}}\simeq\frac{(\Delta z_{\rm subl}^{2}+z_{\rm
diff}^{2})^{1/2}}{H_{\rm g}}$ $\displaystyle\simeq\left(\frac{\Delta z_{\rm
subl}^{2}}{H_{\rm g}^{2}}+\frac{2}{3}\frac{\alpha_{{\rm D}z}}{\alpha_{\rm
acc}}\frac{\Delta x_{\rm subl}/H_{\rm g}}{H_{\rm g}/r}\right)^{1/2}.$ (34)

Furthermore, during $t_{\rm drift}$, the particles are radially mixed in the
range of $x_{\rm diff}\simeq\sqrt{D_{r}t_{\rm drift}}\sim(\alpha_{{\rm
D}r}/\alpha_{{\rm D}z})^{1/2}z_{\rm diff}$. In general, the dust scale height
increases with the distance from the snow line. If $x_{\rm diff}\gg\Delta
x_{\rm subl}$, the radial diffusion brings the high $z$ particles back to the
injection region and raises the value of $h_{\rm d/g}$ there. To take this
effect into account, we found through comparison with the Monte Carlo
simulation results that a good-fit analytical formula is obtained by
multiplying the above $h_{{\rm d/g}}$ by the factor of $[1+(x_{\rm
diff}/\Delta x_{\rm subl})^{2}]$. When $\alpha_{{\rm D}r}$ is comparable to
$\alpha_{\rm acc}$, the silicate scale height is radially smoothed out and the
effect of radial mixing for the scale height becomes weak. To take this into
account, we replace $\alpha_{{\rm D}r}/\alpha_{\rm acc}$ in $(x_{\rm
diff}/\Delta x_{\rm subl})^{2}$ by $(\alpha_{{\rm D}r}/\alpha_{\rm
acc})/[1+(C_{r,\rm diff}\,\alpha_{{\rm D}r}/\alpha_{\rm acc})^{2}]$ in the
above equation. Comparing the formula with the numerical results (Figs. 7 and
8), we empirically adopt $C_{r,\rm diff}=10$. With these modifications,

$\displaystyle h_{\rm d/g,*}\simeq$ $\displaystyle\left(\frac{\Delta z_{\rm
subl}^{2}}{H_{\rm g}^{2}}+\frac{2}{3}\frac{\alpha_{{\rm D}z}}{\alpha_{\rm
acc}}\frac{\Delta x_{\rm subl}/H_{\rm g}}{H_{\rm g}/r}\right)^{1/2}$
$\displaystyle\times\left[1+\frac{2}{3}\frac{(\alpha_{{\rm D}r}/\alpha_{\rm
acc})}{1+(C_{r,\rm diff}\,\alpha_{{\rm D}r}/\,\alpha_{\rm
acc})^{2}}\frac{1}{(H_{\rm g}/r)(\Delta x_{\rm subl}/H_{\rm g})}\right].$ (35)

Because $h_{\rm d/g}$ cannot exceed the equilibrium value $h_{\rm
d/g,eq}\simeq(1+\tau_{\rm s,d}/\alpha_{{\rm D}z})^{-1/2}$, the final form is

$h_{\rm d/g}\simeq(h_{\rm d/g,eq}^{-1}+h_{\rm d/g,*}^{-1})^{-1}.$ (36)

As the lower panels in Figure 7 show, the analytical formula reproduces the
envelope curves of $H_{\rm d}/H_{\rm p}$ obtained by the numerical simulations
for the non pile-up cases. We note that $h_{\rm d/g,0}$ and $(\rho_{\rm
d}/\rho_{\rm g})_{0}$ depend on the ratios of $\alpha_{{\rm D}r}/\alpha_{\rm
acc}$ and $\alpha_{{\rm D}z}/\alpha_{\rm acc}$, but not on their absolute
values.

Equation (35) is the analytical formula for $h_{\rm d/g}$ at $x\sim 0$. This
formula can be extrapolated to any $x$ by replacing $\Delta x_{\rm subl}$ with
$\max(\Delta x_{\rm subl},|x|,\epsilon)$, where $\epsilon$ is added to avoid
the divergence in the case of extremely small $\Delta x_{\rm subl}$; we adopt
$\epsilon=0.01\,H_{\rm g}$. In Figs. 4 and 5, the scale height distribution
obtained by the Monte Carlo simulations are well reproduced by the analytical
formula (the blue dashed curves) except for the runaway pile-up cases where
the back-reaction modulates the scale height as we discussed in this
subsection.

###  4.4 Silicate particle pile-up

####  4.4.1 Analytical formula for the silicate pile-up

We can also derive an analytical formula for $(\rho_{\rm d}/\rho_{\rm g})_{0}$
for the non runaway pile-up cases and for the phase before the runaway pile-up
proceeds. If the estimated $(\rho_{\rm d}/\rho_{\rm g})_{0}$ exceeds $\sim 1$,
the actual value of $\rho_{\rm d}/\rho_{\rm g}$ should deviate from the
estimate and increase with time.

Basically, $(\rho_{\rm d}/\rho_{\rm g})_{0}$ is predicted by Eq. (15) with the
scale hight formula, Eqs. (36) and (35). A possible detailed correction is to
include the effect of outward diffusion flux ($F_{\rm d,D}$) of silicate
particles beyond the snow line that can be found in the snapshots with
$\alpha_{{\rm D}r}=\alpha_{\rm acc}=10^{-2}$ in Figs. 4 and 5. As shown in
Appendix B, calibrating with the results of Monte Carlo simulations in Figs. 7
and 8, the effect can be done with $f_{\rm d/p}F_{\rm p/g}$ replaced by
$F_{\rm d,net}=f_{\rm d/p}F_{\rm p/g}-F_{\rm d,D}$, where

$F_{\rm d,net}=f_{\rm d/p}F_{\rm p/g}\left[1-\left(1+\frac{15}{2}\frac{H_{\rm
g}}{r}\frac{\alpha_{\rm acc}}{\alpha_{{\rm D}r}}\right)^{-1}\right].$ (37)

Accordingly, Eq. (15) is reduced to

$\displaystyle\left(\frac{\rho_{\rm d}}{\rho_{\rm g}}\right)_{0}=\frac{F_{\rm
d,net}}{h_{\rm d/g}-F_{\rm d,net}}.$ (38)

This formula reproduces the numerical results in Figs. 7 and 8, as long as
$\rho_{\rm d}/\rho_{\rm g}\ \lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\ 1$.
Note, however, that if the effect of sticking to icy pebbles by silicate
particles that diffuse beyond the snow line is considered, most of the
silicate grains beyond the snow line may eventually come back with the
pebbles, so that $F_{\rm d,net}$ would become more similar to $f_{\rm
d/p}F_{\rm p/g}$.

Figure 8: Same as Figure 7 except for $\alpha_{{\rm D}z}/\alpha_{{\rm D}r}$
relations and the value of $\Delta x_{\rm subl}$. The maximum $\rho_{\rm
d}/\rho_{\rm g}$ and $h_{\rm d/g}$ $(=H_{\rm d}/H_{\rm g})$ there for
different $F_{\rm p/g}$ are given as functions of (a) $\alpha_{{\rm
D}r}=\alpha_{{\rm D}z}$ with $\Delta x_{\rm subl}=0.1\,H_{\rm g}$, which is
the same as Figure 7b for a comparison, (b) $\alpha_{{\rm D}r}=\alpha_{{\rm
D}z}$ with $\Delta x_{\rm subl}=H_{\rm g}$, (c) $\alpha_{{\rm D}z}$ with
$\alpha_{{\rm D}r}=10^{-4}$ and $\Delta x_{\rm subl}=0.1\,H_{\rm g}$, and (d)
$\alpha_{{\rm D}r}$ with $\alpha_{{\rm D}z}=10^{-4}$ and $\Delta x_{\rm
subl}=0.1\,H_{\rm g}$. For all the cases, $K=1$ is adopted. In the upper and
bottom panels, the red dashed curves are the analytical formulas given by Eqs.
(35) and (67) with corresponding parameter values.

We have derived the analytical formulas for the silicate particle scale
height, Eqs. (35) and (36), and those for the pile-up, Eqs. (37) and (38). The
formulas make clear the intrinsic physics of the pile-up of the silicate
particles released from the drifting icy pebble by their sublimation, which is
a complicated process. Here we study the parameter dependences using the
analytical formulas and the simulations with $K=1$.

The formulas imply that the results are scaled by $\alpha_{{\rm
D}z}/\alpha_{\rm acc}$ and $\alpha_{{\rm D}r}/\alpha_{\rm acc}$, for given
$\Delta z_{\rm subl}$, $\Delta x_{\rm subl}$, and $H_{\rm g}/r$. We also
carried out the Monte Carlo simulations with $\alpha_{\rm acc}=10^{-3}$ and
confirmed that the plots in Figure 7 are almost the same as long as we use
$\alpha_{{\rm D}r}/\alpha_{\rm acc}$ ($=\alpha_{{\rm D}z}/\alpha_{\rm acc}$)
as the horizontal coordinate.

The sublimation width, $\Delta x_{\rm subl}$, is an independent parameter. The
results with $\Delta x_{\rm subl}=H_{\rm g}$ (Fig. 8b), are compared with
those with $\Delta x_{\rm subl}=0.1H_{\rm g}$ (Fig. 8a, which are identical to
Fig. 7b). In Fig. 8b, the injected particles are distributed in a 10 times
broader region than in the case of Fig. 8a. As a result, the vertical mixing
is more effective ($t_{\rm drift}$ is longer) and the pile-up (the maximum
$\rho_{\rm d}/\rho_{\rm g}$) is lower.

So far, the results with $\alpha_{{\rm D}r}=\alpha_{{\rm D}z}$ have been
shown. In Figure 8c, $\alpha_{{\rm D}r}$ is fixed and $\alpha_{{\rm D}z}$ is
independently changed. In Figure 8d, $\alpha_{{\rm D}z}$ is fixed. The
envelope curves for $H_{\rm d}/H_{\rm g}$ ($=h_{\rm d/g}$) are reproduced by
the analytical formula also in these cases, which strongly suggests that the
correction factors with the vertical and radial diffusion are physically
correct.

####  4.4.2 The condition for the silicate runaway pile-up

Figure 9: The contour maps of the maximum $\rho_{\rm d}/\rho_{\rm g}$ on the
plane of $\alpha_{{\rm D}z}/\alpha_{\rm acc}(=\alpha_{{\rm D}r}/\alpha_{\rm
acc})$ and $F_{\rm p/g}$ with (a) $\Delta x_{\rm subl}=0.1\,H_{\rm g}$, and
(b) $\Delta x_{\rm subl}$ given by Eq. (44). In both cases, $\alpha_{\rm
acc}=10^{-2}$ and $K=1$ are adopted. The left panels are the results of Monte
Carlo simulations at $t=15000\Omega^{-1}$ and the right ones are the results
using analytical formulas given by Eq. (67) with Eqs. (35), (36), (68), and
the corresponding $\Delta x_{\rm subl}$. The Monte Carlo simulations were
performed with a grid size of $\Delta\log_{10}(\alpha_{{\rm D}z}/\alpha_{\rm
acc})=\Delta\log_{10}F_{\rm p/g}=0.2$. We used more refined grids to make the
maps by the analytical formulas. The color scales are based on
$\log_{10}(\rho_{\rm d}/\rho_{\rm g})$.

The analytical formulas reproduce the numerical results for broad range of
multi-dimensional parameters. Figure 9a shows the contour maps of the maximum
$\rho_{\rm d}/\rho_{\rm g}$ on the plane of $\alpha_{{\rm D}z}/\alpha_{\rm
acc}(=\alpha_{{\rm D}r}/\alpha_{\rm acc})$ and $F_{\rm p/g}$ at $t=1.5\times
10^{4}\Omega^{-1}$ for $\Delta x_{\rm subl}=0.1\,H_{\rm g}$. In Figure 9b,
$\Delta x_{\rm subl}$ is given by Eq. (44). In both cases, $\alpha_{\rm
acc}=10^{-2}$ and $K=1$. The left and right panels are the results of Monte
Carlo simulations and the analytical results, respectively. The red regions
represent the parameter ranges of runaway pile-up of silicate particles. As
expected, the analytical formulas reproduce the numerical results.

For these simulations, $\alpha_{\rm acc}=10^{-2}$ and $\alpha_{{\rm
D}z}(=\alpha_{{\rm D}r})$ is in the range of $[10^{-6},10^{-2}]$. The
analytical formula to predict $(\rho_{\rm d}/\rho_{\rm g})_{0}$ is Eq. (37)
with Eqs. (38), (35), and (36). They are functions of $\alpha_{{\rm
D}r}/\alpha_{\rm acc}$ and $\alpha_{{\rm D}z}/\alpha_{\rm acc}$ for given
$\Delta x_{\rm subl}$, as long as the term $\Delta z_{\rm subl}^{2}$ in Eq.
(35) is negligible, which is satisfied in the parameter range we consider as
shown in Figs. 7 and 8, We also carried out simulations with $\alpha_{\rm
acc}=10^{-3}$ and $\alpha_{{\rm D}z}(=\alpha_{{\rm D}r})$ in the range of
$[10^{-7},10^{-3}]$ to find that the contour map is almost identical,
confirming that the results are scaled by $\alpha_{{\rm D}z}/\alpha_{\rm
acc}(=\alpha_{{\rm D}r}/\alpha_{\rm acc})$.

Using the formulas, we derive the boundary of the runaway pile-up on the plane
of $\alpha_{{\rm D}z}/\alpha_{\rm acc}$ $(=\alpha_{{\rm D}r}/\alpha_{\rm
acc})$ and $F_{\rm p/g}$. Equation (38) shows that the boundary is given by
$F_{\rm d,net}\simeq h_{\rm d/g}$. As shown in Figure 9a and b, the boundary
is located at $\alpha_{{\rm D}z}(=\alpha_{{\rm D}r})<3\times 10^{-4}$ for
$F_{\rm p/g}<1$. In this parameter range, Figure 7 shows $h_{\rm d/g}\simeq
h_{\rm d/g,*}$. Since $\alpha_{{\rm D}r}/C_{r,\rm diff}\,\alpha_{\rm acc}\ll
1$ for $\alpha_{\rm acc}\sim 10^{-2}$, $F_{\rm d,net}\simeq f_{\rm d/p}F_{\rm
p/g}$. Therefore the boundary is approximately given by

$\displaystyle F_{\rm p/g}$ $\displaystyle\simeq f_{\rm d/p}^{-1}h_{\rm
d/g,*}$ $\displaystyle\simeq f_{\rm
d/p}^{-1}\left(\frac{2}{3}\frac{\alpha_{{\rm D}z}}{\alpha_{\rm
acc}}\frac{\Delta x_{\rm subl}/H_{\rm g}}{H_{\rm
g}/r}\right)^{1/2}\left[1+\frac{2}{3}\frac{\alpha_{{\rm D}r}}{\alpha_{\rm
acc}}\frac{1}{(H_{\rm g}/r)(\Delta x_{\rm subl}/H_{\rm g})}\right].$ (39)

Adopting $H_{\rm g}/r=0.04$ and $f_{\rm d/p}=1/2$ as nominal values,

$\displaystyle F_{\rm p/g}$ $\displaystyle\simeq$ $\displaystyle
0.13\left(\frac{f_{\rm d/p}}{1/2}\right)^{-1}\left(\frac{\alpha_{{\rm
D}z}/\alpha_{\rm acc}}{10^{-2}}\right)^{1/2}\left(\frac{H_{\rm
g}/r}{0.04}\right)^{1/2}\left(\frac{\Delta x_{\rm subl}}{0.1\,H_{\rm
g}}\right)^{1/2}$ (42)
$\displaystyle\times\left[1+\frac{5}{3}\left(\frac{\alpha_{{\rm
D}r}/\alpha_{\rm acc}}{10^{-2}}\right)\left(\frac{H_{\rm
g}/r}{0.04}\right)^{-1}\left(\frac{\Delta x_{\rm subl}}{0.1\,H_{\rm
g}}\right)^{-1}\right]$ $\displaystyle\propto$
$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\left(\frac{\alpha_{{\rm
D}z}}{\alpha_{\rm acc}}\right)^{1/2}\left(\frac{\alpha_{{\rm D}r}}{\alpha_{\rm
acc}}\right)\Delta x_{\rm subl}^{-1}&\displaystyle[\frac{\alpha_{{\rm
D}z}}{\alpha_{\rm acc}}>0.6\times 10^{-2}],\\\
\displaystyle\left(\frac{\alpha_{{\rm D}z}}{\alpha_{\rm
acc}}\right)^{1/2}\Delta x_{\rm subl}^{1/2}&\displaystyle[\frac{\alpha_{{\rm
D}z}}{\alpha_{\rm acc}}<0.6\times 10^{-2}].\end{array}\right.$

This equation explicitly shows that the pile-up condition is scaled by
$\alpha_{{\rm D}r}/\alpha_{\rm acc}$ and $\alpha_{{\rm D}z}/\alpha_{\rm acc}$.

In Eq. (42), the sublimation width, $\Delta x_{\rm subl}$, is a free
parameter. We have used $\Delta x=0.1\,H_{\rm g}$ as a nominal parameter. The
sublimation width is determined by the sublimation and drift rates of icy
pebbles. The sublimation rate is given by the partial pressure of water vapor
and the disk temperature. The water vapor partial pressure is affected by
radial diffusion. In Paper II, we numerically solve the sublimation width
taking account of the icy pebble size evolution and the water vapor partial
pressure in a turbulent accretion disk. Fitting the numerical results, the
sublimation width for $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}$ and $\alpha_{\rm
acc}=10^{-2}$ is roughly approximated as

$\displaystyle\log_{10}\left(\frac{\Delta x_{\rm subl}}{H_{\rm g}}\right)$
$\displaystyle\simeq\frac{X_{+}-X_{-}}{2}{\rm
erf}\left[3\log_{10}\left(\frac{\alpha_{{\rm D}r}/\alpha_{\rm
acc}}{0.08}\right)\right]+\frac{X_{+}+X_{-}}{2},$ $\displaystyle
X_{-}=\log_{10}(\Delta x_{\rm subl}/H_{\rm g})_{-}=\log_{10}2;$ (43)
$\displaystyle X_{+}=\log_{10}(\Delta x_{\rm subl}/H_{\rm
g})_{+}=\log_{10}0.1,$ (44)

where $(\Delta x_{\rm subl}/H_{\rm g})_{-}=2$ and $(\Delta x_{\rm subl}/H_{\rm
g})_{+}=0.1$ are the lower $\alpha_{{\rm D}r}/\alpha_{\rm acc}$ limit where
advection dominates and higher one where radial diffusion dominats,
respectively (Paper II). The error function smoothly connects the two limits.
Figures 9 show that the threshold value of $F_{\rm p/g}$ is several times
higher than in Fig. 9a for $\alpha_{{\rm D}z}/\alpha_{\rm acc}\sim 10^{-4}$,
while it is similar for $\alpha_{{\rm D}z}/\alpha_{\rm acc}\sim 10^{-2}$. In
the former case, $F_{\rm p/g}\propto\Delta x_{\rm subl}^{1/2}$ (Eq. 42), so
that $F_{\rm p/g}$ must be larger in Fig. 9b than in Fig. 9a by a factor of
$(2/0.1)^{1/2}\sim 4.5$ (Eq. (44)), while the dependence on $\Delta x_{\rm
subl}$ must be slightly weaker in the latter case (Eq. (44)). Thereby, the
silicate particle runaway pile-up condition with Eq. (44) is roughly given by

$F_{\rm p/g}\ \lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\
\left(\frac{\alpha_{{\rm D}z}/\alpha_{\rm acc}}{3\times
10^{-2}}\right)^{1/2},$ (45)

which agrees with the results in Fig. 9.

Figure 9 show the results with $\alpha_{{\rm D}r}=\alpha_{{\rm D}z}$. MHD
simulations showed $\alpha_{{\rm D}r}<\alpha_{{\rm D}z}$ (Zhu et al. 2015;
Yang et al. 2018). We also performed simulations with $\alpha_{{\rm
D}r}=0.1\,\alpha_{{\rm D}z}$. As we can predict from Eq. (42), the runaway
pile-up region is extended to a several times larger value of $\alpha_{{\rm
D}z}/\alpha_{\rm acc}$ at $F_{\rm p/g}\sim 1$ in the case of $\Delta x_{\rm
subl}=0.1\,H_{\rm g}$, while the runaway pile-up region is almost the same for
smaller values of $F_{\rm p/g}$. Because $\alpha_{{\rm D}z}/\alpha_{\rm acc}$
is very uncertain, we surveyed a rather broad parameter space and leave the
estimate of $\alpha_{{\rm D}z}/\alpha_{\rm acc}$ for future study.

####  4.4.3 Comparison with Hyodo et al. (2019)

Hyodo et al. (2019) performed 1D advection-diffusion grid code simulations
with $\alpha_{\rm acc}=10^{-2}$ and $\alpha_{{\rm D}r}=\alpha_{{\rm
D}z}=10^{-3},3\times 10^{-3},10^{-2}$ to find that the silicate runaway pile-
up occurs for $\alpha_{{\rm D}z}=\alpha_{{\rm D}r}\approx 10^{-1}\alpha_{\rm
acc}$ and $F_{\rm p/g}\geq 0.3$ in the case of $K=1$. However, Figure 9b shows
that the silicate runaway pile-up region is restricted to a smaller parameter
range, i.e., to values of $\alpha_{{\rm D}z}/\alpha_{\rm acc}=\alpha_{{\rm
D}r}/\alpha_{\rm acc}$ which are one order of magnitude smaller for a given
$F_{\rm p/g}$ value. As already mentioned, this is because Hyodo et al. (2019)
assumed that the silicate scale height is $\sim\Delta z_{\rm subl}$ at the
snow line and it is gradually increased by the vertical turbulent stirring as
the particles drift inward whereas the Monte Carlo simulations performed here
show that the silicate scale height is larger, due to the effect of a coupled
radial and vertical diffusion in a finite-size sublimation width. As a result,
smaller $\alpha_{{\rm D}z}$ and $\alpha_{{\rm D}r}$ are required for a runway
pile-up.

##  5 Conclusion and Discussions

The runaway pile-up of silicate particles released from sublimating icy
pebbles that pass through the snow line is a potential mechanism to form rock-
rich planetesimals. Ida & Guillot (2016) showed that the back-reaction
(inertia) of silicate particles to gas drag can lead to a runaway pile-up of
dust particles. They provided a simple criterion for this to occur, as a
function of both the pebble-to-gas mass-flux ratio and the silicate particle
scale-height. Schoonenberg & Ormel (2017) found instead that a runaway pile-up
of ice-rich pebbles would be possible and that of silicate particles does not
occur. However they did not include the back-reaction of dust particles, and
thus found no pile-up of dust particles. Hyodo et al. (2019) performed 1D
diffusion-advection grid code simulations including turbulent diffusion and
the back-reaction to radial drift and diffusion for both icy pebble and
silicate particles. They allowed the local turbulent diffusion governing
radial and vertical diffusion to differ from the turbulent viscosity
controlling the gas mass flux in the disk. They found that a runaway pile-up
of either dust or pebbles could be achieved, depending on the values of
turbulent diffusion and of the pebble-to-gas mass-flux ratio. However, they
had to approximate the calculation of the silicate particle scale-height based
on simple arguments. In this work, we revisited this issue.

We have developed a new 2D ($r$-$z$) Monte Carlo code to simulate the pile-up
of small silicate particles released from sublimating pebbles in a turbulent
protoplanetary disk, taking account of the back-reactions to the drift
velocity and the diffusion of silicate particles. From the simulation results,
we have derived semi-analytical formulas for the maximum silicate-to-gas
density ratio near the injection region and the silicate scale height there as
a function of the pebble mass flux, $\alpha_{\rm acc}$, $\alpha_{{\rm D}r}$,
and $\alpha_{{\rm D}z}$. Using the derived formulas, we determined the
detailed condition for the silicate runaway pile-up. We found that the
silicate particle scale height is larger than the Hyodo et al. (2019)’s
estimate due to the coupled effect of radial and vertical diffusion, and that
for the silicate particle runaway pile-up to occur, $\alpha_{{\rm
D}z}/\alpha_{\rm acc}\ \lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\
10^{-2}\times$ (pebble-to-gas mass flux)2 is required, which is more
restrictive than Hyodo et al. (2019)’s result ($\alpha_{{\rm D}z}/\alpha_{\rm
acc}\ \lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\ 10^{-1}$ and $F_{\rm p/g}\
\lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\ 0.3$). To clarify if the
condition is actually satisfied, detailed non-ideal MHD simulations are needed
to evaluate $\alpha_{{\rm D}z}$, $\alpha_{{\rm D}r}$, and $\alpha_{\rm acc}$.

Thus far, we have not included the pile-up of icy pebbles. This would occur
upstream, i.e., beyond the snow line, and could thus suppress the supply of
pebbles and thus dust grains, affecting the pile-up of dust-rich
planetesimals. It is also important to understand when and where runaway pile-
ups of icy pebbles and silicate particles may occur in the course of disk
evolution. We will investigate these issues in Paper II.

One would want to extend the arguments developed in this work to other ice
lines, such as for NH3 and CO2. The pile-up process proposed here occurs only
if the Stokes numbers of the particles before and after the sublimation of a
volatile component ($\tau_{\rm s,0}$ and $\tau_{\rm s,1}$, respectively)
satisfy $\tau_{\rm s,0}\gg\alpha_{\rm acc}\ \lower
3.0pt\hbox{${\buildrel>\over{\sim}}$}\ \tau_{\rm s,1}$. It is unlikely that
this condition is satisfied at the ice lines of volatile elements other than
H2O. This pile-up process is therefore expected to be effective only for the
water snow line.

###### Acknowledgements.

We thank Vardan Elbakyan for providing his simulation data and Chao-Chin Yang
for helpful comments. SI was supported by MEXT Kakenhi 18H05438. TG was
partially supported by a JSPS Long Term Fellowship at the University of Tokyo.
RH was supported by JSPS Kakenhi JP17J01269 and 18K13600). RH also
acknowledges JAXA’s International Top Young program. SO was supported by JSPS
Kakenhi 19K03926 and 20H01948. ANY was supported by NASA Astrophysics Theory
Grant NNX17AK59G and NSF grant AST-1616929.

## References

  * Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
  * Appelgren et al. (2020) Appelgren, J., Lambrechts, M., & Johansen, A. 2020, A&A, 638, A156
  * Armitage et al. (2013) Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152
  * Blum & Wurm (2000) Blum, J. & Wurm, G. 2000, Icarus, 143, 138
  * Brauer et al. (2008) Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Cho et al. (2003) Cho, J., Lazarian, A., & Vishniac, E. T. 2003, MHD Turbulence: Scaling Laws and Astrophysical Implications, ed. E. Falgarone & T. Passot, Vol. 614, 56–98
  * Ciesla (2010) Ciesla, F. J. 2010, ApJ, 723, 514
  * Ciesla (2011) Ciesla, F. J. 2011, ApJ, 740, 9
  * Ciesla & Cuzzi (2006) Ciesla, F. J. & Cuzzi, J. N. 2006, Icarus, 181, 178
  * Drażkowska & Alibert (2017) Drażkowska, J. & Alibert, Y. 2017, A&A, 608, A92
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Elbakyan et al. (2020) Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V., & Vorobyov, E. I. 2020, A&A, 637, A5
  * Fromang & Papaloizou (2006) Fromang, S. & Papaloizou, J. 2006, A&A, 452, 751
  * Garaud (2007) Garaud, P. 2007, ApJ, 671, 2091
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, arXiv e-prints, arXiv:2001.10000
  * Gonzalez et al. (2017) Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984
  * Gundlach et al. (2018) Gundlach, B., Schmidt, K. P., Kreuzig, C., et al. 2018, MNRAS, 479, 1273
  * Hasegawa et al. (2017) Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31
  * Hyodo et al. (2019) Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90
  * Hyodo et al. (2020) Hyodo, R., Guillot, T., Ida, S., Okuzumi, S., & Youdin, A. N. 2020, submitted (Paper II)
  * Hyodo et al. (2020) Hyodo, R., Ida, S., & Guillot, T. 2020, submitted
  * Ida & Guillot (2016) Ida, S. & Guillot, T. 2016, A&A, 596, L3
  * Ida et al. (2016) Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72
  * Ida & Lin (2008) Ida, S. & Lin, D. N. C. 2008, ApJ, 685, 584
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75
  * Kanagawa et al. (2017) Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142
  * Kimura et al. (2015) Kimura, H., Wada, K., Senshu, H., & Kobayashi, H. 2015, ApJ, 812, 67
  * Kretke & Lin (2007) Kretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55
  * Krijt et al. (2016) Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20
  * Kunitomo et al. (2018) Kunitomo, M., Guillot, T., Ida, S., & Takeuchi, T. 2018, A&A, 618, A132
  * Laibe & Price (2014) Laibe, G. & Price, D. J. 2014, MNRAS, 440, 2136
  * Lambrechts et al. (2014) Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35
  * Lin & Youdin (2017) Lin, M.-K. & Youdin, A. N. 2017, ApJ, 849, 129
  * Morbidelli et al. (2015) Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418
  * Mousis et al. (2019) Mousis, O., Ronnet, T., & Lunine, J. I. 2019, ApJ, 875, 9
  * Musiolik & Wurm (2019) Musiolik, G. & Wurm, G. 2019, ApJ, 873, 58
  * Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
  * Ros & Johansen (2013) Ros, K. & Johansen, A. 2013, A&A, 552, A137
  * Saito & Sirono (2011) Saito, E. & Sirono, S.-i. 2011, ApJ, 728, 20
  * Sato et al. (2016) Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
  * Schoonenberg & Ormel (2017) Schoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21
  * Steinpilz et al. (2019) Steinpilz, T., Teiser, J., & Wurm, G. 2019, ApJ, 874, 60
  * Stevenson & Lunine (1988) Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146
  * Suzuki et al. (2016) Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74
  * Takeuchi & Lin (2005) Takeuchi, T. & Lin, D. N. C. 2005, ApJ, 623, 482
  * Vericel & Gonzalez (2019) Vericel, A. & Gonzalez, J. F. 2019, in SF2A-2019: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. P. Di Matteo, O. Creevey, A. Crida, G. Kordopatis, J. Malzac, J. B. Marquette, M. N’Diaye, & O. Venot, Di
  * Wada et al. (2013) Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62
  * Wada et al. (2011) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Weidling et al. (2012) Weidling, R., Güttler, C., & Blum, J. 2012, Icarus, 218, 688
  * 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
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57
  * Zsom et al. (2011) Zsom, A., Sándor, Z., & Dullemond, C. P. 2011, A&A, 527, A10

##  Appendix A The effects of back-reaction to gas motion

The silicate dust particle and gas radial velocities taking account of the
inertia of the particles (“back-reaction”) at $|z|<H_{\rm d}$ are given
respectively by (Schoonenberg & Ormel 2017)

$\displaystyle\varv_{r}$ $\displaystyle=$
$\displaystyle-\Lambda^{2}\frac{2\tau_{\rm s,d}}{1+\Lambda^{2}\tau_{\rm
s,d}^{2}}\eta\varv_{\rm K}+\Lambda\frac{1}{1+\Lambda^{2}\tau_{\rm
s,d}^{2}}u_{\nu},$ (46) $\displaystyle u_{r}$ $\displaystyle=$ $\displaystyle
Z\Lambda^{2}\frac{2\tau_{\rm s,d}}{1+\Lambda^{2}\tau_{\rm
s,d}^{2}}\eta\varv_{\rm K}+\Lambda\frac{1+\Lambda\tau_{\rm
s,d}^{2}}{1+\Lambda^{2}\tau_{\rm s,d}^{2}}u_{\nu},$ (47)

where $\tau_{\rm s,d}$ is Stokes number of the particles, $Z=\rho_{\rm
d}/\rho_{\rm g}$, $\Lambda=\rho_{\rm g}/(\rho_{\rm g}+\rho_{\rm d})=1/(1+Z)$,
$u_{\nu}$ is an unperturbed disk gas accretion velocity given by

$\displaystyle u_{\nu}\simeq-\frac{3\nu}{2r}\simeq-\frac{3\alpha_{\rm
acc}}{2C_{\eta}}\,\eta\varv_{\rm K},$ (48)

and $C_{\eta}=\eta/(H_{\rm g}/r)^{2}$. When we consider the icy pebbles, the
subscript “d” is replaced by “p.” Equations (46) and (47) are rewritten as

$\displaystyle\varv_{r}$
$\displaystyle\simeq-\frac{\Lambda}{1+\Lambda^{2}\tau_{\rm
s,d}^{2}}\left(2\Lambda\,\tau_{\rm s,d}+\frac{3}{2C_{\eta}}\alpha_{\rm
acc}\right)\eta\varv_{\rm K},$ (49) $\displaystyle u_{r}$
$\displaystyle\simeq-\frac{\Lambda}{1+\Lambda^{2}\tau_{\rm
s,d}^{2}}\left(-2Z\Lambda\,\tau_{\rm
s,d}+\frac{3}{2C_{\eta}}(1+\Lambda\tau_{\rm s,d}^{2})\alpha_{\rm
acc}\right)\eta\varv_{\rm K}.$ (50)

We adopt a two-layer model: a dust-rich midplane layer with scale height
$H_{\rm d}$ and a dust-poor upper layer. To evaluate gas surface density, we
use the vertically averaged $u_{r}$, such that

$\displaystyle Z_{\Sigma}=\frac{\Sigma_{\rm d}}{\Sigma_{\rm
g}}\simeq\frac{\bar{u}_{r}}{\varv_{r}}\frac{f_{\rm d/p}\dot{M}_{\rm
p}}{\dot{M}_{\rm g}}=\frac{\bar{u}_{r}}{\varv_{r}}f_{\rm d/p}F_{\rm p/g},$
(51)

where

$\displaystyle\bar{u}_{r}$ $\displaystyle\simeq\frac{u_{r}|_{z=0}\,H_{\rm
d}+u_{\nu}(H_{\rm g}-H_{\rm d})}{H_{\rm g}}=h_{\rm
d/g}\,u_{r}|_{z=0}+(1-h_{\rm d/g})u_{\nu}.$ (52)

Substituting Eq. (50) into this equation, we obtain

$\displaystyle\bar{u}_{r}$
$\displaystyle\simeq\frac{\Lambda}{1+\Lambda^{2}\tau_{\rm s,d}^{2}}$
$\displaystyle\times\left[\left(\frac{2\tau_{\rm s,d}}{1+Z}+\frac{3\alpha_{\rm
acc}}{2C_{\eta}}\right)Z\,h_{\rm d/g}-\frac{3\alpha_{\rm
acc}}{2C_{\eta}}(1+Z+\Lambda\tau_{\rm s,d}^{2})\right]\eta\varv_{\rm K}.$ (53)

We note that even the vertically averaged gas motion is outward for

$\frac{\alpha_{\rm acc}}{\tau_{\rm s,d}}<\frac{4C_{\eta}}{3}h_{\rm
d/g}\frac{Z}{(1+Z)[1+Z(1-h_{\rm d/g})]},$ (54)

where we assumed $\tau_{\rm s,d}^{2}\ll 1$. For $Z\ \lower
3.0pt\hbox{${\buildrel<\over{\sim}}$}\ 1$ and $C_{\eta}=11/8$, this condition
is reduced to

$\frac{\alpha_{\rm acc}}{\tau_{\rm s,d}}\ \lower
3.0pt\hbox{${\buildrel<\over{\sim}}$}\ \frac{11}{6}h_{\rm d/g}\,Z.$ (55)

In the case of silicate dust particles, since $\tau_{\rm s,d}\ll\alpha_{\rm
acc}$, this condition is not satisfied and the gas flow is always inward.

For pebbles, the outward flow condition is the same as Eq. (55) with “d”
replaced by “p”. The condition is less restrictive than that for silicate
particles, however, it is not satisfied within the parameter range in this
paper, as shown below. We usually consider the cases of $\tau_{\rm
s,p}>\alpha_{\rm acc}$, and $h_{\rm p/g}\simeq(\alpha_{{\rm D}z}/\tau_{\rm
s,p})^{1/2}$. In this case, the modified Eq. (55) is $Z\ \lower
3.0pt\hbox{${\buildrel>\over{\sim}}$}\ (6/11)(\alpha_{\rm acc}/\alpha_{{\rm
D}z})h_{\rm p/g}$. If this equation is satisfied, the (vertically averaged)
gas flow is outward beyond the snow line. However, in our case, $Z$,
$\alpha_{\rm acc}$, and $\alpha_{{\rm D}z}$ are not independent. Equation (62)
is $Z\simeq(1/2)(3/4C_{\eta})(\alpha_{{\rm D}z}/\alpha_{\rm acc})h_{\rm
p/g}F_{\rm p/g}$ for $Z<1$. Thus, the outward flow condition in the icy pebble
region is $F_{\rm p/g}\ \lower 3.0pt\hbox{${\buildrel>\over{\sim}}$}\
2(\alpha_{\rm acc}/\alpha_{{\rm D}z})^{2}$. It is out of the parameter range
that we cover in this paper. In the parameter regimes we cover, $\bar{u}_{r}$
is always negative (inward flow) in the regions of both silicate particles and
icy pebbles.

Substituting Eqs. (53) and (46) into Eq. (51),

$\displaystyle Z\,h_{\rm d/g}$ $\displaystyle\simeq
Z_{\Sigma}\simeq\frac{\bar{u}_{r}}{\varv_{r}}f_{\rm d/p}F_{\rm p/g}$ (56)
$\displaystyle\simeq\frac{-2\frac{Z}{1+Z}\,h_{\rm d/g}\tau_{\rm
s,d}+\frac{3}{2C_{\eta}}\alpha_{\rm acc}\left[1+Z(1-h_{\rm
d/g})+\frac{\tau_{\rm s,d}^{2}}{1+Z}\right]}{2\Lambda\,\tau_{\rm
s,d}+\frac{3}{2C_{\eta}}\alpha_{\rm acc}}f_{\rm d/p}F_{\rm p/g}.$ (57)

Now we adopt the approximation only appropriate for the silicate dust,
$\tau_{\rm s,d}\ll\alpha_{\rm acc}\ll 1$. With this approximation, Eq. (57) is
reduced to

$\displaystyle Z\,h_{\rm d/g}$ $\displaystyle\simeq\left[1+Z(1-h_{\rm
d/g})\right]f_{\rm d/p}F_{\rm p/g},$ (58)

which is solved in terms of $Z$ as

$\displaystyle Z$ $\displaystyle\simeq\frac{f_{\rm d/p}F_{\rm p/g}}{h_{\rm
d/g}-(1-h_{\rm d/g})f_{\rm d/p}F_{\rm p/g}}.$ (59)

As we discussed in Section 2, the pile-up of silicate particles is radially
local and the corresponding local surface density variation of $\Sigma_{\rm
g}$ could be smoothed out. In that case, it may be better to use
$u_{r}=u_{\nu}$ than $\bar{u}_{r}$ here, and we obtain

$\displaystyle Z$ $\displaystyle\simeq\frac{f_{\rm d/p}F_{\rm p/g}}{h_{\rm
d/g}-f_{\rm d/p}F_{\rm p/g}}.$ (60)

Because $h_{\rm d/g}\ll 1$ for the pile-up case, the pileup condition differs
only slightly between Eqs. (59) and (60).

From Eqs. (59) and (60), the runaway pile-up conditions of silicate particles
are given by

$\displaystyle F_{\rm p/g}$ $\displaystyle>\frac{f_{\rm d/p}^{-1}}{1-h_{\rm
d/g}}h_{\rm d/g}\>{\rm or}\>>f_{\rm d/p}^{-1}h_{\rm d/g}$ $\displaystyle[\rm
silicate\>particles],$ (61)

This equation implies that the particle scale height is the key parameter for
the runaway pile-up: the pile-up is favored for a smaller scale height. In
Section 4.3, we discuss the effect of radial and vertical diffusion on the
silicate particle scale height $h_{\rm d/g}$, because it is not in an
equilibrium state.

On the other hand, for pebbles, $\alpha_{\rm acc}\ll\tau_{\rm s,p}\ll 1$ is an
appropriate approximation. From $\alpha_{\rm acc}\ll\tau_{\rm s,p}$, we can
also assume $h_{\rm p/g}\ll 1$. With these approximations, Eq. (57) with “d”
replaced by “p” and $f_{\rm d/p}$ by 1 is reduced to

$\displaystyle Z\,h_{\rm p/g}$
$\displaystyle\simeq\frac{3}{4C_{\eta}}\frac{\alpha_{\rm acc}}{\tau_{\rm
s,p}}(1+Z)^{2}F_{\rm p/g}.$ (62)

This quadratic equation does not have a solution if

$\displaystyle\left(1-\frac{2C_{\eta}}{3}\frac{\tau_{\rm s,p}}{\alpha_{\rm
acc}}\frac{h_{\rm p/g}}{F_{\rm p/g}}\right)^{2}-1<0.$ (63)

Therefore, the runaway pile-up condition for pebbles is

$\displaystyle F_{\rm p/g}$ $\displaystyle>\frac{C_{\eta}}{3}\frac{\tau_{\rm
s,p}}{\alpha_{\rm acc}}h_{\rm p/g}.$ (64)

We will discuss the pebble runaway pile-up condition with the effect of the
vertical stirring by Kelvin-Hermholtz Instability in more details in Paper II.

##  Appendix B Consequences of outward diffusion over the snow line

In Eq. (15), the radial diffusion effect of silicate particles is not
included. The (scaled) radial diffusion flux of silicate particles is

$F_{\rm d,D}\sim 2\pi r\,\alpha_{{\rm D}r}H_{\rm
g}^{2}\Omega\frac{d\Sigma_{\rm d}}{dr}\dot{M}^{-1}_{\rm g}\sim
2\pi\,C_{F}\,r\,\alpha_{{\rm D}r}H_{\rm g}^{2}\Omega\frac{\Sigma_{\rm
d}}{H_{\rm g}}\dot{M}^{-1}_{\rm g},$ (65)

where $C_{F}$ is an unknown numerical factor. We calibrate $C_{F}$ by the
results of Monte Carlo simulations in Figs. 7 and 8 as follows. The (scaled)
net silicate dust mass flux, $F_{\rm d,net}\sim f_{\rm d/p}F_{\rm p/g}-F_{\rm
d,D}$, in a steady state should satisfy $F_{\rm d,net}\sim 3\pi\Sigma_{\rm
d}\alpha_{\rm acc}H_{\rm g}^{2}\Omega\dot{M}^{-1}_{\rm g}$. Substituting Eq.
(65) into this equation, we obtain

$\displaystyle F_{\rm d,D}$ $\displaystyle\sim f_{\rm d/p}{F}_{\rm
p/g}\times\left[1+\frac{3}{2\,C_{F}}\frac{H_{\rm g}}{r}\frac{\alpha_{\rm
acc}}{\alpha_{{\rm D}r}}\right]^{-1}.$ (66)

Accordingly, $(\rho_{\rm d}/\rho_{\rm g})_{0}$ is given by Eq. (15) with
$f_{\rm d/p}F_{\rm p/g}$ replaced by $F_{\rm d,net}=f_{\rm d/p}F_{\rm
p/g}-F_{\rm d,D}$, as

$\displaystyle\left(\frac{\rho_{\rm d}}{\rho_{\rm g}}\right)_{0}=\frac{F_{\rm
d,net}}{h_{\rm d/g}-F_{\rm d,net}},$ (67)

where

$F_{\rm d,net}=f_{\rm d/p}F_{\rm
p/g}\left\\{1-\left[1+\frac{3}{2\,C_{F}}\frac{H_{\rm g}}{r}\frac{\alpha_{\rm
acc}}{\alpha_{{\rm D}r}}\right]^{-1}\right\\}.$ (68)

Figures 7 and 8 show that the analytical results by Eq. (67) with $C_{F}=0.2$
and they fit the numerical results as long as $\rho_{\rm d}/\rho_{\rm g}\
\lower 3.0pt\hbox{${\buildrel<\over{\sim}}$}\ 1$. If the effect of sticking to
icy pebbles by silicate particles that diffuse beyond the snow line is
considered, most of the silicate grains beyond the snow line may eventually
come back with the pebbles, so that $F_{\rm d,net}$ would become more similar
to $f_{\rm d/p}F_{\rm p/g}$.