Free Access
Issue
A&A
Volume 620, December 2018
Article Number A134
Number of page(s) 18
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201834047
Published online 10 December 2018

© ESO 2018

1 Introduction

Rocky planets and the cores of gas giants form from micrometer-sized dust grains in gaseous disks around young stars. It is generally accepted that an intermediate stage on the way from small dust grains in protoplanetary disks to full-sized planets is the formation of approximately kilometer-sized planetesimals, which mark the transition to a gravitation-dominated growth phase (Safronov 1969; Pollack et al. 1996; Benz 2000). Planetesimal formation is still an active field of research since theories face several problems. First, typical relative velocities between approximately centimeter-sized particles are often too large for coagulation, such that particles fragment or bounce off each other upon collision, rather than stick to form even larger particles (Blum & Wurm 2000; Zsom et al. 2010). Second, particles that are large enough to be aerodynamically decoupled from the gas disk (“pebbles”) lose angular momentum and drift toward the central star (Whipple 1972; Weidenschilling 1977): planetesimals need to form before the solids are lost due to this process.

A promising method to form planetesimals directly from small particles – without the need for growth through all intermediate sizes – is by streaming instability: drifting pebbles clump together and can collapse under their own gravity to form planetesimals (Youdin & Goodman 2005; Johansen et al. 2007; Johansen & Youdin 2007). For streaming instability to operate, however, a solids-to-gas ratio that is enhanced with respect to the canonical value of 1% is required (Johansen et al. 2009; Bai & Stone 2010; Carrera et al. 2015; Yang et al. 2017). In classical planet formation theory, planetesimals are assumed to form throughout the entire disk, but various studies propose that an enhanced solids-to-gas ratio forms preferentially at specific locations instead, for example in the inner disk (Dra̧żkowska et al. 2016), in the outer disk at late times (Carrera et al. 2017), or near the water snowline, where water transitions from the vapour phase to the solid phase (Saito & Sirono 2011; Ros & Johansen 2013; Ida & Guillot 2016; Armitage et al. 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). Detections by the Atacama Large Millimeter/submillimeter Array (ALMA) of substructures in protoplanetary disks such as axisymmetric gaps and rings (ALMA Partnership 2015; Akiyama et al. 2015; Nomura et al. 2016; Isella et al. 2016; Fedele et al. 2018) may support the view that planetesimals form at specific distances from the star. Ormel et al. (2017) present a complete formation scenario for the TRAPPIST-1 planets that is based on the formation of planetesimals at a single location: the water snowline. However, this scenario still lacks a thorough numerical model.

In this work we present a new, versatile model for planetesimal formation. We follow the growth and drift of pebbles in protoplanetary disks and include planetesimal formation via the streaming instability, partly based on the results of Schoonenberg & Ormel (2017; hereafter: SO17). In SO17 we presented a local model where we investigated whether water diffusion and condensation could lead to conditions conducive to streaming instability outside the water snowline, and how this depends on certain parameters such as the turbulence strength and the particle size. In contrast to SO17, in the current work we consider the entire protoplanetary disk. The pebble mass flux and particle sizes are not treated as input parameters as in SO17, but follow self-consistently from the simulation and evolve in time. Another difference with SO17 is that in this paper, the planetesimal formation process is followed (mass is removed from pebbles as planetesimals form), whereas in SO17 we focused only on the conditions for planetesimal formation. Our model makes use of the Lagrangian method where super-particles represent groups of particles with identical physical properties. In the context of dust evolution in protoplanetary disks, this method waspioneered by several studies (e.g., Laibe et al. 2008; Krijt et al. 2016; Gonzalez et al. 2017). In contrast to Eulerian approaches to dust evolution (e.g., Birnstiel et al. 2012; Okuzumi et al. 2012), this method is mesh-free: physical quantities are computed at the locations of the super-particles, not at the locations of grid cells, and individual particle characteristics such as water content can easily be followed as particles move through the disk.

Exoplanet data show that low-mass stars are efficient at forming super-Earths (Mulders et al. 2015; Mulders 2018). Differences in the occurrence rates between stars of different masses could (partly) originate in differences in planetesimal formation efficiencies. In this paper we apply our model to testing the efficiency of icy and dry planetesimal formation as a function of stellar mass. We find that generally, low-mass stars convert a larger fraction of the solids reservoir in their disks to planetesimals than high-mass stars, but other disk parameters play important roles as well. We also discuss the effects that pebble accretion, the viscous evolution of the gas disk, changing the metallicity and fragmentation threshold have on planetesimal formation.

We describe the different components of our model in Sects. 24 and provide a summary of the model in Sect. 5. The results are discussed in Sects. 67 and we discuss our findings in Sect. 8. Our main conclusions are listed in Sect. 9.

2 Gas disk model

2.1 Surface density profile

In our standard disk model, we adopt for simplicity a steady-state gas surface density profile Σgas for a given gas accretion rate gas and dimensionless turbulence parameter α (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974): (1)

where the viscosity ν is related to α as follows: (2)

with cs the sound speed and Ω the Keplerian orbital frequency. The gas moves inward at a speed |vgas| = 3ν∕2r where r is the radial distance from the star. The surface density Σgas for our fiducial model (Table 1) is plotted by the black line in Fig. 1.

We will discuss the validity of the assumption of constant gas in the context of our model, and investigate the effects of relaxing this assumption in Sect. 7.2.

Table 1

Input parameters, their symbols, and their fiducial values.

thumbnail Fig. 1

Gas disk profiles for our fiducial model (Table 1). The black line corresponds to the gas surface density profile Σgas (Eq. (1); the solid blue line corresponds to the temperature profile T, which receives contributions from viscous heating (dashed blue line) and stellar irradiation (dotted blue line; Eq. (5)).

Open with DEXTER

2.2 Temperature profile

We consider two mechanisms that heat the protoplanetary disk: viscous heating and stellar irradiation. Viscous heating is only important in the innermost region of the disk and leads to a radial dependence of Tviscr−3∕4 (Frank et al. 2002). Stellar irradiation results in a temperature profile in the outer disk that goes as Tirrr−1∕2 (Kenyon & Hartmann 1987). These general temperature profiles can only be directly used for the disk midplane temperature (of interest in this work) if the vertical optical depth were radially constant, which is probably not the case. A more sophisticated model would calculate the temperature from the dust properties self-consistently. However, for the purposes of the current work, the simple power laws suffice. For our fiducial model, we use the following viscous and irradiated midplane temperature profiles:

We take the global temperature profile T(r) to be the (smoothed) maximum of the temperature profiles Tvisc and Tirr: (5)

which then defines the isothermal sound speed profile cs: (6)

where kB is the Boltzmann constant and μ is the mean molecular weight of the gas, for which we take a value of 2.34 times the proton mass appropriate for a solar metallicity gas. The temperatureprofiles T, Tvisc, and Tirr for our fiducial model (Table 1) are plotted by the blue lines in Fig. 1.

We assume that the disk is vertically isothermal, such that the disk scale height Hgas is given by (7)

3 Treatment of solids: Lagrangian smooth-particle model

In this work we adopt a Lagrangian method to solve for the growth and the radial movement of the solid particles in the disk. A clear advantage of the Lagrangian method is that it is naturally suited to follow characteristics of individual particles, such as their composition and porosity, as they grow and move through the disk. It is not possible to follow all particles in the disk, however. We therefore use super-particles: groups of particles with the same physical properties. In this work we consider two different classes of super-particles: dust/pebble super-particles1, and planetesimal super-particles (in future work, more super-particle classes may be added to the model). Pebble super-particles can move through the disk, whereas planetesimal super-particles are inert. We assume a mono-disperse particle size distribution at each point in space and time, such that a pebble super-particle represents particles that all have the same individual particle mass.

The characteristics of the pebble super-particles that we follow are location, total mass, individual particle mass, water mass fraction.

The planetesimal super-particles in our model are characterized by location (fixed), total mass, and water mass fraction.

The total mass of a planetesimal superparticle can change because of streaming instability and pebble accretion: when planetesimal formation and/or pebble accretion occurs, mass is transferred from the pebble super-particles to the planetesimal super-particles.

Formally, we end up with a system of ordinary differential equations (ODEs) that govern the time evolution of the super-particles, of the following form: (8)

where i indicates the super-particle number and j corresponds to their properties: radial position r, individual particle mass mp (for pebble super-particles only), total super-particle mass M. The right-hand side of Eq. (8) is a vector Y that consists of the planetesimal formation rates, drift velocities, and particle growth rates, which depend on the properties of the super-particles X. We solve Eq. (8) making use of the Python package scipy.integrate.ode.

3.1 Smooth-particle method

In Krijt et al. (2016), the solids surface density and its radial derivative at the location of each super-particle is calculated using a “tripod” method: each super-particle consists of three “legs” that can move closer to each other (higher surface density) or further apart (lower surface density). A consequence of this approach is that the super-particles evolve independentof each other. In contrast, in our method all super-particles are connected. As in SPH methods (e.g. Laibe et al. 2008), we approximate the solids surface density at each super-particle location using a weighting kernel W, that accounts for the contribution of neighbouring super-particles as function of their mass and distance. Therefore, we simultaneously solve for the evolution of all solids in the disk. As in Krijt et al. (2016) but unlike Laibe et al. (2008), our model is one-dimensional; we only deal with the radial dimension r. The vertical dimension of the disk is taken into account by means of the gas and solids scale heights, and the model is symmetric in the azimuthal dimension.

The value of a quantity F at each particle location x is kernel-approximated by (9)

where x and x′ are vectors and W(xx′, h) is the kernel, for which we take (Hicks & Liebrock 2000; Liu & Liu 2003): (10)

where h is the smoothing length, Δ x = |xx′| is the absolute distance between the location of super-particle x′ and the location of the super-particle-of-interest x. The prefactor ensures that the kernel is normalized. Particles that are separated from x by a distance larger than the smoothing length (|xx′| > h) do not influence the value F(x); in other words, the kernel W is compact.

Since the values of F(x) are only defined at discrete locations (the super-particle locations), we turn the integral in Eq. (9) into a sum over all simulated super-particles. For the surface density Σ at the location xi of super-particle i, we then find (11)

where Mj is the mass of supporting super-particle j and hi is the smoothing length of super-particle i. We treat the smoothing length hi as a variable, demanding that at each super-particle location and time hi takes on a value such that there are five neighbouring super-particles (including super-particle i itself) in the support group contributing to the density at location xi (for the simulations presented in this paper, changing the number of neighbours to three or seven did not change the results). The advantage of a variable smoothing length is that the code can adapt to regions of high density as well as regions of low density. Because of the compactness of the kernel W, for each super-particle we only have to sum over the super-particles in its support group.

3.2 Boundary treatment and initial conditions

The condition we set on the inner boundary of the simulated disk is a constant solids surface density gradient. The innermost particles – which lack particles interior to them to fill their support group – get assigned a surface density value in agreement with the surface density gradient in the inner region.

At the outer boundary, we use an exponentially cut-off initial solid surface density profile, such that particles close to the outer boundary barely grow and drift over the lifetime of the disk, to prevent any unwanted outer boundary effect.

At the start of the disk evolution, a fraction of the total disk mass is in the form of dust grains (ice + silicates). We take this fraction, also called the metallicity Z, to be initially constant throughout the disk with a sharp exponential cutoff at the outer disk radius rout. Therefore, the solids surface density profile Σsolids initially follows: (12)

where the sharp exponential cut-off is just to ensure that no unwanted numerical effects occur at the outer boundary. The initial size of solid particles is set to 0.1 μm. We assume that in the outer disk, the water ice mass fraction fice of particles is 50% (fice,out = 0.5) (Lodders 2003; Morbidelli et al. 2015). Interior to the water snowline, the solids consist of pure silicate (fice,in = 0). The location of the water snowline and the transition between fice,out and fice,in will be discussed in Sect. 3.5.

The initial mass of a super-particle follows from the initial solids surface density profile Σsolids(r) and the initial number of super-particles N. We discretize the radial distance to the central star r into N annuli with edges ri between the inner edge of the disk rin and the outer edge rout. The kth super-particle then gets assigned a mass (corresponding to the total mass in the kth annulus) and an initial position xk = (rk+1rk)∕2. If the initial solids surface density profile is proportional to r−1solids(r) = Σ0r−1 where Σ0 is constant) and the spacing between the annuli is linear, each super-particle has the same initial mass m = 2πΔ0 with Δ r the annulus width. Alternatively, the initial particle locations can be chosen such that the resolution varies across the disk (for example, closely spaced particles in the inner disk and particles that are further apart in the outer disk); in that case the initial super-particle masses differ. In this work we have used a logarithmic spacing between the initial positions of the dust/pebble super-particles, so that super-particles in the outer disk have larger total masses than super-particles in the inner disk.

