EDP Sciences
Free Access
Issue
A&A
Volume 594, October 2016
Article Number A105
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201628983
Published online 19 October 2016

© ESO 2016

1. Introduction

Planets seem to be omnipresent in our Galaxy with most stars orbited by one or more of them (Cassan et al. 2012). However, we are a long way from a complete understanding of how these planets form. Despite recent progress in planet formation research, there is no model that can explain the growth of planets beginning with micron-sized dust grains. Most significantly, there is a disconnection between the models of the early stages of planet formation that are dealing with dust growth, and the late stages which follow the final accretion of planetary systems starting with planetesimals and embryos. The late stage models usually assume that planetesimals form rapidly with a smooth radial distribution (Armitage 2013). On the other hand, the early stage models typically end without producing any aggregates larger than cm-sized (see e.g. Testi et al. 2014). This is because of so-called growth barriers that result from the collisional physics of dust aggregates and loss of solids because of their radial drift. The dust aggregates tend to bounce and fragment at the impact speeds predicted for the protoplanetary disc environment (Zsom et al. 2010; Birnstiel et al. 2011). On the other hand, the radial drift timescale is shorter than the growth timescale for aggregates approaching centimetre sizes and the initial aggregates are removed from the disc before any larger bodies can form (Brauer et al. 2008; Birnstiel et al. 2009). We give a concise overview of the typical predictions for dust evolution in the following section.

This paper is organised as follows. We give a brief introduction to the problem of growth barriers in Sect. 1.1 and highlight the streaming instability as a possible solution in Sect. 1.2. In Sect. 2, we explain the numerical approach that we employ to investigate dust evolution and planetesimal formation in a viscously evolving protoplanetary disc. We present results of our models in Sect. 3. Finally, we discuss limitations of our approach in Sect. 4 and conclude the work in Sect. 5. For the readers convenience, the symbols used throughout this paper are summarised in Table B.1.

thumbnail Fig. 1

Overview of the growth barriers in our fiducial disc model. The red line corresponds to a Stokes number of unity. The yellow region marks the fragmentation barrier where the impact speeds would be higher then the threshold value of 10 m s-1. The grey region shows the size scale for which the drift timescale is shorter than the growth timescale, meaning that the aggregates would be removed from a given location faster than the growth could replenish them: this is the drift barrier. The light grey vertical lines are tracks of test particles evolving by growth and drift in a steady-state disc. Finally, the blue nearly horizontal line shows the size scale corresponding to the Stokes number of 10-2, which is a minimum size required for planetesimal formation via the streaming instability.

Open with DEXTER

1.1. Growth barriers

Dust evolution is driven by its interaction with the sub-Keplerian gas disc. The interaction between a dust aggregate and the surrounding gas can be conveniently parametrised with the so-called Stokes number (1)where ts is the timescale over which the aggregate adjusts its velocity to the velocity of surrounding gas, and ΩK is Keplerian frequency. In this way, the Stokes number is the ratio between how quickly the aggregate reacts to the gas and the orbital period. Particles with St ≪ 1 are tightly coupled to gas and St ≫ 1 indicates decoupled solids. For our needs, it is convenient to rewrite Eq. (1) as (2)where a is the radius and ρ is the internal density of the dust aggregate, and Σg is the surface density of the gas. Equation (2) applies to small and compact solids, where the radius a is smaller than the mean-free path in gas. This is the so-called Epstein drag regime, where particles are well-coupled to the gas, and this regime is valid for all models presented in this paper.

Figure 1 gives an overview of global dust evolution in a protoplanetary disc with total mass of 0.1M, in terms of distance to the central star and dust aggregate size. The size of aggregates corresponding to the Stokes number of unity, which is when they are affected by the interaction with the gas the most, is marked with the red nearly horizontal line. The blue line below shows a minimum size that is required for an efficient streaming instability (corresponding to St = 10-2, see Sect. 1.2). The outer part of the protoplanetary disc is dominated by the radial drift barrier (grey triangle-shaped region), where the growth timescale (3)is longer than the radial drift timescale (4)where ρd is the density of dust, Δv is the impact velocity for collisions between dust aggregates, r is the radial distance to the star, and vr,d is the radial drift velocity. We assume that the growth preferentially happens between equal-sized grains that have settled to the midplane when deriving Eq. (3). This means that the dust density may be written as (5)where Σd is the dust surface density and Hd is the scale height of the dust, which depends on the turbulence strength αt and the Stokes number of grains St. In the case of a single-sized population (6)where is the scale height of the gas.

We note that the growth timescale τgrowth given by Eq. (3) simplifies to (7)under the assumptions that the collisions are driven by turbulence, in which case the impact speed (Ormel & Cuzzi 2007), the dust grains are in the Epstein regime (see Eq. (2)), and the dust density in the midplane is given by Eq. (5). We emphasise that, because Δvαt1/2 and ρdαt− 1/2, the dependence of the dust growth timescale on the turbulence strength parameter αt cancels itself out.

This initial growth stage is significant particularly at large orbital distances, where the growth is much slower than in the inner parts of the protoplanetary disc. This results in a delayed delivery of pebbles from the outer disc to its inner regions. The same effect was also described by Lambrechts & Johansen (2014), who called it a “pebble formation front”.

The radial drift limits the maximum size of grains that are available in the outer parts of the disc. The pebbles that grow far away from the star are then shifted to the inner disc, where they undergo fragmentation during high-speed collisions. We assume that the fragmentation happens when the collision speed exceeds a threshold value, which we set to vth = 10 m s-1 for Fig. 1.

