Planetesimal formation by the streaming instability in a photoevaporating disk

TL;DR

This study presents a comprehensive global model of protoplanetary disk evolution that incorporates streaming instability and photoevaporation, demonstrating efficient planetesimal formation especially beyond 100 au, but highlighting the need for additional processes for inner disk planetesimal formation.

Contribution

It introduces a detailed, global simulation model combining streaming instability with photoevaporation, dust evolution, and turbulence, advancing understanding of planetesimal formation across the disk.

Findings

Massive planetesimal belts form beyond 100 au.

Inner disk forms up to 8 M⊕ of planetesimals within 3 au.

Photoevaporation influences the timing and location of planetesimal formation.

Abstract

Recent years have seen growing interest in the streaming instability as a candidate mechanism to produce planetesimals. However, these investigations have been limited to small-scale simulations. We now present the results of a global protoplanetary disk evolution model that incorporates planetesimal formation by the streaming instability, along with viscous accretion, photoevaporation by EUV, FUV, and X-ray photons, dust evolution, the water ice line, and stratified turbulence. Our simulations produce massive (60-130 $M_\oplus$) planetesimal belts beyond 100 au and up to $\sim 20 M_\oplus$ of planetesimals in the middle regions (3-100 au). Our most comprehensive model forms 8 $M_\oplus$ of planetesimals inside 3 au, where they can give rise to terrestrial planets. The planetesimal mass formed in the inner disk depends critically on the timing of the formation of an inner cavity in the disk by high-energy photons. Our results show that the combination of photoevaporation and the streaming instability are efficient at converting the solid component of protoplanetary disks into planetesimals. Our model, however, does not form enough early planetesimals in the inner and middle regions of the disk to give rise to giant planets and super-Earths with gaseous envelopes. Additional processes such as particle pileups and mass loss driven by MHD winds may be needed to drive the formation of early planetesimal generations in the planet forming regions of protoplanetary disks.

Figures

Tables

Table 1. Table 1: Model parameters

Model	SI	Ice line	Bouncing	${\alpha }_t$	$u_{\mathrm{frag}}$
NP	No	No	No	${10}^{-2}$	10 $m{s}^{-1}$
SI	Yes	No	No	${10}^{-2}$	10 $m{s}^{-1}$
I	Yes	Yes	No	${10}^{-2}$	10 $m{s}^{-1}$
IU	Yes	Yes	No	${10}^{-2}$	1 $m{s}^{-1}$
IB	Yes	Yes	Yes	${10}^{-2}$	10 $m{s}^{-1}$
IA	Yes	Yes	No	${10}^{-4}$	10 $m{s}^{-1}$
IAB	Yes	Yes	Yes	${10}^{-4}$	10 $m{s}^{-1}$

Table 2. Table 2: Simulation results

		Planetesimal mass in radial bins ($M_{\oplus }$)
Model	Summary	Total	1 au	1-3 au	3-10 au	10-30 au	30-100 au	$>$100 au
SI	Streaming instability only	76.34	0.01	0.27	4.03	4.16	7.98	59.89
I	SI + water ice line	122.62	0.02	0.32	6.78	9.03	12.78	93.70
IU	I + reduced $u_{\mathrm{frag}}$	153.08	0.01	0.17	5.09	7.01	11.92	128.87
IB	I + bouncing barrier	122.98	0.00	0.08	5.26	10.84	13.12	93.67
IA	I + reduced ${\alpha }_t$	74.52	0.08	0.15	0.16	0.05	0.16	73.92
IAB	IAB + bouncing barrier	83.58	1.03	7.35	1.62	0.44	0.18	72.95

Equations

$$t_{\nu }=\frac{1}{{\alpha }_{\mathrm{v}}\Omega }{\left(\frac{H}{R}\right)}^{-2},$$

$$t_{\mathrm{s}}=\frac{{\rho }_{\mathrm{s}}a}{\rho c_{\mathrm{s}}}\sqrt{\frac{\pi }{8}},$$

$$\Delta u_{\mathrm{coll}}\sim c_{\mathrm{s}}\sqrt{{\alpha }_{\mathrm{t}}\mathrm{St}}.$$

$${\mathrm{St}}_{\mathrm{frag}}\sim \frac{u_{\mathrm{frag}}^2}{{\alpha }_{\mathrm{t}}{c}_{\mathrm{s}}^2},$$

$$u_{\mathrm{drift}}=\frac{-2\Delta v}{\mathrm{St}+{\mathrm{St}}^{-1}}$$

$$\Delta u_{\mathrm{drift}}=\frac{c_{\mathrm{s}}^2}{u_{\mathrm{kep}}}\left|\frac{\partial \ln P}{\partial \ln r}\right||{\mathrm{St}}_1-{\mathrm{St}}_2|.$$

