Planetesimal formation starts at the snow line

Planetesimal formation starts at the snow line

Key Words.:
accretion, accretion disks – stars: circumstellar matter – protoplanetary disks – planet and satellites: formation – methods: numerical


Context:Planetesimal formation stage represents a major gap in our understanding of the planet formation process. The late-stage planet accretion models typically make arbitrary assumptions about planetesimals and pebbles distribution while the dust evolution models predict that planetesimal formation is only possible at some orbital distances.

Aims:We want to test the importance of water snow line for triggering formation of the first planetesimals during the gas-rich phase of protoplanetary disk, when cores of giant planets have to form.

Methods:We connect prescriptions for gas disk evolution, dust growth and fragmentation, water ice evaporation and recondensation, as well as transport of both solids and water vapor, and planetesimal formation via streaming instability into a single, one-dimensional model for protoplanetary disk evolution.

Results:We find that processes taking place around the snow line facilitate planetesimal formation in two ways. First, due to the change of sticking properties between wet and dry aggregates, there is a ”traffic jam” inside of the snow line that slows down the fall of solids onto the star. Second, ice evaporation and outward diffusion of water followed by its recondensation increases the abundance of icy pebbles that trigger planetesimal formation via streaming instability just outside of the snow line.

Conclusions:Planetesimal formation is hindered by growth barriers and radial drift and thus requires particular conditions to take place. Snow line is a favorable location where planetesimal formation is possible for a wide range of conditions, but still not in every protoplanetary disk model. This process is particularly promoted in large, cool disks with low intrinsic turbulence and increased initial dust-to-gas ratio.

1 Introduction

Our understanding of planet formation is severely limited by the fact that we cannot explain the connection between its early and late stages. As a consequence, models that deal with the late-stage planet accretion, when a planetary embryo grows to its final size and structure, typically use the same input for the radial distribution of gas and solids as the early-stage models dealing with dust growth and planetesimal formation. However, the latter models show that growing large bodies is not easy. This is because of the growth barriers: the dust growth is inhibited at centimeter sizes and some particular conditions are needed for the formation of larger, gravitationally bound planetesimals and planetary embryos.

Probably the most widely accepted planetesimal formation scenario at the moment is the streaming instability (Youdin & Goodman, 2005; Johansen et al., 2007). For sufficiently large pebbles and increased metallicity (Bai & Stone, 2010; Dra̧żkowska & Dullemond, 2014; Carrera et al., 2015), streaming instability leads to formation of dense filaments that become gravitationally unstable and collapse to planetesimals. This scenario allows us to bypass the growth barriers and form gravitationally bound object directly from pebbles.

The streaming instability is typically simulated in local boxes due to the high computational cost of hydrodynamical simulations (Johansen et al., 2011; Kowalik et al., 2013; Simon et al., 2016). Thus, the initial conditions are already set up in a way that streaming instability occurs. However these conditions do not necessarily happen in a realistic disk that starts its evolution with dust-to-gas ratio on the order of 1%, which gets depleted because of removal of solids by the radial drift (Birnstiel et al., 2010; Hughes & Armitage, 2012; Krijt et al., 2016b). Pebble pile-ups may be necessary to allow for planet formation in the gas-rich phase of protoplanetary disk (Dra̧żkowska et al., 2016; Gonzalez et al., 2017), while disk dispersal via photoevaporation may allow for late planetesimal formation (Carrera et al., 2017). Both processes may be needed to explain the existence of different planet types and debris belts, for instance in the Solar System. In this paper, we focus on the former mechanism, with the aim of triggering planetesimal formation early in the evolution of gas disk in order to allow sufficient time for the formation of gas-rich planets.

The great dichotomy of the Solar System used to be commonly attributed to the jump change of condition around the snow line (Stevenson & Lunine, 1988; Wuchterl et al., 2000; Morbidelli et al., 2015), with water ice greatly enhancing the abundance of solids outside of this point. Such a rapid change of conditions may contribute to a pressure bump build-up that would halt the radial drift of solids thus facilitating planet formation (Kretke & Lin, 2007; Brauer et al., 2008; Dra̧żkowska et al., 2013). It was also proposed that the icy dust aggregates can pile-up (Cuzzi & Zahnle, 2004) or even significantly grow thanks to water vapor recondensation at the snow line (Ros & Johansen, 2013; Wang, 2015). Laboratory and numerical experiments dealing with collisional properties of dust aggregates concluded that the icy aggregates are significantly more sticky than silicate grains (Wada et al., 2011; Gundlach et al., 2011) and can thus grow to larger sizes before they fragment, or even grow directly to planetesimal sizes if they are sufficiently porous (Okuzumi et al., 2012; Kataoka et al., 2013).

In this paper, we analyze how the presence of the snow line could trigger formation of the first gravitationally bound planetesimals at the very beginning of planet formation. We start our simulations with a smooth protoplanetary disk and let it evolve taking into account dust growth to pebbles, their drift and fragmentation, as well as ice evaporation and recondensation. In order to demonstrate the universality of our findings, we test our scenario in three different protoplanetary disk models. We conclude that the water component has an immense effect on the growth and redistribution of solids and leads to a pile-up of icy pebbles and planetesimal formation via streaming instability just outside of the snow line. Significant contribution to this pile-up comes from the change of sticking properties between icy and dry aggregates, an effect that was previously included in some of the models (Birnstiel et al., 2010; Banzatti et al., 2015; Estrada et al., 2016; Cridland et al., 2017) but not discussed directly in the context of planetesimal formation.

This paper is organized as follows. We describe our numerical modeling approach and typical initial conditions in Sect. 2 and present typical results and their dependence on input parameters in Sect. 3. We discuss the differences between our work and other published results as well as the implications of our findings in Sect. 4, and finally summarize our work in Sect. 5.

2 Model

