© ESO, 2017

1. Introduction

Planets form in gaseous disks around young stars. Initially, micron-sized dust grains are present in such a disk. It is generally thought that planets form from planetesimals of the order of a kilometer in size, which are large enough for gravity to play an important role (Safronov 1969; Pollack et al. 1996; Benz 2000). However, the first step in the planet formation process, the coagulation of dust grains into planetesimals, is still an unsolved problem. Although micron-sized grains can quite easily coagulate to mm-sized particles (Dominik & Tielens 1997; Ormel et al. 2007), typical relative velocities quickly become so large that particles fragment or bounce rather than stick (Blum & Wurm 2000; Zsom et al. 2010). The typical size at which coagulation is no longer feasible is referred to as the fragmentation or bouncing barrier. Observations at radio wavelengths indeed infer a large reservoir of these mm-to-cm size particles (pebbles) (e.g., Testi et al. 2014).

Even if it were possible to coagulate uninterruptedly from micron-sizes to large bodies, the growth process would face another challenge: that of radial drift. As soon as particles reach sizes at which they aerodynamically start to decouple from the gas, they lose angular momentum and spiral towards the star (Whipple 1972; Weidenschilling 1977). Radial drift is a fast process and hence poses a challenge to planetesimal formation theories, as planetesimals should form before the solids are lost to the star.

A promising mechanism to quickly form planetesimals directly from pebbles is streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Johansen & Youdin 2007). Streaming instability leads to clumping of pebbles that subsequently become gravitationally unstable and form planetesimals. In order for streaming instabilities to occur, the metallicity of the disk has to be locally enhanced above the typical value of 1% (Johansen et al. 2009; Bai & Stone 2010; Dra¸żkowska & Dullemond 2014). However, if one considers a smooth disk (i.e., without any radial pressure bumps or special locations), solids are flushed to the star in an inside-out manner due to dust growth and radial drift (e.g., Birnstiel et al. 2010, 2012; Krijt et al. 2016b; Sato et al. 2016)¹. Radial drift lowers the metallicity of the disk and – at first sight – works against creating conditions favorable to streaming instability.

A special location that can affect the distribution of solids in the disk is the snowline. The snowline is the radius from the central star interior to which all water is in the form of vapor. Recent observations of structures in protoplanetary disks (for example HL Tau; ALMA Partnership et al. 2015) have been proposed to be associated with icelines (not necessarily water), e.g., Banzatti et al. (2015), Okuzumi et al. (2016), Nomura et al. (2016).

The viability of a pile-up of solids at the snowline has been a topic of interest for considerable time. Several mechanisms have been proposed, such as the creation of a pressure maximum (e.g., Brauer et al. 2008) or a deadzone (e.g., Kretke & Lin 2007). Another pile-up mechanism is water diffusion and condensation: water vapor in the inner disk can diffuse outward across the snowline and condense, thereby increasing the solids reservoir just outside the snowline. Works that have investigated the role of water diffusion and condensation on the growth and distribution of solids include Stevenson & Lunine (1988), Cuzzi & Zahnle (2004), Ciesla & Cuzzi (2006), Ros & Johansen (2013), Estrada et al. (2016), Armitage et al. (2016) and, recently, Krijt et al. (2016a), which focuses on the atmospheric snowline. In this paper, we investigate the role of the radial snowline to enhance the solids-to-gas ratio.

The works most relevant to our problem are Stevenson & Lunine (1988) and Ros & Johansen (2013). Stevenson & Lunine (1988) investigated the role of outward diffusion of water vapor across the snowline. In their work, a “cold-finger” effect was considered: water vapor from the inner disk condenses on (pre-existing) solids beyond the snowline. In this way, they found a steep increase in the solid density beyond the snowline and hypothesized that this could accelerate the formation of Jupiter. However, dynamical aspects of the disk, such as radial drift of solids and accretion of gas, were not taken into account; except for the vapor, their model is essentially static. More recently, Ros & Johansen (2013) also simulated pebble growth by condensation, adopting a Monte Carlo method. In their simulation, the purely icy particles grow by condensation, until the point that they are “lost” by evaporation (modeled as an instantaneous process). Because their work does not feature silicate “seeds” that are immune to evaporation, the remaining ice particles can only grow in size. Consequently, a steady-state is never reached, also not because the simulation domain does not connect to the rest of the disk: there is no inflow of small pebbles or removal of vapor-rich gas.

In this work we consider a model that emphasizes the dynamic nature of the disk, characterized by the accretion rates of pebbles and gas. Our model of the iceline therefore includes inflow of pebbles (at the outer boundary) and outflow of vapor (at the inner boundary) – a design that results in a steady-state. Our goal is to investigate whether-or-not this more dynamic design will result in a solids enhancement that could trigger streaming instabilities. In particular, we will investigate how the solids-to-gas ratio will depend on the properties of the disk: the accretion rates (of solids and gas), the turbulence strength, and the aerodynamical properties of the drifting pebbles. We will also investigate which side of the snowline the solids-to-gas ratio is boosted most: interior or exterior.

We use a characteristic particle method, assuming a single pebble size at each time and location in the disk (Ormel 2014; Sato et al. 2016). We numerically solve the partial differential equations that govern the transport and condensation/evaporation processes for solids and vapor. We consider two model designs for the composition of pebbles: “single seed” (all silicates in one core) and “many seeds” (many tiny silicates grains). Our model also accounts for feedback of solids on the gas (a.k.a. collective effects) and feedback of water vapor on the scale height of the disk. A semi-analytic model is presented that approximates the numerical results, enabling us to carry out parameter searches.

We start out with a description of our model in Sect. 2. In Sect. 4 we discuss numerical results for a fiducial set of model parameters in detail. We first describe the time-dependent numerical results in Sect. 4.1, after which we focus on the steady-state solution in Sect. 4.2. In Sect. 5 we summarize the semi-analytical approximate model (which is described in more detail in Appendix C) and present our results concerning streaming instability conditions. We discuss our results in Sect. 6 and present our main conclusions in Sect. 7. A list of frequently used symbols can be found in Table A.1.

2. Model

2.1. Disk model

Throughout this paper, we use the α-accretion model (Shakura & Sunyaev 1973) to model disks. We fix the mass of the central star M_⋆ = M_⊙. We fix the temperature profile T(r) to: $\mathit{T}\mathrm{=}\mathrm{150}{\left(\frac{\mathit{r}}{\mathrm{3}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{K}$(1)where r is the radial distance from the star. This temperature power-law profile corresponds to a passively-irradiated disk and a solar-mass star (Kenyon & Hartmann 1987; Armitage et al. 2016); for simplicity, we neglect viscous heating. The isothermal sound speed c_s is related to the temperature as: ${\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }}}$(2)with k_B the Boltzmann constant and μ the molecular weight of the gas, which is 2.34 times the proton mass for a typical solar metallicity gas. The disk scale height H_gas is given by: ${\mathit{H}}_{\mathrm{gas}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}}{\mathrm{\Omega }}\mathrm{\approx }\mathrm{0.033}{\left(\frac{\mathit{r}}{\mathrm{1}\mathrm{au}}\right)}^{\mathrm{5}\mathit{/}\mathrm{4}}\mathrm{au}$(3)where Ω is the Keplerian orbital frequency. The gas accretes to the central star at a rate Ṁ_gas. Given α and Ṁ_gas, the steady-state gas surface density is (Lynden-Bell & Pringle 1974): ${\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{gas}}{\mathrm{2}\mathit{\pi r}{\mathit{v}}_{\mathrm{gas}}}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{gas}}{\mathrm{3}\mathit{\pi \nu }}$(4)where v_gas = 3ν/ 2r is the radial speed of the accreted gas, with ν the turbulent viscosity, which is given by: $\mathit{\nu }\mathrm{=}\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}{\mathit{H}}_{\mathrm{gas}}$(5)where α is the dimensionless turbulence parameter. As a result of our choices for the temperature and scale height profiles, Σ_gas(r) is a power-law with index − 1. We take the turbulent gas diffusivity D_gas equal to the turbulent viscosity ν.

2.2. Radial motion of ice, silicates and vapor

In this study, we are interested in the pebble-sized fraction of the solids – we do not consider the smaller-sized dust. Pebbles are characterized by the fact that they are aerodynamically partly decoupled from the gas disk. The gas disk is partly pressure-supported, and therefore rotates at a sub-Keplerian velocity. The pebbles do not feel the pressure gradient, and therefore tend to move toward Keplerian orbital velocities. This means that pebbles are moving through a more slowly rotating gas disk, thereby losing angular momentum through gas drag. The loss of angular momentum results in pebbles spiralling (“drifting”) towards the central star. The extent to which pebbles drift depends on the stopping time t_stop, which measures how strongly the pebbles are coupled to the gas – it is the time after which any initial momentum is lost due to gas friction.

Concerning the stopping times of pebbles we can distinguish between two regimes: in the Stokes regime, the particle size is larger than the mean-free path of the gas molecules l_mfp, and the stopping time is calculated in a fluid description; in the Epstein regime, the particle size is smaller than the mean-free path of the gas molecules, and a particle description is needed instead. The mean-free path l_mfp is given by: ${\mathit{l}}_{\mathrm{mfp}}\mathrm{=}\frac{\mathit{\mu }}{\sqrt{\mathrm{2}}{\mathit{\rho }}_{\mathrm{gas}}{\mathit{\sigma }}_{\mathrm{mol}}}$(6)where σ_mol is the molecular collision cross-section, and ρ_gas is the gas density. We take σ_mol = 2 × 10^-15cm² as the collisional cross-section of molecular hydrogen (Chapman & Cowling 1970). The stopping time is given by: ${\mathit{t}}_{\mathrm{stop}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{(}\mathrm{Epstein:}{\mathit{s}}_{\mathrm{p}}\mathit{<}\frac{\mathrm{9}}{\mathrm{4}}{\mathit{l}}_{\mathrm{mfp}}\mathrm{)}\\ & \mathrm{(}\mathrm{Stokes:}{\mathit{s}}_{\mathrm{p}}\mathit{>}\frac{\mathrm{9}}{\mathrm{4}}{\mathit{l}}_{\mathrm{mfp}}\mathrm{)}\end{array}$(7)where s_p is the particle radius, ρ_•,p is the particle internal density, and v_th is the thermal velocity of the gas molecules, defined as ${\mathit{v}}_{\mathrm{th}}\mathrm{=}\sqrt{\mathrm{8}\mathit{/}\mathit{\pi }}{\mathit{c}}_{\mathrm{s}}$. The dimensionless stopping time (sometimes referred to as the Stokes number) is τ_S = t_stopΩ.

The radial drift velocity of pebbles v_peb is given by (Weidenschilling 1977; Nakagawa et al. 1986): ${\mathit{v}}_{\mathrm{peb}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{gas}}\mathrm{+}\mathrm{2}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}{\mathit{\tau }}_{\mathrm{S}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}_{\mathrm{S}}^{\mathrm{2}}}$(8)where ηv_K is the magnitude of the azimuthal motion of the gas disk below the Keplerian velocity v_K: $\mathit{\eta }{\mathit{v}}_{\mathrm{K}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{K}}}\frac{\mathit{\partial }\log \mathit{P}}{\mathit{\partial }\log \mathit{r}}$(9)with P the gas pressure. Under normal disk conditions, the solids-to-gas ratio is about 1 in 100, and it is reasonable to neglect the back-reaction of the solids on the gas. In this paper we are interested in situations where the solids-to-gas ratio is boosted, however, and therefore we study the effects of the back-reaction of solids on the gas, or “collective effects”, as well. Collective effects have been calculated by Nakagawa et al. (1986). However, their formula was derived in the inviscid limit. In Appendix B we have calculated the drift velocity, accounting for both the backreaction and viscous forces. When including the back-reaction, we replace Eq. (8) by Eq. (B.7).

