© ESO 2021

1 Introduction

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

The role of ice lines, and particularly of the water snow line, has long been recognized as potentially important: the sublimation of pebbles across the snow line leads to an outward diffusion of vapor and recondensation, thus increasing the pebble surface density (e.g., Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). Recondensation of outwardly diffused water vapor onto pebbles also induces local pebble growth (Ros & Johansen 2013; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). The sublimation of icy pebbles changes the ionization rate across the snow line to form a local pressure bump, and the bump could stop the pebble drift there, resulting in their pile-up under favorable disk conditions (Kretke & Lin 2007; Brauer et al. 2008; Ida & Lin 2008). These mechanisms potentially lead to the formation of ice-rich planetesimals¹. However, in the present work, our aim is to show that the sublimation of icy pebbles can also result in the formation of silicate-rich (rocky) planetesimals inside the snow line.

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

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

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

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

2 Previous studies and an analytical solution

2.1 Setting

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

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

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

Fig. 1 Histogram of the time-dependent values of Fp/g ≡Ṁp∕Ṁg, the ratio of the pebble to gas mass flux at the inner boundary of three protoplanetary disks (models L, M and S, respectively) from Elbakyan et al. (2020). For this setting, irregularities in the radial and azimuthal structures of the disks yield a highly time-dependent Fp/g value.

2.2 A comparison of previous studies

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

Schoonenberg & Ormel (2017) performed detailed 1D diffusion-advection simulations of gas and icy pebbles with an Eulerian grid code to study pile-ups of the icy pebbles in an accretion disk, taking account of the back reaction to their drift velocity, sublimation process, and both radial and vertical turbulent diffusion. Near the snow line, re-condensation of water vapor diffused from the region inside the snow line enhances the ice surface density. For icy pebbles with Stokes numbers of 10⁻² −1, SI can occur for the solid-to-gas surface density ratio Z_Σ ≳ 0.02−0.03 (Johansen et al. 2009; Carrera et al. 2015; Yang et al. 2017). They found that the ice surface density is enhanced by a factor of several by the recycling process in the case of relatively vigorous turbulence (α ~ 10⁻³−10⁻²) and a relatively high pebble-to-gas mass accretion rate (F_p/g ~ 0.8).

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

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

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

2.3 Analytical derivation

As we described in Sect. 1, we assume that icy pebbles consist of a large number of small silicate particles covered by icy mantles. When the pebbles drift inward and pass the snow line, the icy mantles are sublimated, and consequently the small silicate particles are released. We set the Stokes number of icy pebbles to be τ_s,p = 0.1, which is relevant to the Epstein regime (Okuzumi et al. 2012; Ida et al. 2016), and that of silicate particles is τ_s,d ≪ 1, where the Stokes number is defined by the stopping time due to gas drag t_stop and Keplerian frequency Ω_K as τ_s= t_stopΩ_K. The notations we use here are summarized in Table 1.

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

Firstly, we derived the vertically averaged solid-to-gas mass ratio (metallicity), inside the snow line. In this paper, the subscripts, “g”,“p”, and “d” represent the disk’s H-He gas, icy pebbles, and silicate dustparticles, respectively. In the steady state, the surface density of the silicate dust particles Σ_d and the disk gas Σ_g are given by $$\begin{array}{ll}{\Sigma }_{\text{d}}\hfill & =f_{\text{d/p}}{\dot{M}}_{\text{p}}/2\pi r{v}_r,\hfill \\ {\Sigma }_{\text{g}}\hfill & ={\dot{M}}_{\text{g}}/2\pi r{u}_r,\hfill \end{array}$$

where v_r and u_r are the drift velocity of the silicate particles and the gas accretion velocity, and f_d/p is the silicate dust mass fraction in the icy pebbles. We adopted f_d/p = 0.5 as a nominal value. The vertically averaged metallicity is $$Z_{\Sigma }\equiv \frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}=f_{\text{d/p}}{F}_{\text{p/g}}\frac{u_r}{v_r},$$(3)

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

We also used the local solid-to-gas ratio defined by $$Z\equiv \frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}},$$(5)

where ρ_g and ρ_d are the local densities of gas and silicate dust particles.

The radial drift velocity of the particles and disk gas is given by Ida & Guillot (2016) and Schoonenberg & Ormel (2017): $$v_r=-{\Lambda }^2\frac{2{\tau }_{\text{s},\text{d}}}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\eta v_{\text{K}}-\Lambda \frac{1}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}{u}_{\nu },$$(6)

where v_K is Keplerian velocity, τ_s,d is the Stokes number of the particles, and η is the degree of deviation of the gas rotation angular velocity (Ω) from Keplerian one (Ω_K), given by $$\begin{array}{l}\eta \equiv \frac{{\Omega }_{\text{K}}-\Omega }{{\Omega }_{\text{K}}}=C_{\eta }{\left(\frac{H_{\text{g}}}{r}\right)}^2,\\ C_{\eta }\equiv \frac{1}{2}\frac{\partial \ln P}{\partial \ln r},\end{array}$$

where C_η = 1.3−1.4, depending on the disk structure (Ida et al. 2016). We used C_η = 11∕8 for Σ_g ∝ r⁻¹ and T ∝ r^−1∕2 in this paper. In the above equation, u_ν is an unperturbed disk gas accretion (advection) velocity given by $$u_{\nu }\simeq \frac{3{\nu }_{\text{acc}}}{2r}\simeq \frac{3{\alpha }_{\text{acc}}{H}_{\text{g}}^2\Omega }{2r}\simeq \frac{3{\alpha }_{\text{acc}}}{2}{\left(\frac{H_{\text{g}}}{r}\right)}^2{v}_{\text{K}},$$(9)

where ν_acc is the effective viscosity for the angular momentum transfer of the disk gas, H_g is the disk gas scale height, and u_ν is defined to be positive for inward flow. The effect of back reaction from the piled-up particles is represented by Λ, which is defined by $$\Lambda \equiv \frac{{\rho }_{\text{g}}}{{\rho }_{\text{g}}+{\rho }_{\text{d}}}=\frac{1}{1+Z}.$$(10)

Because ${\rho }_{\text{d}}\simeq {\Sigma }_{\text{d}}/\sqrt{2\pi }{H}_{\text{d}}$, where H_d is the scale hight of the silicate dust particles, and ${\rho }_{\text{g}}\simeq {\Sigma }_{\text{g}}/\sqrt{2\pi }{H}_{\text{g}}$, Eq. (10) for z < H_d is approximated for the particle motion by $${\Lambda }^{-1}=1+\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\simeq 1+Z_{\Sigma }\frac{H_{\text{g}}}{H_{\text{d}}}\simeq 1+h_{\text{d/g}}^{-1}{Z}_{\Sigma }.$$(11)

We define the scale height ratio: $$h_{\text{d/g}}\equiv H_{\text{d}}/H_{\text{g}}.$$(12)

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

Because τ_s,d ≪ 1, silicate particle motions are strongly coupled with disk gas, but the back reaction to the drift velocity of the silicate particles must be be taken into account. Equation (6) with τ_s,d ≪ 1 implies v_r ≃ Λu_ν. Substituting u_r ≃ u_ν, v_r ≃ Λu_ν, and Eq. (11) into Eq. (3), we obtain $$Z_{\Sigma }\simeq \frac{f_{\text{d/p}}}{\Lambda }{F}_{\text{p/g}}\simeq f_{\text{d/p}}(1+h_{\text{d/g}}^{-1}{Z}_{\Sigma })F_{\text{p/g}},$$(13)

which issolved as (Ida & Guillot 2016) $$Z_{\Sigma }=\frac{f_{\text{d/p}}{F}_{\text{p/g}}}{1-h_{\text{d/g}}^{-1}{f}_{\text{d/p}}{F}_{\text{p/g}}}.$$(14)

The midplane Z = ρ_d∕ρ_g is given by $$Z=\frac{Z_{\Sigma }}{h_{\text{d/g}}}=\frac{f_{\text{d/p}}{F}_{\text{p/g}}}{h_{\text{d/g}}-f_{\text{d/p}}{F}_{\text{p/g}}}.$$(15)

If we consider the back reaction to gas motion, the divergence condition is only slightly modified to (see Appendix A) $$Z=\frac{f_{\text{d/p}}{F}_{\text{p/g}}}{h_{\text{d/g}}-(1-h_{\text{p/g}})f_{\text{d/p}}{F}_{\text{p/g}}}.$$(16)

The metallicity Z diverges, which would lead to rocky planetesimal formation, for $$\begin{array}{ll}{F}_{\text{p/g}}\hfill & >\frac{h_{\text{d/g}}}{f_{\text{d/p}}}[\text{w/o the back reaction to gas}],\hfill \\ F_{\text{p/g}}\hfill & >\frac{h_{\text{d/g}}}{(1-h_{\text{p/g}})f_{\text{d/p}}}[\text{w/ the back reaction to gas}].\hfill \end{array}$$

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