3.3 Resampling

During the evolution of the disk, the spacing between adjacent super-particles can become larger than desired. This can, for example, occur when the resolution needs to be high in the inner disk or around a special location (e.g., the snowline), whereas in the outer disk the particles are initially further apart. Due to radial drift, closely spaced particles from the inner disk get accreted to the star and less closely spaced particles from the outer disk enter the inner disk, such that the resolution in the inner disk decreases. Ifat any point the resolution becomes too low, we initiate a “resampling” algorithm. In our model we resample when the separation between particles at a certain point becomes 20% larger than the initial particle separation at that point. The number andlocations of super-particles are reset to the initial configuration, and the characteristics of the new super-particle population are sampled from the super-particle population right before the resampling process. The total mass of each new super-particle is extracted from the cumulative mass distribution right before “resampling”, in order to ensure mass conservation. Figure 2 shows a cartoon of the concept of resampling.

thumbnail Fig. 2

Procedure cartoon of the resampling concept. When the separation between neighbouring pebble super-particles becomes toolarge, the number and configuration of particles is reset such that the resolution is again equal to that of the initial set-up. The total mass in pebble super-particles is conserved and the characteristics (such as water ice fraction, indicated bycolor here) are sampled from the situation before resampling.

Open with DEXTER

3.4 Particle growth and radial drift

In this section we describe how we treat dust evolution (particle growth and radial drift) in our model. In Appendix A we test our model against results from the literature.

The dust population at the locations of the super-particles is represented by a monodisperse dust size distribution (Krijt et al. 2016; Sato et al. 2016). This means that the dust particles represented by one super-particle all have the same size sp and mass mp. Within each super-particle then, particles can coagulate, increasing the individual particle mass mp of that super-particle. The mass growth rate dmp∕dt assuming perfect sticking is (Krijt et al. 2016) (13)

with Hsolids the solids scale height, σcol the collisional cross section (which we take equal to the geometrical cross section), and vrel the relative velocity between particles. As long as the particles are small and well-coupled to the gas, Hsolids = Hgas. Once particles have grown to sizes where they start to aerodynamically decouple from the gas disk, their scale height is reduced with respect to that of the gas due to vertical settling, and is given by (Youdin & Lithwick 2007) (14)

where τ is the dimensionless stopping time, which measures the degree of coupling between the particle and the gas. We distinguish between two regimes of τ: in the Stokes regime, the particle size is larger than the mean-free path of the gas molecules lmfp, and the stopping time is calculated in a fluid description; in the Epstein regime, the particle size is smaller than the mean-free path of the gas molecules, and a particle description is required instead (Weidenschilling 1977). The mean-free path lmfp is given by (15)

where ρgas is the gas volume density at the disk midplane () and σmol is the molecular collisional cross-section, for which we take σmol = 2 × 10−15 cm2 as appropriate for hydrogen (Chapman & Cowling 1970). The dimensionless stopping time is given by (16)

where sp is the particle radius, ρ∙,p is the particle internal density, and vth is the thermal velocity of the gas molecules, defined as . The particle radius and internal density are determined from the particle mass mp and water mass fraction fice:

where ρ∙,ice = 1 g cm−3 is the density of a pure ice particle and ρ∙,sil = 3 g cm−3 the density of a pure silicate particle.

For therelative velocity between individual particles within a super-particle vrel, we consider the relative velocity due to turbulence, Brownian motion, and radial and azimuthal drift, in the same way as Krijt et al. (2016). The total relative velocity is given by the maximum of these four contributions. We find that in most cases, the turbulent relative velocity vrel,turb dominates. In the so-called intermediate regime – where the turnover time scale of the smallest eddies ( where ReT is the turbulence Reynolds number) is much smaller than the stopping time of the particles (Ormel & Cuzzi 2007), which is the regime of importance in our model – vrel,turb is given by (19)

Therefore, the larger the particles grow, the more violently they impact each other.

3.4.1 Fragmentation

When the relative velocity vrel between particles becomes larger than a certain fragmentation threshold vfrag, particles start to fragment rather than coagulate when they collide. Concerning vfrag we distinguish between icy particles and dry particles. The surface energy of an aggregate containing water ice is larger than that of a silicate one, and therefore icy particles are more sticky than silicate particles (Chokshi et al. 1993). Motivated by laboratory experiments, for the fragmentation threshold velocity of silicate particles, we take vfrag,sil = 3 m s−1 (Blum & Münch 1993); and for icy particles, we take vfrag,ice = 10 m s−1 (Sirono 1999; Gundlach & Blum 2015). The fragmentation threshold vfrag for particles with ice fraction fice is then given by (20)

where we simply coupled the two fragmentation velocities according to the ice mass fraction fice (similar to Lorek et al. 2016, who interpolate the fragmentation velocities using the fractional abundance of icy monomers). Our assumption that the fragmentation velocity increases linearly with fice is arbitrary but unimportant, since fice changes rapidly across the water snowline (Dra̧żkowska & Alibert 2017).

If vrel > vfrag, we instantaneously set the particle mass mp to the fragmentation mass limit mfrag. In Fig. 3 we have plotted mfrag as a function of semi-major axis r for our fiducial model. The assumption that fragmentation is instantaneous is justified by the fact that mfrag increases for particles that are drifting from the outer parts of the disk inward within the drift-limited region. In this region, drift is faster than growth/fragmentation (Birnstiel et al. 2012). In the viscosity-dominated inner disk, mfrag decreases for inward-drifting particles (see Fig. 3). In this region therefore, inward-drifting particles that have a mass that is limited to mfrag are constantly fragmenting to smaller masses on their way to the star. However, the viscosity-dominated region is fragmentation-limited, meaning that growth/fragmentation are faster than drift (Birnstiel et al. 2012). Therefore, as long as the region where the temperature is dominated by viscous heating is in the fragmentation-limited regime, treating fragmentation as an instantaneous process is justified throughout the disk.

thumbnail Fig. 3

Fragmentation mass limit mfrag (black line) for our fiducial model with initial ice mass fraction profile fice (blue line). The snowline is located at around 2.4 au. Because of the different stickiness between icy and silicate particles, the fragmentation velocity, and therefore the fragmentation mass limit, varies rapidly across the snowline.

Open with DEXTER

3.4.2 Radial drift

The radial velocity of solid particles vp (taking into account the back-reaction of the solids onto the gas) is given by (Nakagawa et al. 1986): (21)

where vgas is negative (the gas is moving inward as well), ξ is the midplane solids-to-gasvolume density ratio, and ηvK is the difference between the azimuthal motion of the gas disk and the Keplerian velocity vK: (22)

with P the gas pressure, given by (23)

with ρgas the midplane volume density of gas.

We compare the coagulation and radial drift components of our model with well-known literature results in Appendix A.

3.5 Location of the snowline and the pebble ice fraction