Water vapor (denoted by the subscript Z) mixes with the gas and gets accreted to the central star at the same velocity as the gas.

2.3. Evaporation and condensation

The rate of change in mass (dm_p/ dt) of a spherical icy pebble is given by (Lichtenegger & Komle 1991; Ciesla & Cuzzi 2006; Ros & Johansen 2013): $\frac{\mathrm{d}{\mathit{m}}_{\mathrm{p}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{s}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{v}}_{\mathrm{th}\mathit{,Z}}{\mathit{\rho }}_{\mathit{Z}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{P}}_{\mathrm{eq}}}{{\mathit{P}}_{\mathit{Z}}}\right)$(10)where s_p is the radius of the pebble, v_th,Z is the thermal velocity of vapor particles, ρ_Z is the vapor density, P_Z is the vapor pressure, and P_eq is the saturated or equilibrium pressure, which is given by the Clausius-Clapeyron equation: ${\mathit{P}}_{\mathrm{eq}}\mathrm{=}{\mathit{P}}_{\mathrm{eq}\mathit{,}\mathrm{0}}{\mathrm{e}}^{\mathrm{-}{\mathit{T}}_{\mathit{a}}\mathit{/}\mathit{T}}$(11)where T_a and P_eq,0 are constants depending on the species. For water, T_a = 6062 K and P_eq,0 = 1.14 × 10¹³gcm^-1s^-2 (Lichtenegger & Komle 1991).

Assuming spherical pebbles and one typical particle mass m_p at each location and time in the disk (i.e., the particle size distribution is narrow), and using the ideal gas law to write P_Z in terms of the (midplane) vapor density ρ_Z, we can rewrite Eq. (10) into a change in the ice surface density Σ_ice due to condensation and evaporation: $\mathrm{\Sigma ̇}\mathrm{ice}\mathit{,}\mathrm{C}\mathit{/}\mathrm{E}\mathrm{=}{\mathit{R}}_{\mathrm{c}}{\mathrm{\Sigma }}_{\mathrm{ice}}{\mathrm{\Sigma }}_{\mathit{Z}}\mathrm{-}{\mathit{R}}_{\mathrm{e}}{\mathrm{\Sigma }}_{\mathrm{ice}}$(12)where the dot denotes the time derivative. In Eq. (12) the vapor surface density Σ_Z is related to the midplane vapor density as ${\mathit{\rho }}_{\mathit{Z}}\mathrm{=}{\mathrm{\Sigma }}_{\mathit{Z}}\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}$ and R_c and R_e are the condensation rate and evaporation rate, respectively, defined as: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{c}}& \mathrm{=}& \mathrm{8}\sqrt{\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{\mu }}_{\mathit{Z}}}}\frac{{\mathit{s}}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{p}}{\mathit{H}}_{\mathrm{gas}}}\\ {\mathit{R}}_{\mathrm{e}}& \mathrm{=}& \mathrm{8}\sqrt{\mathrm{2}\mathit{\pi }}\frac{{\mathit{s}}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{p}}}\sqrt{\frac{{\mathit{\mu }}_{\mathit{Z}}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}}{\mathit{P}}_{\mathrm{eq}}\end{array}$with k_B the Boltzmann constant, and μ_Z = 18m_H the mean molecular weight of water vapor, where m_H is the proton mass. Evaporating ice turns into vapor, and condensing vapor turns into ice, and therefore we find for the vapor source terms: $\mathrm{\Sigma ̇}\mathit{Z,}\mathrm{C}\mathit{/}\mathrm{E}\mathrm{=}\mathrm{-}\mathrm{\Sigma ̇}\mathrm{ice}\mathit{,}\mathrm{C}\mathit{/}\mathrm{E}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathrm{c}}{\mathrm{\Sigma }}_{\mathrm{ice}}{\mathrm{\Sigma }}_{\mathit{Z}}\mathrm{+}{\mathit{R}}_{\mathrm{e}}{\mathrm{\Sigma }}_{\mathrm{ice}}\mathit{.}$(15)

2.4. Internal structure and composition of pebbles

We assume that in the outer disk, half of the mass in pebbles consists of water ice, and the other half consists of silicates (Lodders 2003; Morbidelli et al. 2015). That is, the dust-fraction in pebbles in the outer disk ζ = 0.5. Adopting ρ_•,ice = 1gcm^-3 as the density of a pure ice particle and ρ_•,sil = 3gcm^-3 as the density of a pure silicate particle, we find that in the outer disk, the internal density of a pebble is ρ_•,p = 1.5gcm^-3. Generally, the internal density of an ice/silicate pebble is given by: ${\mathit{\rho }}_{\mathrm{•}\mathit{,}\mathrm{p}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{•}\mathit{,}\mathrm{ice}}{\mathit{\rho }}_{\mathrm{•}\mathit{,}\mathrm{sil}}\mathrm{(}{\mathit{m}}_{\mathrm{sil}}\mathrm{+}{\mathit{m}}_{\mathrm{ice}}\mathrm{)}}{{\mathit{m}}_{\mathrm{sil}}{\mathit{\rho }}_{\mathrm{•}\mathit{,}\mathrm{ice}}\mathrm{+}{\mathit{m}}_{\mathrm{ice}}{\mathit{\rho }}_{\mathrm{•}\mathit{,}\mathrm{sil}}}$(16)where m_ice is the mass in ice within the pebble and m_sil is the mass in silicate within the pebble. Equation (10) is valid for a pure ice particle, while we use it for icy pebbles with (a) silicate core(s). We checked that using m_ice instead of m_p in the condensation and evaporation rates does not change the results significantly.

We consider two possibilities for the internal structure of icy pebbles (Fig. 1).

• Single-seed model. In this model we view pebbles in the outer diskas balls of ice with a single core made of silicates. When the ice insuch a pebble evaporates, the silicate core is left behind. In thisframework, the number of pebbles does not change whencrossing the snowline.

• Many-seeds model. In this model we view icy pebbles as lots of small silicate particles that are “glued” together by ice. The small silicate particles are released by evaporation, and free silicates can stick onto icy pebbles. We do not “resolve” the internal structure of the pebbles, and only trace the total amount of silicates “locked” into icy pebbles. Following Saito & Sirono (2011), we assume that the released silicate particles are of micron size. Because they are so small, the silicates do not undergo any aerodynamic drift and follow the radial motion of the gas.

2.5. Transport equations

In this subsection we provide the systems of transport equations corresponding to both pebble composition models presented in Sect. 2.4. Our model is 1D; we integrate over the vertical dimension of the disk and follow column (surface) densities. Midplane densities are calculated assuming a Gaussian distribution in height with scaleheights H_gas (for the gas and the vapor) and H_peb (for the pebbles), respectively.

2.5.1. Single-seed pebble model

In the single-seed model, we follow the vapor surface density Σ_Z, the ice surface density Σ_ice, and the number density of particles N_p. The equations governing the time evolution of these quantities are: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{ice}}}{\mathit{\partial t}}& \mathrm{+}& \mathrm{\nabla }\mathrm{·}{\mathrm{ℳ}}_{\mathrm{ice}}\mathrm{=}\mathrm{\Sigma ̇}\mathrm{ice}\mathit{,}\mathrm{C}\mathit{/}\mathrm{E}\\ \frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathit{Z}}}{\mathit{\partial t}}& \mathrm{+}& \mathrm{\nabla }\mathrm{·}{\mathrm{ℳ}}_{\mathit{Z}}\mathrm{=}\mathrm{\Sigma ̇}\mathit{Z,}\mathrm{C}\mathit{/}\mathrm{E}\\ \frac{\mathit{\partial }{\mathit{N}}_{\mathrm{p}}}{\mathit{\partial t}}& \mathrm{+}& \mathrm{\nabla }\mathrm{·}{\mathrm{𝒩}}_{\mathrm{p}}\mathrm{=}\mathrm{0}\end{array}$where the source terms and are defined in Eqs. (12) and (15), ℳ_ice and ℳ_Z are the ice and vapor mass flux, respectively, and is the number flux of solid particles. These fluxes are given by:

$\begin{array}{ccc}{\mathrm{ℳ}}_{\mathrm{ice}}& \mathrm{=}& \mathrm{-}{\mathrm{\Sigma }}_{\mathrm{ice}}{\mathit{v}}_{\mathrm{peb}}\mathrm{-}{\mathit{D}}_{\mathrm{p}}{\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{\nabla }\frac{{\mathrm{\Sigma }}_{\mathrm{ice}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\\ {\mathrm{ℳ}}_{\mathit{Z}}& \mathrm{=}& \mathrm{-}{\mathrm{\Sigma }}_{\mathit{Z}}{\mathit{v}}_{\mathrm{gas}}\mathrm{-}{\mathit{D}}_{\mathrm{gas}}{\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{\nabla }\frac{{\mathrm{\Sigma }}_{\mathit{Z}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\\ {\mathrm{𝒩}}_{\mathrm{p}}& \mathrm{=}& \mathrm{-}{\mathit{N}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{peb}}\mathrm{-}{\mathit{D}}_{\mathrm{p}}{\mathit{N}}_{\mathrm{gas}}\mathrm{\nabla }\frac{{\mathit{N}}_{\mathrm{p}}}{{\mathit{N}}_{\mathrm{gas}}}\\ & & \end{array}$where v_peb is the drift velocity of the pebbles given by Eq. (8) or by Eq. (B.7) when including collective effects, N_gas is the number of gas particles, and D_p is the particle diffusivity, which is related to D_gas through the Schmidt number Sc (Youdin & Lithwick 2007): ${\mathit{D}}_{\mathrm{p}}\mathrm{=}{\mathit{D}}_{\mathrm{gas}}\mathrm{\times }\mathrm{Sc}\mathrm{=}\frac{{\mathit{D}}_{\mathrm{gas}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}_{\mathrm{S}}^{\mathrm{2}}}\mathrm{·}$(23)Since in our models the stopping times are much smaller than unity, we set D_p equal to D_gas.

Note that a negative flux means that material is transported inwards, towards the star. Under the single-size approximation, the typical pebble mass m_p is given by: ${\mathit{m}}_{\mathrm{p}}\mathrm{=}{\mathit{m}}_{\mathrm{core}}\mathrm{+}\frac{{\mathrm{\Sigma }}_{\mathrm{ice}}}{{\mathit{N}}_{\mathrm{p}}}$(24)with m_core the mass of the bare silicate cores of the pebbles.

2.5.2. Many-seeds pebble model

We implement the many-seeds model as follows. We now consider two additional populations of solids: dirt and free silicates. Dirt refers to the silicates that are captured in icy pebbles, whereas free silicates are small silicate grains without any ice coating. The transport equations for Σ_ice, Σ_Z, and N_p are the same as in the single-seed pebble model (Eqs. (17)–(19)). The transport equations for the dirt surface density Σ_dirt and free silicates surface density Σ_sil are given by: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{dirt}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}{\mathrm{ℳ}}_{\mathrm{dirt}}& \mathrm{=}& \mathrm{-}\max \left[\mathrm{0}\mathit{,}\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathrm{-}{\mathit{R}}_{\mathrm{c}}{\mathrm{\Sigma }}_{\mathit{Z}}\mathrm{)}\right]{\mathrm{\Sigma }}_{\mathrm{dirt}}\mathrm{+}{\mathit{R}}_{\mathrm{s}}{\mathrm{\Sigma }}_{\mathrm{sil}}\\ \frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{sil}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}{\mathrm{ℳ}}_{\mathrm{sil}}& \mathrm{=}& \max \left[\mathrm{0}\mathit{,}\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathrm{-}{\mathit{R}}_{\mathrm{c}}{\mathrm{\Sigma }}_{\mathit{Z}}\mathrm{)}\right]{\mathrm{\Sigma }}_{\mathrm{dirt}}\mathrm{-}{\mathit{R}}_{\mathrm{s}}{\mathrm{\Sigma }}_{\mathrm{sil}}\\ & & \end{array}$where ℳ_dirt and ℳ_sil are the dirt and free silicate mass fluxes, respectively. Free silicates are released from the pebbles at a rate proportional to the decrease of ice surface density due to evaporation and to the amount of dirt that can evaporate from within the pebbles: (R_e − R_cΣ_vap)Σ_dirt, but only when there is net evaporation of ice. Silicate grains stick to icy pebbles at a collision/sticking rate R_s: ${\mathit{R}}_{\mathrm{s}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{ice}}\mathrm{+}{\mathrm{\Sigma }}_{\mathrm{dirt}}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{sil}}{\mathit{m}}_{\mathrm{p}}}\mathrm{·}\mathrm{\Delta }{\mathit{v}}_{\mathrm{sil}\mathit{,}\mathrm{peb}}\mathrm{·}\mathit{\pi }{\mathit{s}}_{\mathrm{p}}^{\mathrm{2}}$(27)where Δv_sil,peb is the relative velocity between the pebbles and the silicates and H_sil is the scale height of the free silicates. Since the inward drift velocity of the pebbles is much larger than the outward diffusion velocity of the silicates, we take Δv_sil,peb ≈ v_peb².

The dirt and silicate mass fluxes, ℳ_dirt and ℳ_sil, are given by:

$\begin{array}{ccc}& & {\mathrm{ℳ}}_{\mathrm{dirt}}\mathrm{=}\mathrm{-}{\mathrm{\Sigma }}_{\mathrm{dirt}}{\mathit{v}}_{\mathrm{peb}}\mathrm{-}{\mathit{D}}_{\mathrm{p}}{\mathrm{\Sigma }}_{\mathrm{gas}}\mathrm{\nabla }\frac{{\mathrm{\Sigma }}_{\mathrm{dirt}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\\ & & \end{array}$where v_peb is again the radial drift velocity of the pebbles, and silicates are advected with the gas.

2.6. Effects of variable μ

When icy pebbles evaporate, the gas becomes enriched in water vapor. This means that the mean molecular weight μ of the gas increases: $\frac{\mathrm{1}}{\mathit{\mu }}\mathrm{=}\frac{{\mathit{f}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}\mathrm{+}\frac{{\mathit{f}}_{\mathrm{gas}}}{{\mathit{\mu }}_{\mathrm{gas}}}$(30)where f_H2O is the mass fraction of water molecules, μ_H2O = 18m_H is the mean molecular weight of water, and f_gas is the mass fraction of gas molecules (without the vapor contribution), with μ_gas = 2.34m_H.

Since the sound speed c_s is proportional to $\sqrt{\mathrm{1}\mathit{/}\mathit{\mu }}$, the sound speed decreases as more water vapor is added. The disk scale height H_gas is proportional to c_s and will therefore decrease as well: the disk becomes more vertically compact with increasing water vapor. An increase in μ across the water evaporation front has an effect on the pebble drift velocity v_peb (Eq. (8)) in this region, since v_peb is proportional to 1 /μ through its dependence on c_s, and sensitive to variations in the gas pressure P through ∂log P/∂log r. The total gas pressure (hydrogen/helium plus water) is given by: $\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{gas}}\mathrm{+}{\mathit{P}}_{\mathit{Z}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{gas}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{/}\mathrm{He}}}\mathrm{+}\frac{{\mathrm{\Sigma }}_{\mathit{Z}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}\mathrm{·}$(31)We take the viscosity ν to be independent of variations in μ. Therefore, Σ_gas = Ṁ_gas/ 3πν is independent of μ as well. However, the H₂O-enriched gas is characterized by a reduced scaleheight, since ${\mathit{H}}_{\mathrm{gas}}\mathrm{\propto }\sqrt{\mathrm{1}\mathit{/}\mathit{\mu }}$, increasing the midplane gas pressure P_gas as ${\mathit{P}}_{\mathrm{gas}}\mathrm{\propto }\sqrt{\mathit{\mu }}$. Evaporation hence increases the pressure gradient and the headwind of the gas through Eq. (9): pebbles tend to drift faster. Finally, the drift velocity v_peb also depends on the stopping time τ_S, which in turn depends on the gas density and the thermal velocity – both quantities that depend on μ. We take these dependencies into account as well.

We now define two models that we study in the following sections. In the simple model, we do not include collective effects and the variability of the mean molecular weight μ. In the simple model, the sound speed and scale height are independent of variations in μ, ηv_K (Eq. (9)) is constant and the midplane gas density is μ-independent: ${\mathit{\rho }}_{\mathrm{gas}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{gas}}\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}$. In the complete model, we do take collective effects and variations of μ into account, as described above.

2.7. Input parameters

The four key input parameters that can be varied in our model are:

• The gas accretion rate Ṁ_gas.

• The ratio between the pebble accretion rate and the gas accretion rate, denoted by the symbol ℱ_{s / g}, and defined as: ${\mathrm{ℱ}}_{\mathrm{s}\mathit{/}\mathrm{g}}\mathrm{\equiv }\mathit{Ṁ}\mathrm{peb}\mathit{/}\mathit{Ṁ}\mathrm{gas}\mathit{.}$(32)

• The dimensionless turbulence parameter α.

• The initial stopping time of pebbles just outside the snowline (at 3 au) denoted by the symbol τ₃. This parameter is simply a proxy for the mass of the incoming pebbles.

Together, these input parameters define the “zero-model” of the disk: the gas and solids distribution in the case of no condensation and evaporation. The gas surface density profile Σ_gas is defined through Eq. (4), where the turbulent viscosity ν is determined by α through Eq. (5). With Σ_gas and the stopping time normalization τ₃ at 3 au, we can find τ(r) at any other radius r. We can then also determine the pebbles surface density Σ_peb for the zero-model: it is given by ℱ_{s / g}Ṁ_gas/ 2πrv_peb(r), where the pebble velocity v_peb(r) depends on τ(r) and is given by Eq. (8) (without backreaction) or Eq. (B.7) (including backreaction).

2.8. Model assumptions

In constructing our model we employed several assumptions, which we discuss below.

• We do not include a particle size distribution, but rather considerone typical pebble size at each time and location in the disk. In themany-seeds model, we do distinguish between the population ofsmall silicate grains and large icy pebbles, however. We comeback to the importance of condensation onto small grains ratherthan onto pebbles in the discussion (Sect. 6).

• Our model does not account for coagulation among pebbles: the enhancement of the solids-to-gas ratio is solely due to diffusion and condensation. We remain agnostic as to what happens to the pebbles in terms of coagulation and fragmentation until they reach the evaporation front. In the many-seeds model, we do not take into account coagulation of silicate particles. We will come back to this point in the discussion (Sect. 6).

• Our model is 1D, and has a straight vertical snowline. In reality, the vertical snowline is curved rather than straight, because the temperature varies with height above the midplane (e.g., Oka et al. 2011). However, Ros & Johansen (2013) have shown that particle growth due to diffusion and subsequent condensation of water across the vertical snowline is negligible compared to the same process across the radial snowline.

• Similarly, we assume that the scale height of the vapor – and, in the many-seeds implementation, of the free silicates – is always equal to the scale height of the disk; in other words, when an icy pebble evaporates, the released water vapor immediately distributes itself over the same vertical extent as the background gas, whereas it originates from the icy pebbles that reside in a thinner disk, characterized by a scale height H_peb: ${\mathit{H}}_{\mathrm{peb}}\mathrm{=}{\mathit{H}}_{\mathrm{gas}}\sqrt{\frac{\mathit{\alpha }}{{\mathit{\tau }}_{\mathrm{S}}\mathrm{+}\mathit{\alpha }}}\mathrm{·}$(33)The assumption of rapid vertical mixing extends to the calculation of the mean molecular weight of the gas in the complete model. Also there, it is assumed that the water vapor is distributed evenly throughout the gas at all times. In Sect. 4.2, we verify that the assumption of taking the vapor scale height equal to the disk scale height rather than the pebbles scale height, does not have a large effect on the results.

• We assume a steady-state gas distribution, that is, Ṁ_gas and T remain constant. This assumption is justified since the solids-to-gas enhancement takes place on a radial scale Δr that is small compared to the size of the disk. In other words, the local viscous evolution timescale (Δr²/ν), on which the effect happens, is much smaller than the global viscous evolution timescale of the disk, on which the gas accretion rate and global temperature profile change. The temperature profile (Eq. (1)) corresponds to a passively-irradiated disk. However, in case of viscous heating T(r) will still be described by a power-law (Armitage 2015). We have checked the effect of adopting a viscous heating temperature profile (T(r) ∝ r^{− 3 / 4}) in our numerical simulation. The key consequence is that the location of the snowline is shifted in radius; all other model results are very similar.

• We adopt an α-viscosity model, with constant α (the value of which we allow to vary between very high and very low values), independent of the solids-to-gas ratio. However, if turbulence is driven by the magneto-rotational instability, the local turbulence strength at the midplane depends on the local grain size and abundance. So-called “dead zones”, where turbulence is strongly reduced, could occur at locations in the disk with a sharp increase in the solids-to-gas ratio (Kretke & Lin 2007). Since our results depend strongly on the amount of turbulence near the snowline, the creation of dead zones would change our results.

• Similarly, we do not consider the interplay between (small) grains and the local temperature: we do not solve for the temperature profile self-consistently, but take it as constant, given by Eq. (1). However, in the many-seeds model, a large amount of small silicate grains is released in the evaporation front. This means that locally the opacity could increase. Conceivably, this could lead to a steeper temperature profile in the region near the snowline and therefore to a more narrow evaporation front. Including these effects requires a radiation transfer model, which is beyond the scope of this work.

• In the many-seeds implementation, we do not include condensation onto the small, free silicate particles. We also do not model the porosity of pebbles. We come back to these issues in Sect. 6.3.

3. Numerical implementation

We make use of the open source partial differential equations solver FiPy (Guyer et al. 2009) to solve the transport equations numerically. The coupled systems of equations are solved on a cylindrical 1D grid that ranges from 1 to 5 au.

We implement boundary conditions as follows. At the outer boundary r₂ = 5 au, the pebble mass flux ℱ_{s / g}Ṁ_gas is fixed³. We convert the input parameter τ₃ to a physical pebble size s_p,start using Eq. (7), and since our model does not include coagulation, we set the pebble size (which can be converted to a pebble mass m_p,start) at the outer boundary r₂ to this value. For the single-seed model, the mass of the bare silicate core is then given by m_core = ζm_p,start = m_p,start/ 2. At the inner boundary r₁ = 1 au, we implement outflow boundary conditions: solid particles and water vapor are accreted to the central star at a constant rate. We also demand that Σ_ice = 0 at r₁ and Σ_Z = 0 at r₂, for obvious reasons. Additionally, in the many-seeds implementation we have Σ_dirt = Σ_ice = 0 at r₁ and Σ_sil = Σ_Z = 0 at r₂.

3.1. Time-dependent

In our time-dependent method, we adopt small time steps (Δt ~ 100yr) and a small value for the allowed residuals, to ensure mass conservation by demanding the residuals⁴ to be small. This method requires a lot of computational power because the system of partial differential equations is numerically very stiff, but the advantage is that we find not only the steady-state solution to the equations, but are able to follow the solution in time. This allows us to get a measure of the time needed for the solution to settle into a steady-state. We present the time-dependent results for the “simple, single-seed” model with fiducial input parameters in Sect. 4.1.

3.2. Time-independent

Our time-independent method solves for the steady-state solution directly, by taking very large time-steps (in total 40 time steps that increase to 10⁷yr) and decreasing the residuals threshold after each run, taking the solution from the previous run as the new input. We continue this process until we reach convergence. We present the steady-state results for the “simple, single-seed”, “simple, many-seeds”, “complete, single-seed”, and “complete, many-seeds” models with fiducial input parameters in Sect. 4.2.

4. Results – fiducial model

In this section we introduce a standard set of parameters which we then use to illustrate different aspects of the model. The fiducial model is defined by the parameters listed in Table 1. We first discuss the time-dependent results for the simple, single-seed, fiducial model in Sect. 4.1. We then go on to discuss the steady-state results for the other variants of the fiducial model in Sect. 4.2.

4.1. Time-dependent solution

In Fig. 2 we plot several quantities as function of radial distance from the star r, at four points in time. In Fig. 2a, the surface densities of ice (solid lines) and vapor (dashed lines) are plotted. Initially (t = 0), the ice and vapor surface density are zero, Σ_ice = Σ_Z = 0. Icy pebbles are drifting in from the outer disk (right) to the inner disk (left). At the point where the icy component of the pebbles evaporate. At the snowline r_snow ≈ 2.1 au all the H₂O is in the gas phase. Since the vapor is advected with the gas velocity v_gas, which is much smaller than the pebble velocity v_peb, the vapor surface density (dashed line) quickly exceeds that of the ice (solid line). The vapor density increases until the steady-state vapor distribution Σ_Z,a has been reached. Σ_Z,a is given by⁵: ${\mathrm{\Sigma }}_{\mathit{Z,}\mathrm{a}}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{gas}}{\mathrm{3}\mathit{\pi \nu }}$(34)and is plotted by the dotted curve. Meanwhile, due to the outward diffusion of water vapor across the snowline and subsequent condensation onto the inward-drifting pebbles, the ice surface density just exterior to the snowline increases. Diffusion of vapor therefore results in a distinct bump in the ice profile and we denote the location corresponding to maximum Σ_ice the peak radius r_peak. A similar maximum is seen in the solids-to-gas ratio (Fig. 2b), which rises to ~0.45. After t = 2 × 10⁵ yr a steady state has been reached.

