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

HAL
open science

# Planetesimal formation around the snow line

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

► To cite this version:

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

HAL Id: hal-03456293

https://hal.science/hal-03456293

Submitted on 19 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
## II. Dust or pebbles?

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

Received 12 November 2020 / Accepted 11 December 2020

ABSTRACT

*Context.* Forming planetesimals is a major challenge in our current understanding of planet formation. Around the snow line, icy pebbles and silicate dust may locally pile up and form icy and rocky planetesimals via a streaming instability and/or gravitational instability. The scale heights of both pebbles and silicate dust released from sublimating pebbles are critical parameters that regulate the midplane concentrations of solids.

*Aims.* Here, using a realistic description of the scale height of silicate dust and that of pebbles, we wish to understand disk conditions for which a local runaway pile-up of solids (silicate dust or icy pebbles) occurs inside or outside the snow line.

*Methods.* We performed 1D diffusion-advection simulations that include the back-reaction (the inertia) to radial drift and diffusion of icy pebbles and silicate dust, ice sublimation, the release of silicate dust, and their recycling through the recondensation and sticking onto pebbles outside the snow line. We used a realistic description of the scale height of silicate dust obtained from a companion paper and that of pebbles including the effects of a Kelvin-Helmholtz instability. We study the dependence of solid pile-up on distinct effective viscous parameters for turbulent diffusions in the radial and vertical directions ($α_{Dr}$ and $α_{Dz}$) and for the gas accretion to the star ($α_{acc}$) as well as that on the pebble-to-gas mass flux ($F_{p/g}$).

*Results.* Using both analytical and numerical approaches, we derive the sublimation width of drifting icy pebbles which is a critical parameter to characterize the pile-up of silicate dust and pebbles around the snow line. We identify a parameter space (in the $F_{p/g} - α_{acc} - α_{Dz} (= α_{Dr})$ space) where pebbles no longer drift inward to reach the snow line due to the back-reaction that slows down the radial velocity of pebbles (we call this the “no-drift” region). We show that the pile-up of solids around the snow line occurs in a broader range of parameters for $α_{acc} = 10^{-3}$ than for $α_{acc} = 10^{-2}$. Above a critical $F_{p/g}$ value, the runaway pile-up of silicate dust inside the snow line is favored for $α_{Dr}/α_{acc} \ll 1$, while that of pebbles outside the snow line is favored for $α_{Dr}/α_{acc} \sim 1$. Our results imply that a distinct evolutionary path in the $α_{acc} - α_{Dr} - α_{Dz} - F_{p/g}$ space could produce a diversity of outcomes in terms of planetesimal formation around the snow line.

**Key words.** accretion, accretion disks – planets and satellites: formation – planet-disk interactions – protoplanetary disks

**1. Introduction**

Forming planetesimals in protoplanetary disks is a major challenge in our current understanding of planet formation. Streaming instability (SI; e.g., Youdin & Goodman 2005) and gravitational instability (GI; e.g., Goldreich & Ward 1973) could be prominent candidate mechanisms to directly form planetesimals from small particles, although the detailed conditions and applicability are still a matter of debate.

The water snow line may be a favorable location where solids selectively pile up (Fig. 1). During the passage of icy pebbles drifting through the water snow line, icy pebbles sublimate and silicate dust is ejected<fnref n="1" /> (Saito & Sirono 2011; Morbidelli et al. 2015). The formation of icy planetesimals by the streaming instability outside the snow line could be triggered via the local

enhancement of pebble spatial density through the recondensation of diffused water vapor and sticking of silicate dust onto pebbles beyond the snow line (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019). Inside the snow line, rocky planetesimals might be preferentially formed by the gravitational instability of piled up silicate dust because the drift velocity suddenly drops there as the size significantly decreases from sublimating pebbles (Ida & Guillot 2016; Hyodo et al. 2019).

Recent 1D numerical simulations aimed to study the pile-up of solids around the water snow line (Saito & Sirono 2011; Morbidelli et al. 2015; Estrada et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Charnoz et al. 2019; Hyodo et al. 2019; Gárate et al. 2020). Although these different studies considered similar settings of pebbles and silicate dust, they neglected some of the following essential physical processes to regulate the midplane solid-to-gas ratio: these are, for example, (1) gas-dust friction for solids (pebbles and/or silicate dust), that is, the back-reaction to their drift velocity and diffusive motion in the radial and vertical

<footnote n="1">¹ ISAS/JAXA, Sagamihara, Kanagawa, Japan
e-mail: hyodo@elsi.jp</footnote>

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

<footnote n="3">³ Earth-Life Science Institute, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8550, 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">¹ A critical collision velocity for rebound and fragmentation is conventionally expected to be ~10 times lower for silicate than for ice (e.g., Blum & Wurm 2000; Wada et al. 2011), although recent studies showed updated sticking properties of silicate (e.g., Kimura et al. 2015; Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019).</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Summary of solid pile-up around the snow line. Icy pebbles, which uniformly contain micron-sized silicate dust formed at the outer disk, drift inward due to the gas drag. Pebbles sublimate through the passage of the snow line (Sect. 3) and silicate dust is released together with the water vapor. Silicate dust is well coupled to the gas and the radial drift velocity significantly drops from that of pebbles. This causes so-called “traffic-jam” effect and the dust piles up inside the snow line. A fraction of water vapor and silicate dust diffuses from inside to outside the snow line recycles onto pebbles via the recondensation and sticking, causing a local pile-up of pebbles just outside the snow line. Midplane concentrations (i.e., midplane solid-to-gas ratio) of pebbles and silicate dust strongly depend on their scale heights. Near the snow line, the scale height of silicate dust is the minimum as the dust is released from pebbles whose scale height is much smaller than that of the gas. This causes the maximum midplane concentration of silicate dust just inside the snow line. Significant midplane concentrations of silicate dust and/or icy pebbles would cause gravitational instability and/or streaming instability, forming rocky and/or icy planetesimals, respectively (Sect. 5).</figcaption></figure>

directions, and (2) a realistic prescription of the scale height of silicate dust as well as that of pebbles.

Back-reactions that slow down drift velocities of pebbles and silicate dust as their pile-up proceeds (hereafter, Drift-BKR) are essential for pile-up as $Σ_{p(ord)} ∝ 1/v_{p(ord)}$ (a subscript of “p (or d)” describes physical parameters of either pebbles or those of silicate dust), where $Σ_{p(ord)}$ and $v_{p(ord)}$ are the surface density and the radial velocity of pebbles or silicate dust, respectively (Ida & Guillot 2016). Back-reaction to the diffusive motion (hereafter, Diff-BKR) can trigger runaway pile-up of solids as their diffusivities are progressively weakened as pile-up proceeds (Hyodo et al. 2019, see Eqs. (21) and (22) below).

The scale height of silicate dust $H_d$ is another critical consideration on the solid pile-up. The midplane spatial density of silicate dust is $\propto 1/H_d$ and thus the midplane solid-to-gas ratio is directly affected by $H_d$. The efficiency of the recycling (sticking) of diffused silicate dust onto pebbles outside the snow line is $\propto 1/H_d$ as the number density of colliding silicate dust is $\propto 1/H_d$. All of the previous studies of 1D calculations used a too simple prescription for $H_d$ and their results overestimated or underestimated the pile-up of solids, depending on their assumed $H_d$ and the disk conditions (see more details in Appendix A).

In our companion paper (Ida et al. 2021; hereafter Paper I), we performed 2D (radial-vertical) Monte Carlo simulations (a Lagrange method) of silicate dust that was released around the snow line from drifting icy pebbles under ice sublimation. Based on these Monte Carlo simulations and using semi-analytical arguments, Paper I derived a formula for the evolving scale height of the silicate dust for a given disk parameters as a function of the distance to the snow line (Eq. (36) below). The derived scale height of silicate dust is a function of the radial width where the dust is released from sublimating pebbles (i.e., sublimation width of pebbles; see Eq. (36)). However, Paper I did not include the physical processes of the sublimation. Here, our 1D code is suitable for including the pressure-dependent sublimation around the snow line and we study the sublimation width of the drifting pebbles (Sect. 3).

Paper I also studied conditions for the pile-up of silicate dust just inside the snow line and identified disk parameters where runaway pile-up of silicate dust occurs (i.e., rocky planetesimal

formation). However, Paper I did not include recycling of water vapor and of silicate dust outside the snow line due to their diffusive motion and hence they did not consider pile-up process of pebbles outside the snow line.

The scale height of pebbles is a function of the Stokes number (larger pebbles tend to reside in a vertically thinner disk midplane layer). However, for a sufficiently thin pebble-dominated midplane layer may interact with the upper/lower gas-dominated layer, inducing a Kelvin-Helmholtz (KH) instability (Sekiya 1998; Chiang 2008), which was not considered before. Here, we also consider the effects of a KH instability for a more realistic picture of the scale height of pebbles and their midplane concentration (Sect. 4.1).

In reality, pile-ups of silicate dust and/or pebbles around the snow line are the consequence of combinations of ice sublimation, dust release, and their recycling onto pebbles together with their complex radial and vertical motions (Fig. 1). In this work, we perform new 1D simulations that include both Drift-BKR and Diff-BKR (the same as Hyodo et al. 2019) and that adopt, for the first time, a realistic scale height of silicate dust (derived in Paper I) and that of pebbles (derived in Sect. 4.1). We discuss favorable conditions for runaway pile-ups of pebbles and/or silicate dust by investigating a much broader disk parameter space than previous studies.

In Sect. 2, we explain our numerical methods and models. In Sect. 3, using 1D numerical simulations and analytical arguments, we derive sublimation width of drifting icy pebbles, which is a critical parameter to regulate a local solid pile-up around the snow line. In Sect. 4, we derive the scale height of pebbles considering the effects of a KH instability, and we also derive a critical $F_{p/g}$ above which pebble drift is stopped due to the back-reaction onto the gas that slows down the radial velocity of pebbles. In Sect. 5, we show our overall results of 1D simulations that combined with a realistic scale heights of silicate dust and pebbles. In Sect. 6, we summarize our paper.

## 2. Numerical methods and models

Here, we describe numerical methods and models. In this paper, we use three distinct non-dimensional parameters to describe the
gas-solid evolutions in protoplanetary disk: $\alpha_{acc}$ which describes the global efficiency of angular momentum transport of gas (i.e., corresponds to the gas accretion rate) via turbulent viscosity $v_{acc}$, $\alpha_{Dr}$ which describes the radial turbulence strength (i.e., corresponds to the radial velocity dispersion), and $\alpha_{Dz}$ which describes the vertical turbulence strength (i.e., corresponds to the vertical velocity dispersion) and which is used to describe the scale height of pebbles<fnref n="2" />.

The classical $\alpha$-disk model (Shakura & Sunyaev 1973) has adopted fully turbulent disks, i.e., $\alpha_{acc} \approx \alpha_{Dr} \approx \alpha_{Dz}$, whereas recent theoretical developments suggest that the vertically averaged efficiency of the disk angular momentum transport, i.e., $\alpha_{acc}$, could be different from local radial and vertical turbulence strengths $\alpha_{Dr}$ and $\alpha_{Dz}$ (e.g., Armitage et al. 2013; Liu et al. 2019). The gas accretion can be dominated by a process other than radial turbulent diffusion such as angular momentum transport by large-scale magnetic fields (e.g., Bai et al. 2016). Magnetohydrodynamic (MHD) simulations have shown that a region of ionization that is too low for the magneto-rotational instability (MRI) to operate (“dead zone”; Gammie 1996) might ubiquitously exist in the inner part of disk midplane and only surface layers are magnetically active ($\alpha_{Dr}, \alpha_{Dz} = 10^{-5}-10^{-3}$ within a dead zone; e.g., Gressel et al. 2015; Simon et al. 2015; Bai & Stone 2013; Bai et al. 2016; Mori et al. 2017). These results indicate $\alpha_{Dr}, \alpha_{Dz} < \alpha_{acc}$ (see also Zhu et al. 2015; Hasegawa et al. 2017; Yang et al. 2018) and we study a wide parameter range where $\alpha_{Dr} (= \alpha_{Dz}) \le \alpha_{acc}$.

## 2.1. Structure of the gas disk

The surface density of the gas at a radial distance $r$ is expressed as a function of the gas accretion rate $\dot{M}_g$ and the gas radial velocity $v_g$ (i.e., the classical $\alpha$-accretion disk model; Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974) as

$$ \Sigma_g = \frac{\dot{M}_g}{2\pi r v_g}. \qquad (1) $$

The isothermal sound speed of the gas is written as

$$ c_s = \sqrt{\frac{k_B T}{\mu_g m_{proton}}}, \qquad (2) $$

where $k_B$ is the Boltzmann constant, $\mu_g$ is the mean molecular weight of the gas ($\mu_g = 2.34$), $m_{proton}$ is the proton mass, and $T$ is the temperature of the disk gas. In this paper, the temperature profile of the disk is fixed as

$$ T(r) = T^* \left( \frac{r}{3.0 \text{ au}} \right)^{-\beta}, \qquad (3) $$

where we use $T^* = 150$ K and $\beta = 1/2$. Then, the scale height of the gas is given as

$$ H_g = \frac{c_s}{\Omega_K} \approx 0.033 \text{ au} \times \left( \frac{r}{1.0 \text{ au}} \right)^{5/4}, \qquad (4) $$

where $\Omega_K$ is the Keplerian orbital frequency. Disk gas pressure at the midplane $P_g$ is given as

$$ P_g = \frac{\Sigma_g}{\sqrt{2\pi H_g} \mu_g m_{proton}}. \qquad (5) $$

The gas rotates at sub-Keplerian speed. The degree of deviation of the gas rotation frequency from that of Keplerian $\eta$ is given by

$$ \eta \equiv \frac{\Omega_K - \Omega}{\Omega_K} = -\frac{1}{2} \frac{\partial \ln P_g}{\partial \ln r} \left( \frac{H_g}{r} \right)^2 = C_{\eta} \left( \frac{H_g}{r} \right)^2, \qquad (6) $$