$${\mathrm{St}}_{\mathrm{df}}=\frac{u_{\mathrm{frag}}{u}_{\mathrm{kep}}}{c_{\mathrm{s}}^2(1-N)}{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1},$$

$${\mathrm{St}}_{\mathrm{drift}}\sim \frac{{\Sigma }_{\mathrm{p}}}{{\Sigma }_{\mathrm{g}}}\frac{u_{\mathrm{kep}}^2}{c_{\mathrm{s}}^2}{\left|\frac{\partial \ln P}{\partial \ln r}\right|}^{-1}.$$

$$\frac{\partial {\Sigma }_{\mathrm{g}}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}\left(\nu {\Sigma }_{\mathrm{g}}\sqrt{r}\right)\right]-{\dot{\Sigma }}_{\mathrm{pe}}(r,t),$$

$$\nu ={\alpha }_{\mathrm{v}}{c}_{\mathrm{s}}^2/{\Omega }_{\mathrm{K}},$$

$$L_{\mathrm{FUV},\mathrm{acc}}\approx 0.04\left(\frac{0.8G{M}_{\star }\dot{M}}{R_{\star }}\right),$$

$$\frac{\partial {\Sigma }_{\mathrm{p}}^i}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left[r{\Sigma }_{\mathrm{g}}{D}^i\frac{\partial }{\partial r}\frac{{\Sigma }_{\mathrm{p}}^i}{{\Sigma }_{\mathrm{g}}}-r{\Sigma }_{\mathrm{p}}^i{u}^i\right],$$

$$a_1=\min (a_{\mathrm{frag}},a_{\mathrm{df}},a_{\mathrm{drift}}).$$

$$\frac{{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{g}}}=Z\sqrt{\frac{\mathrm{St}+\delta }{\delta }},$$

$${\alpha }_{\mathrm{t}}$$

$$Z$$

$$Z>Z_2\equiv 10^{-1.86+0.3{\left(0.98+{\log }_{10}\mathrm{St}\right)}^2}.$$

$$Z$$

$$\mathrm{St}$$

$$P_{\mathrm{sat}}=2.53\times 10^{13}\exp \left(\frac{-6070\mathrm{K}}{T}\right)\mathrm{dyne}{\mathrm{cm}}^{-2}$$

$$P_{{\mathrm{H}}_2\mathrm{O}}=n_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{k}}_{\mathrm{B}}T=n_{{\mathrm{H}}_2}{X}_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{k}}_{\mathrm{B}}T,$$

$$n_{{\mathrm{H}}_2}>\frac{2.53\times 10^{13}\mathrm{dyne}{\mathrm{cm}}^{-2}}{X_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{k}}_{\mathrm{B}}T}\exp \left(\frac{-6070\mathrm{K}}{T}\right).$$

$${\log }_{10}\left(\frac{m_{\mathrm{b}}}{\mathrm{g}}\right)=-2.08{\log }_{10}\left(\frac{\Delta u_{\mathrm{coll}}}{\mathrm{cm}{\mathrm{s}}^{-1}}\right)-7.33$$

$${\rho }_{\mathrm{R}}=\frac{9{\Omega }^2}{4\pi G}=\frac{9}{4\pi }\frac{M_{\star }}{r^3},$$

$${\rho }_{\mathrm{p},\mathrm{SI}}=ϵ{\rho }_{\mathrm{g}},$$

$${\rho }_{\mathrm{g},\mathrm{crit}}={\rho }_{\mathrm{g},2.5Myr}\frac{Z_{2.5Myr}}{Z_{\mathrm{crit}}},$$

$${\rho }_{\mathrm{p},\mathrm{SI}}=10^4{\rho }_{\mathrm{g},\mathrm{crit}}$$

$$⟨a⟩=\sum _{i=1}^N\frac{a_i{\Sigma }_{\mathrm{p},i}}{{\Sigma }_{\mathrm{p}}},$$

Full text

accretion disks — planetary systems: protoplanetary disks — planets and satellites: formation — planets and satellites: terrestrial planets

\noaffiliation

1 Introduction