The enhancement of the pebbles-to-gas ratio is reflected by the typical pebble mass m_p and density ρ_•,p, plotted in Figs. 2c and d, respectively. Pebbles at the inner boundary are half as massive as pebbles at the outer boundary, because at the inner boundary all the ice has evaporated off the pebbles, and the dust fraction of pebbles ζ = 0.5. The typical pebble mass just outside the snowline increases over time, indicating that water vapor condenses onto the pebbles. This is also shown by the behavior of the typical pebble density ρ_•,p. At the outer boundary, ρ_•,p = 1.5gcm^-3, corresponding to pebbles containing as much water ice as silicates in mass; at the inner boundary, ρ_•,p = 3gcm^-3, corresponding to pure silicate pebbles; whereas just outside the snowline, ρ_•,p decreases over time to values close to unity – corresponding to pure water ice pebbles. A decrease in particle internal density outside the snowline due to water condensation was also found by Estrada et al. (2016).

In Fig. 3 we plot the growth rate of the peak in Σ_ice outside the snowline as function of time t. At the start of the simulation, ice is not yet present at the location of the peak. After the icy pebbles have reached the location of the eventual peak, at t ~ 2 × 10⁴yr, the rate of growth decreases exponentially – the peak continues to grow, but at an ever slower rate, while the system is converging to its steady-state solution. The e-folding growth timescale is ~3 × 10⁴yr, which is approximately equal to the gas advection timescale (r_peak − r_snow) /v_gas; the time it takes for the vapor to traverse the peak and contribute to the steep vapor concentration gradient along which outward diffusion can take place. We will come back to the peak formation timescale in Sect. 5. The time in which 90% of the peak forms is ~10⁵yr. Note that we have assumed that the incoming pebble mass flux does not change over time (ℱ_{s / g} is constant at the outer boundary r₂). If the pebble flux would decrease on a timescale shorter than the time it takes to form the peak, a steady-state would not exist. If the pebble flux would decrease on a timescale longer than it takes to form the peak, we would reach a quasi-steady-state solution; in that case our model is fully applicable.

4.2. Steady-state solution

The time-dependent solution described in Sect. 4.1 converges to the steady-state solution. In this section we take a closer look at the steady-state solution for the fiducial parameters, for the model designs discussed in Sect. 2.

4.2.1. Single-seed versus many-seeds

In Fig. 4 we compare the simple single-seed model (upper two panels) with the simple many-seeds model (lower two panels). The left panels show the surface densities of ice, vapor, and in the many-seeds case, also of locked silicates (dirt) and free silicates. The ice profiles are very similar in both models, but in the many-seeds case there is an additional peak in the dirt surface density. This is because the free silicates behave like vapor: they are released throughout the evaporation front, and follow the gas profile interior to the snowline. Turbulent diffusion not only causes vapor to be transported outward and condense onto the pebbles, but leads to the same effect for the small silicates – they diffuse outward and stick onto the pebbles, resulting in a peak in the dirt surface density distribution.

In the right panels we plot pebbles-to-gas midplane density ratios ρ_peb/ρ_gas and column density ratios Σ_peb/ Σ_gas. In the many-seeds model, the pebbles-to-gas ratio is the sum of the dirt-to-gas ratio and the ice-to-gas ratio. In both models, a clear peak is present just outside the snowline. The midplane density ratio is larger than the column density ratio due to settling of pebbles. In the many-seeds case, we also plot the free silicates-to-gas midplane density ratio ρ_sil/ρ_gas. This ratio is constant interior to the snowline because the silicates are perfectly coupled to the gas. In the many-seeds model, small silicate grains diffuse outward across the snowline, whereas in the single-seed model, there is no significant outward diffusion of the larger silicate pebbles due to the absence of a steep concentration gradient. Therefore, the peak in the pebbles-to-gas ratio outside the snowline is larger in the many-seeds model than in the single-seed model.

4.2.2. Simple versus complete models

Figure 5 shows the results for the complete model with the single-seed pebble design (upper two panels, a and b) and with the many-seeds pebble design (lower two panels, c and d). In Fig. 5a (single-seed) and Fig. 5c (many-seeds), we plot the surface densities of ice, vapor, and in the many-seeds case the surface densities of locked silicates (dirt) and free silicates (sil). Dashed lines denote the simple model results of Figs. 4a and c, for comparison. The grey lines correspond to the mean molecular weight of the gas. From the outer disk to the inner disk, the mean molecular weight increases from 2.34 m_H to 3.11 m_H across the snowline. Clearly, the peak in the ice surface density is higher and broader for the complete model than for the simple model. This is because collective effects (which reduce the pebble drift velocity) outweigh the effect of the enhanced gas pressure gradient in the evaporation front (which enhances the pebble drift velocity). In the complete model, therefore, pebbles effectively have a smaller radial velocity than in the simple model, and therefore the resulting ice peak is higher and broader. This is also the case in the many-seeds implementation. The width of the peak strengthens the justification of the vertical mixing assumption (see below).

In Figs. 5b and d we compare the midplane pebbles-to-gas ratios between the simple and the complete model. Collective effects help to boost the pebbles-to-gas ratio, for the reason outlined above.

4.2.3. Validation of vertical mixing assumption

As discussed in Sect. 2.8, we assume rapid vertical mixing of vapor. However, it could be argued that (part of) the released vapor recondenses onto icy pebbles before relaxing to the same vertical distribution as the gas, thereby boosting the condensation rate. This assumption also extends to the small silicate grains in the many-seeds model. Similarly, one could argue that the silicate grains can collide with a pebble before they can diffuse to higher vertical layers (e.g., Krijt & Ciesla 2016). In order to investigate what the difference is in terms of the steady-state ice distribution between two extreme cases, we take the simple, single-seed, fiducial model and run it once with H_Z = H_gas and once with H_Z = H_peb. We compare the steady-state outcomes of the two cases in Fig. 6. Even though the location of the snowline and of the pebbles-to-gas peak differs a bit between the two cases, the height of the peak in the midplane pebbles-to-gas ratio increases only by ~0.3% if one takes H_Z = H_peb instead of H_Z = H_gas. Similarly, assuming a smaller scale height for the silicate grains in the evaporation region would probably not significantly change the height of the peak in the “locked” silicates, since the silicate grains behave like vapor. In the rest of this paper, we assume H_Z = H_sil = H_gas.

5. Streaming instability conditions

Carrera et al. (2015) have investigated for what values of the solids-to-gas column density ratio streaming instability occurs, as function of the stopping time of the solid particles. They found that, typically, the solids-to-gas column density ratio should exceed ~2% for stopping times ~5 × 10^-2. However, these results should be interpreted with caution, since Carrera et al. (2015) only accounted for self-driven (Kelvin-Helmholtz) turbulence. We do include global turbulence, and therefore cannot simply use their conditions for streaming instability. A more robust condition for streaming instability in the presence of global turbulence is probably a midplane solids-to-gas ratio that reaches order unity (Johansen & Youdin 2007).