Schoonenberg & Ormel (2017) neglected the back reaction to the silicate dust drift velocity and assumed Λ =1 in Eq. (13). In this case, Eq. (13) is reduced to Z_Σ ≃ f_d/pF_p/g, which showsno divergence. This demonstrates that the back reaction to the silicate particle radial velocity plays an essential role in the occurrence of the runaway pile-up.

In the runaway pile-up condition, $F_{\text{p/g}}>f_{\text{d/p}}^{-1}{h}_{\text{d/g}}$, the estimation of H_d is very important. The scale height of pebbles (H_p) is given with their Stokes number, τ_s,p, and the vertical mixing parameter α_Dz as (Dubrulle et al. 1995; Youdin & Lithwick 2007) $$h_{\text{p/g}}=\frac{H_{\text{p}}}{H_{\text{g}}}={\left(1+\frac{{\tau }_{\text{s},\text{p}}}{{\alpha }_{\text{D}z}}\right)}^{-1/2}.$$(19)

We discuss the increase of H_p by Kelvin-Helmholtz instability for the vertical shear due to pile-up of particles in Paper II.

The scale height of silicate dust particles (H_d) must be the same as H_p just after the release from the icy pebbles. Because τ_s,d ≪ τ_s,p, H_d is increased by the vertical turbulent stirring afterward. Ida & Guillot (2016) assumed that H_d ≃ H_p (h_d/g = h_p/g) until the runaway pile-up develops. For h_d/g ~ 0.1 (τ_s,peb∕α_Dz ~ 10²) and f_d/p ~ 0.5, the optimistic estimate by Eq. (15) predicts that the runaway pile-up occurs for F_p/g > 0.2, which would be available for relatively early phase of disk evolution.

On the other hand, Schoonenberg & Ormel (2017) assumed immediate vertical stirring by turbulence, that is, H_d ≃ H_g (h_d/g ≃ 1). Because they neglected the back reaction to v_r, they never found the runway pile-up in their simulations. Even if we apply Eq. (15), which was derived with the back reaction, for the case of h_d/g ≃ 1 and f_d/p ~ 0.5, the runaway pile-up condition is F_p/g > 2. The high value of F_p/g may not be completely ruled out, but it is not easily established (Ida & Guillot 2016, and see Fig. 1 in this paper).

Hyodo et al. (2019) treated the effects of the vertical and radial diffusion more carefully. They considered the vertical stirring characterized by α_Dz to derive a model for h_d/g as a function of x = r − r_snow. They also introduced the back reaction to the diffusion coefficients D_r and D_z by the pile-up of the particles, such that $$\begin{array}{ll}{D}_r\hfill & ={\alpha }_{\text{D}r}{H}_{\text{g}}^2\Omega {\Lambda }^K,\hfill \\ D_z\hfill & ={\alpha }_{\text{D}z}{H}_{\text{g}}^2\Omega {\Lambda }^K.\hfill \end{array}$$

A simple argument that $D_z,D_r\propto H_{\text{g}}^2\propto c_s^2\propto P/({\rho }_{\text{g}}+{\rho }_{\text{d}})\propto \Lambda$ (Laibe & Price 2014; Lin & Youdin 2017) suggests K = 1. Another argument is based on Kolmogorov theory (e.g., Cho et al. 2003). In the eddy cascade, the energy transfer rate, $({\rho }_{\text{g}}+{\rho }_{\text{d}})v_{\mathcal{l}}^2/t_{\text{eddy},}$ is independent of ℓ, where ℓ, v_ℓ, and t_eddy are the eddy size, velocity, and turnover time ~ℓ∕v_ℓ. Because D_z, D_r ~ ℓ₀v_ℓ0 where ℓ₀ and v_ℓ0 are the eddy size and velocity controlling the turbulence, $D_z,D_r\propto {\mathcal{l}}_0{[{\mathcal{l}}_0/({\rho }_{\text{g}}+{\rho }_{\text{d}})]}^{1/3}\propto {\Lambda }^{1/3}$. This argument thus implies K = 1∕3. This should be a lower limit because an increase in the solid-to-gas is expected to suppress turbulence itself. Given these uncertainties, Hyodo et al. (2019) performed runs with K = 0, 1, and 2. They found that the runaway pile-up is caused by inclusion of the back reaction to v_r, and that the condition for the runaway pile-up depends on the relations among α_acc, α_{D r}, and α_{D z}, while Ida & Guillot (2016) considered the limits of α_Dr, α_Dz → 0 and their condition, Eq. (15), does not include α_acc. The condition derived by Hyodo et al. (2019) is more severe than the prediction by Ida & Guillot (2016). The runaway pile-up occurs for F_p/g > 0.2−0.3, and α_{D r} = α_Dz ≃ 10⁻³ in the case of K = 1 and 2 with α_acc = 10⁻², while it does not occur with K = 0 as long as it is in the parameter space of α_Dr = α_Dz = 10⁻³−10⁻² and F_p/g < 0.6. In this paper, we show this back-reaction effect on the diffusion coefficients more clearly.

2.4 Stokes number of silicate particles

The models discussed thus far and the model that we present are based on the assumption that many silicate dust particles released by the individual icy pebbles remain small compared to pebbles (i.e., with a Stokes number that is much smaller than unity). In this case, the radial drift is much faster for icy pebbles than for silicate particles. Because the gas accretion velocity is given by $u_{\nu }~(3/2){\alpha }_{\text{acc}}{(H_{\text{g}}/r)}^2{v}_{\text{K}}$ (Eq. (9)) and the drift speed relative to the gas is given by $v_r~{\tau }_{\text{s},\text{d}}{(H_{\text{g}}/r)}^2{v}_{\text{K}}$ (Eq. (6)), the total drift velocity is dominated by u_ν if τ_s,d ≲ α_acc. Therefore, the assumption here is expressed as τ_s,d ≲ α_acc < τ_s,p.

In Sect. 2.3, we assume that τ_s,p ~ 0.1, α_acc ~ 10⁻³−10⁻², and τ_s,d ≲ α_acc. Saito & Sirono (2011) assumed that many μm-sized silicate particles are embedded in pebbles. Morbidelli et al. (2015) considered chondrule-sized (~mm) particles. Because the Stokes number is proportional to the particle size in the Epstein regime and its square in the Stokes regime, it is reasonable to assume that τ_s,d < α_acc. This model corresponds to the “many-seeds model” of Schoonenberg & Ormel (2017). In the Monte Carlo simulations here, we adopt τ_s,d = 10⁻⁵.

In the silicate particle pile-up region, the particles may quickly grow up to sizes determined by the threshold collision velocity for fragmentation/rebound. The collision velocity is set by the velocity dispersion induced by the turbulence, $~{(3{({\alpha }_{\text{D}r}^2+{\alpha }_{\text{D}z}^2)}^{1/2}{\tau }_{\text{s},\text{d}})}^{1/2}{c}_s,$ or by the drift velocity difference between the particles ~ 2τ_s,d η v_K (e.g., Sato et al. 2016). The collision velocity between silicate particles is therefore ${({\tau }_{\text{s},\text{d}}/{\tau }_{\text{s},\text{p}})}^{-(0.5-1)}$ times that between icy pebbles. If the threshold collision velocity is ten times lower for the silicate particles than for the icy particles (Blum & Wurm 2000; Zsom et al. 2010, 2011; Wada et al. 2011, 2013; Weidling et al. 2012) and τ_p,d ~ 0.1 is close to the fragmentation/rebound barrier for ice, the silicate particles can grow only up to τ_s,d ~ (10⁻¹–10⁻²)τ_s,p ~ 10⁻³–10⁻² in size. Therefore, for α_acc = 10⁻³–10⁻², the drift velocity of the silicate particles is comparable to u_ν and the pile-up would not be significantly affected by the value of τ_s,d.

We note that this conventional view has been challenged recently (Kimura et al. 2015; Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019). A higher fragmentation threshold for silicates would open the possibility for dust to coagulate and grow. We hence performed runs with larger values of τ_s,d to find that the results are not affected for τ_s,d ≲ α_acc and that the runaway pile-up region is removed from lower F_p/g regions as τ_s,d becomes larger.

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

3 The Monte Carlo approach

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

We considered a radially local region near the snow line and neglected the r-dependence of the Keplerian frequency (Ω) and the disk’s gas scale height (H_g). As discussed in Sect. 2.4, we adopted the conventional view to set τ_s,p = 0.1 and τ_s,d = 10⁻⁵ in our simulations.

3.1 Injection of particles