where $\Omega$ is the orbital frequency of the gas and we define

$$ C_{\eta} \equiv -\frac{1}{2} \frac{\partial \ln P_g}{\partial \ln r}, \qquad (7) $$

which depends on the gas structure. In this paper, $C_{\eta} = 11/8$ for $\Sigma_g \propto r^{-1}$ and $T \propto r^{-1/2}$.

## 2.2. Radial drifts of gas, pebbles, and silicate dust

The gas radial velocity including the effects of gas-solid friction is given as (Ida & Guillot 2016; Schoonenberg & Ormel 2017)

$$ v_{g,BK} = \frac{\Lambda}{1 + \Lambda^2 \tau_s^2} \left( 2Z\Lambda\tau_s\eta v_K + (1 + \Lambda\tau_s^2) v_{g,v} \right), \qquad (8) $$

where $Z$ is the midplane solid-to-gas ratio of pebbles or silicate dust, $\Lambda$ characterizes the strength of the back-reaction due to pile-up of either pebbles or silicate dust, $\tau_s$ is the Stokes number of pebbles or silicate dust, and $v_{g,v}$ is an unperturbed disk gas accretion velocity. $Z$ is defined as

$$ Z \equiv \frac{\rho_{p \ (\text{or } d})}{\rho_g}, \qquad (9) $$

where $\rho_g = \Sigma_g / \sqrt{2\pi H_g}$ and $\rho_{p \ (\text{or } d)} = \Sigma_{p \ (\text{or } d)} / \sqrt{2\pi H_{p \ (\text{or } d)}}$ are the midplane spatial density of the gas and that of pebbles (or silicate dust), respectively. $H_p$ and $H_d$ are the scale heights of pebbles and silicate dust, respectively. $\Lambda$ is defined as

$$ \Lambda \equiv \frac{\rho_g}{\rho_g + \rho_{p \ (\text{or } d)}} = \frac{1}{1 + Z}. \qquad (10) $$

Using $H_p$ or $H_d$, $\Lambda$ is rewritten as

$$ \Lambda^{-1} = 1 + \frac{\rho_{p \ (\text{or } d)}}{\rho_g} \qquad (11) $$

$$ \approx 1 + \left( \frac{H_g}{H_{p \ (\text{or } d)}} \right) \left( \frac{\Sigma_{p \ (\text{or } d)}}{\Sigma_g} \right) = 1 + h_{p/g \ (\text{or } d/g)}^{-1} Z_{\Sigma}, \qquad (12) $$

where $h_{p/g \ (\text{or } d/g)} \equiv H_{p \ (\text{or } d)}/H_g$. $Z_{\Sigma} \equiv \Sigma_{p \ (\text{or } d)}/\Sigma_g$ is the vertically averaged metallicities of pebbles or silicate dust. Using the dimensionless effective viscous parameter $\alpha_{acc}$ and the effective viscosity $v_{acc} = \alpha_{acc} c_s^2 \Omega_K^{-1}$, $v_{g,v}$ is given as

$$ v_{g,v} = -\frac{3v_{acc}}{r} \frac{\partial \ln(r^{1/2}v_{acc}\Sigma_g)}{\partial \ln r} = -\frac{3v_{acc}}{2r} (1 + 2Q), \qquad (13) $$

where $Q = \frac{\partial \ln(v_{acc}\Sigma_g)}{\partial \ln r}$ (Desch et al. 2017). For a steady-state disk, the gas accretion velocity is $v_{g,v}^* = -\frac{3v_{acc}}{2r}$ and is rewritten as

$$ v_{g,v}^* = -\frac{3v_{acc}}{2r} = -\frac{3\alpha_{acc} H_g^2 \Omega_K}{2r} = -\frac{3\alpha_{acc}}{2} \left(\frac{H_g}{r}\right)^2 v_K \qquad (14) $$

$$ \approx -\frac{3}{2} \alpha_{acc} \eta v_K \left( -\frac{1}{2} \frac{\partial \ln P_g}{\partial \ln r} \right)^{-1} = -\frac{3}{2} \alpha_{acc} \eta v_K C_{\eta}^{-1}, \qquad (15) $$

where $v_K$ is the Keplerian velocity.

<footnote n="2">² In Hyodo et al. (2019), $\alpha_{tur}$ is used to describe $\alpha_{Dr}$ and $\alpha_{Dz}$.</footnote>
In this work, we adopt a vertical two-layer model to describe the radial velocity of the gas: a pebble/dust-rich midplane layer with scale height $H_{p(or d)}$ and a pebble/dust-poor upper layer. The vertically averaged radial velocity of the gas $v_g$ is given as

$$v_{g,p(or d)} = \frac{v_{g,BK}H_{p(or d)} + (H_g - H_{p(or d)})v_{g,v}}{H_g} \quad (16)$$

We note that the above simple two-layer model well reproduces the results of Kanagawa et al. (2017) where the approximation of the Gaussian distribution is used.

The radial velocities of pebbles and silicate dust including the effects of gas-solid friction – drift back-reaction (Drift-BKR) – are given as (Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019)

$$v_p = - \frac{\Lambda}{1 + \Lambda^2 \tau_{s,p}^2} \left( 2\tau_{s,p} \Lambda \eta v_K - v_{g,v} \right) \quad (17)$$

$$v_d = - \frac{\Lambda}{1 + \Lambda^2 \tau_{s,d}^2} \left( 2\tau_{s,d} \Lambda \eta v_K - v_{g,v} \right), \quad (18)$$

where $\tau_{s,p}$ and $\tau_{s,d}$ are the Stokes number of pebbles and silicate dust, respectively.

## 2.3. Diffusion of gas and solids

The radial and vertical diffusivities of gas and vapor ($D_{Dr,g}$ and $D_{Dz,g}$) are given as

$$D_{Dr,g} = v_{Dr} = \alpha_{Dr} c_s^2 \Omega_K^{-1} \quad (19)$$

$$D_{Dz,g} = v_{Dz} = \alpha_{Dz} c_s^2 \Omega_K^{-1}, \quad (20)$$

where $v_{Dr}$ and $v_{Dz}$ are the turbulent viscosities that regulate radial and vertical diffusions in association with dimensionless turbulence parameters $\alpha_{Dr}$ and $\alpha_{Dz}$, respectively.

The nature of solid diffusion in the gas is partial coupling with the gas eddies (Youdin & Lithwick 2007). Hyodo et al. (2019) considered that the diffusivity of pebbles (and that of silicate dust) is reduced owing to the gas-solid friction as pile-up of solids proceeds – diffusion back-reaction (Diff-BKR) – and radial and vertical diffusivities of pebbles and silicate dust are given as

$$D_{Dr,p(or d)} = D_{Dr,g} \times \frac{\Lambda^K}{1 + \tau_{s,p(or d)}^2} \quad (21)$$

$$D_{Dz,p(or d)} = D_{Dz,g} \times \frac{\Lambda^K}{1 + \tau_{s,p(or d)}^2}, \quad (22)$$

where K is the coefficient. K = 1 could be applied because the diffusivity is $D_g \propto c_s^2 \propto 1/\rho_g$ under a constant P assumption. Using the effective gas density, $\rho_{g,eff} = \rho_g + \rho_{p(or d)}$, the effective sound velocity of a mixture of gas and small particles is given as

$$c_{s,eff}^2 = c_{s,K=0}^2 \left( \frac{\rho_g}{\rho_g + \rho_{p(or d)}} \right) = c_{s,K=0}^2 \Lambda, \quad (23)$$

where $c_{s,K=0}$ is the sound velocity when K = 0 (Eq. (2)). However, from an energy dissipation argument, the lower limit is K = 1/3. We note that the result of solid pile-up around the snow line does not significantly depend on the choice of K as long as the back reaction to diffusion is included (i.e., K > 0). In this work, we use K = 1 as a representative case.

## 2.4. Equations of radial transport of gas, vapor, and solids

In this paper, we solve the radial motions of gas, pebbles, silicate dust, and water vapor. Pebbles are modeled as a homogenous mixture of water ice and micron-sized silicate dust (see also Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019). Due to the sublimation of icy pebbles near the snow line, water vapor and silicate dust are produced.

The governing equations of the surface density of the gas $\Sigma_g$, that of pebbles $\Sigma_p$, that of silicate dust $\Sigma_d$, that of water vapor $\Sigma_{vap}$, and the number density of pebbles $N_p$ are given as (Desch et al. 2017; Hyodo et al. 2019);

$$\frac{\partial \Sigma_g}{\partial t} = -\frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_g v_g \right) \quad (24)$$

$$\frac{\partial \Sigma_p}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_p v_p - r D_{Dr,p} \Sigma_g \frac{\partial}{\partial r} \left( \frac{\Sigma_p}{\Sigma_g} \right) \right) = \dot{\Sigma}_p \quad (25)$$

$$\frac{\partial \Sigma_d}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_d v_d - r D_{Dr,d} \Sigma_g \frac{\partial}{\partial r} \left( \frac{\Sigma_d}{\Sigma_g} \right) \right) = \dot{\Sigma}_d \quad (26)$$

$$\frac{\partial \Sigma_{vap}}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_{vap} v_{g,v} - r D_{Dr,g} \Sigma_g \frac{\partial}{\partial r} \left( \frac{\Sigma_{vap}}{\Sigma_g} \right) \right) = \dot{\Sigma}_{vap} \quad (27)$$

$$\frac{\partial N_p}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left( r N_p v_p - r D_{Dr,p} N_g \frac{\partial}{\partial r} \left( \frac{N_p}{N_g} \right) \right) = 0. \quad (28)$$

The right-hand sides of Eqs. (25)–(27) are due to sublimation of icy pebbles, the recondensation of water vapor, and sticking of silicate dust onto pebbles (see Hyodo et al. 2019). The evolution of $N_p$ needs to be solved because the size of pebbles changes during sublimation and condensation, and the mass of pebbles is given as $m_p = \Sigma_p/N_p$ under the single-size approximation.

The decrease rate of the pebble mass $m_p$ is given by (Lichtenegger & Komle 1991; Ros et al. 2019)

$$\frac{dm_p}{dt} = 4\pi r_p^2 v_{th}(\rho_{vap} - \rho_{sat}), \quad (29)$$

where $r_p$ is the particle physical radius, $\rho_{vap}$ and $\rho_{sat}$ are vapor and saturation spatial densities, and $v_{th}$ is the averaged normal component of the velocity passing through the particle surface. The averaged normal velocity of only outgoing particles with Maxwell distribution is

$$v_{th} = \frac{c_s}{\sqrt{2\pi}} = \frac{1}{\sqrt{2\pi}} \left( \frac{k_B T}{\mu m_{proton}} \right)^{1/2}, \quad (30)$$

where $\mu$ is the mean molecular weight. Another definition of the thermal velocity is the averaged value of the magnitude of three-dimensional velocity, which is given by

$$v_{th,3D} = \sqrt{\frac{8}{\pi}} c_s. \quad (31)$$

Because this definition of $v_{th}$ does not include the effect of oblique ejection/collisions, the cross section must be considered rather than the surface density, that is,

$$\frac{dm_p}{dt} = \pi r_p^2 v_{th,3D} (\rho_{vap} - \rho_{sat}). \quad (32)$$

Both Eqs. (29) and (55) give the same $dm_p/dt$. We note that Schoonenberg & Ormel (2017) and Hyodo et al. (2019) used $v_{th,3D}$ for Eq. (29), which caused overestimation of the sublimation and recombination rates by a factor of 4.
## 2.5. Scale heights of pebbles and silicate dust

The scale heights of pebbles and silicate dust describe the degree of concentration of solids in the disk midplane ($\rho_{p(or d)} \propto H_{p(or d)}^{-1}$).

In the steady-state, the scale height of pebbles is regulated by the vertical turbulent stirring $H_{p,tur}$ (Dubrulle et al. 1995; Youdin & Lithwick 2007; Okuzumi et al. 2012) as<fnref n="3" />

$$H_{p,tur} = \left( 1 + \frac{\tau_{s,p}}{\alpha_{Dz} (1 + Z_p)^{-K}} \right)^{-1/2} H_g. \quad (33)$$

However, for small $\alpha_{Dz}$ (i.e., for a small $H_{p,tur}$), a vertical shear Kelvin-Helmholtz (KH) instability would prevent the scale height of pebbles from being a smaller value. The scale height of pebbles regulated by a KH instability $H_{p,KH}$ is given (see the derivation in Sect. 4.1) as

$$H_{p,KH} \approx R i^{1/2} \frac{Z^{1/2}}{(1+Z)^{3/2}} C_{\eta} \left(\frac{H_g}{r}\right) H_g, \quad (34)$$

where $Ri = 0.5$ is used here as a critical value for a KH instability to operate. Thus, the scale height of pebbles $H_p$ is given as

$$H_p = \max \{H_{p,tur}, H_{p,KH}\}. \quad (35)$$

Using a Monte Carlo simulations, Paper I derived the scale height of silicate dust $H_d$ as a function of the scaled distance to the snow line $\Delta\tilde{x}_{snow} = (r - r_{snow})/H_g$ and the scaled sublimation width $\Delta\tilde{x}_{subl} = \Delta x_{subl}/H_g$ (Sect. 3) as

$$H_d \approx (h_{d/g,0}^{-1} + h_{d/g,*}^{-1})^{-1} H_g, \quad (36)$$

where

$$h_{d/g,0} = \left(1 + \frac{\tau_{s,d}}{\alpha_{Dz}}\right)^{-1/2}, \quad (37)$$

and

$$h_{d/g,*} \approx \left( h_{p/g,0}^2 + \frac{2}{3} \frac{\alpha_{Dz}}{\alpha_{acc}} \frac{\Delta \tilde{x}}{H_g/r} \right)^{1/2} \times \left( 1 + \frac{2}{3} \frac{\alpha_{Dr}/\alpha_{acc}}{1 + (C_{\text{diff}}\alpha_{Dr}/\alpha_{acc})^2} \frac{1}{(H_g/r)(\Delta \tilde{x} + \epsilon)} \right), \quad (38)$$

