Free Access
Issue
A&A
Volume 572, December 2014
Article Number A78
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201424809
Published online 01 December 2014

© ESO, 2014

1. Introduction

The number of known exoplanets is steadily rising, and it already exceeds 18001. Indirect methods indicate that every star in our Galaxy has at least one planet on average (Cassan et al. 2012). Planet formation is ubiquitous and results in a variety of planetary system architectures. However, the details remain a mystery, as they are impossible to follow with direct observations, which only cover the very first stages of planet formation, when protoplanetary disks consist mainly of gas and small dust particles, and the result in the form of the exoplanet population.

Planets form in disks surrounding young stars. The solid material initially consists only of μm-sized grains, which are already present in the interstellar medium (Pagani et al. 2010). Analytical and numerical models of dust evolution aim to explain the growth of the primordial μm-sized grains through around 40 orders of magnitude in mass to planets that are greater than 1000 km in size. However, the growth already encounters serious obstacles at the beginning of the size range, and these are referred to as growth barriers. There are, in principle, two kinds of growth barriers: those resulting from collisional physics of dust aggregates, and those resulting from radial drift timescale. The first type are known as the bouncing and fragmentation barriers. They arise when the impact speeds are too high to allow aggregate sticking. An overview of collisional physics of dust aggregates is given by Güttler et al. (2010). Implementing a complex model derived from laboratory work into a dust coagulation code, Zsom et al. (2010) found that silicate grain growth is inhibited by aggregate compaction and bouncing already at millimeter sizes. Icy particles, which can exist outside the snow line, are considered to be more sticky and presumably avoid the bouncing behavior (Wada et al. 2011). However, even without bouncing, impact velocities reaching a few tens m s-1 are too high to allow grain growth and instead lead to fragmentation (Blum & Münch 1993; Blum & Wurm 2008; Meru et al. 2013b). At even smaller sizes, for <mm-sized aggregates, there was another sticking barrier found, called the charge barrier, which results from the electrostatic charge of small aggregates (Okuzumi 2009).

Barriers of the second kind are connected to the radial drift, which is caused by the sub-Keplerian rotation of pressure supported gas disk. The dust grains interact with the gas via aerodynamic drag, lose their angular momentum, and drift toward the central star. For roughly decimeter-sized grains, the drift timescale is typically shorter than the growth timescale. Thus, such grains are removed from a given location in the disk before they can produce any larger aggregates (Weidenschilling 1977; Nakagawa et al. 1986; Brauer et al. 2007, 2008a).

There are several concepts that facilitate overcoming the growth barriers and producing bodies which are held together by self-gravity, for example planetesimals with minimum size of 100 m (Benz & Asphaug 1999). Such concepts are, among others, pressure bumps (Whipple 1972; Kretke & Lin 2007; Brauer et al. 2008b; Pinilla et al. 2012; Drążkowska et al. 2013), dust accumulation in vortices (Barge & Sommeria 1995; Bracco et al. 1999; Klahr & Bodenheimer 2006; Lyra et al. 2009; Meheut et al. 2010), sweep-up growth scenario (Windmark et al. 2012a,b; Garaud et al. 2013; Meru et al. 2013a; Drążkowska et al. 2014), ice condensation (Cuzzi & Zahnle 2004; Ros & Johansen 2013), and ultra-porous grain growth (Okuzumi et al. 2012; Kataoka et al. 2013a). A review of the current understanding of planetesimal formation can be found in Johansen et al. (2014). In this paper, we focus on planetesimal formation by gravitational collapse of dense dust clumps produced by two-fluid instability known as the streaming instability.

One of the early scenarios of planetesimal formation was proposed by Goldreich & Ward (1973). This model relied on dust settling toward the midplane of the protoplanetary disk and the formation of a thin, unstable disk of solids. This thin disk would then fragment because of the gravitational instability and the fragments would collapse directly into planetesimals. However, it was found that formation of such a thin midplane layer is not possible, because as soon as the local dust-to-gas ratio in the midplane exceeds unity, the shear instabilities occur (Weidenschilling 1980, 1995; Cuzzi et al. 1993), in particular the Kelvin-Helmholtz instability (Johansen et al. 2006; Barranco 2009).

Formation of planetesimals by gravitational instability is still possible if a strong clumping of dust is present. The complicated interactions between gas and dust may lead to the development of a powerful instability, called streaming instability, that causes strong inhomogeneities in the dust density (Goodman & Pindor 2000; Youdin & Goodman 2005). The instability is most efficient for high dust-to-gas ratios and large particles, with stopping times corresponding to the orbital timescale. Youdin & Johansen (2007) and Johansen et al. (2007) confirmed that this instability is able to produce very dense dust clumps and leads to rapid planetesimal formation. This result was later confirmed by high resolution numerical simulations (Johansen et al. 2011), as well as by other authors using different methods and codes (Balsara et al. 2009; Tilley et al. 2010; Bai & Stone 2010a; Jacquet et al. 2011; Kowalik et al. 2013).

The numerical models of the streaming instability are generally computationally expensive and do not allow us to perform wide parameter studies or use extended simulation domains. They are also typically initialized with an arbitrary particle size distribution, often with even just a single particle species. The aggregates sizes used in the simulations are usually large, corresponding to Stokes numbers of 10-2−1, whereas in a primordial protoplanetary nebula we expect μm-sized grains, which correspond to Stokes numbers of 10-6. Larger grains may be produced by coagulation, but their maximum size is limited by the growth barriers at Stokes numbers in the range 10-4−10-1. Choosing only a narrow size distribution and assuming that all of the dust particles are large leads to a kind of super-critical initial condition and the instability is triggered very quickly. We expect that gradual dust growth leads to a quite wide and highly problem-dependent size distribution, meaning that only a fraction of the largest grains will be able to participate in clumping (Bai & Stone 2010a). These particles can already form planetesimals, while other particles are still growing and gradually refill the population of big grains.

In addition to the abundant large grains, streaming instability also requires that the vertically integrated dust-to-gas ratio is super-solar and the local dust-to-gas ratio is higher than unity, which means that it can only happen in a dense midplane layer. Such a midplane layer can be relatively easily formed in a dead zone, where no turbulent diffusion is present. In order to investigate the interplay of dust coagulation and planetesimal formation in the streaming instability triggered in a dead zone of protoplanetary disks, we build a semi-analytical model of the latter effect into our dust coagulation code based on the Monte Carlo algorithm and the representative particle approach (Zsom & Dullemond 2008; Drążkowska et al. 2013). We base our model of the streaming instability on the work of Bai & Stone (2010a,c). We describe the numerical model in Sect. 2. Basing our work on this approach, we present simple estimates in Sect. 3 and results of the full numerical models in Sect. 4. We offer an analytical formula that explains these results in Sect. 4.3. We discuss the limitations of our approach in Sect. 5 and summarize our work in Sect. 6.

2. Numerical methods

2.1. Dust evolution

We use the Monte Carlo method and the representative particle approach first described by Zsom & Dullemond (2008) to model dust evolution. We developed this method by adding an adaptive grid routine, which enables multi-dimensional models. The code was presented by Drążkowska et al. (2013). For this paper, we focus on local models, only including the vertical dimension. We have already tested our code with such setups (see Sect. 4 of Drążkowska et al. 2013) and we found a good agreement with the results of Dullemond & Dominik (2005) and Zsom et al. (2011). We have also noticed that the adaptive grid routine allows us to obtain very high resolution of the midplane layer and thus to follow dust evolution in dead zones well.

The representative particle approach relies on the assumption that the evolution of many dust particles can be resolved by following the evolution of only a limited number of so called representative particles. One representative particle represents a swarm of identical physical particles. In our code, all the swarms have equal mass, thus the mass of one swarm Mswarm is an adequate fraction of the total dust mass MtotMswarm=Mtot/Nswarms,\begin{equation} \label{mswarm} M_{\rm{swarm}} = M_{\rm{tot}}/N_{\rm{swarms}}, \end{equation}(1)where Nswarms is the number of representative particles, or swarms, present in our simulation (for discussion of the approach restrictions, see Drążkowska et al. 2014). The collisions between dust particles are modeled with a Monte Carlo algorithm, described in detail in Zsom & Dullemond (2008) and Drążkowska et al. (2013).

