Planetesimal Formation by the Streaming Instability in a Photoevaporating Disk

Many current models seek to produce an over-density of dust grains dense enough to self-gravitate and undergo gravitational collapse (e.g., Safronov 1969; Goldreich & Ward 1973; Johansen et al. 2014). These over-densities may be produced by various hydrodynamic processes including long-lived gaseous vortices (Barge & Sommeria 1995; Meheut et al. 2012), pressure bumps (Johansen et al. 2009, 2011; Simon et al. 2012), or through the run-away accumulation of particles by radial drift, known as the streaming instability (e.g., Youdin & Goodman 2005; Johansen et al. 2007; Bai & Stone 2010). The streaming instability has drawn recent attention because it has proven effective at reaching very high particle densities and forming planetesimals on very short (orbital) timescales. In this work, we consider the formation of planetesimals via the streaming instability.

The way that a solid grain responds to aerodynamic drag is determined by its Stokes number, $\,\mathrm{St}={t}_{{\rm{s}}}\,{{\rm{\Omega }}}_{{\rm{K}}}$, where t_s is the stopping time and Ω_K is the Keplerian frequency. For particles in the Epstein regime (particles smaller than the mean-free path of gas) the stopping time is given by

$$\begin{eqnarray}&&{t}_{{\rm{s}}}=\displaystyle \frac{{\rho }_{{\rm{s}}}a}{\rho \,{c}_{{\rm{s}}}}\sqrt{\displaystyle \frac{\pi }{8}},\end{eqnarray} \tag{ 2 }$$

where ρ_s is the material density, ρ is the local gas density, a is the grain size, and c_s is the local sound speed (we have verified that the particles in our simulations are always in the Epstein regime). Carrera et al. (2015) have shown that the streaming instability can be effective at producing particle clumps for a wide range of Stokes numbers, but grains smaller than St = 0.1 require increasingly large densities for the streaming instability to occur. Their measurement of the threshold metallicity as a function of the Stokes number will be used in this work to incorporate planetesimal formation into a global model of an evolving protoplanetary disk. Our prescription for forming planetesimals will be discussed further in Section 3.2.

The size of the grains is determined by the balance between dust coagulation, collisional fragmentation, and radial drift (e.g., Birnstiel et al. 2012; Estrada et al. 2016). Silicate grains also experience a bouncing barrier (Zsom et al. 2010). In a standard α-disk (Shakura & Sunyaev 1976), turbulence generates viscosity, which is parametrized by a parameter ${\alpha }_{{\rm{v}}}$. In the simplest models, turbulence is assumed to be uniform and the turbulence parameter ${\alpha }_{{\rm{t}}}={\alpha }_{{\rm{v}}}$ (from Equation (1)). The collision speed between small particles is

$$\begin{eqnarray}&&{\rm{\Delta }}{u}_{\mathrm{coll}}\sim {c}_{{\rm{s}}}\sqrt{{\alpha }_{{\rm{t}}}\,\mathrm{St}}.\end{eqnarray} \tag{ 3 }$$

Here, we make a distinction between the midplane turbulence α_t and the global ${\alpha }_{{\rm{v}}}$—which sets the viscous evolution of the disk—since α_t and ${\alpha }_{{\rm{v}}}$ may differ in a real disk. This formula (Equation (3)) is valid for St roughly between 10⁻³ and 1 (e.g., Weidenschilling 1984; Ormel et al. 2008). This means that, in the fragmentation-limited regime, the Stokes number of the largest particles is independent of the gas density,

$$\begin{eqnarray}&&\,{\mathrm{St}}_{\mathrm{frag}}\sim \displaystyle \frac{{u}_{\mathrm{frag}}^{2}}{{\alpha }_{{\rm{t}}}\,{c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 4 }$$

where u_frag is the fragmentation speed (Birnstiel et al. 2009). In other words, as the disk evolves and the gas density drops, any increase in Stokes number (Equation (2)) immediately leads to greater collision speeds (Equation (3)) so that the particles fragment until the Stokes number once again reaches St_frag (Equation (4)). As a canonical example, u_frag = 10 m s⁻¹, c_s = 800 m s⁻¹, and ${\alpha }_{{\rm{v}}}={\alpha }_{{\rm{t}}}={10}^{-2}$ would give St_frag ∼ 0.01, independent of the gas density. The u_frag = 10 m s⁻¹ value is consistent with ice–dust aggregates (Wada et al. 2009; Gundlach & Blum 2015). If most of the mass in the dust is in grains close to this maximum size, the streaming instability requires a dust-to-gas ratio of at least Σ_p/Σ_g ∼ 0.03 (Carrera et al. 2015).

Given that the dust-to-gas ratio in the ISM is closer to ∼0.01 (e.g., Leroy et al. 2013), the streaming instability requires processes that either concentrate dust locally, or that preferentially remove the gas component of the disk while leaving the dust behind. We now discuss some mechanisms that can accomplish the required dust concentration.

Inside a protoplanetary disk, solid particles experience a head wind as the gas component orbits at a sub-Keplerian speed. The aerodynamic drag leads to a loss of angular momentum, causing the solids to drift toward the star. The drift velocity is given by

$$\begin{eqnarray}&&{u}_{\mathrm{drift}}=\displaystyle \frac{-2{\rm{\Delta }}v}{\,\mathrm{St}+\,{\mathrm{St}}^{-1}},\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}v={u}_{\mathrm{kep}}-{u}_{\mathrm{gas}}$ is the head wind speed (Weidenschilling 1977). Particles that are either very small or very large are relatively unaffected by radial drift. The fastest drift speed occurs for St = 1. In some cases, the collision speed Δu_coll is dominated by differential drift. The relative drift speed for small particles ( St ≪ 1) is given by the expression

$$\begin{eqnarray}&&{\rm{\Delta }}{u}_{\mathrm{drift}}=\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{u}_{\mathrm{kep}}}\left|\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}\right|| \,{\mathrm{St}}_{1}-\,{\mathrm{St}}_{2}| .\end{eqnarray} \tag{ 6 }$$

Following the notation of Birnstiel et al. (2012), we write ${\mathrm{St}}_{2}=N\,{\mathrm{St}}_{1}\lt \,{\mathrm{St}}_{1}$. Setting Δu_drift = u_frag one gets the drift-driven fragmentation limit,

$$\begin{eqnarray}&&\,{\mathrm{St}}_{\mathrm{df}}=\displaystyle \frac{{u}_{\mathrm{frag}}{u}_{\mathrm{kep}}}{{c}_{{\rm{s}}}^{2}(1-N)}{\left|\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}\right|}^{-1},\end{eqnarray} \tag{ 7 }$$