To model the release from drifting icy pebbles due to sublimation, at each time step δt, we randomly injected a new silicate super-particle with mass m with a radially uniform distribution near the snow line in the − 0.5Δx_subl < x < 0.5Δx_subl range, and a Gaussian distribution of z of the root mean square Δz_subl as $$\begin{array}{ll}x\hfill & =0.5R\Delta x_{\text{subl}},\hfill \\ z\hfill & =\sqrt{2}{\text{erf}}^{-1}(|R|)\frac{R}{\left|R\right|}\Delta z_{\text{subl}},\hfill \end{array}$$

where $R$ is a random number in the range of [−1, 1], Δx_subl is the characteristic sublimation radial length, and Δz_subl corresponds to the scale height of theincoming icy pebbles, H_p, given by Eq. (19). We adopted Δx_subl = 0.1H_g as a nominal case. In Sect. 4.4.2, we also consider a more complicated function that fits the numerical result of the grid code simulation obtained by Paper II.

Here, we distinguish the effective viscosity parameters for advection, radial mixing, and vertical mixing, denoted by α_acc, α_{D r}, and α_{D z}, respectively. According to different values of α_Dz, we used the consistent value of the scale height of the injected silicate particles given by Eq. (19). Although Hasegawa et al. (2017) suggested α_Dr, α_Dz ~ 0.1 α_acc, the relations among α_Dr, α_Dz and α_acc are not clear (e.g., Armitage et al. 2013). Therefore, we surveyed broad parameter ranges of α_{D r}, α_Dz and α_acc; we mostly show the results with α_acc ≥ α_Dr = α_Dz. In some cases, the results with α_Dr≠α_Dz (α_acc ≥ α_Dr, α_Dz) are also shown.

3.2 Advection and diffusion

At each time step δt, we changed the r–z locations of the particles that were injected before, following Ciesla (2010, 2011), by $$\begin{array}{ll}\delta x\hfill & =v_{\text{adv},r}\delta t+R{(6D_r\delta t)}^{1/2},\hfill \\ \delta z\hfill & =v_{\text{adv},z}\delta t+R{(6D_z\delta t)}^{1/2},\hfill \end{array}$$

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

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

The radial drift (advection) velocity v_adv,r is the same as v_r in Eq. (6). The vertical advection velocity is (Ciesla 2010) $$v_{\text{adv},z}=-({\alpha }_{\text{D}z}+{\tau }_{\text{s},\text{d}})\Omega z.$$(26)

The first term (∝ α_Dz) comes from the effect of diffusion smoothing the concentration (ρ_d∕ρ_g) but not ρ_d itself.

3.3 Inclusion of back reaction effects

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

We utilized the “super-particle” method, where one super-particle represents the mass of a large number of particles while it suffers the same specific drag force as the individual particles. We injected one super-particle at every time step δt = Ω⁻¹. For given Ṁ_p, we can determine the individual mass of the super-particles as $$m=f_{\text{d/p}}{\dot{M}}_{\text{p}}\delta t=f_{\text{d/p}}{\dot{M}}_{\text{p}}{\Omega }^{-1},$$(27)

where f_d/p is a silicate fraction in migratingpebbles (f_d/p = 0.5 in the nominal case). To locally average the silicate dust mass density, we used linear grids, Δx and Δz. In the nominal cases, we adopted Δx = 0.1 H_g. We used Δz ~ Δz_subl near the midplane and larger Δz for upper regions to keep statistically enough number of particles in the individual grids.

The silicate particle mass density at a grid of [x, x + Δx] and [|z|, |z| + Δz] is $${\rho }_{\text{d}}=\frac{m\Delta N_{x,z}}{2\pi r\Delta x\times 2\Delta z}=\frac{H_{\text{g}}}{r}\frac{\Delta N_{x,z}{f}_{\text{d/p}}{\dot{M}}_{\text{p}}}{4\pi (\Delta x/H_{\text{g}})(\Delta z/H_{\text{g}})H_{\text{g}}^3\Omega },$$(28)

where ΔN_x,z is a total number of the particles in the grid. We simply assumed the vertical isothermal hydrodynamical equilibrium. Then, the gas density is given by $$\begin{array}{ll}{\rho }_{\text{g}}\hfill & =\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{H}_{\text{g}}}\exp \left(-\frac{z^2}{2H_{\text{g}}^2}\right)\hfill \\ \hfill & =\frac{{\dot{M}}_{\text{g}}}{\sqrt{2\pi }3\pi {\alpha }_{\text{acc}}{H}_{\text{g}}^3\Omega }\exp \left(-\frac{z^2}{2H_{\text{g}}^2}\right).\hfill \end{array}$$(29)

The silicate dust-to-gas ratio at (x, z) is $$\begin{array}{ll}Z\hfill & =\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}=\frac{3\sqrt{2\pi }}{4}\frac{H_{\text{g}}}{r}{\alpha }_{\text{acc}}\frac{f_{\text{d/p}}{\dot{M}}_{\text{p}}}{{\dot{M}}_{\text{g}}}\frac{\Delta N_{x,z}}{(\Delta x/H_{\text{g}})(\Delta z/H_{\text{g}})}\exp \left(\frac{z^2}{2H_{\text{g}}^2}\right)\hfill \\ \hfill & =3.75\times {10}^{-2}\left(\frac{H_{\text{g}}/r}{0.04}\right)\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-2}}\right)\left(\frac{f_{\text{d/p}}}{0.5}\right)F_{\text{p/g}}\hfill \\ \hfill & \times {\left(\frac{\Delta x}{0.1H_{\text{g}}}\right)}^{-1}{\left(\frac{\Delta z}{0.1H_{\text{g}}}\right)}^{-1}\Delta N_{x,z}\exp \left(\frac{z^2}{2H_{\text{g}}^2}\right).\hfill \end{array}$$(30)

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

In the simulations, we normalized length by H_g, including the range of particle injection, Δx_subl and Δz_subl, and the grid sizes, Δx and Δz. The time was scaled by Ω⁻¹, that is, each super-particle was injected at every unit’s scaled time. Therefore, our results can be applied for any location of r_snow, as long as H_g∕r is the same.

In the simulations here, we adopted f_d/p = 0.5 and H_g∕r = 0.04 at r = r_snow, as nominal values. We do not need to specify the values of H_g and Ṁ_g or the value of r_snow, because only the scaled values, h_d/g(= H_d∕H_g), h_p/g(= H_p∕H_g), H_g∕r and F_p/g(= Ṁ_p∕Ṁ_g), are used in the simulations. For the viscous α-parameters, we show that only the ratios α_Dz∕α_acc and α_Dr∕α_acc are important.

4 Results

4.1 Codecheck

We tested our code by comparing the simulated particle scale heights with the existing analytical argument, as Ciesla (2010) already did, and by reproducing the simple analytical result by Ida & Guillot (2016) and the numerical result without the back reaction to silicate particles by Schoonenberg & Ormel (2017). We also offer a detailed comparison with the results obtained by an updated code of Hyodo et al. (2019) in Paper II. While Hyodo et al. (2019) included both nonzero α_{D r} and α_{D z}, they assumed the evolution of H_d analytically estimated only from the vertical stirring with α_Dz. In our simulation, the evolution of H_d was self-consistently calculated, and we find that the radial mixing with α_{D r} is also important for H_d. On the other hand, the sublimation width Δx_subl is self-consistently calculated in Hyodo et al. (2019), while it needs to be assumed in the simulation here.

4.1.1 Confirmation of the theoretically predicted particle scale height

Figure 2 shows the distribution of particles and their scale heights (the local root mean square of z of the particles) in a steady state at t = 10⁴Ω⁻¹ (after 10⁴ particles are injected), obtained by our Monte Carlo simulations for (a) α_{D z} = 10⁻⁵, (b) 10⁻³, and (c) 10⁻¹. For these runs, we neglected the back reactions to v_r (Λ = 1) and α_Dz (K = 0). To highlight the effect of vertical stirring, we set α_Dr = 0. The other parameters are the particles’ Stokes numbers τ_s = 10⁻³ and α_acc = 10⁻². The equilibrium scale height predicted from Eq. (19) with τ_s,p replaced by τ_s = 10⁻³ are (a) 0.1H_g, (b) 0.71H_g, and (c) 1.0H_g, respectively,which are shown via the green dashed lines in the plots. As the particles drift inward, the local root mean square of z (representedby the blue solid curves) asymptotically approaches the theoretical values in each panel. Thus, we confirm that the analytically derived Eq. (19) is reproduced by our Monte Carlo simulations.

Fig. 2 Snapshots of particle distribution (red dots) on the x-z plane at time t = 104Ω−1, and the root mean square of z (blue solid curve) as a function of the radial direction x (centeredon the snow line location and in units of the gas pressure scale height Hg). We plot ∣z∣, while z takes either positive or negative values. We set αacc = 10−2, αD r = 0 and τs = 10−3. The vertical mixing parameter in each panel is (a) αDz = 10−5, (b) αD z = 10−3, and (c) αDz = 10−1. The analytical estimates of the equilibrium particle scale height for each αDz are shown by the dashed green lines. The back reactions are not included for vr (Λ = 1) and αD z (K = 0). We injected the particles in the range of x = [−0.05, 0.05]Hg and z = [−0.1, 0.1]Hg.