In the original code, the particles were influenced by the gas but the gas did not change its state. Here, we assume that the particles are active, i.e., we enable the back-reaction on gas. The effects of the back-reaction are implemented in a semi-analytical way. First, a very crude approach of implementing the turbulence triggered by the streaming instability was already made in the Drążkowska et al. (2013). However, that approach did not include the dust clumping and the possibility of planetesimal formation, which we include now. As we do not directly model the hydrodynamics of the gas, we cannot capture the two-fluid interactions explicitly. We model dust growth and settling, and include the effects of back-reaction based on results of the direct numerical simulations that investigated the dynamics of dust grains in the midplane of protoplanetary disks presented by Bai & Stone (2010a,c), and the analytical model by Takeuchi et al. (2012). We describe this approach in Sect. 2.2.

The dust evolution in a protoplanetary disk is driven by its interactions with gas. We assume that the gas disk is described by the minimum mass solar nebula (MMSN) model proposed by Hayashi (1981), where the surface density follows Σg=1700×(rAU)-1.5gcm-2,\begin{equation} \Sigma_{\rm g}=1700\times\left(\frac{r}{\rm AU}\right)^{-1.5} {\rm g~cm}^{-2}, \end{equation}(2)where r is radial distance to the central star of mass 1 M. The surface density of dust is parametrized by metallicity2Z such that Σd = Z × Σg. The vertical structure of the gas is described by the local hydrostatic equilibrium and thus the gas density follows ρg(z)=Σg2πHgexp(z22Hg2),\begin{equation} \rho_{\rm g}(z)=\frac{\Sigma_{\rm g}}{\sqrt{2\pi}H_{\rm g}}\exp\left(\frac{-z^2}{2H_{\rm g}^2}\right), \end{equation}(3)where z is the distance to the midplane and Hg = cs/ ΩK is the pressure scale height of gas defined by the sound speed cs and orbital frequency ΩK. We assume an isothermal disk with temperature T=280×(rAU)-0.5K.\begin{equation} T = 280\times\left(\frac{r}{\rm AU}\right)^{-0.5} {\rm K}. \end{equation}(4)To investigate dust particles dynamics, it is convenient to use the Stokes number defined as St=tsΩK,\begin{equation} {\it St}=t_{\rm{s}} \Omega_{\rm{K}}, \end{equation}(5)where ts is the stopping time of the dust particle. The stopping time of a particle determines the timescale that the particle needs to adjust its velocity to the velocity of the surrounding gas. The exact expression that we use to compute the ts depends on the particle radius, and we use the formulas given by Weidenschilling (1977). The Stokes number can be used as a particle-gas coupling strength indicator. The particles with St ≪ 1 are completely coupled and follow the motion of the gas, whereas the particles with St ≫ 1 are fully decoupled from the gas. The particles that have St ≈ 1, sometimes called pebbles, are marginally coupled to the gas and suffer the most from this interaction by acquiring high drift and impact speeds. However, these are also the grains that can trigger the streaming instability and form planetesimals, which is what we are investigating in this paper.

In a laminar disk, dust settles because of the vertical component of the star gravity and gas drag. We use the vertical velocity that implements the damped oscillations of large grains, such that the sedimentation rate is consistent with the results of Carballido et al. (2011): vz=zΩKSt1+St2·\begin{equation} \label{vsett} v_{{z}} = -z \Omega_{\rm{K}} \frac{\it St}{1+{\it St}^2}\cdot \end{equation}(6)If turbulence is present, we use the α prescription (Shakura & Sunyaev 1973), and the turbulent mixing of dust is implemented as random kicks (Ciesla 2010). The height of the dust layer is regulated by the turbulence strength and settling. Including orbital oscillations of grains with St> 1, we get the scale height of the dust layer (Carballido et al. 2006, 2011; Youdin & Lithwick 2007) Hd=Hgαα+St·\begin{equation} \label{hdust} H_{\rm{d}} = H_{\rm{g}} \sqrt{\frac{\alpha}{\alpha+{\it St}}}\cdot \end{equation}(7)The collisions between dust aggregates are driven by the Brownian motion, relative radial, transverse and vertical drift, and turbulence. Although we do not include the radial drift of particles explicitly, we keep the drift speeds as a source of impact velocities, to correctly capture the growth rates and collisional physics.

2.2. Streaming instability

In a laminar disk, all dust grains would settle to the midplane and form a very thin and gravitationally unstable midplane layer (Nakagawa et al. 1981). However, the shear between such a layer and gas triggers turbulence that maintains its height at a level that typically does not allow direct gravitational collapse. The strength of turbulence triggered by the midplane instability was recently investigated by Takeuchi et al. (2012). They estimated the turbulence efficiency from the energy supplied by the radial drift of dust. Irrespectively of the turbulence mechanism, which can be both the Kelvin-Helmholz and streaming instability, they find that for particles with St< 1 the resulting turbulence strength can be parametrized with α=[(C1CeffηZ)2/3+(C2CeffηZ-1)-2]-1St,\begin{equation} \label{alpha} \alpha = \left[\left({C}_1{C}_{\rm eff}\eta Z\right)^{-2/3}+\left({C}_2 {C}_{\rm eff}\eta Z^{-1}\right)^{-2}\right]^{-1}{\it St}, \end{equation}(8)where C1 = 1, C2 = 1.6, and Ceff = 0.19. As the derivation of this equation was based on the assumption that all dust grains have an equal size corresponding to one Stokes number St, in our implementation we take \hbox{${{\it St} = \bar{\it St}}$}, a mass-weighted average Stokes number of all the particles. Parameter η is related to the gas midplane pressure Pg gradient η=12ρgrΩK2Pg∂r·\begin{equation} \label{eta} \eta = \frac{1}{2\rho_{\rm g}r\Omega_{\rm K}^2} \frac{\partial P_{\rm g}}{\partial r}\cdot \end{equation}(9)It is useful to parametrize the pressure gradient as Π=|η|vKcs,\begin{equation} \label{pi} \Pi = \frac{|\eta| v_{\rm K}}{c_{\rm s}}, \end{equation}(10)the ratio of maximum radial drift speed | η | vK, with vK being the orbital velocity, and the isothermal sound speed cs.

Bai & Stone (2010a,c, henceforth BS10) performed local 2D and 3D numerical simulations of particles and gas dynamics in the midplane of a protoplanetary disk using the Athena code, including gas as a fluid and dust as superparticles (Bai & Stone 2010b). Both the aerodynamic coupling and particles to gas feedback were included. Magnetic forces and self-gravity were ignored. They focused on a laminar disk, where the turbulence is initially not present. They regulated the aggregate sizes by varying the minimum and maximum Stokes number of particles and by assuming that each logarithmic particle size bin is represented by the same amount of mass. They found that the settling of particles triggers the streaming instability and the related turbulence maintained the particles height before the Kelvin-Helmholtz instability could emerge. The strength of this turbulence is comparable with the results of Takeuchi et al. (2012) represented by Eq. (8).

The most interesting property of the streaming instability in the planet formation context is its ability to concentrate dust particles in dense clumps. Bai & Stone (2010a,c) performed a parameter study, varying dust sizes, the dust-to-gas ratio, and the pressure gradient, in order to define threshold conditions for the strong clumping, over the Roche density, which could lead to planetesimal formation. The grain size distributions used by BS10 had the minimum Stokes number of 10-4 and the maximum of 1. Analyzing the results of their runs, in particular the dense clumps composition, they concluded that the smallest particles needed to trigger such a strong clumping correspond to a critical Stokes number of Stcrit = 10-2 and that a super-solar dust abundance is required. They found that a higher vertically integrated dust-to-gas ratio Ztot lowers the threshold abundance of large grains, as a high dust-to-gas ratio reduces the turbulence and makes it easier to form a dense midplane layer. They also found that the lower the radial pressure gradient, the more easily the streaming instability can be triggered, consistent with findings of Johansen et al. (2007).

