© 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. 2–4 and provide a summary of the model in Sect. 5. The results are discussed in Sects. 6–7 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 c_s the sound speed and Ω the Keplerian orbital frequency. The gas moves inward at a speed |v_gas| = 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.

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

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 T_visc∝ r^−3∕4 (Frank et al. 2002). Stellar irradiation results in a temperature profile in the outer disk that goes as T_irr ∝ r^−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 T_visc and T_irr: (5)

which then defines the isothermal sound speed profile c_s: (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, T_visc, and T_irr 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 H_gas 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-particles¹, 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 m_p (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(x − x′, h) is the kernel, for which we take (Hicks & Liebrock 2000; Liu & Liu 2003): (10)

where h is the smoothing length, Δ x = |x − x′| 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 (|x − x′| > 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 x_i of super-particle i, we then find (11)

where M_j is the mass of supporting super-particle j and h_i is the smoothing length of super-particle i. We treat the smoothing length h_i as a variable, demanding that at each super-particle location and time h_i 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 x_i (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 r_out. 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 f_ice of particles is 50% (f_ice,out = 0.5) (Lodders 2003; Morbidelli et al. 2015). Interior to the water snowline, the solids consist of pure silicate (f_ice,in = 0). The location of the water snowline and the transition between f_ice,out and f_ice,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 r_i between the inner edge of the disk r_in and the outer edge r_out. The kth super-particle then gets assigned a mass (corresponding to the total mass in the kth annulus) and an initial position x_k = (r_k+1 − r_k)∕2. If the initial solids surface density profile is proportional to r⁻¹ (Σ_solids(r) = Σ₀r⁻¹ where Σ₀ is constant) and the spacing between the annuli is linear, each super-particle has the same initial mass m = 2πΔrΣ₀ 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.

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.

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 s_p and mass m_p. Within each super-particle then, particles can coagulate, increasing the individual particle mass m_p of that super-particle. The mass growth rate dm_p∕dt assuming perfect sticking is (Krijt et al. 2016) (13)

with H_solids the solids scale height, σ_col the collisional cross section (which we take equal to the geometrical cross section), and v_rel the relative velocity between particles. As long as the particles are small and well-coupled to the gas, H_solids = H_gas. 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 l_mfp, 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 l_mfp 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⁻¹⁵ cm² as appropriate for hydrogen (Chapman & Cowling 1970). The dimensionless stopping time is given by (16)

where s_p is the particle radius, ρ_∙,p is the particle internal density, and v_th is the thermal velocity of the gas molecules, defined as . The particle radius and internal density are determined from the particle mass m_p and water mass fraction f_ice:

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

For therelative velocity between individual particles within a super-particle v_rel, 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 v_rel,turb dominates. In the so-called intermediate regime – where the turnover time scale of the smallest eddies ( where Re_T 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 – v_rel,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 v_rel between particles becomes larger than a certain fragmentation threshold v_frag, particles start to fragment rather than coagulate when they collide. Concerning v_frag 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 v_frag,sil = 3 m s⁻¹ (Blum & Münch 1993); and for icy particles, we take v_frag,ice = 10 m s⁻¹ (Sirono 1999; Gundlach & Blum 2015). The fragmentation threshold v_frag for particles with ice fraction f_ice is then given by (20)

where we simply coupled the two fragmentation velocities according to the ice mass fraction f_ice (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 f_ice is arbitrary but unimportant, since f_ice changes rapidly across the water snowline (Dra̧żkowska & Alibert 2017).

If v_rel > v_frag, we instantaneously set the particle mass m_p to the fragmentation mass limit m_frag. In Fig. 3 we have plotted m_frag as a function of semi-major axis r for our fiducial model. The assumption that fragmentation is instantaneous is justified by the fact that m_frag 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, m_frag decreases for inward-drifting particles (see Fig. 3). In this region therefore, inward-drifting particles that have a mass that is limited to m_frag 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.

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.

3.4.2 Radial drift

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

where v_gas is negative (the gas is moving inward as well), ξ is the midplane solids-to-gasvolume density ratio, and ηv_K is the difference between the azimuthal motion of the gas disk and the Keplerian velocity v_K: (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 r_snow as the radius interior to which the water vapour pressure P_vap drops below the saturated (equilibrium) pressure P_eq. Outside r_snow, the water vapour distribution always follows the saturated profile, which is given by the Clausius–Clapeyron equation: (24)

where T_a and P_eq,0 are constants depending on the species. For water, T_a = 6062 K and P_eq,0 = 1.14 × 10¹³ g cm⁻¹ s⁻² (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 SO17². 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 r_snow 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 time³. 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 f_ice,outΣ_solids(r, t = 0).

We define the ice fraction of pebbles f_ice 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 f_ice 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 f_ice for our fiducial model (Table 1).

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

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⁻³ (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 N_pltsml 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 dM_pl,SI,direct∕dt: (30)

where M_peb,SI is the mass of the pebble super-particles that meet the requirement for SI, and t_settle = 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 times⁴. 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: f_ice,pltsml = f_ice,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 t_peak,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 r_peak, and the snowline location r_snow. The width of the peak W_peak = r_peak − r_snow 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 dM_{pl,SI,snowline}∕dt on a peak formation timescale t_peak,buildup = W_peak∕|v_gas|. The value ofdM_{pl,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 planetesimals⁵. The mass of the planetesimal super-particle closest to the peak location r_peak (which is given by the semi-analytic model) is increased at a rate dM_{pl,SI,snowline}∕dt.

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

where the fraction f(r_i) < 1 depends on the distance between r_i and the peak location⁶.

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 − f_ice,out) + f_ice]; 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 H₂ 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.

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.

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 M_disk to 0.04 M_⊙. Together with the fixed values for α and Ṁ_gas, this constrains the outer radius of the disk r_out, which for the fiducial model is 55 au.

Figure 6 shows the particle mass m_p 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 × 10⁵ 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 < 10⁵ yr, the region just outside the snowline is still fragmentation-dominated (solids surface density ∝ r^−3∕2; Birnstiel et al. 2012). At t = 2 × 10⁵ yr, the region outside the snowline has become drift-dominated, which is shown by an r⁻¹ surface density power law.

Figure 8 shows the pebble mass flux that arrives at the snowline r_snow 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 10³ 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 r_snow; dashed line) drops below Ṁ_peb,crit, streaming instability shuts off and the pebble mass flux reaching r_snow is again equal to the pebble mass flux reaching the planetesimals annulus outside r_snow. This happens at around 3.2 × 10⁴ years. The total yield in planetesimals is about ten Earth masses.

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

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

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.

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 T_visc 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 T_irr 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 M_disk or the outer disk radius r_out 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 M_disk to 4% M_⋆, and the value of r_out follows. In the second one, we take r_out = 150 au, from which the value of M_disk follows. We omit disks that would have M_disk > 10% M_⋆ or M_disk < 1% M_⋆, as well as disks for which r_out > 500 au or r_out < 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 v_frag,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⁻³. 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,crit∕Ṁ_gas 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 r_snow,out is the location just outside the snowline, Σ_peb,snow,out is the pebble surface density at r_snow,out, and v_p is the drift velocity of pebbles at r_snow,out. Making use of the relation Σ_gas ∝Ṁ_gasα⁻¹, 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 α⁻¹ (Birnstiel et al. 2012). The drift limit on τ is proportionalto with V_K the Kepler velocity and c_s 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 (Σ_gas ∝Ṁ_gasα⁻¹) 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 r₁ is the characteristic radius of thedisk and ν₁ is the viscosity at r₁ (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). The evolution of the solids content of the disk proceeds faster than t_visc, 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⁻², where the viscous temperature T_visc 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 r_out 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 v_frag,icy from 10 to 60 m s⁻¹ (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 r₁ = 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⁻¹ 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⁻⁴, 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.

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 .

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.

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 v_gas 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 v_gas~ ν∕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⁻³, and Dra̧żkowska & Alibert (2017) did not allow for streaming instability by particles with stopping times smaller than 10⁻².

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.

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.

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 A₃ = 0.39 is a numerical constant (Ormel & Liu 2018), q_p 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.

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.

7.5 Increasing the fragmentation threshold for icy particles

Up to now, we have fixed the fragmentation threshold velocity for icy particles at v_frag,icy = 10 m s⁻¹. 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⁻¹ to 60 m s⁻¹. 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 10³ yr, after which the critical pebble mass flux decreases. At around 6 × 10³ 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).

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.

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 CO₂ (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 ~ 10⁵ 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 r_out, 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 H₂O 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 r₁ = 1 au if r < r_out 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 v_dr (the expression for which we copy from Youdin & Shu 2002) for a constant particle size: v_dr ∝ r^d. For our temperature profile T ∝ r^−0.5 and gas surface density profile Σ_gas ∝ r^−3∕2 we find d = 3∕2. r_i (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.

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.

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

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 D_gas = 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 < D_gas. 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 r_snow, surface density to Σ_snow and velocities to D_gas∕r_snow: (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 H₂O vapor at the snowline r_snow – the pebble drift velocity is normalized by D_gas∕r_snow and the radius r is normalized by r_snow. 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 + b_gas = 3∕2 (because , the total mass flux in ice, is directed inwards) in an α-disk. Applying these, we obtain (B.2)

where x = r∕r_snow − 1 as in SO17 and a_eq is a constant that enters the expression for the equilibrium vapor density of H₂ O beyond the snowline. Getting rid of ~ notation and putting D_gas = r_snow = 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

References

All Tables

All Figures

	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)).
In the text

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.

In the text

	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.
In the text

	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).
In the text

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.

In the text

	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).
In the text

	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).
In the text

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.

In the text

	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 .
In the text

	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.
In the text

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.

In the text

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.

In the text

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.

In the text

	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.
In the text

	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).
In the text