We implement one-dimensional protoplanetary disk model where we follow the radial distribution of solids and formation of planetesimals over one million years. We focus on the gas-rich phase of protoplanetary disk, before photoevaporation is efficient, thus we either include only viscous evolution or implement a static gas disk. The initial dust grains size is 1 m at every orbital distance. We assume that the dust may be composed of water ice and silicates, and that the initial ice mass fraction outside of the snow line is 50%. At the beginning of each simulation, the ice and water vapor are distributed across the disk such that the vertically integrated water-to-gas ratio is uniform and equal to 0.5 of the total metallicity. Solid ice is only present outside of the snow line, and water vapor is present inside of the snow line. The refractory dust component is present both outside and inside of the snow line. We follow dust growth to pebbles, fragmentation, and radial drift as well as ice evaporation and recondensation of water vapor.

Solids and water vapor evolution is governed by their interactions with the sub-Keplerian, turbulent gas. We treat the gas disk as a background for solids evolution and plug-in different gas disk models, as described in the following section.

2.1 Gas disk models

To test versatility of the planetesimal formation scenario discussed in this paper, we will use three different protoplanetary disk models. All of them have total gas mass of  M within 100 AU distance to the central star of  M. Their other properties are described in subsequent paragraphs. Comparison of the basic properties of these disks: surface density, temperature, and deviation from the Keplerian rotation, is displayed in Fig. 1.

Figure 1: Comparison of the initial conditions for our protoplanetary disk models: the power-law disk (black), the non-irradiated disk (blue), and the irradiated disk (red). The panels show, from top to bottom, gas surface density , temperature , and the difference between gas and Keplerian rotation, which is equal to the maximum drift speed of dust pebbles . Each of the disks has the total mass of 0.1 .

Power-law disk

Very simple protoplanetary disk models are commonly used in planet formation research. Here, we adopt one of them, with gas surface density profile set as the straightforward function of the distance to the central star


and temperature profile is fixed to


For simplicity, as the model is already basic, we keep this disk static, however we include effects of gas accretion velocity when calculating radial drift of solids. This velocity is estimated as


where turbulent gas diffusivity is equal to gas viscosity


which is calculated based on the standard -accretion model (Shakura & Sunyaev, 1973) with being the dimensionless parameter describing the efficiency of angular momentum transfer. We typically use in this paper. The is the Keplerian orbital frequency and the gas sound speed is calculated as


where is the Boltzmann constant and is the mean molecular weight of gas. When calculating , we include the contribution of water vapor:


which increases the molecular weight and thus decreases the sound speed in the inner part of the protoplanetary disk. We take  m and  m, with m denoting the proton mass. The effect of variable is visible in the bottom panel of Fig. 1 as a small jump of the maximum drift speed (equal to the difference between the gas and Keplerian rotation) around the snow line location. The maximum radial drift speed is roughly constant in the power-law disk model, with  m s.

Non-irradiated disk

We compare the simple static power-law disk model to a more complex, viscously evolving circumstellar disk model described in Alibert et al. (2005, 2013), which was also used in our previous work presented in Dra̧żkowska et al. (2016). In this model, the vertical structure of the disk is first computed for every distance to the star by solving the hydrostatic equation, the energy conservation equation, and the radiative diffusion equation (energy is assumed to be transported by radiation). The solution of the vertical structure equations gives the thermodynamical quantities (temperature, pressure, density, but also disk scale height) as a function of the gas surface density. The same calculation also gives the vertically averaged viscosity, which is again computed in the framework of the formalism. We finally solve the diffusion equation, the viscosity being the one derived from the vertical structure calculation, in order to compute the time evolution of the gas surface density.

As can be seen in the middle panel of Fig. 1, this disk is very cold in its outer part, as the temperature drops down to 10 K outside of 10 AU. Inside of the 10 AU, the temperature profile is steeper than the one implemented in the power-law disk and reaches over 1000 K at the inner edge of the disk. This steeper temperature profile leads to a flatter surface density in the inner part of the disk, while it is very similar to the power-law disk in the outer part. In the inner part of the disk, where the temperature profile is steep, the maximum radial drift speed is higher than for the power-law disk, but it drops in the outer disk, forming a wide minimum around 10 AU (see the bottom panel of Fig. 1).

Irradiated disk

To test the effect of stellar irradiation, we implement a simple analytical model proposed by Bitsch et al. (2015) (B15 in the following), which was designed to fit 2-D radiative hydrodynamic simulations of protoplanetary disks. In this model the disk evolves with time and the accretion rate decreases as follows (Hartmann et al., 1998)


where is the evolution time and is related to the viscosity and the gas surface density via


where the scale height of the disk, and the Keplerian frequency. In this disk model, we use following B15. We note that is only used in B15 as a heating parameter and not to evolve the disk viscously as it is in the case of the non-irradiated disk model. For a given , all quantities of the disk can be derived, except for the temperature profile. In order to compute the latter, we use the formulas presented in Appendix A of B15.

Once the temperature is determined, it can be linked to the disk aspect ratio via


where is the gravitational constant, the location in the disk, is the gas constant and is the mean molecular weight.

As can be seen in Fig. 1, the temperature in this disk is significantly higher than in the non-irradiated disk, particularly in its outer part, thus the surface density flattening is also more pronounced. Also, the maximum drift speed is higher than in the non-irradiated disk, but still lower than the power-law disk in the outer parts.

2.2 Evolution of solids

We follow the solids surface density evolution by solving the advection-diffusion equation


where is the mass weighed average radial velocity of solids, which encompasses information about their size.

Dust aggregates sizes are modeled using method based on the two-population algorithm proposed by Birnstiel et al. (2012). The logic behind this algorithm is to reduce computational intensity by not solving for the dust coagulation directly but rather predict its outcome based on more complete models. Thus, the dust size distribution is set in each radial grid cell depending on the dominating process: coagulation-fragmentation equilibrium or radial drift. We waive the description of further details of this model as they are thoroughly outlined in the original work as well as in our previous paper (Dra̧żkowska et al., 2016). Importantly, the outcome of this algorithm is the size, or Stokes number, distribution at every radial distance. The Stokes number informs us about the interaction between the solids grain and gas, and is connected to the grain size :


