<figure><img src="image_1.png" /></figure>

HAL
open science

# Planetesimal formation around the snow line

Tristan Guillot, Shigeru Ida, Ryuki Hyodo, Satoshi Okuzumi, Andrew Youdin

► To cite this version:

Tristan Guillot, Shigeru Ida, Ryuki Hyodo, Satoshi Okuzumi, Andrew Youdin. Planetesimal formation around the snow line. Astronomy and Astrophysics - A&A, 2021, 646, pp.A13. 10.1051/0004-6361/202039705. hal-03456292

HAL Id: hal-03456292

https://hal.science/hal-03456292

Submitted on 18 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Planetesimal formation around the snow line

## I. Monte Carlo simulations of silicate dust pile-up in a turbulent disk

Shigeru Ida<fnref n="1" />, Tristan Guillot<fnref n="2" />, Ryuki Hyodo<fnref n="3" />, Satoshi Okuzumi<fnref n="4" />, and Andrew N. Youdin<fnref n="5" />

Received 17 October 2020 / Accepted 24 November 2020

## 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 the 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 developed 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 calculated automatically 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 particle-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 $\lesssim [(\alpha_{Dz}/\alpha_{acc})/3 \times 10^{-2}]^{1/2}$, where $\alpha_{Dz}$ and $\alpha_{acc}$ are the $\alpha$ parameters for vertical turbulent diffusion and disk gas accretion.

**Key words.** planets and satellites: formation – protoplanetary disks

## 1. Introduction

In order for planets to form, micron-sized grains must grow to planetesimal size (i.e., kilometers or more) before they are lost into the central star by gas drag. One possibility is through streaming instabilities (SI), a mechanism to concentrate small solid particles at 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, 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 recondensation, thus increasing the pebble surface density (e.g., Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). Recondensation of outwardly diffused water vapor onto pebbles also induces local pebble growth (Ros & Johansen 2013; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). The

sublimation of icy pebbles changes the ionization rate across the snow line to form a local pressure bump, and the bump could stop the pebble drift there, resulting 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 planetesimals<fnref n="1" />. However, 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 disk gas (Sect. 2.4). Because incoming icy pebbles drift by gas drag quickly, 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 understanding the formation of planetesimals, the composition

<footnote n="1">¹ Earth-Life Science Institute, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8550, Japan
e-mail: ida@elsi.jp</footnote>

<footnote n="2">² Laboratoire J.-L. Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, 06304 Nice, France</footnote>

<footnote n="3">³ ISAS/JAXA, Sagamihara, Kanagawa, Japan</footnote>

<footnote n="4">⁴ Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8551, Japan</footnote>

<footnote n="5">⁵ Steward Observatory/The Lunar and Planetary Laboratory, University of Arizona, Tucson, Arizona 85721, USA</footnote>

<footnote n="1">¹ Our subsequent papers, Hyodo et al. (2021a, referred to as Paper II) and Hyodo et al. (2021b), investigate pile-ups of both icy pebble and silicate dust particles.</footnote>
of planets, and ultimately to accounting for the global ice-to-rock ratio in the Solar System (see Kunitomo et al. 2018).

This pile-up process has been addressed in past literature (Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dążkowska & Alibert 2017; Hyodo et al. 2019). However, they lead to different, seemingly contradictory, results. In order to examine the problem from a new angle, we developed two tools. The 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. We then used 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. The 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 in 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 focused on the 2D Monte Carlo code and its results. In a companion paper (Paper II), we apply the results to the diffusion-advection simulation, allowing us to include the formation of both silicate-rich and water-rich planetesimals. We intend to consider the case of complex protoplanetary disk models in a third paper.

The present article is organized as follows. In Sect. 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 Sect. 3, we describe our Monte Carlo simulation code. In Sect. 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 Sect. 5.

# 2. Previous studies and an analytical solution

## 2.1. Setting