As in SO17, we define the location of the snowline rsnow as the radius interior to which the water vapour pressure Pvap drops below the saturated (equilibrium) pressure Peq. Outside rsnow, the water vapour distribution always follows the saturated profile, which is given by the Clausius–Clapeyron equation: (24)

where Ta and Peq,0 are constants depending on the species. For water, Ta = 6062 K and Peq,0 = 1.14 × 1013 g cm−1 s−2 (Lichtenegger & Komle 1991). From Eq. (24) we find the saturated water vapour surface density profile Σeq: (25)

where is the molecular weight of water. However, in contrast to SO17, in this work we do not follow the water vapour distribution directly. Instead, we calculate the expected water vapour surface density profile ΣZ,a from the water mass flux that is delivered by icy pebbles : (26)

which is equal to the steady-state advection-only vapor surface density profile from SO172. We note that ΣZ,a is not the physical water vapour surface density profile, but rather the steady-state surface density profile if all the water in the disk is in the form of vapour. Under our assumptions, the physical water vapour profile would be given by the minimum of ΣZ,a and Σeq.

The location of the snowline rsnow can now be calculated by setting ΣZ,a equal to Σeq (SO17): (27)

Therefore, the location of the snowline depends on the pebble mass flux, which changes in time3. In the simulations of SO17, the water mass flux was fixed and therefore the snowline did not move. In the current work, the water mass flux follows from the simulation: we determine ΣZ,a from Eq. (26) by measuring just outside the snowline (at the innermost position where pebbles are icy) during the simulation. At the start of the simulation, the initial vapour surface density profile ΣZ,a is given by fice,outΣsolids(r, t = 0).

We define the ice fraction of pebbles fice as follows: (28)

where we have smoothed the maximum of ΣZ,a and Σeq in order to avoid a sharp transition in the ice fraction profile. Using Eq. (28) is a convenient way to calculate fice without having to model it directly by following the water vapour distribution. In Fig. 4 we show the initial Σeq and ΣZ,a profiles and the corresponding ice fraction profile fice for our fiducial model (Table 1).

thumbnail Fig. 4

Initial pebble ice fraction (blue solid line) as function of semi-major axis r for the fiducial model. The saturated water vapour profile Σeq is plotted by the black dashed line, and ΣZ,a (see text for more details) is given by the black solid line. The snowline is located where Σeq crosses ΣZ,a, in this case at rsnow ~ 2.4 au (blue dotted line).

Open with DEXTER

3.6 Diffusion of particles

Our model does not directly account for radial diffusion of solid particles. However, the semi-analytic model that we employ for planetesimal formation outside the snowline (Sect. 4.2) does include turbulent radial diffusion of water vapour and of solid particles. Also, the radial velocity of solid particles (Eq. (21)) includes the gas velocity component, meaning that if particles aresmall and well-coupled to the gas, they are moving toward the star at the gas accretion velocity. Therefore, for high values of α, small particles in the inner disk can be in the so-called mixing regime (Birnstiel et al. 2012), even if particle diffusion is not taken into account. Vertical diffusion of particles is accounted for by means of an α- and τ-dependent particle scale height (Eq. (14)).

4 Planetesimal formation by streaming instability

A promisingplanetesimal formation mechanism is streaming instability. Streaming instability occurs in the presence of radial drift, and leads to clumping of pebbles. These clumps can become dense enough to collapse under their own gravity and form planetesimals (Youdin & Goodman 2005; Johansen et al. 2007; Johansen & Youdin 2007). For streaming instability to be triggered, however, the solids-to-gas ratio needs to be locally enhanced above the typical, expected value of 1% (Johansen et al. 2009; Bai & Stone 2010; Carrera et al. 2015; Yang et al. 2017). The condition for streaming instability that we adopt in this work is a midplane solids-to-gas mass density ratio (ρsolidsρgas) exceeding unity (Dra̧żkowska & Dullemond 2014), independent of the dimensionless stopping time of particles τ. Several studies have found that the streaming instability threshold depends on τ, with increasingly higher metallicities Z ≡ Σsolids∕Σgas needed to trigger streaming instability for smaller values of τ (e.g., Carrera et al. 2015; Yang et al. 2017). However, those works considered laminar disks, whereas in our simulations the disks are intrinsically turbulent, so we cannot simply adopt their results. We find that the requirement of a midplane solids-to-gas ratio exceeding unity is always more constraining than the metallicity constraints found by Carrera et al. (2015) and Yang et al. (2017). The behavior of the streaming instability threshold on the metallicity with changing τ (increasing with decreasing τ) is also captured by our choice of threshold, because the midplane solids-to-gas ratio decreases with decreasing τ for a given metallicity: (29)

where in the last step we made use of Eq. (14). The smallest particles that participate in planetesimal formation in our simulations have τ ~ 10−3 (Sect. 7.3), which according to Yang et al. (2017) should indeed be able to trigger the streaming instability.

In this work, we are specifically (but not exclusively) interested in planetesimal formation around the water snowline. An enhanced solids-to-gas ratio can form interior to the snowline due to a traffic-jam effect caused by the variation of the fragmentation velocity across the snowline (Saito & Sirono 2011; Ida & Guillot 2016), as well as outside the snowline due to the effect of water vapour diffusion and condensation (Cuzzi & Zahnle 2004; Armitage et al. 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). In our model, we account for the effect of outward diffusion and re-condensation of water vapour by means of a recipe for planetesimal formation outside the snowline that we distill from SO17, since in this work we do not model the transport and condensation of water vapour directly.

Planetesimals are treated in a similar way as the dust particles: they are represented by super-particles containing planetesimals with the same physical characteristics. Planetesimal super-particles are characterized by a position, a total mass and a ice mass fraction. In contrast to pebble super-particles, the location of planetesimal super-particles is fixed; they are treated as sink particles. At the beginning of the simulation, a number Npltsml planetesimal super-particles are initiated at a given location, with zero mass. If the solids-to-gas ratio at a certain radius exceeds unity (either directly or inferred from the “snowline recipe”), mass is transferred from the pebble super-particles to the nearest planetesimal super-particle.

4.1 Direct streaming instability

Interior to the snowline, the solids-to-gas ratio can exceed unity “directly”, due to a traffic-jam effect caused by the change in fragmentation velocity across the snowline. If this is the case, mass is transferred from pebbles to planetesimals at a rate dMpl,SI,direct∕dt: (30)

where Mpeb,SI is the mass of the pebble super-particles that meet the requirement for SI, and tsettle = 1∕τΩ is the timescale on which SI filaments form (Yang et al. 2017), which equals the vertical settling timescale. The ϵSI parameter is an efficiency parameter, that is considered a free parameter in some studies (e.g., Dra̧żkowska & Dullemond 2014). In our model, we take ϵSI = 0.1, and we come back to the relevance of this parameter in Sect. 7.3.

4.2 Outside (but near) the snowline

A local peak in the solids-to-gas ratio that exceeds unity is formed outside the snowline if the pebble mass flux peb,outside outside the snowline (that is delivered from the outer disk) is larger than a certain critical value peb,crit, which depends on α, gas, and on the dimensionless stopping time τ of particles outside the snowline (SO17). For a given disk model with constant α and gas, we employ the semi-analytic model presented in SO17 to calculate peb,crit for different stopping times4. We assume the “many-seeds” scenario in this calculation, which leads to an ice fraction of planetesimals that is similar to that of pebbles in the outer disk: fice,pltsml = fice,out = 0.5 (Table 1). This is because not only water vapour is transported across the snowline, but small silicate grains diffuse outward as well and stick to the inward-drifting pebbles (SO17). The tabulated values of peb,crit are then used to decide whether planetesimal formation takes place outside the snowline during the simulation.

SO17 showed that the timescale tpeak,buildup on which the solids enhancement outside the snowline forms is related to the time it takes for water vapour to traverse the distance between the location where the solids-to-gas ratio peaks rpeak, and the snowline location rsnow. The width of the peak Wpeak = rpeakrsnow depends on the disk parameters and is another outcome of the semi-analytic model.

The incoming pebble mass flux outside the snowline peb,outside first increases as particles grow to pebble sizes. At some point it reaches a peak value after which it slowly decreases for the greater part of the lifetime of the disk (see also Lambrechts & Johansen 2014). As soon as peb,outside exceeds the critical mass flux, we increase the planetesimal formation rate outside the snowline from zero to dMpl,SI,snowlinedt on a peak formation timescale tpeak,buildup = Wpeak∕|vgas|. The value ofdMpl,SI,snowline∕dt is given by (31)

where we assume that the peak in the solids-to-gas ratio is continuously sustained, and only the excess incoming pebble mass flux is transformed to planetesimals5. The mass of the planetesimal super-particle closest to the peak location rpeak (which is given by the semi-analytic model) is increased at a rate dMpl,SI,snowlinedt.

We distribute the pebble mass loss rate over all pebble super-particles using a Gaussian kernel centered on the peak location rpeak. The mass loss rate of pebble super-particle i is then given by (32)

where the fraction f(ri) < 1 depends on the distance between ri and the peak location6.

Throughout the simulation, we make sure to measure peb,outside outside the region where pebbles are converted to planetesimals.

The total mass that is lost from the pebbles (disregarding the term [(1 − fice,out) + fice]; see footnote 6) is of course the exact negative of Eq. (31), such that when streaming instability outside the snowline is going on (and the solids enhancement outside the snowline has saturated), the pebble mass flux that reaches the actual snowline is peb,crit. A consequenceof this is that the location of the snowline is fixed during planetesimal formation outside the snowline (Sect. 3.5).

Once the pebble flux has dropped below the critical value peb,crit, a peak in the midplane solids-to-gas density ratio that exceeds unity cannot be sustained outside the snowline any longer, and the pebble massflux that reaches the snowline is again the same as the pebble mass flux outside the snowline, because just outside the snowlineno pebbles are converted to planetesimals anymore.

5 Model summary