In this section we investigate under what disk conditions, the midplane solids-to-gas ratio near the snowline gets most enhanced due to the effect of water diffusion and condensation. To this end, we have constructed a semi-analytic approximate model to be able to quickly test for many different values of the input parameters (Sect. 2.7), within the single-seed implementation (Sect. 2.4). Before presenting the results, we will provide a short summary of the semi-analytical model, which is discussed in greater detail in Appendix C.

5.1. Semi-analytical model: Summary

Our semi-analytical model is an analytic prescription to find the approximate location, height and width of the ice surface density at the location of the ice peak, r_peak, of the steady-state (time-independent) solution (see Sect. 3.2 for our numerical method to find the steady-state solution to the transport equations). For simplicity, we have tailored the semi-analytical model towards the “single-seed” design. (An extension to also include the “many-seeds” design will be considered in a future work.)

The model works by iteratively solving a set of analytical equations, until convergence is reached. The key parameter is the ratio between the pebble velocity at r_peak and the normalized diffusivity at this location, ϵ = v_peb(r_peak) /D(r_peak)r_peak. Initially, we start with an estimate for ϵ, which can, for example, be taken from the zero-model solution (Sect. 2.7). For a certain ϵ the semi-analytical model then updates the surface density at the peak (Σ_peak) as well as r_peak (or rather the difference between r_ice and r_peak). With these updates we obtain a new value of ϵ and the procedure can be repeated until the fractional change in the parameters has become sufficiently low.

We emphasize that the semi-analytical model gives approximate solutions (correct at the ~20% level), but that they are very useful since they come at almost zero computational costs. It is therefore ideal to be applied to the parameter searches that we consider in this section. In Appendix C the model is described in detail and compared to runs of the numerical steady-state model.

5.2. Enhancement of the ice surface density

We first look at the relative effect of enhancing the ice surface density near the snowline: for each disk model, we normalize the resulting peak ice surface density by the ice surface density at the peak radius in case there is only advection: ${\mathit{f}}_{\mathrm{\Sigma }\mathit{,}\mathrm{peak}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{ice}\mathit{,}\mathrm{peak}}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\zeta }\mathrm{)}{\mathrm{ℱ}}_{\mathrm{s}\mathit{/}\mathrm{g}}\mathit{Ṁ}\mathrm{gas}\mathit{/}\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{peak}}{\mathit{v}}_{\mathrm{peb}}}$(35)where v_peb is calculated in the absence of condensation and evaporation. As before, we take ζ = 0.5. In Fig. 7a we plot contours of constant f_Σ,peak as function of α and ℱ_{s / g}, where Ṁ_gas is fixed at 10^-8M_⊙yr^-1 and τ₃ is fixed at 0.03. Clearly, the higher the solids-to-gas accretion rate ℱ_{s / g}, the higher f_Σ,peak – at fixed α the normalized ice peak increases from low ℱ_{s / g}-values to high ℱ_{s / g}-values, due to collective effects playing an increasingly important role. The plot also shows that disks with α-values of about 10^-3 are best at enhancing the ice surface density with respect to the advection-only expected ice surface density. This is because higher values of α lead to a broader and relatively lower peak, while lower α-values, indicating larger gas densities, move pebbles into the Stokes drag regime. The fact that the relative height of the ice peak becomes smaller when pebbles enter the Stokes regime, is explained by the pebble size-dependency of the drift velocity. In the Epstein regime the drift velocity depends linearly on the pebble size, whereas in the Stokes regime the drift velocity is proportional to the square of the pebble size. This means that when pebbles enter the Stokes regime, their stopping time – and hence their drift velocity – increases rapidly. The effect of this is that f_Σ,peak decreases rapidly. Therefore, the transition from the Epstein regime to the Stokes drag regime is characterized by closely-spaced contours.

Figure 7b shows contours of constant f_Σ,peak, as function of α and Ṁ_gas. ℱ_{s / g} is fixed at 0.4 and τ₃ is fixed at 0.03. The cyan lines correspond to a Toomre parameter (Q_T) of unity at 10 au (1 au). Q_T is defined as: ${\mathit{Q}}_{\mathrm{T}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}\mathrm{\Omega }}{\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{gas}}}$(36)where G is the gravitational constant. A disk is gravitationally unstable if Q_T< 1. In Fig. 7b, models in the parameter space below the cyan line have Q_T< 1 at 10 au (1 au). In Fig. 7b densely-spaced contours again indicate the transition between Epstein and Stokes drag.

The enhancement of the ice surface density near the snowline is relatively modest, but enrichment factors of a few can still be important for triggering streaming instability (Johansen et al. 2009). Besides, Fig. 7 corresponds to our single-seed model, which is conservative. Larger enhancements can be achieved in the framework of the many-seeds model (see Fig. 4). We can imitate the many-seeds model in the analytical model of Appendix C, crudely, when we let ζ → 0, which assumes that micron-sized silicate grains behave like vapor. In that case we see f_Σ,peak increase by another factor of two (results not plotted). Therefore, total enhancements of f_Σ,peak ~ 5–10 are feasible, which significantly boost the likelihood of triggering streaming instability.

5.3. Solids-to-gas ratios

Let us now look at the absolute results, in terms of the midplane solids-to-gas ratio. In Fig. 8a, black lines denote contours of constant peak midplane solids-to-gas ratios, as function of α and ℱ_{s / g}. We again fix Ṁ_gas = 10^-8M_⊙yr^-1 and τ₃ = 0.03. The orange contours correspond to the total amount of solids in units of Earth masses that is needed to form the eventual peak. This quantity (M_{solids,needed}) is calculated as: ${\mathit{M}}_{\mathrm{solids}\mathit{,}\mathrm{needed}}\mathrm{~}\frac{{\mathrm{ℱ}}_{\mathrm{s}\mathit{/}\mathrm{g}}\mathit{Ṁ}\mathrm{gas}\mathrm{(}{\mathit{r}}_{\mathrm{peak}}\mathrm{-}{\mathit{r}}_{\mathrm{snow}}\mathrm{)}}{{\mathit{v}}_{\mathrm{gas}}}$(37)where (r_peak − r_snow) /v_gas is the typical timescale on which the peak forms, as discussed in Sect. 4.1. The value of M_{solids,needed} is similar to the amount of solids that gets “locked up” in the peak⁶. Figure 8a shows that the higher the value of ℱ_{s / g}, the easier streaming instability can be triggered, as one may expect. It also shows that in order to reach midplane solids-to-gas ratios of ~10%, a solids-to-gas accretion rate ℱ_{s / g} of the same order is required. This also holds for other values of the gas accretion rate Ṁ_gas and initial pebble size τ₃.

Figure 8b is the same as Fig. 7b, except that the black contours now denote constant values of the peak midplane solids-to-gas ratio. Just like in Fig. 8a, orange contours denote constant M_{solids,needed}, but now with Ṁ_gas as the parameter on the x-axis. The region with closely-spaced contours in Fig. 8a and the oblique regions in the otherwise fairly horizontal contours in Fig. 8b are again due to the transition between the Epstein and Stokes drag regimes.

Up to this point we have kept the parameter τ₃ constant at its fiducial value of 0.03. Recently, Yang et al. (2016) have shown that the streaming instability mechanism can also work for small particles (τ ~ 10^-3). In Fig. 9 we now vary τ₃ and plot the resulting peak midplane solids-to-gas ratio, whilst fixing Ṁ_gas at 10^-8M_⊙yr^-1, ℱ_{s / g} at 0.4, and α at 3 × 10^-3. Different regimes in the plot are denoted by Roman numbers I–IV. In region I, the particles are well-coupled to the gas. The solids-to-gas ratio flattens out towards low τ₃-values, because radial drift and settling do not play a role for these particles. The boundary between regime I and regime II corresponds to the value of α. The highest solids-to-gas ratios are reached for particles with stopping times of the same order as α or smaller. For higher values of τ, the larger drift velocity (proportional to τ) starts to dominate over the larger degree of settling (proportional to $\sqrt{\mathit{\tau }}$). Therefore, in regime II, the solids-to-gas ratio decreases with increasing τ. The transition from regime I and regime II corresponds to the τ₃-value for which ϵ is unity. The ϵ parameter (introduced in Sect. 5.1) increases with increasing τ₃, and the peak surface density decreases at a faster rate with increasing ϵ for ϵ> 1 than for ϵ< 1 (Eq. (C.13)). The physical explanation is that the solution changes from diffusion-dominated (ϵ< 1) to drift-dominated (ϵ> 1). Therefore, the solids-to-gas ratio decreases more steeply in regime III than in regime II. The final transition, from regime III to regime IV, is caused by the transition from the Epstein regime to the Stokes regime. For τ₃ values of ~1, the enhancement effect becomes minimal and the peak solids-to-gas ratio flattens out.

From Figs. 8a and b we learn that in the context of the α-viscosity disk model, the largest solids-to-gas ratios near the snowline are reached for high α-values. This means that our diffusion model does not require a very quiescent disk to achieve high solids-to-gas ratio, and also works in the early disk phase when α and Ṁ_gas are presumably high. Our results show that for high α and high Ṁ_gas, tens of Earth masses are locked up in pebbles in an annulus outside the snowline with a width of the order of 1 au.

6. Discussion

6.1. Conditions for streaming instability

Our results show that the ice surface density can be enhanced by a factor ~3–5 outside the snowline, due to the effect of water diffusion and subsequent condensation onto icy pebbles (adopting ζ< 0.5 or the many-seeds model increases the enhancement by another factor ~2). This differs from the results of Stevenson & Lunine (1988), who found that the same effect could lead to an enhancement by as much as 75. However, they considered a static system, whereas our model accounts for radial drift of pebbles and gas accretion onto the star. Even though we found a modest enhancement, we showed that the pebble density at the midplane can reach tens of percent of the gas density, under the condition that the solids-to-gas accretion rate ℱ_{s / g} is of order 0.1.

But what are realistic values for ℱ_{s / g}? Ida et al. (2016) have calculated the pebble mass flux Ṁ_peb in a viscously-evolving disk, taking into account dust growth and radial drift. They found that after 10⁶ yr, ℱ_{s / g} ~ 0.3 for a turbulence parameter α = 10^-3 and a gas accretion rate Ṁ_gas = 10^-8M_⊙yr^-1. Since this result suggests ℱ_{s / g}-values of order 0.1 are plausible, water condensation can plausibly trigger streaming instability near the snowline. According to Ida et al. (2016), Ṁ_peb is inversely proportional to α. This implies that even though our results show that for high α-values, water diffusion and condensation leads to the largest solids-to-gas ratio at the snowline, a smaller pebble mass flux is expected, which reduces the effect. The pebble mass flux is, however, dependent on time t as Ṁ_peb ∝ t^{− 1 / 3}, meaning that the pebble mass flux earlier in the disk is higher than the fiducial 30% of the gas accretion rate quoted before. Because higher α-values lead to larger solids-to-gas ratios but to lower pebble mass fluxes, we conclude that intermediate α-values (~10^-3) are most favorable to trigger streaming instabilities. The fact that the diffusion-driven model of this work already produces enhancement of solids at intermediate α-values and high Ṁ_gas, implies that planetesimal formation does not have to await the arrival of quiescent disk conditions. Provided a sufficiently large pebble flux, our results imply that planetesimals can form early.