The BS10 work suggests a critical total metallicity as a criterion for strong clumping and they find its value to be in the range 0.020.07 for different values of the pressure gradient and different sizes of dust grains. They assumed a flat mass distribution in logarithmic size bins in their models. In our runs, the size distribution is an outcome of the Monte Carlo dust coagulation modeling. In order to use the BS10 results in our code, we build a model relying on splitting the dust mass distribution in two parts: particles larger and smaller than the size corresponding to the critical value of the Stokes number Stcrit = 10-2. In our model, we calculate the metallicity Z(St> 10-2) taking into account only the large grains, with Stokes number higher than Stcrit = 10-2, Z(St>10-2)=St>10-2Σd(St)Σg,\begin{equation} \label{zcritdef} Z({\it St}>10^{-2}) = \sum_{{\it St}>10^{-2}}\frac{\Sigma_{\rm d}({\it St})}{\Sigma_{\rm g}}, \end{equation}(11)where Σd(St) is dust surface density contributed by particles corresponding to a Stokes number St. If the abundance of large grains Z(St> 10-2) is higher than a critical value Zcrit, we assume that clumping over the Roche density happens and planetesimal formation is possible.

The work of BS10 suggests that the value of Zcrit depends on the total metallicity Ztot and the radial pressure gradient, which is parametrized with Π. We assume the dependence to be Zcrit=a×Ztot+b×Π+c.\begin{equation} \label{zcrit} Z_{\rm crit} = {a} \times Z_{\rm tot} + {b} \times \Pi + {c}. \end{equation}(12)In order to determine values of the parameters a, b, and c, we analyze the results presented by Bai & Stone (2010c). For each set of runs sharing the same size distribution and pressure support, we check what is the minimum total metallicity that allows planetesimal formation. Then, we calculate the metallicity contributed by the large grains in these runs (Zcrit) and compare it to the total metallicity (Ztot). By numerical fitting to the data, we find a = −0.88, b = 0.912714, and c = 0.004125. We present the fit and data extracted from Bai & Stone (2010c) in Fig. 1. With these values, we find that the absolute minimum metallicity that triggers strong clumping, for the standard value of Π = 0.05, is Zcrit = 0.026 (assuming that all the particles are large), which is significantly higher than the standard solar value of 0.01. This result is consistent with the findings of Johansen et al. (2009b).

thumbnail Fig. 1

Critical metallicity of large aggregates necessary to trigger planetesimal formation as a function of total metallicity. The points present data obtained in direct numerical simulations by Bai & Stone (2010c) for three different values of the Π parameter. We note that Bai & Stone (2010c) use the symbol Zcrit for the critical total metallicity they find for each run, whereas we use it for the metallicity of grains larger than Stcrit = 10-2. The lines show our fit to the data (Eq. (12)). The shaded region is where the metallicity of large grains would be higher than the total one, which is not possible.

Our planetesimal formation algorithm works as follows. At every advection time step of the Monte Carlo coagulation calculation, we check the mass distribution of grains produced by the interplay of settling and coagulation, and calculate the metallicity of large grains Z(St> 10-2). If the Z(St> 10-2) is higher than the threshold value of Zcrit, planetesimal formation may happen. We also check whether the dust-to-gas ratio in the midplane exceeds unity, as this is a general condition for the streaming instability to be triggered. If both the conditions are fulfilled, we remove some amount of our largest aggregates, which corresponds to an assumed mass of clump Mclump, from the dust coagulation code and refer to them as planetesimals. The evolution of the newly formed planetesimals is not included in the current version of the code.

The BS10 models did not include self-gravity, thus the formation of planetesimals was not followed. In our models, the quantity of dust grains removed in one planetesimal-forming event is estimated with the mass of a collapsing clump Mclump=ρd0×Hd3,\begin{equation} \label{mclump} M_{\rm clump}=\rho_{\rm d}^{0} \times H_{\rm{d}}^3, \end{equation}(13)where ρd0\hbox{$\rho_{\rm d}^{0}$} is the dust density in the midplane and Hd is the vertical scale-height of the dust. For MMSN at 5 AU, we get Mclump ≈ 1022 g, which corresponds to 100 km-sized planetesimals with the internal density of 1 g cm-3, consistent with constraints from the asteroid belt (Morbidelli et al. 2009). Equation (13) presents a very crude order-of-magnitude estimate based on the typical height of the dust layer, which results from the interplay of settling and the turbulence driven by the streaming instability (Eqs. (7) and (8)). In reality, the streaming instability forms elongated filaments that fragment to form planetesimals. Recently, Yang & Johansen (2014) estimated the width of the filaments from direct numerical simulations to be on the order of 10-2 × Hg, which is consistent with our Hd estimate (for St = 10-2 and α = 10-6). In their simulations, the filaments include dust of density ~102×ρg0\hbox{$10^2\times \rho_{\rm g}^{0}$}, where ρg0\hbox{$\rho_{\rm g}^{0}$} is the gas density in the midplane. As in our models ρd0ρg0\hbox{$\rho_{\rm d}^{0}\ga\rho_{\rm g}^{0}$}, the Mclump estimated by Yang & Johansen (2014) would be by 2 orders of magnitude higher than calculated from Eq. (13). Nevertheless, we find that the final outcome of our models does not depend on the exact value of the Mclump as long as there are no resolution problems coming into play.

To avoid the resolution difficulties, we take care that the mass represented by one swarm Mswarm (Eq. (1)) is lower than the mass of the collapsing clump Mclump, meaning that we remove more than one representative particle to account for collapse of one clump. Thanks to this, the amount of dust removed in one step is not too high. At the same time, we need the Mswarm to be higher than the maximum mass of aggregate that can be produced with coagulation, which can be estimated thanks to Eq. (14). The representative particles approach only works when the physical particles are less massive than the mass of one swarm, because of the assumption that one swarm represents multiple physical particles.

3. Preliminary estimates

We investigate whether the dust coagulation can produce aggregates that are large enough to trigger strong clumping in the streaming instability, which can then lead to planetesimal formation, and model the planetesimal formation when it is possible. We present results of our numerical simulations in Sect. 4, but first we motivate our choice of parameter space with simple estimates.

The size of dust grains that can be obtained by coagulation is limited by numerous effects, such as collisional physics and radial drift. On the other hand, existence of the large grains is crucial for the streaming instability and subsequent planetesimal formation.

The maximum size of aggregates depends mainly on the critical velocity above which particles do not stick. Fragmentation velocities of silicate aggregates are measured in laboratory experiments to be in the range of a few tens of cm s-1 to a few m s-1, while bouncing collisions already happen at velocities of a few cm s-1 (Güttler et al. 2010; Seizinger & Kley 2013; Kelling et al. 2014). The collision outcome is known to depend strongly on the porosity, and porous grains may grow even at velocities of 30 m s-1 (Wada et al. 2011, 2013; Meru et al. 2013b). However, numerical models including the porosity have shown that the porous silicate grains will be collisionally compacted and the growth will be halted by bouncing (Zsom et al. 2010). The dust grains consisting of ices are considered to be significantly more sticky and resistant to compaction (Okuzumi et al. 2012; Kataoka et al. 2013b; Aumatell & Wurm 2014). However, as the laboratory experiments involving ices are more challenging than for silicate grains, there is still no detailed collision model for such grains. In the molecular dynamics models, ice grains are found to be very porous and able to grow even at velocities of 50 m s-1 (Wada et al. 2009).

To account for the difference in growth efficiency of silicates and ices in our simple model, we assume a different critical velocity for growth inside and outside the snow line. For the silicate particles present inside the snow line (located at 3 AU in this model) we assume an impact velocity for bouncing/fragmentation of vin = 10 cm s-1 and for icy particles vout = 10 m s-1.