To summarize the various effects that we take into account in our model, Fig. 5 shows a schematic of two possible “lifelines” (trajectories through space and time) of pebble super-particles that start in the icy outer disk and eventually convert partof their mass to planetesimals. The first super-particle (denoted by a “1” in the cartoon) begins its journey justoutside the snowline, representing small icy particles of 0.1 micron (Table 1). Once these particles have grown to pebble sizes, the super-particle starts to drift inward (Eq. (21)). At the time when it reaches the water snowline, the conditions for streaming instability outside the snowline are not yet met (the pebble mass flux does not exceed the critical mass flux yet; Eq. (31)). Therefore, the super-particle does not lose any mass to planetesimals outside the snowline. As the super-particle crosses the snowline, its H2 O ice evaporates (Eq. (28)), and the physical particles it represents fragment down to smaller sizes during the snowline crossing because of the lower fragmentation threshold for silicate particles (Eq. (20)). The inner disk is fragmentation-limited, which leads to a steeper radial dependence of the solids surface density than of the gas surface density (Birnstiel et al. 2012), and therefore to an increasing solids-to-gas ratio with decreasing distance to the star. In this cartoon the solids-to-gas ratio reaches unity close to the star, resulting in dry planetesimal formation (Eq. (30)).

The second super-particle (indicated by a “2” in the cartoon) starts out further away from the star, also representing 0.1 micrometer-sized icy particles. These particles grow to pebble-sizes at a slower pace than the particles of super-particle 1, because the growth timescale (Eq. (13)) is longer for larger distances from the star (because the solids density goes down with increasing semi-major axis). At the time when the super-particle reaches the snowline by radial drift, the conditions for streaming instability outside the snowline are met, and part of its mass is converted to planetesimals before it crosses the snowline (Eq. (31)). There, the water ice content of the super-particle evaporates and the physical particles that it represents fragment to smaller sizes.

During the evolution of the solids, the position of the snowline changes. At first, it is pushed closer to the star as the pebble mass flux increases, after which it slowly recedes as the pebble mass flux decreases (Eq. (27)). When streaming instability outside the snowline starts (indicated with “SI outside snowline starts”), the pebble mass flux reaching the snowline is reduced because part of it is converted to planetesimals outside the snowline. When the pebble mass flux outside the snowline falls below the critical pebble mass flux for streaming instability, planetesimal formation outside the snowline stalls, and the snowline resumes to recede as the pebble mass flux delivered from the outer disk continues to decrease.

thumbnail Fig. 5

Cartoon of two possible trajectories through space and time of pebble super-particles (1 and 2). Both super-particles start out in the outer disk with a small typical particle mass. The physical particles represented by the super-particles grow until they reach a size where radial drift becomes faster than growth, and they start to drift inward. Planetesimal formation occurs outside the snowline as well as in the dry inner disk. See text for a more elaborate description.

Open with DEXTER

6 Fiducial model results

In this section we present results for our fiducial model, where we consider a protoplanetary disk around a solar-mass star. The parameter values for the fiducial model can be found in Table 1. We constrain the total disk mass Mdisk to 0.04 M. Together with the fixed values for α and gas, this constrains the outer radius of the disk rout, which for the fiducial model is 55 au.

Figure 6 shows the particle mass mp as function of semi-major axis r at different time points. The red line indicates the fragmentation limit on the mass (see Sect. 3.4.1) at the last time point plotted (5 × 105 years). After a few hundred years, particles interior to the snowline have grown to their fragmentation limit. In the outer part of the disk (≳10 au) the fragmentation limit is never reached; there the drift timescale is shorter than the growth timescale. The boundary between the fragmentation-limited and the drift-limited regions gets closer to the star in time, as the solids surface density goes down.

In Fig. 7 we plot the solids surface density profiles at different points in time. After a few thousand years, a “traffic jam”materializes interior to the snowline. This pile-up then spreads throughout the inner disk because an icy pebble mass flux from the outer disk keeps delivering material to the snowline region. At some point, as the pile-up is extending through the inner disk, the pebble mass flux that is delivered interior to the snowline starts to decrease, and the solids surface density interior to the snowline goes down as well. The decrease of the pebble mass flux is exacerbated by the streaming instability, because of which pebbles are converted to planetesimals outside the snowline. This is made clear in Fig. 8. The surface density profiles in Fig. 7 also show the transition from the fragmentation-dominated phase to the drift-dominated phase outside the snowline. For t < 105 yr, the region just outside the snowline is still fragmentation-dominated (solids surface density ∝ r−3∕2; Birnstiel et al. 2012). At t = 2 × 105 yr, the region outside the snowline has become drift-dominated, which is shown by an r−1 surface density power law.

Figure 8 shows the pebble mass flux that arrives at the snowline rsnow as a function of time (solid black line), as well as the pebble mass flux that arrives outside the planetesimal annulus (dashed black line) and the streaming instability threshold on the pebble mass flux (dotted line). The total mass in planetesimals is plotted in blue. The pebble mass flux first increases as particles grow (Fig. 6). The critical pebble mass flux needed for streaming instability varies with the stopping time of particles outside the snowline, and reaches a constant value when particles outside the snowline have reached their fragmentation limit, at about 103 years. Whenthe pebble mass flux exceeds the critical pebble mass flux, streaming instability operates and pebbles are being converted toplanetesimals at a rate given by Eq. (31). The pebble mass flux reaching the snowline is decreased to the critical pebble mass flux on the peak formation timescale (Sect. 4.2). The pebble mass flux that is still reaching the snowlineafter the enhancement has saturated, during the planetesimal formation phase, is equal to the critical pebble mass flux peb,crit that is needed to sustain the high solids-to-gas ratio outside the snowline (Sect. 4.2). Whenthe pebble mass flux that is delivered to the streaming instability region (outside rsnow; dashed line) drops below peb,crit, streaming instability shuts off and the pebble mass flux reaching rsnow is again equal to the pebble mass flux reaching the planetesimals annulus outside rsnow. This happens at around 3.2 × 104 years. The total yield in planetesimals is about ten Earth masses.

thumbnail Fig. 6

Particle mass mp as function of semi-major axis r at different points in time for the fiducial model (1a). The red line corresponds to the fragmentation limit on the mass (see Fig. 3) for the last time point plotted (t = 5 × 105 yr).

Open with DEXTER
thumbnail Fig. 7

Solids surface density profiles Σ at different points in time for the fiducial model (1a). The blue dots correspond to the mass of planetesimal super-particles pltsmls in Earth masses at the last time point plotted (t = 5 × 105 yr) (this is also the final total mass in planetesimals, because the planetesimal formation phase ends at ~ 3.2 × 104 yr).

Open with DEXTER
thumbnail Fig. 8

Pebble mass flux peb that arrives at the snowline rsnow in Earth masses per year, as a function of time t (solid black line), for the fiducial model (1a). The dashed black line corresponds to the pebble mass flux peb that arrives at the planetesimal annulus outside rsnow. During the planetesimal formation phase, the pebble mass flux at the snowline is constant and equal to the pebble mass flux delivered by the “peak” outside the snowline peb,crit (dotted line): the excess flux is converted to planetesimals. The total mass in planetesimals Mpltsmls is plotted inblue.

Open with DEXTER

7 Parameter study results

We investigate planetesimal formation in disks around stars of different masses. We scale the gas accretion rate gas with stellar mass M as , consistent with observations (Manara et al. 2017; Mulders et al. 2017). The viscous temperature Tvisc scales with , where M is the stellar mass (Frank et al. 2002). It then follows that . We use the following mass–luminosity relation: , with which the irradiation temperature Tirr scales linearly with M.

In our steady-state viscous gas disk model, the parameters gas and α determine the gas surface density profile Σgas (Eq. (1)). One then needs to either fix the total disk mass Mdisk or the outer disk radius rout to find the other. In constructing disk models around different stellar masses for our parameter study, we consider two scenarios. In the first one, we constrain Mdisk to 4% M, and the value of rout follows. In the second one, we take rout = 150 au, from which the value of Mdisk follows. We omit disks that would have Mdisk > 10% M or Mdisk < 1% M, as well as disks for which rout > 500 au or rout < 30 au.

In Table C.1 we list the parameter values for the different disk models, as well as the results (total masses in icy and dry planetesimals, and the time period during which streaming instability outside the snowline is operational (if it is operational at all)). We run models with stellar masses M of 0.1 M, 0.5 M, 1 M and 2 M (model IDs starting with 2, 3, 1, and 4, respectively). We also included models where the gas accretion rate gas varies with time; where we increasedthe metallicity Z; where we include pebble accretion; and where we increase the fragmentation velocity for icy particles vfrag,ice. These models are discussed in Sects. 7.2, 7.3, 7.4, and 7.5, respectively.

7.1 Varying stellar mass

Let us first compare the results of models 1a, 2b, 3a, and 4a to take a look at the stellar mass dependence for a fixed value of α =10−3. We find that models 1a, 2b, and 3a (corresponding to 1, 0.1, and 0.5 solar-mass stars, respectively) convert similar percentages of their initial solids reservoir to planetesimals (~ 10%). Model 4a (M = 2M) does not form any planetesimals. The critical pebbles-to-gas mass flux ratio peb,critgas that is required for streaming instability depends only weakly on the values of α and gas. We find however that the ratio of the actual pebbles-to-gas mass flux ratio just outside the snowline does depend on the stellar mass, through the snowline location. We can write as (35)

where here rsnow,out is the location just outside the snowline, Σpeb,snow,out is the pebble surface density at rsnow,out, and vp is the drift velocity of pebbles at rsnow,out. Making use of the relation Σgasgasα−1, and neglecting effects such as the time-evolution of the solids disk and the peak formation timescale, to zeroth order, we can approximate Eq. (35) as follows: (36)