4.1.2 Reproduction of the results of Schoonenberg & Ormel (2017) and Ida & Guillot (2016)

Ida & Guillot (2016) adopted the following setup: (1) the turbulent diffusion is neglected for simplicity, and (2) the silicate particle pile-up is fast enough that H_d is not increased from Δz_subl(= H_p) by the vertical stirring. To mimic the setup of Ida & Guillot (2016), we set the following in the Monte Carlo simulations: (1) we adopted α_Dr = 0, K = 0, and Δz_subl = 0.03H_g; (2) to keep H_d = Δz_subl = 0.03H_g as an equilibrium, we adopted α_Dz = 10⁻⁵ and τ_s,d = 10⁻² (Eq. (19)); (3) the dust particle drift speed is set to be v_r = −Λu_ν, where the artificially enlarged τ_s,d is not reflected (Eq. (6)), and (4) for gas, u_r = −u_ν is always assumed. For comparison, v_r = −u_ν (without the back reaction to v_r) is also examined.

In the rest of the paper, we always set τ_s,p = 10⁻¹ and τ_s,d = 10⁻⁵, and we consistently calculated H_p(= Δz_subl) by Eq. (19) and the evolution of H_d from Δz_subl by the vertical and radial diffusion, with given α_Dz and α_Dr. The above artificial setting is only used in the particular runs in Fig. 3.

Figure 3 shows the time evolution of the maximum value of ρ_d ∕ρ_g (the value on the midplane near the snow line) for different values of F_p/g, in the case of α_acc = 10⁻², Δx_subl = 0.1H_g, and Δz_subl = 0.03H_g. In panel a, the back reaction to v_r is not included (Λ = 1) in the Monte Carlo calculations,that is, v_r = u_r = −u_ν, which corresponds to the result of Schoonenberg & Ormel (2017), except that they assumed h_d/g = 1. The maximum ρ_d∕ρ_g quickly reaches the equilibrium values that are given by Eq. (3), with v_r = u_r = −u_ν and Z_Σ = Z h_d/g as $$\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\simeq f_{\text{d/p}}{h}_{\text{d/g}}^{-1}{F}_{\text{p/g}}\simeq 16.7F_{\text{p/g}}.$$(31)

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

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

A comparison between Figs. 3a and b clearly shows that the back reaction to v_r of the silicate particles plays an essential role in the occurrence of the runaway pile-up. The increase in the pile-up slows down the drift velocity, and, accordingly, increases the pile-up itself. Because it is certain that the back reaction to v_r exists, we include it in the results in the rest of the paper.

If the supply rate of silicate particles (f_d/pṀ_p) exceeds a threshold value, the mode of the pile-up becomes runaway. The threshold value must be regulated by the local value of ρ_d ∕ρ_g, because the back reaction becomes effective once ρ_d∕ρ_g exceeds unity. While Ṁ_p, more exactly F_p/g, determines Z_Σ = Σ_d∕Σ_g, the local value of Z = ρ_d∕ρ_g is regulated byH_d and Z_Σ (Z = Z_Σ∕h_d/g) near the snow line. The silicate dust particlescale height H_d is regulated bythe pebble scale height H_p (= Δz_subl in our simulation) and turbulent mixing parmeters, α_Dz and α_Dr, as shown in Sect. 4.3.

Fig. 3 Time evolution of the maximum ρd∕ρg. The back reaction to vr is included in the results in panel b, but not in panel a. The light blue, magenta, blue, green, and red curves are the results with Fp/g = 0.025, 0.05, 0.1, 0.2, 0.4, and 0.8, respectively. The other parameters are fixed as αacc = 10−2 and Δxsubl = 0.1Hg and Δzsubl = 0.03Hg. To mimic the settings of Ida & Guillot (2016), we artificiallly set αDr = 0, αD z = 10−5 and τs,d = 10−2 (however, the relatively large τs,d for the silicate dust is not reflected to the dust drift speed). The latter two parameters maintain Hd = Δzsubl. The circles in panel a represent the analytical solution given by Eq. (31). That in panel b is the analytical solution for non-divergent case (Fp/g = 0.025) given by Eq. (15).

Fig. 4 Snapshots of silicate particles at t = 104 Ω−1 (red dots), and the root mean square of z (blue solid curves) at each grid of x. The dashed blue curves are given by Eq. (35) with Δxsubl replaced by max(Δxsubl, |x|) (see thediscussion at the end of Sect. 4.3.2). The magenta curve represents log10 (ρd∕ρg) near the midplane. We set αacc = 10−2, τs,d = 10−5, and fd/p = 0.5, and Ṁg = 10−8M⊙ yr−1. The nine panels adopt different values of Fp/g and αD z = αDr, as indicated in each panel.

4.2 Typical results

We first show typical results in our simulations for different F_p/g and α_{D z} (= α_Dr). As we show below, the radial diffusion also contributes to H_d. The back reaction to v_r is included. In these results, Δx_subl = 0.1H_g is adopted, as well as α_acc = 10⁻².

Figure 4 shows the snapshots of silicate particles at 10⁴ Ω⁻¹ for the cases without the back reaction to diffusion coefficients (K = 0). In the right panels with α_Dz = α_Dr = 10⁻² (α_{D z}∕τ_s,d = 10³), the vertical mixing is fast enough to realize the upper limit of H_d (H_d ≃ H_g) for all of F_p/g = 0.1, 0.3, and 1. Even for a high-mass flux of pebbles (F_p/g = 1.0), no pile-up of the silicate particles is found for α_Dz = α_Dr = 10⁻². In the case of α_Dz = α_Dr = 10⁻⁶, on the other hand, the mixing is so weak that H_d is not increased to the equilibrium value [$\simeq {({\alpha }_{\text{D}r}/{\tau }_{\text{s},\text{d}})}^{1/2}{H}_{\text{g}}~0.3H_{\text{g}}$] in the range of x > −2 H_g. Nevertheless, H_d near the injection point (x ~ 0) is much higherthan Δz_subl = H_p = 0.003H_g. This means that H_d is regulated by α_Dr as well as byα_Dz, which is discussed in Sect. 4.3.

Figure 5 show the results with the back reaction to diffusion coefficients (K = 1). In general, the runaway pile-up is more pronounced for smaller α_Dz = α_Dr and larger F_p/g. The magenta curves represent log₁₀(ρ_d∕ρ_g) near the midplane. The midplane ρ_d∕ρ_g usually takes the maximum value, at the injection (sublimation) region, x ~ 0, which is denoted by ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$. In the panels with ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0\lesssim 1$ (the panels with α_Dz = 10⁻² and those with α_Dz = 10⁻⁴ and F_p/g = 0.3 and 0.1), the results are similar to those in Fig. 4. If ρ_d∕ρ_g ≳several in Fig. 4, however, the pile-up is much more pronounced in this case, because the radial and vertical diffusions become much weaker as the pile-up proceeds.

Figure 6 shows the time evolution of ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ and H_d∕H_p at the maximum ρ_d∕ρ_g locations (h_d/g,0), for different F_p/g. In panel a, the back reaction to diffusion is not included (K = 0). As Figs. 4 and 5 suggest, the transient α_{D z} (= α_Dr) values between the runaway pile-up and non-pile-up cases are α_Dz ≃ 10⁻⁵ for K = 0 and α_{D z} ≃ 10⁻⁴ for K = 1, so we plot the results with α_Dz = 10⁻⁵ in panel a and α_Dz = 10⁻⁴ in panel b. These figures show that transition to the runaway pile-up is much clearer in the case of K = 1.

This difference comes from the H_d evolution. The comparison with the result of Ida & Guillot (2016) in Sect. 2 suggests that H_d is directly related to a threshold value of F_p/g for the runaway pile-up: a lower H_d leads to a lower threshold value of F_p/g. However, Fig. 6 shows that H_d increases as ρ_d∕ρ_g increases in the case of K = 0, while it decreases for K = 1. Therefore, the pile-up suffers negative feedback for K = 0, but it suffers positive feedback for K = 1, resulting in a much clearer transition to the runaway pile-up.

In both cases, the back reaction to v_r is included. As the pile-up proceeds, v_r decreases, which enhances the pile-up. However, we also included the vertical stirring. As v_r decreases, the silicate particles are stirred up for a longer period of time before they leave the sublimation region, resulting in higher H_d near that location. In panel a, H_d near the sublimation region quickly asymptotes to an equilibrium value for F_p/g ≤ 0.2 (lower ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ cases), and it increases with time for F_p/g = 0.8 (the runaway pile-up cases). In panel b, in the runaway pile-up cases (F_p/g = 0.4 and 0.8), H_d initially increases with time. However, after ρ_d∕ρ_g exceeds ~1, the diffusion coefficients become so small that the settling of the particles overwhelms the vertical stirring, and H_d decreases. Thus, the different responses to the pile-up between K = 0 and K = 1 result in different H_d and ρ_d∕ρ_g evolution in the pile-up cases. As we discussed in Sect. 3.2, K = 1 (with the diffusion back reaction) is more realistic.