where $C_{\text{diff}} = 10$, $\Delta\tilde{x} = \max\{\Delta\tilde{x}_{subl}, \Delta\tilde{x}_{snow}\}$, and $\epsilon = 0.01$ is a softening parameter that prevents unphysical peak for very small $\Delta\tilde{x}_{subl}$ at the snow line ($\Delta\tilde{x}_{snow} = 0$), respectively (Paper I). The scale height of pebbles is given at the snow line (with $\tau_{s,p} = 0.1$) as

$$h_{p/g,0} = \max \left\{ \frac{H_{p,tur}|_{K=0}}{H_g}, \frac{H_{p,KH}|_{Z=0.5}}{H_g} \right\}. \quad (39)$$

Previous works adopted different simplifying prescriptions for $H_d$ (e.g., Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019), leading to very different results for the resultant pile-up of silicate dust and pebbles (Appendix A). Equations (36)–(38) provide a more realistic treatment based on the 2D calculation from Paper I.

## 2.6. Numerical settings

Following Paper I, the ratio of the mass flux of pebbles $\dot{M}_p$ to that of the gas $\dot{M}_g$ is a parameter and we define it as

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

We set $\dot{M}_g = 10^{-8} M_{\odot} \text{yr}^{-1}$. The inner and outer boundaries of our 1D simulations are set to $r_{\text{in}} = 0.1$ au and $r_{\text{out}} = 5.0$ au, respectively. We set the initial constant size of pebbles (radius of $r_p$) with $\tau_{s,p}(r) = 0.1$ at $r=4$ au before sublimation takes place (Okuzumi et al. 2012; Ida & Guillot 2016, Appendix B). The mass fraction of silicate in icy pebbles is set to $f_{d/p} = 0.5$ and we modeled that silicate dust uniformly embedded in icy pebbles is micron-sized small grains (Schoonenberg & Ormel 2017; Hyodo et al. 2019; Ida et al. 2021). The surface density of pebbles at the outer boundary is set by considering Drift-BKR and Diff-BKR (see Sect. 4.2 and Eq. (86)), which is fixed during the simulations. We consider sublimation/condensation of ice as well as release/recycling of silicate dust around the snow line (Hyodo et al. 2019). 1D simulations are initially run without back-reactions to reach a steady-state (~$10^5$ yr) and then we turn on the back-reactions (Drift-BKR and Diff-BKR with $K=1$) to see the further evolution (up to $5 \times 10^5$ yr).

# 3. Sublimation of drifting icy pebbles

A description of $H_d$ is critical to understand the resultant pile-up of silicate dust inside the snow line. Also, $H_d$ describes the efficiency of recycling of silicate dust onto icy pebbles and pile-up of pebbles just outside the snow line because its efficiency is $\propto H_d^{-1}$ (Hyodo et al. 2019) (see Appendix A for a detailed comparison of different models of $H_d$). Paper I derived a semi-analytical expression of $H_d$ that depends on the sublimation width $\Delta x_{subl}$ (see Eq. (36)). Although the global distribution of water within protoplanetary disks has been a topic of research (e.g., Ciesla & Cuzzi 2006), the detailed study of the radial width where drifting pebbles sublimate around the snow line (i.e., sublimation width) has not been conducted yet.

Below, we discuss the radial width of sublimation of drifting icy pebbles $\Delta x_{subl}$ in a protoplanetary disk considering the local temperature-pressure environment. Here, we consider the cases of pure icy pebbles (i.e., $f_{d/p} = 0$ and $\tau_{s,p} = 0.1$ at 4 au). We neglect recycling processes of the recondensation of water vapor as well as sticking of silicate dust onto pebbles outside the snow line via their diffusive processes, so that we can purely focus on the sublimation process.

In the 1D simulations, diffusion radially mixes different-sized sublimating pebbles near the snow line, which comes from solving $N_d$ evolution (i.e., the effective size of pebbles at a given radial distance is calculated). In this section, to suppress these mixing effects and to see the size evolution of a single drifting pebble, we set $\alpha_{Dr}/\alpha_{acc} = 10^{-3}$ for pebbles and silicate dust, while we change $\alpha_{Dr}/\alpha_{acc}$ for the evolution of vapor. The arguments in this section are applicable to the 1:1 rock-ice-mixed pebbles (see Sect. 2.6) because the size and density of pebbles are expected to only weakly change with slightly different compositions.

The sublimation of drifting icy pebbles takes place as long as $P_{\text{vap}}(r) < P_{\text{sat}}(r)$ where $P_{\text{vap}}(r)$ is the water vapor pressure and $P_{\text{sat}}$ is the saturation vapor pressure (Fig. 2). $P_{\text{vap}} > P_{\text{sat}}$ indicates

<footnote n="3">³ Previous studies (Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019) did not consider the effects of Diff-BKR on the scale height of pebbles (i.e., $K=0$ is assumed for Eq. (33)).</footnote>
Advection-dominated regime ($\alpha_{Dr}/\alpha_{acc} \ll 1$)

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

Diffusion-dominated regime ($\alpha_{Dr}/\alpha_{acc} \approx 1$)

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

**Fig. 2.** Schematic illustration of the sublimation of drifting icy pebbles. In the case of the advection-dominated regime (**left**; Sect. 3.1), gas accretion onto the central star dominates the radial transport of vapor. The loss of vapor due to the inward advection of the vapor then balances the vapor supply from sublimating icy pebbles. Thus, drifting pebbles sublimate outside the snow line (i.e., $r > r_{snow}$). In the case of the diffusion-dominated regime (**right**; Sect. 3.2), water vapor produced inside the snow line efficiently diffuses outside the snow line, leading to a super-saturation state ($P_{vap}(r) > P_{sat}(r)$). The sublimation of pebbles can therefore take place only inside the snow line (i.e., $r < r_{snow}$).

a super-saturation state. The saturation vapor pressure of water is given as

$$P_{\text{sat}}(r) = P_0 \exp\left(-\frac{T_0}{T(r)}\right), \qquad (41)$$

where $P_0 = 1.14 \times 10^{13}$ g cm$^{-1}$ s$^{-2}$ and $T_0 = 6062$ K (Lichtenegger & Komle 1991), which can be rewritten as

$$\log_{10}\left(\frac{P_{\text{sat}}(r)}{\text{dyn cm}^{-2}}\right) \approx 13.06 - \frac{2632 \text{ K}}{T(r)}. \qquad (42)$$

As shown in Fig. 2, the sublimation width of drifting pebbles can be divided into two regimes depending on $\alpha_{Dr}/\alpha_{acc}$: the advection-dominated regime (Sect. 3.1; $\Delta x_{subl} \approx 2H_g$ for $\alpha_{Dr} \ll \alpha_{acc}$) and the diffusion-dominated regime (Sect. 3.2; $\Delta x_{subl} \approx 0.2H_g$ for $\alpha_{Dr} \sim \alpha_{acc}$). Figure 3 shows the results of 1D simulations as a function of $\alpha_{Dr}/\alpha_{acc}$ (the top panel being a radial distance of the snow line and the bottom panel being the sublimation width). Here, the snow line $r_{snow}$ is defined by the innermost radial distance where $P_{vap}(r) = P_{sat}(r)$. The sublimation width $\Delta x_{subl}$ is defined by the radial width where drifting pebbles sublimate by 16.7 wt.% to 88.3 wt.% from the initial value. The advection-dominated and diffusion-dominated regimes are clearly identified by their different $\Delta x_{subl}$ values in Fig. 3. We explain the details of these different regimes below.

As reference values, $\tau_{s,p} \sim 0.06$, $r_p \sim 2$ cm at $r_{snow} \sim 2.4$ au for $F_{p/g} = 0.1$ and $\alpha_{acc} = 10^{-2}$, while $\tau_{s,p} \sim 0.06$, $r_p \sim 20$ cm at $r_{snow} \sim 2.1$ au for $F_{p/g} = 0.1$ and $\alpha_{acc} = 10^{-3}$. Here, the drag relations correspond to those of the Epstein regime. Below, a notation of “snow” indicates values at the snow line.

### 3.1. Advection-dominated regime

In the advection-dominated regime ($\alpha_{Dr} \ll \alpha_{acc}$), the gas accretion onto the central star dominates the radial transport of vapor. Outside the snow line, the loss of vapor due to the inward advection of the gas balances the supply of vapor from inward drifting

icy pebbles. Thus, the icy pebbles locally sublimate to satisfy $P_{sat}(r) \approx P_{vap}(r)$ outside the snow line ($r > r_{snow}$; see Fig. 2). Neglecting the effects of the radial diffusion implies that

$$P_{\text{sat}}(r) = P_{\text{vap}}(r), \qquad (43)$$

which is rewritten assuming a rapid vertical mixing with the background gas (the molecular weight $\mu_{vap}$) as

$$P_0 \exp\left(-\frac{T_0}{T(r)}\right) = \frac{k_B T(r)}{\mu_{\text{vap}} m_{\text{proton}}} \frac{\Sigma_{\text{vap}}(r)}{\sqrt{2\pi H_g(r)}}, \qquad (44)$$

where $T(r) \propto r^{-\beta}$, $H_g \propto r^{(3-\beta)/2}$, and $\nu_{acc}(r) \propto \alpha_{acc} r^{(3-2\beta)/2}$.

For $r \ge r_{snow}$, the surface density of the water vapor that satisfies the local saturation state is given as

$$\Sigma_{\text{vap}}(r) = \frac{\sqrt{2\pi H_g(r)\mu_{\text{vap}}m_{\text{proton}}}}{k_B T(r)} P_0 \exp\left(-\frac{T_0}{T(r)}\right). \qquad (45)$$

Using the surface density of the water vapor and temperature at the snow line (i.e., $\Sigma_{vap,snow} \equiv \Sigma_{vap}(r_{snow})$ and $T_{snow} \equiv T(r_{snow})$, respectively), and Eq. (45), $P_0$ is removed and $\Sigma_{vap}(r)$ for $r \ge r_{snow}$ is given as

$$\begin{aligned} \Sigma_{\text{vap}}(r \ge r_{\text{snow}}) &= \Sigma_{\text{vap,snow}} \frac{T_{\text{snow}}}{T(r)} \frac{H_g(r)}{H_g(r_{\text{snow}})} \exp\left(-\frac{T_0}{T(r)} + \frac{T_0}{T_{\text{snow}}}\right) \\ &= \Sigma_{\text{vap,snow}} \left(\frac{r}{r_{\text{snow}}}\right)^{\beta} \left(\frac{r}{r_{\text{snow}}}\right)^{\frac{3-\beta}{2}} \exp\left[-T_0 \left(\frac{1}{T(r)} - \frac{1}{T_{\text{snow}}}\right)\right] \\ &= \Sigma_{\text{vap,snow}} \left(\frac{r}{r_{\text{snow}}}\right)^{\frac{3+\beta}{2}} \exp\left[\frac{T_0}{T_{\text{snow}}} \left(1 - \left(\frac{r}{r_{\text{snow}}}\right)^{\beta}\right)\right]. \end{aligned} \qquad (46)$$

Inside the snow line (i.e., $r \le r_{snow}$), $\Sigma_{vap}(r)$ is written in terms of the gas mass flux $\dot{M}_g$ and the effective viscosity $\nu_{acc}$ as

$$\Sigma_{\text{vap}}(r \le r_{\text{snow}}) = \frac{(1 - f_{d/p}) F_{p/g} \dot{M}_g}{3\pi \nu_{\text{acc}}(r)} = \frac{(1 - f_{d/p}) F_{p/g} \dot{M}_g}{3\pi \alpha_{\text{acc}} H_g^2(r) \Omega_K(r)}. \qquad (47)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Radial location of the snow line (top panel) and sublimation width of drifting pebbles $\Delta x_{subl}$ (bottom panel) as a function of $\alpha_{Dr}/\alpha_{acc}$. Here, $f_{d/p} = 0$ (i.e., pure icy pebbles). Circles ($\alpha_{acc} = 10^{-2}$) and squares ($\alpha_{acc} = 10^{-3}$) represent the results of 1D simulations. Blue, black and red colors indicate those for $F_{p/g} = 0.03$, 0.1 and 0.3, respectively. Analytically derived locations of the snow lines (solving Eq. (44)) are shown by the gray lines (top panel). The light-blue line in the bottom panel shows the analytical $\Delta x_{subl}$ in the case of the advection-dominated regime (Eq. (54)) with $r_{snow} = 2.4$ au and $\beta = 0.5$). The light-green line in the bottom panel shows the analytical $\Delta x_{subl}$ in the case of the diffusion-dominated regime (Eq. (68)) with $\tau_{s,p} = 0.06$, $r_p = 2.3$ cm, $r_{snow} = 2.4$ au). The gray dashed line in the bottom panel shows the analytically derived boundary between the advection- and diffusion-dominated regimes (Eq. (73)) with $r = 2.4$ au and $\beta = 0.5$). A fitting function is shown by a black line (Eq. (74)).</figcaption></figure>

The above equation shows a good match with our 1D simulations for $\alpha_{Dr} \ll \alpha_{acc}$, and thus the advection dominates the radial mass transfer of vapor for such a low-diffusivity case. As discussed above, the local vapor pressure equals the saturation pressure as

$$P_{\text{vap}}(r) = P_0 \exp\left(-\frac{T_0}{T(r)}\right). \qquad (48)$$

Using the vapor pressure at the snow line $P_{\text{vap,snow}} \equiv P(r_{\text{snow}})$ and removing $P_0$, the equation is rewritten as

$$\begin{aligned} \frac{P_{\text{vap}}(r)}{P_{\text{vap,snow}}} &= \exp\left[-\frac{T_0}{T(r)} + \frac{T_0}{T_{\text{snow}}}\right] \\ &= \exp\left[\frac{T_0}{T_{\text{snow}}} \left(1 - \left(\frac{r}{r_{\text{snow}}}\right)^{\beta}\right)\right]. \end{aligned} \qquad (49)$$