The evolution of dust in the outer disc is determined by the interplay between growth and drift. Because the timescale of drift is shorter than what is needed to grow to centimetre sizes at a few tens of AU from the central star, this outer region is gradually depleted on a timescale of a few Myr. However, as may be seen in Fig. 1, the radial drift barrier does not stretch all the way to the inner edge of the disc. In the inner disc, the collisional timescale is shorter than the drift timescale. If growth happens even at high impact speeds, for example in the case of very porous dust aggregates, planetesimal formation via direct growth may be possible (Okuzumi et al. 2012; Kataoka et al. 2013). Even for compact grains, interplay between the radial drift and dust growth can lead to a pile-up of solids in the inner regions of the disc (Birnstiel et al. 2012; Laibe et al. 2012; Pinte & Laibe 2014). This pile-up is required by the planetesimal formation models which include the streaming instability (Johansen et al. 2009; Bai & Stone 2010a). At the same time, the streaming instability requires the presence of relatively large pebbles, with sizes corresponding to the Stokes number of St > 10-2 (Dra¸żkowska & Dullemond 2014). This is not necessarily the case because of the bouncing and fragmentation barriers. However, as can be seen in Fig. 1, there is a region around 110 AU where the maximum size of grains that can be reached because of fragmentation is above the St = 10-2 line and the radial drift barrier is not efficient. In this paper, we check whether the redistribution of solids driven by the radial drift and growth in a realistic viscous disc can lead to planetesimal formation via the streaming instability.

1.2. Streaming instability

Owing to the growth barriers described in the previous section, direct growth from micron to km-sizes seems unlikely. Streaming instability that is able to produce planetesimals directly out of cm-sized pebbles is a good solution to the planetesimal formation issue (Johansen et al. 2007).

However, planetesimal formation via the streaming instability requires enhancement of the dust-to-gas ratio by a factor of a few over the standard solar value of 10-2 (Johansen et al. 2009; Bai & Stone 2010a). There are different scenarios that modify disc structure and introduce pressure bumps that make it possible to obtain such enhancements: the zonal flows (Johansen et al. 2011; Dittrich et al. 2013), vortices (Raettig et al. 2015; Surville et al. 2016), dead zone edges (Lyra et al. 2009), and planet-disc interactions (Ayliffe et al. 2012). In the models presented in this paper, we check whether the streaming instability can work thanks to a dust pile-up induced in the inner disc by the combination of growth and drift. This scenario does not require any ad-hoc assumptions about the disc structure or pre-existing planets.

Until now the streaming instability was only modelled in local or quasi-global simulations. We use a global 1D semi-analytical protoplanetary disc model together with a prescription for streaming instability extracted from the local hydrodynamic simulations, similar to Dra¸żkowska & Dullemond (2014) for a local case. As a consequence, we are able to identify regions of the disc in which planetesimal formation happens. We perform an extended parameter study to investigate how the planetesimal formation is affected by the disc mass, metallicity, and stickiness of dust aggregates. The modelling methods that we use are described in Sect. 2.

2. Methods

In our models, we take into account the viscous evolution of protoplanetary discs, dust growth and drift, as well as planetesimal formation via the streaming instability. We model how the surface density of gas, dust, and planetesimals evolves over timescales of up to 10 Myr. Our approach combines semi-analytical models of viscous protoplanetary disc evolution that follows the work of Alibert et al. (2005, 2013), dust evolution model based on the approach proposed by Birnstiel et al. (2012), and planetesimal formation via the streaming instability similar to Dra¸żkowska & Dullemond (2014). The following subsections describe each of these components.

2.1. Disc model

Our gas disc model is the same as that used in Alibert et al. (2005, 2013) and is computed using the method originally derived by Papaloizou & Terquem (1999). Here, we give a general description of the model and refer the reader to the above-mentioned papers for further details.

The evolution of the disc is computed using a 1+1D model in a two-step process. In the first step, we compute the vertical structure of the disc for each distance to the star by solving the equation of hydrostatic equilibrium, energy conservation (taking into account viscous dissipation), and radiative diffusion. For this calculation, the opacity is assumed to be that of Bell & Lin (1994), which is representative of the opacity in the interstellar medium. In the calculations present in the paper, we did not include the effect of stellar irradiation. The viscosity is assumed to follow the well-known Shakura-Sunyaev parametrisation (Shakura & Sunyaev 1973), with an αt parameter equal to 10-3. For higher αt values, it is very hard to obtain a dense midplane layer of dust that is necessary to trigger the streaming instability (see Sect. 2.3). The solution of the vertical structure equation gives the distribution of temperature, pressure and radiative flux along the vertical axis for each distance to the star and as a function of the gas surface density, as well as the mass-weighted mean viscosity. This latter quantity is finally used, in the second step, to solve the radial diffusion equation that provides the evolution of the gas surface density as a function of time. This diffusion equation was modified to include far-ultraviolet (FUV) photoevaporation following the prescription of Veras & Armitage (2004) with the mass-loss rate of wind = 10-6M yr-1. With this prescription, photoevaporation is significant particularly at large orbital distances where the gas is less bound, and the viscous spreading of the gas disc is suppressed.

Figure 2 highlights the evolution of gas in our fiducial model. The initial mass of this disc is 0.1M, which decreases to 0.02M after 10 Myr. The disc cools with time, which modifies the impact speeds of dust aggregates, as well as the pressure structure of the disc and thus the radial drift. The headwind ηvK, which determines the maximum radial drift speed, is set by the gas pressure Pg gradient (8)where ρg is the gas density and ΩK is the Keplerian frequency. The pressure gradient is negative and the negative values of velocity translate into inward drift. The maximum drift speed decreases during the evolution. This facilitates the retention of solids and growth of the dust pile-up in the inner part of the disc (see Sect. 3).