Fig. 5 Same as Fig. 4, except for K = 1.

Fig. 6 Time evolution of ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ and the dust scale height at the maximum ρd∕ρg locations scaled by the initial one (= Hp). The back reaction to diffusion is included with K = 1 in panel b, but not (K = 0) in panel a, while the back reaction to vr is included in both panels a and b. The light blue, magenta, blue, green, and red curves are the results with Fp/g = 0.05, 0.1, 0.2, 0.4 and 0.8, respectively. The other parameters are fixed as αacc = 10−2, Δxsubl= 0.1, and Δzsubl = 0.03. The diffusion parameters are αDz = αDr = 10−5 in panel a and αDz = αDr = 10−4 in panel b.

4.3 Silicate particle scale height

4.3.1 Monte Carlo simulation results

Figure 7 shows ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ and hd/g, 0 at t = 15 000 Ω⁻¹ as a function of α_Dz = α_Dr for different F_p/g. This figure more clearly shows that the asymptotic equilibrium values of hd/g, 0 are independentof F_p/g in the non-runaway pile-up cases, which is suggested by Fig. 6. The envelope curves in the plots for (a) K = 0 and (b) K = 1 match the analytical formula derived in the next subsubsection (Eqs. (3) and (35) with Eq. (36)). The equilibrium values are generally larger than the scale height of the injected silicate particles, Δz_subl, due to the effects of vertical and radial diffusion. In panels a and b, the values of hd/g, 0 that deviate from the envelope curves correspond to the runaway pile-up; hd/g, 0 deviates to higher values for K = 0 and lower values for K = 1. If we calculate on longer timescales, the deviations increase.

4.3.2 Analytical formula for the silicate scale height

The red dashed curves in Fig. 7 are analytical formulas for the non pile-up cases, which are derived as follows. The timescale to drift by Δx_subl is $$\begin{array}{ll}{t}_{\text{drift}}\hfill & \simeq \frac{\Delta x_{\text{subl}}}{v_r}\hfill \\ \hfill & \simeq \frac{\Delta x_{\text{subl}}}{(3/2)\Lambda {\alpha }_{\text{acc}}{(H_{\text{g}}/r)}^2{v}_{\text{K}}}\simeq \frac{\Delta x_{\text{subl}}/H_{\text{g}}}{(3/2)\Lambda {\alpha }_{\text{acc}}(H_{\text{g}}/r)\Omega }.\hfill \end{array}$$(32)

As the pile-up proceeds (Λ = 1∕(1 + Z) decreases), the particle drift becomes slower, and, accordingly, t_drift becomes longer. The vertical diffusion length during t_drift is $$\begin{array}{ll}{z}_{\text{diff}}\simeq \sqrt{D_z{t}_{\text{drift}}}\hfill & \simeq {\left(\frac{2}{3}{\alpha }_{D,z}{\Lambda }^K{H}_{\text{g}}^2\Omega \times \frac{\Delta x_{\text{subl}}/H_{\text{g}}}{\Lambda {\alpha }_{\text{acc}}(H_{\text{g}}/r)\Omega }\right)}^{1/2}\hfill \\ \hfill & \simeq {\left(\frac{2}{3}{\Lambda }^{K-1}\frac{\Delta x_{\text{subl}}/H_{\text{g}}}{H_{\text{g}}/r}\frac{{\alpha }_{D,z}}{{\alpha }_{\text{acc}}}\right)}^{1/2}{H}_{\text{g}}.\hfill \end{array}$$(33)

This expression suggests that z_diff becomes higher as the pile-up proceeds (Λ decreases) for the case of K = 0. Hereafter, we derive the equilibrium value of H_d∕H_g for the non-pile-up cases with Λ ≃ 1 (the envelope curves in Fig. 7). If Z ≫ 1 for the equilibrium value of H_d∕H_g, the runaway pile-up occurs, so Z keeps increasing and the equilibrium value of H_d∕H_g no longer exists.

If only vertical diffusion is considered, the silicate particle scale height scaled by H_g is predicted as $$\begin{array}{ll}{h}_{\text{d/g}}=\frac{H_{\text{d}}}{H_{\text{g}}}\hfill & \simeq \frac{\max \left(\Delta z_{\text{subl}},z_{\text{diff}}\right)}{H_{\text{g}}}\simeq \frac{{(\Delta z_{\text{subl}}^2+z_{\text{diff}}^2)}^{1/2}}{H_{\text{g}}}\hfill \\ \hfill & \simeq {\left(\frac{\Delta z_{\text{subl}}^2}{H_{\text{g}}^2}+\frac{2}{3}\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}\frac{\Delta x_{\text{subl}}/H_{\text{g}}}{H_{\text{g}}/r}\right)}^{1/2}.\hfill \end{array}$$(34)

Furthermore, during t_drift, the particles are radially mixed in the range of $x_{\text{diff}}\simeq \sqrt{D_r{t}_{\text{drift}}}~{({\alpha }_{\text{D}r}/{\alpha }_{\text{D}z})}^{1/2}{z}_{\text{diff}}$. In general, the dust scale height increases with distance from the snow line. If x_diff ≫ Δx_subl, the radial diffusion brings the high-z particles back to the injection region and raises the value of h_d/g there. To take this effect into account, we found through comparison with the Monte Carlo simulation results that a good-fit analytical formula is obtained by multiplying the above h_d/g by a factor of $[1+{(x_{\text{diff}}/\Delta x_{\text{subl}})}^2]$. When α_{D r} is comparable to α_acc, the silicate scale height is radially smoothed out, and the effect of radial mixing for the scale height becomes weak. To take this into account, we replaced α_Dr∕α_acc in ${(x_{\text{diff}}/\Delta x_{\text{subl}})}^2$ by $({\alpha }_{\text{D}r}/{\alpha }_{\text{acc}})/[1+{(C_{r,\text{diff}}{\alpha }_{\text{D}r}/{\alpha }_{\text{acc}})}^2]$ in the above equation. Comparing the formula with the numerical results (Figs. 7 and 8), we empirically adopted C_r,diff = 10. With these modifications, $$\begin{array}{ll}{h}_{\text{d/g},\text{*}}\simeq \hfill & {\left(\frac{\Delta z_{\text{subl}}^2}{H_{\text{g}}^2}+\frac{2}{3}\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}\frac{\Delta x_{\text{subl}}/H_{\text{g}}}{H_{\text{g}}/r}\right)}^{1/2}\hfill \\ \hfill & \times \left[1+\frac{2}{3}\frac{({\alpha }_{\text{D}r}/{\alpha }_{\text{acc}})}{1+{(C_{r,\text{diff}}{\alpha }_{\text{D}r}/{\alpha }_{\text{acc}})}^2}\frac{1}{(H_{\text{g}}/r)(\Delta x_{\text{subl}}/H_{\text{g}})}\right].\hfill \end{array}$$(35)

Because h_d/g cannot exceed the equilibrium value $h_{\text{d/g},\text{eq}}\simeq {(1+{\tau }_{\text{s},\text{d}}/{\alpha }_{\text{D}z})}^{-1/2}$, the final form is $$h_{\text{d/g}}\simeq {(h_{\text{d/g},\text{eq}}^{-1}+h_{\text{d/g},\text{*}}^{-1})}^{-1}.$$(36)

As the lower panels in Fig. 7 show, the analytical formula reproduces the envelope curves of H_d ∕H_p obtained by the numerical simulations for the non-pile-up cases. We note that hd/g, 0 and ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ depend on the ratios of α_Dr∕α_acc and α_Dz∕α_acc, but not on their absolute values.

Equation (35) is the analytical formula for h_d/g at x ~ 0. This formula can be extrapolated to any x by replacing Δx_subl with max(Δx_subl, |x|, ϵ), where ϵ is added to avoid the divergence in the case of extremely small Δx_subl; we adopted ϵ = 0.01 H_g. In Figs. 4 and 5, the scale height distribution obtained by the Monte Carlo simulations are well reproduced by the analytical formula (the blue dashed curves), except for the runaway pile-up cases, where the back reaction modulates the scale height as we discussed in this subsection.

Fig. 7 Maximum ρd∕ρg and Hd ∕Hg at the maximum ρd∕ρg locations at t = 15000 Ω−1 as a functionof αDz = αDr for different Fp/g. The back reaction to diffusion is included with K = 1 in panel b, but not (K = 0) in panel a. The other parameters are fixed as αacc = 10−2, Δxsubl= 0.1Hg, and ΔzsublHp given by Eq. (19). Upper and bottom panels: the red dashed curves are the analytical formulas given by Eqs. (3) and (35) with Eq. (36). Lower panels: the green dashed lines represent Δzsubl ∕Hg.