where P is the gas pressure and r is the radial coordinate (Birnstiel et al. 2012). This quantity depends weakly on r, and $\,{\mathrm{St}}_{\mathrm{df}}\approx 0.19$ for most of the disk. In some cases, particles drift before they can reach the fragmentation limit. The maximum Stokes number that grains can reach in a drift-dominated region is

$$\begin{eqnarray}&&\,{\mathrm{St}}_{\mathrm{drift}}\sim \displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\displaystyle \frac{{u}_{\mathrm{kep}}^{2}}{{c}_{{\rm{s}}}^{2}}{\left|\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}\right|}^{-1}.\end{eqnarray} \tag{ 8 }$$

In the outer disk, where u_kep is small, St_drift <  St_frag and dust growth is drift-limited. However, for our model parameters, fragmentation by differential drift is the more stringent condition.

The water ice line may also play a role in transferring solids into the inner disk. Beyond the ice line, solids would be in the form of ice–silicate aggregates. When they drift across the ice line, the ice component sublimates, leaving behind the silicates. If silicates have a lower fragmentation speed than ice, or if their growth is limited by a bouncing barrier (Zsom et al. 2010), they would be unable to grow to their former size to start drifting again. This process could result in an additional accumulation of dust grains in the vicinity of the water ice line (Sirono 2011a). However, the viscous accretion speed of gas and particles can limit the degree of pile-up.

Another way to increase the dust-to-gas ratio and trigger the streaming instability is to preferentially remove the gas component while leaving the solids in the disk. A well studied mechanism that can accomplish this is photoevaporation. It occurs when energetic photons, primarily in the form of FUV, EUV, and X-rays, heat gas to the point that it becomes gravitationally unbound from the star.

EUV heating occurs through the direct ionization of hydrogen atoms. The absorption cross-section at the Lyman limit is large, and drops quickly for more energetic photons. The ionized layer is heated to T ∼ 10⁴ K (c_s ∼ 10 km s⁻¹), which can launch a flow at ∼1 au for a star of mass 1 M_⊙. When the photoevporative mass loss rate exceeds the local accretion rate, a gap forms in the disk, and viscous accretion depletes the inner disk to form a hole. The mass loss rate due to EUV depends almost exclusively on the flux at ∼1 au. However, the intrinsic EUV luminosity is poorly determined and the flux incident on the disk is very uncertain because the accretion column, as well as any jets or winds close to the star are extremely optically thick to EUV (Alexander et al. 2014). Once the disk forms an optically thin cavity, direct irradiation onto the inner edge produces an order-of-magnitude increase in the photoevaporation rate and the remaining disk dissipates rapidly (Alexander et al. 2006a, 2006b).

X-ray heating occurs indirectly, through the K-shell ionization of heavy elements (mainly O, C, and Fe), where ejected electrons heat the hydrogen gas. T Tauri stars are bright X-ray sources with a median luminosity of L_X ∼ 10³⁰ erg s⁻¹. The photon energy ranges from ∼0.1–10 keV and gas temperatures from a few 100 K to over 10⁴ K. The softer X-ray photons heat the gas to higher temperatures (T ∼ 1–2 × 10⁴ K) but are also easily absorbed at small column densities. The mass loss profile is wider than EUV, with a peak around 3 au (Alexander et al. 2014). While the X-ray photon flux is better determined than that of EUV, the heating process is much more complex. The main uncertainties come from imprecise knowledge of the soft X-ray flux incident on the disk and assumptions about disk chemistry, which affect gas cooling calculations (Alexander et al. 2014).

FUV photoevaporation is by far the most complex. Heating mainly occurs indirectly through the emission of photo-electrons from small dust grains and polycyclic aromatic hydrocarbon molecules, which in turn collide with hydrogen gas. Because dust grains provide the opacity, FUV heating depends not only on the FUV flux, but also on the abundance of small dust grains. The cooling further depends on gas chemistry, similar to the X-ray case, and FUV mass loss rates are subject to the same uncertainties.

Young stars have active chromospheres and hence powerful high-energy photon fluxes across the X-ray, EUV and FUV bands. In addition, accretion shocks due to flows onto the stellar surface at early epochs of star formation can generate a high FUV flux (e.g., Gullbring et al. 2000; Ingleby et al. 2011; Yang et al. 2012). Gorti et al. (2009) found that this accretion-dominated FUV flux can cause photoevaporation during very early stages of disk evolution. It is unclear whether the accretion component has substantial contributions to EUV and X-rays, and if these photons can penetrate the accretion columns and irradiate the disk (Güdel et al. 2007). Gorti et al. (2009) argue that EUV and soft X-rays are strongly absorbed when the disk accretion rates are high and only become relevant at late stages of evolution, and found that X-rays can be important if they have a predominantly soft spectrum.

Mass loss due to FUV photoevaporation occurs primarily in the inner disk (∼3–10 au) and in the outer disk (≳50 au). FUV photoevaporation in the outer disk is particularly significant for disk evolution because this is where most of the disk mass resides and where temperatures as low as 100 K can unbind the hydrogen gas (Gorti et al. 2015). Thus, FUV photoevaporation truncates the outer edge of the disk even at early epochs. As the disk dissipates and accretion drops, the FUV luminosity is also mainly chromospheric and becomes more comparable to EUV and X-rays. As noted earlier, once an optically thin cavity forms, FUV, EUV and X-rays push the cavity outward even as FUV continues to bring the outer edge inward. The two evaporation fronts meet in the middle, and eventually completely disperse the disk (Gorti et al. 2015).

The dust evolution model that we presented in this section was derived by Birnstiel et al. (2012) under an assumption that the gas is a dynamically dominant medium. Photoevaporation could bring the dust-to-gas ratio above unity, so that the assumptions of the model no longer apply. If the dust-to-gas ratio were to move past unity, we would expect the solid mass to dampen the gas turbulence, which would lead to lower collision speeds and larger particles. However, we expect that the streaming instability prevents the dust-to-gas ratio from ever reaching unity in the first place. We will see in Section 3.2 that the streaming instability has a critical dust-to-gas ratio that is less than unity, and above that value any excess dust efficiently converted into planetesimals.

Although we do not include disk winds in our base model, we will discuss them briefly here. Disk winds are driven by a combination of thermal heating and magnetocentrifugal acceleration along large-scale magnetic fields going through the disk. As such, these winds exert a torque on the disk so they carry a significant portion of the disk's angular momentum and drive much of the accretion. The process was first proposed by Blandford & Payne (1982), but see also the review by Turner et al. (2014) and the recent model by Bai et al. (2016).

Similar to photoevaporation, disk winds should preferentially remove gas rather than solids. Bai (2016) has noted that mass loss predominantly occurs in the outer disk and that the corresponding increase in Σ_p/Σ_g may promote planetesimal formation by the streaming instability.