As we place our models in a dead zone, we do not consider turbulence to be a source of impact velocities. The impact velocities are thus determined by relative drift, which can be parametrized by η (Eq. (9)). Besides the bouncing and fragmentation, the maximum size of grains can also be restricted by removal of material with radial drift. However, we find that for the enhanced dust abundance, which is a prerequisite for an efficient streaming instability, the growth rate is enhanced as well, and the drift barrier does not influence the maximum size of grains, as this would only occur for particles larger than the size limited by fragmentation. We neglect the removal of material by the radial drift in this paper.

The maximum Stokes number of aggregates resulting from the collisions driven by relative drift can be estimated as (Birnstiel et al. 2012) Stmaxvf|η|vK·\begin{equation} \label{maxst} {\it St}_{\rm max}\approx\frac{v_{\rm f}}{|\eta| v_{\rm K}}\cdot \end{equation}(14)We plot the size of aggregates corresponding to the Stmax and the critical value of Stcrit = 10-2 in the laminar MMSN disk in Fig. 2. In other words, Fig. 2 shows where the dust growth can produce aggregates that are large enough to trigger the streaming instability. We find that to obtain grains with St> 10-2, the velocity vf has to be typically higher than ~1 m s-1. Thus, obtaining the particles of St> 10-2 is very hard inside the snow line, where the growth of silicate particles is halted by bouncing. However, it should be relatively easy in regions where solid ice can exist. The snow line location is fixed at 3 AU in the simple model presented in this section. In a realistic disk, the snow line migrates with time (Davis 2005; Min et al. 2011; Martin & Livio 2012; Bitsch et al. 2014). If the snow line moves inward, the region where planetesimals can form extends.

thumbnail Fig. 2

Comparison of maximum dust particle size produced by coagulation (red shaded region) and minimum size needed for planetesimal formation in a dead zone (blue crosshatched region). We assume that particle growth is limited by the relative drift velocities, ignoring turbulence, which corresponds to the dead zone. We find that the maximum particle size exceeds the size corresponding to the St = 10-2 only for the ices that can exist beyond the snow line, where the presence of ice makes the dust particles more sticky.

thumbnail Fig. 3

Time evolution of the dust size distribution in our fiducial run with and without the planetesimal formation by the streaming instability enabled. The vertical dashed line indicates grain size corresponding to a Stokes number of 10-2. Larger grains can form planetesimals and be removed from the simulation.

4. Results

We use our dust evolution code together with the planetesimal formation prescription described in Sect. 2 to model dust coagulation and planetesimal formation in a dead zone of the MMSN disk. Motivated by the estimates presented in the previous section, which show that the streaming instability can only form planetesimals outside the snow line, we locate our numerical models at 5 AU, where the cores of giant planets in the solar system were presumably formed. We assume that the dust aggregates have internal density of ρp = 1 g cm-3 and we treat them as compact spheres. The dust grains have an initial size of 1 μm and are distributed such that the dust-to-gas ratio is constant within the whole vertical range. We let the grains stick, fragment, settle down toward the midplane, and be stirred by turbulence triggered by the streaming instability when the dust-to-gas ratio in the midplane reaches unity. Fragmentation occurs for collisions with impact speeds higher than a critical value vf. All our runs have vertical resolution of 100 grid cells and we place 400 representative particles in each cell. For this level of resolution, we only find minor differences between runs started with identical parameters but different random seeds, as is usually practiced for Monte Carlo methods.

4.1. Fiducial run

For our fiducial run we choose Π = 0.05, corresponding to a pressure gradient slightly reduced with respect to the nominal MMSN model, where Π ≈ 0.08 at 5 AU. However, this value matches a more realistic disk model presented by Chiang & Youdin (2010). In general, the value of Π increases with the radial distance from the star. In the MMSN model Π ≈ 0.055 × (r/ AU)1/4, whereas in the Chiang & Youdin (2010) model Π ≈ 0.036 × (r/ AU)2/7. We start the run with the vertically integrated dust-to-gas ratio of Z = 0.05, a factor of five higher than the usual solar metallicity. For the impact velocity above which the particles fragment, we take vf = 1000 cm s-1. A setup of this kind could correspond to a pressure trap induced by a long-lived zonal flow (Dittrich et al. 2013).

Figure 3 shows the evolution of grain size distribution in the fiducial run. For comparison, we also show the evolution of the same run but without the planetesimal formation algorithm enabled. A typical sedimentation-driven coagulation scenario happens: the equal-sized particles initially grow slowly thanks to the Brownian motion. Then, the particles in the upper layers start to grow much faster than those in the midplane, because the settling velocities increase with height (Eq. (6)). The largest aggregates sweep-up smaller particles while they settle down and thus further increase their settling velocity, resulting in formation of a bimodal size distribution at ~600 yrs. Then, the particles encounter the fragmentation barrier and a coagulation-fragmentation equilibrium develops, leading to a power-law-like size distribution. The slope of the distribution depends on mass distribution of fragments implemented. We implement the fragment distribution n(m) dmm− 11/6 dm, corresponding to the MRN distribution (Mathis et al. 1977). We find that the coagulation-fragmentation equilibrium in the dead zone, where collisions are mainly driven by the systematic drift, leads to the size distribution n(a)·a·mdlogaa1/2dloga,\begin{equation} \label{sizedistr} n(a)\cdot a\cdot m\ {\rm d}{\log}\,a \propto a^{1/2}\ {\rm d}{\log}\,a, \end{equation}(15)visible in the bottom panel of Fig. 3.

The evolution proceeds identically in both cases: with and without planetesimal formation, until ~1000 yrs. After this time, both the conditions described in Sect. 2.2 are fulfilled and some particles are removed to account for the planetesimal formation. Removing the large grains slows down the collisional evolution, which can be seen in the fourth panel of Fig. 3, because the missing large grains possess very high interaction rates and would normally participate in many collisions. However, this effect is quickly smeared out by the development of the coagulation-fragmentation equilibrium. The final difference between the size distributions is not very pronounced.

thumbnail Fig. 4

Time evolution of our fiducial run: a) total planetesimals and dust mass; b) ratio of the surface density of particles larger than St = 10-2 to the gas density compared to the threshold value Zcrit; and c) dust-to-gas ratio in the midplane compared to the threshold value of unity. The dashed vertical line corresponds to the time when the first planetesimals are formed. It can be seen that although the metallicity of large grains was already higher than Zcrit at t ~ 700 years b), the planetesimal formation only starts at t ~ 990 years a) because the large grains have not settled to the midplane yet c).

Figure 4 presents an extended overview of the fiducial run time evolution. Panel a) presents how the dust and planetesimals abundance changes with time. Planetesimal formation starts at t ~ 990 years and lasts ~6 × 104 years. Planetesimal formation efficiency decreases over time. Panels b) and c) correspond to the two conditions that have to be fulfilled simultaneously to allow planetesimal formation: high metallicity in the form of large grains (Z(St> 10-2) >Zcrit) and the midplane dust-to-gas ratio higher than unity. In the presented case, the first condition is already fulfilled at t ~ 700 years, but the settling takes an additional ~300 years which delays the onset of planetesimal formation. The large quantity of large grains already present results in a planetesimal formation outburst: around 40% of the final number of planetesimals is formed in the first planetesimal-forming step. After this, the turbulence generated by the streaming instability stirs the midplane layer below the threshold dust density again. Interplay of growth and settling leads to refilling the reservoir of large grains. Each time a clump is removed, the total metallicity Ztot decreases and the threshold metallicity Zcrit increases (Eq. (12)). After ~6 × 104 years, Zcrit is so high that the coagulation cannot produce enough large grains and planetesimal formation is no longer possible.

The fact that the sedimentation takes longer than growth up to the required sizes is not a general rule. In some of the other runs, presented in the following section, we observe the opposite relation, i.e., that the growth takes longer than settling. This occurs because both timescales are comparable in the dead zone case. Thanks to this, the sedimentation-driven coagulation may happen, where the growth and settling timescales for the rain-out particles are equal.

It is worth noting that the modeling of sedimentation-driven coagulation requires the vertical structure to be included. The sedimentation-driven coagulation determines the initial timescale of growth. The process we describe cannot be modeled in detail using vertically integrated algorithms, as these assume that equilibrium between vertical settling and mixing is reached on a timescale shorter than the growth timescale. We find that planetesimal formation starts before this equilibrium is reached.