4.4 Silicate particle pile-up

4.4.1 Analytical formula for the silicate pile-up

We can also derive an analytical formula for ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ (the midplane particle-to-gas ratio at x ~ 0) for the non-runaway pile-up cases and for the phase before the runaway pile-up proceeds. If the estimated ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ exceeds ~ 1, the actual value of ρ_d∕ρ_g should deviate from the estimate and increase with time.

Basically, ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ is predicted by Eq. (15) with the scale hight formula, Eqs. (36) and (35). A possible detailed correction is to include the effect of outward diffusion flux (F_d,D) of silicate particles beyond the snow line that can be found in the snapshots with α_{D r} = α_acc = 10⁻² in Figs. 4 and 5. As shown in Appendix B, calibrating with the results of Monte Carlo simulations in Figs. 7 and 8, the effect can be achieved with f_d/pF_p/g being replaced by F_d,net = f_d/pF_p/g − F_d,D, where $$F_{\text{d},\text{net}}=f_{\text{d/p}}{F}_{\text{p/g}}\left[1-{\left(1+\frac{15}{2}\frac{H_{\text{g}}}{r}\frac{{\alpha }_{\text{acc}}}{{\alpha }_{\text{D}r}}\right)}^{-1}\right].$$(37)

Accordingly, Eq. (15) is reduced to $${\left(\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\right)}_0=\frac{F_{\text{d},\text{net}}}{h_{\text{d/g}}-F_{\text{d},\text{net}}}.$$(38)

This formula reproduces the numerical results in Figs. 7 and 8, as long as ρ_d ∕ρ_g ≲ 1. We note, however, that if the effect of silicate particles that diffuse beyond the snow line sticking to icy pebbles is considered, most of the silicate grains beyond the snow line may eventually come back with the pebbles, so F_d,net would become more similarto f_d/pF_p/g.

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

The formulas imply that the results are scaled by α_Dz∕α_acc and α_Dr∕α_acc for given Δz_subl, Δx_subl, and H_g ∕r. We also carried out the Monte Carlosimulations with α_acc = 10⁻³ and confirmed that the plots in Fig. 7 do not change as long as we use α_{D r} ∕α_acc (= α_Dz∕α_acc) as the horizontal coordinate.

The sublimation width, Δx_subl, is an independent parameter. The results with Δx_subl = H_g (Fig. 8b), are compared to those with Δx_subl = 0.1H_g (Fig. 8a, which are identical to Fig. 7b). In Fig. 8b, the injected particles are distributed in region ten-times broader than in the case of Fig. 8a. As a result, the vertical mixing is more effective (t_drift is longer) and the pile-up (the maximum ρ_d∕ρ_g) is lower.

So far, the results with α_Dr = α_Dz have been shown. In Fig. 8c, α_Dr is fixed and α_Dz is independently changed. In Fig. 8d, α_Dz is fixed. The envelope curves for H_d∕H_g (= h_d/g) are also reproduced by the analytical formula in these cases, which strongly suggests that the correction factors with the vertical and radial diffusion are physically correct.

Fig. 8 Same as Fig. 7, except for αDz∕αDr relations and the value of Δxsubl. The maximum ρd∕ρg and hd/g (= Hd ∕Hg) for different Fp/g are given as functions of (a) αD r = αDz with Δxsubl = 0.1 Hg, which is the same as Fig. 7b, for comparison, (b) αD r = αDz with Δxsubl = Hg, (c) αD z with αDr = 10−4 and Δxsubl = 0.1 Hg, and (d) αD r with αDz = 10−4 and Δxsubl = 0.1 Hg. For all the cases, K = 1 is adopted. Upper and bottom panels: the red dashed curves are the analytical formulas given by Eqs. (35) and (3) with corresponding parameter values.

4.4.2 The condition for the silicate runaway pile-up

The analytical formulas reproduce the numerical results for a broad range of multidimensional parameters. Figure 9a shows the contour maps of the maximum ρ_d∕ρ_g on the planeof α_Dz∕α_acc(= α_Dr∕α_acc), and F_p/g at t = 1.5 × 10⁴ Ω⁻¹ for Δx_subl = 0.1 H_g. In Fig. 9b, Δx_subl is given by Eq. (41). In both cases, α_acc = 10⁻² and K =1. The left and right panels show the results of Monte Carlo simulations and the analytical results, respectively. The red regions represent the parameter ranges of runaway pile-up of silicate particles. As expected, the analytical formulas reproduce the numerical results.

For these simulations, α_acc = 10⁻² and α_Dz(= α_Dr) are in the range of [10⁻⁶, 10⁻²]. The analytical formula to predict ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$ is Eq. (37) with Eqs. (38), (35), and (36). They are functions of α_{D r} ∕α_acc and α_{D z} ∕α_acc for given Δx_subl, as long as the term $\Delta z_{\text{subl}}^2$ in Eq. (35) is negligible, which is satisfied in the parameter range we considered, as shown inFigs. 7 and 8. We also carried out simulations with α_acc = 10⁻³ and α_{D z} (= α_Dr) in the range of [10⁻⁷, 10⁻³] to find thatthe contour map is almost identical, confirming that the results are scaled by α_{D z} ∕α_acc(= α_Dr∕α_acc).

Using the formulas, we derive the boundary of the runaway pile-up on the plane of α_{D z} ∕α_acc (= α_Dr∕α_acc) and F_p/g. Equation (38) shows that the boundary is given by F_d,net ≃ h_d/g. As shown in Figs. 9a and b, the boundary is located at α_{D z} (= α_Dr) < 3 × 10⁻⁴ for F_p/g < 1. In this parameter range, Fig. 7 shows h_d/g ≃ h_d/g,*. Since α_Dr∕C_r,diff α_acc ≪ 1 for α_acc ~ 10⁻², F_d,net ≃ f_d/pF_p/g. Therefore, the boundary is approximately given by $$\begin{array}{ll}{F}_{\text{p/g}}\hfill & \simeq f_{\text{d/p}}^{-1}{h}_{\text{d/g},\text{*}}\hfill \\ \hfill & \simeq f_{\text{d/p}}^{-1}{\left(\frac{2}{3}\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}\frac{\Delta x_{\text{subl}}/H_{\text{g}}}{H_{\text{g}}/r}\right)}^{1/2}\left[1+\frac{2}{3}\frac{{\alpha }_{\text{D}r}}{{\alpha }_{\text{acc}}}\frac{1}{(H_{\text{g}}/r)(\Delta x_{\text{subl}}/H_{\text{g}})}\right].\hfill \end{array}$$(39)

Adopting H_g∕r = 0.04 and f_d/p = 1∕2 as nominal values, $$\begin{array}{lll}{F}_{\text{p/g}}\hfill & \simeq \hfill & 0.13{\left(\frac{f_{\text{d/p}}}{1/2}\right)}^{-1}{\left(\frac{{\alpha }_{\text{D}z}/{\alpha }_{\text{acc}}}{{10}^{-2}}\right)}^{1/2}{\left(\frac{H_{\text{g}}/r}{0.04}\right)}^{1/2}{\left(\frac{\Delta x_{\text{subl}}}{0.1H_{\text{g}}}\right)}^{1/2}\hfill \\ \hfill & \hfill & \times \left[1+\frac{5}{3}\left(\frac{{\alpha }_{\text{D}r}/{\alpha }_{\text{acc}}}{{10}^{-2}}\right){\left(\frac{H_{\text{g}}/r}{0.04}\right)}^{-1}{\left(\frac{\Delta x_{\text{subl}}}{0.1H_{\text{g}}}\right)}^{-1}\right]\hfill \\ \hfill & \propto \hfill & \{\begin{array}{ll}{\left(\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}\right)}^{1/2}\left(\frac{{\alpha }_{\text{D}r}}{{\alpha }_{\text{acc}}}\right)\Delta x_{\text{subl}}^{-1}\hfill & [\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}>0.6\times {10}^{-2}],\hfill \\ {\left(\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}\right)}^{1/2}\Delta x_{\text{subl}}^{1/2}\hfill & [\frac{{\alpha }_{\text{D}z}}{{\alpha }_{\text{acc}}}<0.6\times {10}^{-2}].\hfill \end{array}\hfill \end{array}$$(40)

This equation explicitly shows that the pile-up condition is scaled by α_{D r} ∕α_acc and α_{D z} ∕α_acc.

In Eq. (40), the sublimation width, Δx_subl, is a free parameter. We used Δx = 0.1 H_g as a nominal parameter. The sublimation width was determined by the sublimation and drift rates of icy pebbles. The sublimation rate is given by the partial pressure of water vapor and the disk temperature. The water vapor partial pressure is affected by radial diffusion. In Paper II, we numerically solve the sublimation width, taking account of the icy pebble size evolution and the water vapor partial pressure in a turbulent accretion disk. Fitting the numerical results, the sublimation width for α_Dz = α_Dr and α_acc = 10⁻² is roughly approximated as