6.2. Interior or exterior?

In Saito & Sirono (2011), it is argued that the evaporation of icy pebbles leads to a large enhancement of small silicate grains interior to the snowline. In their model, micron-sized silicate grains are instantaneously released at the snowline, and then follow the radial motion of the gas. Since the gas moves at a much lower velocity than the drift velocity of the icy pebbles, the silicate grains pile-up and become gravitationally unstable. A similar mechanism for the formation of planetesimals in the inner disk has recently been presented by Ida & Guillot (2016).

We confirm that the solids-to-gas ratio interior to the snowline is boosted due to the slow radial motion of the small silicates (see Fig. 4c), but find smaller enhancements than Saito & Sirono (2011) and Ida & Guillot (2016). These works assumed that the evaporated silicate grains remained in the same vertical layer as the icy pebbles, which had settled to the midplane due to their larger size. This is a natural assumption if one adopts instantaneous evaporation. However, in our model silicate grains are not instantaneously released from the pebbles at a single location, but are instead released across the ice evaporation front – the region between the snowline and the peak radius. With the exception of very low α-values, the width of the ice peak is (much) larger than the gas scale height. Consequently, the timescale for vertical mixing is smaller than the timescale on which the silicate grains traverse the evaporation front, justifying our assumption that silicate grains (and water vapor) are distributed over H_gas. Allowing for the vertical mixing of silicate grains is why we found smaller solids-to-gas ratios interior to the snowline than Saito & Sirono (2011) and Ida & Guillot (2016).

In this work we did not model the porosity of pebbles, but assumed pebbles are compact spheres. For homogeneous porous spheres, the results are the same as for compact spheres with the same stopping time, if they are in the Epstein regime. This is because the evaporation and condensation rates (Eqs. (13), (14)) depend on ${\mathit{s}}_{\mathrm{p}}^{\mathrm{2}}\mathit{/}{\mathit{m}}_{\mathrm{p}}$ (the surface to mass ratio of a pebble), which is effectively the Epstein stopping time (Eq. (7)). For fractal aggregates, however, the evaporation rate could conceivably be increased due to a larger total surface area available for evaporation. In that case, the location of the snowline (r_snow) and of the peak (r_peak) are pushed to larger radial distances, because drifting pebbles start to evaporate earlier. The shape of the solids distribution profile and hence the width of the evaporation front, however, remain similar. We have confirmed this behavior by running our simulation with an evaporation rate that has been arbitrarily increased by a factor 100 compared to Eq. (13). Therefore, a fractal interior structure of the pebbles would not much change the extent of the region in which pebbles evaporate. We note however that in the many-seeds model, a large increase in opacity due to small silicate grains might lead to a steeper temperature profile than assumed here, and therefore to a more narrow evaporation front.

Even though we found smaller enhancements interior to the snowline, adopting the many-seeds model leads to larger enhancements outside the snowline compared to the single-seed model. This is because outward radial diffusion and subsequent sticking of silicate grains to icy pebbles adds to the ice peak and therefore enforces the enhancement of the solids-to-gas ratio outside the snowline (see Fig. 5).

6.3. Condensation onto small grains

In our single-seed model, we adopted a single-size assumption for the solid particles, applicable to a drift-limited scenario, in which particles grow until they decouple from the gas and start to drift inward (see, e.g., Birnstiel et al. 2010; Krijt et al. 2016b). Because the solids density is larger close to the star than in the outer disk, particles close to the star start to drift inward earlier than particles in the outer disk. This leads to an inside-out clearing of the solids component of the disk. In this case, one would not expect small particles to be present in significant amounts near the snowline at the time when icy pebbles approach the snowline, drifting in from the outer disk. Therefore, in the drift-limited case the single-size approximation for pebbles is valid.

The single-size assumption does not hold for a fragmentation-limited size distribution, however. In that case, pebbles dominate the total mass of solids, but small particles dominate the total surface area (see, e.g., Birnstiel et al. 2012). Then, condensation of water vapor will predominantly happen onto small particles, instead of onto pebbles as considered in this work. Also, in our many-seeds model, we do not allow for condensation onto the small grains that get released from icy pebbles upon evaporation.

Within the single-seed framework, we can mimic the situation where condensation happens onto small grains rather than onto pebbles by adopting a very small typical stopping time just outside the snowline (τ₃). This results in approximately the same or even larger peaks, as illustrated by Fig. 9. In that sense, assuming condensation occurs solely onto pebbles, rather than onto small grains, is a conservative choice. However, a more correct approach would allow for condensation onto small grains in the many-seeds model, in which icy pebbles break up into small grains upon evaporation, or, alternatively, within a fragmentation-limited scenario by adopting a full size distribution. Such improvements will be considered in future work.

6.4. Relevance of coagulation

The many-seeds model ignores coagulation among the small silicate grains, as well as the possibility that water vapor condenses on to silicates. Because there are so many of them, the average silicate grain size would not increase significantly due to water condensation. In that case, the assumption that they are well-coupled to the gas would still hold. However, growth of silicates might also occur due to coagulation. Growth increases the drift velocity, thereby reducing the effect of the outward-diffusion and recondensation of the silicate component. Consequently, coagulation among silicates would bring the many-seeds model closer to the single-seed model due to a decrease in the enhancement of “locked” silicates outside the snowline. A similar trend is expected when the silicates that are released upon evaporation are larger. According to the experimental work of Aumatell & Wurm (2011), an icy aggregate of 1 cm in size breaks up into about 200 sub-aggregates when it evaporates. In any case, we expect that a more realistic “break-up” or coagulation model will lead to results that lie in between the results for our many-seeds and single-seed model designs.

In our work, we have also neglected coagulation among pebbles. On the one hand, one could think that icy pebbles on their way to the snowline would grow even more because of coagulation – at least as long their relative velocities stay low enough. On the other hand, according to Sirono (2011) and Okuzumi et al. (2016), the fragmentation threshold velocity for icy pebbles near the snowline is reduced due to sintering. This might motivate a smaller average icy pebble size (a lower τ₃ in our nomenclature) and therefore a smaller radial drift velocity, which would further increase the solid surface density exterior to the snowline (see Fig. 9).

6.5. Observational implications

We showed that intermediate-to-high α-values are most favorable for triggering streaming instability near the snowline. The onset of streaming instability manifests itself through the formation of a peak in the solids-to-gas ratio that grows in height (increasing Σ_ice) but also in width – the latter because collective effects cause pebbles to move slower. Even though more solids are required for high α and high Ṁ_gas than for low α and low Ṁ_gas, tens of Earth masses are sufficient to form the ice peak. This also means that tens of Earth masses can be stored in an annulus outside the snowline, especially when the system is characterized by high pebble and gas accretion rates. Such annuli of solid enhancement can be observable with facilities as ALMA. In the context of the snowline, the ring structure seen in the TW Hya system with a width of about 1 au (Andrews et al. 2016) might very well correspond to a realisation of the water-diffusion effect in a highly turbulent disk with a high gas accretion rate.

The fact that the planets in the inner solar system are water-poor (Marty 2012; McCubbin et al. 2012), even though the snowline should have migrated to 1 au before the disk was depleted, was pointed out as the “snowline problem” by Oka et al. (2011). The early formation of a protoplanet near the snowline might provide a solution to this problem, as proposed recently by Morbidelli et al. (2016). Once a protoplanet near the snowline has formed it could halt the inward-drifting pebble flux, either by accreting the pebbles (Guillot et al. 2014), or by trapping them in pressure maxima created by the newborn planet (Zhu et al. 2014; Lambrechts et al. 2014). In the meantime, the water vapor gets accreted to the star and the snowline location migrates inward (e.g., Oka et al. 2011). Consequently, an early formation of a protoplanet near the snowline might explain the lack of water in the inner solar system. The water diffusion/condensation effect discussed in this paper provides a way to form planetesimals near the snowline in an early stage of the disk (i.e., when the snowline is still located outside of the current water-poor region of the solar system). After planetesimals of a certain size have formed, they can subsequently accrete pebbles and quickly grow to larger bodies (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014; Bitsch et al. 2015; Visser & Ormel 2016). Therefore, our water diffusion/condensation model provides the first key step to realize the early formation of a protoplanet near the snowline.

7. Conclusions

Our main findings can be summarized as follows.

• 1.

  Water diffusion and condensation near the snowline can resultin an enhancement of a factor several in the ice surface density.Because of the dynamic setup of our model where water vapor iscarried away with the accreting gas, this enhancement is muchless than found by Stevenson & Lunine (1988).Nevertheless, the boost in the solids-to-gas ratio can still triggerstreaming instability, provided a large enough pebble flux.

• 2.

  The peak in the ice surface density is not located at the snowline, but exterior to it (r_peak>r_snow). With larger incoming pebble flux, the peak becomes broader, because the back-reaction of the solids on the gas reduces the pebble drift velocity. Depending on the turbulence strength, the width of the peak can be as much as ~1 au.

• 3.

  Such broad peaks can contain tens of Earth masses in pebbles, appearing as bright annuli at radio wavelengths.

• 4.

  The release of many micron-sized silicate grains upon evaporation of icy pebbles produces a peak in “locked” silicates exterior to the snowline, due to diffusion and sweep-up of the silicates. This peak adds to the ice peak and therefore enforces the enhancement of the solids-to-gas ratio outside the snowline.

• 5.

  Interior to the snowline, the solids-to-gas ratio is boosted if many micron-sized silicate grains are released during evaporation, because they cause a “traffic jam” effect. However, because silicate grains mix with the background gas before they cross the snowline, the solids-to-gas ratio enhancement just interior to r_snow is limited.

• 6.

  In the context of a viscous disk model, the ratio between the accretion rate of solids and gas needs to be of order ~0.1 in order to reach solids-to-gas midplane ratios of the order of tens percent. The mechanism operates best at intermediate (~10^-3) α-values. Therefore, planetesimals can form at an early time in the evolution of the disk.

---

Acknowledgments

D.S. and C.W.O. are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We would like to thank Jonathan Guyer and Daniel Wheeler for their elaborate answers to FiPy-related questions. We would also like to thank Carsten Dominik, Shigeru Ida, Satoshi Okuzumi, Anders Johansen, Daniel Carrera and the other participants of the Lorentz workshop “New directions in Planet Formation” (Leiden, July 2016) for stimulating discussions. We are indebted to Sebastiaan Krijt, Jeff Cuzzi, and the anonymous referee for constructive feedback on the manuscript.

References

Appendix A: List of symbols

A list of frequently-used symbols is given in Table A.1.

Appendix B: Expression collective effects

Expressions for the aerodynamic drift of single particles have been derived by Whipple (1972) and Weidenschilling (1977) for non-accreting disks. Later, Nakagawa et al. (1986) accounted for the backreaction forces that the particles collectively exert on the gas. The gas motion is then outwards while the total dust+gas mass flux becomes zero. These expressions can be generalized to account for a size distribution of particles (Tanaka et al. 2005; Estrada & Cuzzi 2008; Bai & Stone 2010).