where is the aggregate internal density. Eq. (11) is derived under an assumption that the solids are in the Epstein drag regime, which means that their sizes do not exceed the mean free path in the gas, which is true for all our models. The internal density of aggregates is calculated based on their composition:


where and are mass of silicate and mass of ice contained within the aggregate. A pure water ice aggregate would have  g cm and a silicate aggregate  g cm.

When calculating the advection speed of solids , we take into account both radial drift caused by the interaction with the sub-Keplerian gas and by gas accretion flow:


where is the maximum radial drift speed driven by the radial gas pressure gradient:


where is the Keplerian velocity of gas and is pressure calculated taking into account contributions both from nebular hydrogen and helium gas and water vapor. is the mass weighted average Stokes number of solids at given radial distance, which is calculated from the size distribution returned by the above-mentioned algorithm. The gas accretion velocity is calculated as in Eq. (3). We take into account the collective drift effect, which means that the drift velocity decreases as the solids-to-gas ratio increases. As most of the pebbles are settled to the midplane, we implement the midplane solids-to-gas ratio in Eq. (13). This equation is equivalent to the one used by Ida & Guillot (2016) and Schoonenberg & Ormel (2017).

Fragmentation threshold

Evolution of solids is dominated by the radial drift and by fragmentation that may stop the growth as the impact speeds increase with the Stokes number of grains (for ). The maximum aggregate size that can be obtained before fragmentation kicks in is sensitively dependent on the fragmentation threshold velocity (see Birnstiel et al., 2012):


so the choice of value is in fact very important to the model outcome. is the midplane turbulence strength parameter that regulates impact speeds of pebbles and their settling. We purposely distinguish from , the efficiency of angular momentum transport via turbulent viscosity that is used in the gas disk models. This is motivated by the fact that in many recent protoplanetary disk models the angular momentum transfer is not necessarily driven by turbulence anymore, and even if it is, a quiescent midplane layer is often formed (Dzyurkevich et al., 2013; Turner et al., 2014; Bai, 2016). In most of our runs we assume , but we discuss the impact of this value in Sect. 3.3.1.

Laboratory experiments estimated threshold velocities of around 1 m s for the onset of fragmentation of silicate dust aggregates (see e.g. Güttler et al., 2010). It is commonly accepted that aggregates containing water ice fragment at higher velocities, as their surface energies are about 10 times higher than those of silicates (Wada et al., 2011; Gundlach et al., 2011; Aumatell & Wurm, 2014; Gundlach & Blum, 2015). Thus, we set the fragmentation threshold velocity according to aggregates composition. For dry aggregates we set  m s and for aggregates containing more than 1% of water ice  m s. This threshold amount of ice above which we consider our aggregates to be more sticky is arbitrary, but we tested that the exact value does not make much difference to the results of our models as the ice-to-dust ratio drops very rapidly across the snow line.

Evaporation and recondensation

We take water ice evaporation and recondensation of water vapor into account following the treatment proposed by Ciesla & Cuzzi (2006). We trace the transport of solid ice that is incorporated into aggregates and thus migrates through the disk much faster than the gas, and the water vapor that moves at the same speed as the gas. To follow the evolution of water vapor surface density , we solve the following transport equation:


where is the gas velocity. Eq. (16) is analogous to Eq. (10) used to follow the transport of solids.

We assume that all the dust grains at a given orbital distance (i.e. in a given radial bin) have the same composition (i.e. ice mass ratio), which means that water is added to and removed from the dust with constant during recondensation and evaporation. This is consistent with instantaneous redistribution of ice component because of coagulation and fragmentation, which happens when the collisional timescale


(see Birnstiel et al., 2012) is shorter than the radial drift timescale


In our runs, this is indeed true around the snow line, where nominal timescales are  years and  years. When the pebble pile-up is formed around the snow line, the vertically integrated dust-to-gas ratio increases making the growth timescale even shorter, and the radial drift speed decreases (because of the collective drift effect) making the drift timescale even longer.

At every timestep, we calculate the equilibrium pressure, which is given by the Clausius-Clapeyron equation:


where the constants  g cm s and  K are taken from Lichtenegger & Komle (1991). We compare the value of to the water vapor pressure


where is the water vapor surface density and is the gas scale-height.

If , evaporation takes place and the surface density of ice decreases by


where and are the average size and mass of pebbles aggregates respectively, and is the time-step. The material removed from solid ice phase is added to the vapor reservoir.

If , the water vapor condenses onto grains and the surface density of ice is increased by


which essentially means that all the excess vapor is added to the solid phase, such that the vapor pressure drops to the equilibrium pressure. The surface density added to the ice is subsequently removed from the vapor supply.

In the original work of Ciesla & Cuzzi (2006) there was a distinction between populations of dust and migrators (pebbles). Evaporation was considered to happen both from dust and migrators, while condensation was only happening on dust, thus it is assumed to be instantaneous. In our model, we do not explicitly make this distinction between dust and pebbles, but we assume that there is a continuous size distribution in every radial cell, which around the snow line is set by coagulation-fragmentation equilibrium. However, we do not model the evaporation and condensation on every size bin but take into account the surface-weighted average size . Condensation on grains larger than micron-sized should take some time because they have smaller surface area available. In reality however, vapor is mostly condensing on the smallest grains that have the most surface area available (see e.g. Stammler et al., 2017). Since the snow line region is in the fragmentation dominated regime (see Fig. 4), the small grains should be constantly replenished, so we keep the assumption of instantaneous recondensation. In practice, the same assumption might have been made for evaporation, as it is very fast (a centimeter sized pebble crossing the snow line would lose its ice content within  year) in all of the runs presented in this paper.