Fig. 9 Contour maps of the maximum ρd∕ρg on the plane of αDz∕αacc(= αDr∕αacc) and Fp/g with (a) Δxsubl= 0.1 Hg, and (b) Δxsubl given by Eq. (41). In both cases, αacc = 10−2 and K = 1 are adopted. Left panels: results of Monte Carlo simulations at t = 15 000 Ω−1, and right ones: results after using analytical formulas given by Eq. (3) with Eqs. (35), (36), (4), and the corresponding Δxsubl. The Monte Carlo simulations were performed with a grid size of Δlog10(αDz∕αacc) = Δlog10Fp/g = 0.2. We used more refined grids to make the maps by the analytical formulas. The color scales are based on log10 (ρd∕ρg).

$$\begin{array}{l}{\log }_{10}\left(\frac{\Delta x_{\text{subl}}}{H_{\text{g}}}\right)\simeq \frac{X_{+}-X_{-}}{2}\text{erf}\left[3{\log }_{10}\left(\frac{{\alpha }_{\text{D}r}/{\alpha }_{\text{acc}}}{0.08}\right)\right]+\frac{X_{+}+X_{-}}{2},\\ X_{-}={\log }_{10}{(\Delta x_{\text{subl}}/H_{\text{g}})}_{-}={\log }_{10}2;\\ X_{+}={\log }_{10}{(\Delta x_{\text{subl}}/H_{\text{g}})}_{+}={\log }_{10}0.1,\end{array}$$(41)

where ${(\Delta x_{\text{subl}}/H_{\text{g}})}_{-}=2$ and ${(\Delta x_{\text{subl}}/H_{\text{g}})}_{+}=0.1$ are the lower α_Dr∕α_acc limit at which advection dominates, and the higher one at which radial diffusion dominates, respectively (Paper II). The error function smoothly connects the two limits. Figures 9a and b show that the threshold value of F_p/g is several times higher than in Fig. 9a for α_Dz∕α_acc ~ 10⁻⁴, while it is similar for α_Dz∕α_acc ~ 10⁻². In the former case, $F_{\text{p/g}}\propto \Delta x_{\text{subl}}^{1/2}$ (Eq. (40)), so F_p/g must be larger in Fig. 9b than in Fig. 9a by a factor of ${(2/0.1)}^{1/2}~4.5$ (Eq. (41)), while the dependence on Δx_subl must be slightly weaker in the latter case (Eq. (41)). Thereby, the silicate particle runaway pile-up condition with Eq. (41) is roughly given by $$F_{\text{p/g}}\gtrsim {\left(\frac{{\alpha }_{\text{D}z}/{\alpha }_{\text{acc}}}{3\times {10}^{-2}}\right)}^{1/2},$$(42)

which agrees with the results in Fig. 9.

Figure 9 shows the results with α_Dr = α_Dz. As mentioned in Sect. 1, magneto-hydrodynamic simulations gave α_{D r} < α_Dz (Zhu et al. 2015; Yang et al. 2018). We also performed simulations with α_{D r} = 0.1 α_Dz. As we can predict from Eq. (40), the runaway pile-up region is extended to a value several times larger of α_{D z} ∕α_acc at F_p/g ~ 1 in the case of Δx_subl = 0.1 H_g, while the runaway pile-up region is almost the same for smaller values of F_p/g. Because α_{D z} ∕α_acc is very uncertain, we surveyed a rather broad parameter space and leave the estimate of α_{D z} ∕α_acc for future study.

4.4.3 Comparison with Hyodo et al. (2019)

Hyodo et al. (2019) performed 1D advection-diffusion grid code simulations with α_acc = 10⁻² and α_{D r} = α_Dz = 10⁻³, 3 × 10⁻³, 10⁻² to find that the silicate runaway pile-up occurs for α_Dz = α_Dr ≈ 10⁻¹α_acc and F_p/g ≥ 0.3 in the case of K = 1. However, Fig. 9b shows that the silicate runaway pile-up region is restricted to a smaller parameter range, that is, to values of α_{D z} ∕α_acc = α_Dr∕α_acc that are one order of magnitude smaller for a given F_p/g value. As already mentioned, this is because Hyodo et al. (2019) assumed that the silicate scale height is ~ Δz_subl at the snow line and it is gradually increased by the vertical turbulent stirring as the particles drift inward, whereas the Monte Carlo simulations performed here show that the silicate scale height is larger, due to the effect ofa coupled radial and vertical diffusion in a finite-sized sublimation width. As a result, smaller α_{D z} and α_{D r} are required for a runway pile-up.

5 Conclusion and discussions

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

We developed a new 2D (r-z) Monte Carlo code to simulate the pile-up of small silicate particles released from sublimating pebbles in a turbulent protoplanetary disk, taking account of the back reactions to the drift velocity and the diffusion of silicate particles. From the simulation results, we derived semi-analytical formulas for the maximum silicate-to-gas density ratio near the injection region and the silicate scale height there as a function of the pebble mass flux, α_acc, α_{D r}, and α_{D z}. Using the derived formulas, we determined the detailed condition for the silicate runaway pile-up. We found that the silicate particle scale height is larger than the estimate of Hyodo et al. (2019), due to the coupled effect of radial and vertical diffusion, and that for the silicate particle runaway pile-up to occur, α_{D z} ∕α_acc ≲ 10⁻²× (pebble-to-gas mass flux)² is required, which is more restrictive than the result of Hyodo et al. (2019) (α_{D z}∕α_acc ≲ 10⁻¹ and F_p/g ≳ 0.3). To clarify if the condition is actually satisfied, detailed nonideal MHD simulations are needed to evaluate α_{D z}, α_{D r}, and α_acc.

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

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

Acknowledgements

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

Appendix A The effects of back reaction to gas motion

The silicate dust particle and gas radial velocities taking account of the inertia of the particles (“back reaction”) at |z| < H_d are given, respectively, by (Schoonenberg & Ormel 2017) $$\begin{array}{lll}{v}_r\hfill & =\hfill & -{\Lambda }^2\frac{2{\tau }_{\text{s},\text{d}}}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\eta v_{\text{K}}+\Lambda \frac{1}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}{u}_{\nu },\hfill \\ u_r\hfill & =\hfill & Z{\Lambda }^2\frac{2{\tau }_{\text{s},\text{d}}}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\eta v_{\text{K}}+\Lambda \frac{1+\Lambda {\tau }_{\text{s},\text{d}}^2}{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}{u}_{\nu },\hfill \end{array}$$

where τ_s,d is the Stokes number of the particles, Z = ρ_d∕ρ_g, Λ = ρ_g∕(ρ_g + ρ_d) = 1∕(1 + Z), u_ν is an unperturbed disk gas accretion velocity given by $$u_{\nu }\simeq -\frac{3\nu }{2r}\simeq -\frac{3{\alpha }_{\text{acc}}}{2C_{\eta }}\eta v_{\text{K}},$$(A.3)