In (steady) accretion disks gas moves inwards at a rate that is set by the viscosity. This modifies the single particle drift expressions. In particular, small particles, which show negligible radial drift, nevertheless move inwards as they are carried by the gas. While for the single particle case the modification of the drift expressions for accreting disks is straightforward, expressions accounting for both collective effects and viscosity have, to the best of our knowledge, not yet been presented. We derive those here.

The equations of motions in a frame co-rotating with the local Keplerian period read: $\begin{array}{ccc}\frac{\mathrm{D}{v}}{\mathrm{D}\mathit{t}}& \mathrm{=}& \mathrm{-}\frac{{v}\mathrm{-}{u}}{{\mathit{t}}_{\mathrm{stop}}}\mathrm{-}\mathrm{2}{{\Omega }}_{\mathrm{K}}\mathrm{\times }{v}\mathrm{+}{{F}}_{\mathrm{Euler}\mathrm{-}\mathrm{dust}}\\ \frac{\mathrm{D}{u}}{\mathrm{D}\mathit{t}}& \mathrm{=}& \frac{{\mathit{\rho }}_{\mathit{p}}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{{v}\mathrm{-}{u}}{{\mathit{t}}_{\mathrm{stop}}}\mathrm{-}\mathrm{2}{{\Omega }}_{\mathrm{K}}\mathrm{\times }{u}\mathrm{+}{{F}}_{\mathrm{Euler}\mathrm{-}\mathrm{gas}}\mathrm{+}{{F}}_{\mathrm{\nabla P}}\mathrm{+}{{F}}_{\mathit{\nu }}\end{array}$where v is the particle velocity and u the gas velocity, D / Dt is the material derivative, the first term on the RHS is the drag force, the second term is the Coriolis force, and F_Euler = − (dΩ_K/ dt) × r is the Euler force that arises due to the radial motion. It therefore appears in the radial equations. In addition, the equation of motion for the gas includes a hydrostatic correction due to a radial pressure gradient (F_∇P) and viscous forces (F_ν). Following convention, we write the pressure gradient in terms of a nondimensional η: ${{F}}_{\mathrm{\nabla }\mathit{P}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{\equiv }\mathrm{2}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}{\mathrm{\Omega }}_{\mathrm{K}}{{e}}_{\mathit{r}}$(B.2)and use the thin disk approximation (consider column densities Σ_gas instead of ρ) for the viscous force: ${{F}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\mathrm{\nabla }\mathrm{·}{T}$(B.3)where T is the viscous stress tensor. In cylindrical coordinates, the only relevant term is: ${\mathit{T}}_{\mathit{r\phi }}\mathrm{=}\mathit{\nu }{\mathrm{\Sigma }}_{\mathrm{gas}}\mathit{r}\frac{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{K}}}{\mathrm{d}\mathit{r}}\mathit{,}$(B.4)with which the viscous force becomes ${{F}}_{\mathit{\nu }}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\mathrm{\nabla }\mathrm{·}{T}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\mathit{\nu }{\mathrm{\Sigma }}_{\mathrm{gas}}\mathit{r}\frac{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{K}}}{\mathrm{d}\mathit{r}}\right){e}\phi$(B.5)where ν is the kinematic viscosity. In the α-disk model we have that νΣ_gas is constant and therefore ${{F}}_{\mathit{\nu }}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\nu }{\mathrm{\Omega }}_{\mathrm{K}}}{\mathrm{4}\mathit{r}}{{e}}_{\mathit{\phi }}\mathit{.}$(B.6)In order to solve for the drift velocities, we put the accelerations on the LHS of Eq. (B.1) zero, Du/ Dt = Dv/ Dt = 0. This is justified, because a change in the drift velocities occurs on a timescale of ~r/v_r, which is always much longer than the stopping time. Inserting Eqs. (B.2) and (B.6) into (B.1) gives a system of four equations and four unknowns (the velocities). Its solution is: $\left(\begin{array}{c}\\ {\mathit{v}}_{\mathit{r}}\\ {\mathit{v}}_{\mathit{\phi }}\\ {\mathit{u}}_{\mathit{r}}\\ {\mathit{u}}_{\mathit{\phi }}\end{array}\right)\mathrm{=}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\xi }{\mathrm{)}}^{\mathrm{2}}}\left(\begin{array}{c}\\ \mathrm{-}\mathrm{2}\mathit{\eta }{\mathit{\tau }}_{\mathrm{s}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\xi }\mathrm{)}\\ \frac{\mathrm{3}}{\mathrm{4}}\begin{array}{c}\mathrm{˜}\\ \mathit{\nu }\end{array}{\mathit{\tau }}_{\mathrm{s}}\mathrm{-}\mathit{\eta }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\xi }\mathrm{)}\\ \mathrm{2}\mathit{\eta }{\mathit{\tau }}_{\mathrm{s}}\mathit{\xi }\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\tau }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathit{\xi }\mathrm{)}\\ \mathrm{-}\frac{\mathrm{3}}{\mathrm{4}}\begin{array}{c}\mathrm{˜}\\ \mathit{\nu }\end{array}{\mathit{\tau }}_{\mathrm{s}}\mathit{\xi }\mathrm{-}\mathit{\eta }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\tau }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathit{\xi }\mathrm{)}\end{array}\right)\mathit{v}K$(B.7)where τ_s = t_stopΩ_K, ξ = ρ_d/ρ_g, and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\nu }\end{array}\mathrm{=}\mathit{\nu }\mathit{/}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{K}}$. Note that in the α-prescription $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\nu }\end{array}\mathrm{~}\mathit{\alpha \eta }$, that is, for α ≪ 1 the modification due to viscosity is minor.

Appendix C: Semi-analytical model

In this section we aim for an analytic prescription of the key characteristics of the steady-state ice and vapor distribution for any given disc model – not necessarily with the same gas and temperature profiles as assumed in this study. We will derive expressions for the peak of the ice surface density and its location, which enable us to determine whether the conditions for streaming instability beyond the snowline can be reached.

Appendix C.1: Location of the snowline

The location of the snowline, r_snow, is the outer-most radius where the entire incoming ice flux can thermodynamically exist in the gaseous state. Its location is found by equating the equilibrium (saturation) pressure to the vapor pressure of H₂O corresponding to the incoming icy pebble flux, that is, P_eq = P_Z,a, where P_Z,a is the steady-state vapor pressure; i.e., the vapor pressure obtained when all the H₂O is in the gas phase. Therefore, the snowline r_snow is given by the radius r where ${\mathit{P}}_{\mathrm{eq}\mathit{,}\mathrm{0}}\exp \left[\mathrm{-}\frac{{\mathit{T}}_{\mathrm{a}}}{\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right]\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }}\frac{{\mathrm{\Sigma }}_{\mathit{Z,}\mathrm{a}}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{\mu }}\frac{\mathit{Ṁ}\mathrm{ice}\mathit{/}\mathrm{3}\mathit{\pi \nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}$(C.1)see Eqs. (11) and (34). Given Ṁ_ice, α and a disk model, we can therefore readily solve (numerically) for r_snow.

A caveat of Eq. (C.1) is that the role of transport processes are neglected. These could potentially introduce disequilibrium corrections. However, we find that the times to achieve the equilibrium (saturation) pressure is typically much shorter than the transport timescales (by diffusion or drift).

Appendix C.2: The equilibrium density

Generalizing the above arguments to locations beyond r_snow, we approximate the vapor pressure with the saturation pressure – an approximation that is again justified as long as the timescales for evaporation (and condensation) are sufficiently short (in other words: the super-saturation level is small). Similarly, we can define a vapor surface density (Σ_eq) corresponding to P_eq. Re-arranging Eq. (C.1) in terms of Σ we obtain this equilibrium surface density: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{eq}}& \mathrm{\equiv }& \frac{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{gas}}\mathit{\mu }}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{P}}_{\mathrm{eq}}\\ & \mathrm{=}& \end{array}$(C.2)where Σ_snow and T_snow are the vapor surface density and the temperature at the location of the snowline. In Eq. (C.2) we have assumed power-laws for the temperature, gas surface density, and the gas scaleheight. In our case we have p = 7 / 4 and q = 1 / 2 (see Sect. 2.1).

Taylor-expanding the (weakly-varying) power-laws involved in the expression of Σ_eq around r_snow gives⁷: ${\mathrm{\Sigma }}_{\mathrm{eq}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\approx }{\mathrm{\Sigma }}_{\mathrm{snow}}\left(\mathrm{1}\mathrm{+}\mathit{px}\right)\exp \left[\mathrm{-}\frac{\mathit{q}{\mathit{T}}_{\mathrm{a}}\mathit{x}}{{\mathit{T}}_{\mathrm{snow}}}\right]\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{snow}}\left(\mathrm{1}\mathrm{+}\mathit{px}\right){\mathrm{e}}^{\mathrm{-}{\mathit{a}}_{\mathrm{eq}}\mathit{x}}$(C.3)where x = (r − r_snow) /r_snow is the fractional distance beyond the snowline and ${\mathit{a}}_{\mathrm{eq}}\mathrm{=}\frac{\mathit{q}{\mathit{T}}_{\mathrm{a}}}{{\mathit{T}}_{\mathrm{snow}}}\mathrm{\approx }\mathrm{15}\left(\frac{\mathit{q}}{\mathrm{0.5}}\right){\left(\frac{{\mathit{T}}_{\mathrm{snow}}}{\mathrm{200}\mathrm{K}}\right)}^{-1}$(C.4)is a dimensionless number that weakly depends on the location of the snowline.

In Fig. C.1 the respective vapor density profiles have been plotted. Parameters correspond to those of the default model (Table 1). The simulated result of Sect. 4 is plotted by the thick gray curve. It is characterized by a sharp dent at r_snow. Exterior to the snowline Σ_Z follows the equilibrium profile (solid black curve), while interior to r_snow, it is limited by the imposed mass flux (solid dashed curve). The approximate expressions derived in Eq. (C.3), valid for r>r_snow, are shown by the blue curves. We obtain r_snow = 2.09 au, Σ_snow = 92 g cm^-2 and a_eq = 16.9. The solid blue curve approximates Σ_eq very closely. The dashed blue curves gives a further approximation to Σ_eq, obtained by putting p = 0. This slightly underestimates Σ_Z but nevertheless gives a reasonable approximation that we will use below.