where τ is the dimensionless stopping time of pebbles just outside the snowline (where particles are icy) and p is the power-law index of the pebbles surface density. In case the region outside the snowline is fragmentation-limited, as is the case in our simulations (except for model 1l; see Sect. 7.5); p = −3∕2. In the drift-limitedcase p = −1. Assuming that the temperature at the snowline is independent of stellar mass, the fragmentation limit on τ is proportional to α−1 (Birnstiel et al. 2012). The drift limit on τ is proportionalto with VK the Kepler velocity and cs the sound speed (Birnstiel et al. 2012), which boils down to if we again assume a snowline temperature independent of stellar mass. We then find if the region outside the snowline is drift-limited and for the fragmentation-limited case. In both cases, the pebbles-to-gas mass flux ratio that is important for planetesimal formation becomes smaller for more massive stars that have their snowline further away. Because the snowline is located further away in model 4a than in models 1a, 2b, and 3c, the critical pebble mass flux is (just) not reached in this disk and no planetesimals are formed.

The negative relation between and α (partly due to the set-up of the α-viscosity gas disks (Σgasgasα−1) as was already found in SO17), is clearly seen in Table C.1. For a given stellar mass, the lower α, the higher the total pebbles-to-planetesimal conversion efficiency. Model 4a has the lowest α value among all runs with M = 2 M, and comes closest to forming planetesimals: the pebble mass fluxes in the other disks in batch 4, which have higher α values, are even lower.

We have to keep in mind though that the zeroth-order dependencies on the stellar mass discussed above can be overshadowed by the effects of the gas disk. We have already discussed the effect of the value of α (which by no means needs to be constant throughout the disk as assumed in this work); but the transition from the Epstein to the Stokes drag regime; the compactness of the disk; deviations from the viscously-relaxed gas profile; the time when the fragmentation-limited region outside the snowline becomes drift-limited, etcetera, could play important roles as well. For example, one might notice that model 3a has a higher pebbles-to-planetesimals conversion efficiency than model 2b, which is not expected from the snowline-dependency argument alone. However, disk 3a is much more compact than disk 2b. This means that disk 2b has a lot of solids material in the outer disk that is not used for streaming instability, because at the time when pebbles that originate in the outer disk reach the snowline, the planetesimal formation phase has already ended. This reduces the conversion efficiency in model 2b. The compactness of a disk can also be a limiting factor to the total amount of planetesimals that can form. Model disk 3d is the same as model disk 3a except that it is larger. Model 3d forms more planetesimals than model 3a because a high enough pebble mass flux is sustained for a longer period of time (but again, the total conversion efficiency is lower). Therefore, even though the pebbles-to-gas mass flux ratio outside the snowline depends on the snowline location, we cannot directly translate this dependency to a simple scaling law for the pebbles-to-planetesimals conversion efficiency. Such secondary effects are precisely the reason why fast and versatile planetesimal formation models such as the one presented here are necessary.

7.2 Time-dependent accretion rate

Up to now, we have kept gas constant in our model runs. In reality, the gas accretion rate gas decreases on a viscous evolution timescale , where r1 is the characteristic radius of thedisk and ν1 is the viscosity at r1 (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). The evolution of the solids content of the disk proceeds faster than tvisc, and therefore a constant Σgas profile is justified. However, the location of the water snowline, which is an important quantity in our model, is highly sensitive to the temperature structure in the disk. The temperature in the inner disk, which is dominated by viscous heating, depends on gas. Therefore, a decrease in gas – however small – could lead to an inward movement of the water snowline during the planetesimal formation phase (Garaud & Lin 2007; Oka et al. 2011; Martin & Livio 2012). The viscous evolution timescale is shorter for larger α. Therefore, to investigate how viscous evolution affects our results, we run a model with a high α-value of 3 × 10−2, where the viscous temperature Tvisc goes down with decreasing gas (Sect. 2.2). The gas surface density profile Σgas is adjusted accordingly, but we do not account for viscous spreading of the disk (we keep rout fixed). This is not a problem, because it would not change the situation around the snowline and the inner disk. Because of the high α-value, we also increase the fragmentation velocity for icy particles vfrag,icy from 10 to 60 m s−1 (see also Sect. 7.5), to get non-negligible drift velocities. The time-dependent accretion rate gas(t) is given by (37)

where we take for the characteristic radius of the disk r1 = 50 au.

In Fig. 9 we plot the resulting temperature profile at different points in time, as well as the corresponding snowline locations for a given, arbitrary ice flux of .

We find that in this model no planetesimals form. We plot the solids surface density profiles at different time points in Fig. 10. Because of the large value for α, silicate particles interior to the snowline move with the gas velocity (the radial drift component to their velocity is negligible compared to the gas speed). Therefore, the solids surface density profile inside the snowline has the same power-law index as the gas profile. Outside the snowline, particles are fragmentation-limited, again due to the high α-value. This leads to an r−3∕2 profile outside the snowline. Even further away from the star, particles are drift-limited and the solids profile evolves toward an r−1 profile. If we would have chosen a smaller value for α in this time-dependent gas accretion model, such that planetesimals would form, viscous evolution would be slower and the position of the snowline would not change significantly during the planetesimal formation phase. Therefore, we conclude that viscous evolution does not play an important role during the planetesimal formation phase. Recently, Dra̧żkowska & Dullemond (2018) concluded that for α ≳ 10−4, the build-up stage of the disk (when material is still falling onto the disk) is also not important for planetesimal formation. Of course, changes in the disk conditions due to viscous evolution would matter for the subsequent stages of planetesimal–planetesimal mergers and embryo migration, which occur after planetesimal formation and may take place on timescales longer than the viscous evolution timescale.

thumbnail Fig. 9

Temperature T as function of semi-major axis r at different points in time, for model 1g (Table C.1), where the accretion rate starts out at gas = 5 × 10−8 M yr−1 and decreases on the viscous evolution timescale (see text for more details). Time is denoted in years, and the vertical dashed lines indicate the corresponding snowline locations for a constant, arbitrary ice flux .

Open with DEXTER
thumbnail Fig. 10

Solids surface density profiles Σ at different points in time, for model 1g (Table C.1). In this model the gas accretion rate is time-dependent and α = 3 × 10−2. No planetesimals have formed.

Open with DEXTER

7.3 Higher metallicity: icy and dry planetesimal formation

Figure 11 shows results for model run 1i, where we increased the metallicity Z from 1 to 2%. Planetesimals form both outside and inside the snowline, but the icy planetesimals dominate the total mass in planetesimals. During the streaming instability phase outside the snowline, the pebble mass flux that reaches the snowline is equal to the critical pebble mass flux (Sect. 4.2), just as in the fiducial model. Before and after the icy planetesimal formation phase, however, the pebble mass flux delivered to the snowline is higher in this model than in the fiducial model because of the increased metallicity. Interior to the snowline the pebbles pile up after their ice has evaporated and they have fragmented to smaller sizes. The solids-to-gas ratio just interior to the snowline is not large enough to trigger streaming instability. However, as the pile up extends through the inner disk, dry planetesimals are formed closer to the star. This is because the surface density profile in the inner disk tends to an r−3∕2 power law: particle sizes are limited by fragmentation and the radial velocity is not dominated by the gas accretion velocity (as is the case in Fig. 10, where α and hence the gas accretion velocity vgas are higher). An r−3∕2 solids surface density power law is steeper than the gas surface density profile, leading to an increasing solids-to-gas ratio with decreasing distance to the star. The result is dry planetesimal formation that peaks close to the star. This effect was also described inDra̧żkowska et al. (2016). We note that in our model, turbulent radial diffusion of particles is not accounted for (though it is implicitly included in our treatment of planetesimal formation outside the snowline; see Sect. 3.6), and might oppose a particle pile up in the inner disk. However, small particles in the inner disk move inward with the gas accretion velocity, and therefore small particles are removed from (and delivered to) the inner pile-up region on a timescale that is of the same order as the radial diffusion timescale (because vgas~ νr; Eq. (2)).

In the work of Dra̧żkowska & Alibert (2017) no dry planetesimals were formed interior to the snowline even at high metallicities. This difference with our results is due to the fact that the dimensionless stopping times τ of particles involved in dry planetesimal formation in our simulations are ~ 10−3, and Dra̧żkowska & Alibert (2017) did not allow for streaming instability by particles with stopping times smaller than 10−2.

The ratio of the total amount of icy planetesimals that forms in model 1i (Z = 0.02) compared to the amount that forms in model 1a (Z = 0.01) is much more than a factor 2. This is because in model 1a, the pebble mass flux does not exceed the pebble mass flux threshold by a lot (see Fig. 8). Therefore, if we increase Z by a factor 2, the “excess” pebble mass flux becomes larger by a factor much larger than 2.

In model runs 1j and 1k, we increase the “direct” planetesimal formation efficiency ϵSI to 0.1 and 1.0, respectively (see Sect. 4.1). We find that the amount of dry planetesimals that forms close to the star increases by a factor ~ 3 between model 1i and 1k, suggesting that the planetesimal formation efficiency is mildly dependent on the streaming instability efficiency.

thumbnail Fig. 11

Results for model 1i where Z = 0.02 and planetesimal formation occurs not only outside the snowline, but also interior to it. The masses of the icy planetesimal super-particles are plotted in blue; the masses of the water-poor planetesimal super-particles are plotted in brown. A massive annulus of planetesimals forms outside the snowline. Interior to the snowline, a traffic-jam effect leads to a high solids-to-gasratio in the inner disk such that streaming instability is triggered by dry pebbles as well.

Open with DEXTER

7.4 Including pebble accretion