4.2. Parameter study

In order to check how our results depend on the choice of parameters, we vary the fragmentation velocity vf, the pressure gradient Π, and the total vertically integrated metallicity Ztot. We perform 24 runs in total. Table 1 gives an overview of the parameter values used in our simulations and their results in terms of how much dust is turned into planetesimals and how long we have to wait for the planetesimal formation to be triggered. Our fiducial run is referred to as P5Z5V1000 in this table. Different panels of Fig. 5 show how the planetesimal formation changes when one of the parameters is varied from its fiducial value. We discuss the parameter dependencies further in this section.

Table 1

Overview of parameters used and results obtained in different runs.

thumbnail Fig. 5

Planetesimal formation in our runs depends on parameters: a) radial pressure gradient parametrized by Π (Eq. (10)); b) total metallicity Ztot; c) fragmentation velocity vf. The red line in each panel corresponds to our fiducial run P5Z5V1000.

4.2.1. Dependence on the pressure gradient

Both Johansen et al. (2007) and Bai & Stone (2010c) observed that with a lower pressure gradient it is easier to trigger the planetesimal formation. This effect may seem counterintuitive, as the streaming instability is driven by the existence of the pressure gradient and the sub-Keplerian rotation of gas. However, the radial pressure gradient provides free energy that triggers the turbulence; therefore, a lower pressure gradient weakens the turbulence and promotes particle settling to higher densities (see Eq. (8)). In our model, the radial pressure gradient is parametrized by Π (Eq. (10)). Based on the BS10 results, the critical metallicity for planetesimal formation, Zcrit, is dependent on Π such that a lower pressure gradient promotes planetesimal formation (Eq. (12)). What is more, with a lower Π, the collision speeds between particles are also lower and thus the maximum particle size and the corresponding Stokes number (Stmax) are higher. This is because, as seen in Eq. (14), the maximum Stokes number we can obtain is inversely proportional to the maximum radial drift velocity, which is proportional to the radial pressure gradient. With the two effects combined, the effect of the radial pressure gradient is very strong. As can be seen in panel a) of Fig. 5, lowering the pressure gradient by 60% results in a greater than 80% increase in the number of planetesimals that can be formed. On the other hand, by increasing the gradient by 60%, we get a dramatic decrease in the number of planetesimals produced: from 40% to only 5% of the initial mass of dust.

4.2.2. Dependence on the metallicity

Increasing the total amount of dust available has a positive effect on the number of planetesimals that are formed. First of all, higher total metallicity directly translates into a larger quantity of large grains (Z(St> 10-2)) that can participate in the strong clumping that can lead to planetesimal formation. As seen in Eq. (12), higher total dust abundance lowers the critical vertically integrated dust-to-gas ratio Zcrit because the higher dust abundance reduces the turbulence triggered by the instability (Eq. (8)) and thus makes the clumping easier. The higher dust-to-gas ratio in the midplane also leads to a faster rotation of the gas and reduction of the difference between the gas and Keplerian rotation (Bai & Stone 2010a). This has a similar effect to reducing of the radial pressure gradient, which facilities the planetesimal formation, as discussed above. Finally, the higher metallicity speeds up the growth and the planetesimal formation is able to start earlier, as the growth timescale scales inversely with the vertically integrated dust-to-gas ratio (Birnstiel et al. 2012), τgrowth(ZtotvK)-1,\begin{equation} \tau_{\rm growth}\propto \left(Z_{\rm tot}v_{\rm K}\right)^{-1}\!, \end{equation}(16)where vK is the Keplerian velocity. We observe that with metallicity Ztot> 0.03, the growth proceeds so quickly that it is the settling timescale which determines the onset of planetesimal formation. This effect takes place in our fiducial run P5Z5V1000. With these combined positive effects, metallicity strongly influences the number of planetesimals formed. However, the impact of metallicity is nonlinear: increasing it by 40% results in >40% increase in the number of planetesimals formed, but decreasing it by 40% with respect to the fiducial value results in almost no planetesimals being formed (panel b of Fig. 5). The strong dependence on the vertically integrated dust-to-gas ratio is consistent with the results of direct numerical simulations (Johansen et al. 2009b; Bai & Stone 2010a).

4.2.3. Dependence on the fragmentation velocity

Dependence on the fragmentation velocity is the most straightforward. The value of vf only influences the maximum size of particles we are able to obtain by coagulation. The maximum Stokes number Stmax increases proportionally to the fragmentation velocity, as can be seen in Eq. (14). The higher the fragmentation velocity, the more grains will have a Stokes number above the critical value (Stcrit> 10-2), and it is easier to reach the condition of the high metallicity of large grains (Z(St> 10-2) >Zcrit). We present the dependence of our results on the fragmentation velocity in panel c) of Fig. 5. For the chosen values, the vf does not influence the outcome as strongly as the other two parameters. However, as we show in Sect. 3, with even lower fragmentation velocities we would not be able to produce any grains larger than the critical size. As most of the mass is found in large grains, which is connected to the fragmentation law we chose, we conclude that the fragmentation velocity (once it has permited growth of grains with a Stokes number St> 10-2) does not have a very strong influence on the number of planetesimals that can be formed.

4.3. Analytical model

We presented results of our numerical simulations following interplay between coagulation and planetesimal formation under different conditions in the disk. In this section, we construct a relatively simple formula that estimates the results and compute the efficiency of planetesimal formation in terms of the pressure gradient, initial vertically integrated dust-to-gas ratio, and fragmentation velocity.

As we describe in the previous sections, the most important property that determines whether planetesimal formation is possible is the quantity of dust grains with Stokes number higher than the critical value Stcrit = 10-2. If the Stokes number of the largest grains that can be produced by coagulation is Stmax>Stcrit, the quantity of grains of interesting sizes can be estimated taking into account the dust mass distribution. The mass distribution we obtain when the coagulation and fragmentation are in equilibrium is described with Eq. (15). We emphasize that this particular slope is a result of the chosen distribution of fragments, which was described in Sect. 4.1.

Taking into account that mSt3 (in the Epstein drag regime, which applies to the small grains that we obtain in our simulations) and dlog  m = dm/m, we can rewrite the size distribution (Eq. (15)) in terms of the Stokes number n(St)·mdStSt1/2dSt.\begin{equation} \label{stokesdistr} n({\it St})\cdot m\ {\rm d}{\it St} \propto {\it St}^{-1/2}\ {\rm d}{\it St}. \end{equation}(17)The relative amount of dust above a critical Stokes number Stcrit is Z(St>10-2)Ztot=StcritStmaxn(St)mdStStminStmaxn(St)mdSt=StmaxStcritStmaxStmin·\begin{equation} \label{ratio} \frac{Z({\it St}>10^{-2})}{Z_{\rm tot}} = \frac{\int_{{\it St}_{\rm crit}}^{{\it St}_{\rm max}} n({\it St})m\ {\rm d}{\it St}}{\int_{{\it St}_{\rm min}}^{{\it St}_{\rm max}} n({\it St})m\ {\rm d}{\it St}} = \frac{\sqrt{{\it St}_{\rm max}}-\sqrt{{\it St}_{\rm crit}}}{\sqrt{{\it St}_{\rm max}}-\sqrt{{\it St}_{\rm min}}}\cdot \end{equation}(18)The minimum Stokes number Stmin is in principle determined by the size of monomers, and for μm-sized grains Stmin ≈ 10-6. The maximum Stokes number Stmax is determined by the fragmentation velocity vf and the impact speeds. In the dead zone case, collisions are mainly driven by radial, azimuthal and vertical drift and the maximum size of grains is estimated accurately by Eq. (14). We assume the critical Stokes number of Stcrit = 10-2, as discussed in Sect. 2.2.