Planetesimal formation via streaming instability

With our one-dimensional model, we cannot resolve streaming instability that would locally condense pebbles into dense filaments, which would then gradually collapse to form planetesimals. To include the possibility of planetesimal formation via this process, we use the same approach as we did in Dra̧żkowska et al. (2016). At every time-step and in every radial bin, we check if the midplane density of pebbles exceeds unity. With the turbulence level parameter that we use in this paper, this condition is always stronger than the critical metallicity conditions proposed for laminar disk by Dra̧żkowska & Dullemond (2014) and Carrera et al. (2015) (recently Carrera et al. (2017) arrived at the same conclusion).

If the criterion for planetesimal formation via streaming instability is fulfilled, namely , part of the surface density of pebbles is transferred to planetesimals:


with the efficiency of , that is motivated by numerical models presented by Simon et al. (2016) and for which Dra̧żkowska et al. (2016) found convergence of the amount of planetesimal formed, i.e. the amount of planetesimal would not change significantly for higher values.

2.3 Model assumptions

Our algorithm is limited by several assumptions, which we list here for clarity.

  • We consider one-dimensional, locally isothermal disk models and focus on the evolution of their midplane, where pebbles and planetesimals reside. Thus, we only consider the radial snow line and neglect effects associated with the atmospheric snow line.

  • Grain sizes are set by either the coagulation-fragmentation or the growth-drift equilibrium. We do not consider the impact of evaporation and condensation on aggregate sizes. This is equivalent with assuming that the ice added during recondensation is quickly redistributed by coagulation and fragmentation. Any aggregates that would increase their size over the maximum are immediately fragmented or removed by radial drift.

  • During one time-step and at a given orbital distance, all the aggregates are either in evaporation or in condensation mode. We do not consider grain curvature effects that could switch between those effects from grain to grain.

  • We treat all grains as compact spheres and neglect effects of porosity.

  • We assume that vertical structure of solids is always in equilibrium between settling and turbulent mixing, which leads to the dust scale-height derived by Dubrulle et al. (1995)


    and that the water vapor is instantly mixed up to gas scale-height , even though it is released by pebbles with . Recently, Krijt et al. (2016a) showed that the Eq. (24) breaks at high dust-to-gas ratios, when the collisional evolution is faster than the vertical redistribution.

  • We assume that minimum pebble size necessary to trigger planetesimal formation via streaming instability corresponds to following Bai & Stone (2010) and Dra̧żkowska & Dullemond (2014). More recently, Carrera et al. (2015) and Yang et al. (2016) suggested that streaming instability is also possible for smaller grains, albeit at higher metallicity. However, we find that for our assumed turbulence strength, the smaller grains do not settle the the midplane efficiently (see Eq. (24)), so their contribution to planetesimal formation would not be possible anyway in our models.

  • The structure of our gas disks is independent on solids evolution. This looses validity for high dust-to-gas ratio, when , which happens when the pile-up forms in our models. Recently, Gonzalez et al. (2017) and Kanagawa et al. (2017) showed that including effects of dust backreaction on gas disk promotes formation and sustainability of dust pile-ups.

3 Results

3.1 Traffic jam inside and pebble pile-up outside of the snow line

Figure 2: This sketch explains processes facilitating formation of the snow line pile-up: a) In the initial condition, ice presence increases solids density outside of the snow line; b) Coagulation is more efficient for aggregates that incorporate water ice. Thus, solids grow to larger sizes and drift faster outside of the snow line. The quick drift results in an efficient delivery of the embedded refractory material, which do not drift rapidly, causing ”traffic jam” and increase of dust concentration in the inner disk; c) The outward diffusion and recondensation of water vapor locally enhances abundance of solids just outside of the snow line, contributing to the pile-up of icy pebbles.

The common outcome of all the runs, independently of the underlying gas disk model, is that the highest solids-to-gas ratios are obtained in the region directly outside of the snow line. A general pattern for formation of this snow line pile-up is explained in the Fig. 2.

The dust-to-gas ratio in the inner parts of the disk is enhanced because of the general pattern of dust transport by radial drift that shifts mass inwards. As explained by Birnstiel et al. (2012), if the maximum size of dust grains is regulated by fragmentation, the surface density of solids becomes proportional to . As the gas surface density is shallower, this redistribution leads to depletion of the outer disk and increase of solids-to-gas ratio in the inner disk. As demonstrated by Dra̧żkowska et al. (2016), the magnitude of this increase may already be sufficient to trigger planetesimal formation at the inner edge of the disk at a timescale of  years. However, in the present work this picture is complicated by the difference in fragmentation speeds of aggregates outside and inside of the snow line.

Refractory aggregates fragment at lower impact velocities and thus reach sizes that are two orders of magnitude smaller than the icy aggregates (see Eq. (15) and Fig. 4). As the drift velocity decreases with decreasing Stokes number, dust is retained in the inner disk causing sort of ”traffic jam” effect. The dry dust aggregates inside of the snow line are small and well-coupled to the gas, so they undergo significant diffusion and do not form any local pile-up, unlike the icy pebbles outside of the snow line. The outward diffusion of water vapor and subsequent recondensation causes further increase in surface density of the icy pebbles just outside of the snow line. These pebbles are large enough to trigger streaming instability and form planetesimals. The exact location, extent, and mass of resulting planetesimal annulus depends on applied disk model, as discussed in Sect. 3.3.

3.2 Fiducial simulation

Figure 3: Panels show time evolution of surface density of gas, solids, planetesimals, water ice, and water vapor in the irradiated disk model with initial dust-to-gas ratio of . In the bottom panel, surface density of solids in the classical minimum mass solar nebula model (MMSN) of Weidenschilling (1977) is displayed for reference. Corresponding animation is available at̃oannad/snowline.mp4.