Using $\Delta x^* = (r - r_{\text{snow}})/r_{\text{snow}} \ll 1$ near the snow line (and $\Delta x^* > 0$), the above equation is approximated as follows:

$$\begin{aligned} \frac{P_{\text{vap}}(r)}{P_{\text{vap,snow}}} &= \exp\left[\frac{T_0}{T_{\text{snow}}} \left(1 - \left(\frac{r - r_{\text{snow}}}{r_{\text{snow}}} + 1\right)^{\beta}\right)\right] \\ &= \exp\left[\frac{T_0}{T_{\text{snow}}} \left(1 - (1 + \Delta x^*)^{\beta}\right)\right] \end{aligned} \qquad (50)$$

$$\approx \exp\left[\frac{T_0}{T_{\text{snow}}} \left(1 - (1 + \beta \Delta x^*)\right)\right] \qquad (51)$$

$$= \exp\left(-\frac{\beta T_0}{T_{\text{snow}}} \Delta x^*\right), \qquad (52)$$

where $P_{\text{vap}}(r)/P_{\text{vap,snow}} = 0$ is the innermost radius where the sublimation takes place and $P_{\text{vap}}(r)/P_{\text{vap,snow}} = 1$ is the radial location where all ice sublimates (in reality, a very small fraction of pebbles could sublimate inside the snow line). Here, the sublimation width $\Delta x_{subl}$ is equivalent to the radial width of $P_{\text{vap}}(r)/P_{\text{vap,snow}} = 0.167$ to 0.833, which corresponds to the change in mass of pebbles ($\dot{m}_p \propto \dot{\Sigma}_{\text{vap}} \propto \dot{P}_{\text{vap}}$). Using Eq. (52), $\Delta x_{subl}^* = \Delta x_{subl}/r_{\text{snow}}$ is given by

$$\begin{aligned} \Delta x_{subl}^* &= \Delta x_{(P_{\text{vap}}/P_{\text{vap,snow}} = 0.167)}^* - \Delta x_{(P_{\text{vap}}/P_{\text{vap,snow}} = 0.883)}^* \\ &= -\frac{T_{\text{snow}}}{\beta T_0} \left(\ln 0.167 - \ln 0.883\right) \end{aligned} \qquad (53)$$

$$\approx 0.093 \times \left(\frac{T_{\text{snow}}}{170 \text{ K}}\right) \left(\frac{T_0}{6062 \text{ K}}\right)^{-1} \left(\frac{\beta}{0.5}\right)^{-1}, \qquad (54)$$

where $T_{\text{snow}}$ depends on $\alpha_{acc}$ and $F_{p/g}$ (Fig. 3). For $r_{\text{snow}} \approx 2.4$ au with $\beta = 0.5$ ($T_{\text{snow}} \approx 170$ K), $\Delta x_{subl} = \Delta x_{subl}^* \times r_{\text{snow}} \approx 0.22$ au ($\approx 2.2H_g$), which is consistent with the results of 1D simulations (the light-blue line in Fig. 3 bottom panel). Both 1D simulations and analytical arguments show that the sublimation width in the advection-dominated regime is very weakly dependent on $F_{p/g}$ and $\alpha_{Dr}/\alpha_{acc}$. Equation (54) also indicates that the sublimation width in the advection-dominated regime is independent of the Stokes number of pebbles and initial pebble physical size, which are supported by 1D simulations (Fig. 3; $r_p \sim 2$ cm and $\sim 20$ cm for $\alpha_{acc} = 10^{-2}$ and $10^{-3}$, respectively). These are because the sublimation width in the advection-dominated regime is regulated by $P_{\text{sat}}$ and $T_{\text{snow}}$.

## 3.2. Diffusion-dominated regime

In the diffusion-dominated regime ($\alpha_{Dr} \sim \alpha_{acc}$), the local sublimation of drifting pebbles does not take place outside of the snow line ($r > r_{\text{snow}}$; see Fig. 2) because the water vapor produced in the inner region efficiently diffuses outward and the local pressure by the diffused vapor becomes larger than the local saturation vapor pressure, being a super-saturation state (i.e., $P_{\text{vap}}(r) > P_{\text{sat}}(r)$ for $r > r_{\text{snow}}$). Thus, the sublimation only takes place after pebbles pass through the snow line where $P_{\text{vap}}(r) < P_{\text{sat}}$ (i.e., $r < r_{\text{snow}}$).

The solar composition H/He gas sound velocity is given by $c_s = \sqrt{k_B T/2.34 m_{\text{proton}}^*} \approx 10^5 (T/280 \text{ K})^{1/2}$ cm s$^{-1}$. The water vapor sound velocity with the molecular weight $\mu = 18$ is given by $c_{s,v} = \sqrt{k_B T/\mu m_{\text{proton}}^*} \approx 0.36 c_s \approx 0.28 \times 10^4 (T/170 \text{ K})^{1/2}$ cm s$^{-1}$. Because the saturation pressure and vapor pressure are respectively given by $P_{\text{sat}} = c_{s,v}^2 \rho_{\text{sat}}$ and $P_{\text{vap}} = c_{s,v}^2 \rho_{\text{vap}}$ and $m_p = (4\pi/3)\rho_{\text{bulk}} r_p^3$, where $\rho_{\text{bulk}} \sim 1$ g cm$^{-3}$)
is the bulk density of the particle, the rate of change in particle size is given as (see Eqs. (29) and (55))

$$ \frac{dr_p}{dt} = - \frac{1}{\sqrt{2\pi} c_{s,v}} \frac{P_{\text{sat}} - P_{\text{vap}}}{\rho_{\text{bulk}}} \quad (55) $$

For the disk density similar to MMSN, the snow line corresponds to $T = T_{\text{snow}} \approx 170$ K. Using Eq. (44), the vapor pressure with $170$ K at $r = r_{\text{snow}}$ is $P_{\text{vap,snow}} [\text{dyn cm}^{-2}] = 10^{13.06-2632/170} \approx 10^{-2.4} \approx 3.8 \times 10^{-3}$. We assume that $\Sigma \propto r^{-\beta_S}$, $T \propto r^{-\beta_T}$, $P_{\text{vap}} \propto r^{-\beta_P}$.

We consider the case where the vapor pressure $P_{\text{vap}}$ is proportional to disk gas pressure $P_{\text{H/He}}$. In that case, $\beta_P = \beta_S + \beta_T/2 + 3/2$ for ideal gas. We note that $P_{\text{vap}} \propto P_{\text{H/He}}$ is established in the limit of efficient radial diffusion (i.e., diffusion-dominated regime). For a nominal case, $\beta_S = 1$, $\beta_T = 1/2$, and $\beta_P = 11/4$. Inside the snow line ($r < r_{\text{snow}}$), we can write

$$ P_{\text{sat}} \approx 10^{13.06 - \frac{2632 K}{r_{\text{snow}}(r/r_{\text{snow}})^{-\beta_T}}} = P_{\text{vap,snow}} \times 10^{15.5[1-(r/r_{\text{snow}})^{\beta_T}]}, \quad (56) $$

$$ P_{\text{vap}} \approx P_{\text{vap,snow}} (r/r_{\text{snow}})^{-\beta_P}, \quad (57) $$

where $T_{\text{snow}} \approx 170$ K was used. At $r \sim r_{\text{snow}}$, Eq. (55) is reduced to

$$ \begin{aligned} \frac{dr_p}{dt} &\approx \frac{P_{\text{vap}}}{\sqrt{2\pi} c_{s,v} \rho_{\text{bulk}}} \left(1 - \frac{P_{\text{sat}}}{P_{\text{vap}}}\right) \\ &\approx C_P \left(\frac{r}{r_{\text{snow}}}\right)^{-\beta_P + \beta_T/2} \left[1 - 10^{15.5(1-(r/r_{\text{snow}})^{\beta_T})}\left(\frac{r}{r_{\text{snow}}}\right)^{\beta_P}\right] c_s, \end{aligned} \quad (58) $$

where

$$ \begin{aligned} C_P &\equiv \frac{P_{\text{vap,snow}}}{\sqrt{2\pi} c_{s,v} c_s \rho_{\text{bulk}}} = \frac{P_{\text{vap,snow}}}{\sqrt{2\pi} \times 0.36 c_s^2 \rho_{\text{bulk}}} \\ &= 0.70 \times 10^{-12} \left( \frac{c_{s,snow}}{0.78 \times 10^5 \text{ cm s}^{-1}} \right)^{-2} \left( \frac{\rho_{\text{bulk}}}{1 \text{ g cm}^{-3}} \right)^{-1}. \end{aligned} \quad (59) $$

Defining $\Delta\tilde{r} \equiv (r - r_{\text{snow}})/r_{\text{snow}}$,

$$ \begin{aligned} \frac{dr_p}{dt} &\approx C_P \left(\frac{r}{r_{\text{snow}}}\right)^{-\beta_P + \beta_T/2} \left[1 - 10^{-15.5\beta_T \Delta\tilde{r}} (1 + \beta_P \Delta\tilde{r})\right] c_s \\ &\approx C_P (35.7\beta_T - \beta_P) \Delta\tilde{r} c_s. \end{aligned} \quad (60) $$

For the nominal parameters, $35.7\beta_T - \beta_P \approx 15.1$. The radial drift of a sublimating pebble is given by

$$ \frac{dr}{dt} = - \Lambda^2 \frac{2\tau_{s,p}}{1 + \Lambda^2 \tau_{s,p}^2} \eta v_K + \Lambda \frac{1}{1 + \Lambda^2 \tau_{s,p}^2} v_{g,v}^* \quad (61) $$

In the case of $\Lambda \approx 1$ and $\tau_{s,p}^2 \ll 1$,

$$ \frac{d\Delta\tilde{r}}{dt} \approx -(\beta_P \tau_{s,p} + (3/2)\alpha_{\text{acc}}) \left(\frac{H_g}{r}\right) c_s r_{\text{snow}}^{-1} \quad (62) $$

From Eqs. (60) and (62) with $\tilde{r}_p = r_p/r_{p,snow}$,

$$ \frac{d\tilde{r}_p}{d\Delta\tilde{r}} \approx -C_P \frac{35.7\beta_T - \beta_P}{\beta_P \tau_{s,p} + (3/2)\alpha_{\text{acc}}} \frac{r_{\text{snow}}}{r_{p,snow}} \left(\frac{H_g}{r}\right)_{\text{snow}}^{-1} \Delta\tilde{r} \quad (63) $$

$$ \begin{aligned} &= -0.98 \times 10^4 \left( \frac{H_g}{0.04} \right)_{\text{snow}}^{-1} \\ &\quad \times \left( \frac{\rho_{\text{bulk}}}{1 \text{ g cm}^{-3}} \right)^{-1} \left( \frac{r_{p,snow}}{1 \text{ cm}} \right)^{-1} \left( \frac{r_{\text{snow}}}{2.5 \text{ au}} \right) \frac{1}{\beta_P \tau_{s,p} + (3/2)\alpha_{\text{acc}}} \Delta\tilde{r}. \end{aligned} \quad (64) $$

We further set $\tau_{s,p} = \tau_{s,p,snow} \tilde{r}_p^\zeta$ (where $\zeta = 2$ for Stokes drag and $\zeta = 1$ for Epstein drag). Then,

$$ \frac{d\tilde{r}_p}{d\Delta\tilde{r}} \approx \begin{cases} -3.6 \times 10^4 C_R \left(\frac{\tau_{s,p,snow}}{0.1}\right)^{-1} \left(\frac{r_{p,snow}}{1 \text{ cm}}\right)^{-1} \tilde{r}_p^{-\zeta} \Delta\tilde{r} \\ \left[ \text{for } \tilde{r}_p > \tilde{r}_{p,\text{tr}} \equiv \left( \frac{\alpha_{\text{acc}}}{1.83\tau_{s,p,snow}} \right)^{1/\zeta} \right] \\ -6.6 \times 10^5 C_R \left(\frac{\alpha_{\text{acc}}}{10^{-2}}\right)^{-1} \left(\frac{r_{p,snow}}{1 \text{ cm}}\right)^{-1} \Delta\tilde{r} \\ \text{[otherwise]} \end{cases} \quad (65) $$

$$ C_R \equiv \left( \frac{H_g/r}{0.04} \right)_{\text{snow}}^{-1} \left( \frac{\rho_{\text{bulk}}}{1 \text{ g cm}^{-3}} \right)^{-1} \left( \frac{r_{\text{snow}}}{2.5 \text{ au}} \right). \quad (66) $$

We integrate the upper equation from $\tilde{r}_p = 1$ as

$$ \begin{aligned} \tilde{r}_p(\Delta\tilde{r}) &\approx \left[1 - 3.6 \times 10^4 C_R \left(\frac{\zeta + 1}{2}\right) \left(\frac{\tau_{s,p,snow}}{0.1}\right)^{-1} \left(\frac{r_{p,snow}}{1 \text{ cm}}\right)^{-1} \Delta\tilde{r}^2\right]^{1/(\zeta+1)}. \end{aligned} \quad (67) $$

If we define the sublimation width $\Delta x_{\text{subl}}$ by the radial separation between $\tilde{r}_p^3 = \xi_1$ and $\tilde{r}_p^3 = \xi_2$ (as $m_p \propto \tilde{r}_p^3$),

$$ \begin{aligned} \Delta x_{\text{subl}} &\approx \left[ \left(1 - \xi_1^{(\zeta+1)/3}\right)^{1/2} - \left(1 - \xi_2^{(\zeta+1)/3}\right)^{1/2} \right] \\ &\quad \times 0.13 C_R^{-1/2} \left(\frac{\zeta + 1}{2}\right)^{-1/2} \left(\frac{H_g}{0.04}\right)^{-1} \left(\frac{\tau_{s,p,snow}}{0.1}\right)^{1/2} \left(\frac{r_{p,snow}}{1 \text{ cm}}\right)^{1/2} H_g. \end{aligned} \quad (68) $$