When small planetesimals have merged to larger bodies – a process that we do not model in this work – they can start to grow more efficiently by accreting pebbles (Visser & Ormel 2016). In our current models, we do not have information about the individual planetesimal sizes, because for simplicity we treat planetesimals as super-particle “sinks” that are only characterized by a position and a total mass. However, if the planetesimals are residing in a narrow annulus, we can assume that the pebble accretion efficiency ϵPA is proportional to the total mass in planetesimals, without having any information on the individual planetesimals. This assumption is only valid in the three-dimensional regime (Ormel 2017), when the planetesimals are residing in a vertical layer that is less extended than that of the pebbles, and only when planetesimal eccentricities and inclinations are small. The three-dimensional pebble accretion rate is given by (Liu & Ormel 2018): (38)

where A3 = 0.39 is a numerical constant (Ormel & Liu 2018), qp is the mass ratio between the planetesimals and the star, and peb is the pebble flux just outside the planetesimal annulus.

We takeour fiducial model, in which a narrow annulus of planetesimals forms outside the snowline (Sect. 6), and now take pebble accretion into account. We assume that the pebble mass flux available to streaming instability outside the snowline peb,SI equals the mass flux that remains after the pebbles that are accreted by the already existing planetesimals have been removed: (39)

The total mass loss rate of pebble super-particles due to pebble accretion is the negative of Eq. (38) and the individual pebble super-particle mass loss rates are calculated according to Eq. (33).

The resulting pebble mass fluxes at the snowline and outside the planetesimal annulus, as well as the total mass in planetesimals, are plotted as a function of time in Fig. 12. The planetesimal formation phase, which is characterized by the pebble mass flux that reaches the snowline being constant, is much shorter than in the fiducial model (Fig. 8). This is because pebble accretion reduces the pebble mass flux available for streaming instability peb,SI, which quickly becomes smaller than the critical pebble mass flux peb,crit. Therefore, even though the pebble mass flux delivered from the outer disk is still considerably larger than peb,crit, planetesimal formation stops and pebble accretion takes over. Planetesimal formation is a self-limiting process: the more planetesimals are formed, the larger the pebble accretion efficiency becomes, which eventually prohibits the formation of new planetesimals. The final total mass in planetesimals is 89.9 Earth masses, which is much larger than in the case without pebble accretion.

Two key effects will alter these results. First, planetesimals can self-excite each other to highly eccentric and inclined orbits when the individual planetesimals are large (Levison et al. 2015), or when large embryos emerge as part of a runaway growth process (Kokubo & Ida 1998). In that case pebble accretion on the smaller planetesimals will terminate, because the relative motions between pebbles and planetesimals become too large (Liu & Ormel 2018). Second, the eventual planetary embryos that emerge as the result of the combined planetesimal and pebble accretion are likely to migrate out of their birth zone due to type-I migration (Tanaka et al. 2002). Both effects will suppress the amount of planetesimals formed compared to the above estimate. A thorough analysis of the efficacy of this process is underway (Liu et al. 2018). Another effect of a large embryo crossing the snowline after the planetesimal formation phase has terminated could be that the pebble mass flux reaching the snowline becomes large enough for streaming instability again, such that a new generation of planetesimals forms.

thumbnail Fig. 12

Results for model 1 h, which is the fiducial model including pebble accretion. The pebble mass flux arriving at the snowline is plotted by the solid black line; the pebble mass flux outside the planetesimal annulus is plotted by the dashed black line. The dotted black line corresponds to the streaming instability threshold pebble mass flux. The total mass in planetesimals including pebble accretion is given by the blue line. A little after 104 yr, the pebble accretion efficiency is 100% and no pebbles reach the snowline anymore.

Open with DEXTER

7.5 Increasing the fragmentation threshold for icy particles

Up to now, we have fixed the fragmentation threshold velocity for icy particles at vfrag,icy = 10 m s−1. However, some molecular dynamics simulations suggest icy particles fragment only at higher velocities (Dominik & Tielens 1997; Wada et al. 2009, 2013). We therefore explore if and how the results change if we increase the icy fragmentation threshold from 10 m s−1 to 60 m s−1. Because of the higher fragmentation threshold, the sizes of particles outside the snowline are now limited by radial drift rather than fragmentation. In Fig. 13 the pebble mass fluxes peb outside the snowline (dashed black line) and at the snowline (solid black line), as well as the critical pebble mass flux required for triggering streaming instability, are plotted as a function of time. First, as the stopping time of pebbles outside the snowline increases due to particle growth, the pebble mass flux increases and the critical pebble mass flux (which depends on the stopping time) increases as well. In this case particles outside the snowline reach dimensionless stopping times τ ~ 0.3 at about 103 yr, after which the critical pebble mass flux decreases. At around 6 × 103 yr the stopping time at the snowline has reached a maximum of about 0.6, after which it gradually become smaller again because of depletion of solids mass outside the snowline due to radial drift (see also Lambrechts & Johansen 2014), and the threshold for streaming instability on the pebble mass flux increases again. The total mass in planetesimals as a function of time is plotted in blue. The final total planetesimal mass is 17.1 Earth masses, which is larger than the planetesimals yield of the fiducial model. This is caused by the fact that the pebble mass flux outside the snowline reaches larger values in this model than in the fiducial model due to the higher stopping times, as well as by the fact that the pebble mass flux required for streaming instability is for a large period of time reduced compared to the fiducial model, again because of the higher stopping times. All in all this leads to more excess pebble mass flux outside the snowline that can be converted to planetesimals. Note also that the shape of the pebble mass flux as a function of time is different in the drift-limited case (Fig. 13) than in the fragmentation-limited case (Fig. 8).

thumbnail Fig. 13

Results for model 1l, which is the fiducial model but with icy fragmentation threshold vfrag,icy = 60 m s−1. In solid black: the pebble mass flux peb that reachesthe snowline; in dashed black: the pebble mass flux outside the snowline that is available for planetesimal formation outside the snowline; and in dotted black the critical pebble mass flux needed for streaming instability outside the snowline, all as a function of time. The total mass in planetesimals is plotted in blue.

Open with DEXTER

8 Discussion

We find that lower-mass stars are more efficient at forming planetesimals than higher-mass stars, that is, disks around low-mass stars tend to convert a larger fraction of their initial solids content to planetesimals. This is because in our model, planetesimal formation works by virtue of a high flux of pebbles through the disk. The pebbles-to-gas mass flux ratio, which dictates planetesimal formation outside the snowline, depends inversely on snowline location, which for lower mass stars is closer in than for higher-mass stars that have hotter disks. The pebbles-to-planetesimals conversion efficiency also depends on other disk parameters that may vary across the stellar mass range, however, such as the outer disk radius, the turbulence parameter α, and the metallicity. A spread in disk properties (e.g., Ansdell et al. 2017, 2018; Mulders et al. 2017; Manara et al. 2017; Tazzari et al. 2017) might obscure any observable correlation between stellar mass and planetesimal formation efficiency.

For a given stellar mass and gas accretion rate, the total amount of planetesimals that forms depends sensitively on the value of the turbulence parameter α. This is partly due to the nature of the viscous gas disk model: for a given gas accretion rate, the gas surface density is larger for lower α. In our model this means that for a given stellar mass and metallicity, a lower α leads to a larger pebble mass flux. This conclusion was already made in Schoonenberg & Ormel (2017). Another quantity that proved to be important in our models is the fragmentation velocity threshold, because it determines the drift velocity of pebbles and hence the pebble flux. Therefore, more stringent laboratory constraints on the fragmentation threshold in the snowline region (including snowline-specific effects such as sintering), as well as constraints on particle sizes around observed snowlines (such as around FU Ori objects; Banzatti et al. 2015; Cieza et al. 2016; Schoonenberg et al. 2017), would help narrow down the parameter space for our models. Additionally, recent studies predict observable signatures of drifting pebbles in the gas-phase abundances of CO and CO2 (Booth et al. 2017; Bosman et al. 2018; Krijt et al. 2018), which could constrain the pebble mass flux observationally.

Because the crucial factor in our planetesimal formation model is a high pebble flux, planetesimals form in the early stages of disk evolution (in the first ~ 105 yr) when the pebble mass flux is still high. A side effect of fast planetesimal formation is the well-known “radial drift” problem: depletion of the solids content of the disk on short timescales is in disagreement with observations, which show that small solid particles are present in the outer regions of protoplanetary disks at later ages. A solution to this problem could be to take into account a particle size distribution. In this work, we assume a mono-disperse particle size distribution at each point in space and time, and thereby only focus on the large particles. Large particles dominate the mass, but small particles that do not drift significantly might be present at all times. This could also be the clue to late planetesimal formation: while in our model only early-stage planetesimal formation is possible due to the depletion of solids over time, other works have suggested models in which late planetesimal formation occurs in the outer disk at a late evolutionary disk stage due to gas depletion (Carrera et al. 2017). Also, the early formation ofa protoplanet could halt radial drift (Kobayashi et al. 2012; Morbidelli et al. 2016). Finally, another possible solution to the radial drift problem in the context of our model is choosing a larger outer disk radius rout, which leadsto longer drift timescales in the outer regions of the disk.

The Lagrangian smooth-particle model presented in this paper can be used for different research directions that we have not explored yetin this paper. Characteristics of individual (groups of) particles can be naturally followed during their evolution in time and through the disk. In this work, we have focused on the water content, but any other compositional information could be included as well. For example, the D/H ratio of water depends on the local conditions in the evolving protoplanetary disk (Yang et al. 2013), and therefore measurements of the D/H content of bodies in the solar system might provide information on where they formed (Hallis et al. 2015; Hallis 2017). An interesting direction for further research using our model would therefore be to focus on theD/H ratio of water in planetesimals by keeping track of where and when in the disk their water is incorporated. One could also think of different disk structures than the ones we have considered in this paper. Axisymmetric features such as rings and gaps seem to be ubiquitous in the gas and dust of observed disks (e.g., ALMA Partnership 2015; Akiyama et al. 2015; Andrews et al. 2016; Isella et al. 2016; Nomura et al. 2016; Fedele et al. 2018). In our model scenario, planetesimal formation occurs at earlier stages than the ages of the observed disks. However, in principle our model could deal with any disk structure. It would be interesting to investigate what the consequences would be for planetesimal formation if one assumes a disk with radial pressure bumps rather than a smooth diskas we did in the current work. Such structures could also alleviate the radial drift problem (Pinilla et al. 2012).