Figure 3 shows the evolution of the gas and dust surface density, including ice and water vapor as well as the forming planetesimals in the irradiated disk model with initial dust-to-gas ratio of , which we will refer to as the fiducial simulation later in the paper.

As can be seen in the uppermost panel of Fig. 3, presence of the solid ice increases the initial surface density in the outer part of the disk. The predominant effect shaping the evolution of solids is their redistribution driven by growth and radial drift, which shifts mass inwards causing depletion of the outer parts and solids-to-gas ratio increase in the inner parts of the disk. Initially, the evolution outside of the snow line is dominated by fragmentation of the icy aggregates and thus the dust surface density is evolving toward , as discussed in previous section. After  years of evolution, as the outer disk gets depleted, it becomes dominated by the radial drift and the surface density profile at the outer edge becomes more shallow again, as it evolves to (Birnstiel et al., 2012).

This picture is complicated here by ice evaporation and recondensation and different sticking properties of wet and dry aggregates. Keeping the grains inside of the snow line small, slows down their removal and keeps the dust-to-gas ratio in the inner disk high. The surface density of dust inside of the snow line keeps the same profile as the gas surface density for most of the time, because the aggregates are so small that their evolution is dominated by gas viscosity rather than radial drift (see Fig. 4). Ice evaporation leads to a jump in the surface density at the snow line. Recondensation of water vapor additionally bumps pebbles-to-gas ratio just outside of the snow line. The combined action of the the ”traffic jam” inside and pile-up outside of the snow line, which spreads outwards thanks to the collective drift effect (as the drift velocity decreases with increasing solids-to-gas ratio, see Eq. (13)), leads to the conditions allowing for planetesimal formation via streaming instability.

The planetesimals start to appear after  years of evolution and their formation last for another  years. During this time, the snow line, which marks the inner edge of planetesimal formation region, moves inwards as the disk cools down. At the same time, the planetesimal formation region slightly spreads outwards because of the collective drift effect. About 20 M of planetesimals are produced in this model. After  years, the inward flux of pebbles is not sufficient to supply the pile-up anymore and the surface density of solids quickly drops, terminating the planetesimal formation phase.

Pebble sizes

Figure 4: Maximum pebble size () as a function of the radial distance after  years of evolution obtained in our fiducial run. In the outer part of the disk, this size is limited by the radial drift (), while in the inner part of the disk it is limited by fragmentation. Aggregates that are large enough to participate in planetesimal formation via streaming instability (, purple solid line) are only present outside of the snow line. Inside of the snow line, aggregates size corresponds to (gray solid line), which means their transport is dominated by viscosity rather than radial drift.

Figure 4 shows pebble sizes obtained from our simplified growth and fragmentation treatment after  years of evolution (during the period of planetesimal formation) in our fiducial simulation. These effective sizes are very similar in other runs.

The dust growth follows pattern described by Birnstiel et al. (2012) and Dra̧żkowska et al. (2016), with the inner disk being dominated by fragmentation driven by turbulence and the outer disk being gradually depleted by the radial drift before the particles have time to grow to the fragmentation limit. The modification of fragmentation velocity described in Sect. 2.2.1 introduces the rapid change of pebble size around the snow line. The icy pebbles outside of the snow line grow to sizes of several centimeters, corresponding to Stokes numbers of , which allows for planetesimal formation. The dry aggregates inside of the snow line only grow to sub-millimeter sizes, corresponding to , which means that they are well-coupled to the gas and follow its viscous evolution, or are in so-called mixing regime, as discussed by Birnstiel et al. (2012). Thus, these small grains drift at much lower speed than the icy pebbles, which contributes to enhancement of dust-to-gas ratio inside of the snow line and the retention of icy pebble pile-up outside of the snow line thanks to the outward diffusion of small grains.

Which is more important: traffic jam or recondensation?

Figure 5: Vertically integrated solids-to-gas ratio around the snow line after  yrs (dotted lines) and  yrs (solid lines) of evolution for our fiducial run (black lines) and analogous runs without the water recondensation (purple lines) and without the difference between sticking properties of wet and dry aggregates (gray lines).

An obvious question arises: which of the two processes facilitating planetesimal formation at the snow line is more important? To answer this dilemma, we performed additional models with setup identical as in the fiducial run, but with one of the mechanisms switched off. Fig. 5 compares solids-to-gas ratio around the snow line obtained in our fiducial run, the same run without recondensation of water vapor (but still including ice evaporation), and analogous run with fragmentation velocity for both wet and dry aggregates being equal to  m s, designed to get rid of the ”traffic jam” effect inside of the snow line.

As can be seen in Fig. 5, the difference between sticking properties of wet and dry aggregates is the dominant process facilitating formation of the solids-to-gas ratio enhancement at the snow line at every evolutionary stage of the disk. Switching the recondensation off decreases the amplitude of the solids-to-gas ratio enhancement at the snow line only by . On the other hand, letting the dry aggregates grow to the same sizes as the icy pebbles leads to much more dramatic decrease of the solids-to-gas ratio in the bump. If the dry aggregates grow to pebble sizes, the solids-to-gas ratio falls more rapidly across the snow line and stays at lower levels inside of it. The surface density outside of the snow line reaches lower values and the pile-up vanishes more quickly in this case.

3.3 Planetesimal formation