Model NP is a reproduction of the disk evolution model used by Gorti et al. (2015) and is described here for completeness. At its core, it is the radial 1D time-dependent viscous evolution model as described by Lynden-Bell & Pringle (1974), with an additional sink term that describes the mass loss due to photoevaporation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \,{{\rm{\Sigma }}}_{{\rm{g}}}}{\partial t}=\displaystyle \frac{3}{r}\displaystyle \frac{\partial }{\partial r}\left[\sqrt{r}\displaystyle \frac{\partial }{\partial r}(\nu \,{{\rm{\Sigma }}}_{{\rm{g}}}\sqrt{r})\right]-{\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}(r,t),\end{eqnarray} \tag{ 9 }$$

where r is the radial coordinate, ν is the viscosity, Σ_g is the gas surface density, and ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ is the instantaneous photoevaporation rate. The viscosity term ν follows the usual α-model of Shakura & Sunyaev (1976),

$$\begin{eqnarray}&&\nu ={\alpha }_{{\rm{v}}}{c}_{{\rm{s}}}^{2}/{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 10 }$$

where c_s is the local sound speed and Ω_K is the Keplerian frequency. As discussed earlier, the viscous ${\alpha }_{{\rm{v}}}$ may differ from the midplane turbulence α_t. We use the instantaneous accretion rate onto the star to compute the accretion component of the FUV flux, which we add to the chromospheric flux. We assume that the accretion shock has a blackbody spectrum with a temperature of 9000 K (Calvet & Gullbring 1998) and compute the FUV spectrum accordingly. At that temperature, around 4% of the accretion luminosity is in the FUV band (91.2–200 nm), so we estimate

$$\begin{eqnarray}&&{L}_{\mathrm{FUV},\mathrm{acc}}\approx 0.04\left(\displaystyle \frac{0.8\,G\,{M}_{\star }\,\dot{M}}{{R}_{\star }}\right),\end{eqnarray} \tag{ 11 }$$

where M_⋆ and R_⋆ are the stellar mass and radius, and $\dot{M}$ is the instantaneous accretion rate obtained from the solution to Equation (9). We assume that the chromospheric luminosity in FUV and EUV is of the same magnitude as the X-ray luminosity, which is L_X ∼ 10³⁰ erg s⁻¹ for M_⋆ = 1 M_⊙ (Flaccomio et al. 2003).

To compute ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$, at each time step, we build a 1+1D model to resolve the vertical structure of the disk. Photons from the star irradiate the disk surface to heat the dust and gas. Thermal and hydrostatic equilibrium yield the disk structure including the gas temperature T and number density of gas and dust as a function of position (r, z). Following Gorti et al. (2015), we use the analytical approximations of Adams et al. (2004) to compute the evaporation rate ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ at each point (r, z). We consider the photoevaporative flow to be launched from the height z, where ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ is greatest. We refer the reader to Gorti et al. (2015) for further details.

The dust surface density is evolved independently of the gas component. The dust component is divided into 10 logarithmically spaced size bins ranging from 5 nm to 1 cm. The initial total dust-to-gas ratio is 0.01 across the entire disk. The time evolution of each dust bin ${{\rm{\Sigma }}}_{{\rm{p}}}^{i}$ with particle size aⁱ is given by the rate of radial mass advection and diffusion. The diffusive flux itself is proportional to the gradient in the mass density ratio of solids to gas. This gives the relation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{p}}}^{i}}{\partial t}=\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left[r\,{{\rm{\Sigma }}}_{{\rm{g}}}{D}^{i}\displaystyle \frac{\partial }{\partial r}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}}}^{i}}{{{\rm{\Sigma }}}_{{\rm{g}}}}-r{{\rm{\Sigma }}}_{{\rm{p}}}^{i}{u}^{i}\right],\end{eqnarray} \tag{ 12 }$$

where uⁱ is the total radial velocity of the dust due to gas drag and the radial pressure gradient, ${D}^{i}=\nu /(1+\,{\mathrm{St}}_{i}^{2})$ is the dust diffusivity, and St is the Stokes number. The evolution of grain sizes follows the prescription in Section 5.2 of Birnstiel et al. (2011), with the effects of radial drift added from Birnstiel et al. (2012). At every given r and time, grains could drift before they grow (Equation (8)), and drift could result in fragmentation if drift velocities are higher than the threshold (Equation (7)). Following the notation of Birnstiel et al. (2012), we call these limits a_drift and a_df.

We then compute

$$\begin{eqnarray}&&{a}_{1}=\min ({a}_{\mathrm{frag}},{a}_{\mathrm{df}},{a}_{\mathrm{drift}}).\end{eqnarray} \tag{ 13 }$$

Once the maximum particle size a₁ is computed, the dust mass is distributed across the size bins following the prescription of Birnstiel et al. (2011). All of our models follow the same basic prescription. Most models use fragmentation speed u_frag = 10 m s⁻¹ and α_t = 10⁻², but Model IU uses u_frag = 1 m s⁻¹, and Models IA and IAB use α_t = 10⁻⁴. In addition, Models IB and IAB add a maximum size limit due to the bouncing barrier to Equation (13).

In Model SI and later, we instantaneously convert dust into planetesimals whenever the conditions for the streaming instability are met. This is a reasonable approximation because the streaming instability and gravitational collapse of particle clouds occurs on the timescale of just a few orbits (e.g., Johansen et al. 2007).

We do not track the evolution of planetesimals after they form, but we record the planetesimal mass produced. To trigger the streaming instability, we have two requirements. The first requirement is that ${\rho }_{{\rm{p}}}/{\rho }_{{\rm{g}}}\gt 1$ at the midplane, as this threshold marks the transition to high analytical growth rates by the streaming instability (Youdin & Goodman 2005). The value of ρ_p/ρ_g is given by the particle size and diffusion coefficient δ,

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{{\rm{p}}}}{{\rho }_{{\rm{g}}}}=Z\sqrt{\displaystyle \frac{\,\mathrm{St}+\delta }{\delta }},\end{eqnarray} \tag{ 14 }$$

where Z = Σ_p/Σ_g. In the case of turbulence driven by the magnetorotational instability, δ ≈ α_t, but the value may be different if, for example, accretion is driven by disk winds (e.g., Johansen et al. 2014; Johansen & Klahr 2005; Bai & Stone 2013). The condition ρ_p/ρ_g > 1 can be thought of as a minimum Z or a maximum value for the midplane turbulence α_t,