The formation of planetesimals from dust is clearly only the first piece of the planet formation puzzle. The next step is how planetesimals grow into protoplanets and beyond, for which we need an N-body model that follows the interactions between planetesimals, their growthby mutual collisions and pebble accretion, and their migration through the disk. The combination of the Lagrangian planetesimal formation model presented in this paper with such an N-body code will allow us to model the formation of entire planetary systems from A to Z while automatically following composition, which is a research direction we are planning to pursue. Our first target is the H2O fraction of the TRAPPIST-1 system (Schoonenberg et al., in prep.).

9 Conclusions

Our main findings can be summarized as follows:

  • 1.

    Planetesimals form early, when the pebble mass flux is still high. In the context of our model, in this early formation phase the migration of the snowline due to a decreasing gas accretion rate is not important.

  • 2.

    Planetesimals form preferentially in a narrow annulus outside the snowline. Even if rocky planetesimals are formed interior to the snowline, icy planetesimals dominate the total planetesimal mass.

  • 3.

    Pebble accretion leads to self-limiting planetesimal formation (the more planetesimals form, the higher the pebble accretion efficiency and the smaller the pebble mass flux available to form new planetesimals). However, taking into account migration and planetesimal–planetesimal scattering will change this picture.

  • 4.

    The planetesimal formation efficiency depends on the location of the water snowline. The cooler the disk (the closer-in the water snowline), the higher the pebbles-to-planetesimals conversion efficiency. Therefore, in general, we find that low-mass stars are better at producing planetesimals than high-mass stars. However, this result also depends on other factors such as the compactness of disks, and may change for gas disks that deviate from a steady-state α-viscosity disk where the gas surface density depends on a constant value of α as we assumed in this work.

Acknowledgements

D.S. and C.W.O are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). S.K. acknowledges support from NASA through Hubble Fellowship grant HST-HF2-51394 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. D.S. would like to thank Carsten Dominik for useful discussions and Shigeru Ida for pointing out a typographical error in Schoonenberg & Ormel (2017). We thank Beibei Liu, Joanna Dra̧żkowska, and the anonymous referee for constructive feedback on the manuscript.

Appendix A Validation

We test our Lagrangian smooth-particle model of dust growth and radial drift against analytical results from the literature.

A.1 Drift-only

First, we check if our model reproduces the analytical predictions of Youdin & Shu (2002) for the surface density evolution of a single particle size (0.5 cm) without growth. We take a gas disk with surface density profile and temperature profile . We consider an initial solids surface density profile Σ(r, 0) that has a cutoff at outer radius with r1 = 1 au if r < rout and Σ(r, 0) = 0 otherwise. Following the derivation in Sect. 4.1 of Youdin & Shu (2002), we then find that Σ evolves with time t as (A.1)

where d is defined through the radial dependence of the drift velocity vdr (the expression for which we copy from Youdin & Shu 2002) for a constant particle size: vdrrd. For our temperature profile Tr−0.5 and gas surface density profile Σgasr−3∕2 we find d = 3∕2. ri (r, t) is the initial location of a particle that ends up at radius r at time t, which is given by (A.2)

In Fig. A.1, we compare Eq. (A.1) with our numerical results at different points in time, demonstrating that our model reproduces the analytical results presented by Youdin & Shu (2002).

A.2 Including simple growth

We next turn to Lambrechts & Johansen (2014, hereafter: LJ14) to compare the results of our model including particle growth with their analytical results. We take the same disk model as LJ14 and plot the resulting surface densities at different points in time in Fig. A.2. LJ14 introduced the term “pebble front” for the distance from the star up to which particles have grown to pebble-sizes and are drifting inwards. This distance, calculated by their Eq. (1), is plotted by the vertical dashed lines in Fig. A.2. Rather than a razor-sharp transition between the pebble region to the dust region as in the analytical work of LJ14, our numerical results show a smooth connection between the two regions, the radius of which is in agreement with LJ14. The pebble front radius that follows from our simulation progresses a little bit faster in time than the pebble front radius in LJ14, which is due to the fact that in their work, a constant value of ξ (the number of growth e-foldings needed to grow to pebble-sizes), was used regardless of semi-major axis, whereas in our model it is not. The slope of the surface density interior to the pebble front falls off as r−3∕4, as was also found by LJ14.

thumbnail Fig. A.1

Solids surface density Σ as function of radial distance from the star r, plotted at different time points. Solid lines correspond to the analytical predictions from Youdin & Shu (2002, Eq. (A.1)), and scattered points correspond to our model results. The physical size of particles is 0.5 cm.

Open with DEXTER
thumbnail Fig. A.2

Reproduction of Fig. 1 of Lambrechts & Johansen (2014) with our model. Gray lines correspond to solids surface densities at different points in time (labels denote time in years). Gas surface density profiles are plotted in blue for the same time points. The vertical dashed lines indicate the radial extent of the pebble front defined in Eq. (1) of Lambrechts & Johansen (2014).

Open with DEXTER

Appendix B Updatedsemi-analytic model for solids enhancement outside the snowline

Appendix C of Schoonenberg & Ormel (2017) contains several minor errors or unclarified approximations:

  • the term was omitted in the top line of Eq. (C.10);

  • in Eq. (C.10), the normalized of Eq. (C.7) was substituted. However, the difference between the particle and gas diffusivities was not accounted for;

  • Eq. (C.11) also implicitly assumes Dgas = Dp.

As a result the semi-analytic model will give an erroneous result for particles reaching dimensionless stopping time approaching unity (at the location of the ice peak), since in that case Dp < Dgas. For particles τ ≪ 1, however, the model outlined in SO17 remains correct.

In this appendix, we briefly re-derive the key expressions.

Starting from Eq. (C.8) of SO17, we normalize lengths to rsnow, surface density to Σsnow and velocities to Dgasrsnow: (B.1)

The left-hand side represents the mass flux of the ice; ~ denote non-dimensional quantities. The ice surface density (Σ) is normalized by the quantity Σsnow – the surface density in H2O vapor at the snowline rsnow – the pebble drift velocity is normalized by Dgasrsnow and the radius r is normalized by rsnow. The first term on the right-hand side (RHS) is given by Eq. (C.7) of SO17. The second term on the RHS is equals + bgas = 3∕2 (because , the total mass flux in ice, is directed inwards) in an α-disk. Applying these, we obtain (B.2)

where x = rrsnow − 1 as in SO17 and aeq is a constant that enters the expression for the equilibrium vapor density of H2 O beyond the snowline. Getting rid of ~ notation and putting Dgas = rsnow = 1: (B.3)

An analytical solution to this differential equation exists: (B.4)

Of particular importance is the location where this function peaks (Σ′ = 0). The maximum corresponding to Σ(x) is located at (B.5)

and the corresponding value of the (normalized) surface density is (B.6)

which follows from Eq. (B.3).

Appendix C Additional table

Table C.1.

Parameter values and results for different model runs.

References

  1. Akiyama, E., Muto, T., Kusakabe, N., et al. 2015, ApJ, 802, L17 [NASA ADS] [CrossRef] [Google Scholar]
  2. ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
  3. Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
  6. Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bai, X.-N., & Stone, J. M. 2010, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
  8. Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15 [NASA ADS] [CrossRef] [Google Scholar]
  9. Benz, W. 2000, Space Sci. Rev., 92, 279 [NASA ADS] [CrossRef] [Google Scholar]
  10. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Blum, J., & Münch, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
  12. Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
  13. Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bosman, A. D., Tielens, A. G. G. M., & van Dishoeck, E. F. 2018, A&A, 611, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [NASA ADS] [CrossRef] [Google Scholar]
  17. Chapman, S., & Cowling, T. 1970, The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge Mathematical Library (Cambridge: Cambridge University Press) [Google Scholar]
  18. Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, 258 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490 [NASA ADS] [CrossRef] [Google Scholar]
  21. Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
  22. Dra̧żkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Dra̧żkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Dra̧żkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn. (Cambridge: Cambridge University Press), 398 [Google Scholar]
  28. Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [NASA ADS] [CrossRef] [Google Scholar]
  29. Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
  30. Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hallis, L. J. 2017, Phil. Trans. R. Soc. London, Ser. A, 375, 20150390 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hallis, L. J., Huss, G. R., Nagashima, K., et al. 2015, Science, 350, 795 [NASA ADS] [CrossRef] [Google Scholar]
  33. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ida, S., & Guillot, T. 2016, A&A, 596, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
  36. Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
  37. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  38. Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
  39. Kenyon, S. J.,& Hartmann, L. 1987, ApJ, 323, 714 [NASA ADS] [CrossRef] [Google Scholar]
  40. Kobayashi, H., Ormel, C. W., & Ida, S. 2012, ApJ, 756, 70 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
  42. Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hicks, D. L., & Liebrock, L. 2000, Appl. Math. Comput., 112, 63 [Google Scholar]
  45. Laibe, G., Gonzalez, J.-F., Fouchet, L., & Maddison, S. T. 2008, A&A, 487, 265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322 [NASA ADS] [CrossRef] [Google Scholar]
  48. Lichtenegger, H. I. M., & Komle, N. I. 1991, Icarus, 90, 319 [NASA ADS] [CrossRef] [Google Scholar]
  49. Liu, G., & Liu, M. 2003, Smoothed Particle Hydrodynamics: A Meshfree Particle Method (Singapore: World Scientific) [Google Scholar]
  50. Liu, B., & Ormel, C. W. 2018, A&A, 615, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Liu, B., Ormel, C. W., & Johansen, A. 2018, A&A, submitted [Google Scholar]
  52. Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
  53. Lorek, S., Gundlach, B., Lacerda, P., & Blum, J. 2016, A&A, 587, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  55. Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Martin, R. G., & Livio, M. 2012, MNRAS, 425, L6 [NASA ADS] [CrossRef] [Google Scholar]
  57. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
  58. Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 [NASA ADS] [CrossRef] [Google Scholar]
  59. Mulders, G. D. 2018, Handbook of Exoplanets, eds. H. Deeg & J. Belmonte (Cham: Springer) [Google Scholar]
  60. Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130 [NASA ADS] [CrossRef] [Google Scholar]
  61. Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31 [NASA ADS] [CrossRef] [Google Scholar]
  62. Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
  63. Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJ, 819, L7 [NASA ADS] [CrossRef] [Google Scholar]
  64. Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
  65. Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
  66. Ormel, C. W. 2017, Astrophys. Space Sci. Lib., 445, 197 [NASA ADS] [CrossRef] [Google Scholar]
  67. Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Ormel, C. W., & Liu, B. 2018, A&A, 615, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  72. Ros, K., & Johansen, A. 2013, A&A, 552, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Safronov, V. S. 1969, in Evolution of the Protoplanetary Cloud and Formation of Earth and the Planets, ed. V. S. Safronov (Moscow: Nauka. Transl. 1972 NASA Tech. F-677) [Google Scholar]
  74. Saito, E., & Sirono, S.-i. 2011, ApJ, 728, 20 [NASA ADS] [CrossRef] [Google Scholar]
  75. Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Schoonenberg,D., Okuzumi, S., & Ormel, C. W. 2017, A&A, 605, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  79. Sirono, S. 1999, A&A, 347, 720 [NASA ADS] [Google Scholar]
  80. Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
  81. Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Visser, R. G., & Ormel, C. W. 2016, A&A, 586, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
  84. Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  86. Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (Slona: Almqvist och Wiksell), 211 [Google Scholar]
  87. Yang, L., Ciesla, F. J., & Alexander, C. M. O. 2013, Icarus, 226, 256 [NASA ADS] [CrossRef] [Google Scholar]
  88. Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
  90. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
  91. Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
  92. 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]