Assuming Σ_Z = Σ_eq for the vapor density beyond r_snow, we write for the associated mass flux ℳ_Z:$\begin{array}{ccc}{\mathrm{ℳ}}_{\mathit{Z}}\mathrm{\approx }{\mathrm{ℳ}}_{\mathrm{eq}}& \mathrm{=}& \mathrm{-}{\mathit{v}}_{\mathrm{gas}}{\mathrm{\Sigma }}_{\mathrm{eq}}\mathrm{-}{\mathit{D}}_{\mathrm{gas}}{\mathrm{\Sigma }}_{\mathrm{gas}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left(\frac{{\mathrm{\Sigma }}_{\mathrm{eq}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\right)\\ & \mathrm{\approx }& \mathrm{-}{\mathit{v}}_{\mathrm{gas}}{\mathrm{\Sigma }}_{\mathrm{eq}}\mathrm{-}{\mathit{D}}_{\mathrm{gas}}\frac{\mathrm{d}{\mathrm{\Sigma }}_{\mathrm{eq}}}{\mathrm{d}\mathit{r}}\end{array}$(C.5)where we neglected the curvature of the disk, i.e., we assumed that the scales over which Σ_eq and Σ_gas change are much smaller than r. This is justified as long as x ≪ 1. In a similar vein, we assume that quantites as D_gas and v_gas do not vary. With Eq. (C.3) for Σ_eq we obtain: $\begin{array}{ccc}\frac{{\mathrm{ℳ}}_{\mathrm{eq}}}{{\mathit{D}}_{\mathrm{gas}}\mathit{/}{\mathit{r}}_{\mathrm{snow}}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{v}}_{\mathrm{gas}}{\mathit{r}}_{\mathrm{snow}}}{{\mathit{D}}_{\mathrm{gas}}}{\mathrm{\Sigma }}_{\mathrm{eq}}\mathrm{-}\frac{\mathrm{d}{\mathrm{\Sigma }}_{\mathrm{eq}}}{\mathrm{d}\mathit{x}}\\ & \mathrm{=}& {\mathrm{\Sigma }}_{\mathrm{snow}}\mathrm{(}\mathrm{-}{\mathit{b}}_{\mathrm{gas}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{px}\mathrm{)}\mathrm{+}{\mathit{a}}_{\mathrm{eq}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{px}\mathrm{)}\mathrm{-}{\mathit{p}}^{\mathrm{)}}{\mathrm{e}}^{\mathrm{-}{\mathit{a}}_{\mathrm{eq}}\mathit{x}}\\ & \mathrm{\approx }& {\mathrm{\Sigma }}_{\mathrm{snow}}\mathrm{(}{\mathit{a}}_{\mathrm{eq}}\mathrm{-}{\mathit{b}}_{\mathrm{gas}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}{\mathit{a}}_{\mathrm{eq}}\mathit{x}}\end{array}$where b_gas = v_gasr_snow/D_gas is the dimensionless gas velocity and in the last step p = 0. In a viscously-evolving disk v_gas = 3ν/ 2r and hence b_gas = 1.5.

Appendix C.3: The ice profile

With the equilibrium profile as the approximation to the vapor density, we solve for the surface density of the ice by invoking conservation of mass. In steady-state the total H₂O mass flux ℳ_tot,ice = ℳ_ice + ℳ_Z is constant. Interior to the snowline, ℳ_Z = ℳ_tot,ice while for r ≫ r_snowℳ_Z ≪ ℳ_ice ≈ ℳ_tot,ice. But near the snowline the mass fluxes change rapidly. Here, the steep gradient in the vapor density causes an outward (positive) mass flux, which increases the inwardly-directed icy pebble flux: $\mathrm{-}{\mathrm{ℳ}}_{\mathrm{ice}}\mathrm{\equiv }{\mathrm{\Sigma }}_{\mathrm{ice}}{\mathit{v}}_{\mathrm{ice}}\mathrm{+}{\mathit{D}}_{\mathrm{p}}\frac{\mathrm{d}{\mathrm{\Sigma }}_{\mathrm{ice}}}{\mathrm{d}\mathit{r}}\mathrm{=}{\mathrm{ℳ}}_{\mathit{Z}}\mathrm{-}{\mathrm{ℳ}}_{\mathrm{tot}\mathit{,}\mathrm{ice}}\mathit{.}$(C.8)(Note that just outside r_snow, ℳ_ice and ℳ_tot,ice are negative and ℳ_Z is positive.) With ℳ_eq for ℳ_Z Eq. (C.8), in terms of x, reads $\frac{\mathrm{d}{\mathrm{\Sigma }}_{\mathrm{ice}}}{\mathrm{d}\mathit{x}}\mathrm{+}{\mathit{b}}_{\mathrm{peb}}{\mathrm{\Sigma }}_{\mathrm{ice}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{snow}}}{{\mathit{D}}_{\mathrm{p}}}\mathrm{(}{\mathrm{ℳ}}_{\mathrm{eq}}\mathrm{-}{\mathrm{ℳ}}_{\mathrm{tot}}\mathrm{)}$(C.9)where b_peb = r_snowv_peb/D_p is the ratio between the speed at which the pebbles drift in and the radial velocity of the gas. Usually, pebbles outpace the gas and b_peb ≫ 1.

To solve Eq. (C.9) we will assume that b_peb and other parameters are constant in x. Assuming constant b_peb is clearly an approximation, since the pebbles may (i) acquire thick icy mantles on their approach to the snowline but (ii) lose all their ice once their are very close to r_snow. Hence, the aerodynamical properties of the pebbles and therefore v_peb are expected to vary considerably. Disregarding these concerns, momentarily, the solution to Eq. (C.9) reads: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{ice}}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}{\mathit{b}}_{\mathrm{peb}}\mathit{x}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathrm{e}}^{{\mathit{b}}_{\mathrm{peb}}{\mathit{x}}^{\mathrm{\prime }}}\frac{{\mathit{r}}_{\mathrm{snow}}{\mathrm{ℳ}}_{\mathrm{eq}}\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}}{{\mathit{D}}_{\mathrm{p}}}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\\ & \mathrm{=}& \end{array}$(C.10)where we used Eq. (C.6) for ℳ_eq, put p = 0⁸ and defined ε = b_peb/a_eq, ε_gas = b_gas/a_eq and c_ℳ = − ℳ_tot,icer_snow/ Σ_snowD_p(a_eq − b_gas) to keep the notation concise. Typically, ε_gas ≈ 1.5 /a_eq and c_ℳ are small numbers. In a steady-state viscous disk Σ_snow ≈ Ṁ_ice/ 2πrv_gas and ℳ_tot,ice = − Ṁ_ice/ 2πr_snow. Hence in a viscous disk c_ℳ is related to ε_gas: ${\mathit{c}}_{\mathrm{ℳ}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{gas}}{\mathit{r}}_{\mathrm{snow}}}{{\mathit{D}}_{\mathrm{p}}\mathrm{(}{\mathit{a}}_{\mathrm{eq}}\mathrm{-}{\mathit{b}}_{\mathrm{gas}}\mathrm{)}}\mathrm{=}\frac{{\mathit{\epsilon }}_{\mathrm{gas}}}{\mathrm{1}\mathrm{-}{\mathit{\epsilon }}_{\mathrm{gas}}}\mathrm{\approx }\mathrm{0.1}\mathit{.}$(C.11)In Fig. C.2 we plot Σ_ice for the parameters of the default model (Table 1). The simulated data are given by the thick grey curve while the analytical profiles are in blue – the more precise approximation (with p = 7 / 4; solid) and the p = 0 case in dashed. The analytical profiles fit the simulation data very well. The key parameter for the analytical profiles is the velocity of the pebbles v_peb. Here, we have chosen v_peb to be the pebble velocity at the ice peak r_peak, v_peb ≈ 2.4 m s^-1, which provides a decent fit. The fit is somewhat sensitive to the choice of v_peb; for example, adopting v_peb(r_peak) with the initial value for m_ice (i.e., without deposition of vapor) would overestimate the numerical curve by 20%.

Appendix C.4: Peak values and final tuning

Equation (C.10) has a maximum at (x,Σ_ice) = (x_peak,Σ_peak) where ${\mathit{x}}_{\mathrm{peak}}\mathrm{=}\frac{\log \left[\mathit{\epsilon }\mathrm{+}{\mathit{c}}_{\mathrm{ℳ}}\mathrm{(}\mathit{\epsilon }\mathrm{-}\mathrm{1}\mathrm{)}\right]}{{\mathit{a}}_{\mathrm{eq}}\mathrm{(}\mathit{\epsilon }\mathrm{-}\mathrm{1}\mathrm{)}}$(C.12)or r_peak = (1 + x_peak)r_snow in dimensional units, and ${\mathrm{\Sigma }}_{\mathrm{peak}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\epsilon }}_{\mathrm{gas}}\mathrm{)}\frac{{\mathit{c}}_{\mathrm{ℳ}}\mathrm{+}{\left[{\mathit{c}}_{\mathrm{ℳ}}\mathrm{(}\mathit{\epsilon }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}\mathit{\epsilon }\right]}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\epsilon }\mathrm{)}}}{\mathit{\epsilon }}{\mathrm{\Sigma }}_{\mathrm{snow}}\mathit{.}$(C.13)Since c_ℳ and ε_gas are fixed, the latter expression only depends on ε. Lower ε – meaning: a smaller pebble velocity at r_peak – increases both the width (x_peak) and the magnitude of the ice peak, resulting in a more pronounced, “fatter” ice bump.

All that remains is to find an expression for the pebble velocity at the ice peak v_peb, which in turn depends on the amount of ice those pebbles have accreted. To obtain m_peb = m_core + m_ice we invoke conservation of the pebble number flux ${\mathrm{𝒩}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{sil}\mathit{/}{\mathit{m}}_{\mathrm{core}}}{\mathrm{2}\mathit{\pi r}{\mathit{v}}_{\mathrm{peb}}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{ice}}}{{\mathit{m}}_{\mathrm{ice}}}\mathrm{·}$(C.14)Hence, we can obtain m_ice from Σ_peak, calculate the aerodynamical properties of the pebbles at the peak (i.e., their stopping time) and obtain the pebble velocity v_peb.

With Eq. (C.14) we have a closed system of expressions to obtain the relevant quantities at the ice peak. These can best be computed by an iterative scheme, e.g.,

• 1.

  From (an initial guess for) ε obtain Σ_peak and the position of the ice peak r_peak from the ice flux conservation model, Eq. (C.13) .

• From Σ_peak, obtain the ice mass m_ice at the ice peak from the pebble conservation law, Eq. (C.14).

• 3.

  From m_ice and the gas properties at r = r_peak, compute the stopping time of the pebbles at r_peak. Update the pebble velocity v_peb and its normalized variant ε.

This scheme can easily be extended with collective and mean molecular weight effects (i.e., the complete model). Finally, in order to obtain a better match to our numerical results we adopt two ad-hoc empirical “fixes”. First we reduce Eq. (C.13) by 10%, which accounts for the fact that the p = 0 approximation tends to overestimate the peak. Secondly, we slightly increase ε when it drops below unity. This fix approximately accounts for curvature effects (not included in the model) that become apparent when the ice peak becomes broad (at low ε). Hence we take ε = max(ε ∗ ,ε ∗ ^0.8) where ϵ ∗ = (v_peb,peakr_peak/a₂D_peak) is the uncorrected value.

In Table C.1 a comparison of the model with the numerically-obtained peak parameters is given for a number of model parameters. Generally, the agreement is very good; the error in Σ_peak is at most 20%. Given its crudeness, the analytical model, however, does not always produce a perfect match – especially not when it comes to the profile. For example the low turbulence runs (α = 3 × 10^-4), which imply high gas densities when Ṁ_gas is kept the same, still show a good match to Σ_peak, but the profile corresponding to Eq. (C.10) is very different from the numerical model. The reason for this is that pebbles enter the Stokes drag regime when approaching the ice peak, causing a sharp increase in the stopping time (and ε). The stopping time at the ice peak is in that case not representative for other radii.

All Tables

All Figures