Most young stars are surrounded by a circumstellar disk of gas and dust with a mass between 0.01% and 10% that of the central star (Andrews & Williams, 2005). The disk has a characteristic lifetime of 2-5 Myr, with a significant scatter due to variations in stellar mass and local environment (e.g. Hartmann et al., 1998; Haisch et al., 2001; Mamajek, 2009). Thanks to radial velocity surveys and the results of the Kepler mission it is now established that planets between the mass of Earth and Neptune (“super-Earths”) are among the most abundant in the Galaxy (Howard et al., 2010; Borucki et al., 2011; Batalha et al., 2013). From the bulk density of these planets we know that many of them must possess large H/He envelopes (e.g. Adams et al., 2008; Lopez & Fortney, 2014) which must have been accreted onto a solid core before the disk dissipated. In the pebble accretion model, the solid core is produced by the accretion of cm-size “pebbles” onto a seed planetesimal (e.g. Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Bitsch et al., 2015; Levison et al., 2015). After the disk dissipates, the region close to the star is thought to be populated by planetesimals and Mars-sized planetary embryos (Kokubo & Ida, 1996; Johansen et al., 2015). These bodies are the seed material for the formation of terrestrial planets over the next 100 Myr (e.g. Chambers & Wetherill, 1998; Raymond et al., 2006). This has two implications for planetesimals:

• •

Some planetesimals should form early to allow sufficient time for the formation of super-Earths and giant planet cores before the disk dissipates. These cores experience disk migration, so they must have formed at greater distances than where we see them today. For example, a Jupiter-mass planet at 3 au may come from a seed that formed in the $20-30$ au region (Bitsch et al., 2015).

• •

Some planetesimals must form dry — inside the water ice line — in order to reproduce the very low water content of the terrestrial planets in the solar system. The water present in the terrestrial planets likely comes from a late accretion of icy planetesimals scattered inward (Raymond et al., 2004).

In other words, a complete model of planetesimal formation should be able to produce planetesimals beyond the water ice line early, and it must also be able to produce dry planetesimals in the terrestrial zone at some point (possibly late). Furthermore, exoplanet observations suggest a wide range of super-Earth compositions, which may indicate a diversity in their formation pathways and perhaps also times (e.g. Lopez & Fortney, 2014).

The aim of this paper is to address the formation of planetesimals in a protoplanetary disk model that includes a state-of-the-art description of photoevaporation (based on the recent model of Gorti et al., 2015). Photoevaporation may be key to forming planetesimals by the streaming instability, because pebble-sized particles can only concentrate in regions of elevated metallicity compared to the nominal 1% value of the solar photosphere (Johansen et al., 2009; Bai & Stone, 2010; Carrera et al., 2015). In a complementary work, Drazkowska et al. (2016) also modelled the evolution of dust in an evolving disk. The most salient difference between our work and theirs is that we use a more sophisticated model for photoevaporation. Drazkowska et al. (2016) assume a fixed mass-loss rate that is insensitive to changes in accretion, or the radial and size distribution of dust. Our model is described in sections 2.4 and 3.1. In a real protoplanetary disk, the FUV flux is dominated by photons released in the accretion process, and the resulting FUV heating depends in a complicated way on grain size and distribution.

Another important difference between our work and Drazkowska et al. (2016) is that we use the turbulent viscosity motivated by observation. In a turbulent disk with viscosity $\nu=\,\alpha_{\rm v}\,c_{\rm s}H$ (Shakura & Sunyaev, 1976) where $\,c_{\rm s}$ is the local sound speed and $\,\alpha_{\rm v}$ is the dimension-less turbulence parameter, the timescale for viscous evolution is

[TABLE]

where $\Omega$ is the Keplerian frequency and $H/R\sim 0.1$ is the disk aspect ratio (e.g. Alexander et al., 2014). For disk lifetimes of a few Myr, we obtain $\,\alpha_{\rm v}\sim 10^{-2}$, and that is the value we use (see also Hartmann et al., 1998; Alexander et al., 2014). In contrast, Drazkowska et al. (2016) set $\,\alpha_{\rm v}$ to $10^{-3}$ to facilitate planetesimal formation. We also allow smaller particles to participate in the streaming instability, following the results of Carrera et al. (2015); however this is a minor difference. One interesting feature of Drazkowska et al. (2016) is that the radial drift of solids decreases with increasing dust-to-gas ratio, as they include the back-reaction of the dust on the gas in their models. This facilitates the pile-up of solids in the inner disk, and may play a role in planetesimal formation. In an upcoming paper we will combine our more sophisticated models for photoevaporation and streaming instability with the more sophisticated radial drift model of Drazkowska et al. (2016).

This paper is organized as follows. In section 2 we introduce main concepts pertaining to our disk evolution model. In section 3 we explain how the model is implemented and our key assumptions. Our results are presented in section 4, where we show that planetesimal formation is a natural by-product of the evolution of a protoplanetary disk. We show that planetesimal formation in the outer disk is a natural outcome of disk evolution. We also show that standard prescriptions for disk processes can lead to several Earth masses of planetesimals forming in the terrestrial zone. Then in section 5 we discuss various processes that we could not model adequately in our simulations, and speculate on how they might affect our results. Finally, we summarize our results and draw conclusions in section 6.

2 Main Theoretical Concepts