1

We use the term “dust” for solid particles that are well-coupled to the gas, and the term “pebbles” for particles that have a non-negligible radial drift velocity (Stokes numbers of ~ 10−3−101; see Sect. 3.4). Throughout the paper, the term “pebble super-particle” is used instead of “dust/pebble super-particle”.

2

There is a typo in Eq. (34) of SO17 for the steady-state advection-only vapor surface density profile. That equation should be the same as the current Eq. (26).

3

If we take the accretion rate gas to be a decreasing function of time and the snowline is located in the viscosity-dominated temperature region, the decrease in accretion heating is an additional effect that pushes the snowline inward (Garaud & Lin 2007; Oka et al. 2011; Martin & Livio 2012). This is discussed in Sect. 7.2.

4

The semi-analytic model of SO17 is not valid for dimensionless stopping times τ > 1. We have updated it so that it can deal with τ > 1 values as well. We present this update in Appendix B.

5

One could imagine an “episodic” nature of streaming instability rather than the “continuous” approach we use here. In the episodic scenario, once streaming instability outside the snowline has been triggered, the annulus of pebbles outside the snowline is emptied and need to be refilled by inward-drifting pebbles before streaming instability is triggered again, leading to episodic bursts of planetesimal formation. If the width of the peak Wpeak is much larger than ηr (the distance a drifting pebble traverses during one vertical settling timescale) and all pebbles in the peak participate in planetesimal formation, the episodic description may be better than the continuous one. However, we checked that the timescales involved in emptying and refilling the snowline region with pebbles are much shorter than the planetesimal formation phase, which justifies the continuous approach.

6

To be specific, Eq. (32) reads: (33)

where Npeb is the number of pebble super-particles and WGauss is given by (34)

The term [(1 − fice,out) + fice] in Eq. (33) takes into account the fact that in our model the location of the solids enhancement outside the snowline rpeak can be at a distance from the star where fice < fout, again because we do not follow the transport of water vapour explicitly in this work.

All Tables

Table 1

Input parameters, their symbols, and their fiducial values.

Table C.1.

Parameter values and results for different model runs.

All Figures

thumbnail Fig. 1

Gas disk profiles for our fiducial model (Table 1). The black line corresponds to the gas surface density profile Σgas (Eq. (1); the solid blue line corresponds to the temperature profile T, which receives contributions from viscous heating (dashed blue line) and stellar irradiation (dotted blue line; Eq. (5)).

Open with DEXTER
In the text
thumbnail Fig. 2

Procedure cartoon of the resampling concept. When the separation between neighbouring pebble super-particles becomes toolarge, the number and configuration of particles is reset such that the resolution is again equal to that of the initial set-up. The total mass in pebble super-particles is conserved and the characteristics (such as water ice fraction, indicated bycolor here) are sampled from the situation before resampling.

Open with DEXTER
In the text
thumbnail Fig. 3

Fragmentation mass limit mfrag (black line) for our fiducial model with initial ice mass fraction profile fice (blue line). The snowline is located at around 2.4 au. Because of the different stickiness between icy and silicate particles, the fragmentation velocity, and therefore the fragmentation mass limit, varies rapidly across the snowline.

Open with DEXTER
In the text
thumbnail Fig. 4

Initial pebble ice fraction (blue solid line) as function of semi-major axis r for the fiducial model. The saturated water vapour profile Σeq is plotted by the black dashed line, and ΣZ,a (see text for more details) is given by the black solid line. The snowline is located where Σeq crosses ΣZ,a, in this case at rsnow ~ 2.4 au (blue dotted line).

Open with DEXTER
In the text
thumbnail Fig. 5

Cartoon of two possible trajectories through space and time of pebble super-particles (1 and 2). Both super-particles start out in the outer disk with a small typical particle mass. The physical particles represented by the super-particles grow until they reach a size where radial drift becomes faster than growth, and they start to drift inward. Planetesimal formation occurs outside the snowline as well as in the dry inner disk. See text for a more elaborate description.

Open with DEXTER
In the text
thumbnail Fig. 6

Particle mass mp as function of semi-major axis r at different points in time for the fiducial model (1a). The red line corresponds to the fragmentation limit on the mass (see Fig. 3) for the last time point plotted (t = 5 × 105 yr).

Open with DEXTER
In the text
thumbnail Fig. 7

Solids surface density profiles Σ at different points in time for the fiducial model (1a). The blue dots correspond to the mass of planetesimal super-particles pltsmls in Earth masses at the last time point plotted (t = 5 × 105 yr) (this is also the final total mass in planetesimals, because the planetesimal formation phase ends at ~ 3.2 × 104 yr).

Open with DEXTER
In the text
thumbnail Fig. 8

Pebble mass flux peb that arrives at the snowline rsnow in Earth masses per year, as a function of time t (solid black line), for the fiducial model (1a). The dashed black line corresponds to the pebble mass flux peb that arrives at the planetesimal annulus outside rsnow. During the planetesimal formation phase, the pebble mass flux at the snowline is constant and equal to the pebble mass flux delivered by the “peak” outside the snowline peb,crit (dotted line): the excess flux is converted to planetesimals. The total mass in planetesimals Mpltsmls is plotted inblue.

Open with DEXTER
In the text
thumbnail Fig. 9

Temperature T as function of semi-major axis r at different points in time, for model 1g (Table C.1), where the accretion rate starts out at gas = 5 × 10−8 M yr−1 and decreases on the viscous evolution timescale (see text for more details). Time is denoted in years, and the vertical dashed lines indicate the corresponding snowline locations for a constant, arbitrary ice flux .

Open with DEXTER
In the text
thumbnail Fig. 10

Solids surface density profiles Σ at different points in time, for model 1g (Table C.1). In this model the gas accretion rate is time-dependent and α = 3 × 10−2. No planetesimals have formed.

Open with DEXTER
In the text
thumbnail Fig. 11

Results for model 1i where Z = 0.02 and planetesimal formation occurs not only outside the snowline, but also interior to it. The masses of the icy planetesimal super-particles are plotted in blue; the masses of the water-poor planetesimal super-particles are plotted in brown. A massive annulus of planetesimals forms outside the snowline. Interior to the snowline, a traffic-jam effect leads to a high solids-to-gasratio in the inner disk such that streaming instability is triggered by dry pebbles as well.

Open with DEXTER
In the text
thumbnail Fig. 12

Results for model 1 h, which is the fiducial model including pebble accretion. The pebble mass flux arriving at the snowline is plotted by the solid black line; the pebble mass flux outside the planetesimal annulus is plotted by the dashed black line. The dotted black line corresponds to the streaming instability threshold pebble mass flux. The total mass in planetesimals including pebble accretion is given by the blue line. A little after 104 yr, the pebble accretion efficiency is 100% and no pebbles reach the snowline anymore.

Open with DEXTER
In the text
thumbnail Fig. 13

Results for model 1l, which is the fiducial model but with icy fragmentation threshold vfrag,icy = 60 m s−1. In solid black: the pebble mass flux peb that reachesthe snowline; in dashed black: the pebble mass flux outside the snowline that is available for planetesimal formation outside the snowline; and in dotted black the critical pebble mass flux needed for streaming instability outside the snowline, all as a function of time. The total mass in planetesimals is plotted in blue.

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

Solids surface density Σ as function of radial distance from the star r, plotted at different time points. Solid lines correspond to the analytical predictions from Youdin & Shu (2002, Eq. (A.1)), and scattered points correspond to our model results. The physical size of particles is 0.5 cm.

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

Reproduction of Fig. 1 of Lambrechts & Johansen (2014) with our model. Gray lines correspond to solids surface densities at different points in time (labels denote time in years). Gas surface density profiles are plotted in blue for the same time points. The vertical dashed lines indicate the radial extent of the pebble front defined in Eq. (1) of Lambrechts & Johansen (2014).

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.