2.2. Dust evolution

Models that include both dust advection and growth, as pioneered by Weidenschilling (1997), are computationally demanding and their results may be difficult to interpret owing to their complexity. We focus on simplified models, only including the necessary physics in a semi-analytical manner since we want to gain an overview of the dust evolution process and understand its connection to the planetesimal formation stage.

Instead of directly modelling collisions between dust grains, we prescribe their growth and fragmentation using a prescription similar to the one proposed by Birnstiel et al. (2012). We describe the main aspects of this approach here and refer the interested readers to the full paper.

thumbnail Fig. 2

Time evolution of the gas surface density Σg, gas temperature, and the headwind speed ηvK in our fiducial disc model. The disc mass is initially 0.1M and it decreases to 0.02M after 10 Myrs. The disc is shallower then the Minimum Mass Solar Nebula (MMSN, red line) and thus most of its mass resides in its outer parts. The headwind, which determines the maximum radial drift velocity, is not constant throughout the disc, as opposed to the MMSN model, and it gradually decreases during the evolution, facilitating the retention of solids.

Open with DEXTER

The approach is based on solving the advection-diffusion equation for dust surface density Σd(9)where r is the distance to the central star, Σg is the surface density of gas, and Dg is gas diffusivity. The dust advection speed is a mass-weighted average (see Eq. (A.9)). The advection velocity corresponding to each Stokes number is described by Eq. (A.3). This consists of two components: one corresponding to the difference between gas and Keplerian rotation and one related to the radial velocity of gas. The radial drift speed depends on the maximum size of dust aggregates and on the size distribution. The size distribution is regulated by the processes that limits growth at given location: fragmentation or radial drift (see Sect. 1.1).

If the maximum size of dust aggregates is limited by the radial drift, then we assume that the size distribution is very narrow and focussed around (10)where vK is the Keplerian velocity, cs is the sound speed, and fd = 0.55 is a model parameter calibrated by Birnstiel et al. (2012). In the inner part of the disc, the size distribution is regulated by the equilibrium between coagulation and fragmentation (Birnstiel et al. 2011) and it takes the form of a power law described by Eq. (A.8). In this case, we assume that the maximum size is (11)where vth is the fragmentation threshold velocity and the calibration factor ff = 0.37. In all the models presented in this paper, fragmentation is driven by the turbulent velocities rather than differential drift. If the αt parameter would be very low, i.e. in a disc with dead zone, then Eq. (11) would have to be modified to take the differential drift into account. The sizes indicated by Eqs. (10) and (11) correspond to the lower boundaries of the growth-barrier regions depicted in Fig. 1.

We assume that at the beginning of the evolution all the dust grains have a size of a0, independent of orbital distance. For this paper, we adopt a0 = 1μm. The dust grows until it reaches the size limited by fragmentation or radial drift. The growth is prescribed such that (12)where the first condition accounts for the initial growth phase and the growth timescale τgrowth is described by Eq. (7).

As it was demonstrated by Birnstiel et al. (2012), this type of simplified approach is able to fit the results of full dust evolution simulations surprisingly well, while greatly reducing the computational cost.

One major difference between the code of Birnstiel et al. (2012) and ours is that we take into account the so-called collective effects, or the effects of backreaction of dust on gas, while calculating the radial drift velocity . Bai & Stone (2010a) noticed that the radial drift is modified when clumping occurs in their streaming instability models. The rate of radial drift decreases in a way that depends both on the local dust-to-gas ratio and on the size distribution of grains. This is because the different dust species are indirectly coupled owing to their interaction with gas. We implement these effects using equations derived by Tanaka et al. (2005), which are also found in Okuzumi et al. (2012). Details of our implementation are presented in Appendix A.

This modification of radial drift significantly strengthens the effects of dust overdensities. If there is a location where the dust-to-gas ratio is enhanced, the drift slows down and a pile-up of material arises. We demonstrate the importance of the backreaction on the global dust retention in Fig. 5 (the black dashed line versus the grey solid line).

In contrast to Birnstiel et al. (2012), who employed an implicit integration scheme, we perform explicit integration of Eq. (9) using a total variation diminishing scheme and employing the input provided by the gas evolution module. We limit the time step to (13)where dr is the width of radial grid cell, and C< 1 to assure numerical accuracy. Explicit integration is necessary to add planetesimal formation into this model, as is described in the following section.

2.3. Planetesimal formation

As mentioned in Sect. 1.2, we account for planetesimal formation via the streaming instability using a semi-analytic model similar to Dra¸żkowska & Dullemond (2014). However, they focussed on local models and studying the connection between dust growth and planetesimal formation. The planetesimals were formed out of the reservoir of pebbles that has been replenished by growth and the loss of pebbles because of radial drift was neglected. In this paper, we extend this approach to a global disc model, including both dust growth (albeit in a simplified fashion as discussed in the previous section) and the radial drift. The radial drift has both positive and negative effects on planetesimal formation: it makes planetesimal formation possible by creating pile-ups of dust in the inner disc and then it delivers the pebbles from the outer disc, extending the planetesimal formation period. However, it also gradually removes solids from the disc, restricting the amount of planetesimals that can be formed.