Equation (68) (with $\zeta = 1$, $\tau_{s,p} = 0.06$, $r_p = 2.3$ cm, $r_{\text{snow}} = 2.4$ au, $\xi_1 = 0.167$, and $\xi_2 = 0.883$) is shown by the light-green line in the bottom panel of Fig. 3. The analytical arguments above are generally in accordance with the 1D numerical results. The sublimation width in the diffusion-dominated regime depends on $F_{p/g}$, in that, a larger $F_{p/g}$ leads to a smaller $\Delta x_{\text{subl}}$. This is because the location of the snow line becomes closer to the star for a larger $F_{p/g}$, which leads to a larger $P_{\text{sat}} - P_{\text{vap}}$ (see Eq. (55)). In the case of $\alpha_{\text{acc}} = 10^{-3}$ (squares in Fig. 3), the snow line is located closer to the star (a larger $P_{\text{sat}} - P_{\text{vap}}$) than that of $\alpha_{\text{acc}} = 10^{-2}$, which makes $\Delta x_{\text{subl}}$ smaller, whereas the size of pebbles is 10 times larger than that of $\alpha_{\text{acc}} = 10^{-2}$, which makes $\Delta x_{\text{subl}}$ larger.

### 3.3. Boundary between advection/diffusion-dominated regimes

The boundary between two regimes – the advection-dominated regime and the diffusion-dominated regime – can be evaluated by considering the mass fluxes of advection and diffusion. In the steady-state for an unperturbed case (no back-reaction), the
rate of change in the surface density of the vapor (Eq. (27)) is
described as:

$$
\begin{align}
& \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_{\text{vap}} v_{\text{g},v} - r D_{\text{Dr,g}} \Sigma_g \frac{\partial}{\partial r} \left( \frac{\Sigma_{\text{vap}}}{\Sigma_g} \right) \right) \nonumber \\
&= \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_{\text{vap}} v_{\text{g},v} - D_{\text{Dr,g}} \Sigma_{\text{vap}} \frac{\partial \ln \Sigma_{\text{vap}}}{\partial \ln r} \right) \nonumber \\
&= \frac{1}{r} \frac{\partial}{\partial r} \left( r \Sigma_{\text{vap}} v_{\text{g},v} - D_{\text{Dr,g}} \Sigma_{\text{vap}} \frac{\partial \ln P_{\text{vap}}}{\partial \ln r} \right), \tag{69}
\end{align}
$$

where we assume that the change in the water vapor is more rapid
than the other quantities (i.e., we only consider the local change
in the water vapor and ∂P<sub>vap</sub>/∂r ≈ ∂Σ<sub>vap</sub>/∂r). As discussed in the
previous subsection (Sect. 3.1), P<sub>vap</sub> satisfies the local saturation
vapor (Eq. (48)) and

$$
\frac{\partial \ln P_{\text{vap}}}{\partial \ln r} = -\beta \frac{T_0}{T(r)}. \qquad (70)
$$

From Eq. (69), the diffusion dominates over the advection when

$$
r \Sigma_{\text{vap}} v_{\text{g},v} - D_{\text{Dr,g}} \Sigma_{\text{vap}} \frac{\partial \ln P_{\text{vap}}}{\partial \ln r} > 0, \quad (71)
$$

which is rewritten as

$$
c_s^2 \Omega_K^{-1} \Sigma_{\text{vap}} \left( -\frac{3}{2} \alpha_{\text{acc}} + \beta \frac{T_0}{T(r)} \alpha_{\text{Dr}} \right) > 0. \quad (72)
$$

Therefore, the diffusion dominates over the advection when

$$
\frac{\alpha_{\text{Dr}}}{\alpha_{\text{acc}}} > \frac{3}{2\beta} \frac{T(r)}{T_0} \approx 0.08 \times \left(\frac{\beta}{0.5}\right)^{-1} \left(\frac{T_0}{6062 \text{ K}}\right)^{-1} \left(\frac{T(r)}{170 \text{ K}}\right). \quad (73)
$$

The analytical prediction of Eq. (73) (with r = 2.4 au for β = 0.5) is shown by a gray dashed line in Fig. 3 and it shows a good consistency with 1D numerical simulations.

3.4. A quick summary of sublimation width of drifting pebbles

As discussed above, the sublimation width of icy drifting pebbles
is categorized by two regimes – the advection-dominated regime
(Sect. 3.1) and the diffusion-dominated regime (Sect. 3.2) –
which is a function of α<sub>Dr</sub>/α<sub>acc</sub>. In the lower and higher α<sub>Dr</sub>/α<sub>acc</sub>
regions, advection and radial diffusion are dominated, respec-
tively. Here, we provide a fitting function that smoothly connects
the two regimes as

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

where $X_-$ = $\log_{10}$(2) and $X_+$ = $\log_{10}$(0.1).

Equation (74) is shown by a black line in Fig. 3 and is used to
characterize $H_d(\Delta x_{\text{subl}})$ in Sect. 5.

3.5. Effects of sintering on the sublimation of pebbles

When sintering effects are considered (e.g., Saito & Sirono 2011;
Okuzumi et al. 2016), a sudden destruction of pebbles may occur
during the sublimation and the resultant sublimation width may
become smaller than those considered above (Sects. 3.1 and 3.2),
where a gradual decrease in size is expected. Pebbles would sub-
lim ate by a certain fraction to produce enough vapor that triggers

a sudden destruction of pebbles due to sintering effects. Thus,
the sublimation width may be reduced by some fraction from
the above arguments when sintering is taken into account. This
would lead to a higher midplane concentration of silicate dust
(smaller Δx<sub>subl</sub> leads to a smaller H<sub>d</sub>). Details of the effects
of sintering on the sublimation width of drifting pebbles are
beyond the scope of this paper. We leave this matter for future
investigations.

4. Scale height and maximum flux of pebbles

4.1. Pebble scale height regulated by KH instability

When the midplane concentration of pebbles/dust increases,
a vertical shear in the gas velocity between the pebble-
concentrated layer and the upper layer induces vorticity at the
interface between the two fluids, i.e., Kelvin-Helmholtz (KH)
instability (e.g., Sekiya 1998; Youdin & Shu 2002; Chiang
2008).

The azimuthal component of gas velocity (Schoonenberg &
Ormel 2017) is

$$
v_{g,\phi} = - \Lambda \frac{1 + \Lambda \tau_s^2}{1 + \Lambda^2 \tau_s^2} \eta v_K - Z \frac{\frac{1}{2} \Lambda^2 \tau_s}{1 + \Lambda^2 \tau_s^2} v_{g,v}^* \quad (75)
$$

$$
\approx - \frac{\Lambda}{1 + \Lambda^2 \tau_s^2} \left( (1 + \Lambda \tau_s^2) - \frac{3 Z \tau_s}{4 (1 + Z)} \alpha_{\text{acc}} C_{\eta}^{-1} \right) \eta v_K. \quad (76)
$$

Because $\tau_s \ll 1$, $\alpha_{\text{acc}} \ll 1$, and $|v_{g,\phi}| \gg |v_{g,r}|$, $v_{g,\phi}$ is approximated as

$$
v_{g,\phi} \approx -\Lambda \eta v_K = \frac{1}{1+Z} \eta v_K. \qquad (77)
$$

The Richardson number for pebbles at vertical height z is

$$
Ri = - \frac{g \left( \frac{d(\rho_g + \rho_p)}{dz} \right)}{\left( \rho_g + \rho_p \right) \left( \frac{dv_{g,\phi}}{dz} \right)^2}
$$

$$
\approx \frac{\Omega_K^2 z \times (\rho_g H_g^{-2} + \rho_p H_p^{-2}) z}{(\rho_g + \rho_p) \left( \frac{Z}{(1+Z)^3} \eta v_K z H_p^{-2} (1 - h_{p/g}^2) \right)^2} \quad (78)
$$

$$
\approx \frac{(1 + Z h_{p/g}^{-2}) \Omega_K^2 (z/H_g)^2}{\frac{Z^2}{(1+Z)^3} \eta^2 \Omega_K^2 (rz/H_p^2)^2 (1 - h_{p/g}^2)^2} \quad (79)
$$

$$
\approx \frac{(1 + Z)^3}{Z^2} C_{\eta}^{-2} \left( \frac{H_g}{r} \right)^{-2} (1 + Z h_{p/g}^{-2}) \left( \frac{h_{p/g}^2}{1 - h_{p/g}^2} \right)^2, \quad (80)
$$

where assumptions are made for ρ<sub>g</sub>(z) ∝ exp(-z<sup>2</sup>/2H<sub>g</sub><sup>2</sup>), ρ<sub>p</sub>(z) ∝
exp(-z<sup>2</sup>/2H<sub>p</sub><sup>2</sup>), and g ≈ Ω<sub>K</sub><sup>2</sup>z. We define C as

$$
C \equiv Ri \frac{Z^2}{(1+Z)^3} C_{\eta}^2 \left(\frac{H_g}{r}\right)^2 . \qquad (81)
$$

Defining $X \equiv h_{p/g}^2$ and using $C$ and $X$, we can rewrite Eq. (80) as

$$
C = \frac{(X + Z)X}{(1 - X)^2}, \qquad (82)
$$

which is equivalently written as

$$
(1 - C) X^2 + (Z + 2C) X - C = 0. \qquad (83)
$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Midplane pebble-to-gas density ratio ($Z = \rho_p/\rho_g$) as a function of $\alpha_{Dz}$ and $F_{p/g}$ (at $r=5$ au and $\tau_{s,p}=0.125$). Top panels: cases where the back-reaction onto the gas motion is neglected (i.e., $v_g = v_{g,v}$). Bottom panels: cases where the back-reaction to the gas motion is included, following Eq. (16). The color contours are obtained by directly solving Eq. (86). Left and right two panels: cases where $\alpha_{acc} = 10^{-3}$ and $\alpha_{acc} = 10^{-2}$ with/without a KH instability ($Ri=0.5$), respectively. The red-colored regions indicate where a steady-state solution is not found (i.e., $\rho_p/\rho_g > 100$ and it keeps increasing; the “no-drift (ND)” region). The white (Eq. (91)) and black (Eq. (95) with $Ri=0.5$) dashed lines are analytically derived critical $F_{p/g,\text{crit}}$ above which no steady-state are found (the “no-drift (ND)” region). Analytical predictions (dashed lines) well predict critical boundaries of the ND region. Including a KH instability reduces the ND region as the minimum scale height of pebbles is regulated by a KH instability. Due to the back-reaction to the gas motion that slows down the radial velocity of the gas, the midplane concentration of pebbles as well as the ND region (red-colored region) are shifted toward slightly larger $F_{p/g}$ values.</figcaption></figure>

As $Z \gg C$ by definition and $X \ll 1$, the solution is $X \approx C/Z$, that is, the scale height of pebbles regulated by a KH instability is written as

$$h_{p/g,KH} \approx Ri^{1/2} \frac{Z^{1/2}}{(1+Z)^{3/2}} C_{\eta} \left(\frac{H_g}{r}\right) \quad (84)$$

$$\approx 0.039 \left(\frac{Ri}{0.5}\right)^{1/2} \frac{Z^{1/2}}{(1+Z)^{3/2}} \left(\frac{C_{\eta}}{11/8}\right) \left(\frac{H_g/r}{0.04}\right), \quad (85)$$

where Eq. (85) has a peak value at $Z=1/2$. We note that Eq. (85) with $Ri=0.5$ and $Z=0.5$ ($h_{p/g} \equiv H_p/H_g$) is equivalent, within a factor order unity, to a classical description of $H_p \propto \eta r$ (e.g., Chiang & Youdin 2010). Although $Ri$ has a dependence on $Z$ and $Z > 1$ could lead to SI and/or cusp (Sekiya 1998; Youdin & Shu 2002), the average $Ri \sim 0.75$ is reported for $Z < 1$ (or $Z_{\Sigma} < 0.01$) (Gerbig et al. 2020). Thus, our simple prescription is valid for $Z < 0$ (this is where we use the above prescription to derive a critical $H_p$). Further detailed discussion with numerical simulations is beyond the scope of this paper.

We note also that we estimated $h_{p/g}$, assuming a Gaussian distribution for heights of the particles, while the precise height distribution function has a cusp for $Z \gtrsim 1$ and the cusp becomes sharper as $Z$ increases (Sekiya 1998; Chiang 2008). The height of the cusp ($z_{\text{max}}$) at which the particle distribution is truncated becomes constant in the limit of $Z \gg 1$ (Chiang 2008), while Eq. (85) shows $h_{p/g,KH} \propto Z^{-1}$ for $Z \gg 1$. On the other hand, the precise distribution is more similar to the Gaussian distribution for $Z \ll 1$, Actually, in this case, $z_{\text{max}}$ and $H_p$ given by Eq. (85) are similar (only differ by $\sqrt{2}$). Because we define $H_p$ as a root mean square of heights of the particle distribution, the derivation assuming a Gaussian distribution is more appropriate than the estimation of $H_p = z_{\text{max}}$.

## 4.2. “No-drift (ND)” region

Drift-BKR reduces the radial drift velocity of pebbles as a local concentration of pebbles is elevated (see also Gonzalez et al. 2017). The midplane concentration of pebbles increases with increasing $F_{p/g}$ and decreasing $H_p$.

For sufficiently large $F_{p/g}$ and small $H_p$ (i.e., small $\alpha_{Dz}$) values, the drift velocity of pebbles can decrease in a self-induced manner due to Drift-BKR, eventually leading to a complete stop of the radial motion of pebbles. In such a case, pebbles no longer reach the snow line from the outer region of the disk. We call this parameter regime as the “no-drift (ND)” region. As shown below, we explain that the ND region appears at an arbitrary choice of the radial distance, i.e., irrespective of the vicinity of the snow line.

Below, we discuss the concentration of icy pebbles in the disk midplane considering Drift-BKR onto pebbles and gas. We numerically and analytically derive the degree of the midplane concentration of pebbles. The results are also used for our initial conditions at the outer boundary of 1D simulations. We identify the ND region (within the $F_{p/g} - \alpha_{acc} - \alpha_{Dz}$ space).