and $C_{\eta }=\eta /{(H_{\text{g}}/r)}^2$. When we consider the icy pebbles, the subscript “d” is replaced by “p.” Equations (A.1) and (A.2) are rewritten as $$\begin{array}{ll}{v}_r\hfill & \simeq -\frac{\Lambda }{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\left(2\Lambda {\tau }_{\text{s},\text{d}}+\frac{3}{2C_{\eta }}{\alpha }_{\text{acc}}\right)\eta v_{\text{K}},\hfill \\ u_r\hfill & \simeq -\frac{\Lambda }{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\left(-2Z\Lambda {\tau }_{\text{s},\text{d}}+\frac{3}{2C_{\eta }}(1+\Lambda {\tau }_{\text{s},\text{d}}^2){\alpha }_{\text{acc}}\right)\eta v_{\text{K}}.\hfill \end{array}$$

We adopted a two-layer model: a dust-rich midplane layer with scale height H_d and a dust-poor upper layer. To evaluate gas surface density, we used the vertically averaged u_r, such that $$Z_{\Sigma }=\frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}\simeq \frac{{\overline{u}}_r}{v_r}\frac{f_{\text{d/p}}{\dot{M}}_{\text{p}}}{{\dot{M}}_{\text{g}}}=\frac{{\overline{u}}_r}{v_r}{f}_{\text{d/p}}{F}_{\text{p/g}},$$(A.6)

where $$\begin{array}{ll}{\overline{u}}_r\hfill & \simeq \frac{u_r{|}_{z=0}{H}_{\text{d}}+u_{\nu }(H_{\text{g}}-H_{\text{d}})}{H_{\text{g}}}=h_{\text{d/g}}{u}_r{|}_{z=0}+(1-h_{\text{d/g}})u_{\nu }.\hfill \end{array}$$(A.7)

SubstitutingEq. (A.5) into this equation, we obtain $$\begin{array}{ll}{\overline{u}}_r\hfill & \simeq \frac{\Lambda }{1+{\Lambda }^2{\tau }_{\text{s},\text{d}}^2}\hfill \\ \hfill & \times \left[\left(\frac{2{\tau }_{\text{s},\text{d}}}{1+Z}+\frac{3{\alpha }_{\text{acc}}}{2C_{\eta }}\right)Z{h}_{\text{d/g}}-\frac{3{\alpha }_{\text{acc}}}{2C_{\eta }}(1+Z+\Lambda {\tau }_{\text{s},\text{d}}^2)\right]\eta v_{\text{K}}.\hfill \end{array}$$(A.8)

We note that even the vertically averaged gas motion is outward for $$\frac{{\alpha }_{\text{acc}}}{{\tau }_{\text{s},\text{d}}}<\frac{4C_{\eta }}{3}{h}_{\text{d/g}}\frac{Z}{(1+Z)[1+Z(1-h_{\text{d/g}})]},$$(A.9)

where we assumed ${\tau }_{\text{s},\text{d}}^2\ll 1$. For Z ≲ 1 and C_η = 11∕8, this condition is reduced to $$\frac{{\alpha }_{\text{acc}}}{{\tau }_{\text{s},\text{d}}}\lesssim \frac{11}{6}{h}_{\text{d/g}}Z.$$(A.10)

In the case of silicate dust particles, since τ_s,d ≪ α_acc, this condition is not satisfied and the gas flow is always inward.

For pebbles, the outward flow condition is the same as Eq. (A.10) with “d” replaced by “p”. The condition is less restrictive than that for silicate particles; however, it is not satisfied within the parameter range in this paper, as shown below. We usually consider the cases of τ_s,p > α_acc and $h_{\text{p/g}}\simeq {({\alpha }_{\text{D}z}/{\tau }_{\text{s},\text{p}})}^{1/2}$. In this case, the modified Eq. (A.10) is Z ≳ (6∕11)(α_acc∕α_Dz)h_p/g. If this equation is satisfied, the (vertically averaged) gas flow is outward beyond the snow line. However, in our case, Z, α_acc, and α_{D z} are not independent. Equation (A.16) is Z ≃ (1∕2)(3∕4C_η)(α_Dz∕α_acc)h_p/gF_p/g for Z < 1. Thus, the outward flow condition in the icy pebble region is $F_{\text{p/g}}\gtrsim 2{({\alpha }_{\text{acc}}/{\alpha }_{\text{D}z})}^2$. It is out of the parameter range that we cover in this paper. In the parameter regimes we cover, ū_r is always negative (inward flow) in the regions of both silicate particles and icy pebbles.

Substituting Eqs. (A.8) and (A.1) into Eq. (A.6), $$\begin{array}{ll}Z{h}_{\text{d/g}}\hfill & \simeq Z_{\Sigma }\simeq \frac{{\overline{u}}_r}{v_r}{f}_{\text{d/p}}{F}_{\text{p/g}}\hfill \\ \hfill & \simeq \frac{-2\frac{Z}{1+Z}{h}_{\text{d/g}}{\tau }_{\text{s},\text{d}}+\frac{3}{2C_{\eta }}{\alpha }_{\text{acc}}\left[1+Z(1-h_{\text{d/g}})+\frac{{\tau }_{\text{s},\text{d}}^2}{1+Z}\right]}{2\Lambda {\tau }_{\text{s},\text{d}}+\frac{3}{2C_{\eta }}{\alpha }_{\text{acc}}}{f}_{\text{d/p}}{F}_{\text{p/g}}.\hfill \end{array}$$(A.11)

Now we adopt the approximation only appropriate for the silicate dust, τ_s,d ≪ α_acc ≪ 1. With this approximation, Eq. (A.11) is reduced to $$\begin{array}{ll}Z{h}_{\text{d/g}}\hfill & \simeq \left[1+Z(1-h_{\text{d/g}})\right]f_{\text{d/p}}{F}_{\text{p/g}},\hfill \end{array}$$(A.12)

which is solved in terms of Z as $$\begin{array}{ll}Z\hfill & \simeq \frac{f_{\text{d/p}}{F}_{\text{p/g}}}{h_{\text{d/g}}-(1-h_{\text{d/g}})f_{\text{d/p}}{F}_{\text{p/g}}}.\hfill \end{array}$$(A.13)

As we discuss in Sect. 2, the pile-up of silicate particles is radially local and the corresponding local surface density variation of Σ_g could be smoothed out. In that case, it may be better to use u_r = u_ν than ū_r here, and we obtain $$\begin{array}{ll}Z\hfill & \simeq \frac{f_{\text{d/p}}{F}_{\text{p/g}}}{h_{\text{d/g}}-f_{\text{d/p}}{F}_{\text{p/g}}}.\hfill \end{array}$$(A.14)

Because h_d/g ≪ 1 for the pile-up case, the pileup condition differs only slightly between Eqs. (A.13) and (A.14).

From Eqs. (A.13) and (A.14), the runaway pile-up conditions of silicate particles are given by $$\begin{array}{lll}{F}_{\text{p/g}}\hfill & >\frac{f_{\text{d/p}}^{-1}}{1-h_{\text{d/g}}}{h}_{\text{d/g}}\text{or}>f_{\text{d/p}}^{-1}{h}_{\text{d/g}}\hfill & [\text{silicate}\text{particles}].\hfill \end{array}$$(A.15)

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

On the other hand, for pebbles, α_acc ≪ τ_s,p ≪ 1 is an appropriate approximation. From α_acc ≪ τ_s,p, we can also assume h_p/g ≪ 1. With these approximations, Eq. (A.11) with “d” replaced by “p” and f_d/p by 1 is reduced to $$\begin{array}{ll}Z{h}_{\text{p/g}}\hfill & \simeq \frac{3}{4C_{\eta }}\frac{{\alpha }_{\text{acc}}}{{\tau }_{\text{s},\text{p}}}{(1+Z)}^2{F}_{\text{p/g}}.\hfill \end{array}$$(A.16)

This quadratic equation does not have a solution if $${\left(1-\frac{2C_{\eta }}{3}\frac{{\tau }_{\text{s},\text{p}}}{{\alpha }_{\text{acc}}}\frac{h_{\text{p/g}}}{F_{\text{p/g}}}\right)}^2-1<0.$$(A.17)

Therefore, the runaway pile-up condition for pebbles is $$\begin{array}{ll}{F}_{\text{p/g}}\hfill & >\frac{C_{\eta }}{3}\frac{{\tau }_{\text{s},\text{p}}}{{\alpha }_{\text{acc}}}{h}_{\text{p/g}}.\hfill \end{array}$$(A.18)

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

Appendix B Consequences of outward diffusion over the snow line

In Eq. (15), the radial diffusion effect of silicate particles is not included. The (scaled) radial diffusion flux of silicate particles is $$F_{\text{d},\text{D}}~2\pi r{\alpha }_{\text{D}r}{H}_{\text{g}}^2\Omega \frac{\text{d}{\Sigma }_{\text{d}}}{\text{d}r}{\dot{M}}_{\text{g}}^{-1}~2\pi C_{F}r{\alpha }_{\text{D}r}{H}_{\text{g}}^2\Omega \frac{{\Sigma }_{\text{d}}}{H_{\text{g}}}{\dot{M}}_{\text{g}}^{-1},$$(B.1)

where C_F is an unknown numerical factor. We calibrated C_F by the results of Monte Carlo simulations in Figs. 7 and 8 as follows. The (scaled) net silicate dust mass flux, F_d,net ~ f_d/pF_p/g − F_d,D, in a steady state, should satisfy $F_{\text{d},\text{net}}~3\pi {\Sigma }_{\text{d}}{\alpha }_{\text{acc}}{H}_{\text{g}}^2\Omega {\dot{M}}_{\text{g}}^{-1}$. SubstitutingEq. (B.1) into this equation, we obtain $$\begin{array}{ll}{F}_{\text{d},\text{D}}\hfill & ~f_{\text{d/p}}{F}_{\text{p/g}}\times {\left[1+\frac{3}{2C_F}\frac{H_{\text{g}}}{r}\frac{{\alpha }_{\text{acc}}}{{\alpha }_{\text{D}r}}\right]}^{-1}.\hfill \end{array}$$(B.2)

Accordingly, the midplane particle-to-gas ratio at x ~ 0, ${({\rho }_{\text{d}}/{\rho }_{\text{g}})}_0$, is given by Eq. (15) with f_d/pF_p/g replaced by F_d,net = f_d/pF_p/g − F_d,D, as $${\left(\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\right)}_0=\frac{F_{\text{d},\text{net}}}{h_{\text{d/g}}-F_{\text{d},\text{net}}},$$(B.3)