Another major difference between our approach and that of Dra¸żkowska & Dullemond (2014) is that they focussed on dead zones, where turbulence was only present when triggered by the midplane instability caused by dust sedimentation. In this paper, we focus on a viscous disc model, where external turbulence is present and it is parametrised by αt = 10-3. This causes a change to the dominant condition determining the possibility of planetesimal formation. Dra¸żkowska & Dullemond (2014) defined two conditions that have to be simultaneously fulfilled to allow for planetesimal formation:

  • the vertically integrated dust-to-gas ratio of pebbles with sizes corresponding to St> 10-2 exceeds a critical value Zcrit;

  • the dust-to-gas ratio in the midplane exceeds unity.

The value of Zcrit was calibrated on the hydrodynamical simulations results presented by Bai & Stone (2010a,b) who considered a laminar disc. Similar results were presented recently by Carrera et al. (2015), who also included dependence of the Zcrit on the Stokes number of particles. Here, with the external turbulence present, the vertically integrated dust-to-gas ratio Σd/ Σg has to be much higher than Zcrit used in previous work to form a dense midplane layer of solids. In other words, the second condition pointed out by Dra¸żkowska & Dullemond (2014) is always more restrictive than the first one. This is why we do not use the Zcrit, but rather focus on the midplane dust-to-gas ratio in this paper. Following Eq. (6), the relation between the vertically integrated and the midplane dust-to-gas ratio in the case of a single dust size is (14)where ρd and ρg are, respectively, the dust and the gas midplane density.

The planetesimal formation algorithm we implement works as follows. In each timestep and in each grid cell, we check whether the condition for planetesimal formation is fulfilled, namely that the density of sufficiently large pebbles in the midplane is comparable to the gas density (15)If this is the case, we transfer part of the surface density of dust into planetesimals, such that (16)where TK is the orbital period and ζ is the planetesimal formation efficiency, that is the ratio of sufficiently large pebbles that are turned to planetesimals within one orbital period. We test different values of ζ in Sect. 3.5.

3. Results

3.1. Fiducial model

For our reference run, we choose a protoplanetary disc of mass Mdisc = 0.1M, global dust-to-dust ratio of Z = 10-2, and turbulence strength parametrised by αt = 10-3. We set the fragmentation threshold to vth = 10 m s-1, independent of orbital distance, and the planetesimal formation efficiency parameter ζ = 10-4. As our initial condition, we assume a constant Σd/ Σg = Z for the radial distance between 0.3 AU and 100 AU and Σd = 0 otherwise.

Figure 3 shows the evolution of dust sizes restricted by growth barriers during our fiducial run. The combined action of dust growth and drift leads to the gradual loss of solids from the disc. Because of this, the importance of drift barriers is modified and the radial drift barrier strengthens with time (the grey triangle-shaped region moves towards the star). The fragmentation barriers slightly weakens because the gas temperature drops (see Fig. 2). The representative size of dust follows the restrictions imposed by growth barriers, except for the beginning of the simulation, when it is determined by dust growth.

thumbnail Fig. 3

Time evolution of the growth barriers. The plotting style corresponds to Fig. 1: the yellow region corresponds to the fragmentation barrier while the grey region indicates the drift regime. The red and blue horizontal lines denote sizes corresponding to the Stokes number of unity and 10-2, respectively. The purple line shows the maximum size of dust aggregates, which generally follows the growth barriers, except for the earliest stages of evolution, when it is restricted by the growth timescale. The radial drift barrier gains significance as the dust-to-gas ratio decreases during the evolution.

Open with DEXTER

thumbnail Fig. 4

Time evolution of a) the vertically integrated dust-to-gas ratio and b) surface density of planetesimals. The outer parts of the disc are depleted of dust because they are drift dominated. At the same time, there is a temporary pile-up of dust in the fragmentation-dominated inner disc, which enables planetesimal formation via the streaming instability. The dashed line in panel b) corresponds to the surface density of solids in the MMSN model.

Open with DEXTER

thumbnail Fig. 5

Retention of solids in our fiducial model and in a model without backreaction of solids on the gas flow. Initially all the solids are in the form of dust (dotted line at 1). Solids are lost primarily because of the radial drift. If the radial drift is independent of the dust-to-gas ratio (“without backreaction”, black dashed line), 99% of solids are lost within 2 Myr and no planetesimals are formed. In the model presented in this paper that includes the backreaction, the radial drift is slowed down when dust piles up. This allows for planetesimal formation (red solid line) which keeps around 17% of the initial mass of solids (light grey solid line) after the rest of the disc disperses (dark grey solid line).

Open with DEXTER

In its outer, drift-dominated part, the dust density gradually decreases, but the interplay between dust growth and drift leads to the formation of a temporary pile-up in the inner disc (upper panel of Fig. 4). The formation of this pile-up is aided by the disc dispersal. As can be seen in the bottom panel of Fig. 2, the maximum drift speed decreases with time in the inner part while increasing in the outer. This means that less and less solids flow out through the inner edge of the disc, which we set to 0.3 AU for dust, based on the evaporation temperature of silicates, while more and more dust arrives to the inner part from the outside.

As mentioned in Sect. 2.2, the radial drift velocity in our models includes the so-called collective effects or, in other words, the effects of backreaction of dust on gas, that becomes important when dust-to-gas ratio increases. The drift slows down when the dust-to-gas ratio increases. Because of this, the initial pile-up becomes stronger, finally leading to conditions that are sufficient for triggering planetesimal formation via the streaming instability (see Sect. 2.3) after ~3 × 10-5 yrs (lower panel of Fig. 4).