Disk1 Z R2 3 R4 R5
0.01 100
0.02 100 3.54 1.98 2.34
P-L 0.03 100 24.28 1.93 2.99
0.04 100 62.76 1.89 3.90
0.05 100 115.10 1.87 4.90
0.03 100 217.13 1.66 5.19
P-L 0.03 100 128.79 1.81 3.44
0.03 100
0.03 100
0.03 60 8.90 1.89 2.45
P-L 0.03 80 16.73 1.91 2.75
0.03 140 38.52 1.94 3.50
0.03 200 49.50 1.98 4.10
0.01 100 119.80 1.03 2.20
0.02 100 371.37 0.98 4.08
N-IRR 0.03 100 663.18 0.96 5.50
0.04 100 928.35 0.96 6.83
0.05 100 1130.9 0.95 8.00
0.01 100
0.02 100 3.32 2.12 4.33
IRR 0.03 100 22.50 1.64 7.16
0.04 100 63.31 1.25 10.37
0.05 100 129.93 1.06 16.99
Table 1: Details on the planetesimal annuli: their final mass and radial extent (R and R), formed in runs with different underlying protoplanetary disk models, global solids-to-gas ratio , turbulence strength , and initial disk size R.

Planetesimal formation is triggered in the pile-up arising outside of the snow line in many of the runs that we performed. Table 1 summarizes the information about mass and extent of the planetesimals annulus formed in models with different underlying gas parameters, initial dust-to-gas ratio, intrinsic turbulence level and disk extent. We notice that it is significantly easier to trigger planetesimal formation when the snow line is closer to the central star. The initial metallicity of is sufficient for planetesimal formation in the non-irradiated disk model, which is the coldest (see the middle panel of Fig. 1), while the power-law and the irradiated disks need at least . This is because close-in pile-up formation is aided by the process described by Dra̧żkowska et al. (2016): the surface density of solids in the fragmentation dominated regime naturally evolves to , a profile steeper than the gas disk, which leads to additional enhancement of dust-to-gas ratio in the inner disk. The further away the border between fragmentation dominated and drift dominated regions (see Fig. 4) is, the stronger enhancement we obtain. The maximum enhancement formed by this process would fall at the inner edge of the disk (although it is diffused in our models because the small aggregates inside of the snow line falling into the mixing regime). For this reason, the annuli formed closer to the central star tend to be more massive.

The inner edge of the planetesimal formation zone depends primarily on the underlying disk model, and slightly changes with the assumed metallicity. This is because higher abundance of solids translates into higher flux of pebbles coming to the snow line, which deliver more water vapor and thus increase the vapor pressure that determines the snow line location (see Sect. 2.2.2). At the same time, the higher the initial metallicity, the wider planetesimal annulus we get. The broadening of planetesimal annulus with higher incoming pebble flux is caused by the collective drift effect. The radial drift slows down near the solids-to-gas ratio peak and the more pebbles are delivered to this region, the wider this peak becomes.

We stress that the collective drift effect is a critical component of our model without which obtaining the significant pile-up and planetesimal formation is nearly impossible. Recently, Schoonenberg & Ormel (2017) arrived at the same conclusion.

Impact of the turbulence strength

Figure 6: Planetesimals annuli obtained in the power-law disk models with different turbulence strength parameter (upper panel) and initial disk extent (bottom panel). Black solid line corresponds to the same model in both panels.

All the models presented in this paper up to this point were performed assuming that the fragmentation and settling of solids is regulated by turbulence with . However this value is rather vague, as it is very challenging to estimate the turbulence strength from the observational data. Recent estimates range from for the outer parts of the disk around HD 163296 (Flaherty et al., 2015) to in the outer parts of the TW Hya disk (Teague et al., 2016).

To test the impact of turbulence strength, we perform a suite of models where we vary the parameter value and keep all the other parameters constant. For this purpose, we use the setup with static, power-law disk and initial dust-to-gas ratio of , and vary the between and . As can be seen in the upper panel of Fig. 6, lower values of the parameters facilitate planetesimal formation. The lower the turbulence strength, the wider and more massive the resulting planetesimal annulus becomes. We find that no planetesimal formation is possible for significantly higher than our fiducial (see Table 1). This is because the higher value reduces size of pebbles that can grow (see Eq. (15)) and decreases the possibility of their settling (see Eq. (24)), both of the factors counteract the possibility of obtaining conditions necessary to trigger the streaming instability, namely the minimum size of pebbles corresponding to and the midplane pebbles-to-gas ratio exceeding unity (see Sect. 2.2.3).

Impact of the initial disk size

The initial size of protoplanetary disk is rather uncertain. Observational constraints place the outer edge of the disk anywhere between 60 AU and several hundreds AU (Andrews et al., 2009, 2010). Thus, we decided to test the impact that initial disk extent has on the planetesimal formation.

The bottom panel of Fig. 6 shows the impact that initial distribution of material has. We used the power-law disk model with the initial dust-to-gas ratio of and turbulence strength of again and tested how the results change when the extent of this disk is different from fiducial 100 AU (however keeping the total mass of the disk constant). The variation in the resulting planetesimal formation is not quite as pronounced as when varying the turbulence strength, but the larger the disk, the more massive and more extended the planetesimals annulus becomes. This is because larger disks provide long lasting supply of pebbles, as their growth takes longer at larger orbital distances. In other words, increasing the disk size (while keeping its mass unchanged) shifts more solid mass to its outer regions, and this reservoir can be then used to form more planetesimals.

Pebbles and planetesimals composition

Figure 7: Panels show time evolution of the ice fraction of planetesimals (red) and pebbles (gray) for the irradiated disk model with initial solids-to-gas ratio of 0.03. The initial ice mass fraction (outside of the snow line) of 0.5 is marked with the dashed line.