The concentration of pebbles at the midplane is written as

$$Z \equiv \frac{\rho_p}{\rho_g} = \frac{\sum_p h_{p/g}^{-1}}{\sum_g h_{p/g}^{-1}} = \frac{v_g}{v_p} h_{p/g}^{-1} F_{p/g}, \quad (86)$$

where $v_g$, $v_p$, and $h_{p/g}$ are functions of $\Lambda(Z)$ and $Z = \rho_p/\rho_g$ (Eqs. (16), (17), (35), and (10)). We numerically find a solution (i.e., solve for $Z$ with $\tau_{s,p}=0.1$) for Eq. (86) for a given $F_{p/g}$ and $\alpha_{Dz}$. The color contours in Fig. 4 show the numerically obtained $\rho_p/\rho_g$ values. The red-colored regions correspond to the parameters for which no steady-state solution is found (numerically, this corresponds to $\rho_p/\rho_g \gg 100$), defined as the ND region for which pebbles no longer reach the snow line. We found that the boundaries of the “no-drift” regimes are characterized by a
critical value of $\rho_p/\rho_g = 1$ (i.e., $Z_{\text{cri}} = 1$; see the color contours in Fig. 4). The dashed lines in Fig. 4 are analytically derived critical $F_{\text{p,g,crit}}$ above which no steady-state solution is found for $\rho_p/\rho_g$ (i.e., no solution for Z). The analytical estimations well predict the direct solutions (the color contours). Below, we discuss that $F_{\text{p,g,crit}}$ can be divided into two different regimes.

The first regime is when $H_{\text{p,tur}} > H_{\text{p,KH,max}}$ where $H_{\text{p,KH,max}}$ is the maximum scale height of pebbles regulated by a KH-instability (Sect. 4.1). In this case, using Eqs. (16) and (17), the vertically averaged metallicity of pebbles is given as

$$Z_{\Sigma} = \frac{\Sigma_{\text{p}}}{\Sigma_{\text{g}}} = F_{\text{p/g}} \times \frac{v_{\text{g}}}{v_{\text{p}}} \qquad (87) \\ \approx F_{\text{p/g}}$$

$$\times \frac{\left(\frac{\Lambda}{1+\Lambda^2 \tau_{s,p}^2}\right) \left[ \left(\frac{2\tau_{s,p}}{1+Z} + \frac{3\alpha_{\text{acc}}}{2C_{\eta}}\right) Z h_{\text{p/g}} - \frac{3\alpha_{\text{acc}}}{2C_{\eta}} \left(\frac{1+\Lambda^2 \tau_{s,p}^2}{\Lambda}\right) \right] \eta v_{\text{K}}}{-\left(\frac{\Lambda}{1+\Lambda^2 \tau_{s,p}^2}\right) \left(2\Lambda \tau_{s,p} + \frac{3\alpha_{\text{acc}}}{2C_{\eta}}\right) \eta v_{\text{K}}} \qquad (88)$$

$$\approx \frac{3\alpha_{\text{acc}}F_{\text{p/g}}}{4\tau_{s,p}C_{\eta}\Lambda^2} = \frac{3\alpha_{\text{acc}}F_{\text{p/g}}}{4\tau_{s,p}C_{\eta}} (1 + h_{\text{p/g}}^{-1}Z_{\Sigma})^2, \qquad (89)$$

where approximation is made for $\alpha_{\text{acc}} \ll \tau_{s,p} \ll 1$, $Z \ll 1$, $\Lambda \approx 1$ (because $\Lambda = 1 \to 0$ as pile-up proceeds), and we use Eq. (11). Solving this equation gives

$$Z_{\Sigma} = \frac{(1 - 2ab) \pm \sqrt{1 - 4ab}}{2ab^2}, \qquad (90)$$

where $a = \frac{3\alpha_{\text{acc}}F_{\text{p/g}}}{4\tau_{s,p}C_{\eta}}$ and $b = h_{\text{p/g}}^{-1}$, respectively. $Z_{\Sigma}$ has a real solution when $1 - 4ab > 0$ and a critical $F_{\text{p/g,crit1}}$ for the first regime is given as

$$F_{\text{p/g,crit1}} = \frac{C_{\eta}}{3} \frac{\tau_{s,p}}{\alpha_{\text{acc}}} h_{\text{p/g}} \approx \frac{(\alpha_{\text{Dz}} \tau_{s,p})^{1/2} C_{\eta}}{3\alpha_{\text{acc}}} \qquad (91) \\ \approx 0.15 \times \left(\frac{\alpha_{\text{acc}}}{10^{-2}}\right)^{-1} \left(\frac{\alpha_{\text{Dz}}}{10^{-4}}\right)^{1/2} \left(\frac{\tau_{s,p}}{0.1}\right)^{1/2} \left(\frac{C_{\eta}}{11/8}\right), \qquad (92) \\ \text{[for } H_{\text{p,tur}} > H_{\text{p,KH,max}} \text{]}$$

where $h_{\text{p/g}} \approx (\tau_{s,p}/(\alpha_{\text{Dz}}(1+Z)^{-K}))^{-1/2} \approx (\tau_{s,p}/\alpha_{\text{Dz}})^{-1/2}$ for $Z \ll 1$ and $\alpha_{\text{Dz}} \ll \tau_{s,p}$ (Eq. (35)). The white dashed line in Fig. 4 shows Eq. (91) and it generally shows a good consistency with the direct solutions of Eq. (86) (i.e., color contours in Fig. 4)$^4$.

The second regime appears when a KH instability plays a role, which corresponds to $H_{\text{p,tur}} < H_{\text{p,KH,max}}$, where $H_{\text{p,KH,max}}$ is the maximum scale height of pebbles regulated by a KH instability (see Sect. 4.1). Thus, in this case, the scale height of pebbles is described by $H_{\text{p,KH}}$, which is independent of $\alpha_{\text{Dz}}$. As $Z = \rho_p/\rho_g$ increases from zero value, $H_{\text{p,KH}}$ initially increases and reaches its maximum $H_{\text{p,KH,max}}$ at $Z = 1/2$ (see Sect. 4.1). As Z further increases, $H_{\text{p,KH}}$ decreases. Using $Z_{\text{cri}} = 1$ as a critical value, the critical $F_{\text{p/g}}$ in the second regime $F_{\text{p/g,crit2}}$ adopts $H_{\text{p,KH}}^{Z=1}$ ($H_{\text{p,KH}}$ with $Z_{\text{cri}} = 1$). $Z_{\text{cri}}$ is given as

$$Z_{\text{cri}} = \frac{\rho_{\text{p}}}{\rho_{\text{g}}} \left|_{H_{\text{p,KH}}^{Z=1}} = \frac{\Sigma_{\text{p}}}{\Sigma_{\text{g}}} \left(h_{\text{p/g,KH}}^{Z=1}\right)^{-1} \qquad (93)$$

$$\approx \frac{3\alpha_{\text{acc}}F_{\text{p/g,crit2}}}{4\tau_{s,p}C_{\eta}(\Lambda_{\text{KH}}^{Z=1})^2} (h_{\text{p/g,KH}}^{Z=1})^{-1}, \qquad (94)$$

where $h_{\text{p/g,KH}}^{Z=1} \equiv H_{\text{p,KH}}^{Z=1}/H_g$ and approximation is made for $\alpha_{\text{acc}} \ll \tau_{s,p} \ll 1$ and $\Lambda \approx 1$. $\Lambda_{\text{KH}}^{Z=1}$ is $\Lambda$ using $H_{\text{p,KH}}^{Z=1}$. Here, $Z = 1$ and $\Lambda_{\text{KH}}^{Z=1} = 1/2$. Thus, the critical $F_{\text{p/g}}$ for the second regime $F_{\text{p/g,crit2}}$ with $C_{\eta} = 11/8$ for $T \propto r^{-1/2}$ is given as

$$F_{\text{p/g,crit2}} = \frac{C_{\eta}}{3} \left(\frac{\tau_{s,p}}{\alpha_{\text{acc}}}\right) h_{\text{p/g,KH}}^{Z=1} = \frac{11}{24} \left(\frac{\tau_{s,p}}{\alpha_{\text{acc}}}\right) h_{\text{p/g,KH}}^{Z=1} \qquad (95)$$

$$\approx 0.06 \times \left(\frac{\alpha_{\text{acc}}}{10^{-2}}\right)^{-1} \left(\frac{\tau_{s,p}}{0.1}\right) \left(\frac{Ri}{0.5}\right)^{1/2} \left(\frac{H_g}{0.04}\right), \qquad (96) \\ \text{[for } H_{\text{p,tur}} < H_{\text{p,KH,max}} \text{]}$$

which is independent on $\alpha_{\text{Dz}}$. Equation (95) is plotted by the black dashed line in Fig. 4 and it shows a good consistency with the direct solution of Eq. (86). Including a KH instability (i.e., the second regime and Eq. (95)) prevents the scale height of pebbles becomes further smaller as $\alpha_{\text{Dz}}$ becomes small, reducing the parameter range of the “no-drift” region (see Fig. 4 for the cases with/without a KH instability).

The above arguments are irrespective of the vicinity of the snow line and the ND mode could occur at an arbitrary radial location to the central star. The application of the “no-drift” mechanism in protoplanetary disks is discussed in Hyodo et al. (2021).

## 5. Pile-ups around the snow line: dust or pebbles?

In this section, we finally show the overall results of our full 1D simulations that include sublimation of ice, the release of silicate dust, and the recycling processes (water vapor as well as silicate dust via the recondensation and sticking onto pebbles outside the snow line – see Fig. 1). We adopt a realistic scale height of silicate dust obtained from Paper I. We also consider the effects of a KH instability for the scale height of pebbles (Sect. 4.1), which was not included in the previous works. We include Drift-BKR and Diff-BKR onto the motions of pebbles and silicate dust. We choose $K = 1$ for Diff-BKR. We change $F_{\text{p/g}}$, $\alpha_{\text{acc}}$, and $\alpha_{\text{Dr}} = \alpha_{\text{Dz}}$ as parameters. For simplicity and to have the same setting as Paper I, we fix the gas structure as an unperturbed background disk (i.e., $v_g = v_{g,v}$).

Figure 5 shows the maximum midplane solid-to-gas ratio around the water snow line in the $F_{\text{p/g}}-\alpha_{\text{acc}}-\alpha_{\text{Dz}}(=\alpha_{\text{Dr}})$ space. The background color contours correspond to the results obtained for silicate dust inside the snow line in Paper I. The red-colored area indicates runaway pile-up of silicate dust (i.e., no steady-state and an indication of planetesimal formation via direct gravitational collapse). The dependence on $\Delta x_{\text{sub}}$ is discussed in Appendix C.

Points in Fig. 5 represent the results of our 1D simulations. The left and right numbers in the parenthesis represent the maximum $\rho_p/\rho_g$ (pebbles) and $\rho_d/\rho_g$ (silicate dust), respectively. When a runaway pile-up occurs in our 1D simulations, it is labeled either by “Peb-RP” for pebbles or by “Sil-RP” for silicate dust. The numbers below the parenthesis in Fig. 5 show the rock fraction in pebbles. For a combination of a small $\alpha_{\text{Dz}}(=\alpha_{\text{Dr}})/\alpha_{\text{acc}}$ and a large $F_{\text{p/g}}$ as an initial setting (as we set pebbles initially exist only at the outer boundary), pebbles do not drift from the outer boundary as a consequence of continuous slowing down of the radial drift velocity of pebbles due to Drift-BKR (we label

<footnote>$^4$ Even when the back-reaction to the gas motion is neglected (i.e., $v_g = v_{g,v}$ as Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019), $F_{\text{p/g,crit1}}$ is approximated to $\sim \frac{C_{\eta}}{3} \frac{\tau_{s,p}}{\alpha_{\text{acc}}} h_{\text{p/g}}$ (i.e., replace the denominator of Eq. (88) with $v_{g,v}$ and approximation is made for $\alpha_{\text{acc}} \ll \tau_{s,p} \ll 1$, $Z \ll 1$, $\Lambda \approx 1$).</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 5.** Zones of pile-up of pebbles and/or silicate dust around the snow line as a function of $F_{p/g}$ and $\alpha_{Dz} (= \alpha_{Dr})/\alpha_{acc}$ (case of a variable $\Delta x_{subl}$; Eq. (74) in Sect. 3.4). The numbers in the parenthesis represent the maximum $\rho_p/\rho_g$ (left) and $\rho_d/\rho_g$ (right) obtained from our full 1D simulations, respectively. The numbers below the parenthesis represent the rock fraction in the surface density of pebbles at the maximum $\rho_p/\rho_g$ (when a runaway pile-up of silicate dust occurs, the fraction is not shown because the system would gravitationally collapse to form 100% rocky planetesimals). The runaway pile-up of pebbles (labeled as “Peb-RP”) or that of silicate dust (labeled as “Sil-RP”) occurs around the snow line, depending on the combinations of $F_{p/g}$ and $\alpha_{Dz} (= \alpha_{Dr})$. The color contour shows an analytical prediction of the maximum $\rho_d/\rho_g$ just inside the snow line obtained from Paper I. The black dashed line indicates analytically derived critical $F_{p/g}$ above which the “no-drift (ND)” takes place (labeled by “ND” for the results of 1D simulations in the panels) at the outer boundary ($r = 5$ au and $\tau_{s,p} = 0.125$; Sect. 4.2). Drift-BKR and Diff-BKR onto the motions of pebbles and silicate dust are included with $K = 1$. Left and right panels: cases of $\alpha_{acc} = 10^{-2}$ and $\alpha_{acc} = 10^{-3}$, respectively. In the case of $\alpha_{acc} = 10^{-2}$ the runaway pile-up of silicate dust is inhibited because the required parameter regime is overlapped by the “no-drift” region (i.e., pebbles do not reach the snow line). In the case of $\alpha_{acc} = 10^{-3}$, the runaway pile-ups of both pebbles and silicate dust occur.

these parameters by “ND” as the “no-drift” region in Fig. 5 – see also Sect. 4.2).

