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/00046361/201834047  
Published online  10 December 2018 
A Lagrangian model for dust evolution in protoplanetary disks: formation of wet and dry planetesimals at different stellar masses
^{1}
Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Science Park 904,
1090 GE Amsterdam,
The Netherlands
email: d.schoonenberg@uva.nl
^{2}
Department of the Geophysical Sciences, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637, USA
Received:
8
August
2018
Accepted:
3
October
2018
We introduce a new Lagrangian smoothparticle method to model the growth and drift of pebbles in protoplanetary disks. The Lagrangian nature of the model makes it especially suited to following characteristics of individual (groups of) particles, such as their composition. In this work we focus on the water content of solid particles. Planetesimal formation via streaming instability is taken into account, partly based on previous results on streaming instability outside the water snowline that were presented in a recent publication. We validated our model by reproducing earlier results from the literature and apply our model to steadystate viscous gas disks (with constant gas accretion rate) around stars with different masses. We also present various other models where we explore the effects of pebble accretion, the fragmentation velocity threshold, the global metallicity of the disk, and a timedependent gas accretion rate. We find that planetesimals preferentially form in a local annulus outside the water snowline, at early times in the lifetime of the disk (≲10^{5} yr), when the pebble mass fluxes are high enough to trigger the streaming instability. During this first phase in the planet formation process, the snowline location hardly changes due to slow viscous evolution, and we conclude that assuming a constant gas accretion rate is justified in this first stage. The efficiency of converting the solids reservoir of the disk to planetesimals depends on the location of the water snowline. Cooler disks with a closerin water snowline are more efficient at producing planetesimals than hotter disks where the water snowline is located further away from the star. Therefore, lowmass stars tend to form planetesimals more efficiently, but any correlation may be overshadowed by the spread in disk properties.
Key words: planets and satellites: formation / protoplanetary disks / methods: numerical
© ESO 2018
1 Introduction
Rocky planets and the cores of gas giants form from micrometersized 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 fullsized planets is the formation of approximately kilometersized planetesimals, which mark the transition to a gravitationdominated 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 centimetersized 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 solidstogas 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 solidstogas 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 TRAPPIST1 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 selfconsistently 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 superparticles 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 meshfree: physical quantities are computed at the locations of the superparticles, 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 lowmass stars are efficient at forming superEarths (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, lowmass stars convert a larger fraction of the solids reservoir in their disks to planetesimals than highmass 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 steadystate gas surface density profile Σ_{gas} for a given gas accretion rate Ṁ_{gas} and dimensionless turbulence parameter α (Shakura & Sunyaev 1973; LyndenBell & 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.
Input parameters, their symbols, and their fiducial values.
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 selfconsistently. 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 smoothparticle 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 superparticles: groups of particles with the same physical properties. In this work we consider two different classes of superparticles: dust/pebble superparticles^{1}, and planetesimal superparticles (in future work, more superparticle classes may be added to the model). Pebble superparticles can move through the disk, whereas planetesimal superparticles are inert. We assume a monodisperse particle size distribution at each point in space and time, such that a pebble superparticle represents particles that all have the same individual particle mass.
The characteristics of the pebble superparticles that we follow are location, total mass, individual particle mass, water mass fraction.
The planetesimal superparticles 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 superparticles to the planetesimal superparticles.
Formally, we end up with a system of ordinary differential equations (ODEs) that govern the time evolution of the superparticles, of the following form: (8)
where i indicates the superparticle number and j corresponds to their properties: radial position r, individual particle mass m_{p} (for pebble superparticles only), total superparticle mass M. The righthand 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 superparticles X. We solve Eq. (8) making use of the Python package scipy.integrate.ode.
3.1 Smoothparticle method
In Krijt et al. (2016), the solids surface density and its radial derivative at the location of each superparticle is calculated using a “tripod” method: each superparticle 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 superparticles evolve independentof each other. In contrast, in our method all superparticles are connected. As in SPH methods (e.g. Laibe et al. 2008), we approximate the solids surface density at each superparticle location using a weighting kernel W, that accounts for the contribution of neighbouring superparticles 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 onedimensional; 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 kernelapproximated 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 superparticle x′ and the location of the superparticleofinterest 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 superparticle locations), we turn the integral in Eq. (9) into a sum over all simulated superparticles. For the surface density Σ at the location x_{i} of superparticle i, we then find (11)
where M_{j} is the mass of supporting superparticle j and h_{i} is the smoothing length of superparticle i. We treat the smoothing length h_{i} as a variable, demanding that at each superparticle location and time h_{i} takes on a value such that there are five neighbouring superparticles (including superparticle 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 superparticle we only have to sum over the superparticles 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 cutoff 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 cutoff 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 superparticle follows from the initial solids surface density profile Σ_{solids}(r) and the initial number of superparticles 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 superparticle 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^{−1} (Σ_{solids}(r) = Σ_{0}r^{−1} where Σ_{0} is constant) and the spacing between the annuli is linear, each superparticle has the same initial mass m = 2πΔrΣ_{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 superparticle masses differ. In this work we have used a logarithmic spacing between the initial positions of the dust/pebble superparticles, so that superparticles in the outer disk have larger total masses than superparticles in the inner disk.
3.3 Resampling
During the evolution of the disk, the spacing between adjacent superparticles 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 superparticles are reset to the initial configuration, and the characteristics of the new superparticle population are sampled from the superparticle population right before the resampling process. The total mass of each new superparticle 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 superparticles becomes toolarge, the number and configuration of particles is reset such that the resolution is again equal to that of the initial setup. The total mass in pebble superparticles 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 superparticles 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 superparticle all have the same size s_{p} and mass m_{p}. Within each superparticle then, particles can coagulate, increasing the individual particle mass m_{p} of that superparticle. 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 wellcoupled 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 meanfree 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 meanfree path of the gas molecules, and a particle description is required instead (Weidenschilling 1977). The meanfree path l_{mfp} is given by (15)
where ρ_{gas} is the gas volume density at the disk midplane () and σ_{mol} is the molecular collisional crosssection, for which we take σ_{mol} = 2 × 10^{−15} cm^{2} 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^{−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 superparticle 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 socalled 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^{−1} (Blum & Münch 1993); and for icy particles, we take v_{frag,ice} = 10 m s^{−1} (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 semimajor 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 driftlimited region. In this region, drift is faster than growth/fragmentation (Birnstiel et al. 2012). In the viscositydominated inner disk, m_{frag} decreases for inwarddrifting particles (see Fig. 3). In this region therefore, inwarddrifting 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 viscositydominated region is fragmentationlimited, 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 fragmentationlimited regime, treating fragmentation as an instantaneous process is justified throughout the disk.
Fig. 3 Fragmentation mass limit m_{frag} (black line) for our fiducial model with initial ice mass fraction profile f_{ice} (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 backreaction 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 solidstogasvolume 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 wellknown 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^{13} 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 steadystate advectiononly vapor surface density profile from SO17^{2}. We note that Σ_{Z,a} is not the physical water vapour surface density profile, but rather the steadystate 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^{3}. 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 semimajor 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 r_{snow} ~ 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 semianalytic 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 wellcoupled 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 socalled 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 solidstogas 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 solidstogas 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 solidstogas 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 solidstogas 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 solidstogas ratio can form interior to the snowline due to a trafficjam 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 recondensation 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 superparticles containing planetesimals with the same physical characteristics. Planetesimal superparticles are characterized by a position, a total mass and a ice mass fraction. In contrast to pebble superparticles, the location of planetesimal superparticles is fixed; they are treated as sink particles. At the beginning of the simulation, a number N_{pltsml} planetesimal superparticles are initiated at a given location, with zero mass. If the solidstogas ratio at a certain radius exceeds unity (either directly or inferred from the “snowline recipe”), mass is transferred from the pebble superparticles to the nearest planetesimal superparticle.
4.1 Direct streaming instability
Interior to the snowline, the solidstogas ratio can exceed unity “directly”, due to a trafficjam 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 superparticles 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 solidstogas 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 semianalytic model presented in SO17 to calculate Ṁ_{peb,crit} for different stopping times^{4}. We assume the “manyseeds” 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 inwarddrifting 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 solidstogas 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 semianalytic 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 solidstogas ratio is continuously sustained, and only the excess incoming pebble mass flux is transformed to planetesimals^{5}. The mass of the planetesimal superparticle closest to the peak location r_{peak} (which is given by the semianalytic model) is increased at a rate dM_{pl,SI,snowline}∕dt.
We distribute the pebble mass loss rate over all pebble superparticles using a Gaussian kernel centered on the peak location r_{peak}. The mass loss rate of pebble superparticle i is then given by (32)
where the fraction f(r_{i}) < 1 depends on the distance between r_{i} and the peak location^{6}.
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 solidstogas 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 superparticles that start in the icy outer disk and eventually convert partof their mass to planetesimals. The first superparticle (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 superparticle 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 superparticle does not lose any mass to planetesimals outside the snowline. As the superparticle crosses the snowline, its H_{2} 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 fragmentationlimited, 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 solidstogas ratio with decreasing distance to the star. In this cartoon the solidstogas ratio reaches unity close to the star, resulting in dry planetesimal formation (Eq. (30)).
The second superparticle (indicated by a “2” in the cartoon) starts out further away from the star, also representing 0.1 micrometersized icy particles. These particles grow to pebblesizes at a slower pace than the particles of superparticle 1, because the growth timescale (Eq. (13)) is longer for larger distances from the star (because the solids density goes down with increasing semimajor axis). At the time when the superparticle 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 superparticle 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 superparticles (1 and 2). Both superparticles start out in the outer disk with a small typical particle mass. The physical particles represented by the superparticles 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 solarmass 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 semimajor 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^{5} 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 fragmentationlimited and the driftlimited 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 pileup 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 pileup 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 fragmentationdominated phase to the driftdominated phase outside the snowline. For t < 10^{5} yr, the region just outside the snowline is still fragmentationdominated (solids surface density ∝ r^{−3∕2}; Birnstiel et al. 2012). At t = 2 × 10^{5} yr, the region outside the snowline has become driftdominated, which is shown by an r^{−1} 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^{3} 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 solidstogas 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^{4} years. The total yield in planetesimals is about ten Earth masses.
Fig. 6 Particle mass m_{p} as function of semimajor 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 × 10^{5} 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 superparticles Ṁ_{pltsmls} in Earth masses at the last time point plotted (t = 5 × 10^{5} yr) (this is also the final total mass in planetesimals, because the planetesimal formation phase ends at ~ 3.2 × 10^{4} yr). 
Fig. 8 Pebble mass flux Ṁ_{peb} that arrives at the snowline r_{snow} 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 r_{snow}. 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 M_{pltsmls} 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 steadystate 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^{−3}. We find that models 1a, 2b, and 3a (corresponding to 1, 0.1, and 0.5 solarmass stars, respectively) convert similar percentages of their initial solids reservoir to planetesimals (~ 10%). Model 4a (M_{⋆} = 2M_{⊙}) does not form any planetesimals. The critical pebblestogas 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 pebblestogas 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}α^{−1}, and neglecting effects such as the timeevolution 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 powerlaw index of the pebbles surface density. In case the region outside the snowline is fragmentationlimited, as is the case in our simulations (except for model 1l; see Sect. 7.5); p = −3∕2. In the driftlimitedcase 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 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 driftlimited and for the fragmentationlimited case. In both cases, the pebblestogas 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 setup of the αviscosity gas disks (Σ_{gas} ∝Ṁ_{gas}α^{−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 pebblestoplanetesimal 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 zerothorder 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 viscouslyrelaxed gas profile; the time when the fragmentationlimited region outside the snowline becomes driftlimited, etcetera, could play important roles as well. For example, one might notice that model 3a has a higher pebblestoplanetesimals conversion efficiency than model 2b, which is not expected from the snowlinedependency 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 pebblestogas mass flux ratio outside the snowline depends on the snowline location, we cannot directly translate this dependency to a simple scaling law for the pebblestoplanetesimals 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 Timedependent 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_{1} is the characteristic radius of thedisk and ν_{1} is the viscosity at r_{1} (LyndenBell & 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^{−2}, 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^{−1} (see also Sect. 7.5), to get nonnegligible drift velocities. The timedependent accretion rate Ṁ_{gas}(t) is given by (37)
where we take for the characteristic radius of the disk r_{1} = 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 powerlaw index as the gas profile. Outside the snowline, particles are fragmentationlimited, 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 driftlimited and the solids profile evolves toward an r^{−1} profile. If we would have chosen a smaller value for α in this timedependent 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 buildup 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 semimajor 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 timedependent 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 solidstogas 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 solidstogas 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 pileup 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^{−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.
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 superparticles are plotted in blue; the masses of the waterpoor planetesimal superparticles are plotted in brown. A massive annulus of planetesimals forms outside the snowline. Interior to the snowline, a trafficjam effect leads to a high solidstogasratio 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 superparticle “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 threedimensional 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 threedimensional pebble accretion rate is given by (Liu & Ormel 2018): (38)
where A_{3} = 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 superparticles due to pebble accretion is the negative of Eq. (38) and the individual pebble superparticle 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 selflimiting 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 selfexcite 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 typeI 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 10^{4} 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^{−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 10^{3} yr, after which the critical pebble mass flux decreases. At around 6 × 10^{3} 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 driftlimited case (Fig. 13) than in the fragmentationlimited case (Fig. 8).
Fig. 13 Results for model 1l, which is the fiducial model but with icy fragmentation threshold v_{frag,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 lowermass stars are more efficient at forming planetesimals than highermass stars, that is, disks around lowmass 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 pebblestogas 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 highermass stars that have hotter disks. The pebblestoplanetesimals 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 snowlinespecific 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 gasphase abundances of CO and CO_{2} (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^{5} yr) when the pebble mass flux is still high. A side effect of fast planetesimal formation is the wellknown “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 monodisperse 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 earlystage 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 smoothparticle 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 Nbody 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 Nbody 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_{2}O fraction of the TRAPPIST1 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 selflimiting 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 closerin the water snowline), the higher the pebblestoplanetesimals conversion efficiency. Therefore, in general, we find that lowmass stars are better at producing planetesimals than highmass 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 steadystate α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 HSTHF251394 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS526555. 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 smoothparticle model of dust growth and radial drift against analytical results from the literature.
A.1 Driftonly
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} = 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 pebblesizes 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 razorsharp 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 efoldings needed to grow to pebblesizes), was used regardless of semimajor 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 Updatedsemianalytic 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 semianalytic 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 rederive 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 lefthand side represents the mass flux of the ice; ~ denote nondimensional quantities. The ice surface density (Σ) is normalized by the quantity Σ_{snow} – the surface density in H_{2}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 righthand 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_{2} 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
Parameter values and results for different model runs.
References
 Akiyama, E., Muto, T., Kusakabe, N., et al. 2015, ApJ, 802, L17 [NASA ADS] [CrossRef] [Google Scholar]
 ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Benz, W. 2000, Space Sci. Rev., 92, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Münch, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994 [NASA ADS] [CrossRef] [Google Scholar]
 Bosman, A. D., Tielens, A. G. G. M., & van Dishoeck, E. F. 2018, A&A, 611, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S., & Cowling, T. 1970, The Mathematical Theory of Nonuniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge Mathematical Library (Cambridge: Cambridge University Press) [Google Scholar]
 Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [Google Scholar]
 Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Dra̧żkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn. (Cambridge: Cambridge University Press), 398 [Google Scholar]
 Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [Google Scholar]
 Gonzalez, J.F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Hallis, L. J. 2017, Phil. Trans. R. Soc. London, Ser. A, 375, 20150390 [NASA ADS] [CrossRef] [Google Scholar]
 Hallis, L. J., Huss, G. R., Nagashima, K., et al. 2015, Science, 350, 795 [NASA ADS] [CrossRef] [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S., & Guillot, T. 2016, A&A, 596, L3 [Google Scholar]
 Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low, M.M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J.,& Hartmann, L. 1987, ApJ, 323, 714 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, H., Ormel, C. W., & Ida, S. 2012, ApJ, 756, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Hicks, D. L., & Liebrock, L. 2000, Appl. Math. Comput., 112, 63 [Google Scholar]
 Laibe, G., Gonzalez, J.F., Fouchet, L., & Maddison, S. T. 2008, A&A, 487, 265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322 [CrossRef] [Google Scholar]
 Lichtenegger, H. I. M., & Komle, N. I. 1991, Icarus, 90, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, G., & Liu, M. 2003, Smoothed Particle Hydrodynamics: A Meshfree Particle Method (Singapore: World Scientific) [Google Scholar]
 Liu, B., & Ormel, C. W. 2018, A&A, 615, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, B., Ormel, C. W., & Johansen, A. 2018, A&A, submitted [Google Scholar]
 Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Lorek, S., Gundlach, B., Lacerda, P., & Blum, J. 2016, A&A, 587, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martin, R. G., & Livio, M. 2012, MNRAS, 425, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Mulders, G. D. 2018, Handbook of Exoplanets, eds. H. Deeg & J. Belmonte (Cham: Springer) [Google Scholar]
 Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2016, ApJ, 819, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W. 2017, Astrophys. Space Sci. Lib., 445, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Liu, B. 2018, A&A, 615, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Ros, K., & Johansen, A. 2013, A&A, 552, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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. F677) [Google Scholar]
 Saito, E., & Sirono, S.i. 2011, ApJ, 728, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg,D., Okuzumi, S., & Ormel, C. W. 2017, A&A, 605, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Sirono, S. 1999, A&A, 347, 720 [NASA ADS] [Google Scholar]
 Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
 Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Visser, R. G., & Ormel, C. W. 2016, A&A, 586, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (Slona: Almqvist och Wiksell), 211 [Google Scholar]
 Yang, L., Ciesla, F. J., & Alexander, C. M. O. 2013, Icarus, 226, 256 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
 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]
We use the term “dust” for solid particles that are wellcoupled to the gas, and the term “pebbles” for particles that have a nonnegligible radial drift velocity (Stokes numbers of ~ 10^{−3}−10^{1}; see Sect. 3.4). Throughout the paper, the term “pebble superparticle” is used instead of “dust/pebble superparticle”.
If we take the accretion rate Ṁ_{gas} to be a decreasing function of time and the snowline is located in the viscositydominated 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.
The semianalytic 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.
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 inwarddrifting pebbles before streaming instability is triggered again, leading to episodic bursts of planetesimal formation. If the width of the peak W_{peak} 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.
To be specific, Eq. (32) reads: (33)
where N_{peb} is the number of pebble superparticles and W_{Gauss} is given by (34)
The term [(1 − f_{ice,out}) + f_{ice}] in Eq. (33) takes into account the fact that in our model the location of the solids enhancement outside the snowline r_{peak} can be at a distance from the star where f_{ice} < f_{out}, again because we do not follow the transport of water vapour explicitly in this work.
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 superparticles becomes toolarge, the number and configuration of particles is reset such that the resolution is again equal to that of the initial setup. The total mass in pebble superparticles 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 m_{frag} (black line) for our fiducial model with initial ice mass fraction profile f_{ice} (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 semimajor 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 r_{snow} ~ 2.4 au (blue dotted line). 

In the text 
Fig. 5 Cartoon of two possible trajectories through space and time of pebble superparticles (1 and 2). Both superparticles start out in the outer disk with a small typical particle mass. The physical particles represented by the superparticles 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 m_{p} as function of semimajor 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 × 10^{5} 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 superparticles Ṁ_{pltsmls} in Earth masses at the last time point plotted (t = 5 × 10^{5} yr) (this is also the final total mass in planetesimals, because the planetesimal formation phase ends at ~ 3.2 × 10^{4} yr). 

In the text 
Fig. 8 Pebble mass flux Ṁ_{peb} that arrives at the snowline r_{snow} 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 r_{snow}. 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 M_{pltsmls} is plotted inblue. 

In the text 
Fig. 9 Temperature T as function of semimajor 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 timedependent 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 superparticles are plotted in blue; the masses of the waterpoor planetesimal superparticles are plotted in brown. A massive annulus of planetesimals forms outside the snowline. Interior to the snowline, a trafficjam effect leads to a high solidstogasratio 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 10^{4} 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 v_{frag,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 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.