The outskirts of protoplanetary disks are a large reservoir of condensed particles, 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”, that is, 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 relative to the gas mass flux can be significantly above the standard value of ~1/100, which is realized for the mass flux of micron-dust dust particles. For example, Appelgren et al. (2020) find that a decrease in the drift rate of pebbles from the outer to the inner region can lead to a ~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 ~20 to 30. The values of $F_{p/g}$, the ratio of the mass flux of pebbles $\dot{M}_p$ to that of the gas $\dot{M}_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_{p/g}$ that is shaped by the perturbations of the disk, as exemplified in Fig. 1. The high values, much beyond the canonical value (~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.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Histogram of the time-dependent values of $F_{p/g} \equiv \dot{M}_p/\dot{M}_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_{p/g}$ value.</figcaption></figure>

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 sub-disk. However, they had assumed a stationary disk without inward-gas accretion through the disk and radial and vertical diffusion. Our goal is to understand, depending on the value of $F_{p/g}$, whether planetesimals may indeed 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 (Sects. 2.3 and 4.1.2). Because Ida & Guillot (2016) assumed that the scale heights of silicate particles are similar to those of incoming icy pebbles (i.e., much smaller than the disk’s 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. 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 numbers of $10^{-2}-1$, SI can occur for the solid-to-gas surface density ratio
$Z_Σ \gtrsim 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_{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’s gas scale height. As a result, they did not find a runaway pile-up of the silicate particles (for details, see Sects. 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 between 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 Sect. 1, we assume that icy pebbles consist 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 the Stokes number of icy pebbles to be $\tau_{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_{s,d} \ll 1$, where the Stokes number is defined by the stopping time due to gas drag $t_{stop}$ and Keplerian frequency $\Omega_K$ as $\tau_s = t_{stop}\Omega_K$. The notations we use here are summarized in Table 1.

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’s gas accretion rate. Actually, the pebble mass flux through the disk ($\dot{M}_p$) is proportional to the disk’s gas accretion rate ($\dot{M}_g$) in the case of the steady accretion before the pebble formation front reaches the disk’s 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}_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

<table><caption>**Table 1.** Notations.</caption><thead><tr><th>Notation</th><th>Definition</th><th>Eq. #</th></tr></thead><tbody><tr><td>$f_{d/p}$</td><td>Silicate dust mass fraction in a pebble beyond the snow line [= 0.5 (nominal)]</td><td></td></tr><tr><td>$F_{p/g}$</td><td>Pebble-to-gas mass flux in the disk (= $\dot{M}_p/\dot{M}_g$)</td><td>(4)</td></tr><tr><td>$Z_Σ$</td><td>Vertically averaged dust-to-gas mass ratio (= $\Sigma_d/\Sigma_g$)</td><td>(3)</td></tr><tr><td>$Z$</td><td>Local dust-to-gas mass ratio (= $\rho_d/\rho_g = Z_Σ/h_{d/g}$)</td><td>(5)</td></tr><tr><td>$\Lambda$</td><td>Parameter for the solid particle inertia [= $\rho_g/(\rho_d + \rho_g) = 1/(1+Z)$]</td><td>(10)</td></tr><tr><td>$C_\eta$</td><td>Power index of pressure gradient [= $-(1/2)(\partial \ln P/\partial \ln r) = 11/8$ (nominal)]</td><td>(8)</td></tr><tr><td>$\alpha_{acc}$<br>$\alpha_{Dr}$<br>$\alpha_{Dz}$</td><td>α-parameter for disk gas accretion<br>α-parameter for radial mixing<br>α-parameter for vertical stirring</td><td>(9)</td></tr><tr><td>$(\rho_d/\rho_g)_0$</td><td>Maximum value of $\rho_d/\rho_g$ in all the grids (usually, at the grid of $(x, z) \sim (0, 0)$)</td><td></td></tr><tr><td>$h_{d/g}$<br>$h_{d/g,0}$<br>$h_{p/g}$</td><td>Dust-to-gas scale height ratio (= $H_d/H_g$)<br>$h_{d/g}$ at $x$ with $(\rho_d/\rho_g)_0$<br>Pebble-to-gas scale height ratio (= $H_p/H_g$)</td><td>(12)<br>(19)</td></tr><tr><td>$\tau_{s,p}$<br>$\tau_{s,d}$</td><td>Stokes number of pebbles [= 0.1 (nominal)]<br>Stokes number of dust [= $10^{-5}$ (nominal)]</td><td></td></tr></tbody></table>

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.

Firstly, we derived 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’s H-He gas, icy pebbles, and silicate dust particles, respectively. In the steady state, the surface density of the silicate dust particles $\Sigma_d$ and the disk gas $\Sigma_g$ are given by

$$ \Sigma_d = f_{d/p} \dot{M}_p / 2\pi r v_r, \quad (1) $$

$$ \Sigma_g = \dot{M}_g / 2\pi r u_r, \quad (2) $$

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

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

where

$$ F_{p/g} \equiv \frac{\dot{M}_p}{\dot{M}_g}. \quad (4) $$

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

$$ Z \equiv \frac{\rho_d}{\rho_g}, \quad (5) $$

where $\rho_g$ and $\rho_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):

$$v_r = -\Lambda^2 \frac{2\tau_{s,d}}{1 + \Lambda^2 \tau_{s,d}^2} \eta v_K - \Lambda \frac{1}{1 + \Lambda^2 \tau_{s,d}^2} u_v, \quad (6)$$

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

$$\eta \equiv \frac{\Omega_K - \Omega}{\Omega_K} = C_\eta \left(\frac{H_g}{r}\right)^2, \quad (7)$$

$$C_\eta \equiv \frac{1}{2} \frac{\partial \ln P}{\partial \ln r}, \quad (8)$$

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

$$u_v \approx \frac{3v_{\text{acc}}}{2r} \approx \frac{3\alpha_{\text{acc}}H_g^2\Omega}{2r} \approx \frac{3\alpha_{\text{acc}}}{2}\left(\frac{H_g}{r}\right)^2 v_K, \quad (9)$$

where $v_{\text{acc}}$ is the effective viscosity for the angular momentum transfer of the disk gas, $H_g$ is the disk gas scale height, and $u_v$ is defined to be positive for inward flow. The effect of back reaction from the piled-up particles is represented by $\Lambda$, which is defined by

$$\Lambda \equiv \frac{\rho_g}{\rho_g + \rho_d} = \frac{1}{1 + Z}. \quad (10)$$

Because $\rho_d \approx \Sigma_d / \sqrt{2\pi H_d}$, where $H_d$ is the scale height of the silicate dust particles, and $\rho_g \approx \Sigma_g / \sqrt{2\pi H_g}$, Eq. (10) for $z < H_d$ is approximated for the particle motion by

$$\Lambda^{-1} = 1 + \frac{\rho_d}{\rho_g} \approx 1 + Z_{\Sigma} \frac{H_g}{H_d} \approx 1 + h_{d/g}^{-1} Z_{\Sigma}. \quad (11)$$

We define the scale height ratio:

$$h_{d/g} \equiv H_d/H_g. \quad (12)$$

In the case of silicate particle pile-up due to the sublimation of icy pebbles, the pile-up is radially localized just inside of the snow line. If we apply Eq. (A.2), the gas accretion velocity ($u_r$) is locally slowed down from $u_v$ by the back reaction to the gas motion. If the steady accretion with constant $\alpha_{\text{acc}}$ is assumed, the gas surface density given by Eq. (2) is increased locally. However, in reality, the assumption for vertically constant $\alpha_{\text{acc}}$ may be broken down in that case, or some local instability may smooth out the local gas concentration. Therefore, we assumed that $u_r \approx u_v$ 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 below and in Appendix A.

Because $\tau_{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 taken into account. Equation (6) with $\tau_{s,d} \ll 1$ implies $v_r \approx \Lambda u_v$. Substituting $u_r \approx u_v$, $v_r \approx \Lambda u_v$, and Eq. (11) into Eq. (3), we obtain

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

which is solved as (Ida & Guillot 2016)

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

The midplane $Z = \rho_d/\rho_g$ is given by

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

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

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

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

$$F_{p/g} > \frac{h_{d/g}}{f_{d/p}} \quad [\text{w/o the back reaction to gas}], \quad (17)$$

$$F_{p/g} > \frac{h_{d/g}}{(1 - h_{p/g}) f_{d/p}} \quad [\text{w/ the back reaction to gas}]. \quad (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} \approx f_{d/p} F_{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_{p/g} > f_{d/p}^{-1} h_{d/g}$, the estimation of $H_d$ is very important. The scale height of pebbles ($H_p$) is given with their Stokes number, $\tau_{s,p}$, and the vertical mixing parameter $\alpha_{Dz}$ as (Dubrulle et al. 1995; Youdin & Lithwick 2007)

$$h_{p/g} = \frac{H_p}{H_g} = \left(1 + \frac{\tau_{s,p}}{\alpha_{Dz}}\right)^{-1/2}. \quad (19)$$

We discuss the increase of $H_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_d$) must be the same as $H_p$ just after the release from the icy pebbles. Because $\tau_{s,d} \ll \tau_{s,p}$, $H_d$ is increased by the vertical turbulent stirring afterward. Ida & Guillot (2016) assumed that $H_d \approx H_p$ ($h_{d/g} = h_{p/g}$) until the runaway pile-up develops. For $h_{d/g} \sim 0.1$ ($\tau_{s,peb}/\alpha_{Dz} \sim 10^2$) and $f_{d/p} \sim 0.5$, the optimistic estimate by Eq. (15) predicts that the runaway pile-up occurs for $F_{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_d \approx H_g$ ($h_{d/g} \approx 1$). Because they neglected the back reaction to $v_r$, they never found the runaway pile-up in their simulations. Even if we apply Eq. (15), which was derived with the back reaction, for the case of $h_{d/g} \approx 1$ and $f_{d/p} \sim 0.5$, the runaway pile-up condition is $F_{p/g} > 2$. The high value of $F_{p/g}$ may not be completely ruled out, but it is not easily 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_{Dz}$ to derive a model for $h_{d/g}$ as a function
of $x = r - r_{\text{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

$$D_r = \alpha_{Dr} H_g^2 \Omega \Lambda^K, \qquad (20)$$

$$D_z = \alpha_{Dz} H_g^2 \Omega \Lambda^K. \qquad (21)$$

A simple argument that $D_z, D_r \propto H_g^2 \propto c_s^2 \propto P/(\rho_g + \rho_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_g + \rho_d) v_l^2 / t_{\text{eddy}}$, is independent of $l$, where $l, v_l$, and $t_{\text{eddy}}$ are the eddy size, velocity, and turnover time $\sim l/v_l$. Because $D_z, D_r \sim \ell_0 v_l \ell_0$ where $\ell_0$ and $v_l$ are the eddy size and velocity controlling the turbulence, $D_z, D_r \propto \ell_0 [\ell_0/(\rho_g + \rho_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 $v_r$, and that the condition for the runaway pile-up depends on the relations among $\alpha_{\text{acc}}$, $\alpha_{Dr}$, and $\alpha_{Dz}$, while Ida & Guillot (2016) considered the limits of $\alpha_{Dr}, \alpha_{Dz} \to 0$ and their condition, Eq. (15), does not include $\alpha_{\text{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_{p/g} > 0.2-0.3$, and $\alpha_{Dr} = \alpha_{Dz} \approx 10^{-3}$ in the case of $K = 1$ and $2$ with $\alpha_{\text{acc}} = 10^{-2}$, while it does not occur with $K = 0$ as long as it is in the parameter space of $\alpha_{Dr} = \alpha_{Dz} = 10^{-3}-10^{-2}$ and $F_{p/g} < 0.6$. In this paper, we 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_v \sim (3/2)\alpha_{\text{acc}}(H_g/r)^2 v_K$ (Eq. (9)) and the drift speed relative to the gas is given by $v_r \sim \tau_{s,d}(H_g/r)^2 v_K$ (Eq. (6)), the total drift velocity is dominated by $u_v$ if $\tau_{s,d} \lesssim \alpha_{\text{acc}}$. Therefore, the assumption here is expressed as $\tau_{s,d} \lesssim \alpha_{\text{acc}} < \tau_{s,p}$.

In Sect. 2.3, we assume that $\tau_{s,p} \sim 0.1$, $\alpha_{\text{acc}} \sim 10^{-3}-10^{-2}$, and $\tau_{s,d} \lesssim \alpha_{\text{acc}}$. Saito & Sirono (2011) assumed that many µm-sized silicate particles are embedded in pebbles. Morbidelli et al. (2015) considered chondrule-sized (~mm) particles. Because the Stokes number is proportional to the particle size in the Epstein regime and its square in the Stokes regime, it is reasonable to assume that $\tau_{s,d} < \alpha_{\text{acc}}$. This model corresponds to the “many-seeds model” of Schoonenberg & Ormel (2017). In the Monte Carlo simulations here, we adopt $\tau_{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_{Dr}^2 + \alpha_{Dz}^2)^{1/2} \tau_{s,d})^{1/2} c_s$, or by the drift velocity difference between the particles $\sim 2\tau_{s,d} \eta v_K$ (e.g., Sato et al. 2016). The collision velocity between silicate particles is therefore $(\tau_{s,d}/\tau_{s,p})^{-(0.5-1)}$ times that between icy pebbles. If the threshold collision velocity is ten times lower for the silicate particles than for the icy particles (Blum & Wurm 2000; Zsom et al. 2010, 2011; Wada et al. 2011, 2013; Weidling et al. 2012) and $\tau_{p,d} \sim 0.1$

is close to the fragmentation/rebound barrier for ice, the silicate particles can grow only up to $\tau_{s,d} \sim (10^{-1}-10^{-2})\tau_{s,p} \sim 10^{-3}-10^{-2}$ in size. Therefore, for $\alpha_{\text{acc}} = 10^{-3}-10^{-2}$, the drift velocity of the silicate particles is comparable to $u_v$ and the pile-up would not be significantly affected by the value of $\tau_{s,d}$.

We note that this conventional view has been challenged recently (Kimura et al. 2015; Gundlach et al. 2018; Musiolik & Wurm 2019; Steinvilz 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_{s,d}$ to find that the results are not affected for $\tau_{s,d} \lesssim \alpha_{\text{acc}}$ and that the runaway pile-up region is removed from lower $F_{p/g}$ regions as $\tau_{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 was explored by Vericel & Gonzalez (2019) (see also Gonzalez et al. 2017). In this case, planetesimal formation may still be possible, but through a slowdown 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 developed 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 $v_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_{\text{snow}}$ is the radial distance from the snow line at $r_{\text{snow}}$ and $z$ is the distance from the midplane.

We considered a radially local region near the snow line and neglected the $r$-dependence of the Keplerian frequency ($\Omega$) and the disk’s gas scale height ($H_g$). As discussed in Sect. 2.4, we adopted the conventional view to set $\tau_{s,p} = 0.1$ and $\tau_{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 injected a new silicate super-particle with mass $m$ with a radially uniform distribution near the snow line in the $-0.5\Delta x_{\text{subl}} < x < 0.5\Delta x_{\text{subl}}$ range, and a Gaussian distribution of $z$ of the root mean square $\Delta z_{\text{subl}}$ as

$$x = 0.5 \mathcal{R} \Delta x_{\text{subl}}, \qquad (22)$$

$$z = \sqrt{2} \text{erf}^{-1}(|\mathcal{R}|) \frac{\mathcal{R}}{|\mathcal{R}|} \Delta z_{\text{subl}}, \qquad (23)$$

where $\mathcal{R}$ is a random number in the range of $[-1, 1]$, $\Delta x_{\text{subl}}$ is the characteristic sublimation radial length, and $\Delta z_{\text{subl}}$ corresponds to the scale height of the incoming icy pebbles, $H_p$, given by Eq. (19). We adopted $\Delta x_{\text{subl}} = 0.1 H_g$ as a nominal case. In Sect. 4.4.2, we also consider 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_{\text{acc}}$, $\alpha_{Dr}$, and $\alpha_{Dz}$, respectively. According to different values of $\alpha_{Dz}$, we used the consistent value of the scale height of the injected silicate particles given by Eq. (19). Although Hasegawa et al. (2017) suggested $\alpha_{Dr}, \alpha_{Dz} \sim 0.1 \alpha_{\text{acc}}$, the relations among
$α_{Dr}$, $α_{Dz}$ and $α_{acc}$ are not clear (e.g., Armitage et al. 2013). Therefore, we surveyed broad parameter ranges of $α_{Dr}$, $α_{Dz}$ and $α_{acc}$; we mostly show the results with $α_{acc} ≥ α_{Dr} = α_{Dz}$. In some cases, the results with $α_{Dr} ≠ α_{Dz}$ ($α_{acc} ≥ α_{Dr}$, $α_{Dz}$) are also shown.

## 3.2. Advection and diffusion

At each time step $δt$, we changed the r-z locations of the particles that were injected before, following Ciesla (2010, 2011), by

$$δx = v_{adv,r}δt + \mathcal{R}(6D_rδt)^{1/2}, \qquad (24)$$

$$δz = v_{adv,z}δt + \mathcal{R}(6D_zδt)^{1/2}, \qquad (25)$$

where the first and second terms on the right-hand side represent advection (drift) and diffusion, respectively. The root mean square of $\mathcal{R}$ is $1/\sqrt{3}$. In our simulations, we adopted $δt = Ω^{-1}$. Fromang & Papaloizou (2006) suggested a turbulent correlation time in a protoplanetary disk is ~0.15$T_K$ ~ $Ω^{-1}$. Ciesla (2011) adopted planar (x-y) random walks, in addition to a vertical one (z), to take account of the global curvature effect. In our case, because we only considered a local radial range, we simply adopted r-z random walks.

We performed 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 the back reaction to diffusion is considered (i.e., $K > 0$). Hereafter, we show the results with $K = 1$ as representative of the cases including a back reaction to diffusion.

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

$$v_{adv,z} = -(α_{Dz} + τ_{s,d})Ωz. \qquad (26)$$

The first term ($\propto α_{Dz}$) comes from the effect of diffusion smoothing the concentration ($ρ_d/ρ_g$) but not $ρ_d$ itself.

## 3.3. Inclusion of back reaction effects

The advection and diffusion scheme was originally developed by Ciesla (2010, 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 utilized 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 injected one super-particle at every time step $δt = Ω^{-1}$. For given $\dot{M}_p$, we can determine the individual mass of the super-particles as

$$m = f_{d/p} \dot{M}_p δt = f_{d/p} \dot{M}_p Ω^{-1}, \qquad (27)$$

where $f_{d/p}$ is a silicate fraction in migrating pebbles ($f_{d/p} = 0.5$ in the nominal case). To locally average the silicate dust mass density, we used linear grids, $Δx$ and $Δz$. In the nominal cases, we adopted $Δx = 0.1 H_g$. We used $Δz ∼ Δz_{subl}$ near the midplane and larger $Δ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 + Δx]$ and $[|z|, |z| + Δz]$ is

$$ρ_d = \frac{m ΔN_{x,z}}{2πrΔx × 2Δz} = \frac{H_g}{r} \frac{ΔN_{x,z} f_{d/p} \dot{M}_p}{4π(Δx/H_g)(Δz/H_g) H_g^3 Ω}, \qquad (28)$$

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

$$ρ_g = \frac{\Sigma_g}{\sqrt{2πH_g}} \exp\left(-\frac{z^2}{2H_g^2}\right) = \frac{\dot{M}_g}{\sqrt{2π} 3πα_{acc}H_g^3 Ω} \exp\left(-\frac{z^2}{2H_g^2}\right). \qquad (29)$$

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

$$Z = \frac{ρ_d}{ρ_g} = \frac{3\sqrt{2π} H_g}{4 r} α_{acc} \frac{f_{d/p} \dot{M}_p}{\dot{M}_g} \frac{ΔN_{x,z}}{(Δx/H_g)(Δz/H_g)} \exp\left(\frac{z^2}{2H_g^2}\right) = 3.75 × 10^{-2} \left(\frac{H_g/r}{0.04}\right) \left(\frac{α_{acc}}{10^{-2}}\right) \left(\frac{f_{d/p}}{0.5}\right) F_{p/g} \times \left(\frac{Δx}{0.1H_g}\right)^{-1} \left(\frac{Δz}{0.1H_g}\right)^{-1} ΔN_{x,z} \exp\left(\frac{z^2}{2H_g^2}\right). \qquad (30)$$

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

In the simulations, we normalized length by $H_g$, including the range of particle injection, $Δx_{subl}$ and $Δz_{subl}$, and the grid sizes, $Δx$ and $Δz$. The time was scaled by $Ω^{-1}$, that is, each super-particle was injected at every unit's scaled time. Therefore, our results can be applied for any location of $r_{snow}$, as long as $H_g/r$ is the same.

In the simulations here, we adopted $f_{d/p} = 0.5$ and $H_g/r = 0.04$ at $r = r_{snow}$, as nominal values. We do not need to specify the values of $H_g$ and $\dot{M}_g$ or the value of $r_{snow}$, because only the scaled values, $h_{d/g} (= H_d/H_g)$, $h_{p/g} (= H_p/H_g)$, $H_g/r$ and $F_{p/g} (= \dot{M}_p/\dot{M}_g)$, are used in the simulations. For the viscous $α$-parameters, we show that only the ratios $α_{Dz}/α_{acc}$ and $α_{Dr}/α_{acc}$ are important.

# 4. Results

## 4.1. Code check

We tested 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 also offer a 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 nonzero $α_{Dr}$ and $α_{Dz}$, they assumed the evolution of $H_d$ analytically estimated only from the vertical stirring with $α_{Dz}$. In our simulation, the evolution of $H_d$ was self-consistently calculated, and we find that the radial mixing with $α_{Dr}$ is also important for $H_d$. On the other hand, the sublimation width $Δx_{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 shows the distribution of particles and their scale heights (the local root mean square of $z$ of the particles) in a steady state at $t = 10^4 Ω^{-1}$ (after $10^4$ particles are injected),
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Snapshots of 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_g$). We plot $|z|$, while $z$ takes either positive or negative values. We set $\alpha_{acc} = 10^{-2}$, $\alpha_{Dr} = 0$ and $\tau_s = 10^{-3}$. The vertical mixing parameter in each panel is (a) $\alpha_{Dz} = 10^{-5}$, (b) $\alpha_{Dz} = 10^{-3}$, and (c) $\alpha_{Dz} = 10^{-1}$. The analytical estimates of the equilibrium particle scale height for each $\alpha_{Dz}$ are shown by the dashed green lines. The back reactions are not included for $v_r$ ($\Lambda = 1$) and $\alpha_{Dz}$ ($K = 0$). We injected the particles in the range of $x = [-0.05, 0.05]H_g$ and $z = [-0.1, 0.1]H_g$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Time evolution of the maximum $\rho_d/\rho_g$. The back reaction to $v_r$ is included in the results in panel b, but not in panel a. The light blue, magenta, blue, green, and red curves are the results with $F_{p/g} = 0.025, 0.05, 0.1, 0.2, 0.4$, and $0.8$, respectively. The other parameters are fixed as $\alpha_{acc} = 10^{-2}$ and $\Delta x_{subl} = 0.1 H_g$ and $\Delta z_{subl} = 0.03 H_g$. To mimic the settings of Ida & Guillot (2016), we artificially set $\alpha_{Dr} = 0$, $\alpha_{Dz} = 10^{-5}$ and $\tau_{s,d} = 10^{-2}$ (however, the relatively large $\tau_{s,d}$ for the silicate dust is not reflected to the dust drift speed). The latter two parameters maintain $H_d = \Delta z_{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_{p/g} = 0.025$) given by Eq. (15).</figcaption></figure>

obtained by our Monte Carlo simulations for (a) $\alpha_{Dz} = 10^{-5}$, (b) $10^{-3}$, and (c) $10^{-1}$. For these runs, we neglected the back reactions to $v_r$ ($\Lambda = 1$) and $\alpha_{Dz}$ ($K = 0$). To highlight the effect of vertical stirring, we set $\alpha_{Dr} = 0$. The other parameters are the particles’ Stokes numbers $\tau_s = 10^{-3}$ and $\alpha_{acc} = 10^{-2}$. The equilibrium scale height predicted from Eq. (19) with $\tau_{s,p}$ replaced by $\tau_s = 10^{-3}$ are (a) $0.1 H_g$, (b) $0.71 H_g$, and (c) $1.0 H_g$, respectively, which are shown via 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 the results of Schoonenberg & Ormel (2017) and Ida & Guillot (2016)

Ida & Guillot (2016) adopted the following setup: (1) the turbulent diffusion is neglected for simplicity, and (2) the silicate particle pile-up is fast enough that $H_d$ is not increased from $\Delta z_{subl} (= F_p)$ by the vertical stirring. To mimic the setup of Ida & Guillot (2016), we set the following in the Monte Carlo simulations: (1) we adopted $\alpha_{Dr} = 0$, $K = 0$, and $\Delta z_{subl} = 0.03 H_g$;

(2) to keep $H_d = \Delta z_{subl} = 0.03 H_g$ as an equilibrium, we adopted $\alpha_{Dz} = 10^{-5}$ and $\tau_{s,d} = 10^{-2}$ (Eq. (19)); (3) the dust particle drift speed is set to be $v_r = -\Lambda u_v$, where the artificially enlarged $\tau_{s,d}$ is not reflected (Eq. (6)), and (4) for gas, $u_r = -u_v$ is always assumed. For comparison, $v_r = -u_v$ (without the back reaction to $v_r$) is also examined.

In the rest of the paper, we always set $\tau_{s,p} = 10^{-1}$ and $\tau_{s,d} = 10^{-5}$, and we consistently calculated $H_p(= \Delta z_{subl})$ by Eq. (19) and the evolution of $H_d$ from $\Delta z_{subl}$ by the vertical and radial diffusion, with given $\alpha_{Dz}$ and $\alpha_{Dr}$. The above artificial setting is only used in the particular runs in Fig. 3.

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

$$ \frac{\rho_d}{\rho_g} \approx f_{d/p} h_{d/g}^{-1} F_{p/g} \approx 16.7 F_{p/g}. \qquad (31) $$

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

In Fig. 3b, the back reaction to $v_r$ is included ($\Lambda < 0$), corresponding to the setup of Ida & Guillot (2016). They predicted that runaway pile-up occurs for $F_{p/g} > h_{d/g} f_{d/p}^{-1}$, which is $F_{p/g} > 0.06$ for the parameters adopted here. This panel shows a clear runaway pile-up for $F_{p/g} \gtrsim 0.1$, which is consistent with the prediction.

A comparison between Figs. 3a and b clearly shows that the back reaction to $v_r$ of the silicate particles plays an essential role
<figure><img src="image_5.png" /><figcaption>**Fig. 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_{subl}$ replaced by $\max(\Delta x_{subl}, |x|)$ (see the discussion at the end of Sect. 4.3.2). The magenta curve represents $\log_{10}(\rho_d/\rho_g)$ near the midplane. We set $\alpha_{acc} = 10^{-2}$, $\tau_{s,d} = 10^{-5}$, and $f_{d/p} = 0.5$, and $\dot{M}_g = 10^{-8} M_{\odot} \text{yr}^{-1}$. The nine panels adopt different values of $F_{p/g}$ and $\alpha_{Dz} = \alpha_{Dr}$, as indicated in each panel.</figcaption></figure>

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 $v_r$ exists, we include it in the results in the rest of the paper.

If the supply rate of silicate particles ($f_{d/p} \dot{M}_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_d/\rho_g$, because the back reaction becomes effective once $\rho_d/\rho_g$ exceeds unity. While $\dot{M}_p$, more exactly $F_{p/g}$, determines $Z_\Sigma = \Sigma_d/\Sigma_g$, the local value of $Z = \rho_d/\rho_g$ is regulated by $H_d$ and $Z_\Sigma$ ($Z = Z_\Sigma/h_{d/g}$) near the snow line. The silicate dust particle scale height $H_d$ is regulated by the pebble scale height $H_p$ ($= \Delta z_{subl}$ in our simulation) and turbulent mixing parameters, $\alpha_{Dz}$ and $\alpha_{Dr}$, as shown in Sect. 4.3.

## 4.2. Typical results

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

Figure 4 shows 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_{Dz} = \alpha_{Dr} = 10^{-2}$ ($\alpha_{Dz}/\tau_{s,d} = 10^3$), the vertical mixing is fast enough to realize the upper limit of $H_d$ ($H_d \approx H_g$) for all of $F_{p/g} = 0.1, 0.3$, and 1. Even for a high-mass flux of pebbles ($F_{p/g} = 1.0$), no pile-up of the silicate particles is found for $\alpha_{Dz} = \alpha_{Dr} = 10^{-2}$. In the case of $\alpha_{Dz} = \alpha_{Dr} = 10^{-6}$, on the other hand, the mixing is so weak that $H_d$ is not increased to the equilibrium value [$\approx (\alpha_{Dr}/\tau_{s,d})^{1/2} H_g \sim 0.3 H_g$] in the range of $x > -2 H_g$. Nevertheless, $H_d$ near the injection point ($x \sim 0$) is much higher than $\Delta z_{subl} = H_p = 0.003 H_g$. This means that $H_d$ is regulated by $\alpha_{Dr}$ as well as by $\alpha_{Dz}$, which is discussed in Sect. 4.3.

Figure 5 shows the results with the back reaction to diffusion coefficients ($K = 1$). In general, the runaway pile-up is more pronounced for smaller $\alpha_{Dz} = \alpha_{Dr}$ and larger $F_{p/g}$. The magenta curves represent $\log_{10}(\rho_d/\rho_g)$ near the midplane. The midplane $\rho_d/\rho_g$ usually takes the maximum value, at the injection (sublimation) region, $x \sim 0$, which is denoted by $(\rho_d/\rho_g)_0$. In the panels with $(\rho_d/\rho_g)_0 \lesssim 1$ (the panels with $\alpha_{Dz} = 10^{-2}$ and those with $\alpha_{Dz} = 10^{-4}$ and $F_{p/g} = 0.3$ and $0.1$), the results are similar to those in Fig. 4. If $\rho_d/\rho_g \gtrsim$ several in Fig. 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 shows the time evolution of $(\rho_d/\rho_g)_0$ and $H_d/H_p$ at the maximum $\rho_d/\rho_g$ locations ($h_{d/g,0}$), for different $F_{p/g}$. In panel a, the back reaction to diffusion is not included ($K = 0$). As Figs. 4 and 5 suggest, the transient $\alpha_{Dz}$ ($= \alpha_{Dr}$) values between the runaway pile-up and non-pile-up cases are $\alpha_{Dz} \approx 10^{-5}$ for $K = 0$ and $\alpha_{Dz} \approx 10^{-4}$ for $K = 1$, so we plot the results with $\alpha_{Dz} = 10^{-5}$ in panel a and $\alpha_{Dz} = 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_d$ evolution. The comparison with the result of Ida & Guillot (2016) in Sect. 2 suggests that $H_d$ is directly related to a threshold value of $F_{p/g}$ for the runaway pile-up: a lower $H_d$ leads to a lower threshold value of $F_{p/g}$. However, Fig. 6 shows that $H_d$ increases as $\rho_d/\rho_g$ increases in the case of $K = 0$, while it decreases for $K = 1$. Therefore, the pile-up suffers negative feedback for $K = 0$, but it suffers positive feedback for $K = 1$, resulting in a much clearer transition to the runaway pile-up.

In both cases, the back reaction to $v_r$ is included. As the pile-up proceeds, $v_r$ decreases, which enhances the pile-up. However, we also included the vertical stirring. As $v_r$ decreases, the silicate particles are stirred up for a longer period of time before they leave the sublimation region, resulting in higher $H_d$.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 4, except for K = 1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Time evolution of $(\rho_d/\rho_g)_0$ and the dust scale height at the maximum $\rho_d/\rho_g$ locations scaled by the initial one (= $H_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 $v_r$ is included in both panels a and b. The light blue, magenta, blue, green, and red curves are the results with $F_{p/g} = 0.05, 0.1, 0.2, 0.4$ and $0.8$, respectively. The other parameters are fixed as $\alpha_{acc} = 10^{-2}$, $\Delta x_{subl} = 0.1$, and $\Delta z_{subl} = 0.03$. The diffusion parameters are $\alpha_{Dz} = \alpha_{Dr} = 10^{-5}$ in panel a and $\alpha_{Dz} = \alpha_{Dr} = 10^{-4}$ in panel b.</figcaption></figure>

near that location. In panel a, $H_d$ near the sublimation region quickly asymptotes to an equilibrium value for $F_{p/g} \le 0.2$ (lower $(\rho_d/\rho_g)_0$ cases), and it increases with time for $F_{p/g} = 0.8$ (the

runaway pile-up cases). In panel b, in the runaway pile-up cases ($F_{p/g} = 0.4$ and $0.8$), $H_d$ initially increases with time. However, after $\rho_d/\rho_g$ exceeds ~1, the diffusion coefficients become so small that the settling of the particles overwhelms the vertical stirring, and $H_d$ decreases. Thus, the different responses to the pile-up between $K = 0$ and $K = 1$ result in different $H_d$ and $\rho_d/\rho_g$ evolution in the pile-up cases. As we discussed in Sect. 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 shows $(\rho_d/\rho_g)_0$ and $h_{d/g,0}$ at $t = 15\,000\,\Omega^{-1}$ as a function of $\alpha_{Dz} = \alpha_{Dr}$ for different $F_{p/g}$. This figure more clearly shows that the asymptotic equilibrium values of $h_{d/g,0}$ are independent of $F_{p/g}$ in the non-runaway pile-up cases, which is suggested by Fig. 6. The envelope curves in the plots for (a) $K = 0$ and (b) $K = 1$ match the analytical formula derived in the next sub-subsection (Eqs. (B.3) and (35) with Eq. (36)). The equilibrium values are generally larger than the scale height of the injected silicate particles, $\Delta z_{subl}$, due to the effects of vertical and radial diffusion. In panels a and b, the values of $h_{d/g,0}$ that deviate from the envelope curves correspond to the runaway pile-up; $h_{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 Fig. 7 are analytical formulas for the non pile-up cases, which are derived as follows. The timescale to drift by $\Delta x_{subl}$ is

$$ t_{\text{drift}} \approx \frac{\Delta x_{\text{subl}}}{v_r} \\ \approx \frac{\Delta x_{\text{subl}}}{(3/2)\Lambda\alpha_{\text{acc}}(H_g/r)^2 v_K} \approx \frac{\Delta x_{\text{subl}}/H_g}{(3/2)\Lambda\alpha_{\text{acc}}(H_g/r)\Omega}. \quad (32) $$
<figure><img src="image_8.png" /><figcaption>Fig. 7. Maximum $\rho_d/\rho_g$ and $H_d/H_g$ at the maximum $\rho_d/\rho_g$ locations at $t = 15000 \, \Omega^{-1}$ as a function of $\alpha_{Dz} = \alpha_{Dr}$ for different $F_{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_{acc} = 10^{-2}$, $\Delta x_{subl} = 0.1 H_g$, and $\Delta z_{subl} H_p$ given by Eq. (19). Upper and bottom panels: the red dashed curves are the analytical formulas given by Eqs. (B.3) and (35) with Eq. (36). Lower panels: the green dashed lines represent $\Delta z_{subl}/H_g$.</figcaption></figure>

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

$$z_{\text{diff}} \approx \sqrt{D_z t_{\text{drift}}} \approx \left( \frac{2}{3} \alpha_{D,z} \Lambda^K H_g^2 \Omega \times \frac{\Delta x_{\text{subl}}/H_g}{\Lambda \alpha_{\text{acc}}(H_g/r)\Omega} \right)^{1/2} \\ \approx \left( \frac{2}{3} \Lambda^{K-1} \frac{\Delta x_{\text{subl}}/H_g}{H_g/r} \frac{\alpha_{D,z}}{\alpha_{\text{acc}}} \right)^{1/2} H_g. \quad (33)$$

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

If only vertical diffusion is considered, the silicate particle scale height scaled by $H_g$ is predicted as

$$h_{d/g} = \frac{H_d}{H_g} \approx \frac{\max(\Delta z_{\text{subl}}, z_{\text{diff}})}{H_g} \approx \frac{(\Delta z_{\text{subl}}^2 + z_{\text{diff}}^2)^{1/2}}{H_g} \\ \approx \left( \frac{\Delta z_{\text{subl}}^2}{H_g^2} + \frac{2}{3} \frac{\alpha_{Dz}}{\alpha_{\text{acc}}} \frac{\Delta x_{\text{subl}}/H_g}{H_g/r} \right)^{1/2}. \quad (34)$$

Furthermore, during $t_{drift}$, the particles are radially mixed in the range of $x_{\text{diff}} \approx \sqrt{D_r t_{\text{drift}}} \sim (\alpha_{Dr}/\alpha_{Dz})^{1/2} z_{\text{diff}}$. In general, the dust scale height increases with distance from the snow line. If $x_{\text{diff}} \gg \Delta x_{\text{subl}}$, the radial diffusion brings the high-$z$ particles back to the injection region and raises the value of $h_{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_{d/g}$ by a factor of $[1 + (x_{\text{diff}}/\Delta x_{\text{subl}})^2]$. When $\alpha_{Dr}$ is comparable to $\alpha_{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 replaced $\alpha_{Dr}/\alpha_{acc}$ in $(x_{\text{diff}}/\Delta x_{\text{subl}})^2$ by $(\alpha_{Dr}/\alpha_{acc})/[1 + (C_{r,\text{diff}} \alpha_{Dr}/\alpha_{acc})^2]$ in the above equation. Comparing the formula with the numerical results (Figs. 7 and 8), we empirically adopted $C_{r,\text{diff}} = 10$. With these modifications,

$$h_{d/g,*} \approx \left( \frac{\Delta z_{\text{subl}}^2}{H_g^2} + \frac{2}{3} \frac{\alpha_{Dz}}{\alpha_{\text{acc}}} \frac{\Delta x_{\text{subl}}/H_g}{H_g/r} \right)^{1/2} \\ \times \left[ 1 + \frac{2}{3} \frac{(\alpha_{Dr}/\alpha_{\text{acc}})}{1 + (C_{r,\text{diff}} \alpha_{Dr}/\alpha_{\text{acc}})^2} \frac{1}{(H_g/r)(\Delta x_{\text{subl}}/H_g)} \right]. \quad (35)$$

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

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

As the lower panels in Fig. 7 show, the analytical formula reproduces the envelope curves of $H_d/H_p$ obtained by the numerical simulations for the non-pile-up cases. We note that $h_{d/g,0}$ and $(\rho_d/\rho_g)_0$ depend on the ratios of $\alpha_{Dr}/\alpha_{acc}$ and $\alpha_{Dz}/\alpha_{acc}$, but not on their absolute values.

Equation (35) is the analytical formula for $h_{d/g}$ at $x \sim 0$. This formula can be extrapolated to any $x$ by replacing $\Delta x_{\text{subl}}$ with $\max(\Delta x_{\text{subl}}, |x|, \epsilon)$, where $\epsilon$ is added to avoid the divergence in the case of extremely small $\Delta x_{\text{subl}}$; we adopted $\epsilon = 0.01 H_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_d/\rho_g)_0$ (the mid-plane particle-to-gas ratio at $x \sim 0$) for the non-runaway pile-up cases and for the phase before the runaway pile-up proceeds. If the estimated $(\rho_d/\rho_g)_0$ exceeds $\sim 1$, the actual value of $\rho_d/\rho_g$ should deviate from the estimate and increase with time.

Basically, $(\rho_d/\rho_g)_0$ is predicted by Eq. (15) with the scale height formula, Eqs. (36) and (35). A possible detailed correction is to include the effect of outward diffusion flux ($F_{d,D}$) of silicate particles beyond the snow line that can be found in the snapshots with $\alpha_{Dr} = \alpha_{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 achieved with $f_{d/p}F_{p/g}$ being replaced by $F_{d,net} = f_{d/p}F_{p/g} - F_{d,D}$, where

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

Accordingly, Eq. (15) is reduced to

$$\left(\frac{\rho_d}{\rho_g}\right)_0 = \frac{F_{d,net}}{h_{d/g} - F_{d,net}}. \quad (38)$$

This formula reproduces the numerical results in Figs. 7 and 8, as long as $\rho_d/\rho_g \lesssim 1$. We note, however, that if the effect of silicate
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 7, except for $\alpha_{Dz}/\alpha_{Dr}$ relations and the value of $\Delta x_{subl}$. The maximum $\rho_d/\rho_g$ and $h_{d/g}$ (= $H_d/H_g$) for different $F_{p/g}$ are given as functions of (a) $\alpha_{Dr} = \alpha_{Dz}$ with $\Delta x_{subl} = 0.1 H_g$, which is the same as Fig. 7b, for comparison, (b) $\alpha_{Dr} = \alpha_{Dz}$ with $\Delta x_{subl} = H_g$, (c) $\alpha_{Dz}$ with $\alpha_{Dr} = 10^{-4}$ and $\Delta x_{subl} = 0.1 H_g$, and (d) $\alpha_{Dr}$ with $\alpha_{Dz} = 10^{-4}$ and $\Delta x_{subl} = 0.1 H_g$. For all the cases, $K = 1$ is adopted. Upper and bottom panels: the red dashed curves are the analytical formulas given by Eqs. (35) and (B.3) with corresponding parameter values.</figcaption></figure>

particles that diffuse beyond the snow line sticking to icy pebbles is considered, most of the silicate grains beyond the snow line may eventually come back with the pebbles, so $F_{d,net}$ would become more similar to $f_{d/p}F_{p/g}$.

We 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 the intrinsic physics of the pile-up of the silicate particles released from the drifting icy pebble by their sublimation clear, 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_{Dz}/\alpha_{acc}$ and $\alpha_{Dr}/\alpha_{acc}$ for given $\Delta z_{subl}$, $\Delta x_{subl}$, and $H_g/r$. We also carried out the Monte Carlo simulations with $\alpha_{acc} = 10^{-3}$ and confirmed that the plots in Fig. 7 do not change as long as we use $\alpha_{Dr}/\alpha_{acc}$ (= $\alpha_{Dz}/\alpha_{acc}$) as the horizontal coordinate.

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

So far, the results with $\alpha_{Dr} = \alpha_{Dz}$ have been shown. In Fig. 8c, $\alpha_{Dr}$ is fixed and $\alpha_{Dz}$ is independently changed. In Fig. 8d, $\alpha_{Dz}$ is fixed. The envelope curves for $H_d/H_g$ (= $h_{d/g}$) are also reproduced by the analytical formula 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

The analytical formulas reproduce the numerical results for a broad range of multidimensional parameters. Figure 9a shows the contour maps of the maximum $\rho_d/\rho_g$ on the plane of $\alpha_{Dz}/\alpha_{acc}$ (= $\alpha_{Dr}/\alpha_{acc}$), and $F_{p/g}$ at $t = 1.5 \times 10^4 \Omega^{-1}$ for $\Delta x_{subl} = 0.1 H_g$. In Fig. 9b, $\Delta x_{subl}$ is given by Eq. (41). In both cases, $\alpha_{acc} = 10^{-2}$ and $K = 1$. The left and right panels show 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_{acc} = 10^{-2}$ and $\alpha_{Dz}$ (= $\alpha_{Dr}$) are in the range of $[10^{-6}, 10^{-2}]$. The analytical formula to predict $(\rho_d/\rho_g)_0$ is Eq. (37) with Eqs. (38), (35), and (36). They are functions of $\alpha_{Dr}/\alpha_{acc}$ and $\alpha_{Dz}/\alpha_{acc}$ for given $\Delta x_{subl}$, as long as the term $\Delta z_{subl}^2$ in Eq. (35) is negligible, which is satisfied in the parameter range we considered, as shown in Figs. 7 and 8. We also carried out simulations with $\alpha_{acc} = 10^{-3}$ and $\alpha_{Dz}$ (= $\alpha_{Dr}$) 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_{Dz}/\alpha_{acc}$ (= $\alpha_{Dr}/\alpha_{acc}$).

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

$$ F_{p/g} \approx f_{d/p}^{-1} h_{d/g,*} \\ \approx f_{d/p}^{-1} \left( \frac{2}{3} \frac{\alpha_{Dz}}{\alpha_{acc}} \frac{\Delta x_{subl}/H_g}{H_g/r} \right)^{1/2} \left[ 1 + \frac{2}{3} \frac{\alpha_{Dr}}{\alpha_{acc}} \frac{1}{(H_g/r)(\Delta x_{subl}/H_g)} \right]. \quad (39) $$

Adopting $H_g/r = 0.04$ and $f_{d/p} = 1/2$ as nominal values,

$$ F_{p/g} \approx 0.13 \left( \frac{f_{d/p}}{1/2} \right)^{-1} \left( \frac{\alpha_{Dz}/\alpha_{acc}}{10^{-2}} \right)^{1/2} \left( \frac{H_g/r}{0.04} \right)^{1/2} \left( \frac{\Delta x_{subl}}{0.1 H_g} \right)^{1/2} \\ \times \left[ 1 + \frac{5}{3} \left( \frac{\alpha_{Dr}/\alpha_{acc}}{10^{-2}} \right) \left( \frac{H_g/r}{0.04} \right)^{-1} \left( \frac{\Delta x_{subl}}{0.1 H_g} \right)^{-1} \right] \\ \propto \left\{ \begin{array}{l} \left( \frac{\alpha_{Dz}}{\alpha_{acc}} \right)^{1/2} \left( \frac{\alpha_{Dr}}{\alpha_{acc}} \right) \Delta x_{subl} \left[ \frac{\alpha_{Dz}}{\alpha_{acc}} > 0.6 \times 10^{-2} \right], \\ \\ \left( \frac{\alpha_{Dz}}{\alpha_{acc}} \right)^{1/2} \Delta x_{subl}^{1/2} \left[ \frac{\alpha_{Dz}}{\alpha_{acc}} < 0.6 \times 10^{-2} \right]. \end{array} \right. \quad (40) $$

This equation explicitly shows that the pile-up condition is scaled by $\alpha_{Dr}/\alpha_{acc}$ and $\alpha_{Dz}/\alpha_{acc}$.

In Eq. (40), the sublimation width, $\Delta x_{subl}$, is a free parameter. We used $\Delta x = 0.1 H_g$ as a nominal parameter. The sublimation width was 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
<figure><img src="image_10.png" /><figcaption>Fig. 9. Contour maps of the maximum $ρ_d/ρ_g$ on the plane of $α_{Dz}/α_{acc}(=α_{Dr}/α_{acc})$ and $F_{p/g}$ with (a) $Δx_{subl} = 0.1 H_g$, and (b) $Δx_{subl}$ given by Eq. (41). In both cases, $α_{acc} = 10^{-2}$ and $K = 1$ are adopted. Left panels: results of Monte Carlo simulations at $t = 15 000 Ω^{-1}$, and right ones: results after using analytical formulas given by Eq. (B.3) with Eqs. (35), (36), (B.4), and the corresponding $Δx_{subl}$. The Monte Carlo simulations were performed with a grid size of $Δlog_{10}(α_{Dz}/α_{acc}) = Δlog_{10} F_{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}(ρ_d/ρ_g)$.</figcaption></figure>

sublimation width for $α_{Dz} = α_{Dr}$ and $α_{acc} = 10^{-2}$ is roughly approximated as

$$ \log_{10}\left(\frac{\Delta x_{subl}}{H_g}\right) \approx \frac{X_+ - X_-}{2}\operatorname{erf}\left[3 \log_{10}\left(\frac{\alpha_{Dr}/\alpha_{acc}}{0.08}\right)\right] + \frac{X_+ + X_-}{2}, $$

$$ X_- = \log_{10}(\Delta x_{subl}/H_g)_- = \log_{10} 2; $$

$$ X_+ = \log_{10}(\Delta x_{subl}/H_g)_+ = \log_{10} 0.1, \qquad (41) $$

where $(\Delta x_{subl}/H_g)_- = 2$ and $(\Delta x_{subl}/H_g)_+ = 0.1$ are the lower $α_{Dr}/α_{acc}$ limit at which advection dominates, and the higher one at which radial diffusion dominates, respectively (Paper II). The error function smoothly connects the two limits. Figures 9a and b show that the threshold value of $F_{p/g}$ is several times higher than in Fig. 9a for $α_{Dz}/α_{acc} \sim 10^{-4}$, while it is similar for $α_{Dz}/α_{acc} \sim 10^{-2}$. In the former case, $F_{p/g} \propto Δx_{subl}^{1/2}$ (Eq. (40)), so $F_{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. (41)), while the dependence on $Δx_{subl}$ must be slightly weaker in the latter case (Eq. (41)). Thereby, the silicate particle runaway pile-up condition with Eq. (41) is roughly given by

$$ F_{p/g} \approx \left( \frac{\alpha_{Dz}/\alpha_{acc}}{3 \times 10^{-2}} \right)^{1/2}, \qquad (42) $$

which agrees with the results in Fig. 9.

Figure 9 shows the results with $α_{Dr} = α_{Dz}$. As mentioned in Sect. 1, magneto-hydrodynamic simulations gave $α_{Dr} < α_{Dz}$ (Zhu et al. 2015; Yang et al. 2018). We also performed simulations with $α_{Dr} = 0.1 α_{Dz}$. As we can predict from Eq. (40), the runaway pile-up region is extended to a value several times larger of $α_{Dz}/α_{acc}$ at $F_{p/g} \sim 1$ in the case of $Δx_{subl} = 0.1 H_g$, while the runaway pile-up region is almost the same for smaller values of $F_{p/g}$. Because $α_{Dz}/α_{acc}$ is very uncertain, we surveyed a rather broad parameter space and leave the estimate of $α_{Dz}/α_{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 $α_{acc} = 10^{-2}$ and $α_{Dr} = α_{Dz} = 10^{-3}, 3 \times$

$10^{-3}, 10^{-2}$ to find that the silicate runaway pile-up occurs for $α_{Dz} = α_{Dr} \approx 10^{-1}α_{acc}$ and $F_{p/g} \ge 0.3$ in the case of $K = 1$. However, Fig. 9b shows that the silicate runaway pile-up region is restricted to a smaller parameter range, that is, to values of $α_{Dz}/α_{acc} = α_{Dr}/α_{acc}$ that are one order of magnitude smaller for a given $F_{p/g}$ value. As already mentioned, this is because Hyodo et al. (2019) assumed that the silicate scale height is $\sim Δz_{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-sized sublimation width. As a result, smaller $α_{Dz}$ and $α_{Dr}$ are required for a runaway 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 effective 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 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 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_{acc}$, $\alpha_{D_r}$, and $\alpha_{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 estimate of Hyodo et al. (2019), due to the coupled effect of radial and vertical diffusion, and that for the silicate particle runaway pile-up to occur, $\alpha_{D_z}/\alpha_{acc} \lesssim 10^{-2} \times (\text{pebble-to-gas mass flux})^2$ is required, which is more restrictive than the result of Hyodo et al. (2019) ($\alpha_{D_z}/\alpha_{acc} \lesssim 10^{-1}$ and $F_{p/g} \gtrsim 0.3$). To clarify if the condition is actually satisfied, detailed nonideal MHD simulations are needed to evaluate $\alpha_{D_z}$, $\alpha_{D_r}$, and $\alpha_{acc}$.

Thus far, we have not included the pile-up of icy pebbles. This would occur upstream, that is, beyond the snow line, and could thus suppress the supply of pebbles and therefore 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 investigate these issues in Paper II.

One would want to extend the arguments developed in this work to other ice lines, such as for NH₃ and CO₂. The pile-up process proposed here only occurs if the Stokes numbers of the particles before and after the sublimation of a volatile component ($\tau_{s,0}$ and $\tau_{s,1}$, respectively) satisfy $\tau_{s,0} \gg \alpha_{acc} \gtrsim \tau_{s,1}$. It is unlikely that this condition is satisfied at the ice lines of volatile elements other than H₂O. This pile-up process is therefore only expected to be effective for the water snow line.

*Acknowledgements.* We thank Vardan Elbakyan for providing his simulation data and Chao-Chin Yang for helpful comments. S.I. was supported by MEXT Kakenhi 18H05438. T.G. was partially supported by a JSPS Long Term Fellowship at the University of Tokyo. R.H. was supported by JSPS Kakenhi JP17J01269 and 18K13600. R.H. also acknowledges JAXA's International Top Young program. S.O. was supported by JSPS Kakenhi 19K03926 and 20H01948. A.N.Y. was supported by NASA Astrophysics Theory Grant NNX17AK59G and NSF grant AST-1616929.

## References

Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phy., 56, 1756

Appelgren, J., Lambrechts, M., & Johansen, A. 2020, A&A, 638, A156

Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14

Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152

Blum, J., & Wurm, G. 2000, Icarus, 143, 138

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43

Cho, J., Lazarian, A., & Vishniac, E. T. 2003, MHD Turbulence: Scaling Laws and Astrophysical Implications, eds. E. Falgarone & T. Passot, Lect. Notes Phys., 614, 56

Ciesla, F. J. 2010, ApJ, 723, 514

Ciesla, F. J. 2011, ApJ, 740, 9

Ciesla, F. J., & Cuzzi, J. N. 2006, Icarus, 181, 178

Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237

Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V., & Vorobyov, E. I. 2020, A&A, 637, A5

Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751

Garaud, P. 2007, ApJ, 671, 2091

Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132

Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Gundlach, B., Schmidt, K. P., Kreuzig, C., et al. 2018, MNRAS, 479, 1273

Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31

Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90

Hyodo, R., Guillot, T., Ida, S., Okuzumi, S., & Youdin, A. N. 2021a, A&A, 646, A14 (Paper II)

Hyodo, R., Ida, S., & Guillot, T. 2021b, A&A, 645, L9

Ida, S., & Guillot, T. 2016, A&A, 596, L3

Ida, S., & Lin, D. N. C. 2008, ApJ, 685, 584

Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72

Johansen, A., Youdin, A., & Mac Low M.-M. 2009, ApJ, 704, L75

Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142

Kimura, H., Wada, K., Senshu, H., & Kobayashi, H. 2015, ApJ, 812, 67

Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55

Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20

Kunitomo, M., Guillot, T., Ida, S., & Takeuchi, T. 2018, A&A, 618, A132

Laibe, G., & Price, D. J. 2014, MNRAS, 440, 2136

Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35

Lin, M.-K., & Youdin, A. N. 2017, ApJ, 849, 129

Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418

Mousis, O., Ronnet, T., & Lunine, J. I. 2019, ApJ, 875, 9

Musiolik, G., & Wurm, G. 2019, ApJ, 873, 58

Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106

Ros, K., & Johansen, A. 2013, A&A, 552, A137

Saito, E., & Sirono, S.-i. 2011, ApJ, 728, 20

Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15

Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21

Steinpilz, T., Teiser, J., & Wurm, G. 2019, ApJ, 874, 60

Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74

Vericel, A., & Gonzalez, J. F. 2019, SF2A-2019: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. P. Di Matteo, O. Creevey, A. Crida, et al.

Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36

Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Weidling, R., Güttler, C., & Blum, J. 2012, Icarus, 218, 688

Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80

Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588

Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81

Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57

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_d$ are given, respectively, by (Schoonenberg & Ormel 2017)

$$v_r = -\Lambda^2 \frac{2\tau_{s,d}}{1 + \Lambda^2 \tau_{s,d}^2} \eta v_K + \Lambda \frac{1}{1 + \Lambda^2 \tau_{s,d}^2} u_v, \quad (A.1)$$

$$u_r = Z\Lambda^2 \frac{2\tau_{s,d}}{1 + \Lambda^2 \tau_{s,d}^2} \eta v_K + \Lambda \frac{1 + \Lambda \tau_{s,d}^2}{1 + \Lambda^2 \tau_{s,d}^2} u_v, \quad (A.2)$$

where $\tau_{s,d}$ is the Stokes number of the particles, $Z = \rho_d/\rho_g$, $\Lambda = \rho_g/(\rho_g + \rho_d) = 1/(1+Z)$, $u_v$ is an unperturbed disk gas accretion velocity given by

$$u_v \approx -\frac{3v}{2r} \approx -\frac{3\alpha_{\text{acc}}}{2C_{\eta}} \eta v_K, \quad (A.3)$$

and $C_\eta = \eta/(H_g/r)^2$. When we consider the icy pebbles, the subscript “d” is replaced by “p.” Equations (A.1) and (A.2) are rewritten as

$$v_r \approx -\frac{\Lambda}{1 + \Lambda^2 \tau_{s,d}^2} \left( 2\Lambda \tau_{s,d} + \frac{3}{2C_{\eta}} \alpha_{\text{acc}} \right) \eta v_K, \quad (A.4)$$

$$u_r \approx -\frac{\Lambda}{1 + \Lambda^2 \tau_{s,d}^2} \left( -2Z\Lambda \tau_{s,d} + \frac{3}{2C_{\eta}}(1 + \Lambda \tau_{s,d}^2)\alpha_{\text{acc}} \right) \eta v_K. \quad (A.5)$$

We adopted a two-layer model: a dust-rich midplane layer with scale height $H_d$ and a dust-poor upper layer. To evaluate gas surface density, we used the vertically averaged $u_r$, such that

$$Z_{\Sigma} = \frac{\Sigma_d}{\Sigma_g} \approx \frac{\bar{u}_r}{v_r} \frac{f_{d/p} \dot{M}_p}{\dot{M}_g} = \frac{\bar{u}_r}{v_r} f_{d/p} F_{p/g}, \quad (A.6)$$

where

$$\bar{u}_r \approx \frac{u_r|_{z=0} H_d + u_v(H_g - H_d)}{H_g} = h_{d/g} u_r|_{z=0} + (1 - h_{d/g}) u_v. \quad (A.7)$$

Substituting Eq. (A.5) into this equation, we obtain

$$\bar{u}_r \approx \frac{\Lambda}{1 + \Lambda^2 \tau_{s,d}^2} \times \left[ \left( \frac{2\tau_{s,d}}{1+Z} + \frac{3\alpha_{\text{acc}}}{2C_{\eta}} \right) Z h_{d/g} - \frac{3\alpha_{\text{acc}}}{2C_{\eta}} (1+Z + \Lambda \tau_{s,d}^2) \right] \eta v_K. \quad (A.8)$$

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

$$\frac{\alpha_{\text{acc}}}{\tau_{s,d}} < \frac{4C_{\eta}}{3} h_{d/g} \frac{Z}{(1+Z)[1+Z(1-h_{d/g})]}, \quad (A.9)$$

where we assumed $\tau_{s,d}^2 \ll 1$. For $Z \lesssim 1$ and $C_\eta = 11/8$, this condition is reduced to

$$\frac{\alpha_{\text{acc}}}{\tau_{s,d}} \lesssim \frac{11}{6} h_{d/g} Z. \quad (A.10)$$

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

For pebbles, the outward flow condition is the same as Eq. (A.10) 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_{s,p} > \alpha_{\text{acc}}$ and $h_{p/g} \approx (\alpha_{Dz}/\tau_{s,p})^{1/2}$. In this case, the modified Eq. (A.10) is $Z \gtrsim (6/11)(\alpha_{\text{acc}}/\alpha_{Dz})h_{p/g}$. If this equation is satisfied, the (vertically averaged) gas flow is outward beyond the snow line. However, in our case, $Z$, $\alpha_{\text{acc}}$, and $\alpha_{Dz}$ are not independent. Equation (A.16) is $Z \approx (1/2)(3/4C_\eta)(\alpha_{Dz}/\alpha_{\text{acc}})h_{p/g}F_{p/g}$ for $Z < 1$. Thus, the outward flow condition in the icy pebble region is $F_{p/g} \gtrsim 2(\alpha_{\text{acc}}/\alpha_{Dz})^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. (A.8) and (A.1) into Eq. (A.6),

$$\begin{aligned} Z h_{d/g} \approx Z_{\Sigma} &\approx \frac{\bar{u}_r}{v_r} f_{d/p} F_{p/g} \\ &\approx \frac{-2 \frac{Z}{1+Z} h_{d/g} \tau_{s,d} + \frac{3}{2C_{\eta}} \alpha_{\text{acc}} \left[ 1 + Z(1 - h_{d/g}) + \frac{\tau_{s,d}^2}{1+Z} \right]}{2\Lambda \tau_{s,d} + \frac{3}{2C_{\eta}} \alpha_{\text{acc}}} f_{d/p} F_{p/g}. \end{aligned} \quad (A.11)$$

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

$$Z h_{d/g} \approx [1 + Z(1 - h_{d/g})] f_{d/p} F_{p/g}, \quad (A.12)$$

which is solved in terms of Z as

$$Z \approx \frac{f_{d/p} F_{p/g}}{h_{d/g} - (1 - h_{d/g}) f_{d/p} F_{p/g}}. \quad (A.13)$$

As we discuss in Sect. 2, the pile-up of silicate particles is radially local and the corresponding local surface density variation of $\Sigma_g$ could be smoothed out. In that case, it may be better to use $u_r = u_v$ than $\bar{u}_r$ here, and we obtain

$$Z \approx \frac{f_{d/p} F_{p/g}}{h_{d/g} - f_{d/p} F_{p/g}}. \quad (A.14)$$

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

From Eqs. (A.13) and (A.14), the runaway pile-up conditions of silicate particles are given by

$$F_{p/g} > \frac{f_{d/p}^{-1}}{1 - h_{d/g}} h_{d/g} \quad \text{or} \quad > f_{d/p}^{-1} h_{d/g} \quad \text{[silicate particles].} \quad (A.15)$$

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 Sect. 4.3, we discuss the effect of radial and vertical diffusion on the silicate particle scale height, $h_{d/g}$, because it is not in an equilibrium state.

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

$$Z h_{p/g} \approx \frac{3}{4C_{\eta}} \frac{\alpha_{\text{acc}}}{\tau_{s,p}} (1 + Z)^2 F_{p/g}. \quad (A.16)$$
This quadratic equation does not have a solution if

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

Therefore, the runaway pile-up condition for pebbles is

$$ F_{p/g} > \frac{C_{\eta}}{3} \frac{\tau_{s,p}}{\alpha_{\text{acc}}} h_{p/g}. \qquad (A.18) $$

We 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_{d,D} \sim 2\pi r \alpha_{Dr} H_g^2 \Omega \frac{d\Sigma_d}{dr} \dot{M}_g^{-1} \sim 2\pi C_F r \alpha_{Dr} H_g^2 \Omega \frac{\Sigma_d}{H_g} \dot{M}_g^{-1}, \quad (B.1) $$

where $C_F$ is an unknown numerical factor. We calibrated $C_F$ by the results of Monte Carlo simulations in Figs. 7 and 8 as follows. The (scaled) net silicate dust mass flux, $F_{d,net} \sim f_{d/p}F_{p/g}$ -

$F_{d,D}$, in a steady state, should satisfy $F_{d,net} \sim 3\pi\Sigma_d\alpha_{acc}H_g^2\Omega\dot{M}_g^{-1}$. Substituting Eq. (B.1) into this equation, we obtain

$$ F_{d,D} \sim f_{d/p} F_{p/g} \times \left[ 1 + \frac{3}{2 C_F} \frac{H_g}{r} \frac{\alpha_{\text{acc}}}{\alpha_{Dr}} \right]^{-1}. \qquad (B.2) $$

Accordingly, the midplane particle-to-gas ratio at $x \sim 0$, $(\rho_d/\rho_g)_0$, is given by Eq. (15) with $f_{d/p}F_{p/g}$ replaced by $F_{d,net} = f_{d/p}F_{p/g} - F_{d,D}$, as

$$ \left(\frac{\rho_d}{\rho_g}\right)_0 = \frac{F_{d,net}}{h_{d/g} - F_{d,net}}, \qquad (B.3) $$

where

$$ F_{d,net} = f_{d/p} F_{p/g} \left\{ 1 - \left[ 1 + \frac{3}{2 C_F} \frac{H_g}{r} \frac{\alpha_{\text{acc}}}{\alpha_{Dr}} \right]^{-1} \right\}. \qquad (B.4) $$

We plotted the analytical estimate of $(\rho_d/\rho_g)_0$ given by Eq. (B.3) with $C_F = 0.2$ in Figs. 7 and 8. They fit the numerical results as long as $(\rho_d/\rho_g)_0 \lesssim 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 $F_{d,net}$ would become more similar to $f_{d/p}F_{p/g}$.