In the following subsections, we discuss the comparison to previous works (Sect. 5.1), the detailed processes of pile-up of silicate dust (Sect. 5.2), that of pebbles (Sect. 5.3), runaway pile-ups of pebbles and silicate dust (Sect. 5.4), rock-to-ice ratio within the pile-up region (Sect. 5.5), and implication for planet formation (Sect. 5.6).

## 5.1. Comparison to previous works

In most previous studies of the pile-up of silicate dust and pebbles around the snow line, $\alpha_{acc} = \alpha_{Dr} = \alpha_{Dz}$ was assumed (Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017). Figure 5 indicates that such parameters do not lead to a runaway pile-up of silicate dust inside the snow line, while a pile-up of pebbles outside the snow line could satisfy a condition for SI to take place for sufficiently large pebble mass flux of $F_{p/g} > 0.8$ (Fig. 5). Previous studies (Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017) assumed $K = 0$ and a runaway pile-up of pebbles was not reported (only a steady-state is reached with $\rho_p/\rho_g > 1$ for $F_{p/g} > 0.8$), while it could take place when $K \neq 0$ with Drift-BKR (see also Hyodo et al. 2019).

## 5.2. Pile-up of silicate dust inside the snow line

As discussed in Paper I, silicate dust can pile up just inside the snow line (Fig. 1; see also Estrada et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019; Gárate et al. 2020). Our 1D simulations include the recycling process via sticking of the diffused dust onto icy pebbles beyond the snow line. It conserves the inward mass flux of silicate dust $F_{d/g}$ across the snow line (all the diffused silicate dust eventually comes back to the snow line with

pebbles), which indicates $F_{d/g} = f_{d/p} \times F_{p/g}$ (here, $f_{d/p} = 0.5$ is the silicate fraction in the original pebbles at the outer boundary) is independent on $\alpha_{acc}$ and $\alpha_{Dr}$. The vertically integrated metallicity of silicate dust can be written as $Z_{\Sigma_d} \propto F_{d/g} \times v_g/v_d$ and is independent on $\alpha_{acc}$ as $v_d \simeq v_g$. The midplane solid-to-gas ratio is written as $\rho_d/\rho_g = Z_{\Sigma_d} h_{d/g}^{-1} (h_{d/g} \equiv H_d/H_g)$. At the snow line, $h_{d/g}$ is regulated by $\alpha_{Dz}/\alpha_{acc}$ (see Eq. (36) for $\alpha_{Dz} = \alpha_{Dr}$). Thus, both $\alpha_{acc} = 10^{-2}$ and $\alpha_{acc} = 10^{-3}$ cases have almost the same $\rho_d/\rho_g$ for the same $\alpha_{Dz}/\alpha_{acc}$ and $F_{p/g}$ when $\rho_d/\rho_g \ll 1$ (the right numbers in the parenthesis in Fig. 5)$^5$.

Our 1D simulations (the right numbers in the parenthesis in Fig. 5) and the results of an analytical formula that is calibrated by the Monte Carlo simulations (a color contour in Fig. 5 obtained from Paper I) show good consistency when pebble pile-up is not significant. Paper I neglected the recycling of water vapor onto pebbles, but it does not significantly affect the pile-up of silicate dust inside the snow line as the mass flux of silicate dust is a critical parameter for the pile-up.

## 5.3. Pile-up of pebbles outside the snow line

In this subsection, we discuss how pebbles pile up just outside the snow line and its dependence on $\alpha_{acc}$. Outside the snow line, the pile-up of icy pebbles is locally enhanced due to two distinct recycling processes: the recondensation of water vapor and sticking of silicate dust (Fig. 1; see also Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019; Gárate et al. 2020). The efficiency of the recycling process (i.e., $\alpha_{Dr}$) is a critical parameter to regulate the pile-up of icy pebbles. Thus, unlike the case of silicate dust (Sect. 5.2), the value $\rho_p/\rho_g$ is different for the same

<footnote>$^5$ When $\rho_d/\rho_g > 1$, the back-reaction changes other quantities, such as the scale height of pebbles (Eq. (35)).</footnote>
$α_{Dr}/α_{acc} = (α_{Dz}/α_{acc})$ and $F_{p/g}$ at a different $α_{acc}$ (see the left and right panels in Fig. 5).

The surface density of pebbles does not strongly depend on $α_{acc}$ as pebbles are only partially coupled to the gas flow ($τ_{s,p} \sim 0.1$; Okuzumi et al. 2012; Ida & Guillot 2016) and the radial velocity of pebbles is mainly dominated by the gas drag (the first term in Eq. (17)), while the surface density of the gas becomes smaller for larger $α_{acc}$ as $Σ_g \propto 1/v_g \propto 1/α_{acc}$. This indicates that $ρ_p/ρ_g$ is $\propto α_{acc}$ when the recycling processes are negligible.

As $F_{p/g}$ increases or/and $α_{Dz}/α_{acc}$ decreases, the pile-up of pebbles tends to increase as the disk metallicity increases ($Σ_p \propto F_{p/g}$) and the scale height of pebbles decreases (Eqs. (35) and (36)). However, when $α_{Dr}/α_{acc} \approx 1$, an efficient recycling of water vapor and dust onto pebbles can change this dependence for some parameters (for example, the cases of $α_{acc} = α_{Dz} = 10^{-3}$). Comparing the maximum $ρ_d/ρ_g$ and $ρ_p/ρ_g$, the pile-up of pebbles is more significant than that of silicate dust when $F_{p/g}$ is sufficiently large and when the diffusion is an efficient process (i.e., $α_{Dr}/α_{acc} \sim 1$ cases) because the recycling of the vapor and the dust onto pebbles efficiently enhances the pile-up of pebbles (Hyodo et al. 2019). The pile-up of silicate dust is more prominent when the diffusivity is weak (i.e., $α_{Dz}/α_{acc} \ll 1$ cases).

## 5.4. Conditions for preferential pile-up of silicate dust versus pebbles

Finally, we aim to understand under which conditions the runaway pile-ups of silicate dust or/and pebbles are favored around the snow line.

The black lines in Fig. 5 indicate where pebbles cannot drift toward the snow line, corresponding to the ND region (see Sect. 4.2; $r=5$ au and $R_i=0.5$ are used here). In the case of $α_{acc} = 10^{-2}$, the red-colored region is located above the black line, which implies that parameters in the $F_{p/g} - α_{acc} - α_{Dz}$ space that would yield the runaway pile-up of silicate dust found in Paper I in fact lie in the forbidden ND region. On the other hand, for $α_{acc} = 10^{-3}$, the pile-up of silicate dust inside the snow line for $α_{Dz} = α_{Dr} < 10^{-5}$ and $0.1 < F_{p/g} < 0.8$ found Paper I is outside of the ND region and confirmed by our simulations.

Generally, runaway pile-ups around the snow line occur over a broader range of parameters for $α_{acc} = 10^{-3}$ than for $α_{acc} = 10^{-2}$. Runaway pile-up of silicate dust is favored for a small $α_{Dz}/α_{acc} < 10^{-3}$ (labeled by “Sil-RP”), while runaway pile-up of pebbles (labeled by “Peb-RP”) only occurs for a limited range of parameters (large $F_{p/g}$ and $α_{Dz}/α_{acc} \sim 1$).

The results in Fig. 5 thus indicate that planetesimal formation around the snow line potentially occurs via a direct gravitational collapse of solids, to form icy or rocky planetesimal inside or outside the snow line, respectively, depending on the parameters of $F_{p/g}$ and $α_{Dr}/α_{acc}$. However, the formation efficiencies of these two kinds of planetesimals strongly depend on the disk conditions. This emphasizes that detailed disk evolution models including a realistic treatment of pebble growth are required to understand planetesimal formation around the snow line.

## 5.5. Rock-to-ice ratio within the pile-up region

Around the snow line, the rock-to-ice mixing ratio can be modified due to sublimation, condensation, and the recycling effects (Hyodo et al. 2019). Our numerical simulations show that pebbles just outside the snow line contain both ice and rock (the numbers below parenthesis in Fig. 5 indicate resulting rock

fractions). The results indicate that $α_{Dr}/α_{acc} = 10^{-1}$ cases tend to most efficiently increase the rock fraction from the original value (the 1:1 initial rock-ice ratio; Sect. 2) for a fixed $F_{p/g}$.

For a higher $α_{Dr}/α_{acc}$, the scale height of silicate dust becomes larger (Eq. (36)), which leads to a less efficient sticking of silicate dust onto icy pebbles (its efficiency $\propto 1/H_d$), leading to a lower rock fraction in pebbles. However, a higher $α_{Dr}/α_{acc}$ yields a more efficient outward diffusion of silicate dust released in the inner region, leading to a higher surface density of silicate dust outside the snow line, hence enhancing the silicate fraction in pebbles (the sticking efficiency $\propto Σ_d$).

Diff-BKR and Drift-BKR also play critical roles as these effects modify the diffusivity and radial motion of solids. As a result of a combination of these different effects, the numerical results show that the rock fraction in pebbles, generally increases from the initial value of 50 wt.% up to ~80 wt.%, depending on $α_{Dr}/α_{acc}$ and $F_{p/g}$ (Fig. 5). When a runaway pile-up of pebbles occurs (labeled by Peb-RP in Fig. 5), the rock fraction becomes ~50 wt.%. This is because the runaway pile-up occurs by accreting pebbles that drift from the outer region.

Thus, forming planetesimals from pebbles (by GI or SI) by runaway pile-up outside the snow line leads to a rock fraction being similar to that of the original pebbles formed in the outer region. The threshold $ρ_p/ρ_g$ for the occurrence of a streaming instability should be close to unity, within an order of magnitude (Carrera et al. 2015; Yang et al. 2017). As shown in Fig. 5, lower values of this threshold could lead to the formation of rock-rich planetesimals even from the pile-up of pebbles outside of the ice line. This is for example the case for $α_{acc} = 10^{-3}$ and $F_{p/g} = 0.3$. This may explain the high dust-to-ice ratio measured in comet 67P (O’Rourke et al. 2020). Planetesimals formed from the pile-up of silicate dust inside of the snow line should be purely rocky.

## 5.6. Implication for planet formation

The structure and evolution of protoplanetary disks remain poorly constrained. We propose with Fig. 6 an example scenario that would lead to the formation of both ice-rich and rock-rich planetesimals as seen in the solar system.

The evolution of the disks is characterized by a decrease of the accretion rate, a decrease of viscous heating and therefore an inward migration of the snow line (Oka et al. 2011). In the same time, $F_{p/g}$ could decrease (Ida et al. 2016). The inner region of the disk may have an MRI-inactive “dead zone” near the disk midplane where turbulence and correspondingly diffusion are weak (Fig. 6; Gammie 1996). The radial boundary between the active zone in the outer disk region ($α_{Dr} = α_{Dz})/α_{acc} \sim 1$) and the dead zone in the inner disk region ($α_{Dr} = α_{Dz})/α_{acc} \ll 1$) may smoothly change, that is, $α_{Dr} = α_{Dz})/α_{acc}$ may gradually become smaller as the distance to the central star becomes smaller (e.g., Bai et al. 2016; Mori et al. 2017).

Our results (see Fig. 5) imply that in the early stage of the disk evolution when the snow line is located in the active outer region, a pile-up of icy pebbles leading to the formation of icy planetesimals by GI or SI may take place. As shown in Fig. 6, late in the evolution, the migration of the snow line into the dead inner region would lead to the formation of rocky planetesimals. In-between, in an intermediate phase, the planetesimal formation around the snow line would be suppressed because of a combination of $F_{p/g}$ and $α_{Dz}/α_{acc}$ that are insufficient to lead to a pile-up of solids in the disk midplane.

As the observed structure and theoretically modeled evolutionary paths of protoplanetary disks are diverse, following a distinct evolutionary path in the $α_{acc}-α_{Dr}-α_{Dz}-F_{p/g}$ space as well
<figure><img src="image_8.png" /></figure>

**Fig. 6.** Illustration showing a possible planet formation scenario. Here we envision a disk containing an MRI-inactive dead zone near the midplane ($\alpha_{Dr} = \alpha_{Dz}/\alpha_{acc} \ll 1$), while the outer disk is active ($\alpha_{Dr} = \alpha_{Dz}/\alpha_{acc} \sim 1$). We also envision an initially high $F_{p/g}$ value that becomes smaller with time, while the snow line moves inward. Our results of planetesimal formation around the snow line imply that the icy pebble pile-up (i.e., icy planetesimals formation by GI and/or SI; panel A) would be favored in the outer region ($\alpha_{Dz}/\alpha_{acc} \sim 1$) in the early stage of the evolution, while the formation of rocky planetesimals by GI of silicate dust would be favored in the later stage when the snow line has reached the dead zone ($\alpha_{Dz}/\alpha_{acc} \ll 1$; panel C). In the intermediate phase ($\alpha_{Dz}/\alpha_{acc} < 1$ but not $\alpha_{Dz}/\alpha_{acc} \ll 1$) and $F_{p/g} < 0.6$, SI/GI of pebbles and GI of silicate dust would not be expected to occur (panel B). A diversity of evolutionary paths for protoplanetary disks would produce a diversity of planetary systems.

as that of the snow line could produce a diversity of the outcomes of the planetesimal formation around the snow line. A preferential formation of planetesimals only around the snow line would produce a localized narrow ring distribution of planetesimals.

**6. Summary**

In this paper, to study the pile-up of solids around the snow line, we incorporated a realistic $H_d$ (obtained from Paper I) and the KH-instability-considered $H_p$ (Sect. 4.1), for the first time, into a 1D diffusion-advection code that includes the back-reactions to radial drift and diffusion of icy pebbles and silicate dust. The code takes into account ice sublimation, the release of silicate dust, and their recycling through the recondensation and sticking onto pebbles. We studied a much wider range of disk parameters (in the $F_{p/g} - \alpha_{acc} - \alpha_{Dr} = \alpha_{Dz}$) space) than previous studies. Our main findings are as follows.