Figure 7 presents the time evolution of ice fraction of pebbles (gray line) and ice fraction of the resulting planetesimals (red line) for the irradiated disk model with initial metallicity of 0.03 (the same model as presented on Figs. 3 and 4). The first panel is plotted at the beginning of planetesimal formation stage and the last panel corresponds to a time shortly after planetesimal formation is terminated. We assume that the composition of planetesimals reflects those of pebbles from which they are forming. We start all our models with dust outside of the snow line consisting of 50% water ice (dashed line). Evaporation removes the solid ice that is delivered to the inner disk by radial drift and turns it to vapor, part of which is able to diffuse outwards and recondense, moderately enhancing the ice content of dust aggregates just outside of the snow line. As explained in Sect. 3.2.2, the main cause for the pile-up of pebbles outside of the snow line is the ”traffic jam” arising in the inner disk and not recondensation, which is reflected here in the low amplitude of the ice fraction enhancement. When planetesimal formation starts after about  years of evolution, the enhancement of pebbles ice fraction goes away completely. This is because the streaming instability turns the icy pebbles from outside of the snow line to planetesimals and thus hinders delivery of water to the evaporation region, and decreases the rate of recondensation. During the planetesimal formation phase, the snow line moves inwards as the disk evolves and thus the ice content of the final planetesimal population decreases more smoothly than those of pebbles, which are evaporating rather quickly. After the planetesimal formation is over ( years of evolution), pebbles composition is set again by the interplay of evaporation and recondensation, causing the mild enhancement of ice fraction outside of the evaporation region.

4 Discussion

4.1 Comparison to published work

Armitage et al. (2016) presented similar idea where pebbles drift radially and pile-up to reach conditions required for planetesimal formation. However, they neglected dust growth and fragmentation, assuming constant pebble size throughout the disk. In such a setup, the Stokes number increases with the radial distance (because the gas surface density drops, see Eq. (11)), making it easier to trigger planetesimal formation further away from the star. At the same time, the particles further away in the disk drift faster (because the radial drift speed depends on the Stokes number, not size, see Eq. (13)), which leads to a pile-up that is harder to generate when particle sizes are decided by fragmentation and radial drift. As derived by Birnstiel et al. (2012), in the fragmentation-dominated regime the steady-state dust surface density is proportional to and in the radial drift regime to . Analogous derivation assuming constant dust size gives , which facilitates obtaining high dust-to-gas ratio in the inner disk more than our models.

Recently, Schoonenberg & Ormel (2017) performed local models focusing on the snow line region. They took into account the inflow of icy pebbles that would be formed in the outer disk and brought to the snow line region by radial drift and the outflow of water vapor carried with gas accretion. They found that water diffusion and recondensation enhances surface density of icy pebbles by a factor of 3-5 outside the snow line that is further increased if the evaporating pebbles release many small refractory ”seeds” that help the pile-up thanks to the ”traffic jam” effect and outward diffusion. With the models presented in this paper, although they are fundamentally different by construction, as Schoonenberg & Ormel (2017) focused on a local box and used particle approach, while in this paper we perform global disk models applying fluid approach to dust dynamics, we find very similar both qualitative and quantitative results. This is very encouraging and proving that a pile-up of pebbles outside the snow line is a robust mechanism that can trigger formation of the first, icy planetesimals.

On the other hand, another recent work by Ida & Guillot (2016) suggested different scenario of planetesimal formation, namely pile-up of dry dust grains released by the icy pebbles inside of the snow line. They find that for a sufficiently high flux of pebbles, the dust-to-gas ratio increases to a point when direct gravitational instability is possible. We do also find an increase in the solids-to-gas ratio inside of the snow line but it is not as significant. The main difference is that we assume that the small grains released from the icy pebbles quickly coagulate until they reach fragmentation limit at about millimeter sizes and that they are vertically mixed by the turbulence, while Ida & Guillot (2016) consider that the grains stay at micron sized and that their scale-height is equal to the scale-height of icy pebbles that released them.

Both Ida & Guillot (2016) and Schoonenberg & Ormel (2017) stressed the importance of sufficiently high pebble mass flux for the possibility of planetesimal formation. In their models, the pebble flux was a free parameter since they did not include dust growth to pebble sizes and their drift self-consistently. With our models, we can measure the pebble flux incoming to the snow line region. For the power-law disk model, which is similar to the models used in the quoted papers, we measure the ratio of pebble mass flux to the gas mass flux during the planetesimal formation stage on the order of 0.5, which Schoonenberg & Ormel (2017) also found sufficient for triggering planetesimal formation via streaming instability outside of the snow line.

4.2 Implications for planet formation

Our results suggest that formation of planetesimals is simpler in the direct vicinity of the water evaporation front. This could naturally explain the fast formation of Jupiter in the Solar System (Kruijer et al., 2017). The surface density obtained in the planetesimal annulus in all of our runs is higher than predicted by the minimum mass solar nebula models (Weidenschilling, 1977; Hayashi, 1981). Starting from a high mass concentration translates into faster growth of planetesimals to planetary embryos and the planetary cores. It is known that enhancement of about 10 times over the solids surface density corresponding to the minimum mass solar nebula is necessary to allow for Jupiter core formation before the gas disk dispersal (Pollack et al., 1996; Ikoma et al., 2000; Kobayashi et al., 2010). The fast formation of a gas giant outside of the snow line could possibly halt the delivery of water to the inner part of planetary system (Morbidelli et al., 2015). If such a barrier is not formed quickly, it may be problematic to explain the low water content of terrestrial planets in the Solar System (Sato et al., 2016).

Our results, including the self-consistent surface densities and composition of planetesimals and pebbles, may be used as an input to models dealing with later stages of planet accretion, particularly discussing the pebble accretion process, when the planetary cores grow by accreting not only planetesimals, but also the leftover pebbles that were not incorporated during the planetesimal formation stage. The radial dependence of sizes and radial flux of pebbles that we can extract from our results are important parameters of the pebble accretion models (Ormel & Klahr, 2010; Lambrechts & Johansen, 2014; Levison et al., 2015; Visser & Ormel, 2016).

