Issue 
A&A
Volume 646, February 2021



Article Number  A14  
Number of page(s)  17  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202039894  
Published online  02 February 2021 
Planetesimal formation around the snow line
II. Dust or pebbles?
^{1}
ISAS/JAXA,
Sagamihara,
Kanagawa, Japan
email: hyodo@elsi.jp
^{2}
Université Côte d’Azur, Laboratoire J.L. Lagrange, CNRS, Observatoire de la Côte d’Azur,
06304 Nice, France
^{3}
EarthLife Science Institute, Tokyo Institute of Technology,
Meguroku,
Tokyo 1528550, Japan
^{4}
Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
Meguroku,
Tokyo 1528551, Japan
^{5}
Steward Observatory/The Lunar and Planetary Laboratory, University of Arizona,
Tucson,
Arizona 85721, USA
Received:
12
November
2020
Accepted:
11
December
2020
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 pileup of solids (silicate dust or icy pebbles) occurs inside or outside the snow line.
Methods. We performed 1D diffusionadvection simulations that include the backreaction (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 KelvinHelmholtz instability. We study the dependence of solid pileup 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 pebbletogas 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 pileup 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 backreaction that slows down the radial velocity of pebbles (we call this the “nodrift” region). We show that the pileup 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 pileup of silicate dust inside the snow line is favored for α_{Dr}∕α_{acc} ≪ 1, while that of pebbles outside the snow line is favored for α_{Dr}∕α_{acc} ~ 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 / planetdisk interactions / protoplanetary disks
© ESO 2021
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^{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; Dra̧ż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 pileup of solids around the water snow line (Saito & Sirono 2011; Morbidelli et al. 2015; Estrada et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dra̧ż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 solidtogas ratio: these are, for example, (1) gasdust friction for solids (pebbles and/or silicate dust), that is, the backreaction to their drift velocity and diffusive motion in the radial and vertical directions, and (2) a realistic prescription of the scale height of silicate dust as well as that of pebbles.
Backreactions that slow down drift velocities of pebbles and silicate dust as their pileup proceeds (hereafter, DriftBKR) are essential for pileup as Σ_{p (or d)} ∝ 1∕v_{p (or d)} (a subscript of “p (or d)” describes physical parameters of either pebbles or those of silicate dust), where Σ_{p (or d)} and v_{p (or d)} are the surface density and the radial velocity of pebbles or silicate dust, respectively (Ida & Guillot 2016). Backreaction to the diffusive motion (hereafter, DiffBKR) can trigger runaway pileup of solids as their diffusivities are progressively weakened as pileup 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 pileup. The midplane spatial density of silicate dust is ∝ 1∕H_{d} and thus the midplane solidtogas ratio is directly affected by H_{d}. The efficiency of the recycling (sticking) of diffused silicate dust onto pebbles outside the snow line is ∝ 1∕H_{d} as the number density of colliding silicate dust is ∝ 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 pileup 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 (radialvertical) 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 semianalytical 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). Thederived 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 pressuredependent sublimation around the snow line and we study the sublimation width of the drifting pebbles (Sect. 3).
Paper I also studied conditions for the pileup of silicate dust just inside the snow line and identified disk parameters where runaway pileup 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 pileup 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 pebbledominated midplane layer may interact with the upper/lower gasdominated layer, inducing a KelvinHelmholtz (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, pileups 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 DriftBKR and DiffBKR (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 pileups 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 pileup 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 backreaction 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.
Fig. 1 Summary of solid pileup around the snow line. Icy pebbles, which uniformly contain micronsized 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 socalled “trafficjam” 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 recyclesonto pebbles via the recondensation and sticking, causing a local pileup of pebbles just outside the snow line. Midplane concentrations (i.e., midplane solidtogas 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). 
2 Numerical methods and models
Here, we describe numerical methods and models. In this paper, we use three distinct nondimensional parameters to describe the gassolid evolutions in protoplanetary disk: α_{acc} which describes the global efficiency of angular momentum transport of gas (i.e., corresponds to the gas accretion rate) via turbulent viscosity ν_{acc}, α_{Dr} which describes the radial turbulence strength (i.e., corresponds to the radial velocity dispersion), and α_{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^{2}.
The classical αdisk model (Shakura & Sunyaev 1973) has adopted fully turbulent disks, i.e., α_{acc} ≃ α_{Dr} ≃ α_{Dz}, whereas recent theoretical developments suggest that the vertically averaged efficiency of the disk angular momentum transport, i.e., α_{acc}, could be different from local radial and vertical turbulence strengths α_{Dr} and α_{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 largescale magnetic fields (e.g., Bai et al. 2016). Magnetohydrodynamic (MHD) simulations have shown that a region of ionization that is too low for the magnetorotational 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 (α_{Dr}, α_{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 α_{Dr}, α_{Dz} < α_{acc} (see also Zhu et al. 2015; Hasegawa et al. 2017; Yang et al. 2018) and we study a wide parameter range where α_{Dr} (= α_{Dz}) ≤ α_{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 Ṁ_{g} and the gas radial velocity v_{g} (i.e., the classical αaccretion disk model; Shakura & Sunyaev 1973; LyndenBell & Pringle 1974) as (1)
The isothermal sound speed of the gas is written as (2)
where k_{B} is the Boltzmann constant, μ_{g} is the mean molecular weight of the gas (μ_{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 (3)
where we use T^{*} = 150 K and β = 1∕2. Then, the scale height of the gas is given as (4)
where Ω_{K} is the Keplerian orbital frequency. Disk gas pressure at the midplane P_{g} is given as (5)
The gas rotates at subKeplerian speed. The degree of deviation of the gas rotation frequency from that of Keplerian η is given by (6)
where Ω is the orbital frequency of the gas and we define (7)
which depends on the gas structure. In this paper, C_{η} = 11∕8 for Σ_{g} ∝ r^{−1} and T ∝ r^{−1∕2}.
2.2 Radial drifts of gas, pebbles, and silicate dust
The gas radial velocity including the effects of gassolid friction is given as (Ida & Guillot 2016; Schoonenberg & Ormel 2017) (8)
where Z is the midplane solidtogas ratio of pebbles or silicate dust, Λ characterizes the strength of the backreaction due to pileup of either pebbles or silicate dust, τ_{s} is the Stokes number of pebbles or silicate dust, and v_{g,ν} is an unperturbeddisk gas accretion velocity. Z is defined as (9)
where and 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. Λ is defined as (10)
Using H_{p} or H_{d}, Λ is rewritten as
where h_{p/g (or d/g)} ≡ H_{p (or d)}∕H_{g}. Z_{Σ} ≡ Σ_{p (or d)}∕Σ_{g} is the vertically averaged metallicities of pebbles or silicate dust. Using the dimensionless effective viscous parameter α_{acc} and the effective viscosity , v_{g,ν} is given as (13)
where (Desch et al. 2017). For a steadystate disk, the gas accretion velocity is and is rewritten as
where v_{K} is the Keplerian velocity.
In this work, we adopt a vertical twolayer model to describe the radial velocity of the gas: a pebble/dustrich midplane layer with scale height H_{p (or d)} and a pebble/dustpoor upper layer. The vertically averaged radial velocity of the gas v_{g} is given as (16)
We note that the above simple twolayer 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 gassolid friction − drift backreaction (DriftBKR) − are given as (Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019)
where τ_{s,p} and τ_{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
where ν_{Dr} and ν_{Dz} are the turbulent viscosities that regulate radial and vertical diffusions in association with dimensionless turbulence parameters α_{Dr} and α_{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 gassolid friction as pileup of solids proceeds − diffusion backreaction (DiffBKR) − and radial and vertical diffusivities of pebbles and silicate dust are given as (21) (22)
where K is the coefficient. K = 1 could be applied because the diffusivity is under a constant P assumption. Using the effective gas density, ρ_{g,eff} = ρ_{g} + ρ_{p (or d)}, the effective sound velocity of a mixture of gas and small particles is given as (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 pileup 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 micronsized 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 Σ_{g}, that of pebbles Σ_{p}, that of silicate dust Σ_{d}, that of water vapor Σ_{vap}, and the number density of pebbles N_{p} are given as (Desch et al. 2017; Hyodo et al. 2019);
The righthand 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} = Σ_{p}∕N_{p} under the singlesize approximation.
The decrease rate of the pebble mass m_{p} is given by (Lichtenegger & Komle 1991; Ros et al. 2019) (29)
where r_{p} is the particle physical radius, ρ_{vap} and ρ_{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 (30)
where μ is the mean molecular weight. Another definition of the thermal velocity is the averaged value of the magnitude of threedimensional velocity, which is given by (31)
Because thisdefinition 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, (32)
Both Eqs. (29) and (55) give the same d m_{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 ().
In the steadystate, 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^{3} (33)
However, forsmall α_{Dz} (i.e., for a small H_{p,tur}), a vertical shear KelvinHelmholtz (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 (seethe derivation in Sect. 4.1) as (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 (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 and the scaled sublimation width (Sect. 3) as (36)
where C_{r.diff} = 10, , and ϵ = 0.01 is a softening parameter that prevents unphysical peak for very small at the snow line (), respectively (Paper I). The scale height of pebbles is given at the snow line (with τ_{s,p} = 0.1) as (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 pileup 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 Ṁ_{p} to that of the gas Ṁ_{g} is a parameter and we define it as (40)
We set Ṁ_{g} = 10^{−8}M_{⊙} yr^{−1}. The inner and outer boundaries of our 1D simulations are set to r_{in} = 0.1 au and r_{out} = 5.0 au, respectively. We set the initial constant size of pebbles (radius of r_{p}) with τ_{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 micronsized 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 DriftBKR and DiffBKR (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 backreactions to reach a steadystate (~ 10^{5} yr) and then we turn on the backreactions (DriftBKR and DiffBKR with K = 1) to see the further evolution (up to 5 × 10^{5} yr).
3 Sublimation of drifting icy pebbles
A description of H_{d} is critical to understand the resultant pileup of silicate dust inside the snow line. Also, H_{d} describes the efficiency of recycling of silicate dust onto icy pebbles and pileup of pebbles just outside the snow line because its efficiency is (Hyodo et al. 2019) (see Appendix A for a detailed comparison of different models of H_{d}). Paper I derived a semianalytical expression of H_{d} that depends on the sublimation width Δ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 Δx_{subl} in a protoplanetary disk considering the local temperaturepressure environment. Here, we consider the cases of pure icy pebbles (i.e., f_{d/p} = 0 and τ_{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 differentsized 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 α_{Dr} ∕α_{acc} = 10^{−3} for pebbles and silicate dust, while we change α_{Dr}∕α_{acc} for the evolution of vapor. The arguments in this section are applicable to the 1:1 rockicemixed 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_{vap}(r) < P_{sat}(r) where P_{vap}(r) is the water vapor pressure and P_{sat} is the saturation vapor pressure (Fig. 2). P_{vap} > P_{sat} indicates a supersaturation state. The saturation vapor pressure of water is given as
Fig. 2 Schematic illustration of the sublimation of drifting icy pebbles. In the case of the advectiondominated 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 diffusiondominated regime (right; Sect. 3.2), water vapor produced inside the snow line efficiently diffuses outside the snow line, leading to a supersaturation 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}). 
where P_{0} = 1.14 × 10^{13} g cm^{−1} s^{−2} and T_{0} = 6062 K (Lichtenegger & Komle 1991), which can be rewritten as (42)
As shown in Fig. 2, the sublimation width of drifting pebbles can be divided into two regimes depending on α_{Dr} ∕α_{acc}: the advectiondominated regime (Sect. 3.1; Δx_{subl} ≃ 2H_{g} for α_{Dr} ≪ α_{acc}) and the diffusiondominated regime (Sect. 3.2; Δx_{subl} ≃ 0.2H_{g} for α_{Dr} ~ α_{acc}). Figure 3 shows the results of 1D simulations as a function of α_{Dr} ∕α_{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 Δ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 advectiondominated and diffusiondominated regimes are clearly identified by their different Δx_{subl} values in Fig. 3. We explain the details of these different regimes below.
As reference values, τ_{s,p} ~ 0.06, r_{p} ~ 2 cm at r_{snow} ~ 2.4 au for F_{p/g} = 0.1 and α_{acc} = 10^{−2}, while τ_{s,p} ~ 0.06, r_{p} ~ 20 cm at r_{snow} ~ 2.1 au for F_{p/g} = 0.1 and α_{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 Advectiondominated regime
In the advectiondominated regime (α_{Dr} ≪ α_{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) ≃ P_{vap}(r) outside the snow line (r > r_{snow}; see Fig. 2). Neglecting the effects of the radial diffusion implies that (43)
which is rewritten assuming a rapid vertical mixing with the background gas (the molecular weight μ_{vap}) as (44)
where T(r) ∝ r^{−β}, H_{g} ∝ r^{(3−β)∕2}, and ν_{acc}(r) ∝ α_{acc}r^{(3−2β)∕2}.
For r ≥ r_{snow}, the surface density of the water vapor that satisfies the local saturation state is given as (45)
Using the surface density of the water vapor and temperature at the snow line (i.e., Σ_{vap,snow} ≡ Σ_{vap}(r_{snow}) and T_{snow} ≡ T(r_{snow}), respectively), and Eq. (45), P_{0} is removed and Σ_{vap}(r) for r ≥ r_{snow} is given as (46)
Inside the snow line (i.e., r ≤ r_{snow}), Σ_{vap}(r) is written in terms of the gas mass flux Ṁ_{g} and the effective viscosity ν_{acc} as (47)
The above equation shows a good match with our 1D simulations for α_{Dr} ≪ α_{acc}, and thus the advection dominates the radial mass transfer of vapor for such a lowdiffusivity case. As discussed above, the local vapor pressure equals the saturation pressure as (48)
Using the vapor pressure at the snow line P_{vap,snow} ≡ P(r_{snow}) and removing P_{0}, the equation is rewritten as (49)
Using Δx^{*} = (r − r_{snow})∕r_{snow} ≪ 1 near the snow line (and Δx^{*} > 0), the above equation is approximated as follows:
where P_{vap}(r)∕P_{vap,snow} = 0 is the innermost radius where the sublimation takes place and P_{vap}(r)∕P_{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 Δx_{subl} is equivalent to the radial width of P_{vap}(r)∕P_{vap,snow} = 0.167 to 0.833, which corresponds to the change in mass of pebbles (). Using Eq. (52), is given by
where T_{snow} depends on α_{acc} and F_{p/g} (Fig. 3). For r_{snow} ≃ 2.4 au with β = 0.5 (T_{snow} ≃ 170 K), au (≃ 2.2H_{g}), which is consistent with the results of 1D simulations (the lightblue line in Fig. 3 bottom panel). Both 1D simulations and analytical arguments show that the sublimation width in the advectiondominated regime is very weakly dependent on F_{p/g} and α_{Dr} ∕α_{acc}. Equation (54) also indicates that the sublimation width in the advectiondominated regime is independent of the Stokes number of pebbles and initial pebble physical size, which are supported by 1D simulations (Fig. 3; r_{p} ~ 2 cm and ~ 20 cm for α_{acc} = 10^{−2} and 10^{−3}, respectively). These are because the sublimation width in the advectiondominated regime is regulated by P_{sat} and T_{snow}.
Fig. 3 Radial location of the snow line (top panel) and sublimation width of drifting pebbles Δx_{subl} (bottom panel) as a function of α_{Dr}∕α_{acc}. Here, f_{d/p} = 0 (i.e., pure icy pebbles). Circles (α_{acc} = 10^{−2}) and squares (α_{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 bythe gray lines (top panel). The lightblue line in the bottom panel shows the analytical Δx_{subl} in the case of the advectiondominated regime (Eq. (54) with r_{snow} = 2.4 au and β = 0.5). The lightgreen line in the bottom panel shows the analytical Δx_{subl} in the case of the diffusiondominated regime (Eq. (68) with τ_{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 betweenthe advection and diffusiondominated regimes (Eq. (73) with r = 2.4 au and β = 0.5). A fitting function is shown by a black line (Eq. (74)). 
3.2 Diffusiondominated regime
In the diffusiondominated regime (α_{Dr} ~ α_{acc}), the local sublimation of drifting pebbles does not take place outside of the snow line (r > r_{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 supersaturation state (i.e., P_{vap}(r) > P_{sat}(r) for r > r_{snow}). Thus, the sublimation only takes place after pebbles pass through the snow line where P_{vap}(r) < P_{sat} (i.e., r < r_{snow}).
The solar composition H/He gas sound velocity is given by cm s^{−1}. The water vapor sound velocity with the molecular weight μ = 18 is given by cm s^{−1}. Because the saturation pressure and vapor pressure are respectively given by and and , where ρ_{bulk} (~ 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)) (55)
For the disk density similar to MMSN, the snow line corresponds to T = T_{snow} ≃ 170 K. Using Eq. (44), the vapor pressure with 170 K at r = r_{snow} is P_{vap,snow} [dyn cm^{−2}] = 10^{13.06−2632∕170} ≃ 10^{−2.4} ≃ 3.8 × 10^{−3}. We assume that Σ ∝ r^{−βΣ}, T ∝ r^{−βT}, P_{vap} ∝ r^{−βP}.
We consider the case where the vapor pressure P_{vap} is proportional to disk gas pressure P_{H/He}. In that case,β_{P} = β_{Σ} + β_{T}∕2 + 3∕2 for ideal gas. We note that P_{vap} ∝ P_{H/He} is established in the limit of efficient radial diffusion (i.e., diffusiondominated regime). For a nominal case, β_{Σ} = 1, β_{T} = 1∕2, and β_{P} = 11∕4. Inside the snow line (r < r_{snow}), we can write
where T_{snow} ≃ 170 K was used. At r ~ r_{snow}, Eq. (55) is reduced to (58)
For the nominal parameters, 35.7β_{T} − β_{P} ≃ 15.1. The radial drift of a sublimating pebble is given by (61)
In the case of Λ ≃ 1 and , (62)
From Eqs. (60) and (62) with , (63) (64)
We further set (where ζ = 2 for Stokes drag and ζ = 1 for Epstein drag). Then,
We integrate the upper equation from as (67)
If we define the sublimation width Δx_{subl} by the radialseparation between and (as ), (68)
Equation (68) (with ζ = 1, τ_{s,p} = 0.06, r_{p} = 2.3 cm, r_{snow} = 2.4 au, ξ_{1} = 0.167, and ξ_{2} = 0.883) is shown by the lightgreen 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 diffusiondominated regime depends on F_{p/g}, in that, a larger F_{p/g} leads to a smaller Δx_{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_{sat}−P_{vap} (see Eq. (55)). In the case of α_{acc} = 10^{−3} (squares in Fig. 3), the snow line is located closer to the star (a larger P_{sat}−P_{vap}) than that of α_{acc} = 10^{−2}, which makes Δx_{subl} smaller, whereas the size of pebbles is 10 times larger than that of α_{acc} = 10^{−2}, which makes Δx_{subl} larger.
3.3 Boundary between advection/diffusiondominated regimes
The boundary between two regimes − the advectiondominated regime and the diffusiondominated regime − can be evaluated by considering the mass fluxes of advection and diffusion. In the steadystate for an unperturbed case (no backreaction), the rate of change in the surface density of the vapor (Eq. (27)) is described as: (69)
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_{vap}∕∂r ≃ ∂Σ_{vap}∕∂r). As discussed in the previous subsection (Sect. 3.1), P_{vap} satisfies the local saturation vapor (Eq. (48)) and (70)
From Eq. (69), the diffusion dominates over the advection when (71)
Therefore, the diffusion dominates over the advection when (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 advectiondominated regime (Sect. 3.1) and the diffusiondominated regime (Sect. 3.2) − which is a function of α_{Dr}∕α_{acc}. In the lower and higher α_{Dr}∕α_{acc} regions, advection and radial diffusion are dominated, respectively. Here, we provide a fitting function that smoothly connects the two regimes as (74)
Equation (74) is shown by a black line in Fig. 3 and is used to characterize H_{d} (Δx_{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 sublimate 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_{subl} leads to a smaller H_{d}). 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 pebbleconcentrated layer and the upper layer induces vorticity at the interface between the two fluids, i.e., KelvinHelmholtz (KH) instability (e.g., Sekiya 1998; Youdin & Shu 2002; Chiang 2008).
The azimuthal component of gas velocity (Schoonenberg & Ormel 2017) is
Because τ_{s} ≪ 1, α_{acc} ≪ 1, and v_{g,ϕ}≫v_{g,r}, v_{g,ϕ} is approximated as (77)
The Richardson number for pebbles at vertical height z is
where assumptions are made for , , and . We define C as (81)
Defining and using C and X, we can rewrite Eq. (80) as (82)
which is equivalently written as (83)
As Z ≫ C by definition and X ≪ 1, the solution is X ≃ C∕Z, that is, the scale height of pebbles regulated by a KH instability is written as
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} ≡ H_{p}∕H_{g}) is equivalent, within a factor order unity, to a classical description of H_{p} ∝ η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 ~ 0.75 is reported for Z < 1 (or Z_{Σ} < 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≳1 and the cusp becomes sharper as Z increases (Sekiya 1998; Chiang 2008). The height of the cusp (z_{max}) at which the particle distribution is truncated becomes constant in the limit of Z ≫ 1 (Chiang 2008), while Eq. (85) shows h_{p/g,KH} ∝ Z^{−1} for Z ≫ 1. On the other hand, the precise distribution is more similar to the Gaussian distribution for Z ≪ 1, Actually, in this case, z_{max} and H_{p} given by Eq. (85) are similar (only differ by ). 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_{max}.
Fig. 4 Midplane pebbletogas density ratio (Z = ρ_{p}∕ρ_{g}) as a functionof α_{Dz} and F_{p/g} (at r = 5 au and τ_{s,p} = 0.125). Top panels: cases where the backreaction onto the gas motion is neglected (i.e., v_{g} = v_{g,ν}). Bottom panels: cases where the backreaction 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 α_{acc} = 10^{−3} and α_{acc} = 10^{−2} with/without a KH instability (Ri = 0.5), respectively. The redcolored regions indicate where a steadystate solution is not found (i.e., ρ_{p} ∕ρ_{g} > 100 and it keeps increasing; the “nodrift (ND)” region). The white (Eq. (91)) and black (Eq. (95) with Ri = 0.5) dashed lines are analytically derived critical F_{p/g,crit} above which no steadystate are found (the “nodrift (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 backreaction to the gas motion that slows down the radial velocity of the gas, the midplane concentration of pebbles as well as the ND region (redcolored region) are shifted toward slightly larger F_{p/g} values. 
4.2 “Nodrift (ND)” region
DriftBKR 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 α_{Dz}) values, the drift velocity of pebbles can decrease in a selfinduced manner due to DriftBKR, 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 “nodrift (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 DriftBKR onto pebbles and gas. We numerically and analytically derive the degree of the midplane concentration of pebbles. The results are also used for ourinitial conditions at the outer boundary of 1D simulations. We identify the ND region (within the F_{p/g} − α_{acc} − α_{Dz} space).
The concentration of pebbles at the midplane is written as (86)
where v_{g}, v_{p}, and h_{p/g} are functions of Λ(Z) and Z = ρ_{p}∕ρ_{g} (Eqs. (16), (17), (35), and (10)). We numerically find a solution (i.e., solve for Z with τ_{s,p} = 0.1) for Eq. (86) for a given F_{p/g} and α_{Dz}. The color contours in Fig. 4 show the numerically obtained ρ_{p} ∕ρ_{g} values. The redcolored regions correspond to the parameters for which no steadystate solution is found (numerically, this corresponds to ρ_{p} ∕ρ_{g} ≫ 100), defined as the ND region for which pebbles no longer reach the snow line. We found that the boundaries of the “nodrift” regimes are characterized by a criticalvalue of ρ_{p}∕ρ_{g} = 1 (i.e., Z_{cri} = 1; see the color contours in Fig. 4). The dashed lines in Fig. 4 are analytically derived critical F_{p/g,crit} above which no steadystate solution is found for ρ_{p}∕ρ_{g} (i.e., no solution for Z). The analytical estimations well predict the direct solutions (the color contours). Below, we discuss that F_{p/g,crit} can be divided into two different regimes.
The first regime is when H_{p,tur} > H_{p,KH,max} where H_{p,KH,max} is the maximum scale height of pebbles regulated by a KHinstability (Sect. 4.1). In this case, using Eqs. (16) and (17), the vertically averaged metallicity of pebbles is given as
where approximation is made for α_{acc} ≪ τ_{s,p} ≪ 1, Z ≪ 1, Λ ≃1 (because Λ = 1 → 0 as pileup proceeds), and we use Eq. (11). Solving this equation gives (90)
where and , respectively. Z_{Σ} has a real solution when 1 − 4ab > 0 and a critical F_{p/g,crit1} for the first regime is given as
where for Z ≪ 1 and α_{Dz} ≪ τ_{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_{p,tur} < H_{p,KH,max}, where H_{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_{p,KH}, which is independentof α_{Dz}. As Z = ρ_{p}∕ρ_{g} increases from zero value, H_{p,KH} initially increases and reaches its maximum H_{p,KH,max} at Z = 1∕2 (see Sect. 4.1). As Z further increases, H_{p,KH} decreases. Using Z_{cri} = 1 as a critical value, the critical F_{p/g} in the second regime F_{p/g,crit2} adopts (H_{p,KH} with Z_{cri} = 1). Z_{cri} is given as
where and approximation is made for α_{acc} ≪ τ_{s,p} ≪ 1 and Λ ≃ 1. is Λ using . Here, Z = 1 and . Thus, the critical F_{p/g} for the second regime F_{p/g,crit2} with C_{η} = 11∕8 for T ∝ r^{−1∕2} is given as
which is independent on α_{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 α_{Dz} becomes small, reducing the parameter range of the “nodrift” 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 “nodrift” mechanism in protoplanetary disks is discussed in Hyodo et al. (2021).
5 Pileups 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 DriftBKR and DiffBKR onto the motions of pebbles and silicate dust. We choose K = 1 for DiffBKR. We change F_{p/g}, α_{acc}, and α_{Dr} = α_{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., ).
Figure 5 shows the maximum midplane solidtogas ratio around the water snow line in the F_{p/g}–α_{acc}–α_{Dz}(= α_{Dr}) space. The background color contours correspond to the results obtained for silicate dust inside the snow line in Paper I. The redcolored area indicates runaway pileup of silicate dust (i.e., no steadystate and an indication of planetesimal formation via direct gravitational collapse). The dependence on Δx_{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 ρ_{p}∕ρ_{g} (pebbles) and ρ_{d}∕ρ_{g} (silicate dust), respectively. When a runaway pileup occurs in our 1D simulations, it is labeled either by “PebRP” for pebbles or by “SilRP” for silicate dust. The numbers below the parenthesis in Fig. 5 show the rock fraction in pebbles. For a combination of a small α_{Dz} (= α_{Dr})/α_{acc} and a large F_{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 DriftBKR (we label these parameters by “ND” as the “nodrift” 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 pileup of silicate dust (Sect. 5.2), that of pebbles (Sect. 5.3), runaway pileups of pebbles and silicate dust (Sect. 5.4), rocktoice ratio within the pileup region (Sect. 5.5), and implication for planet formation (Sect. 5.6).
Fig. 5 Zones of pileup of pebbles and/or silicate dust around the snow line as a function of F_{p/g} and α_{Dz}(= α_{Dr})∕α_{acc} (case of a variable Δx_{subl}; Eq. (74) in Sect. 3.4). The numbers in the parenthesis represent the maximum ρ_{p} ∕ρ_{g} (left) and ρ_{d}∕ρ_{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 ρ_{p} ∕ρ_{g} (when a runaway pileup of silicate dust occurs, the fraction is not shown because the system would gravitationally collapse to form 100% rocky planetesimals). The runaway pileup of pebbles (labeled as “PebRP”) or that of silicate dust (labeled as “SilRP”) occurs around the snow line, depending on the combinations of F_{p/g} and α_{Dz}(= α_{Dr}). The color contour shows an analytical prediction of the maximum ρ_{d}∕ρ_{g} just inside the snow line obtained from Paper I. The black dashed line indicates analytically derived critical F_{p/g} above whichthe “nodrift (ND)” takes place (labeled by “ND” for the results of 1D simulations in the panels) at the outer boundary (r = 5 au and τ_{s,p} = 0.125; Sect. 4.2). DriftBKR and DiffBKR onto the motions of pebbles and silicate dust are included with K = 1. Left and right panels: cases of α_{acc} = 10^{−2} and α_{acc} = 10^{−3}, respectively. In the case of α_{acc} = 10^{−2}, the runaway pileup of silicate dust is inhibited because the required parameter regime is overlapped by the “nodrift” region (i.e., pebbles do not reach the snow line). In the case of α_{acc} = 10^{−3}, the runaway pileups of both pebbles and silicate dust occur. 
5.1 Comparison to previous works
In most previous studies of the pileup of silicate dust and pebbles around the snow line, α_{acc} = α_{Dr} = α_{Dz} was assumed (Dra̧żkowska & Alibert 2017; Schoonenberg & Ormel 2017). Figure 5 indicates that such parameters do not lead to a runaway pileup of silicate dust inside the snow line, while a pileup 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 (Dra̧żkowska & Alibert 2017; Schoonenberg & Ormel 2017) assumed K = 0 and a runaway pileup of pebbles was not reported (only a steadystate is reached with ρ_{p} ∕ρ_{g} > 1 for F_{p/g} > 0.8), while it could take place when K≠0 with DriftBKR (see also Hyodo et al. 2019).
5.2 Pileup 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; Dra̧ż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} × F_{p/g} (here, f_{d/p} = 0.5 is the silicate fraction in the original pebbles at the outer boundary) is independent on α_{acc} and α_{Dr}. The vertically integrated metallicity of silicate dust can be written as and is independent on α_{acc} as v_{d} ≃ v_{g}. The midplane solidtogas ratio is written as (h_{d/g} ≡ H_{d}∕H_{g}). At the snow line, h_{d/g} is regulated by α_{Dz} ∕α_{acc} (see Eq. (36) for α_{Dz} = α_{Dr}). Thus, both α_{acc} = 10^{−2} and α_{acc} = 10^{−3} cases have almost the same ρ_{d}∕ρ_{g} for the sameα_{Dz}∕α_{acc} and F_{p/g} when ρ_{d} ∕ρ_{g} ≪ 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 pileup is not significant. Paper I neglected the recycling of water vapor onto pebbles, but it does not significantly affect the pileup of silicate dust inside the snow line as the mass flux of silicate dust is a critical parameter for the pileup.
5.3 Pileup of pebbles outside the snow line
In this subsection, we discuss how pebbles pile up just outside the snow line and its dependence on α_{acc}. Outside the snow line, the pileup 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; Dra̧żkowska & Alibert 2017; Hyodo et al. 2019; Gárate et al. 2020). The efficiency of the recycling process (i.e., α_{Dr}) is a critical parameter to regulate the pileup of icy pebbles. Thus, unlike the case of silicate dust (Sect. 5.2), the value ρ_{p} ∕ρ_{g} is different for the same α_{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} ~ 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} ∝ 1∕v_{g} ∝ 1∕α_{acc}. This indicates that ρ_{p}∕ρ_{g} is ∝ α_{acc} when the recycling processes are negligible.
As F_{p/g} increases or/and α_{Dz}∕α_{acc} decreases, the pileup of pebbles tends to increase as the disk metallicity increases (Σ_{p}∝ F_{p/g}) and the scale height of pebbles decreases (Eqs. (35) and (36)). However, when α_{Dr} ∕α_{acc} ≃ 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 pileup 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} ~ 1 cases) because the recycling of the vapor and the dust onto pebbles efficiently enhances the pileup of pebbles (Hyodo et al. 2019). The pileup of silicate dust is more prominent when the diffusivity is weak (i.e., α_{Dz} ∕α_{acc} ≪ 1 cases).
5.4 Conditions for preferential pileup of silicate dust versus pebbles
Finally, we aim to understand under which conditions the runaway pileups 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 Ri = 0.5 are used here). In thecase of α_{acc} = 10^{−2}, the redcolored region is located above the black line, which implies that parameters in the F_{p/g} − α_{acc} − α_{Dz} space that would yield the runaway pileup of silicate dust found in Paper I in fact lie in the forbidden ND region. On the other hand, for α_{acc} = 10^{−3}, the pileup 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 pileups around the snow line occur over a broader range of parameters for α_{acc} = 10^{−3} than for α_{acc} = 10^{−2}. Runaway pileup of silicate dust is favored for a small α_{Dz}∕α_{acc} < 10^{−3} (labeled by “SilRP”), while runaway pileup of pebbles (labeled by “PebRP”) only occurs for a limited range of parameters (large F_{p/g} and α_{Dz} ∕α_{acc} ~ 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 Rocktoice ratio within the pileup region
Around the snow line, the rocktoice 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 rockice 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 ∝ 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 ∝ Σ_{d}).
DiffBKR and DriftBKR also play critical roles as these effects modify the diffusivity and radial motion of solids. As a result of a combinationof 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 pileup of pebbles occurs (labeled by PebRP in Fig. 5), the rock fraction becomes ~ 50 wt.%. This is because the runaway pileup occurs by accreting pebbles that drift from the outer region.
Thus, forming planetesimals from pebbles (by GI or SI) by runaway pileup 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 rockrich planetesimals even from the pileup 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 dusttoice ratio measured in comet 67P (O’Rourke et al. 2020). Planetesimals formed from the pileup 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 icerich and rockrich 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 andtherefore 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 MRIinactive “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} ~ 1) and the dead zone in the inner disk region (α_{Dr}(= α_{Dz})∕α_{acc} ≪ 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 pileup 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. Inbetween, 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 pileup 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 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.
Fig. 6 Illustration showing a possible planet formation scenario. Here we envision a disk containing an MRIinactive dead zone near the midplane (α_{Dr}(= α_{Dz})∕α_{acc} ≪ 1), while the outer disk is active (α_{Dr}(= α_{Dz})∕α_{acc} ~ 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 pileup (i.e., icy planetesimals formation by GI and/or SI; panel A) would be favored in the outer region (α_{Dz} ∕α_{acc} ~ 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 (α_{Dz} ∕α_{acc} ≪ 1; panel C). In the intermediate phase (α_{Dz}∕α_{acc} < 1 but not α_{Dz}∕α_{acc} ≪ 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. 
6 Summary
In this paper, to study the pileup of solids around the snow line, we incorporated a realistic H_{d} (obtained from Paper I) and the KHinstabilityconsidered H_{p} (Sect. 4.1), for the first time, into a 1D diffusionadvection code that includes the backreactions 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} − α_{acc} − α_{Dr}(= α_{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} − α_{acc} − α_{Dz} space, in which pebbles cannot reach the snow line by stopping their radial drift due to DriftBKR in a selfinduced manner (the “nodrift” region; Sect. 4.2).
Third, our 1D simulations showed that α_{acc} = 10^{−3} case is more favorable for the pileup of solids around the snow line than that for α_{acc} = 10^{−2}. The “nodrift” regime entirely covers parameter space of the runaway pileup of silicate dust when α_{acc} = 10^{−2} case (left panel in Fig. 5), preventing a runaway pileup of silicate dust inside the snow line. In contrast, in the α_{acc} = 10^{−3} case, large F_{p/g} values are allowed and could lead to a runaway pileup of silicate dust inside the snow line (right panel in Fig. 5). In the case of α_{acc} = 10^{−3}, both silicate dust and pebbles could experience runaway pileups inside and outside the snow line, forming rocky and icy planetesimals, respectively.
Forth, pebbles just outside the snow line is a rockice mixture (Sect. 5.5). When a runaway pileup of pebbles occurs just outside the snow line, the local rocktoice ratio becomes similar to that of pebbles formed at the outer region – i.e., the resultant planetesimal composition would be waterbearing 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 pileups 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. ChaoChin 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 AST1616929.
Appendix A Dependence of dust pileup upon H_{d}
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 1Dradial 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 , which is also used to describe the strength of the backreaction of the gas onto their motion (e.g., Eq. (18)). Therefore, the description of H_{d} is critical for the backreaction and to evaluate the pileup of dust in the midplane. In this section, we demonstrate how different models of H_{d} affect the dust pileup inside the snow line.
Fig. A.1 Scale height of silicate dust in a unit of that of the gas (top panels) and midplane dusttogas ratio (bottom panels) for different models of H_{d} (Case of F_{p/g} = 0.3 and τ_{s,p} = 0.1 at 3 au). Here, DriftBKR and DiffBKR for pebbles are neglected, whereas DriftBKR and DiffBKR (K = 0 or 1) for silicate dust are included. Left two panels: cases of α_{acc} = 10^{−2} and α_{Dr} = α_{Dz} = 10^{−4} with K = 0 or 1, respectively. Right two panels: cases of α_{acc} = 10^{−2} and α_{Dr} = α_{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 pileup when α_{acc} = 10^{−2} and α_{Dr} = α_{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), pileups occur in a runaway fashion for the left two panels, and the timeevolutions are shown for t = 1 × 10^{4} yr, t = 5 × 10^{4} yr, and t = 1.5 × 10^{5} yr. The other cases reach a steadystate within t = 5 × 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}. 
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 (A.1)
where a notation of “snow” indicates values at the snow line. Δt = (r − r_{snow})∕v_{d}. is the diffusion/mixing timescale to the vertical direction at the snow line and is the mean free path of turbulent blobs. Δt and t_{mix} are calculated by using the physical values at the snow line and we use τ_{s,p,snow} = 0.1.
Figure A.1 shows different models of the scale height of silicate dust and their resultant pileups 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 pileup of silicate dust when α_{Dz} ≪ α_{acc}, while it correctly evaluated when α_{Dz} ≃ α_{acc}. Ida & Guillot (2016), H_{d} = H_{p,snow}, overestimated the pileup of silicate dust regardless of the value of α_{Dz}. Hyodo et al. (2019) overestimated pileup of silicate dust near the snow line.
Appendix B Particle size and Stokes number
The Stokes number of a particle τ_{s}(r) depends on the local gas structure and τ_{s}(r) is written byusing the stopping time t_{s} as τ_{s}(r) = t_{s}(r)Ω_{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 τ_{s}(r) is given as
where is the thermal velocity, r_{p} is the particle radius, and ρ_{p} = 1.5 g cm^{−3} is the particle internal density, and is the gas spatial density, respectively. The mean free path of the gas l_{mfp}(r) is given as
Fig. B.1 Particle size as a function of distance to the star. The red, green, and blue lines represent for α_{acc} = 10^{−2}, 10^{−3}, and 10^{−4}, respectively (Eq. (B.3)). The Stokes regime is where , while the Epstein regime is where . Four black and gray lines represent particle size whose τ_{s} = 0.3, 0.1, 0.03, and 0.01 in the Stokesregime and the Epstein regime with α_{acc} = 10^{−2}, respectively. Here, Ṁ_{g} = 10^{−8}M_{⊙} yr^{−1} and are used. 
where σ_{mol} = 2.0 × 10^{−15} 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 α_{acc}.
Rewriting Eqs. (B.1) and (B.2) give particle sizes r_{p} in the Epstein and Stokes regimes for a fixed τ_{s}. Figure B.1 shows critical particle sizes of for different α_{acc} and particle sizes in either Epstein and Stokes regimes for a fixed τ_{s}. Because , , , , and , the particle size in the Epstein regime becomes , that is, r_{p} ∝ r^{−1} for β = 1∕2 (the gray lines in Fig. B.1). Because , the particle size in the Stokes regime becomes , which is independent on α_{acc}, that is, for β = 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 Δ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 (Δx_{subl} = 0.1H_{g}).
As the sublimation width Δ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 α_{Dr} (= α_{Dz})∕α_{acc} (i.e., the advectiondominated regime) as both radial and vertical mixings are ineffective, while large α_{Dr} (= α_{Dz})∕α_{acc} (i.e., the diffusiondominated 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 pileup of silicate dust is reduced for small α_{Dr} (= α_{Dz})∕α_{acc} in the case of a variable Δx_{sub} (Eq. (74)) compared to the case of the fixed Δ_{x} = 0.1H_{g} (compare Fig. C.1 to Fig. 5).
References
 Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 769, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chiang, E. 2008, ApJ, 675, 1549 [CrossRef] [Google Scholar]
 Chiang, E., & Youdin, A. N. 2010, Ann. Rev. Earth Planet. Sci., 38, 493 [Google Scholar]
 Ciesla, F. J., & Cuzzi, J. N. 2006, Icarus, 181, 178 [NASA ADS] [CrossRef] [Google Scholar]
 Desch, S. J., Estrada, P. R., Kalyaan, A., & Cuzzi, J. N. 2017, ApJ, 840, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Dra̧żkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Estrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Gárate, M., Birnstiel, T., Drążkowska, J., & Stammler, S. M. 2020, A&A, 635, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gerbig, K., MurrayClay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91 [CrossRef] [Google Scholar]
 Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Gundlach, B., Schmidt, K. P., Kreuzig, C., et al. 2018, MNRAS, 479, 1273 [NASA ADS] [CrossRef] [Google Scholar]
 Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hyodo, R., Ida, S., & Guillot, T. 2021, A&A, 645, L9 [CrossRef] [EDP Sciences] [Google Scholar]
 Ida, S., & Guillot, T. 2016, A&A, 596, L3 [Google Scholar]
 Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ida, S., Guillot, T., Hyodo, R., Okuzumi, S., & Youdin, A. N. 2021, A&A, 646, A13 [EDP Sciences] [Google Scholar]
 Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Kimura, H., Wada, K., Senshu, H., & Kobayashi, H. 2015, ApJ, 812, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Lichtenegger, H. I. M., & Komle, N. I. 1991, Icarus, 90, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, B., Lambrechts, M., Johansen, A., & Liu, F. 2019, A&A, 632, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Mori, S., Muranushi, T., Okuzumi, S., & Inutsuka, S.i. 2017, ApJ, 849, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Musiolik, G., & Wurm, G. 2019, ApJ, 873, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Momose, M., Sirono, S.i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82 [NASA ADS] [CrossRef] [Google Scholar]
 O’Rourke, L., Heinisch, P., Blum, J., et al. 2020, Nature, 586, 697 [CrossRef] [Google Scholar]
 Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saito, E., & Sirono, S.i. 2011, ApJ, 728, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sekiya, M. 1998, Icarus, 133, 298 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117 [NASA ADS] [CrossRef] [Google Scholar]
 Steinpilz, T., Teiser, J., & Wurm, G. 2019, ApJ, 874, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yang, C.C., Mac Low, M.M., & Johansen, A. 2018, ApJ, 868, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Stone, J. M., & Bai, X.N. 2015, ApJ, 801, 81 [NASA ADS] [CrossRef] [Google Scholar]
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).
In Hyodo et al. (2019), α_{tur} is used to describe α_{Dr} and α_{Dz}.
Previous studies (Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017; Hyodo et al. 2019) did not consider the effects of DiffBKR on the scale height of pebbles (i.e., K = 0 is assumed for Eq. (33)).
Even when the backreaction to the gas motion is neglected (i.e., v_{g} = v_{g,ν} as Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017; Hyodo et al. 2019), F_{p/g, crit1} is approximated to (i.e., replace the denominator of Eq. (88) with v_{g,ν} and approximation is made for α_{acc} ≪ τ_{s,p} ≪ 1, Z ≪ 1, Λ ≃ 1).
When ρ_{d}∕ρ_{g} > 1, the backreaction changes other quantities, such as the scale height of pebbles (Eq. (35)).
All Figures
Fig. 1 Summary of solid pileup around the snow line. Icy pebbles, which uniformly contain micronsized 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 socalled “trafficjam” 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 recyclesonto pebbles via the recondensation and sticking, causing a local pileup of pebbles just outside the snow line. Midplane concentrations (i.e., midplane solidtogas 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). 

In the text 
Fig. 2 Schematic illustration of the sublimation of drifting icy pebbles. In the case of the advectiondominated 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 diffusiondominated regime (right; Sect. 3.2), water vapor produced inside the snow line efficiently diffuses outside the snow line, leading to a supersaturation 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}). 

In the text 
Fig. 3 Radial location of the snow line (top panel) and sublimation width of drifting pebbles Δx_{subl} (bottom panel) as a function of α_{Dr}∕α_{acc}. Here, f_{d/p} = 0 (i.e., pure icy pebbles). Circles (α_{acc} = 10^{−2}) and squares (α_{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 bythe gray lines (top panel). The lightblue line in the bottom panel shows the analytical Δx_{subl} in the case of the advectiondominated regime (Eq. (54) with r_{snow} = 2.4 au and β = 0.5). The lightgreen line in the bottom panel shows the analytical Δx_{subl} in the case of the diffusiondominated regime (Eq. (68) with τ_{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 betweenthe advection and diffusiondominated regimes (Eq. (73) with r = 2.4 au and β = 0.5). A fitting function is shown by a black line (Eq. (74)). 

In the text 
Fig. 4 Midplane pebbletogas density ratio (Z = ρ_{p}∕ρ_{g}) as a functionof α_{Dz} and F_{p/g} (at r = 5 au and τ_{s,p} = 0.125). Top panels: cases where the backreaction onto the gas motion is neglected (i.e., v_{g} = v_{g,ν}). Bottom panels: cases where the backreaction 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 α_{acc} = 10^{−3} and α_{acc} = 10^{−2} with/without a KH instability (Ri = 0.5), respectively. The redcolored regions indicate where a steadystate solution is not found (i.e., ρ_{p} ∕ρ_{g} > 100 and it keeps increasing; the “nodrift (ND)” region). The white (Eq. (91)) and black (Eq. (95) with Ri = 0.5) dashed lines are analytically derived critical F_{p/g,crit} above which no steadystate are found (the “nodrift (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 backreaction to the gas motion that slows down the radial velocity of the gas, the midplane concentration of pebbles as well as the ND region (redcolored region) are shifted toward slightly larger F_{p/g} values. 

In the text 
Fig. 5 Zones of pileup of pebbles and/or silicate dust around the snow line as a function of F_{p/g} and α_{Dz}(= α_{Dr})∕α_{acc} (case of a variable Δx_{subl}; Eq. (74) in Sect. 3.4). The numbers in the parenthesis represent the maximum ρ_{p} ∕ρ_{g} (left) and ρ_{d}∕ρ_{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 ρ_{p} ∕ρ_{g} (when a runaway pileup of silicate dust occurs, the fraction is not shown because the system would gravitationally collapse to form 100% rocky planetesimals). The runaway pileup of pebbles (labeled as “PebRP”) or that of silicate dust (labeled as “SilRP”) occurs around the snow line, depending on the combinations of F_{p/g} and α_{Dz}(= α_{Dr}). The color contour shows an analytical prediction of the maximum ρ_{d}∕ρ_{g} just inside the snow line obtained from Paper I. The black dashed line indicates analytically derived critical F_{p/g} above whichthe “nodrift (ND)” takes place (labeled by “ND” for the results of 1D simulations in the panels) at the outer boundary (r = 5 au and τ_{s,p} = 0.125; Sect. 4.2). DriftBKR and DiffBKR onto the motions of pebbles and silicate dust are included with K = 1. Left and right panels: cases of α_{acc} = 10^{−2} and α_{acc} = 10^{−3}, respectively. In the case of α_{acc} = 10^{−2}, the runaway pileup of silicate dust is inhibited because the required parameter regime is overlapped by the “nodrift” region (i.e., pebbles do not reach the snow line). In the case of α_{acc} = 10^{−3}, the runaway pileups of both pebbles and silicate dust occur. 

In the text 
Fig. 6 Illustration showing a possible planet formation scenario. Here we envision a disk containing an MRIinactive dead zone near the midplane (α_{Dr}(= α_{Dz})∕α_{acc} ≪ 1), while the outer disk is active (α_{Dr}(= α_{Dz})∕α_{acc} ~ 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 pileup (i.e., icy planetesimals formation by GI and/or SI; panel A) would be favored in the outer region (α_{Dz} ∕α_{acc} ~ 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 (α_{Dz} ∕α_{acc} ≪ 1; panel C). In the intermediate phase (α_{Dz}∕α_{acc} < 1 but not α_{Dz}∕α_{acc} ≪ 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. 

In the text 
Fig. A.1 Scale height of silicate dust in a unit of that of the gas (top panels) and midplane dusttogas ratio (bottom panels) for different models of H_{d} (Case of F_{p/g} = 0.3 and τ_{s,p} = 0.1 at 3 au). Here, DriftBKR and DiffBKR for pebbles are neglected, whereas DriftBKR and DiffBKR (K = 0 or 1) for silicate dust are included. Left two panels: cases of α_{acc} = 10^{−2} and α_{Dr} = α_{Dz} = 10^{−4} with K = 0 or 1, respectively. Right two panels: cases of α_{acc} = 10^{−2} and α_{Dr} = α_{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 pileup when α_{acc} = 10^{−2} and α_{Dr} = α_{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), pileups occur in a runaway fashion for the left two panels, and the timeevolutions are shown for t = 1 × 10^{4} yr, t = 5 × 10^{4} yr, and t = 1.5 × 10^{5} yr. The other cases reach a steadystate within t = 5 × 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}. 

In the text 
Fig. B.1 Particle size as a function of distance to the star. The red, green, and blue lines represent for α_{acc} = 10^{−2}, 10^{−3}, and 10^{−4}, respectively (Eq. (B.3)). The Stokes regime is where , while the Epstein regime is where . Four black and gray lines represent particle size whose τ_{s} = 0.3, 0.1, 0.03, and 0.01 in the Stokesregime and the Epstein regime with α_{acc} = 10^{−2}, respectively. Here, Ṁ_{g} = 10^{−8}M_{⊙} yr^{−1} and are used. 

In the text 
Fig. C.1 Same as Fig. 5, but for the case of a constant Δx_{subl} = 0.1H_{g}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.