As the initial disc profile is shallow, corresponding to roughly Σgr-1, at the beginning of the simulation most of the mass of solids is located in the outer part of the disc. As the dust growth timescale is much longer at large orbital radii, these solids are then delivered to the inner disc later, which feeds the zone where planetesimal formation is possible and extends the planetesimal formation period until about 1 Myr of evolution. Afterwards, the inner disc becomes too depleted and planetesimal formation stalls.

The gradual loss of solids, as well as planetesimal formation, are depicted in Fig. 5. In our fiducial model, the total mass of solids (light grey solid line), split between dust (dark grey solid line) and planetesimals (red solid line), never actually decreases to zero and this is because about 17% of the initial dust mass, which corresponds to 60 Earth masses, is turned into planetesimals. The planetesimals do not undergo radial drift. For comparison, we show the evolution of solids in the case when the effects of backreaction on the radial drift speed are not included (black dashed line). No planetesimals are formed and all the solids are lost in this case.

The amount of planetesimals formed, as well as the location of the annulus in which they form, changes with the initial conditions. In the following sections, we describe the dependence of these results on some of the parameters: disc mass Mdisc, initial dust-to-gas ratio Z, fragmentation threshold velocity vth, and planetesimal formation efficiency ζ. All these results are summarised in Fig. 6 and Table 1.

3.2. Dependence on disc mass

Changing the mass of the disc alters the growth barriers (Fig. 1) and thus the dust evolution. In a lower mass disc, the relative importance of the radial drift regime increases, moving the inner pile-up region even closer to the star (panel a in Fig. 6). Fortunately, at the same time the lower mass discs are colder and, because to this, the maximum size of grains in the fragmentation-limited case increases (see Eq. (11)). This makes it possible to grow aggregates to the pebble size and still form planetesimals, although much closer to the star. Surprisingly, we find that the mass of planetesimal annulus relative to the initial mass of solids in the disc does increase with decreasing disc mass (6th column of Table 1). This is because the lower mass discs disperse faster, facilitating the inner pile-up. However, the absolute mass of planetesimals decreases and is as low as six Earth masses for the disc with initial mass of 0.005 M. For the high mass discs, planetesimal formation becomes harder as the fragmentation barrier does not allow the formation of sufficiently large pebbles. No planetesimals are formed in discs more massive than 0.2M.

3.3. Dependence on metallicity

Higher initial dust-to-gas ratios makes planetesimal formation much easier and more efficient. This is not only because there are more solids available, but also because the radial drift regime is pushed outwards as the growth timescale is shorter for higher dust densities (see Eq. (3)). We find that the solar metallicity of Z = 10-2 is the lowest that allows planetesimal formation (while keeping all the other parameters at their default values). The total mass of planetesimals that form increases linearly with the initial mass of solids, reaching 1400 Earth masses for a disc that initially contains 5% solids. At the same time, the outer edge of the planetesimal annulus moves outward, so the planetesimal formation region is broader, reaching 3 AU for Z = 0.05. These results should naturally lead to the well-known metallicity-giant planet occurrence correlation (Udry & Santos 2007).

3.4. Dependence on fragmentation threshold

We find that the aggregates should be able to grow at impact velocities of up to 8 m s-1 to allow planetesimal formation. This is because sufficiently large pebbles need to be produced and the maximum size of grains depends strongly on the fragmentation threshold (see Eq. (11)). On the other hand, we observe that planetesimal formation does not occur if the aggregates are allowed to grow at velocities above vth = 15 m s-1. This is because the aggregates with Stokes numbers close to unity, which are formed with such a high fragmentation threshold, drift very quickly and are not able to form a significant pile-up.

The properties of the planetesimal annulus change with vth. For higher velocities, the annulus moves inwards because of the faster radial drift. The mass of planetesimals that form slightly increases with higher vth values, reaching 69 Earth masses for vth = 11 m s-1, and decreases to 12 Earth masses for vth = 15 m s-1.

Table 1

Summary of parameters used and results obtained in different models.

3.5. Dependence on planetesimal formation efficiency

Equation (16) connects planetesimal formation to evolution of the surface density of pebbles. The parameter ζ determines what fraction of pebbles is turned to planetesimals within one orbital period.

Recent local hydrodynamical simulations of streaming instability presented by Simon et al. (2016) report that planetesimal formation saturates on a timescale of a few tens orbital periods at a level of ~50% of pebbles being turned to planetesimals. This would correspond to ζ ≈ 10-2. However, these models used pebbles corresponding to St = 0.3, order of magnitude larger than considered in this paper. For smaller pebbles, the timescale of streaming instability increases, which forces us to also consider lower ζ values.

We tested a wide range of values for this planetesimal formation efficiency parameter. Naturally, higher ζ values lead to more planetesimals being formed. The width of the planetesimal annulus slightly changes, reaching a minimum for the highest value of ζ = 10-2. Setting ζ to values even higher than 10-2 leads to numerical instabilities since too great a mass of dust is removed in one timestep. The mass of planetesimals that form naturally increases with increasing ζ and saturates at 82 Earth masses for ζ = 10-3. For lower ζ values, the radial drift competes with planetesimal formation in taking the pebbles away from the pile-up region.

thumbnail Fig. 6

Surface density of planetesimals formed in models with different initial a) disc mass; b) dust-to-gas ratio; c) fragmentation threshold; d) planetesimal formation efficiency. All the panels share a common colour map for convenience.

Open with DEXTER

4. Discussion