$$\begin{eqnarray}&&{\alpha }_{{\rm{t}}}\lt \displaystyle \frac{{Z}^{2}\,\mathrm{St}}{1+{Z}^{2}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&Z\gt {Z}_{1}\equiv \sqrt{\displaystyle \frac{{\alpha }_{t}}{\,\mathrm{St}+{\alpha }_{{\rm{t}}}}}.\end{eqnarray} \tag{ 16 }$$

In addition, Carrera et al. (2015) found a minimum threshold in Z for dense particle filaments to form by the streaming instability; this threshold occurs even when α_t is negligible. We computed a least-squares fit for the results of Carrera et al. (2015) and obtained the condition,

$$\begin{eqnarray}&&Z\gt {Z}_{2}\equiv {10}^{-1.86+0.3{(0.98+{\mathrm{log}}_{10}\mathrm{St})}^{2}}.\end{eqnarray} \tag{ 17 }$$

In line with the results of Carrera et al. (2015), we also require St ≥ 0.003. More recent work by Yang et al. (2016) shows that the streaming instability can be effective for even smaller particles. We will incorporate those new results in a future work. Figure 2 shows both criteria together. Let Z₁ and Z₂ be the right-hand sides of Equations (16) and (17), respectively. Together, the conditions for the streaming instability become

$$\begin{eqnarray}&&Z\gt {Z}_{\mathrm{crit}}=\max ({Z}_{1},{Z}_{2}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\,\mathrm{St}\gt 0.003.\end{eqnarray} \tag{ 19 }$$

The simulations conducted by Carrera et al. (2015) always had a uniform particle size. In our simulations, we have to generalize this result to a situation where there is a range of particle size bins. To do this, we start by computing the value of Z for the largest size bin. If Z ≥ Z_crit, we remove 90% of the mass in that bin. If Z < Z_crit, we compute Z for the top two bins, and re-compute Z_crit as if the entire dust mass was in the lower mass bin. We repeat this process until planetesimals form, or we reach the minimum particle size of St = 0.003. This is a conservative choice because it discounts the effect of larger grains, even though they tend to dominate the streaming instability (Bai & Stone 2010).

We assume that planetesimals always form when the streaming instability is active. We believe that this assumption is justified because the streaming instability gathers solids into long-lived stable filaments. If the filaments do not form planetesimals immediately, they will become effective traps for other solid particles drifting radially from the outer disk. Therefore, the filament density will grow on the radial drift timescale. The filament density will drastically increase when photoevaporation dissipates the last remaining gas from the disk. In Section 5.3, we discuss the physics of gravitational collapse in greater detail.

Model I and later include the water ice line. Near the ice line, we must model the sublimation and condensation of water ice onto solid grains. We use a pure thermal desorption/absorption process, meaning that ice forms whenever the partial pressure of water vapor is greater than the saturation pressure. We ignore photodesorption (e.g., from UV) because we are only interested in the formation of ices on the midplane, which is shielded from high-energy photons. The saturation pressure of water is given by

$$\begin{eqnarray}&&{P}_{\mathrm{sat}}=2.53\times {10}^{13}\exp \left(\displaystyle \frac{-6070\,{\rm{K}}}{T}\right)\,\mathrm{dyne}\ \,{\mathrm{cm}}^{-2}\end{eqnarray} \tag{ 20 }$$

where T is the local gas temperature (Leger et al. 1985). The partial pressure of water is

$$\begin{eqnarray}&&{P}_{{{\rm{H}}}_{2}{\rm{O}}}={n}_{{{\rm{H}}}_{2}{\rm{O}}}\,{k}_{{\rm{B}}}T={n}_{{{\rm{H}}}_{2}}{X}_{{{\rm{H}}}_{2}{\rm{O}}}\,{k}_{{\rm{B}}}T,\end{eqnarray} \tag{ 21 }$$

where ${n}_{{{\rm{H}}}_{2}}$ and ${n}_{{{\rm{H}}}_{2}{\rm{O}}}$ are the number densities of hydrogen and water, $\,{k}_{{\rm{B}}}$ is the Boltzmann constant and ${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the water abundance. Ice forms whenever ${P}_{{{\rm{H}}}_{2}{\rm{O}}}\gt {P}_{\mathrm{sat}}$, which simplifies to

$$\begin{eqnarray}&&{n}_{{{\rm{H}}}_{2}}\gt \displaystyle \frac{2.53\times {10}^{13}\ \,\mathrm{dyne}\ \,{\mathrm{cm}}^{-2}}{{X}_{{{\rm{H}}}_{2}{\rm{O}}}\ \,{k}_{{\rm{B}}}T}\exp \left(\displaystyle \frac{-6070\,{\rm{K}}}{T}\right).\end{eqnarray} \tag{ 22 }$$

The run begins with ${X}_{{{\rm{H}}}_{2}{\rm{O}}}=3\times {10}^{-4}$ everywhere, which corresponds to a water mass fraction of 0.005. Note that the CO ice line can also be important for disks, but is ignored here. Our chosen water abundance therefore reflects all the oxygen being bound in water ice and slightly underestimates the total solid mass in the outer disk beyond the putative CO ice line. For simplicity, we further assume that the water is either all vapor or all ice, and we use the midplane temperature with mean number density $\langle {n}_{{{\rm{H}}}_{2}}\rangle =\langle {\rho }_{{{\rm{H}}}_{2}}\rangle /{m}_{{{\rm{H}}}_{2}}$ integrated over one scale height.

When $\langle {n}_{{{\rm{H}}}_{2}}\rangle$ satisfies Equation (22), all the water is added to the dust component uniformly, with constant dm/m across all grain size bins. This differs from the approach of Ros & Johansen (2013), who kept da constant. This latter choice is consistent with pure vapor condensation, but in a disk where dust particles have frequent collisions that result in coagulation and fragmentation, the water component could be quickly distributed across bin sizes.

We keep track separately of Σ_ice and Σ_vapor in order to determine the water fraction in dust or gas. Inside the ice line, Σ_ice is locally released into Σ_vapor, which advects and diffuses like the rest of the disk gas.

Inside the water ice line, the solid component is made entirely of silicate grains, which may experience a bouncing barrier (Güttler et al. 2010; Zsom et al. 2010). Models IB and IAB implement this bouncing barrier. To do this, we manually fit a straight line through the plot in Figure 11 of Güttler et al. (2010). We use the results for collisions between equal-mass compact grains. This is the most conservative choice we can make because the other choices would produce larger grains that more easily participate in the streaming instability. The choice of compact grains is also consistent with our opacity calculation. Our manual fit to the mass of particles at the bouncing barrier m_b, is

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{m}_{{\rm{b}}}}{{\rm{g}}}\right)=-2.08{\mathrm{log}}_{10}\left(\displaystyle \frac{{\rm{\Delta }}{u}_{\mathrm{coll}}}{\mathrm{cm}\,{{\rm{s}}}^{-1}}\right)-7.33.\end{eqnarray} \tag{ 23 }$$

Inside the ice line, we use Equation (23) to compute the corresponding grain size at the bouncing barrier a_b and replace Equation (13) with the minimum of a₁ and a_b.

Most of our models follow a simple α-disk prescription (Shakura & Sunyaev 1976) in which ${\alpha }_{{\rm{t}}}={\alpha }_{{\rm{v}}}={10}^{-2}$ everywhere in the disk. However, on the midplane and in the outer disk, the actual level of turbulence is likely to be significantly lower. This is important for planetesimal formation because solids are mostly on the midplane and the collision speed between solid particles is proportional to $\sqrt{{\alpha }_{{\rm{t}}}}$ (Equation (3)), and because the midplane turbulence sets the amount of sedimentation (set δ = α_t in Equation (14)).

First, from around 0.2–1 to 30–100 au, protoplanetary disks are thought to experience layered accretion, in which the midplane is inside an MRI-inactive "dead zone" that experiences very little turbulence (Gammie 1996; Fromang et al. 2002). Simulations of the magnetorotational instability by Oishi et al. (2007) give an α_t of 10⁻⁵ inside the dead zone. Alternatively, simulations of the vertical shear instability give α_t = 10⁻⁴ on the midplane (Stoll & Kley 2014). While there are no strong theoretical constraints beyond 100 au, ALMA CO observations indicate a 3σ upper limit of α < 10⁻³ for the outer regions (R > 30 au) of the disk around HD 163296 (Flaherty et al. 2015). Recent ALMA observations of the TW Hya disk (Teague et al. 2016) suggest a value of α closer to our adopted default value of 0.01.

Our choice thus far of α_t = 10⁻² in the midplane may therefore be high, and it is possible that true particle collision speeds are in fact lower than calculated. On the other hand, disk lifetimes indicate an effective viscous ${\alpha }_{{\rm{v}}}\sim {10}^{-2}$, and lower values yield disk lifetimes far longer than observed (Gorti et al. 2015). In Models IA and IAB, we separate ${\alpha }_{{\rm{v}}}$ and α_t: ${\alpha }_{{\rm{v}}}$ is used to compute the viscous evolution of the disk, and α_t is used to compute the collision speed between solid particles. We kept the viscous evolution at ${\alpha }_{{\rm{v}}}={10}^{-2}$ as in our previous runs, and we chose a conservative α_t = 10⁻⁴. This lowers collision speeds by a factor of 10 while maintaining the nominal accretion rate of the gas and giving reasonable disk lifetimes.

Figure 3 shows how the formation of planetesimals affects the evolution of a protoplanetary disk. Models NP and SI have the same disk parameters, but Model SI introduces planetesimal formation by the streaming instability as described in Section 3.2. The main consequence of planetesimal formation is the depletion of dust. Model NP leaves behind 70 M_⊕ of dust. This much dust is inconsistent with observations of debris belts (Wyatt 2008). In Model SI, only 0.002 M_⊕ of dust remains because the rest is converted into planetesimals. In reality, once the protoplanetary disk evolves into a debris disk, the remaining dust is quickly removed by either photon pressure or Poynting–Robertson drag, and new dust is re-generated through collisional grinding of surviving planetesimals (Wyatt 2008). We do not attempt to model these processes because the debris phase is beyond the scope of our investigation and not relevant to our goal of understanding planetesimal formation. Therefore, the solid component in Figure 3 stops evolving once the gas dissipates.

Zoom In Zoom Out Reset image size

Figure 4 shows the fate of the dust component in greater detail. Model NP finishes with 23 M_⊕ of 1–100 μm grains, and 47 M_⊕ of sub-μm grains. However, with the inclusion of the streaming instability in Model SI, 76 M_⊕ of dust are converted into planetesimals and only 0.002 M_⊕ of sub-μm grains survive to the end, along with negligible amounts of larger grains.

Zoom In Zoom Out Reset image size

In Model SI, about 24% of the original dust in the disk is converted into planetesimals. However, very few of those planetesimals form in the inner disk. Table 2 and Figure 5 show the final distribution of planetesimals for all the runs. In all cases, most of the planetesimals form beyond 100 au, which indicates the presence of a massive outer debris belt. We note that our debris belt is at larger radii than typically observed (Wyatt 2008). Given that we consider neither migration of planetesimals after formation nor the ensuing stages of planet formation, we cannot meaningfully compare our results to planetary system architectures.

Zoom In Zoom Out Reset image size

Figures 6 and 7 give a general picture of how planetesimal formation occurs in Model SI. Planetesimals always begin to form early (after ∼500,000 years) in the outer disk. Even in the absence of photoevaporation, gas viscously spreads out to larger radii while dust generally drifts in radially. This causes an enhancement in the local dust-to-gas ratio in the outer disk. With photoevaporation, gas is preferentially removed from the outer disk, increasing Z above Z_crit and triggering the streaming instability. FUV-driven photoevaporation begins in the outer disk at early times, and this is where planetesimals first form in all of our models. As the disk evolves, its outer radius decreases due to FUV photoevaporation and planetesimal formation moves inward.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Model SI, there is not much mass loss due to photoevaporation in the inner disk at early epochs. Viscous accretion depletes both gas and dust, and there is no preferential removal of gas at small radii to bring Z above Z_crit. At later epochs, when the surface density of gas decreases, the mass accretion rate eventually falls below the photoevaporation rate and a gap opens in the disk (as shown previously by Alexander & Armitage 2009). The disk dissipates rapidly after this, increasing the solids-to-gas mass ratio globally. Model SI shows a last-minute strong burst of planetesimal formation, which leaves behind very little remaining dust. Models I, IU, and IB look very similar to SI, and show all the same features. Models IA and IAB show different evolution and will be discussed in later sections.

It is important to note that the terrestrial planets of the solar system are water poor, and could not have formed from planetesimals that originated beyond 2–3 au (approximate location of the ice line). Therefore, the total planetesimal mass that forms within 3 au is an important benchmark for any planetesimal formation model. We will refer to the inner 3 au as the terrestrial zone. Model SI does not produce enough planetesimals in this region to form the terrestrial planets in the solar system (Table 2).

Model I (Table 1) has the same model parameters as Model SI, but the addition of a water ice line, which probably plays an important role in planetesimal formation. Beyond the ice line, water contributes to the solid component of the disk—Model I starts with 476 M_⊕ of solids, compared to only 313 M_⊕ for Model SI. Water ice maintains a high fragmentation barrier of u_frag ∼ 10 m s⁻¹ (Wada et al. 2009), so that ice-silicate aggregates can achieve larger Stokes numbers than pure silicates. Finally, water ice may also play an important role in transporting solids to the inner disk because the larger Stokes numbers of ice-silicate aggregates are more susceptible to radial drift (Equation (5)). However, when the aggregate crosses the ice line, the ice component will sublimate and leave behind the silicate dust. Since pure silicate aggregates cannot achieve the same Stokes numbers as the ice-rich aggregates, they experience much lower radial drift, so they would tend to accumulate. The greater solid mass, and the increase in the dust-to-gas ratio could contribute to the formation of planetesimals in the terrestrial zone. We explore this possibility in Model I.

Table 2 shows that Model I does produce significantly more planetesimals than Model SI (123 M_⊕ versus 76 M_⊕) as would be expected given the larger solid mass. However, most of the gains remain in the outer disk. Figure 5 shows that the distribution of planetesimals in Model I is basically the same as Model SI. Water ice alone is not enough to produce a large planetesimal mass in the terrestrial zone.

As noted earlier, the u_frag = 10 m s⁻¹ limit is consistent with the ice–dust aggregates that would be found beyond the ice line (Wada et al. 2009; Gundlach & Blum 2015). However, given the uncertainty in disk dust properties, we next test how our results would change if ices fragment more easily than expected. Model IU is a copy of Model I, with the fragmentation speed reduced to u_frag = 1 m s⁻¹. This change increases the total planetesimal mass (153 M_⊕ versus 123 M_⊕). The reason for this is that easier fragmentation leads to smaller grains that do not drift readily. Table 2 shows that the gain in planetesimals occurs entirely beyond 100 au. Closer to the star, the lack of radial drift actually depletes the planetesimals. In the terrestrial zone (inside 3 au), the planetesimal mass drops by 50%, as more of the solids turn into planetesimals further out.

Figure 5 again shows that the distribution of planetesimals is essentially the same as that of Models SI and I. There are more planetesimals in the outer disk, especially at 1000 au, and fewer planetesimals in the inner disk, inside 3 au. These results show that the fragmentation speed does not drastically change our results, though a low fragmentation barrier may exacerbate the difficulties in forming dry planetesimals in the terrestrial zone.

Model IB is a copy of Model I (both have a water ice line and u_frag = 10 m s⁻¹) with the addition of the bouncing barrier for silicates inside the ice line (Table 1). Our implementation of the bouncing barrier (Section 3.4) is derived from Güttler et al. (2010). The bouncing barrier affects pure silicate grains and keeps them from reaching the fragmentation limit. Figure 8 shows the mass-weighted mean particle size for Models I and IB at three different points in time. There are three key features in this plot.

• 1.  
  In the outer disk, Models I and IB look identical. Inside the ice line (approximately 3 au) the bouncing-limited grains of Model IB are significantly smaller. The transition in particle size is smooth and goes between 3 and ∼10 au because the small bouncing-limited grains can diffuse outward.
• 2.  
  The jump in particle size that occurs in the outer disk (at around 20 au) marks the transition from drift-limited growth (outer disk) to fragmentation-limited growth (middle and inner disk).
• 3.  
  As time passes and the disk dissipates, particles become more strongly coupled to the gas and their collision speeds increase (Equation (3)) so the particles must become smaller.

Zoom In Zoom Out Reset image size

Table 2 and Figure 5 show the effect of the bouncing barrier. Not surprisingly, Model IB looks extremely similar to Model I beyond 3 au. Inside 3 au, the planetesimal mass drops sharply because the smaller grains produced by the bouncing barrier do not easily participate in the streaming instability and instead they accrete onto the star. Model IB also leaves more dust after the gas disk has completely dispersed inside 3 au than Model I (still less than 10⁻³ M_⊕). When the disk is about to disappear (∼3 Myr), the dust mass inside 2 au is completely dominated by sub-micron grains. We note that collisions between planetesimals may generate an additional (debris) dust component, which we do not consider in these evolution models.

All the models presented so far have the turbulence parameter fixed at ${\alpha }_{{\rm{v}}}={\alpha }_{{\rm{t}}}={10}^{-2}$. As we discussed in Section 3.5, turbulence in the midplane may be better represented by low values of α_t = 10⁻⁴. So in Model IA (Table 1), we use two parameters: ${\alpha }_{{\rm{v}}}={10}^{-2}$ to compute the viscous evolution of the disk, and α_t = 10⁻⁴ to compute the collision speeds between solid particles. The fragmentation-limited Stokes number scales as 1/α_t (Equation (4)). In Model IAB, the grain sizes are much larger than in Model I, but note that in IA the particle growth is drift-limited beyond ∼1 au.

Naively, one might expect that larger grains would lead to more planetesimal formation, but that does not seem to be the case. Model IA actually produces the fewest planetesimals (Table 2). This is because of the effects of larger grains on the mass loss rate due to FUV photoevaporation. As discussed in Gorti et al. (2015), changing the cross-sectional grain area leads to two effects: (i) a reduction in FUV attenuation and deeper penetration of photons resulting in a higher flow density, and (ii) reduction in the FUV heating due to small grains resulting in lower flow temperatures. These oppose each other (${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}\propto {\rho }_{{\rm{g}}}\sqrt{T}$) and hence the photoevaporation rate only depends weakly on grain size with the temperature effect being marginally stronger. Larger grain sizes therefore result in slightly lower mass loss rates, and hence longer disk survival times (also see Figure 4 of Gorti et al. 2015). This is seen in Figure 7 where the reduction in the FUV photoevaporation rate is seen to prolong the disk lifetime by 0.82 Myr (based on the time when 99.9% of the disk gas has disappeared). A long disk lifetime is not inherently detrimental to planetesimal formation, but with less FUV photoevaporation, the dust-to-gas ratio cannot easily reach the value required by the streaming instability across most of the disk.

Finally, Model IAB (Table 1) includes the water ice line, the bouncing barrier of Model IB, and the reduced midplane turbulence of Model IA. This is our most complete model. Its most salient feature is that it actually produces a large population of planetesimals in the terrestrial zone. While other models can barely reach 0.4 M_⊕ inside 3 au, Model IAB produces 8.4 M_⊕—easily enough to reproduce the terrestrial planets in the solar system (Table 2). In this section, we investigate how the combination of the bouncing barrier and low midplane turbulence can increase the planetesimal mass inside the terrestrial zone by an order of magnitude compared to Model I.

We begin with a look at the spacetime diagram of Model IAB, shown in Figure 9. For the first few million years, the results look nearly identical to Model IA (shown in Figure 7). However, the reduced midplane turbulence in this case also results in larger grains and slightly longer disk lifetimes. Larger grains also result in more drift and increase the solids-to-gas mass fraction slightly. The addition of the bouncing barrier in the inner disk makes grains smaller in Model IAB and mitigates the loss of dust due to radial drift. The smaller grains also result in higher FUV photoevaporation rates in the inner disk compared to the other models. All these effects combine to produce a different evolutionary sequence in Model IAB. After 2.7 Myr the outer disk suddenly stops producing planetesimals for 60,000 years. Then planetesimals start to form again at a slower rate, and farther out. This break occurs just as the inner disk starts forming planetesimals. In the inner disk, planetesimal formation follows a different pattern from previous models. Planetesimals begin to form inside a narrow band that gradually expands both inward and outward.

Zoom In Zoom Out Reset image size

Figure 10 shows what is happening in Model IAB 100,000 years before the sudden change in the disk. All the other runs are also included for comparison. The figure shows the surface density profile of the gas, dust, and planetesimal components of each model. What sets Model IAB apart from the others is that it starts forming an inner cavity already at 2.6 Myr, while there are still 17 M_⊕ of dust remaining in the disk. As the disk evolves, the bouncing barrier produces small grains and increases FUV photoevaporation in the inner disk. A gap begins to form where the photoevaporation rate exceeds the accretion rate. This decreases the gas surface density in the inner disk. When the inner disk depletes, accretion onto the star decreases sharply (since ${\dot{M}}_{\mathrm{acc}}\propto {{\rm{\Sigma }}}_{{\rm{g}}}$). This has two key consequences.

• 1.  
  When the accretion rate decreases, the accretion-generated FUV flux also drops dramatically. This effectively shuts down FUV photoevaporation in the outer disk. When that happens, the disk begins to expand viscously. At some point, the outer edge becomes distant enough that the lower FUV emission can once again unbind the gas molecules. This causes the hiatus in planetesimal formation that we observe in Figure 9.
• 2.  
  In the inner disk, loss of gas raises the dust-to-gas ratio. Photoevaporation by FUV and X-rays continues and the solid grains eventually reach a dust-to-gas ratio and Stokes number that allow the streaming instability to form planetesimals. Admittedly, the streaming instability has not been studied well in environments with extreme dust-to-gas ratios, but if the streaming instability did not occur, the dust would form an ever denser, ever thinner disk that would eventually become gravitationally unstable in some Toomre-like instability.

Zoom In Zoom Out Reset image size

Figure 9 also gives us a clue of how the final stages of the disk evolution occur. The planetesimals in the spacetime diagram very roughly trace the edges of the disk. After 2.6 Myr, the inner cavity forms and starts to expand outward. After 2.7 Myr, FUV photoevaporation is re-ignited in the outer disk, and the outer edge once gain starts to move inward. As the two disk edges evolve, they leave behind a trail of planetesimals. After 4 Myr the two fronts meet, and the disk is fully dissipated.

Figure 11 helps complete the picture in the inner disk. The top plot of the figure shows the gas density profile in Model IAB after 2.7 Myr—just at the moment when the inner cavity has fully formed. The bottom plot shows the mass distribution of planetesimals at the end of the simulation. This figure makes it clear that the planetesimals form just outside the gap. This is what we would expect if the underlying mechanism is that the gap prevents the accretion of solids onto the star.

Zoom In Zoom Out Reset image size

Lastly, what makes Model IAB produce an early gap while the other models do not? The answer seems to be that the bouncing barrier promotes FUV photoevaporation in the inner disk (by creating very small grains) while in the outer disk, the larger particles inhibit FUV photoevaporation and cause the disk to last longer.

We note that all the results presented here are sensitive to assumptions made on viscous disk evolution. We have chosen a simple α-parametrization of disk viscosity, and moreover assume that α is constant in radius and time as the disk evolves. These assumptions may not be realistic, but a more accurate prescription of disk viscosity is still not available, though considerable progress is being made in this direction (e.g., Simon et al. 2015; Bai et al. 2016). Disk lifetimes due to a combination of viscous evolution and photoevaporation also depend on the choice of the numerical value of α (see Gorti et al. 2015, for a more in-depth discussion), a value of 0.01 best reproduces most observed characteristics, such as the overall disk lifetime, stellar mass accretion rate and accretion luminosity. Viscous evolution sets the basic timescales in the disk (more so at earlier epochs when photoevaporation rates are comparatively low), and all physical processes modeled here are also subject to resultant uncertainties.

We conducted simulations with lower viscosity (${\alpha }_{{\rm{v}}}={10}^{-3}$). As shown in Gorti et al. (2015), the lower accretion luminosity and viscous spreading result in longer disk evolution times. Our key result (that most planetesimals form in the outer disk) remained unchanged. FUV photoevaporation occurs everywhere in the disk. Its presence is most easily seen in the outer disk because the gas removed is not replenished; the dust-to-gas ratio increases and planetesimal formation occurs. Elsewhere, the gas that photoevaporates is replaced by gas that viscously accretes from outside, so the dust-to-gas ratio remains low.

In models I and IA, the pile-up of solids inside the ice line is frustrated by the viscous accretion that drags these solids into the star. This is an artifact of our 1D disk model, which forces the entire gas column to accrete at the same rate, set by ${\alpha }_{{\rm{v}}}={10}^{-2}$. In model IA, we can capture the low collision speeds from the α_t = 10⁻⁴ on the midplane, but the gas component is still forced to evolve with ${\alpha }_{{\rm{v}}}={10}^{-2}$. This means the model greatly overestimates the gas accretion rate on the midplane. Outside the ice line, this is not much of an issue because the radial speed of the solids is dominated by the radial drift relative to the gas (Equation (5)). However, inside the ice line, where the bouncing barrier keeps solids closely coupled to the gas, our models greatly overestimate the radial speed of solids. In fact, a 2D model by Ciesla (2007) shows that there is even a region of outward gas flow around the midplane, and an outward-transport region where the viscous flow pushes solids away from the star.

With a more realistic model that could capture the laminar structure of the disk, we would expect to see a much larger pile-up of solids close to the ice line, and a much larger planetesimal mass in the terrestrial zone.

In this section, we take a closer look at the formal conditions for particle collapse. As noted in the Introduction, the streaming instability belongs to a family of models that aim to produce over-densities of solid particles to the point of initiating gravitational collapse. In the case of the streaming instability, these over-densities take the form of elongated filaments that stretch along the azimuthal direction. The formal (but approximate) condition for gravitational collapse is that the midplane density of solid particles ρ_p must reach the Roche density,

$$\begin{eqnarray}&&{\rho }_{{\rm{R}}}=\displaystyle \frac{9{{\rm{\Omega }}}^{2}}{4\pi G}=\displaystyle \frac{9}{4\pi }\displaystyle \frac{{M}_{\star }}{{r}^{3}},\end{eqnarray} \tag{ 24 }$$

where Ω is the Keplerian frequency, G is the gravitational constant, M_⋆ is the mass of the central star, and r is the orbital distance. When ρ_p > ρ_R, the particle self-gravity overwhelms Keplerian shear, leading to collapse (Goldreich & Ward 1973; Johansen et al. 2014). This means that, even when the streaming instability effectively forms a stable filament of particles, that filament will not produce any planetesimals unless ρ_p > ρ_R. Since the streaming instability is largely a scale-free process, the particle density inside the filaments is proportional to the gas density at the time that the streaming instability occurs,

$$\begin{eqnarray}&&{\rho }_{{\rm{p}},\mathrm{SI}}=\epsilon \ {\rho }_{{\rm{g}}},\end{eqnarray} \tag{ 25 }$$

where ${\rho }_{{\rm{g}}}$ is the gas density at the midplane, and ${\rho }_{{\rm{p}},\mathrm{SI}}$ is the particle density produced by the streaming instability inside a filament. The value of probably depends on the Stokes number, with smaller particles forming less dense filaments (Carrera et al. 2015). As we discussed in Section 3.2, the streaming instability can only occur when ρ_p/ρ_g > 1 on the midplane. The streaming instability then proceeds to concentrate particles into denser filaments, so that $\epsilon ={\rho }_{{\rm{p}},\mathrm{SI}}/{\rho }_{{\rm{g}}}\gg 1$ inside the filaments. However, the value of has only been measured precisely for St = 0.3 particles, which have going locally up to at least 10⁴ (Johansen et al. 2015).

As mentioned in the discussion, the streaming instability occurs whenever Σ_p/Σ_g > Z_crit, where Z_crit is a function of the particle Stokes number. From the point of view of the streaming instability, it makes no difference if Z_crit is achieved by removing gas or by adding solids. However, Equation (25) shows that it is easier to reach ${\rho }_{{\rm{p}},\mathrm{SI}}\gt {\rho }_{{\rm{R}}}$ if the streaming instability can be triggered at a higher gas density.

To proceed with this investigation, we consider a toy model in which we calculate ρ_g, Z, St, and Z_crit( St) at 2.5 Myr (i.e., shortly before the inner disk begins to photoevaporate) and ignore any subsequent change in St or the solid density. In this toy model, St and ρ_p are fixed, and the only way to trigger the streaming instability is to increase Z by removing gas. In other words, the gas component has to be depleted by a factor of Z/Z_crit for the streaming instability to occur,

$$\begin{eqnarray}&&{\rho }_{{\rm{g}},\mathrm{crit}}={\rho }_{{\rm{g}},2.5\mathrm{Myr}}\,\displaystyle \frac{{Z}_{2.5\mathrm{Myr}}}{{Z}_{\mathrm{crit}}},\end{eqnarray} \tag{ 26 }$$

where ρ_g,crit is the (critical) midplane density when the streaming instability begins. We then compute ${\rho }_{{\rm{p}},\mathrm{SI}}$ assuming  = 10⁴,

$$\begin{eqnarray}&&{\rho }_{{\rm{p}},\mathrm{SI}}={10}^{4}\,{\rho }_{{\rm{g}},\mathrm{crit}}.\end{eqnarray} \tag{ 27 }$$

Except for the value of , which is not well understood for the wide range in Stokes numbers considered here, this toy model is conservative. Larger Stokes numbers imply lower Z_crit and larger ${\rho }_{{\rm{p}},\mathrm{SI}}.$

The implications of this toy model are shown in Figure 12. The figure shows a snapshot of three sample runs at 2.5 Myr. The plot on the top-left shows the mass-weighted mean particle size,

$$\begin{eqnarray}&&\langle a\rangle =\displaystyle \sum _{i=1}^{N}\,\displaystyle \frac{{a}_{i}{{\rm{\Sigma }}}_{{\rm{p}},i}}{{{\rm{\Sigma }}}_{{\rm{p}}}},\end{eqnarray} \tag{ 28 }$$

where a_i and ${{\rm{\Sigma }}}_{{\rm{p}},i}$ are the size and surface density of particles in the ith particle bin, and Σ_p is the total particle surface density. The plot on the top-right shows the Stokes number for particles of size $\langle a\rangle$. The dashed line marks the limit St = 0.003, which is the smallest particle size that has been observed to participate in the streaming instability (Carrera et al. 2015). This plot seems to imply that models SI and IB cannot form particle filaments inside 3 au. We remind the reader that the onset of photoevaporation in the inner disk leads to an increase in the particle Stokes numbers inside 3 au. We chose to exclude this fact in this section to keep the toy model simple. In the bottom two plots, we only show results in places where St ≥ 0.003 before photoevaporation.

Zoom In Zoom Out Reset image size

The bottom-left plot of Figure 12 shows the critical dust-to-gas ratio Z_crit assuming a uniform particle size $\langle a\rangle$. Finally, we show ${\rho }_{{\rm{p}},\mathrm{SI}}/{\rho }_{{\rm{R}}}$ in the bottom-right plot, following Equations (24) and (27). We also mark the critical value ${\rho }_{{\rm{p}},\mathrm{SI}}/{\rho }_{{\rm{R}}}=1$, above which the filaments produced by the streaming instability should collapse into planetesimals. The plot shows that model IA, benefiting from large particles, can easily form planetesimals well inside the terrestrial zone. The other models cannot form planetesimals unless the Stokes numbers increase.

One feature of the bouncing barrier that we did not implement is that it is a "soft" barrier—when gas dissipates, the particles do not fragment into smaller ones, so their St increases. As ${\rho }_{{\rm{g}}}\to 0$, St approaches the fragmentation barrier. To conclude this section, we list three more processes that may significantly increase ρ_p inside a filament and have not been adequately explored in the literature.

• 1.  
  Inside the filament, orbits are strongly Keplerian and collision speeds between particles may be much smaller than ${c}_{{\rm{s}}}\sqrt{\alpha \,\mathrm{St}}$ (Equation (3)). If so, the Stokes number will no longer be fixed by Equation (4), and coagulation will proceed to larger Stokes numbers. As particles grow, collision speeds would continue to drop. See Carrera et al. (2015) for an earlier exploration of this concept.
• 2.  
  Even when a stable filament fails to reach the Roche density, ρ_p will continue to grow on the radial drift timescale as solid particles in the outer disk migrate inward and become trapped when they enter the filament. This process will only work if the drift timescale is shorter than the evaporation timescale, which is typically the case.
• 3.  
  Finally, we do not know how a filament produced by the streaming instability will respond to the removal of gas by photoevaporation. At some point, the gas component should not be able to support the full vertical height of the filament and ρ_p should begin to increase. As ${\rho }_{{\rm{g}}}\to 0$, ${\rho }_{{\rm{p}}}\to \infty$ on the midplane and the filament may experience a 1D instability along the lines of Goldreich & Ward (1973), Sekiya (1998), and Youdin & Shu (2002).

The last two mechanisms in particular suggest that, even if ${\rho }_{{\rm{p}},\mathrm{SI}}\lt {\rho }_{{\rm{R}}}$ initially, ρ_p could grow arbitrarily high and eventually reach the Roche density. In other words, once the streaming instability forms a stable filament, planetesimals may be unavoidable. The formation of filaments of low-St pebbles and the evolution of such filaments in response to gas removal are topics that should be explored more in future simulations.