In the interstellar medium most of the solids are in the form of sub-micron dust (e.g. Weingartner & Draine, 2001; Lefèvre et al., 2014). Inside the protoplanetary disk, these grains stick and grow to larger sizes up to at least millimetre sizes in the outer disk (van der Marel et al., 2013). Growth beyond a few millimetres is inhibited by collisional fragmentation (Brauer et al., 2008). The mechanism by which these sub-cm grains are converted into 100-km planetesimals is not well understood.

2.1 The streaming instability

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) time scales. 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, $\,{\rm St}=\,t_{\rm s}\,\Omega_{\rm K}$, where $\,t_{\rm s}$ is the stopping time and $\,\Omega_{\rm 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,

[TABLE]

where $\,\rho_{\rm s}$ is the material density, $\rho$ is the local gas density, $a$ is the grain size, and $\,c_{\rm 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 $\,{\rm 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.

2.2 Dust growth

The size of the grains is determined by the balance between dust coagulation, collisional fragmentation, and radial drift (e.g. Birnstiel et al., 2012). Silicate grains also experience a bouncing barrier (Zsom et al., 2010). In a standard $\alpha$-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

[TABLE]

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

[TABLE]

where $\,u_{\rm 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 $\,{\rm St}_{\rm frag}$ (Equation (4)). As a canonical example, $\,u_{\rm frag}=10\,{\rm m\,s}^{-1}$, $\,c_{\rm s}=800\,{\rm m\,s}^{-1}$, and $\,\alpha_{\rm v}=\,\alpha_{\rm t}=10^{-2}$ would give $\,{\rm St}_{\rm frag}\sim 0.01$, independent of the gas density. The $\,u_{\rm frag}=10\,{\rm m\,s}^{-1}$ 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 $\,\Sigma_{\rm p}/\,\Sigma_{\rm g}\sim 0.03$ (Carrera et al., 2015).

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

2.3 Radial drift and the water ice line

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

[TABLE]

where $\Delta v=u_{\rm kep}-u_{\rm 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 $\Delta u_{\rm coll}$ is dominated by differential drift. The relative drift speed for small particles ($\,{\rm St}\ll 1$) is given by the expression

[TABLE]

Following the notation of Birnstiel et al. (2012), we write $\,{\rm St}_{2}=N\,{\rm St}_{1}<\,{\rm St}_{1}$. Setting $\Delta u_{\rm drift}=u_{\rm frag}$ one gets the drift-driven fragmentation limit,

[TABLE]

where $P$ is the gas pressure and $r$ is the radial coordinate (Birnstiel et al., 2012). This quantity depends weakly on $r$, and $\,{\rm St}_{\rm 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

[TABLE]

In the outer disk, where $u_{\rm kep}$ is small, $\,{\rm St}_{\rm drift}<\,{\rm St}_{\rm 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.

2.4 Photoevaporation

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\sim 10^{4}\,{\rm K}$ ($\,c_{\rm s}\sim 10\,{\rm km\,s}^{-1}$), which can launch a flow at $\sim 1$ au for a star of mass 1 $M_{\odot}$. 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 $\sim 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, b).

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}\sim 10^{30}\,{\rm erg}\;{\rm s}^{-1}$. The photon energy ranges from $\sim 0.1-10$ keV and gas temperatures from a few 100 K to over $10^{4}$ K. The softer X-ray photons heat the gas to higher temperatures ($T\sim 1-2\times 10^{4}\,{\rm 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 (PAH) 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 ($\sim 3-10$ au) and in the outer disk ($\gtrsim 50$au). FUV photoevaporation in the outer disk is particularly significant for disk evolution as 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).

2.5 Dust evolution at high dust-to-gas ratio

The dust evolution model that we presented in this section was derived by Birnstiel et al. (2012) under an assumption that the gas is 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.

2.6 Disk winds

Although we do not include disk winds in our base model, we will discuss them briefly here. Disk winds are 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 $\,\Sigma_{\rm p}/\,\Sigma_{\rm g}$ may promote planetesimal formation by the streaming instability.

3 Methods

All our simulations begin at the end of the Class I phase, after about 0.5 to 1 Myr of evolution. At this stage the central star has acquired most of its final mass and the envelope has dissipated to reveal a gravitationally stable, accreting, protoplanetary disk (Bell et al., 2013; Dunham et al., 2014). The model parameters are all as described in Gorti et al. (2015), as in their fiducial run. We briefly summarize some of these here for completion. The initial disk mass is $0.1M_{\odot}$, distributed in a $1/r$ power law extending to 200au; this profile is gravitationally stable everywhere. The disk is allowed to viscously relax before the simulations are begun, and attains a self-similar (viscous) radial profile (in $\sim 10^{5}$ years). The gas-to-dust ratio is initially 100, and the dust sizes (32 bins, in the range 50Å-1cm) are computed as the disk evolves. The 1-D disk evolution model is solved using a Crank-Nicholson method with an additional criterion imposed by the rate of change of surface density due to photoevaporation ($\Delta t<0.1\Sigma/\dot{\Sigma}$). The temperature and density structure of the disk is determined from the instantaneous surface density profile ($\Sigma(r,t)$) using a 1+1D thermochemical disk model (Gorti et al., 2015).

Figure 1 shows the initial surface density and temperature profile shared by all our runs. The surface density profile is consistent with observed disks (e.g. Williams & Cieza, 2011, figure 2). Note that the disk is has a large radial extent (beyond 1,000 au) but the surface density drops exponentially after 100 au. Because of the low surface density, the disk is gravitationally stable. The Toomre Q parameter reaches a minimum of $Q=6$ at around 72 au, so that the stability criterion ($Q>1$) is met everywhere. The initial dust distribution follows the gas density profile with a dust-to-gas ratio of $Z=0.01$ and a local size distribution following the prescription of Birnstiel et al. (2011).

We performed seven simulations in which we explore how different physical processes impact the formation of planetesimals. Our most complete model is Model IAB. The models are summarized in Table 1. Here we describe the seven models in more detail:

• •

Model NP (“no planetesimals”) is a disk evolution model without planetesimal formation. The model is identical to that of Gorti et al. (2015) and includes viscous accretion, photoevaporation, and dust evolution.

• •

Model SI (“streaming instability”) adds planetesimal formation by the streaming instability to Model NP using the results of Carrera et al. (2015).

• •

Model I (“ice line”) adds the water ice line including condensation, sublimation, and diffusion of water vapor.

• •

Model IU (“I + $\,u_{\rm frag}$”) is a duplicate of Model I with the fragmentation speed reduced to $\,u_{\rm frag}=1\,{\rm m\,s}^{-1}$. This is a test to determine whether the simulation results are sensitive to $\,u_{\rm frag}$.

• •

Model IB (“I + bouncing”) includes the water ice line and the bouncing barrier for silicate grains inside the ice line (Güttler et al., 2010).

• •

Model IA (“I + $\,\alpha_{\rm t}$”) includes the water ice line and a reduced midplane turbulence of $\,\alpha_{\rm t}=10^{-4}$ which is more consistent with a dead zone (Oishi et al., 2007), and is consistent with ALMA observations of outer disks (Flaherty et al., 2015).

• •

Model IAB (“I + $\,\alpha_{\rm t}$ + bouncing”) is our most complete model. It includes the water ice line, the low midplane turbulence, and bouncing barrier.

3.1 Photo-evaporating disk with dust evolution

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 one-dimensional (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 photo-evaporation,

[TABLE]

where $r$ is the radial coordinate, $\nu$ is the viscosity, $\,\Sigma_{\rm g}$ is the gas surface density, and $\dot{\Sigma}_{\rm pe}$ is the instantaneous photo-evaporation rate. The viscosity term $\nu$ follows the usual $\alpha$-model of Shakura & Sunyaev (1976),

[TABLE]

where $\,c_{\rm s}$ is the local sound speed and $\,\Omega_{\rm K}$ is the Keplerian frequency. As discussed earlier, the viscous $\,\alpha_{\rm v}$ may differ from the midplane turbulence $\,\alpha_{\rm 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 black-body spectrum with a temperature of 9,000K (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

[TABLE]

where $M_{\star}$ and $R_{\star}$ 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}\sim 10^{30}\,{\rm erg}\;{\rm s}^{-1}$ for $M_{\star}=1M_{\odot}$ (Flaccomio et al., 2003).

To compute $\dot{\Sigma}_{\rm 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{\Sigma}_{\rm pe}$ at each point $(r,z)$. We consider the photoevaporative flow to be launched from the height $z$ where $\dot{\Sigma}_{\rm 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 $\,\Sigma_{\rm p}^{i}$ with particle size $a^{i}$ 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

[TABLE]

where $u^{i}$ is the total radial velocity of the dust due to gas drag and the radial pressure gradient, $D^{i}=\nu/(1+\,{\rm St}_{i}^{2})$ is the dust diffusivity, and $\,{\rm 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_{\rm drift}$ and $a_{\rm df}$. We then compute

[TABLE]

Once the maximum particle size $a_{1}$ is computed, the dust mass is distributed across the size bins following the prescription of Birnstiel et al. (2011). All our models follow the same basic prescription. Most models use fragmentation speed $\,u_{\rm frag}=10\,{\rm m\,s}^{-1}$ and $\,\alpha_{\rm t}=10^{-2}$, but Model IU uses $\,u_{\rm frag}=1\,{\rm m\,s}^{-1}$, and Models IA and IAB use $\,\alpha_{\rm t}=10^{-4}$. In addition, Models IB and IAB add a maximum size limit due to the bouncing barrier to Equation (13).

3.2 Streaming instability

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}>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 $\,\rho_{\rm p}/\,\rho_{\rm g}$ is given by the particle size and diffusion coefficient $\delta$,

[TABLE]

where $Z=\,\Sigma_{\rm p}/\,\Sigma_{\rm g}$. In the case of turbulence driven by the magneto-rotational instability, $\delta\approx\,\alpha_{\rm 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 $\,\rho_{\rm p}/\,\rho_{\rm g}>1$ can be thought of as a minimum $Z$ or a maximum value for the midplane turbulence $\,\alpha_{\rm t}$,

[TABLE]

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 $\,\alpha_{\rm t}$ is negligible. We computed a least-squares fit for the results of Carrera et al. (2015) and obtained the condition,

[TABLE]

In line with the results of Carrera et al. (2015), we also require $\,{\rm St}\geq 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_{1}$ and $Z_{2}$ be right-hand-side of Equations (16) and (17) respectively. Together, the conditions for the streaming instability become

[TABLE]

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\geq\,Z_{\rm crit}$ , we remove 90% of the mass in that bin. If $Z<\,Z_{\rm crit}$, we compute $Z$ for the top two bins, and re-compute $\,Z_{\rm 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 $\,{\rm 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.

3.3 Water ice line

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

[TABLE]

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

[TABLE]

where $n_{\rm H_{2}}$ and $n_{\rm H_{2}O}$ are the number densities of hydrogen and water, $\,{\rm k_{B}}$ is the Boltzmann constant and $X_{\rm H_{2}O}$ is the water abundance. Ice forms whenever $P_{\rm H_{2}O}>P_{\rm sat}$, which simplifies to

[TABLE]

The run begins with $X_{\rm H_{2}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 d$m/m$ across all grain size bins. This differs from the approach of Ros & Johansen (2013), who kept d$a$ 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 $\Sigma_{\rm ice}$ and $\Sigma_{\rm vapor}$ in order to determine the water fraction in dust or gas. Inside the ice line, $\Sigma_{\rm ice}$ is locally released into $\Sigma_{\rm vapor}$ which advects and diffuses like the rest of the disk gas.

3.4 Bouncing barrier

Inside the water ice line, the solid component is made entirely of silicate grains which may experience a bouncing barrier (Zsom et al., 2010; Güttler 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, as 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_{\rm b}$, is

[TABLE]

Inside the ice line, we use Equation (23) to compute the corresponding grain size at the bouncing barrier $a_{\rm b}$ and replace Equation (13) with the minimum of $a_{1}$ and $a_{\rm b}$.

3.5 Reduced midplane turbulence

Most of our models follow a simple $\alpha$-disk prescription (Shakura & Sunyaev, 1976) in which $\,\alpha_{\rm t}=\,\alpha_{\rm v}=10^{-2}$ everywhere in the disk. But 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 $\delta=\,\alpha_{\rm t}$ in Equation (14)).

First, from around 0.2-1 au 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 magneto-rotational instability by Oishi et al. (2007) give an $\,\alpha_{\rm t}$ of $10^{-5}$ inside the dead zone. Alternatively, simulations of the vertical shear instability give $\,\alpha_{\rm t}=10^{-4}$ on the midplane (Stoll & Kley, 2014). While there are no strong theoretical constraints beyond 100 au, ALMA CO observations indicate a three-sigma upper limit of $\alpha<10^{-3}$ for the outer regions ($R>30$ au) of the disk around around HD 163296 (Flaherty et al., 2015). Recent ALMA observations of the TW Hya disk (Teague et al., 2016) suggest a value of $\alpha$ closer to our adopted default value of 0.01.

Our choice thus far of $\,\alpha_{\rm t}=10^{-2}$ 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 $\,\alpha_{\rm t}$: $\,\alpha_{\rm v}$ is used to compute the viscous evolution of the disk, and $\,\alpha_{\rm 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 $\,\alpha_{\rm t}=10^{-4}$. This lowers collision speeds by a factor 10 while maintaining the nominal accretion rate of the gas and giving reasonable disk lifetimes.

4 Results and discussion

4.1 Planetesimal formation

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 $70M_{\oplus}$ of dust. This much dust is inconsistent with observations of debris belts (Wyatt, 2008). In Model SI, only $0.002M_{\oplus}$ of dust remain, as 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.

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

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.

Figures 6 and 7 give a general picture of how planetesimal formation occurs in Model SI. Planetesimals always begin to form early (after $\sim$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_{\rm 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 our models. As the disk evolves, its outer radius decreases due to FUV photoevaporation and planetesimal formation moves inward.

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_{\rm 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).

4.2 Water ice line

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_{\oplus}$ of solids, compared to only 313 $M_{\oplus}$ for Model SI. Water ice maintains a high fragmentation barrier of $\,u_{\rm frag}\sim 10\,{\rm m\,s}^{-1}$ (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)). But 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_{\oplus}$ vs 76 $M_{\oplus}$) 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.

4.3 Fragmentation barrier

As noted earlier, the $\,u_{\rm frag}=10\,{\rm m\,s}^{-1}$ 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_{\rm frag}=1\,{\rm m\,s}^{-1}$. This change increases the total planetesimal mass (153 $M_{\oplus}$ vs 123 $M_{\oplus}$). 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.

4.4 Bouncing barrier

Model IB is a copy of Model I (both have a water ice line and $\,u_{\rm frag}=10\,{\rm m\,s}^{-1}$) 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,

• •

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 $\sim$10 au because the small bouncing-limited grains can diffuse outward.

• •

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).

• •

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.

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^{-3}M_{\oplus}$). When the disk is about to disappear ($\sim$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.

4.5 Midplane turbulence

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 $\,\alpha_{\rm t}=10^{-4}$. So in Model IA (Table 1) we use two parameters: $\,\alpha_{\rm v}=10^{-2}$ to compute the viscous evolution of the disk, and $\,\alpha_{\rm t}=10^{-4}$ to compute the collision speeds between solid particles. The fragmentation-limited Stokes number scales as $1/\,\alpha_{\rm 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 $\sim$ 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{\Sigma}_{\rm 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.

4.6 Complete model

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.4M_{\oplus}$ inside 3 au, Model IAB produces $8.4M_{\oplus}$ – 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.

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_{\oplus}$ 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}_{\rm acc}\propto\,\Sigma_{\rm g}$). This has two key consequences:

• •

When the accretion rate decreases, the accretion-generated FUV flux also drops dramatically. This effectively shuts down FUV photo-evaporation 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.

• •

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.

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 photo-evaporation 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.

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.

5 Other factors

Planetesimal formation is a complex, multi-physics problem. We have shown that photoevaporation, dust evolution, the water ice line, and the streaming instability all probably play an important role in the formation of planetesimals. In this section we discuss various physical processes that we could not model adequately in our simulations, and we speculate on how these processes might change our results.

5.1 Disk viscosity

We note that all the results presented here are sensitive to assumptions made on viscous disk evolution. We have chosen a simple $\alpha$-parametrization of disk viscosity, and moreover assume that $\alpha$ 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, although considerable progress is being made in this direction (e.g. Bai et al., 2016; Simon et al., 2015). Disk lifetimes due to a combination of viscous evolution and photoevaporation also depend on the choice of the numerical value of $\alpha$ (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 modelled 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.

5.2 Laminar accretion

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 artefact of our one-dimensional 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 $\,\alpha_{\rm t}=10^{-4}$ 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)). But 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.

5.3 Roche density

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 $\,\rho_{\rm p}$ must reach the Roche density,

[TABLE]

where $\Omega$ is the Keplerian frequency, $G$ is the gravitational constant, $M_{\star}$ is the mass of the central star, and $r$ is the orbital distance. When $\,\rho_{\rm p}>\,\rho_{\rm 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 $\,\rho_{\rm p}>\,\rho_{\rm 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,

[TABLE]

where $\,\rho_{\rm g}$ is the gas density at the midplane, and $\rho_{\rm p,SI}$ is the particle density produced by the streaming instability inside a filament. The value of $\epsilon$ 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 $\,\rho_{\rm p}/\,\rho_{\rm g}>1$ on the midplane. The streaming instability then proceeds to concentrate particles into denser filaments, so that $\epsilon=\rho_{\rm p,SI}/\,\rho_{\rm g}\gg 1$ inside the filaments. However, the value of $\epsilon$ has only been measured precisely for $\,{\rm St}=0.3$ particles, which have $\epsilon$ going locally up to at least $10^{4}$ (Johansen et al., 2015).

As mentioned in the discussion, the streaming instability occurs whenever $\,\Sigma_{\rm p}/\,\Sigma_{\rm g}>\,Z_{\rm crit}$, where $\,Z_{\rm crit}$ is a function of the particle Stokes number. From the point of view of the streaming instability, it makes no difference if $\,Z_{\rm crit}$ is achieved by removing gas or by adding solids. However, Equation (25) shows that it is easier to reach $\rho_{\rm p,SI}>\,\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 $\,\rho_{\rm g}$, $Z$, St, and $\,Z_{\rm crit}(\,{\rm 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 $\,\rho_{\rm 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_{\rm crit}$ for the streaming instability to occur,

[TABLE]

where $\rho_{\rm g,crit}$ is the (critical) midplane density when the streaming instability begins. We then compute $\rho_{\rm p,SI}$ assuming $\epsilon=10^{4}$,

[TABLE]

Except for the value of $\epsilon$, 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_{\rm crit}$ and larger $\rho_{\rm p,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,

[TABLE]

where $a_{\rm i}$ and $\Sigma_{{\rm p},i}$ are the size and surface density of particles in the $i^{\rm th}$ particle bin, and $\,\Sigma_{\rm 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 $\,{\rm 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 $\,{\rm St}\geq 0.003$ before photoevaporation.

The bottom-left plot of Figure 12 shows the critical dust-to-gas ratio $\,Z_{\rm crit}$ assuming a uniform particle size $\langle a\rangle$. Finally, we show $\rho_{\rm p,SI}/\,\rho_{\rm R}$ in the bottom-right plot, following Equations (24) and (27). We also mark the critical value $\rho_{\rm p,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 $\,{\rm St}$ increases. As $\,\rho_{\rm g}\rightarrow 0$, $\,{\rm St}$ approaches the fragmentation barrier. To conclude this section, we list three more processes that may significantly increase $\,\rho_{\rm p}$ inside a filament and have not been adequately explored in the literature.

• •

Inside the filament, orbits are strongly Keplerian and collision speeds between particles may be much smaller than $\,c_{\rm s}\sqrt{\alpha\,{\rm 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.

• •

Even when a stable filament fails to reach the Roche density, $\,\rho_{\rm 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.

• •

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 $\,\rho_{\rm p}$ should begin to increase. As $\,\rho_{\rm g}\rightarrow 0$, $\,\rho_{\rm p}\rightarrow\infty$ on the midplane and the filament may experience a 1-dimensional 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,SI}<\,\rho_{\rm R}$ initially, $\,\rho_{\rm 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.

6 Summary and conclusions

In this paper we have presented a global disk evolution models that incorporate photoevaporation and the formation of planetesimals by the streaming instability. Our key result is that planetesimals are a natural by-product of protoplanetary disk evolution. Planetesimal formation does not require a fundamental change in how we understand disk evolution — it merely requires a multi-physics view of planet formation that combines the work done in many areas of protoplanetary disk models and planetesimal formation models over recent years. Concretely, we offer two key findings:

• •

Planetesimals form early in the outer disk as a result of the streaming instability and FUV photo-evaporation. This outcome seems unavoidable, as even the most conservative models produce at least $\sim 60M_{\oplus}$ of planetesimals beyond 100 au. This result is very robust and does not require fine-tuning of any disk parameter. If these planetesimals further assimilate to form planets in the outer disk, their dynamical effects on the protoplanetary disk may be observable with high-sensitivity facilities such as ALMA.

• •

The formation of planetesimals in the terrestrial zone depends mostly on the formation of an inner cavity in the disk. Before the cavity forms, the streaming instability is ineffective and cannot easily produce planetesimals. Once a cavity forms, the solid material can no longer accrete onto the star. As photoevaporation removes the remaining gas, the dust-to-gas ratio continues to increase and the remaining solids are forced to eventually form planetesimals. Our most complete model, Model IAB (which includes the water ice line, bouncing barrier, and low midplane turbulence), produces 8.4 $M_{\oplus}$ of planetesimals in the terrestrial zone. All other models produce at most 0.3 $M_{\oplus}$. What sets Model IAB apart is that it forms a cavity sooner than the others.

All our models produce very few planetesimals in the giant planet zone (3-30 au), and those that do form appear too late and could probably not grow into giant planet cores before gas dissipation (Bitsch et al., 2015). It appears that the seed planetesimals that form giant planets must form through some other type of process, perhaps related to early particle pile-ups near ice lines (Sirono, 2011b, a).

Finally, our results suggest new questions and new opportunities for future research:

• •

How do disk winds affect planetesimal formation by the streaming instability?

• •

How does a particle filament produced by the streaming instability evolve if it does not initially reach the Roche density?

• •

What is the collision speed between particles inside those filaments?

Acknowledgements

D.C. acknowledges Chao-Chin Yang for many helpful discussions and for his insights on the Roche density and the streaming instability. The authors acknowledge the support from the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017, 2014.0048), the Swedish Research Council (grants 2011-3991 and 2014-5775), and the European Research Council Starting Grant 278675-PEBBLE2PLANET. U.G. is supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Co-operation Agreement Notice NNH13ZDA017C issued through the Science Mission Directorate, and a grant from the National Science Foundation (AST1313003). This project made use of NASA HEC supercomputing resources.