Our results show that global redistribution of solids driven by dust growth and radial drift leads to a local pile-up of pebbles in the inner, fragmentation dominated part of the protoplanetary disc. This pile-up generates conditions necessary for planetesimal formation. Because the planetesimal formation may only be triggered under such specific conditions, planetesimals do not form at every orbital distance. This is qualitatively consistent with the latest observations of protoplanetary discs, which suggest that axisymmetric pile-ups of pebbles occur at various distances from the central star (ALMA Partnership et al. 2015; Andrews et al. 2016). In the models presented in this paper, we do not include processes, such as zonal flows that may potentially create pressure bumps in the outer parts of the disc (Flock et al. 2015). Without these bumps, the pebbles pile-up in the inner part of the disc, where their subsequent evolution is driven by fragmentation rather than radial drift.

Planetesimal formation is not very efficient in our models, typically incorporating about 20% of the total solids. However, as this mass is packed into a relatively narrow annulus, the surface density of planetesimals is fairly high, exceeding that predicted by MMSN. The total mass of planetesimals ranges from 1 Earth mass to above 1000 Earth masses in different models. However, we do not know how much mass could actually end up in the final planets – this depends on the efficiency of the late stages of planetary accretion.

Implications of this local, close-in planetesimal formation for the architectures of planetary systems will be the subject of our future work. Our scenario has a promising potential for explaining some of the solar system features. Raymond et al. (2009) show that it is impossible to reproduce the architecture of the inner solar system starting from an initially uniform distribution of planetesimals. Mars analogues produced in these models were always significantly more massive than the actual Mars. The low mass ratio between Mars and Earth suggests that a significant depletion of solids outside of 1 AU existed before the final assembly of the terrestrial planets (Hansen 2009; Izidoro et al. 2014, 2015). One possible explanation for this depletion is that it was sculpted by the migration of Jupiter and Saturn, which is known as the Grand Tack scenario (Walsh et al. 2011). Our results produce initial conditions necessary to reproduce the masses of terrestrial planets in a much more straightforward way.

Another solution to the low mass of the Mars problem is offered by the so-called pebble accretion scenario. Pebbles corresponding to Stokes number of 10-2−1 may be accreted very efficiently onto planetary embryos (Ormel & Klahr 2010; Lambrechts & Johansen 2012). However, the growth efficiency depends on the embryo size and its location in the disc such that there is a steep cutoff between embryos that can and cannot grow (Levison et al. 2015). The rate of the pebble accretion depends also on the flux and size of accreting pebbles. Large pebbles would lead to a very efficient growth of embryos and formation of too many planets in the terrestrial planet region (Kretke & Levison 2014). Recently, Johansen et al. (2015) present more realistic models with pebbles corresponding to Stokes numbers more similar to those suggested by our models. In this case, the pebble accretion contributes to the growth of terrestrial planets at the same rate as planetesimal accretion.

Our models indicate the formation of a significant pile-up of pebbles close to the central star that then allows for planetesimal formation via the streaming instability. These pile-ups in the inner disc were previously reported in the dust evolution models of Youdin & Shu (2002), Youdin & Chiang (2004), and Birnstiel et al. (2012). A global transition from shallow to steep surface density profiles of pebbles was also observed recently in the TW Hydrae disc (Hogerheijde et al. 2016). The universality of pebble pile-ups would support the in situ formation scenario for the tightly packed exoplanetary systems detected by the Kepler mission (Fang & Margot 2012).

In our work, we use a relatively simple 1D model of a viscous protoplanetary disc. However, many of the observed discs do not appear smooth or axisymmetric (Mayama et al. 2012; Rameau et al. 2012; Benisty et al. 2015), and large fraction of discs exhibit an extended inner cavity (Andrews et al. 2011). The peculiar shapes of observed protoplanetary and transition discs may be explained by planet-disc interactions (e.g. Pohl et al. 2015; Pinilla et al. 2016). It is also known that presence of planet-induced gaps in discs facilitates the formation of further planets (Kobayashi et al. 2012; Pinilla et al. 2015). Thus, it is formation of the first planet that is most problematic to explain and this is the problem we have considered, and using a simple disc model is justified by the uncertainty of the disc morphology during the very early stages of evolution.

The shallow gas disc profile we employ is consistent with observations of circumstellar discs (Andrews & Williams 2007). Birnstiel et al. (2012) notice that the steeper surface density profile Σdr− 3/2, consistent with the MMSN, is naturally obtained by the dust component during its evolution in the fragmentation dominated regime, even in a primarily shallow disc, as long as the turbulence level is low. This means that the initially constant dust-to-gas ratio changes and it can be enhanced in the inner regions of a shallow disc. In a steep disc, such as the MMSN with Σgr− 3/2, this type of pile-up does not occur and the inner part of the disc is gradually depleted with the dust-to-gas ratio remaining independent of radial distance. As shown in Fig. 3, the disc becomes more and more drift-dominated during its evolution. In the drift-dominated regime the surface density of dust evolves to Σdr− 3/4, which was also found by Lambrechts & Johansen (2014).

As mentioned in Sect. 2.1, we do not include disc heating by stellar irradiation. The midplane temperature drops to ~10 K at large distances from the star, close to the assumed temperature of interstellar medium. By way of contrast, models including irradiation from the star predict higher temperatures (see, e.g. Bitsch et al. 2015), but these models assume that the accretion rate in the disc is uniform, which is correct close to the star, but questionable at larger distances. In our model, we do not make any assumptions on the time evolution of the accretion rate, but we rather assumed some value for the viscosity (or more precisely the αt parameter) and a value of the mass loss associated with photoevaporation (see Sect. 2.1). We leave testing another disc models for future work.