The planetesimal formation happens as long as the abundance of large grains is higher than critical (Z(St> 10-2) >Zcrit). The threshold abundance Zcrit is dependent on the pressure gradient Π and the total metallicity Ztot, as described with Eq. (12). When planetesimals are produced, the total metallicity decreases and thus the critical metallicity increases. The planetesimal formation stops when the quantity of large grains, corresponding to Z(St> 10-2) cannot reach the threshold value Zcrit anymore. Taking into account the ratio of the mass included in large grains to the total mass resulting from the equilibrium mass distribution (Eq. (18)), we can derive a critical total vertically integrated dust-to-gas ratio Ztot,crit below which planetesimal formation is no longer possible. Using Eqs. (12) and (18), we obtain Ztot,crit=(bΠ+c)×(StmaxStcritStmaxStmina)-1·\begin{equation} Z_{\rm tot,crit} = \left(b\Pi + c\right)\times\left({{ \dfrac{\sqrt{{\it St}_{\rm max}}-\sqrt{{\it St}_{\rm crit}}}{\sqrt{{\it St}_{\rm max}}-\sqrt{{\it St}_{\rm min}}} }-a}\right)^{-1}\cdot \end{equation}(19)The values of a, b, and c are given in Sect. 2.2, and the values of the Stokes numbers are discussed under Eq. (18).

The end of the planetesimal formation phase happens when the total metallicity Ztot has dropped below the critical total metallicity Ztot,crit and the relative number of planetesimals produced Mplts/Mtot can be estimated as MpltsMtot=Ztot,0Ztot,critZtot,0,\begin{equation} \label{mpltslast} \frac{M_{\rm plts}}{M_{\rm tot}} = \frac{Z_{\rm tot,0}-Z_{\rm tot,crit}}{Z_{\rm tot,0}}, \end{equation}(20)where Ztot,0 is the initial total metallicity.

thumbnail Fig. 6

Comparison of the mass in planetesimals as a function of total metallicity obtained in different runs (points) and predicted by our analytical model (lines). The three different colors correspond to different values of the pressure gradient parameter Π. The three lines of the same colors correspond to different values of the fragmentation velocity used. For Π = 0.08 and vf = 100 cm s-1, there is no planetesimal formation predicted, thus the line is not visible.

Figure 6 shows a comparison of the relative planetesimal mass produced, Mplts/Mtot, obtained with Eq. (20) and measured in our simulations (Table 1). The efficiency of planetesimal formation given by Eq. (20) tends to slightly underestimate the results of our simulations. We associate this with the mass distribution peak around maximum grain size visible in the bottom panel of Fig. 3, which is not covered by the simple power law (Eq. (17)) that we used to derive Eq. (20). This peak arises because of the fragmentation barrier: particles that grow mostly in roughly similar-sized collisions, suddenly lack slightly larger collision partners. In order to sustain a steady state, the particles pile up around the maximum size (Birnstiel et al. 2011).

5. Discussion

The work we have presented includes some inevitable assumptions. In this section, we discuss uncertainties contributed by these assumptions as well as possible ways to improve our work.

The collision model used in our work is highly simplified, restricting all the collision physics to one crucial parameter: the fragmentation velocity vf. What is more, we assume that all the dust grains are spherical and compact, with a constant internal density. It is known that the internal structure of grains is important to the coagulation (Ormel et al. 2007; Okuzumi et al. 2009; Meru et al. 2013b) and dynamics (Okuzumi et al. 2012; Hubbard & Ebel 2014). However, we lack reliable models of icy aggregate porosity, because laboratory experiments including the ices are particularly challenging. Our work may be improved in the future, as soon as laboratory data on collisions of icy aggregates are available.

Our streaming instability model relies on the work of BS10 and thus inherits all its uncertainties. The highest uncertainty is probably contributed by the fact that we focus on their 2D simulations results, as due to computational expense reasons, the extended parameter study was performed in 2D only. Bai & Stone (2010a,c) reported that the conditions for the strong clumping are more stringent in 3D than in 2D, thus our model may be a bit too optimistic. On the other hand, the BS10 models do not include the particle’s self-gravity, which might help to obtain strong clumping for lower metallicity or higher pressure gradient. This is because with the self-gravity a dust clump is able to collapse even if its initial density is lower than the Roche density (Michikoshi et al. 2010).

Bai & Stone (2010a,c) reported that their runs saturate within ~50100 orbits. This is the time the dust needs to settle and the streaming instability needs to develop. The settling timescale is self-consistently included in our models, but we neglect the time the streaming instability needs to develop: we assume the streaming instability is able to produce planetesimals immediately after the dust has settled and a sufficient quantity of large grains is present. We motivate this assumption by the fact that the timescale of dust coagulation is typically longer than the timescale to develop the instability. What is more, the timescale for clump collapse measured in simulations that include self-gravity is very short, on the order of one orbit (Johansen et al. 2011).

Bai & Stone (2010a,c) assumed a flat mass distribution in logarithmic bins in the Stokes number, whereas in our models the mass distribution is a result of realistic coagulation-fragmentation cycle and roughly follows the power law described with Eq. (15). The model we built to overcome the problem of incompatible size distributions relies on splitting the dust mass into particles larger and smaller than a size corresponding to the Stcrit = 10-2. We assume that the large particles actively drive the instability and participate in clumps, and fit the dependence of critical abundance of the large aggregates on the total metallicity (Sect. 2.2, Fig. 1). Both the choice of a constant value of Stcrit, and our fits should be rather treated as a first approximation. Given the limited data of the BS10 papers, this is the best we can do. This model may be improved as soon as more detailed parameter studies of the streaming instability become available. What is necessary is a more systematic parameter scan of the 3D streaming instability models, ideally going hand-in-hand with results from coagulation models to ensure that realistic particle size distributions are used in these 3D models.

The BS10 results are restricted to laminar disks, without the turbulence, unlike the work of Johansen et al. (2007). It is known that some source of viscosity is necessary as a mechanism of the angular momentum transport, which enables the disk accretion on observed timescales (Shakura & Sunyaev 1973). The magnetorotational instability (MRI) was considered as a mechanism to drive turbulence that is a source of such viscosity (Balbus & Hawley 1991). Thus, neglecting the turbulence may seem a serious restriction. However, many models of protoplanetary disks suggest that large regions of the disks can be free from the MRI (Gammie 1996; Turner et al. 2007; Dzyurkevich et al. 2013). In this picture, the midplane is quiescent, but the upper layers may still be active. The turbulence from the active layers may stir the midplane, which we neglect in this paper. Recently it was shown that the MRI is inefficient or even completely suppressed in the inner regions of disks when all the non-ideal magnetohydrodynamic effects, such as the ambipolar diffusion and the Hall effect, are taken into account (Bai 2013, 2014; Lesur et al. 2014). The typical extent of this dead zone is from 1 AU to 10 AU, which covers the planetesimal formation region we investigate. Thanks to the existence of a dead zone, not only are the planetesimals able to form, but they also avoid destructive collisions and grow by the runaway growth (Gressel et al. 2011, 2012; Ormel & Okuzumi 2013). The evolution of the planetesimals and their interaction with the remaining dust is not yet included in our model.

The extent of the dead zone is regulated by the quantity of small dust grains, whose large surface-to-mass ratio enables the sweep-up of free electrons from the gas and thus decreases the ionization (Sano et al. 2000; Ilgner & Nelson 2006). In the case of efficient growth of dust grains, the dead zone vanishes (Okuzumi 2009; Okuzumi & Hirose 2012). In our model, we neglect the dead zone evolution, but as we keep a significant number of small grains because of the fragmentation barrier, we consider this a safe approach.

We place our models at 5 AU, relying on our estimates, which show that it is not possible to grow sufficiently large grains out of silicate aggregates when the bouncing barrier is acting (Sect. 3). It was shown that the bouncing barrier can be overcome and large grains can grow when the fragmentation with mass transfer effect is taken into account (Windmark et al. 2012b; Garaud et al. 2013; Drążkowska et al. 2014). However, the impact velocity distribution that is necessary to produce the seed gains in those papers is contributed by the MRI turbulence, which is not present in a dead zone. The MRI turbulence would make it harder to trigger the streaming instability, as it does not allow a very thin midplane layer to form. Thus, the sweep-up growth scenario may be operating and forming planetesimals preferentially in the active zone of protoplanetary disk, whereas the scenario we investigate in this paper works better in the dead zone.