First, we derived the sublimation width of drifting icy pebbles, a critical parameter to regulate the scale height of silicate dust (Sect. 3). We also derived the scale height of pebbles, including the effects of KH instabilities (Sect. 4.1).

Second, using analytical arguments, we identified a parameter regime in the $F_{p/g} - \alpha_{acc} - \alpha_{Dz}$ space, in which pebbles cannot reach the snow line by stopping their radial drift due to Drift-BKR in a self-induced manner (the “no-drift” region; Sect. 4.2).

Third, our 1D simulations showed that $\alpha_{acc} = 10^{-3}$ case is more favorable for the pile-up of solids around the snow line than that for $\alpha_{acc} = 10^{-2}$. The “no-drift” regime entirely covers parameter space of the runaway pile-up of silicate dust when $\alpha_{acc} = 10^{-2}$ case (left panel in Fig. 5), preventing a runaway

pile-up of silicate dust inside the snow line. In contrast, in the $\alpha_{acc} = 10^{-3}$ case, large $F_{p/g}$ values are allowed and could lead to a runaway pile-up of silicate dust inside the snow line (right panel in Fig. 5). In the case of $\alpha_{acc} = 10^{-3}$, both silicate dust and pebbles could experience runaway pile-ups inside and outside the snow line, forming rocky and icy planetesimals, respectively.

Forth, pebbles just outside the snow line is a rock-ice mixture (Sect. 5.5). When a runaway pile-up of pebbles occurs just outside the snow line, the local rock-to-ice ratio becomes similar to that of pebbles formed at the outer region – i.e., the resultant planetesimal composition would be water-bearing materials, if planetesimals ultimately formed from such pebbles. In contrast, almost pure rocky planetesimals formed from silicate dust, if GI occurs inside the snow line.

Lastly, we discussed the implication of the runaway pile-ups around the snow line to planet formation (Sect. 5.6). A diversity of outcomes in terms of planetesimal formation around the snow line would occur for a diversity of protoplanetary disks. Detailed prescriptions of disk evolution and pebble growth are necessary.

*Acknowledgements.* We thank Dr. Chao-Chin Yang for discussions. We thank Dr. Beibei Liu for his constructive comments that improved the manuscript. R.H. was supported by JSPS Kakenhi JP17J01269 and 18K13600. R.H. also acknowledges JAXA's International Top Young program. T.G. was partially supported by a JSPS Long Term Fellowship at the University of Tokyo. S.I. was supported by MEXT Kakenhi 18H05438. 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**

Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14
Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76
R. Hyodo et al.: Planetesimal formation around the snow line. II. Dust or pebbles?

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

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

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

Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50

Chiang, E. 2008, ApJ, 675, 1549

Chiang, E., & Youdin, A. N. 2010, Ann. Rev. Earth Planet. Sci., 38, 493

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

Desch, S. J., Estrada, P. R., Kalyaan, A., & Cuzzi, J. N. 2017, ApJ, 840, 86

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

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

Estrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200

Gammie, C. F. 1996, ApJ, 457, 355

Gárate, M., Birnstiel, T., Drążkowska, J., & Stammler, S. M. 2020, A&A, 635, A149

Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91

Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051

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

Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84

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., Ida, S., & Guillot, T. 2021, A&A, 645, L9

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

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

Ida, S., Guillot, T., Hyodo, R., Okuzumi, S., & Youdin, A. N. 2021, A&A, 646, A13

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

Lichtenegger, H. I. M., & Komle, N. I. 1991, Icarus, 90, 319

Liu, B., Lambrechts, M., Johansen, A., & Liu, F. 2019, A&A, 632, A7

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

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

Mori, S., Muranushi, T., Okuzumi, S., & Inutsuka, S.-i. 2017, ApJ, 849, 86

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

Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141

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

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

O'Rourke, L., Heinisch, P., Blum, J., et al. 2020, Nature, 586, 697

Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65

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

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

Sekiya, M. 1998, Icarus, 133, 298

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117

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

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

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

Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494

Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81
## Appendix A: Dependence of dust pile-up upon $H_d$

<figure><img src="image_9.png" /><figcaption>Fig. A.1. Scale height of silicate dust in a unit of that of the gas (top panels) and midplane dust-to-gas ratio (bottom panels) for different models of $H_d$ (Case of $F_{p/g} = 0.3$ and $\tau_{s,p} = 0.1$ at 3 au). Here, Drift-BKR and Diff-BKR for pebbles are neglected, whereas Drift-BKR and Diff-BKR ($K = 0$ or 1) for silicate dust are included. Left two panels: cases of $\alpha_{acc} = 10^{-2}$ and $\alpha_{Dr} = \alpha_{Dz} = 10^{-4}$ with $K = 0$ or 1, respectively. Right two panels: cases of $\alpha_{acc} = 10^{-2}$ and $\alpha_{Dr} = \alpha_{Dz} = 10^{-2}$ with $K = 0$ or 1, respectively. The gray lines are the case of Ida & Guillot (2016). The blue lines are the case of Schoonenberg & Ormel (2017). The red lines are the case of Hyodo et al. (2019) and the results show runaway pile-up when $\alpha_{acc} = 10^{-2}$ and $\alpha_{Dr} = \alpha_{Dz} = 10^{-4}$ with $K = 1$. The black lines are the case of this work (Eq. (36)). In the cases of the gray lines (i.e., Ida & Guillot 2016), pile-ups occur in a runaway fashion for the left two panels, and the time-evolutions are shown for $t = 1 \times 10^4$ yr, $t = 5 \times 10^4$ yr, and $t = 1.5 \times 10^5$ yr. The other cases reach a steady-state within $t = 5 \times 10^5$ yr. Here, the recycling of water vapor and silicate dust onto pebbles are not included to focus on the different models of $H_d$.</figcaption></figure>

Silicate dust is released during the sublimation of the drifting icy pebbles approaching the snow line. Different works have adopted different models of the scale height of silicate dust $H_d(r)$ in a 1D-radial accretion disk (e.g., Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019). The midplane density of silicate dust inside the snow line is calculated by $\rho_d = \Sigma_d / \sqrt{2\pi H_d}$, which is also used to describe the strength of the back-reaction of the gas onto their motion (e.g., Eq. (18)). Therefore, the description of $H_d$ is critical for the back-reaction and to evaluate the pile-up of dust in the midplane. In this section, we demonstrate how different models of $H_d$ affect the dust pile-up inside the snow line.

Ida & Guillot (2016) neglected the effects of the vertical stirring of the released dust and assumed that the dust keeps the same scale height as that of pebbles, that is, $H_d(r) = H_p(r_{snow})$, where $H_p(r_{snow})$ is the scale height of pebbles at the snow line. In contrast, Schoonenberg & Ormel (2017) assumed that the released silicate dust is instantaneously stirred to the vertical direction and has the same scale height as that of the gas, that is, $H_d(r) = H_g(r)$.

Because silicate dust is released from sublimating pebbles, Hyodo et al. (2019) modeled that the released silicate dust has the same scale height as that of pebbles at the snow line and is gradually stirred to the vertical direction as a function to the distance to the snow line as

$$H_d = \left( 1 + \frac{\tau_{s,p,snow}}{\alpha_{Dz}} \times e^{\frac{-\Delta t}{t_{mix}}} \right)^{-1/2} H_g \quad (A.1)$$

where a notation of “snow” indicates values at the snow line. $\Delta t = |(r - r_{snow})/v_d| \cdot t_{mix} = (H_g/l_{mfp})^2/\Omega_K$ is the diffusion/mixing

timescale to the vertical direction at the snow line and $l_{mfp} = \sqrt{\alpha_{Dz} H_g}$ is the mean free path of turbulent blobs. $\Delta t$ and $t_{mix}$ are calculated by using the physical values at the snow line and we use $\tau_{s,p,snow} = 0.1$.

Figure A.1 shows different models of the scale height of silicate dust and their resultant pile-ups of silicate dust, including a more realistic model used in this work (Eq. (36) which is originally derived in Paper I).

The model of Schoonenberg & Ormel (2017), $H_d = H_g$, underestimated pile-up of silicate dust when $\alpha_{Dz} \ll \alpha_{acc}$, while it correctly evaluated when $\alpha_{Dz} \approx \alpha_{acc}$. Ida & Guillot (2016), $H_d = H_{p,snow}$, overestimated the pile-up of silicate dust regardless of the value of $\alpha_{Dz}$. Hyodo et al. (2019) overestimated pile-up of silicate dust near the snow line.

## Appendix B: Particle size and Stokes number

The Stokes number of a particle $\tau_s(r)$ depends on the local gas structure and $\tau_s(r)$ is written by using the stopping time $t_s$ as $\tau_s(r) = t_s(r)\Omega_K(r)$. The stopping time represents the relaxation timescale of particle momentum through the gas drag. Here, we consider two regimes of the stopping time and $\tau_s(r)$ is given as

$$\tau_s(r) = t_s(r)\Omega_K(r) = \frac{\rho_p r_p}{\rho_g v_{th}} \Omega_K \quad \text{Epstein:} r_p < \frac{9}{4} l_{mfp} \quad (B.1)$$

$$\tau_s(r) = t_s(r)\Omega_K(r) = \frac{4\rho_p r_p^2}{9\rho_g v_{th} l_{mfp}} \Omega_K \quad \text{Stokes:} r_p > \frac{9}{4} l_{mfp} \quad (B.2)$$

where $v_{th}(r) = \sqrt{8/\pi c_s(r)}$ is the thermal velocity, $r_p$ is the particle radius, and $\rho_p = 1.5 \text{ g cm}^{-3}$ is the particle internal density, and
<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Particle size as a function of distance to the star. The red, green, and blue lines represent $r_p = \frac{9}{4}l_{mfp}$ for $\alpha_{acc} = 10^{-2}$, $10^{-3}$, and $10^{-4}$, respectively (Eq. (B.3)). The Stokes regime is where $r_p > \frac{9}{4}l_{mfp}$, while the Epstein regime is where $r_p < \frac{9}{4}l_{mfp}$. Four black and gray lines represent particle size whose $\tau_s = 0.3$, 0.1, 0.03, and 0.01 in the Stokes regime and the Epstein regime with $\alpha_{acc} = 10^{-2}$, respectively. Here, $M_g = 10^{-8}M_\odot$ yr$^{-1}$ and $T(r) = 150 (r/3au)^{-1/2}$ are used.</figcaption></figure>

$\rho_g = \Sigma_g / \sqrt{2\pi H_g}$ is the gas spatial density, respectively. The mean free path of the gas $l_{mfp}(r)$ is given as

$$l_{mfp}(r) = \frac{\mu_g m_{proton}}{\sqrt{2} \rho_g(r) \sigma_{mol}} \qquad (B.3)$$

where $\sigma_{mol} = 2.0 \times 10^{-15} \text{ cm}^2$ is the collisional cross section of the gas molecules. We note that $l_{mfp}$ depends on the local gas structure and is a function of $\alpha_{acc}$.

Rewriting Eqs. (B.1) and (B.2) give particle sizes $r_p$ in the Epstein and Stokes regimes for a fixed $\tau_s$. Figure B.1 shows critical particle sizes of $r_p = \frac{9}{4}l_{mfp}$ for different $\alpha_{acc}$ and particle sizes in either Epstein and Stokes regimes for a fixed $\tau_s$. Because $v_{acc} \propto \alpha_{acc}r^{\frac{3}{2}-\beta}$, $H_g \propto r^{\frac{3-\beta}{2}}$, $\Sigma_g \propto \alpha_{acc}^{-1}r^{\beta-\frac{3}{2}}$, $\rho_g \propto \alpha_{acc}^{-1}r^{\frac{3}{2}-\beta-3}$, and $v_{th} \propto r^{-\frac{\beta}{2}}$, the particle size in the Epstein regime becomes $r_p \propto \alpha_{acc}^{-1}r^{\beta-\frac{3}{2}}$, that is, $r_p \propto r^{-1}$ for $\beta = 1/2$ (the gray lines in Fig. B.1). Because $l_{mfp} \propto 1/\rho_g \propto \alpha_{acc}r^{3-\frac{3}{2}\beta}$, the particle size in the Stokes regime becomes $r_p \propto r^{\frac{3}{4}-\frac{\beta}{4}}$, which is independent on $\alpha_{acc}$, that is, $r_p \propto r^{\frac{5}{8}}$ for $\beta = 1/2$ (the black lines in Fig. B.1).

## Appendix C: Dependence on the sublimation width

Here, we discuss the dependence on the sublimation width $\Delta x_{subl}$. The detailed derivation of the sublimation width is discussed in Sect. 3. In Fig. 5, a realistic variable sublimation width is used (Eq. (74)). In contrast, Fig. C.1 is the same as Fig. 5 but the sublimation width is fixed ($\Delta x_{subl} = 0.1 H_g$).

As the sublimation width $\Delta x_{subl}$ becomes larger, the average scale height of silicate dust becomes larger because the silicate dust released in a wider radial width diffuses vertically and radially mixes (Eq. (36)), leading to a lower concentration of silicate dust in the disk midplane. This effect is more significant for small $\alpha_{Dr} (= \alpha_{Dz})/\alpha_{acc}$ (i.e., the advection-dominated regime) as both radial and vertical mixings are ineffective, while large $\alpha_{Dr} (= \alpha_{Dz})/\alpha_{acc}$ (i.e., the diffusion-dominated regime) case is rarely affected because the mixing effects “erase” the information of the initial distribution of silicate dust. Thus, the parameter space for runaway pile-up of silicate dust is reduced for small $\alpha_{Dr} (= \alpha_{Dz})/\alpha_{acc}$ in the case of a variable $\Delta x_{subl}$ (Eq. (74)) compared to the case of the fixed $\Delta x = 0.1 H_g$ (compare Fig. C.1 to Fig. 5).

<figure><img src="image_11.png" /><figcaption>**Fig. C.1.** Same as Fig. 5, but for the case of a constant $\Delta x_{subl} = 0.1 H_g$.</figcaption></figure>