A caveat of our models is that we do not explicitly account for the impact of dust on the global evolution of the gas disc, which may be important when the local dust-to-gas ratio reaches unity, since the gas rotation would be modified in this case (Taki et al. 2016). However, since we are mostly interested in obtaining the conditions necessary for triggering the planetesimal formation, before the dust-to-gas ratio reaches such a high value, this is a justified approach.

In this paper, we focus on compact growth during which the internal density of aggregates remains constant. Porous growth may significantly influence the fragmentation and radial drift barriers (Okuzumi et al. 2012; Kataoka et al. 2013) is well known. Recently, Krijt et al. (2016) investigated the global evolution of porous aggregates using another simplified method. They find that direct planetesimal growth is hindered by erosion but the sizes of dust aggregates might be sufficient for triggering the streaming instability. However, they did not find significant pile-ups of material as we do, but this is because they only adopted the MMSN model in which the pile-up does not occur as we explained earlier in this section.

We found that fragmentation threshold has to be above 8 m s-1 to allow us for the formation of sufficiently large pebbles in our fiducial disc. These high fragmentation thresholds can be obtained for porous silicate aggregates (Meru et al. 2013), preferably built from small monomers as the fragmentation velocity strongly depends on the monomer size (Wada et al. 2013). Thus, the bouncing that is found in laboratory experiments that consider the silicate grains (Güttler et al. 2010; Kelling et al. 2014), has to be excluded to allow planetesimal formation. This means that planetesimals can only form if grains are relatively sticky. However our models also suggest that the grains cannot grow too large to allow for the pile-up formation (see Sect. 3.4). The so-called stickiness of dust aggregates is strongly connected to their constituent material: icy grains are considered to be significantly more sticky than silicates (Aumatell & Wurm 2014) and thus the fragmentation threshold velocity should depend on the local temperature. We do not include these effects in this study. However, recent models show that if the threshold velocity at which the growth stalls decreases inward to the snow line, a pile-up of solids arises around this location (Banzatti et al. 2015; Estrada et al. 2016), which might be another mechanism for triggering the planetesimal formation. We leave the investigation of this effect for future studies.

5. Conclusions

This paper addresses the connection between dust evolution and planetesimal formation, which represents a major gap in state-of-the-art planet formation models. As dust growth is limited by fragmentation and radial drift, a direct growth to planetesimal sizes seems unlikely. However, the same processes drive a global redistribution of solids and may lead to a pile-up of pebbles that triggers planetesimal formation via the streaming instability.

We show that a narrow ring of planetesimals with a steep surface density profile is naturally produced in the inner part of a shallow protoplanetary disc corresponding to the observed ones. These planetesimals form from the flux of solids originating from the outer disc that are carried by the radial drift. At the same time, the radial drift limits the radial extent of this annulus, as it prevents the dust-to-gas ratio enhancement from spreading further out. On the other hand, the inner edge is caused by lack of sufficiently large pebbles that could trigger the streaming instability. Impact velocities increase towards the inner disc because of higher temperatures and thus the maximum size of grains decreases. The exact properties of planetesimal annulus depend on various parameters, as shown in Fig. 6 and summarised in Table 1.

The narrow annulus of planetesimals around 1 AU enables us to reproduce the unusual masses of terrestrial planets, where Earth and Venus are significantly more massive than Mercury and Mars (Hansen 2009). In some of our models, the planetesimal annulus contains as much as 1000 Earth masses, which could allow for in situ formation of the close-in massive planets that have been detected around many stars (Fressin et al. 2013).

Our results, including the surface density of planetesimals and timing of their formation, may be used as an input to models which investigate the later stages of planet accretion, planet population synthesis, and the internal evolution of asteroids. The size and flux of pebbles that we obtain is important for the pebble accretion models that typically assume unphysically large aggregates and make ad hoc assumptions about their surface density.

Acknowledgments

The authors thank Bertram Bitsch, Kees Dullemond, Barbara Ercolano, Anders Johansen, and Chris Ormel for useful discussions. We also thank the referee for their useful report that helped us to significantly improve this paper. This work has been carried out within the frame of the PlanetS National Center of Competence in Research (NCCR) of the Swiss National Science Foundation. The numerical models were performed on the zBox4 CPU cluster at the Institute for Computational Science, University of Zurich.

References

Appendix A: Radial drift dependence on dust-to-gas ratio and size distribution

Most of the dust evolution models use the radial drift speed derived by Weidenschilling (1977)(A.1)where ηvK is the maximum drift speed that results from the gas pressure gradient and St is the Stokes number of a particle.

Already Nakagawa et al. (1986) noticed that the radial drift velocity should be reduced when the solids density is high. They derived the dependence of the drift speed on the dust-to-gas ratio in the case of monodisperse dust population (A.2)where ϵ is the local dust-to-gas ratio (of a single sized dust species).

However, in a realistic case, when we have dust particles of different sizes, the dust species interact and exchange the angular momentum leading to drift velocity being dependent on the particular size distribution. In this paper, we implement this effect, which turns out to be of great importance to the dust retention and planetesimal formation (see Fig. 5). We use the formulation derived by Tanaka et al. (2005, their Eqs. (14)(18)), which can also be found in Okuzumi et al. (2012; their Eqs. (48)(53)). The drift velocity of particles with Stokes number St is given by (A.3)where the radial velocity of gas is (A.4)and the difference between the azimuthal velocity of gas and the Keplerian velocity is (A.5)The X and Y are dimensionless quantities dependent on the mass distribution where ϵ(m) = ρd(m) /ρg.