Another way of forming sufficiently large aggregates inside the snow line is co-accretion of dust and chondrules, suggested by Ormel et al. (2008) and subsequently investigated in laboratory experiments by Beitz et al. (2012) and Jankowski et al. (2012). In this scenario, chondrules acquire dusty rims that facilitate their sticking and allow the bouncing barrier to be overcome. As shown by Ormel et al. (2008), this scenario requires moderately low turbulence, which would also be favorable for planetesimal formation via the streaming instability.

Our models are local and neglect the radial drift. The radial drift timescale for pebbles of St = 10-2 at 5 AU is on the order of 103 orbits, which is one order of magnitude longer than the coagulation needs to trigger the streaming instability, thus we consider this a safe approach.

We assume that metallicity is enhanced over the standard solar value of 0.01. Such an enhancement could form in a radial drift dominated disk, where the pebbles necessary to trigger the streaming instability are drifting inward and may form pile-ups, as suggested by Youdin & Shu (2002), Youdin & Chiang (2004), Laibe et al. (2012), and Laibe (2014). However, such an enhancement typically happens inside the snow line, where growing sufficiently large particles is suppressed by the non-stickiness of aggregates and the existing large aggregates may be destroyed by ice evaporation and high-speed collisions. A way to overcome this problem may be the reduction of impact speeds in the presence of strong dust-to-gas ratio enhancements (Nakagawa et al. 1986; Bai & Stone 2010a), which is not yet included in our code. The reduction of impact speeds would, in general, help to form larger aggregates and thus to produce more planetesimals.

Another possibility that would justify the enhanced metallicity and the reduced pressure support at the same time is the existence of some kind of pressure bump. However, the formation process and lifetime of such pressure bumps remains uncertain. Pressure bumps caused by zonal flows have been observed in numerical simulations including the MRI turbulence (Johansen et al. 2009a; Dzyurkevich et al. 2010; Uribe et al. 2011), but their lifetimes are up to 50 orbits (Dittrich et al. 2013). A pressure bump arising around the snow line was suggested by Kretke & Lin (2007), but it was recently found to require an unrealistically high viscosity transition (Bitsch et al. 2014). What is more, complicated evaporation and condensation processes have to be taken into account to model the size evolution of dust (Kuroiwa & Sirono 2011) when considering the region near the snow line. Thus, our model is not self consistent in the pressure bump and dust enhancement origin aspect: we start our runs with metallicity already enhanced by a factor of a few with respect to a nominal solar value, as otherwise the planetesimal formation by streaming instability is not possible. We investigate the interplay between vertical settling, coagulation, and planetesimal formation, ignoring the radial drift of dust. Such an enhancement would in reality build up over a timescale determined by the radial drift, even in the presence of a pre-existing pressure bump. We plan to investigate the effects of radial drift on our results in a future work.

The sizes of clumps and planetesimals created by the streaming instability are highly uncertain, as the hydrodynamic simulations have limited resolution and are typically not able to follow the clump collapse with realistic collisional behavior. As described in Sect. 2.2, we assume that the collapsing clumps have identical mass, which is estimated based on the height of the dust layer (Eq. (13)). This is only an order-of-magnitude estimate, which is not necessarily consistent with the recent results of Yang & Johansen (2014), as discussed in Sect. 2.2. However, we find that, although the details of evolution are dependent on the assumed clump mass, the final outcome in terms of total planetesimal mass produced is not.

6. Conclusions

The streaming instability was proposed as an efficient way of overcoming the growth and drift barriers and forming planetesimals. However, strong clumping was proven to require large grains. In this paper, we investigated whether large grains can form in sufficient amount during coagulation under realistic conditions. We developed and implemented a simplified model for planetesimal formation in our dust coagulation code, as described in Sect. 2. Our work is a step toward a unified model for planetesimal formation because we join the dust coagulation modeling with planetesimal formation via streaming instability for the first time.

We find that planetesimal formation by streaming instability is hindered for the silicate aggregates because the bouncing barrier prevents growth to the sizes (Stokes numbers) needed to trigger the streaming instability. It is possible to obtain the minimum size of particles, corresponding to the critical Stokes number (Stcrit = 10-2), for the stickier, icy aggregates, which are present beyond the snow line. If some way can be found to overcome the bouncing barrier, we can also expect the streaming instability to operate inside the snow line. However, the strong clumping may only be triggered when the vertically integrated dust-to-gas ratio is enhanced by a factor of at least three with respect to the solar value of 0.01 and/or the radial pressure gradient is reduced with respect to the standard minimum mass solar nebula model. What is more, a dense midplane layer of solids have to be formed, which is only possible if the turbulence is relatively weak.

We modeled the interplay of dust coagulation and settling and planetesimal formation and performed a parameter study, varying the radial pressure gradient, metallicity, and fragmentation velocity. We investigated how these values influence the amount of planetesimals formed. We proposed a simple explanation of the obtained results, by constructing an analytical expression for the maximum number of dust that can be turned into planetesimals (Sect. 4.3). This model can be used in future projects, for example to constrain initial conditions for planetesimal evolution models or in planet population synthesis codes.


1

1832 known exoplanets according to the Extrasolar Planets Encyclopaedia http://exoplanet.eu on 13 October 2014.

2

In this paper, we use the term metallicity interchangeably with the vertically integrated dust-to-gas ratio.

Acknowledgments