One prominent consequence of the planetesimal formation mechanism we discuss is that the first planetesimals are water-rich (see Fig. 7). As a consequence, planetary cores formed from these planetesimals would also be water-rich. For small mass planets (without a massive gas envelope), the presence of large amounts of water may be detrimental for habitability (see Alibert et al. (2013); Kitzmann et al. (2015), see however Levi et al. (2017) for another view). Planetesimal formation mechanism we describe here could therefore imply that the majority of low mass planets are not habitable. However, the recent models show that low mass, short period planets detected by the Kepler mission should be water-poor (see e.g. Jin & Mordasini, 2017). In the framework of our planetesimal formation model, this imply that other mechanisms, facilitating efficient water loss from existing planetesimals or allowing the formation of dry planetesimals, are at work to prevent the accumulation of water on these short-period planets.

5 Summary

In this paper, we addressed the problem of connection between dust evolution and planetesimal formation, which is still one of the least certain aspects of the planet formation theory. As dust growth is hindered by collisional fragmentation and radial drift, continuous growth from micron to planetesimal sizes appears to be improbable. We propose that the first planetesimals form via streaming instability, in a pile-up of icy pebbles generated outside of the snow line.

The water snow line is a favorable location for planetesimal formation as large icy pebbles efficiently deliver water and embedded refractory material to the inner part of the disk. The water vapor is partially mixed outwards by diffusion and recondenses just outside of the snow line, locally enhancing the solids-to-gas ratio. At the same time, the less sticky, dry aggregates inside of the snow line drift at much lower speed, creating the ”traffic jam” effect and helping to reach high dust concentration which is needed to form planetesimals in the streaming instability scenario.

A relatively compact annulus of icy planetesimal is a common result of our simulations, performed with three diverse protoplanetary disk models and different input parameters. The main condition we find for making this outcome feasible is that the turbulence strength cannot be too high: the corresponding parameter must stay equal to or below . Also, the further away the snow line is located, the higher disk metallicity is needed to allow for planetesimal formation.

On a more general level, this work as well as similar studies (Dra̧żkowska et al., 2016; Carrera et al., 2017; Gonzalez et al., 2017) indicate that dust distribution during and after the planetesimal formation stage is very different from the commonly assumed power-laws, as significant redistribution of solids must take place before the conditions necessary for planetesimal formation happen. In particular, the solids-to-gas ratio is significantly increased at the location where planetesimals form, which should facilitate faster accretion of the final planets.

We thank the anonymous referee for their detailed reports that helped us to improve this paper. This work has been carried out within the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF.


  1. P-L – power law, N-IRR – non-irradiated, IRR - irradiated,
  2. in AU,
  3. in Earth masses
  4. in AU,
  5. in AU,
  6. .


  1. Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109
  2. Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343
  3. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502
  4. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241
  5. Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2
  6. Aumatell, G. & Wurm, G. 2014, MNRAS, 437, 690
  7. Bai, X.-N. 2016, ApJ, 821, 80
  8. Bai, X.-N. & Stone, J. M. 2010, ApJ, 722, 1437
  9. Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15
  10. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79
  11. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  12. Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28
  13. Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1
  14. Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16
  15. Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  16. Ciesla, F. J. & Cuzzi, J. N. 2006, Icarus, 181, 178
  17. Cridland, A. J., Pudritz, R. E., & Birnstiel, T. 2017, MNRAS, 465, 3865
  18. Cuzzi, J. N. & Zahnle, K. J. 2004, ApJ, 614, 490
  19. Dra̧żkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105
  20. Dra̧żkowska, J. & Dullemond, C. P. 2014, A&A, 572, A78
  21. Dra̧żkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37
  22. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  23. Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114
  24. Estrada, P. R., Cuzzi, J. N., & Morgan, D. A. 2016, ApJ, 818, 200
  25. Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99
  26. Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984
  27. Gundlach, B. & Blum, J. 2015, ApJ, 798, 34
  28. Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717
  29. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  30. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
  31. Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  32. Hughes, A. L. H. & Armitage, P. J. 2012, MNRAS, 423, 389
  33. Ida, S. & Guillot, T. 2016, A&A, 596, L3
  34. Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013
  35. Jin, S. & Mordasini, C. 2017, ArXiv e-prints [\eprint[arXiv]1706.00251]
  36. Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62
  37. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
  38. Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142
  39. Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4
  40. Kitzmann, D., Alibert, Y., Godolt, M., et al. 2015, MNRAS, 452, 3752
  41. Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836
  42. Kowalik, K., Hanasz, M., Wóltański, D., & Gawryszczak, A. 2013, MNRAS, 434, 1460
  43. Kretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55
  44. Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016a, ApJ, 833, 285
  45. Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016b, A&A, 586, A20
  46. Kruijer, T. S., Kleine, T., Burkhardt, C., & Budde, G. 2017, in Lunar and Planetary Inst. Technical Report, Vol. 48, Lunar and Planetary Science Conference, 1386
  47. Lambrechts, M. & Johansen, A. 2014, A&A, 572, A107
  48. Levi, A., Sasselov, D., & Podolak, M. 2017, ApJ, 838, 24
  49. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322
  50. Lichtenegger, H. I. M. & Komle, N. I. 1991, Icarus, 90, 319
  51. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418
  52. Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
  53. Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43
  54. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  55. Ros, K. & Johansen, A. 2013, A&A, 552, A137
  56. Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15
  57. Schoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21
  58. Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  59. Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55
  60. Stammler, S. M., Birnstiel, T., Panić, O., Dullemond, C. P., & Dominik, C. 2017, A&A, 600, A140
  61. Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146
  62. Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49
  63. Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
  64. Visser, R. G. & Ormel, C. W. 2016, A&A, 586, A66
  65. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36
  66. Wang, X.-M. 2015, MNRAS, 449, 1084
  67. Weidenschilling, S. J. 1977, Ap&SS, 51, 153
  68. Wuchterl, G., Guillot, T., & Lissauer, J. J. 2000, Protostars and Planets IV, 1081
  69. Yang, C.-C., Johansen, A., & Carrera, D. 2016, ArXiv e-prints [\eprint[arXiv]1611.07014]
  70. Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.