We note that the radial drift velocity described by Eq. (A.3) naturally converges to the nominal radial drift solution of Weidenschilling (1977) (Eq. (A.1)) for low dust-to-gas mass ratios (0.01, the black line in Fig. A.1) and of Nakagawa et al. (1986) (Eq. (A.2)) for single sized grains.

thumbnail Fig. A.1

Midplane radial drift speed dependence on the Stokes number for a simple mass distribution and vertically integrated dust-to-gas ratios of 0.01, 0.03, and 0.1. The inward drift is significantly reduced for the higher dust-to-gas ratio. Outward drift is possible for the smallest grains.

Open with DEXTER

The same process was also described by Bai & Stone (2010a, see their Appendix A), who derived a solution based on the matrix formalism (their Eqs. (A.1)(A.5)). Figure A.1 shows that these two approaches are equivalent and demonstrate how this effect works for vertically integrated dust-to-gas ratios of 0.01, 0.03, and 0.1 (1, 3, and 10 times the solar value). For this plot, we used a power-law size (and Stokes number) distribution that is an effect of the coagulation-fragmentation equilibrium (Birnstiel et al. 2011) (A.8)for seven equally spaced logarithmic bins between St = 10-4 and St = 10-1. We assumed that there is an equilibrium between settling and turbulent mixing with αt = 10-3, resulting in the midplane dust-to-gas ratios described by Eq. (14). The inward drift of the largest particles is suppressed for the increasing dust-to-gas ratio and even the outward drift (vr> 0) is possible for the smallest grains. The average drift velocity depends both on the total dust-to-gas ratio and mass distribution.

We assume the size distribution described with Eq. (A.8) in the fragmentation regime case in our models. The mass weighted average that is used in Eq. (9) is calculated as (A.9)We find that the number of mass bins plays a role in the radial drift velocity calculation. We use 200 bins, since we found the convergence of for this value.

Appendix B: List of symbols used in the paper

In this Appendix, we compile Table B.1 which describes symbols used throughout this paper and provides their typical value when applicable.

Table B.1

List of symbols used in this paper.

All Tables

Table 1

Summary of parameters used and results obtained in different models.

Table B.1

List of symbols used in this paper.

All Figures

thumbnail Fig. 1

Overview of the growth barriers in our fiducial disc model. The red line corresponds to a Stokes number of unity. The yellow region marks the fragmentation barrier where the impact speeds would be higher then the threshold value of 10 m s-1. The grey region shows the size scale for which the drift timescale is shorter than the growth timescale, meaning that the aggregates would be removed from a given location faster than the growth could replenish them: this is the drift barrier. The light grey vertical lines are tracks of test particles evolving by growth and drift in a steady-state disc. Finally, the blue nearly horizontal line shows the size scale corresponding to the Stokes number of 10-2, which is a minimum size required for planetesimal formation via the streaming instability.

Open with DEXTER
In the text
thumbnail Fig. 2

Time evolution of the gas surface density Σg, gas temperature, and the headwind speed ηvK in our fiducial disc model. The disc mass is initially 0.1M and it decreases to 0.02M after 10 Myrs. The disc is shallower then the Minimum Mass Solar Nebula (MMSN, red line) and thus most of its mass resides in its outer parts. The headwind, which determines the maximum radial drift velocity, is not constant throughout the disc, as opposed to the MMSN model, and it gradually decreases during the evolution, facilitating the retention of solids.

Open with DEXTER
In the text
thumbnail Fig. 3

Time evolution of the growth barriers. The plotting style corresponds to Fig. 1: the yellow region corresponds to the fragmentation barrier while the grey region indicates the drift regime. The red and blue horizontal lines denote sizes corresponding to the Stokes number of unity and 10-2, respectively. The purple line shows the maximum size of dust aggregates, which generally follows the growth barriers, except for the earliest stages of evolution, when it is restricted by the growth timescale. The radial drift barrier gains significance as the dust-to-gas ratio decreases during the evolution.

Open with DEXTER
In the text
thumbnail Fig. 4

Time evolution of a) the vertically integrated dust-to-gas ratio and b) surface density of planetesimals. The outer parts of the disc are depleted of dust because they are drift dominated. At the same time, there is a temporary pile-up of dust in the fragmentation-dominated inner disc, which enables planetesimal formation via the streaming instability. The dashed line in panel b) corresponds to the surface density of solids in the MMSN model.

Open with DEXTER
In the text
thumbnail Fig. 5

Retention of solids in our fiducial model and in a model without backreaction of solids on the gas flow. Initially all the solids are in the form of dust (dotted line at 1). Solids are lost primarily because of the radial drift. If the radial drift is independent of the dust-to-gas ratio (“without backreaction”, black dashed line), 99% of solids are lost within 2 Myr and no planetesimals are formed. In the model presented in this paper that includes the backreaction, the radial drift is slowed down when dust piles up. This allows for planetesimal formation (red solid line) which keeps around 17% of the initial mass of solids (light grey solid line) after the rest of the disc disperses (dark grey solid line).

Open with DEXTER
In the text
thumbnail Fig. 6

Surface density of planetesimals formed in models with different initial a) disc mass; b) dust-to-gas ratio; c) fragmentation threshold; d) planetesimal formation efficiency. All the panels share a common colour map for convenience.

Open with DEXTER
In the text
thumbnail Fig. A.1

Midplane radial drift speed dependence on the Stokes number for a simple mass distribution and vertically integrated dust-to-gas ratios of 0.01, 0.03, and 0.1. The inward drift is significantly reduced for the higher dust-to-gas ratio. Outward drift is possible for the smallest grains.

Open with DEXTER
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.