© 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¹ (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 pile-up of solids around the water snow line (Saito & Sirono 2011; Morbidelli et al. 2015; Estrada et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; 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 solid-to-gas ratio: these are, for example, (1) gas-dust friction for solids (pebbles and/or silicate dust), that is, the back-reaction to their drift velocity and diffusive motion in the radial and vertical directions, and (2) a realistic prescription of the scale height of silicate dust as well as that of pebbles.

Back-reactions that slow down drift velocities of pebbles and silicate dust as their pile-up proceeds (hereafter, Drift-BKR) are essential for pile-up as Σ_{p (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). Back-reaction to the diffusive motion (hereafter, Diff-BKR) can trigger runaway pile-up of solids as their diffusivities are progressively weakened as pile-up proceeds (Hyodo et al. 2019, see Eqs. (21) and (22) below).

The scale height of silicate dust H_d is another critical consideration on the solid pile-up. The midplane spatial density of silicate dust is ∝ 1∕H_d and thus the midplane solid-to-gas ratio is directly affected by H_d. The efficiency of the recycling (sticking) of diffused silicate dust onto pebbles outside the snow line is ∝ 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 pile-up of solids, depending on their assumed H_d and the disk conditions (see more details in Appendix A).

In our companion paper (Ida et al. 2021; hereafter Paper I), we performed 2D (radial-vertical) Monte Carlo simulations (a Lagrange method) of silicate dust that was released around the snow line from drifting icy pebbles under ice sublimation. Based on these Monte Carlo simulations and using semi-analytical arguments, Paper I derived a formula for the evolving scale height of the silicate dust for a given disk parameters as a function of the distance to the snow line (Eq. (36) below). 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 pressure-dependent sublimation around the snow line and we study the sublimation width of the drifting pebbles (Sect. 3).

Paper I also studied conditions for the pile-up of silicate dust just inside the snow line and identified disk parameters where runaway pile-up of silicate dust occurs (i.e., rocky planetesimal formation). However, Paper I did not include recycling of water vapor and of silicate dust outside the snow line due to their diffusive motion and hence they did not consider pile-up process of pebbles outside the snow line.

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

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

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

Fig. 1 Summary of solid pile-up around the snow line. Icy pebbles, which uniformly contain micron-sized silicate dust formed at the outer disk, drift inward due to the gas drag. Pebbles sublimate through the passage of the snow line (Sect. 3) and silicate dust is released together with the water vapor. Silicate dust is well coupled to the gas and the radial drift velocity significantly drops from that of pebbles. This causes so-called “traffic-jam” effect and the dust piles up inside the snow line. A fraction of water vapor and silicate dust diffuses from inside to outside the snow line recyclesonto pebbles via the recondensation and sticking, causing a local pile-up of pebbles just outside the snow line. Midplane concentrations (i.e., midplane solid-to-gas ratio) of pebbles and silicate dust strongly depend on their scale heights. Near the snow line, the scale height of silicate dust is the minimum as the dust is released from pebbles whose scale height is much smaller than that of the gas. This causes the maximum midplane concentration of silicate dust just inside the snow line. Significant midplane concentrations of silicate dust and/or icy pebbles would cause gravitational instability and/or streaming instability, forming rocky and/or icy planetesimals, respectively (Sect. 5).

2 Numerical methods and models

Here, we describe numerical methods and models. In this paper, we use three distinct non-dimensional parameters to describe the gas-solid evolutions in protoplanetary disk: α_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².

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 large-scale magnetic fields (e.g., Bai et al. 2016). Magnetohydrodynamic (MHD) simulations have shown that a region of ionization that is too low for the magneto-rotational instability (MRI) to operate (“dead zone”; Gammie 1996) might ubiquitously exist in the inner part of disk midplane and only surface layers are magnetically active (α_Dr, α_Dz = 10⁻⁵−10⁻³ 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; Lynden-Bell & 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 sub-Keplerian 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⁻¹ and T ∝ r^−1∕2.

2.2 Radial drifts of gas, pebbles, and silicate dust

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

where Z is the midplane solid-to-gas ratio of pebbles or silicate dust, Λ characterizes the strength of the back-reaction due to pile-up 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 steady-state disk, the gas accretion velocity is and is rewritten as

where v_K is the Keplerian velocity.

In this work, we adopt a vertical two-layer model to describe the radial velocity of the gas: a pebble/dust-rich midplane layer with scale height H_{p (or d)} and a pebble/dust-poor upper layer. The vertically averaged radial velocity of the gas v_g is given as (16)

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

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

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 gas-solid friction as pile-up of solids proceeds − diffusion back-reaction (Diff-BKR) − and radial and vertical diffusivities of pebbles and silicate dust are given as (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 pile-up around the snow line does not significantly depend on the choice of K as long as the back reaction to diffusion is included (i.e., K > 0). In this work, we use K = 1 as a representative case.

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

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

The governing equations of the surface density of the gas Σ_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 right-hand sides of Eqs. (25)–(27) are due to sublimation of icy pebbles, the recondensation of water vapor, and sticking of silicate dust onto pebbles (see Hyodo et al. 2019). The evolution of N_p needs to be solved because the size of pebbles changes during sublimation and condensation, and the mass of pebbles is given as m_p = Σ_p∕N_p under the single-size approximation.

The decrease rate of the pebble mass m_p is given by (Lichtenegger & Komle 1991; Ros et al. 2019) (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 three-dimensional 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 steady-state, the scale height of pebbles is regulated by the vertical turbulent stirring H_p,tur (Dubrulle et al. 1995; Youdin & Lithwick 2007; Okuzumi et al. 2012) as³ (33)

However, forsmall α_Dz (i.e., for a small H_p,tur), a vertical shear Kelvin-Helmholtz (KH) instability would prevent the scale height of pebbles from being a smaller value. The scale height of pebbles regulated by a KH instability H_{p, KH} is given (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 (37)

and (38)

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 pile-up of silicate dust and pebbles (Appendix A). Equations (36)–(38) provide a more realistic treatment based on the 2D calculation from Paper I.

2.6 Numerical settings

Following Paper I, the ratio of the mass flux of pebbles Ṁ_p to that of the gas Ṁ_g is a parameter and we define it as (40)

We set Ṁ_g = 10⁻⁸M_⊙ yr⁻¹. 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 micron-sized small grains (Schoonenberg & Ormel 2017; Hyodo et al. 2019; Ida et al. 2021). The surface density of pebbles at the outer boundary is set by considering Drift-BKR and Diff-BKR (see Sect. 4.2 and Eq. (86)), which is fixed during the simulations. We consider sublimation/condensation of ice as well as release/recycling of silicate dust around the snow line (Hyodo et al. 2019). 1D simulations are initially run without back-reactions to reach a steady-state (~ 10⁵ yr) and then we turn on the back-reactions (Drift-BKR and Diff-BKR with K = 1) to see the further evolution (up to 5 × 10⁵ yr).

3 Sublimation of drifting icy pebbles

A description of H_d is critical to understand the resultant pile-up of silicate dust inside the snow line. Also, H_d describes the efficiency of recycling of silicate dust onto icy pebbles and pile-up of pebbles just outside the snow line because its efficiency is (Hyodo et al. 2019) (see Appendix A for a detailed comparison of different models of H_d). Paper I derived a semi-analytical expression of H_d that depends on the sublimation width Δ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 temperature-pressure 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 different-sized sublimating pebbles near the snow line, which comes from solving N_d evolution (i.e., the effective size of pebbles at a given radial distance is calculated). In this section, to suppress these mixing effects and to see the size evolution of a single drifting pebble, we set α_Dr ∕α_acc = 10⁻³ 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 rock-ice-mixed pebbles (see Sect. 2.6) because the size and density of pebbles are expected to only weakly change with slightly different compositions.

The sublimation of drifting icy pebbles takes place as long as P_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 super-saturation 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 advection-dominated regime (left; Sect. 3.1), gas accretion onto the central star dominates the radial transport of vapor. The loss of vapor due to the inward advection of the vapor then balances the vapor supply from sublimating icy pebbles. Thus, drifting pebbles sublimate outside the snow line (i.e., r > rsnow). In the case of the diffusion-dominated regime (right; Sect. 3.2), water vapor produced inside the snow line efficiently diffuses outside the snow line, leading to a super-saturation state (Pvap (r) > Psat(r)). The sublimation of pebbles can therefore take place only inside the snow line (i.e., r < rsnow).

(41)

where P₀ = 1.14 × 10¹³ g cm⁻¹ s⁻² and T₀ = 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 advection-dominated regime (Sect. 3.1; Δx_subl ≃ 2H_g for α_Dr ≪ α_acc) and the diffusion-dominated 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 advection-dominated and diffusion-dominated 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⁻², while τ_s,p ~ 0.06, r_p ~ 20 cm at r_snow ~ 2.1 au for F_p/g = 0.1 and α_acc = 10⁻³. Here, the drag relations correspond to those of the Epstein regime. Below, a notation of “snow” indicates values at the snow line.

3.1 Advection-dominated regime

In the advection-dominated regime (α_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) ∝ α_accr^{(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₀ 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 low-diffusivity 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₀, 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 light-blue line in Fig. 3 bottom panel). Both 1D simulations and analytical arguments show that the sublimation width in the advection-dominated regime is very weakly dependent on F_p/g and α_Dr ∕α_acc. Equation (54) also indicates that the sublimation width in the advection-dominated regime is independent of the Stokes number of pebbles and initial pebble physical size, which are supported by 1D simulations (Fig. 3; r_p ~ 2 cm and ~ 20 cm for α_acc = 10⁻² and 10⁻³, respectively). These are because the sublimation width in the advection-dominated 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 Δxsubl (bottom panel) as a function of αDr∕αacc. Here, fd/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 Fp/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 light-blue line in the bottom panel shows the analytical Δxsubl in the case of the advection-dominated regime (Eq. (54) with rsnow = 2.4 au and β = 0.5). The light-green line in the bottom panel shows the analytical Δxsubl in the case of the diffusion-dominated regime (Eq. (68) with τs,p = 0.06, rp = 2.3 cm, rsnow = 2.4 au). The gray dashed line in the bottom panel shows the analytically derived boundary betweenthe advection- and diffusion-dominated regimes (Eq. (73) with r = 2.4 au and β = 0.5). A fitting function is shown by a black line (Eq. (74)).

3.2 Diffusion-dominated regime

In the diffusion-dominated 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 super-saturation 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⁻¹. The water vapor sound velocity with the molecular weight μ = 18 is given by cm s⁻¹. Because the saturation pressure and vapor pressure are respectively given by and and , where ρ_bulk (~ 1 g cm⁻³) 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⁻²] = 10^{13.06−2632∕170} ≃ 10^−2.4 ≃ 3.8 × 10⁻³. 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., diffusion-dominated 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)

where (59)

Defining , (60)

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, ξ₁ = 0.167, and ξ₂ = 0.883) is shown by the light-green line in the bottom panel of Fig. 3. The analytical arguments above are generally in accordance with the 1D numerical results. The sublimation width in the diffusion-dominated regime depends on F_p/g, in that, a larger F_p/g leads to a smaller Δ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⁻³ (squares in Fig. 3), the snow line is located closer to the star (a larger P_sat−P_vap) than that of α_acc = 10⁻², which makes Δx_subl smaller, whereas the size of pebbles is 10 times larger than that of α_acc = 10⁻², which makes Δx_subl larger.

3.3 Boundary between advection/diffusion-dominated regimes

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

which is rewritten as (72)

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 advection-dominated regime (Sect. 3.1) and the diffusion-dominated 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 pebble-concentrated layer and the upper layer induces vorticity at the interface between the two fluids, i.e., Kelvin-Helmholtz (KH) instability (e.g., Sekiya 1998; Youdin & Shu 2002; Chiang 2008).

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

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⁻¹ 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 pebble-to-gas density ratio (Z = ρp∕ρg) as a functionof αDz and Fp/g (at r = 5 au and τs,p = 0.125). Top panels: cases where the back-reaction onto the gas motion is neglected (i.e., vg = vg,ν). Bottom panels: cases where the back-reaction to the gas motion is included, following Eq. (16). The color contours are obtained by directly solving Eq. (86). Left and right two panels: cases where αacc = 10−3 and αacc = 10−2 with/without a KH instability (Ri = 0.5), respectively. The red-colored regions indicate where a steady-state solution is not found (i.e., ρp ∕ρg > 100 and it keeps increasing; the “no-drift (ND)” region). The white (Eq. (91)) and black (Eq. (95) with Ri = 0.5) dashed lines are analytically derived critical Fp/g,crit above which no steady-state are found (the “no-drift (ND)” region). Analytical predictions (dashed lines) well predict critical boundaries of the ND region. Including a KH instability reduces the ND region as the minimum scale height of pebbles is regulated by a KH instability. Due to the back-reaction to the gas motion that slows down the radial velocity of the gas, the midplane concentration of pebbles as well as the ND region (red-colored region) are shifted toward slightly larger Fp/g values.

4.2 “No-drift (ND)” region

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

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

Below, we discuss the concentration of icy pebbles in the disk midplane considering Drift-BKR onto pebbles and gas. We numerically and analytically derive the degree of the midplane concentration of pebbles. The results are also used for 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 red-colored regions correspond to the parameters for which no steady-state 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 “no-drift” 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 steady-state 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 KH-instability (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 pile-up 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)⁴.

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 “no-drift” region (see Fig. 4 for the cases with/without a KH instability).

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

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

In this section, we finally show the overall results of our full 1D simulations that include sublimation of ice, the release of silicate dust, and the recycling processes (water vapor as well as silicate dust via the recondensation and sticking onto pebbles outside the snow line – see Fig. 1). We adopt a realistic scale height of silicate dust obtained from Paper I. We also consider the effects of a KH instability for the scale height of pebbles (Sect. 4.1), which was not included in the previous works. We include Drift-BKR and Diff-BKR onto the motions of pebbles and silicate dust. We choose K = 1 for Diff-BKR. We change F_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 solid-to-gas 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 red-colored area indicates runaway pile-up of silicate dust (i.e., no steady-state and an indication of planetesimal formation via direct gravitational collapse). The dependence on Δ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 pile-up occurs in our 1D simulations, it is labeled either by “Peb-RP” for pebbles or by “Sil-RP” for silicate dust. The numbers below the parenthesis in Fig. 5 show the rock fraction in pebbles. For a combination of a small α_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 Drift-BKR (we label these parameters by “ND” as the “no-drift” region in Fig. 5 – see also Sect. 4.2).

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

Fig. 5 Zones of pile-up of pebbles and/or silicate dust around the snow line as a function of Fp/g and αDz(= αDr)∕αacc (case of a variable Δxsubl; 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 pile-up of silicate dust occurs, the fraction is not shown because the system would gravitationally collapse to form 100% rocky planetesimals). The runaway pile-up of pebbles (labeled as “Peb-RP”) or that of silicate dust (labeled as “Sil-RP”) occurs around the snow line, depending on the combinations of Fp/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 Fp/g above whichthe “no-drift (ND)” takes place (labeled by “ND” for the results of 1D simulations in the panels) at the outer boundary (r = 5 au and τs,p = 0.125; Sect. 4.2). Drift-BKR and Diff-BKR onto the motions of pebbles and silicate dust are included with K = 1. Left and right panels: cases of αacc = 10−2 and αacc = 10−3, respectively. In the case of αacc = 10−2, the runaway pile-up of silicate dust is inhibited because the required parameter regime is overlapped by the “no-drift” region (i.e., pebbles do not reach the snow line). In the case of αacc = 10−3, the runaway pile-ups of both pebbles and silicate dust occur.

5.1 Comparison to previous works

In most previous studies of the pile-up 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 pile-up of silicate dust inside the snow line, while a pile-up of pebbles outside the snow line could satisfy a condition for SI to take place for sufficiently large pebble mass flux of F_p/g > 0.8 (Fig. 5). Previous studies (Dra̧żkowska & Alibert 2017; Schoonenberg & Ormel 2017) assumed K = 0 and a runaway pile-up of pebbles was not reported (only a steady-state is reached with ρ_p ∕ρ_g > 1 for F_p/g > 0.8), while it could take place when K≠0 with Drift-BKR (see also Hyodo et al. 2019).

5.2 Pile-up of silicate dust inside the snow line

As discussed in Paper I, silicate dust can pile up just inside the snow line (Fig. 1; see also Estrada et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; 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 solid-to-gas 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⁻² and α_acc = 10⁻³ 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)⁵.

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

5.3 Pile-up of pebbles outside the snow line

In this subsection, we discuss how pebbles pile up just outside the snow line and its dependence on α_acc. Outside the snow line, the pile-up of icy pebbles is locally enhanced due to two distinct recycling processes: the recondensation of water vapor and sticking of silicate dust (Fig. 1; see also Ida & Guillot 2016; Schoonenberg & Ormel 2017; 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 pile-up 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 pile-up 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⁻³). Comparing the maximum ρ_d∕ρ_g and ρ_p∕ρ_g, the pile-up of pebbles is more significant than that of silicate dust when F_p/g is sufficiently large and when the diffusion is an efficient process (i.e., α_Dr∕α_acc ~ 1 cases) because the recycling of the vapor and the dust onto pebbles efficiently enhances the pile-up of pebbles (Hyodo et al. 2019). The pile-up of silicate dust is more prominent when the diffusivity is weak (i.e., α_Dz ∕α_acc ≪ 1 cases).

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

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

The black lines in Fig. 5 indicate where pebbles cannot drift toward the snow line, corresponding to the ND region (see Sect. 4.2; r = 5 au and Ri = 0.5 are used here). In thecase of α_acc = 10⁻², the red-colored region is located above the black line, which implies that parameters in the F_p/g − α_acc − α_Dz space that would yield the runaway pile-up of silicate dust found in Paper I in fact lie in the forbidden ND region. On the other hand, for α_acc = 10⁻³, the pile-up of silicate dust inside the snow line for α_Dz(= α_Dr) < 10⁻⁵ and 0.1 < F_p/g < 0.8 found Paper I is outside of the ND region and confirmed by our simulations.

Generally, runaway pile-ups around the snow line occur over a broader range of parameters for α_acc = 10⁻³ than for α_acc = 10⁻². Runaway pile-up of silicate dust is favored for a small α_Dz∕α_acc < 10⁻³ (labeled by “Sil-RP”), while runaway pile-up of pebbles (labeled by “Peb-RP”) only occurs for a limited range of parameters (large F_p/g and α_Dz ∕α_acc ~ 1).

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

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

Around the snow line, the rock-to-ice mixing ratio can be modified due to sublimation, condensation, and the recycling effects (Hyodo et al. 2019). Our numerical simulations show that pebbles just outside the snow line contain both ice and rock (the numbers below parenthesis in Fig. 5 indicate resulting rock fractions). The results indicate that α_Dr ∕α_acc = 10⁻¹ cases tend to most efficiently increase the rock fraction from the original value (the 1:1 initial rock-ice ratio; Sect. 2) for a fixed F_p/g.

For a higher α_Dr∕α_acc, the scale height of silicate dust becomes larger (Eq. (36)), which leads to a less efficient sticking of silicate dust onto icy pebbles (its efficiency ∝ 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).

Diff-BKR and Drift-BKR also play critical roles as these effects modify the diffusivity and radial motion of solids. As a result of a 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 pile-up of pebbles occurs (labeled by Peb-RP in Fig. 5), the rock fraction becomes ~ 50 wt.%. This is because the runaway pile-up occurs by accreting pebbles that drift from the outer region.

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

5.6 Implication for planet formation

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

The evolution of the disks is characterized by a decrease of the accretion rate, a decrease of viscous heating 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 MRI-inactive “dead zone” near the disk midplane where turbulence and correspondingly diffusion are weak (Fig. 6; Gammie 1996). The radial boundary between the active zone in the outer disk region (α_Dr(= α_Dz)∕α_acc ~ 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 pile-up of icy pebbles leading to the formation of icy planetesimals by GI or SI may take place. As shown in Fig. 6, late in the evolution, the migration of the snow line into the dead inner region would lead to the formation of rocky planetesimals. In-between, in an intermediate phase, the planetesimal formation around the snow line would be suppressed because of a combination of F_p/g and α_Dz ∕α_acc that are insufficient to lead to a pile-up of solids in the disk midplane.

As the observed structure and theoretically modeled evolutionary paths of protoplanetary disks are diverse, following a distinct evolutionary path in the α_acc–α_Dr–α_Dz–F_p/g space as well 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 MRI-inactive 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 Fp/g value that becomes smaller with time, while the snow line moves inward. Our results of planetesimal formation around the snow line imply that the icy pebble pile-up (i.e., icy planetesimals formation by GI and/or SI; panel A) would be favored in the outer region (α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 Fp/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 pile-up of solids around the snow line, we incorporated a realistic H_d (obtained from Paper I) and the KH-instability-considered H_p (Sect. 4.1), for the first time, into a 1D diffusion-advection code that includes the back-reactions to radial drift and diffusion of icy pebbles and silicate dust. The code takes into account ice sublimation, the release of silicate dust, and their recycling through the recondensation and sticking onto pebbles. We studied a much wider range of disk parameters (in the F_p/g − α_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 Drift-BKR in a self-induced manner (the “no-drift” region; Sect. 4.2).

Third, our 1D simulations showed that α_acc = 10⁻³ case is more favorable for the pile-up of solids around the snow line than that for α_acc = 10⁻². The “no-drift” regime entirely covers parameter space of the runaway pile-up of silicate dust when α_acc = 10⁻² case (left panel in Fig. 5), preventing a runaway pile-up of silicate dust inside the snow line. In contrast, in the α_acc = 10⁻³ case, large F_p/g values are allowed and could lead to a runaway pile-up of silicate dust inside the snow line (right panel in Fig. 5). In the case of α_acc = 10⁻³, both silicate dust and pebbles could experience runaway pile-ups inside and outside the snow line, forming rocky and icy planetesimals, respectively.

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

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

Acknowledgements

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

Appendix A Dependence of dust pile-up 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 1D-radial accretion disk (e.g., Ida & Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019). The midplane density of silicate dust inside the snow line is calculated by , which is also used to describe the strength of the back-reaction of the gas onto their motion (e.g., Eq. (18)). Therefore, the description of H_d is critical for the back-reaction and to evaluate the pile-up of dust in the midplane. In this section, we demonstrate how different models of H_d affect the dust pile-up inside the snow line.

Fig. A.1 Scale height of silicate dust in a unit of that of the gas (top panels) and midplane dust-to-gas ratio (bottom panels) for different models of Hd (Case of Fp/g = 0.3 and τs,p = 0.1 at 3 au). Here, Drift-BKR and Diff-BKR for pebbles are neglected, whereas Drift-BKR and Diff-BKR (K = 0 or 1) for silicate dust are included. Left two panels: cases of α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 pile-up 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), pile-ups occur in a runaway fashion for the left two panels, and the time-evolutions are shown for t = 1 × 104 yr, t = 5 × 104 yr, and t = 1.5 × 105 yr. The other cases reach a steady-state within t = 5 × 105 yr. Here, the recycling of water vapor and silicate dust onto pebbles are not included to focus on the different models of Hd.

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 pile-ups of silicate dust, including a more realistic model used in this work (Eq. (36) which is originally derived in Paper I).

The model of Schoonenberg & Ormel (2017), H_d = H_g, underestimated pile-up of silicate dust when α_Dz ≪ α_acc, while it correctly evaluated when α_Dz ≃ α_acc. Ida & Guillot (2016), H_d = H_p,snow, overestimated the pile-up of silicate dust regardless of the value of α_Dz. Hyodo et al. (2019) overestimated pile-up of silicate dust near the snow line.

Appendix B Particle size and Stokes number

The Stokes number of a particle τ_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⁻³ 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−8M⊙ yr−1 and are used.

(B.3)

where σ_mol = 2.0 × 10⁻¹⁵ cm² 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⁻¹ 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 advection-dominated regime) as both radial and vertical mixings are ineffective, while large α_Dr (= α_Dz)∕α_acc (i.e., the diffusion-dominated regime) case is rarely affected because the mixing effects “erase” the information of the initial distribution of silicate dust. Thus, the parameter space for runaway pile-up of silicate dust is reduced for small α_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).

Fig. C.1 Same as Fig. 5, but for the case of a constant Δxsubl = 0.1Hg.

References

All Figures

Fig. 1 Summary of solid pile-up around the snow line. Icy pebbles, which uniformly contain micron-sized silicate dust formed at the outer disk, drift inward due to the gas drag. Pebbles sublimate through the passage of the snow line (Sect. 3) and silicate dust is released together with the water vapor. Silicate dust is well coupled to the gas and the radial drift velocity significantly drops from that of pebbles. This causes so-called “traffic-jam” effect and the dust piles up inside the snow line. A fraction of water vapor and silicate dust diffuses from inside to outside the snow line recyclesonto pebbles via the recondensation and sticking, causing a local pile-up of pebbles just outside the snow line. Midplane concentrations (i.e., midplane solid-to-gas ratio) of pebbles and silicate dust strongly depend on their scale heights. Near the snow line, the scale height of silicate dust is the minimum as the dust is released from pebbles whose scale height is much smaller than that of the gas. This causes the maximum midplane concentration of silicate dust just inside the snow line. Significant midplane concentrations of silicate dust and/or icy pebbles would cause gravitational instability and/or streaming instability, forming rocky and/or icy planetesimals, respectively (Sect. 5).

In the text

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

In the text

Fig. 3 Radial location of the snow line (top panel) and sublimation width of drifting pebbles Δxsubl (bottom panel) as a function of αDr∕αacc. Here, fd/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 Fp/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 light-blue line in the bottom panel shows the analytical Δxsubl in the case of the advection-dominated regime (Eq. (54) with rsnow = 2.4 au and β = 0.5). The light-green line in the bottom panel shows the analytical Δxsubl in the case of the diffusion-dominated regime (Eq. (68) with τs,p = 0.06, rp = 2.3 cm, rsnow = 2.4 au). The gray dashed line in the bottom panel shows the analytically derived boundary betweenthe advection- and diffusion-dominated 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 pebble-to-gas density ratio (Z = ρp∕ρg) as a functionof αDz and Fp/g (at r = 5 au and τs,p = 0.125). Top panels: cases where the back-reaction onto the gas motion is neglected (i.e., vg = vg,ν). Bottom panels: cases where the back-reaction to the gas motion is included, following Eq. (16). The color contours are obtained by directly solving Eq. (86). Left and right two panels: cases where αacc = 10−3 and αacc = 10−2 with/without a KH instability (Ri = 0.5), respectively. The red-colored regions indicate where a steady-state solution is not found (i.e., ρp ∕ρg > 100 and it keeps increasing; the “no-drift (ND)” region). The white (Eq. (91)) and black (Eq. (95) with Ri = 0.5) dashed lines are analytically derived critical Fp/g,crit above which no steady-state are found (the “no-drift (ND)” region). Analytical predictions (dashed lines) well predict critical boundaries of the ND region. Including a KH instability reduces the ND region as the minimum scale height of pebbles is regulated by a KH instability. Due to the back-reaction to the gas motion that slows down the radial velocity of the gas, the midplane concentration of pebbles as well as the ND region (red-colored region) are shifted toward slightly larger Fp/g values.

In the text

Fig. 5 Zones of pile-up of pebbles and/or silicate dust around the snow line as a function of Fp/g and αDz(= αDr)∕αacc (case of a variable Δxsubl; 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 pile-up of silicate dust occurs, the fraction is not shown because the system would gravitationally collapse to form 100% rocky planetesimals). The runaway pile-up of pebbles (labeled as “Peb-RP”) or that of silicate dust (labeled as “Sil-RP”) occurs around the snow line, depending on the combinations of Fp/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 Fp/g above whichthe “no-drift (ND)” takes place (labeled by “ND” for the results of 1D simulations in the panels) at the outer boundary (r = 5 au and τs,p = 0.125; Sect. 4.2). Drift-BKR and Diff-BKR onto the motions of pebbles and silicate dust are included with K = 1. Left and right panels: cases of αacc = 10−2 and αacc = 10−3, respectively. In the case of αacc = 10−2, the runaway pile-up of silicate dust is inhibited because the required parameter regime is overlapped by the “no-drift” region (i.e., pebbles do not reach the snow line). In the case of αacc = 10−3, the runaway pile-ups 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 MRI-inactive 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 Fp/g value that becomes smaller with time, while the snow line moves inward. Our results of planetesimal formation around the snow line imply that the icy pebble pile-up (i.e., icy planetesimals formation by GI and/or SI; panel A) would be favored in the outer region (α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 Fp/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 dust-to-gas ratio (bottom panels) for different models of Hd (Case of Fp/g = 0.3 and τs,p = 0.1 at 3 au). Here, Drift-BKR and Diff-BKR for pebbles are neglected, whereas Drift-BKR and Diff-BKR (K = 0 or 1) for silicate dust are included. Left two panels: cases of α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 pile-up 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), pile-ups occur in a runaway fashion for the left two panels, and the time-evolutions are shown for t = 1 × 104 yr, t = 5 × 104 yr, and t = 1.5 × 105 yr. The other cases reach a steady-state within t = 5 × 105 yr. Here, the recycling of water vapor and silicate dust onto pebbles are not included to focus on the different models of Hd.

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−8M⊙ yr−1 and are used.

In the text

	Fig. C.1 Same as Fig. 5, but for the case of a constant Δxsubl = 0.1Hg.
In the text