We thank Chris Ormel, Carsten Dominik, and Alessandro Morbidelli for useful discussions and encouragement. We thank Xue-Ning Bai and the referee, Anders Johansen, for their comments that helped us to improve this paper. J.D. would like to acknowledge the use of the computing resources provided by bwGRiD (http://www.bw-grid.de), member of the German D-Grid initiative, funded by the Ministry for Education and Research (Bundesministerium für Bildung und Forschung) and the Ministry for Science, Research and Arts Baden-Wuerttemberg (Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg).

References

  1. Aumatell, G., & Wurm, G. 2014, MNRAS, 437, 690 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bai, X.-N. 2013, ApJ, 772, 96 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bai, X.-N. 2014, ApJ, 791, 137 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bai, X.-N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bai, X.-N., & Stone, J. M. 2010b, ApJS, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bai, X.-N., & Stone, J. M. 2010c, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
  7. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  8. Balsara, D. S., Tilley, D. A., Rettig, T., & Brittain, S. D. 2009, MNRAS, 397, 24 [NASA ADS] [CrossRef] [Google Scholar]
  9. Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
  10. Barranco, J. A. 2009, ApJ, 691, 907 [NASA ADS] [CrossRef] [Google Scholar]
  11. Beitz, E., Güttler, C., Weidling, R., & Blum, J. 2012, Icarus, 218, 701 [NASA ADS] [CrossRef] [Google Scholar]
  12. Benz, W., & Asphaug, E. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef] [Google Scholar]
  13. Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bitsch, B., Morbidelli, A., Lega, E., Kretke, K., & Crida, A. 2014, A&A, 570, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Blum, J., & Münch, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
  17. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, Phys. Fluids, 11, 2280 [NASA ADS] [CrossRef] [Google Scholar]
  19. Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Brauer, F., Dullemond, C. P., & Henning, T. 2008a, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Brauer, F., Henning, T., & Dullemond, C. P. 2008b, A&A, 487, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 1633 [NASA ADS] [CrossRef] [Google Scholar]
  23. Carballido, A., Bai, X.-N., & Cuzzi, J. N. 2011, MNRAS, 415, 93 [NASA ADS] [CrossRef] [Google Scholar]
  24. Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chiang, E., & Youdin, A. N. 2010, Ann. Rev. Earth Planet. Sci., 38, 493 [Google Scholar]
  26. Ciesla, F. J. 2010, ApJ, 723, 514 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102 [NASA ADS] [CrossRef] [Google Scholar]
  29. Davis, S. S. 2005, ApJ, 620, 994 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dittrich, K., Klahr, H., & Johansen, A. 2013, ApJ, 763, 117 [NASA ADS] [CrossRef] [Google Scholar]
  31. Drążkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Drążkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 [NASA ADS] [CrossRef] [Google Scholar]
  36. Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
  37. Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 [NASA ADS] [CrossRef] [Google Scholar]
  38. Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
  39. Goodman, J., & Pindor, B. 2000, Icarus, 148, 537 [NASA ADS] [CrossRef] [Google Scholar]
  40. Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291 [NASA ADS] [CrossRef] [Google Scholar]
  41. Gressel, O., Nelson, R. P., & Turner, N. J. 2012, MNRAS, 422, 1140 [NASA ADS] [CrossRef] [Google Scholar]
  42. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  44. Hubbard, A., & Ebel, D. S. 2014, Icarus, 237, 84 [NASA ADS] [CrossRef] [Google Scholar]
  45. Ilgner, M., & Nelson, R. P. 2006, A&A, 445, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
  47. Jankowski, T., Wurm, G., Kelling, T., et al. 2012, A&A, 542, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219 [NASA ADS] [CrossRef] [Google Scholar]
  49. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  50. Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
  51. Johansen, A., Youdin, A., & Mac Low, M.-M. 2009b, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
  52. Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI (University of Arizona Press), accepted [arXiv:1402.1344] [Google Scholar]
  54. Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013a, A&A, 557, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013b, A&A, 554, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Kelling, T., Wurm, G., & Köster, M. 2014, ApJ, 783, 111 [NASA ADS] [CrossRef] [Google Scholar]
  57. Klahr, H., & Bodenheimer, P. 2006, ApJ, 639, 432 [NASA ADS] [CrossRef] [Google Scholar]
  58. Kowalik, K., Hanasz, M., Wóltański, D., & Gawryszczak, A. 2013, MNRAS, 434, 1460 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kuroiwa, T., & Sirono, S.-I. 2011, ApJ, 739, 18 [NASA ADS] [CrossRef] [Google Scholar]
  61. Laibe, G. 2014, MNRAS, 437, 3037 [NASA ADS] [CrossRef] [Google Scholar]
  62. Laibe, G., Gonzalez, J.-F., & Maddison, S. T. 2012, A&A, 537, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009, A&A, 497, 869 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Martin, R. G., & Livio, M. 2012, MNRAS, 425, L6 [NASA ADS] [CrossRef] [Google Scholar]
  66. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
  67. Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Meru, F., Galvagni, M., & Olczak, C. 2013a, ApJ, 774, L4 [NASA ADS] [CrossRef] [Google Scholar]
  69. Meru, F., Geretshauser, R. J., Schäfer, C., Speith, R., & Kley, W. 2013b, MNRAS, 435, 2371 [NASA ADS] [CrossRef] [Google Scholar]
  70. Michikoshi, S., Kokubo, E., & Inutsuka, S.-I. 2010, ApJ, 719, 1021 [NASA ADS] [CrossRef] [Google Scholar]
  71. Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [NASA ADS] [CrossRef] [Google Scholar]
  72. Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558 [NASA ADS] [CrossRef] [Google Scholar]
  73. Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
  74. Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
  75. Okuzumi, S. 2009, ApJ, 698, 1122 [NASA ADS] [CrossRef] [Google Scholar]
  76. Okuzumi, S., & Hirose, S. 2012, ApJ, 753, L8 [NASA ADS] [CrossRef] [Google Scholar]
  77. Okuzumi, S., Tanaka, H., & Sakagami, M.-A. 2009, ApJ, 707, 1247 [NASA ADS] [CrossRef] [Google Scholar]
  78. Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
  79. Ormel, C. W., & Okuzumi, S. 2013, ApJ, 771, 44 [NASA ADS] [CrossRef] [Google Scholar]
  80. Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Ormel, C. W., Cuzzi, J. N., & Tielens, A. G. G. M. 2008, ApJ, 679, 1588 [NASA ADS] [CrossRef] [Google Scholar]
  82. Pagani, L., Steinacker, J., Bacmann, A., Stutz, A., & Henning, T. 2010, Science, 329, 1622 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  83. Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Ros, K., & Johansen, A. 2013, A&A, 552, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] [Google Scholar]
  86. Seizinger, A., & Kley, W. 2013, A&A, 551, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  88. Takeuchi, T., Muto, T., Okuzumi, S., Ishitsu, N., & Ida, S. 2012, ApJ, 744, 101 [NASA ADS] [CrossRef] [Google Scholar]
  89. Tilley, D. A., Balsara, D. S., Brittain, S. D., & Rettig, T. 2010, MNRAS, 403, 211 [NASA ADS] [CrossRef] [Google Scholar]
  90. Turner, N. J., Sano, T., & Dziourkevitch, N. 2007, ApJ, 659, 729 [NASA ADS] [CrossRef] [Google Scholar]
  91. Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85 [NASA ADS] [CrossRef] [Google Scholar]
  92. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
  93. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
  94. Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  96. Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
  97. Weidenschilling, S. J. 1995, Icarus, 116, 433 [NASA ADS] [CrossRef] [Google Scholar]
  98. Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 [Google Scholar]
  99. Windmark, F., Birnstiel, T., Güttler, C., et al. 2012a, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012b, A&A, 544, L16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Yang, C.-C., & Johansen, A. 2014, ApJ, 792, 86 [NASA ADS] [CrossRef] [Google Scholar]
  102. Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
  103. Youdin, A. N., & Chiang, E. I. 2004, ApJ, 601, 1109 [NASA ADS] [CrossRef] [Google Scholar]
  104. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
  105. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
  106. Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
  107. Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1

Overview of parameters used and results obtained in different runs.

All Figures

thumbnail Fig. 1

Critical metallicity of large aggregates necessary to trigger planetesimal formation as a function of total metallicity. The points present data obtained in direct numerical simulations by Bai & Stone (2010c) for three different values of the Π parameter. We note that Bai & Stone (2010c) use the symbol Zcrit for the critical total metallicity they find for each run, whereas we use it for the metallicity of grains larger than Stcrit = 10-2. The lines show our fit to the data (Eq. (12)). The shaded region is where the metallicity of large grains would be higher than the total one, which is not possible.

In the text
thumbnail Fig. 2

Comparison of maximum dust particle size produced by coagulation (red shaded region) and minimum size needed for planetesimal formation in a dead zone (blue crosshatched region). We assume that particle growth is limited by the relative drift velocities, ignoring turbulence, which corresponds to the dead zone. We find that the maximum particle size exceeds the size corresponding to the St = 10-2 only for the ices that can exist beyond the snow line, where the presence of ice makes the dust particles more sticky.

In the text
thumbnail Fig. 3

Time evolution of the dust size distribution in our fiducial run with and without the planetesimal formation by the streaming instability enabled. The vertical dashed line indicates grain size corresponding to a Stokes number of 10-2. Larger grains can form planetesimals and be removed from the simulation.

In the text
thumbnail Fig. 4

Time evolution of our fiducial run: a) total planetesimals and dust mass; b) ratio of the surface density of particles larger than St = 10-2 to the gas density compared to the threshold value Zcrit; and c) dust-to-gas ratio in the midplane compared to the threshold value of unity. The dashed vertical line corresponds to the time when the first planetesimals are formed. It can be seen that although the metallicity of large grains was already higher than Zcrit at t ~ 700 years b), the planetesimal formation only starts at t ~ 990 years a) because the large grains have not settled to the midplane yet c).

In the text
thumbnail Fig. 5

Planetesimal formation in our runs depends on parameters: a) radial pressure gradient parametrized by Π (Eq. (10)); b) total metallicity Ztot; c) fragmentation velocity vf. The red line in each panel corresponds to our fiducial run P5Z5V1000.

In the text
thumbnail Fig. 6

Comparison of the mass in planetesimals as a function of total metallicity obtained in different runs (points) and predicted by our analytical model (lines). The three different colors correspond to different values of the pressure gradient parameter Π. The three lines of the same colors correspond to different values of the fragmentation velocity used. For Π = 0.08 and vf = 100 cm s-1, there is no planetesimal formation predicted, thus the line is not visible.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.