© S. Charnoz et al. 2019

1 Introduction

The current formation paradigm for planetesimals asserts that they form through streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Bai & Stone 2010). Alternative possibilities also exist, such as gravitational instability just inward of the snow line (Ida & Guillot 2016) or direct growth of very porous aggregates (Arakawa & Nakamoto 2016 and Tatsuuma et al. 2018). However, the conditions for the onset of streaming instability are quite stringent, as they require that most dust particles be close to pebble size (Stokes number S t > 10⁻²) and that the dust-to-gas ratio in the midplane be higher than 1.

Drążkowska et al. (2016), Drążkowska & Alibert (2017), and Drążkowska & Dullemond (2018) have studied this problem using a 1D alpha disk model of an active disk in order to simultaneously track the dust growth and planetesimal formation using a parametric prescription physically motivated by local simulations of streaming instability (Johansen et al. 2007; Bai & Stone 2010). It is found that planetesimal formation in an ice-free disk (i.e., all dust made only of silicates) is very sensitive to the fragmentation velocity considered for dust growth (U_f, which is the dust impact velocity beyond which fragmentation occurs, rather than sticking). For silicate dust U_f ≃ 1 m s⁻¹ (Blum & Wurm 2008) is commonly considered, whereas for water ice U_f ≃ 10 m s⁻¹ (Gundlach & Blum 2015). When water ice is taken into account, it is found that planetesimals form preferentially at the water snow line because of the retro-diffusion of water that recondenses just outward of the snow line, leading to a local enhancement of the solid-to-gas ratio and back-reaction of the gas onto the dust (Drążkowska & Alibert 2017). Schoonenberg & Ormel (2017) found a similar result using a simple and detailed local model. However, as this model is local it does not consider large-scale transport of the dust, and both studies focused on fully active disks, i.e., where the turbulence intensity (quantified by α) is constant everywhere.

It is now believed that the midplane of a protoplanetary disk is not turbulent, at least in the planet forming region (Bai & Stone 2010; Turner et al. 2014 and references therein). Turbulence induced by the magneto-rotational instability (MRI) requires the gas to be ionized. In the disk midplane the ohmic resistivity is high (Sano et al. 2000), supressing MRI, and in the low density regions, the ambipolar diffusion acts to further suppress MRI (Bai & Stone 2010; Dzyurkevich et al. 2013; see Turner et al. 2014 for a review). Recent works now favor disk structures where the midplane is mostly devoid of turbulence (called dead zones) and with a complex vertical structure: an actively accreting upper layer topped by a disk-wind carrying away angular momentum (Suzuki et al. 2010; Bai & Stone 2010; Fromang et al. 2013). In this context planetesimal formation and dust growth may be favored in the midplane, while gas accretion onto the star may be mostly driven in the disk upper layers. Following the above ideas, the present work aims to study planetesimal accretion and dust transport in a vertically stratified disk with a dead midplane and an active upper layer using the 1D approach of Zhu et al. (2010) and Hasegawa & Takeuchi (2015); this allows us to take into account two layers, coupled with the planetesimal formation prescription introduced by Drążkowska et al. 2016.

One original aspect of this work is to simultaneously compute the dynamical and thermal evolution of the disk so that our planetesimal calculation takes place in a time evolving disk, and is not limited to a steady state. Even though steady state approaches provide very important insights into the physics, transient processes must also be studied. In this respect,this work differs from those of Schoonenberg & Ormel (2017) and Drążkowska & Alibert (2017), which only consider a steady state, or where the disk thermal evolution is either fixed or decoupled from the vapor transport. In addition, we consider different species with different condensation temperatures that are of importance for planet formation: highly refractory minerals (like calcium-aluminum-rich inclusions, CAIs), silicate and iron (relevant for terrestrial planets and the core of giant planets), and water ice and CO ice (relevant for the outer solar system). We consider grains whose composition reflects the local solid material composition of the disk. As a consequence grains have a composition that varies continuously with distance to the star. Handling such grains requires a specific numerical method in order to properly account for mass and chemical composition exchanges.

We would like to address the following questions, relevant to disks with a dead zone:

• Where do planetesimals form in a disk with a dead zone?

• What are the key physical parameters?

• What is the composition of planetesimals after their formation?

• What is the pebble-to-planetesimal mass ratio as a function of space and time?

We focus on an already formed protoplanetary disk, i.e., an isolated disk whose initial structure is arbitrarily imposed at the beginning of the simulation, in the spirit of the popular minimum mass solar nebula. We do not consider a forming disk fed by a collapsing molecular cloud, though this is the subject of future work. Drążkowska & Dullemond (2018) recently investigated this question, but without considering the presence of a dead zone.

This paper is organized as follows. In Sect. 2 we describe our method. In Sect. 3 we present the time evolution of disks containing a dead zone and explore two parameters: the disk metallicity and the fragmentation velocity. In Sect. 5 we discuss our results and draw conclusions.

2 Method

The 1D disk model we present here is not really new; it is an extension of the Hueso & Guillot (2005) disk model, with some additional physical modules that take into account multi-species, planetesimal formation, dust growth, dust sublimation, and layering due to the DZ. We detail the different aspects below.

2.1 System

The protoplanetary disk is made of two components: gas and dust. The gas is made of different components: H and He (the dominant ones), H₂O, Si, Fe, refractory species (e.g., aluminium or calcium, designated as REF hereafter), and CO. The solid component of the disk is made of Si-rich species, Fe metal, REF, H2O, and CO. Each species is characterized by a condensation temperature (25, 150, 1500, 1550, 1650 K, for CO, H₂O, Si, Fe, and REF, respectively). Condensation temperatures depend on the local pressure; however, this effect is ignored here for simplicity and to reveal the basic mechanisms. For dust, we assume that all grains are composite and have locally the same composition and the same size distribution, and that they are a mix of the element species described above. As such, the modeling does not give us the possibility to have pure H₂O grains at the same location co-existing with pure CO grains, for example. The central star has a mass M_* that increases with time because of gas accretion. It has an effective temperature of 4000 K and a radius of 3 solar radii, constant over the simulation. The disk is described using a 1D logarithmic grid extending from 0.02 to 500 au with 300 radial cells.

2.2 Transport of gaseous species

Gaseous species are transported in the viscous disk using a classical Sakura-Sunaev α prescription (see, e.g., Hueso & Guillot 2005; Yang & Ciesla 2012; Baillié & Charnoz 2014 and references therein). We recall here the main aspects. We call the total gas surface density σ_g and the local gas velocity V_g. Using the continuity equation, we obtain the equation of gas surface density evolution $$\frac{\partial {\sigma }_{\text{g}}}{\partial t}=\frac{-1}{r}\frac{\partial }{\partial r}(r{\sigma }_{\text{g}}{V}_{\text{g}}),$$(1)

where r is the distance to the star. The term in parentheses on the right-hand side is the mass flux. Assuming that the gas velocity is given by the viscous stress, we have (Yang & Ciesla 2012) $$V_{\text{g}}=\frac{-3r^{-1/2}}{{\sigma }_{\text{g}}}\frac{\partial }{\partial r}(\nu {\sigma }_{\text{g}}{r}^{1/2}),$$(2)

where ν is the local viscosity that is averaged vertically over the disk (because of the presence of a dead zone, the value of ν in the midplane may be different from ν in the upper layers; see below). Minor gaseous species (H-He, REF, SI, Fe, H₂O, CO) are followed through Eq. (3), taking into consideration advection with the gas, and diffusion inside the gas due to turbulence. We call σ_g,i the surface density of a gaseous species i (which may be H-He, H₂O, SI, Fe, REF, or CO) and verifying mass conservation, $${\sigma }_{\text{g}}={\sigma }_{\text{g},\text{H-He}}+{\sigma }_{\text{g},\text{H2O}}+{\sigma }_{\text{g},\text{SI}}+{\sigma }_{\text{g},\text{CO}}+{\sigma }_{\text{g},\text{Fe}}+{\sigma }_{\text{g},\text{REF}},$$(3)

the value of σ_g,i is evolved as $$\frac{\partial {\sigma }_{\text{g},i}}{\partial t}=\frac{-1}{r}\frac{\partial }{\partial r}(r{\sigma }_{\text{g},i}{V}_{\text{g}})+\frac{-1}{r}\frac{\partial }{\partial r}(r{\sigma }_{\text{g}}{V}_{\text{g},i}).$$(4)

Here V_g,i stands for the radial velocity due to turbulent diffusion of gaseous species i. The first term in parentheses on the right-hand side is the advective flux, and the second term is the diffusive flux. Following Yang & Ciesla (2012) we have: $$V_{\text{g},i}=-D_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\sigma }_{\text{g},i}}{{\sigma }_{\text{g}}}\right),$$(5)

with D_g standing for the diffusivity of the gas, taken to be equal to the disk viscosity (Hueso & Guillot 2005 or Fromang & Papaloizou 2006).

2.3 Transport of solid species

The main originality of the present approach is the consideration of composite grains, i.e., grains with a non-homogeneous chemical composition. However, this makes the transport algorithm more complex than usual transport codes because grains composition varies with distance; this means that the chemical composition of the dust must be properly tracked during the system evolution. In particular, self-diffusion will be an important process. Self-diffusion is the process by which, even if the net diffusive flux between two adjacent cells is 0, grains are still exchanged between these two cells because of Brownian motion or turbulence (the net diffusive flux is the difference of the fluxes crossing a frontier in opposite directions). If all the grains were the same size with a constant composition, self-diffusion would have no effect. However, here grain composition and size varies from cell to cell. As a consequence, even in the case of zero net mass exchange between two adjacent cells, Brownian motion and/or turbulent motion induces an exchange of composition between neighboring cells. So it is necessary to explicitlycompute the mass exchange between all neighboring cells in both directions. To know exactly the amount of material transported between neighboring cells, we need to know explicitly the mass flux leaving and entering each cell, on both edges. We detail the transport algorithm, with explicit calculations of material entering and leaving each cell in Appendix B. Below we just recall the governing equations for the dust evolution.

Let σ_d the local total dust surface density, and σ_d,i the differentmaterials surface density in dust form. We have by mass conservation: $${\sigma }_{\text{d}}={\sigma }_{\text{d},\text{H-He}}+{\sigma }_{\text{d},\text{H2O}}+{\sigma }_{\text{d},\text{SI}}+{\sigma }_{\text{d},\text{CO}}+{\sigma }_{\text{d},\text{Fe}}+{\sigma }_{\text{d},\text{REF}}.$$(6)

Diffusion induces a net transport of grains toward locations where dust concentration is smaller. The transport of solid grains is described by the classical advection diffusion dust transport equation (see, e.g., Yang & Ciesla 2012; Birnstiel et al. 2010 their Eq. (21)) $$\frac{\partial {\sigma }_{\text{d}}}{\partial t}=\frac{-1}{r}\frac{\partial }{\partial r}(r{\sigma }_{\text{d}}{V}_{\text{d}})+\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_{\text{d}}{\sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\sigma }_{\text{d}}}{{\sigma }_{\text{g}}}\right)\right).$$(7)

The first derivative on the right-hand side is the advective term, where V_d is the dust radial velocity, given by Eqs. (19) and (20) of Birnstiel et al. (2010). The second term is the diffusive flux and is the difference between the two diffusive fluxes leaving each neighboring cell. The parameter D_d is the diffusivity of the dust, that is $\alpha C_s^2/({\Omega }_{k}Sc)$ with α taken in the dead zone (as most dust resides in the dead zone rather than in the active layer), and Sc is the Schmidt number, which depends on the local Stokes number (Birnstiel et al. 2010). Equation (7) evolves the dust surface density. The dust radial velocity is (Brauer et al. 2008; Dra̧żkowska et al. 2013) $$V_{\text{d}}=\frac{2V_n}{S_t+1/S_t}+\frac{V_{\text{g}}}{1+S_t^2},$$(8)

where V_g is given by Eq. (2), while $$V_n=\frac{1}{2{\Omega }_k{\rho }_g}\frac{\text{d}P}{\text{d}r}$$(9)

(see, e.g., Dra̧żkowska et al. 2013). Now to compute explicitly the surface density of every chemical species transported by the dust, all fluxes of all species between adjacent cells must be known. This is detailed in Appendix B.

Fig. 1 Rosseland mean opacity as a function of the local temperature. The dust composition is given in Table 1.

2.4 Disk viscosity and dead zones

We consider a two-layer model, following the procedure of Zhu et al. (2010): a midplane layer that contains most of the disk mass and that can be turbulent or not (a dead zone), and an upper layer that is always turbulent. The disk viscosity is computed assuming an α prescription, so that $$\nu =\frac{\alpha C_{\text{s}}^2}{{\Omega }_k},$$(10)

where C_s is the local sound velocity (computed using the local gas mean molecular weight derived from the local gas composition; see Sect. 2.2) and Ω_k is the local Keplerian frequency. For a fully active disk, 10⁻³ < α < 10⁻² typically (see, e.g., Balbus & Hawley 1991), and the effective viscosity is defined as in Eq. (10). In the case of a DZ, α is in the range 10⁻⁵ –10⁻³ in the midplane, and could be >10⁻² in the upper layers. Studies of disk ionization (see, e.g., Turner & Sano 2008) show that in general only a column density < 100 Kg m⁻² is subject to ionization, and the remaining mass is not ionized when T < 1000 K (see, e.g., Terquem 2008; Zhu et al. 2010). When T > 1000 K the whole disk is thermally ionized. So in the presence of a dead zone, the effective viscosity is computed as $$\nu =\frac{{\sigma }_{\text{active}}{\nu }_{\text{active}}+{\sigma }_{\text{dead}}{\nu }_{\text{dead}}}{{\sigma }_{\text{g}}},$$(11)

where σ_active and ν_active are the surface density and viscosity of the active layer, and σ_dead and ν_dead are the surface density and viscosity in the dead zone. We note that σ_active + σ_dead = σ_g. Zhu et al. (2010) show that the disk advective flux can be computed by using the vertically averaged viscosity. However, for dust diffusion we take ν_dead as the gas diffusion coefficient because the vast majority of dust resides in the DZ.

Gravitational instability is also taken into account (but only plays a minor role in the present study) by increasing sharply the value of α at the location where the disk becomes gravitationally unstable (see Eqs. (6)–(8) in Zhu et al. 2010).

2.5 Disk temperature and opacity

The disk temperature, which includes viscous and radiative heating, is computed following Eqs. (22)–(25) in Hueso & Guillot (2005). The Rosseland mean opacity is taken from Baillié et al. (2016) and depends on the local temperature (see Fig. 1). Grains are assumed to follow an interstellar distribution (Mathis et al. 1977) and the composition is given in Table 1. It is assumed that the dust-to-gas ratio is 1%. Steep opacity transitions correspond to the sublimation of one of the components. Beyond 1500 K the values are extrapolated from the opacity tables of Bell & Lin (1994).

It should also be noted that the grain composition used in the opacity table is independent of the local dust grain composition in the disk and also independent on the local dust-to-gas ratio, inducing an inconsistency in our approach. Ciesla & Cuzzi (2006) introduced a way to compute the dust opacity self-consistently with the local dust abundance and dust-to-gas ratio. Following their pioneering work, we tried to run simulations modifying the opacity locally scaling with the local dust-to-gas ratio and water abundance, but we found the disk evolution to be very unstable numerically, triggering numerous oscillations. So more work is needed here to investigate this problem and to compute a local opacity that would be self-consistently calculated with chemistry. However, we think that there are already numerous physical processes with non-linear feedback in the present version of the study, and we will deal with these improvements in a future paper.

2.6 Dust sublimation and condensation

As previously noted, the sublimation temperature at which the different dust species (REF, Fe, Si, H₂O, CO) sublimate is fixed (1650, 1550, 1500, 150, 25 K, respectively) and is assumed to be independent of the pressure. This simplification was necessary to make the computation tractable and we do not consider a real chemical network. When a dust species i sublimates, then species i is put in the gas σ_g,i (to conserve mass) and σ_d,i is set to 0. The reverse procedure is used for condensation. Following Ida & Guillot (2016) and Schoonenberg & Ormel (2017) we assume that all grains crossing the snow line inward explode into micron-sized silicate dust grains, while the water ice content transforms into water vapor.

2.7 Dust density, coagulation, and fragmentation

In order to compute the effect of gas drag on dust it is necessary to know the size and density of the local dust, that depends on its composition. We consider that the material density of the different pure species are the following: ρ_REF = 3000 kg m⁻³, ρ_Fe = 7800 kg m⁻³, ρ_Si = 3000 kg m⁻³, ${\rho }_{{\text{H}}_2\text{O}}=$900 kg m⁻³, ρ_CO = 900 kg m⁻³. We also assume that the dust porosity is 50% (p = 0.5). Using these assumptions the average density of a particle located at distance r from the star can be computed from the knowledge of the surface densities of the different solid species (σ_d,i). We obtain: $${\rho }_{\text{part}}=(1-p)\frac{{\sigma }_{\text{d}}}{{\sum }_i^{}\frac{{\sigma }_{\text{d},i}}{{\rho }_i}}.$$(12)

Once the local dust average density is known, it is possible to compute dust growth. Dust growth is taken into account, including coagulation and fragmentation, using the simplified two-population model of Birnstiel et al. (2012) and is also used in Drążkowska et al. (2016) and Drążkowska & Alibert (2017). Birnstiel et al. (2012) propose a physically motivated model that mimics the evolution of the local dust size distribution of dust (see Appendix B in Birnstiel et al. 2012), which only depends on the representative size of the population that dominates the total mass budget (and is close to the maximum size of the population).At each radius in the disk, we keep track of the largest particle radius (a). We start from monomers with radius a = 1 micrometer. The largest particles grow at a rate (Brauer et al. 2008) $$\frac{\text{d}a}{\text{d}t}=\frac{{\rho }_{\text{d}}{\Delta }_u}{{\rho }_{\text{part}}},$$(13)

where ρ_d is the volume density of dust in the disk midplane, ρ_part is the average density of one dust grain, and Δ_u their relative velocity (Eq. (10) of Birnstiel et al. 2012). This growth is limited by three main processes: (1) fragmentation due to turbulent velocity, (2) too rapid drift due to gas drag, and (3) fragmentation due to differential drift velocity. The particle radius increases following Eq. (13) up to the smallest of the three following quantities (Birnstiel et al. 2012):

• the maximum radius due to turbulent fragmentation, a_t:

$$a_{\text{T}}=0.37\frac{2{\sigma }_{\text{g}}{U}_{\text{f}}^2}{3\pi {\rho }_s\alpha C_s^2};$$(14)

• the maximum radius due to drifting (particle drift before they can grow further):

$$a_{\text{d}}=0.55\frac{2{\sigma }_{\text{d}}{V}_k^2}{\gamma \pi {\rho }_s{C}_s^2},$$(15)

where γ is the local pressure gradient $$\gamma =\left|\frac{\text{d}\ln (p)}{\text{d}\ln (r)}\right|;$$(16)

• the maximum particle radius achievable before fragmentation due to differential drag:

$$a_{\text{dd}}=4\frac{{\sigma }_{\text{g}}{U}_{\text{f}}{V}_{\text{k}}}{\gamma \pi {\rho }_s{C}_s^2}.$$(17)

This procedure is an approximation; however, it is a useful and numerically efficient first-order model that captures the main physical ingredients of dust growth, and fits expensive numerical simulations of dust growth quite well (see Birnstiel et al. 2012; Charnoz & Taillifet 2012). This is an essential feature of the current numerical model, as full dust growth codes coupled with full disk evolution models for 1 million years is beyond the current capacity of computers, despite several attempts (see, e.g., Charnoz & Taillifet 2012).

2.8 A simple analytical prescription for planetesimal formation

Drążkowska et al. (2016) introduces a simple physical description to describe the onset of planetesimal formation in an alpha-disk, motivated by numerical simulations of planetesimal formation induced by streaming instability (see Johansen et al. 2007; Bai & Stone 2011). We recall here the main aspect. For the streaming instability two necessary conditions should be met simultaneously. First, dust must be marginally coupled to the gas, and the higher the Stokes number, the faster the instability can take place (Bai & Stone 2011). Following Drążkowska et al. (2016), particles with Stokes number >0.01 participate in the streaming instability. The second criterion requires the dust-to-gas ratio in the midplane to be >1. This quantity, denoted β can be approximated as follows (assuming that sedimentation equilibrates with turbulent diffusion): $$\beta =\frac{{\sigma }_{\text{d}}}{{\sigma }_{\text{g}}}\sqrt{1+st/{\alpha }_{\text{mid}}}.$$(18)

Here α_mid corresponds to α in the disk midplane. This equation is valid for Stokes numbers <1, and overestimates the dust concentration for particle Stokes numbers of approximately 1 (see, e.g., Charnoz et al. 2011), where vertical displacement is controlled by Keplerian oscillation. However, as dust size increases, dust may be incorporated into planetesimals before Eq. (18) becomes invalid. So we keep the above expression. Once the two criteria are met locally (Stokes >0.01 and β > 1), the dust in a radial cell is converted into planetesimals at a rate $$\frac{\text{d}{\sigma }_{\text{p}}}{\text{d}t}=\frac{f{\sigma }_{\text{d}}}{T_{\text{orb}}},$$(19)

where σ_p is the surface density of planetesimals, T_orb is the local orbital period, and f is the mass fraction of dust converted into planetesimals per orbital period. The value of f is not very well constrained and depends on the local turbulence rate, the local pressure gradient, and the dust Stokes number. Following Bai & Stone (2011) and Drążkowska et al. (2016) we set f = 10⁻³ as a fiducial value. Once planetesimals are formed, they no longer evolve and stay at their birth location.

2.9 Fragmentation velocity thresholds

Several works have studied the fragmentation velocity thresholds for grains of different composition using either experimental results, numerical simulations, or empirical laws (Blum & Wurm 2008; Teiser & Wurm 2009; Zsom et al. 2010; Wada et al. 2009; Wada et al. 2013; Meru et al. 2013; Yamamoto et al. 2014; Gundlach & Blum 2015).

These works return a large dispersion of values for fragmentation velocities: between 1 m s⁻¹ (Blum & Wurm 2008) and 30 m s⁻¹ for silicates (Yamamoto et al. 2014), and between 10 m s⁻¹ (Gundlach & Blum (2015) and 50–80 m s⁻¹ (Yamamoto et al. 2014; Wada et al. 2013) for icy-grains.

For simplicity, we explore two values of the fragmentation velocity for silicates: U_f = 1 m s⁻¹ and U_f = 10 m s⁻¹. The value ofU_f is held constant at 10 m s⁻¹ for icy grains. We note that recent papers on CAI formation show that U_f ≃10 m s⁻¹ explains themaximum size of CAIs found in chondrites (Charnoz et al. 2015), suggesting that cold silicate dust (<1000 K) may have a low value of U_f, whereas at higher temperature the silicate dust enters into a plastic regime (>1000 K) that rendersevery collision more dissipative, thus increasing the threshold velocity for fragmentation (e.g., Metzler 2012; Jacquet 2014).

3 Gas composition

The gas starting composition is based on the work of Lodders (2003) and Pignatale et al. (2011). Overall, the Sun’s atomic composition is constrained (Anders & Grevesse 1989), whereas the gas and dust composition of the whole disk depends on the temperature and kinetics of the considered reactions. For example, H₂O, CO, and CH₄ have a particular interconnected chemistry, and their abundances at T < 600 K are unconstrained. At this temperature, equilibrium calculations predict the conversion of CO(g) and part of H2(g) into CH₄(g) and H₂O(g). However, this reaction is kinetically driven and there is a debate on the efficiency of this conversion (Lodders 2003). As a consequence, there is no accepted value for the bulk quantities of these three molecular compounds in the whole disk.

For the present paper (see Table 2) we use the method described in Pignatale et al. (2011), the Sun elemental composition from Asplund et al. (2009), and thermodynamical equilibrium. These values are very close to those derived by Lodders (2003) using the elemental solar composition of Anders & Grevesse (1989). These small differences are brought by the Sun’s different elemental values. We then consider CO as the main carrier of carbon (we do not distinguish between CO at high temperature and CH₄ at low temperature); for Fe we include the metal, the Fe-sulfides, the Fe-silicates, and the Fe-oxides, while in Lodders (2003) they sum Fe metal and Fe-sulfides only.

Table 2 shows that solar metallicity (hereafter Z_⊙) has an equivalent dust-to-gas ratio of about 0.5% for 150 K < T < 1500 K (all species but water and CO are condensed), 1.14% for 25 K < T < 150 K (all species but CO are condensed), and 1.6% for T < 25 K (all species are condensed).

4 Results

We now investigate planetesimal formation in a disk containing a dead zone. We will focus on two important parameters:

• the disk’s metallicity, which controls the amount of dust that can accumulate in the midplane. Two cases are presented, Z = Z_⊙ and Z = 2Z_⊙;

• the fragmentation velocity of silicate grains, U_f, which controls the dust growth. Two cases are presented, U_f = 1 m s⁻¹ and U_f = 10 m s⁻¹.

4.1 Case of low fragmentation velocity for silicate dust (1 ms) and Z = Z_⊙

We ran a simulation containing a dead zone using α = 10⁻⁵ in the midplane and α = 10⁻² in the active layer. We also assumed that U_f = 1 m s⁻¹ for silicate rich dust and 10 m s⁻¹ for water ice rich dust. The initial disk structure is displayed in Fig. 2.

Fig. 2 Initialdisk structure for Z = Z⊙. Shown is the surface density of the different species as a function of distance to the star. The thin black line stands for H-He gas surface density as a function of distance to star, the thick black line stands for planetesimal surface density, the red solid line is for solid silicates, the red dashed line is for silicate vapor, the blue dashed line is for water vapor, the blue solid line is for water ice, the cyan solid line is for CO ice, and the cyan dashed line is for CO vapor. The black dashed line plots the disk temperature at the start of the simulation (right-hand scale).

First stage: dust growth to pebble size

Disk evolution is shown in Figs. 3–5. A dead zone appears from the very beginning (α is displayed in Fig. 4) extending from 0.1 to about 20 au where the vertically averaged α ranges from 5 × 10⁻⁵ up to 10⁻³. Given the low value of α inside the dead zone, impact velocities are low and dust grains can sediment and grow rapidly to pebble size from 1 to 20 au, in 10 to 10⁵ yr (Fig. 5). We see a drop in the Stoke number (Fig. 3b) at the snow line that results from the low fragmentation velocity of silicate rich grains below the snow line (U_f = 1 m s⁻¹) and from the destruction of pebbles into micron-sized silicate grains at the snow line. As a consequence, inward of the snow line particles are maintained at a low value of their Stokes number, whereas outward of the snow line they can grow to values higher than 0.01 (and thus are qualified as pebbles) because of the high value of U_f for icy grains.

Second stage: planetesimal formation due to water recondensation and formation of a dust trap

In orderto form planetesimals, the dust-to-gas ratio in the midplane must be >1. The inward drift of pebbles decreases the solid surface density at r > 10 au, whereas it increases strongly close to 1 au, facilitating dust concentration in the midplane (Fig. 6). The water vapor released by the incoming pebbles at the snow line increases the local water vapor surface density by almost two orders of magnitude (Fig. 3), leading to diffusion of water vapor outside the snow line followed by recondensation into ice (known as the cold finger effect; Stevenson & Lunine 1988). This process increases the dust-to-gas ratio outward of the snow line and facilitates planetesimal formation.

A close look at the planetesimal surface density profile (Fig. 3, in red) reveals that planetesimals form preferentially in two regions where the gas surface density and pressure in the midplane are at a maximum, at and outward of the snow line (Fig. 7). Pressure maxima are known to form dust traps because there; as the pressure gradient vanishes, pebbles do not drift anymore. The two maxima come from feedback between two processes acting on the local viscosity:

• The production of water vapor just inward of the snow line increases locally the mean molecular weight, lowering the sound velocity and local viscosity (Fig. 8) and leading to a local enhancement of surface density (detailed below).

• This process is further amplified by the behavior of α in the dead zone: the vertically averaged value of α is a decreasing function of the surface density (see Eq. (11), we note that σ_dead = σ_g − σ_active; see also Fig. 8b). When the surface density increases, more material is transported through the dead zone than through the active layer. As a consequence the vertically averaged α drops, resulting in the piling up of material. This positive feedback (lower α → higher surface density → lower α) creates a pileup just inward of the snow line.

• As the gas flux is proportional to the gradient of νσ_gr^1∕2, the gas velocity increases at the location of the strong gradient (Figs. 8a and 9). This results in a local drop of gas surface density across the snow line (between the two maxima, around 0.89 au) because of mass conservation. This can be understood qualitatively as follows: the disk viscous evolution tends to smooth the function νσ_g r^1∕2; as ν exhibits radial variations, then σ_g adapts to make νσ_gr^1∕2 smooth (Fig. 8c).

It is known that dust traps can form at pressure maxima (often corresponding to maxima of surface density); however, it would be interesting to check more precisely under what conditions they form. The time evolution of the dust surface density is given by Eq. (7) and we compute analytically its sign. We assume that the gas surface density can be written ${\sigma }_{\text{g}}(r)={\sigma }_{\text{g},0}{(r/r_0)}^p$, the dust surface density ${\sigma }_{\text{d}}(r)={\sigma }_{\text{d},\text{0}}{(r/r_0)}^p$, the sound velocity $C(r)=C_0{(r/r_0)}^q$, and that $\alpha (r)={\alpha }_0{(r/r_0)}^a$, where r₀ is some reference distance. If a dust-trap occurs then dσ_d∕dt > 0, otherwise dσ_d∕dt ≤ 0. Using these expressions we first compute the gas radial velocity and maximum dust drift velocity (V_g and V_n ; see Eqs. (2) and (9), respectively): $$\begin{array}{lll}\hfill & \hfill & V_{\text{g}}=-\frac{3C_0^2{\alpha }_0}{r{\Omega }_k}{\left(\frac{r}{r_0}\right)}^{a+2q}\left(a+p+2q+2\right)\text{},\hfill \\ \hfill & \hfill & \hfill \\ \hfill & \hfill & V_n=\frac{-C_0^2{\left(\frac{r}{r_0}\right)}^{2q}}{r{\Omega }_k}\left(-\frac{p}{2}-\frac{q}{2}+\frac{3}{4}\right).\hfill \end{array}$$(20)

Then V_d is computed using Eq. (8), and inserted into the equation of d σ_d ∕dt (Eq. (7)) where only the advective term is considered while ignoring the diffusive term (because α here is assumed to be small). We obtain the local evolution of the dust surface density : $$\begin{array}{lll}\hfill & \hfill & \text{d}{\sigma }_d/\text{d}t=\frac{C_0^2{\sigma }_{d,0}{\left(\frac{r}{r_0}\right)}^{p+2q}}{4{\Omega }_k{r}^2\left(S_t^2+1\right)}(-4S_t{p}^2-12S_{t}pq-8S_t{q}^2+6S_{t}q\hfill \\ \hfill & \hfill & +9S_{\text{}t}+6{\alpha }_0{\left(\frac{r}{r_0}\right)}^a(2a^2+4ap+8aq+7a+2p^2+8pq\hfill \\ \hfill & \hfill & +7p+8q^2+14q+6)).\hfill \end{array}$$(21)

To clearly identify the effect of viscosity variations, we assume the disk is radiatively heated by the star (q = −1∕4) and also that p = −1 (consistent with our initial conditions). We obtain the following: $$\text{d}{\sigma }_d/\text{d}t=\frac{3C_0^2{\sigma }_{d,0}{\alpha }_0{\left(\frac{r}{r_0}\right)}^{a-\frac{3}{2}}}{2{\Omega }_k{r}^2\left(S_t^2+1\right)}a\left(2a+1\right)\text{}.$$(22)

The value of dσ_d∕dt is then positive for any positive value of a, or a < −1∕2. The snow line is located where α is at a minimum in Fig. 8b corresponding to a dip in solid material surface density. On the left side of the snow line a < −1∕2, while on the right side a > 0 (Fig. 8b). So on both sides of the snow line the conditions are met for an accumulation of dust. Just inward of the minimum of α we have −1∕2 < a ≤ 0, so the dust density decreases (as observed) preventing planetesimal formation. Outward of the minimum, as a > 0 the gas velocity increases, creating a positive gas surface density gradient as gas accelerates. This explains why two sites of planetesimal formation are found on either side of the snow line, where α is minimum. If we now consider variations in σ_∕rmg rather than in α, we get (for a = 0 and q = −1∕4) $$\text{d}{\sigma }_d/\text{d}t=-\frac{C_0^2{\sigma }_{d,0}{\left(\frac{r}{r_0}\right)}^{p-\frac{1}{2}}}{4{\Omega }_k{r}^2\left(S_t^2+1\right)}(p+1)\left(S_t\left(4p-7\right)-6{\alpha }_0\left(2p+3\right)\right).$$(23)

For S_t ≫ α₀ (relevant for pebble-sized dust inside a dead zone), dσ_d∕dt > 0 for − 1 < p < 7∕4. Since at a gas surface density maximum p ≃ 0, an increase in the dust surface density is produced there. This also contributes to explaining why a maximum of planetesimal production is found just inward of the snow line: indeed, it is due to the maximum of surface density.

This process has not been identified in previous studies because gas flux is assumed constant (see, e.g., Schoonenberg & Ormel 2017) and because planetesimal formation in an evolving dead zone, topped by an active layer, has not been studied.

In the planetesimal forming region pebbles represent about 30% of the total solid mass. The mass of planetesimals formed in the system reaches 19 M_⊕ at about 5 × 10⁵ yr (Fig. 10), whereas the total mass of pebbles peaks at 80 M_⊕. So only about 25% of the total amount of pebbles is converted into planetesimals. We also see that at 10⁵ yr, most planetesimals are located outward of the snow line (blue curve), whereas a small fraction (about 3 M_⊕) are inward of the snow line (red curve). After 3 × 10⁵ yr all planetesimals are now inward of the snow line. This occurs because the snow line moves outward in a disk with a DZ, due to the piling up of gas that increases the gas temperature because of viscous heating (Fig. 11). We display in Fig. 12 the location of the DZ as well as the accretion rate as the function of time. The accretion rate increases with time because of progressive heating of the dead zone. Spikes in the accretion rate are sporadic accretion processes, already identified in simulations including a dead zone (Zhu et al. 2010). They are due to the loading of the DZ that becomes gravitationally instable at some point, and then violently releases its material content to the star. This mechanism occurs in layered accretion disks containing dead zones, and has been proposed as a possible explanation to FU Orionis events (Zhu et al. 2010). Since the mechanism occurs here well after planetesimal formation, its description is beyond the scope of the current paper.

Fig. 3 Evolution of the disk with Z = 1 Z⊙ and Uf = 1 m s−1. The thin black line stands for H-He gas surface density as a function of distance to star, the thick black line stands for planetesimal surface density, the red solid line is for solid silicates, the red dashed line is for silicate vapor, the blue dashed line is for water vapor, the blue solid line is for water ice, the cyan solid line is for CO ice, and the cyan dashed line is for CO vapor. Silicate is often superposed with iron because they have the same starting concentrations. The black dashed line plots the Stokes number (right-hand scale).

Fig. 4 Vertically averaged value of α in a disk with a dead zone as a function of distance at different times. The value of α is low in thedead zone, and decreases with increasing gas surface density.

Fig. 5 Highest Stokes number for a particle as a function of time at different locations in the disk. The snow line crosses 1 au at 2 × 105 yr, resultingin a sudden drop of the Stokes number at this moment, due to pebble destruction into small silicate grains during evaporation.

Fig. 6 Dust-to-gas ratio in the disk as a function of distance at T = 10 Kyr. The black line stands for the vertically integrated ratio, and the red line is the midplane ratio.

Fig. 7 Zoom on the region close to the snow line. Top: same legend as in Fig. 3. Bottom: pressure in the disk midpane as a function of distance. The two pressure bumps are visible at 0.86 and 0.93 au. The snow line is located at 0.86 au.

Fig. 8 (a) Gas velocity as a function of distance. (b) Vertically averaged α. (c) Local viscosity. (d) σνr1∕2. The colors designates different times: black = 0.1 Kyrs, blue = 1 Kyr, pink = 10 kyrs, orange = 10 kyrs.

Fig. 9 Evolution of the surface density of gas (H,He) and water vapor (H2O) and vertically averaged value of α just inward of the snow line, as a function of time. Gas surface density (black line, left scale), water vapor surface density (blue line, left scale) and vertically averaged value of α (red, right scale).

Fig. 10 Evolution of the mass of different species in the disk as a function of time. Solid red line: mass of planetesimals inward of the snow line; solid blue line: mass of planetesimals outward of the snow line; dashed red line: mass of pebbles inward of the snow line; dashed blue line: mass of pebbles outward of the snow line; yellow solid line: total mass of solid material; dashed green line: total mass of water vapor.

4.2 Caseof high fragmentation velocity for silicate dust: 10/ms and Z = Z_⊙

Here we explore the effect of setting U_f = 10 m s⁻¹ for silicate-rich dust. We show that dust growth and planetesimal formation are strongly enhanced inward of the snow line. Starting from the same disk used in the previous section, the resulting evolution is displayed in Fig. 14. Compared with Fig. 3, differences clearly emerge. Now particles inward of the snow line are not able to grow up to a Stokes number of about 0.1, large enough to trigger efficient planetesimal formation. This result is in agreement with the simple study presented inAppendix C and with Drążkowska et al. (2016); these studies show that, for a large U_f, planetesimal formation is possible even inward of the snow line. A strong accumulation of silicate dust and planetesimals is observed around 0.1 au, in the form of two peaks, one at 0.07 au and the other at 0.15 au. They correspond to places where the gas surface density is at a maximum due to strong opacity variations (vaporisation of silicates at 0.07 au) and also to the inner edge of the dead zone at 0.15 au (where the temperature drops below 1000 K preventing thermal ionisation). As mentioned in the previous section, an increasing surface density, and sharp α variations induce dust accumulation, as happens here. As pebbles grow quickly they drift inward, and solid material is depleted outward of the inner dust traps and inward of the snow line.

At larger distance, planetesimals form at the snow line, as described in the previous section. Up to 10 M_⊕ planetesimals are formed in the hot inner regions in the first 200 kyrs (Fig. 15). A similar amount of planetesimals forms outward of the snow line (Fig. 15), but on a much shorter timescale (50 kyrs). At 300 kyrs the snow line has shifted outward enough so that the whole planetesimal population is now inward of the snow line. At 2 Myr the entire system contains about 40 M_⊕ of planetesimals inward of the snow line: it decomposes into 15 M_⊕ of planetesimals that formed outward of the snow line (ice-rich) and 35 M_⊕ of planetesimals that formed inward of the snow line.

Planetesimals below 1 au are dominated by silicate and iron material in solar proportions (about 50% each in this simulation; see Fig. 13). Planetesimals that formed just inward of the snow line have silicate, iron, and refractory elements in solar proportion, but no water. Planetesimals that formed outward of the snow line have water, silicate and iron, and refractory elements in solar proportion. A slight increase in H₂O ice is visible for planetesimals just outward of the snow line because of water vapor recondensation.

In conclusion, in a disk with U_f = 10 m s⁻¹ and solar metallicity, a massive formation of planetesimals is observed inside the dead zone, in the inner disk with up to about 30 M_⊕ of planetesimals. While this does not seem compatible with the structure of our solar system, with a total mass inside 2 au of about 2 M_⊕, it may be compatible with some exoplanetary systems with massive populations of super-Earths inward of the snow line, for example the Trappist system (Gillon et al. 2016, 2017). However, it can also be compatible with the formation of a first generation of planetesimals that differentiated and then disrupted forming the oldest iron meteorites (see, e.g., Kruijer et al. 2014).

Fig. 11 Location of snow line as a function of time in the disk.

Fig. 12 Edges of the dead zone, snow-line location, and accretion rate. Top: inner and outer radii of the dead zone (red line) and the snow line (dashed black line). Bottom: accretion rate (solar mass/year) vs. time. Spikes >1 Myrs are FU Orionis-like instabilities and triggers in the innermost region of the dead zone.

Fig. 13 Composition of planetesimals at the end of the simulation for the case Z = Z⊙ and Uf = 1 m s−1. The mass fraction of silicate and iron are often the same and are superimposed in the graph.

4.3 Case of a disk with Z = 2 Z_⊙

We now consider a disk with the same total mass as before, but multiplying the mass of condensible material by 2 (Z = 2 Z_⊙) while keeping U_f = 1 m s⁻¹. The initial dust-to-gas ratio is 1% below the snow line, about 2% from the water snow line to the CO snow line, and about 4% beyond. The disk’s evolution is shown in Fig. 17. In the first 10⁴ yr we see that higher metallicity favors a faster growth of dust (Fig. 17a) and a higher dust-to-gas ratio in the midplane compared to the Z = Z_⊙ case (Fig. 19). This is in agreement with the qualitative analytical study in Appendix C. Planetesimal formation starts at the snow line, as in the Z = Z_⊙ case, and the formation front extends progressively outward up to the outer edge of the dead zone in about 10⁴ yr. This contrasts sharply withthe Z = Z_⊙ case where planetesimals only form in a narrow ring around the snow line. Here we see strong enhancement of planetesimal formation at the snow line, between 0.7 and 2 au (as in the previous cases, due to water vapor recondensation and formation of pressure maxima) and a wide region of planetesimal production from 2 to ≃ 20 au, where planetesimal surface density scales as r^−2.5, which is much steeper than the gas initial surface density. This region is populated by material coming from much larger distances in the form of pebbles that drifted inward and that are incorporated into planetesimals.

5 Evolution of a disk with Z = 2 Z_⊙

The mass of formed planetesimals is also strongly enhanced compared to the Z = Z_⊙ case. Up to 70 M_⊕ of planetesimals are formed outward of the snow line (scaling roughly linearly with the metallicity), whereas only about three Earth masses are formed just inward of the snow line (Fig. 18). Planetesimal formation stops around 10⁵ yr. As in the previouscases, as the snow line moves outward, planetesimals that formed beyond the snow line end up inward of the snow line. At 3 × 10⁵ yr only 10 M_⊕ of planetesimals are still beyond the snow line. They are those that formed earlier between 2 and 20 au. All planetesimals that formed closeto the snow line end up inward of the snow line.

An interesting result is that planetesimals are always more abundant than pebbles (Fig. 20). At 1 au planetesimals are 1000 times more abundant than pebbles. Only at 10 au are there about 2 times more planetesimals than pebbles. Conversely, at 20 au pebbles are 10 times more abundant. This may have implications for giant planet formation as pebble accretion has been invoked as a possible way to form giant planets in a short time (see, e.g., Lambrechts et al. 2014). We see here that when planetesimals form in a disk with a dead zone in the giant planets region, pebbles represent less mass than planetesimals.

If we now set the fragmentation velocity of silicates to U_f = 10 m s⁻¹ while keeping Z = 2Z_⊙, planetesimals form efficiently inward of the snow line and several tens of Earth masses of silicate-rich planetesimals are formed below 1 au, down to 0.1 au at the inner edge of the dead zone, as in the previous case.

6 At what metallicity do planetesimals form in an annulus or over an extended region?

We have seen that in the case of Z = Z_⊙ the vertically integrated dust-to-gas ratio is about 1% beyond the snow line; this results in limited planetesimal formation, taking place in an annulus starting at the snow line and extending over ≃ 1 au. Conversely, in the Z = 2Z_⊙ case, the vertically integrated dust-to-gas ratio is about 2% beyond the snow line, resulting in an extend region of planetesimal formation from the snow line up to about 20 au.

After exploring different values of Z from 1 to 2, we found that the critical value for an extended region of planetesimal formation is 1.2 Z_⊙. This value could be anticipated from examination of the Z = Z_⊙ case where the dust-to-gas ratio in the midplane is very close to (but lower than) 1 (Fig. 6).

Fig. 14 Case of a disk containing a dead zone and with Uf = 10 m s−1 everywhere and Z = Z⊙. Each panel displays the surface density as a function of distance for different types of material. The thick black line stands for H-He gas surface density. Solid and dashed lines respectively stand for the solid and vapor form of the different species. The red line is for silicates, deep blue for H2O, light blue for CO, yellow for refractory minerals, and green for iron (see legend below figure). The black dashed line plots the Stokes number (right-hand scale).

Fig. 15 Mass of planetesimals and pebbles as a function of time in a disk with fragmentation velocity set to Uf = 10 m s−1 for silicate-rich dust.

Fig. 16 Composition of planetesimals in a disk with Uf = 10 m s−1 as a function of distance.

7 Comparison with other works

Drążkowska & Alibert (2017) studied planetesimal formation in a fully active disk. Our results should be compared to the lowest α cases covered in Drążkowska & Alibert 2017, i.e., 10⁻⁴. Similarly, we find that planetesimal formation starts at the snow line in an annulus of about 1 au in width, and that the region of planetesimal formation gets wider with increasing metallicity and decreasing α. However, the presence of a dead zone modifies the results substantially. Unlike Drążkowska & Alibert 2017 the mechanism that triggers planetesimal formation is not the back-reaction of dust onto the gas, but rather the recondensation of water vapor and the creation of pressure maxima around the snow line. The other difference comes from the global evolution of the disk. As the material accumulates in the dead zone, the dead zone region is viscously heated, which shifts the snow line zone outward (see also Zhang & Jin 2015), a feature not seen in Drążkowska & Alibert (2017) that considers a fully active disk where the temperature and gas flow are decoupled, and thus cannot exhibit such a behavior. Despite these differences, interestingly, the global picture remains similar in Drążkowska & Alibert (2017).

Schoonenberg & Ormel (2017) performed local models of dust transport and focused on the region around the snow line with fixed boundary conditions for the pebble and gas fluxes. They find that water vapor recondensation increases the surface density of icy pebbles outward of the snow line by a factor ranging from 3 to 6. We reproduce these features (Appendix A). One major difference is that our pebble flux is derived from the disk evolution, whereas it is imposed as a boundary condition in Schoonenberg & Ormel (2017). We find here a pebble flux that ranges between 10 and 100% of the gas flux in the planetesimal forming region. This value matches the values investigated in Schoonenberg & Ormel (2017), namely 10 to 80%.

Ida & Guillot (2016) suggested a different way that planetesimals can form, namely by the destruction of silicate grains into micron-sized dust at the snow line during evaporation. By comparing the timescale of drifting to the timescale of vertical diffusion, they suggest that the small grains could be concentrated enough in the midplane just after their release (on the vertical scale height of larger pebbles) so that gravitational instability takes over and form planetesimals. Although we take into account the destruction of silicate grains at the snow line, we cannot describe this effect here as we assume that the vertical scale-height of grains is always at a diffusive steady state (in other words, the steady state is reached instantaneously). So we can only confirm that the release of silicates in the form of micron-sized dust at the snow line do increases the surface density of grains. A better modelisation with time dependant vertical diffusion is needed.

Fig. 17 Evolution of the disk with Z = 2Z⊙ and Uf = 1 m s−1. The thin black line stands for H-He gas surface density as a function of distance to star, thick black for planetesimal surface density, red for solid silicates, red dashed for silicate vapor, blue dashed for water vapor, blue for water ice, cyan for CO ice, cyan dashed for CO vapor. The black dashed line plots the Stokes number (right-hand scale).

Fig. 18 Mass of planetesimal and pebbles as a function of time in a disk with fragmentation velocity set to Uf = 1 m s−1 for silicate rich dust and Z = 2 Z⊙.

Fig. 19 Dust-to-gas ratio in the disk as a function of distance at T = 10 Kyr. The black line stands for the vertically integrated ratio, and the red line is the midplane ratio.

Fig. 20 Ratio of the local pebble surface density to planetesimal surface density at different distances from the star.

8 Limitations

We have tried to link a large number of physical effects in the present work; some assumptions and simplifications may affect the overall validity of the presented results:

• The first limitation is our assumption of instantaneous evaporation at the snow line. This choice was made to save computer time and because we use composite grains. To be self-consistent we would have needed a complex expression for the evaporation expression that is material dependent. As a consequence evaporation fronts are “razor thin” here, that is non-physical. This also results in an on-off behavior of vaporisation close to the vaporisation temperature, where small temperature variations can result in big variations of the solid surface density. This assumption will be released and studied in a future work.

• The second significant limitation is that opacity is independent of the dust-to-gas ratio and is independent of changes in the dust composition. This is a common simplification made by the vast majority of works on this subject; however, it could be a major limitation because the opacity depends on the dust-to-gas ratio and on particle sizes, inducing an additional feedback between dust growth and disk temperature. The opacity should decrease as dust grows, thus inducing progressive cooling of the disk, at least in the regions dominated by viscous heating. In an early version of this work, we have attempted to scale the opacity to the local dust-to-gas ratio. This triggered large numerical instabilities. However, we did find, as expected, that the opacity globally decreases as planetesimal form; the viscous heating contribution also decreases.This important point deserves more work.

• In this work we consider a disk that is fully formed and isolated. This assumption should be revised in the future. It has been found that planetesimals form very rapidly, meaning that formation is very sensitive to the initial conditions. While we have no constraint on the size of planetesimals forming in the simulations, the present scenario of very fast accretion seems, at first glance, at odds with the formation times that can be derived for chondritic parent bodies from the absolute U/Pb ages and relative 26Al/27Al ages of their chondrules. Both ages imply that chondrites contain some chondrules formed as late as ≃3 to ≃4 Myr after the start of the solar system (Bollard et al. 2017; Villeneuve et al. 2009). These observations could be reconciled with the present model if we consider that big planetesimals (several hundreds of kilometers) form in two steps, first forming very early planetesimal cores (as modeled here) and then grow at a slower pace by accreting pebbles (analogous to chondrules) that are still significantly present in the disk after a few Ma (see Fig. 10). This would be compatible with recent thermal modelings of the evolution of the parent bodies of magmatic iron meteorites (Kaminski et al. 2017). It would also be consistent with paleomagnetic observations suggesting that the parent body of CV chondrites could have had a differentiated metallic core and a dynamo magnetic field (Weiss & Elkins-Tanton 2013). Drążkowska & Dullemond (2018) investigated a formation using the cloud infall prescription of Hueso & Guillot (2005) in a fully active disk. We plan to lead a similar study in the future.

Regarding the formation of dust traps on either side of the snow line, two objections can be formulated that temper the generality of our results:

• The viscosity gradient found in our simulation could be overestimated because of our simplistic assumption that vaporisation is instantaneous, creating an artificially sharp transition from ice to vapor. Taking into account the water-vapor saturating pressure may naturally lead to a shallower gradient of mean-molecular weight, and thus to a shallower gradient of viscosity. Due to the non-linear feedback between the different processes, it is difficult to anticipate the final result.

• The drop in the average α we find at the snow line results from the 1D modeling of the dead zone. While the dead zone model we use is taken from previous studies (Zhu et al. 2010), the vertical structure of protoplanetary disks is still a matter of debate (Turner et al. 2014), especially when considering the effects of disk winds. Changing the vertical structure of the disk could change the present results.

Finally, our last limitation is the absence of back-reaction in this work. It has been shown that back-reaction becomes an important process when the dust-to-gas ratio gets close to 1 (see, e.g., Tanaka et al. 2005; Okuzumi et al. 2012). It has been shown that gas-dust back-reactions can trigger planetesimal formation (see, e.g., Drążkowska et al. 2016; Drążkowska & Alibert 2017) in disks with higher values of α than considered here. So it could be possible that including dust-gas back-reactions may even further enhance the production of planetesimals.

9 Discussion and conclusion

We have shown that planetesimal formation presents a different picture in a fully active disk compared to a disk with a dead zone. In a fully active disk devoid of dead zone, planetesimal formation seems only possible at the water snow line (Drążkowska & Alibert 2017). Conversely, in a disk containing a dead zone, efficient planetesimal formation can occur, owing to the very low value of α in the disk midplane. In agreement with previous works, we find that the metallicity plays a key role, as does the fragmentation velocity U_f. Our conclusions can be summarized as follows:

• For a disk with one solar metallicity (and U_f = 1 m s⁻¹) an annulus of planetesimals about 1 au in width develops from the snow line (about 1 au). Pebbles represent about 30% of the solid material mass within the planetesimal annulus, and 100% beyond.

• For a disk with Z > 1.2Z_⊙ the planetesimal formation region extends from the snow line (about 1au) to about 20 au. In this case the disk is largely dominated by planetesimals, whereas pebbles represent only a small percentage of the total solid mass budget up to 10 au.

• In all cases the vast majority of planetesimals are ice-rich as they form outward of the snow line. Because the dead zone is filled with gas, its temperature increases and the snow line moves outward by about 1au in 1 Myrs. So, finally, the majority of planetesimals end up inward of the snow line.

For a conservative value of silicate fragmentation velocity (U_f ≃ 1 m s⁻¹) planetesimals do not form significantly inward of the snow line. If U_f = 10 m s⁻¹ massive planetesimal formation occurs inward of the snow line at around 0.1 au, at the inner edge of the dead zone, and close to the silicate condensation front. However, the fragmentation velocities of silicates and other materials are matters of debate.

Planetesimal formation is always very rapid (less than a few 10⁵ yr), which may be at odds with numerous measurements of chondrite ages often longer than 1 Myr (Scott & Krot 2003) but is qualitatively consistent with the rapid depletion of dust found in protoplanetary disks of the σ Orionis cluster (Ansdell et al. 2017).

Future works need to consider nebular infall (Drążkowska & Dullemond 2018) and photo-evaporation (Carrera et al. 2017), which could offer the possibility for several phases of planetesimal formation in the disk.

Acknowledgements

We are indebted to the anonymous referee for the insightful comments that increased the quality of the paper. This work was supported by ANR CRADDLE (ANR-15-CE31-0004) and by IDEX Sorbonne Paris Cité. We acknowledge the financial support from the UnivEarthS Labex program of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). R.H. acknowledges the financial support of JSPS Grants-in-Aid (JP17J01269, 18K13600). We thank J. Drążkowska, Y. Alibert, A. Morbidelli, and S. Okuzumi for very constructive comments.

Appendix A Testing the code

We performed numerous tests of the transport scheme of our code and also individually tested different physical modules. For the gas and dust transport scheme we tried to reproduce the results of Schoonenberg & Ormel (2017), in particular their Figs. 2 and 6a. It was considered an interesting test as it combines the transport of water vapor, icy grains, diffusion, and advection. In order to do this test we artificially created local boundary conditions in our global code to reproduce the gas and pebble fluxes assumed in (Schoonenberg & Ormel (2017; see their Table 1). We also implemented a water saturating vapor pressure procedure (their Eq. (11)). The result of this test is presented in Fig. A.1 and must be compared to Fig. 6A in Schoonenberg & Ormel (2017; see their H_z = H_pebble case). We reproduce similar water and silicate abundances, with the same maximum values and timescale of production, whereas our silicate peak is somewhat more spread, which may be due to our logarithmic grid and to the use of the Lax-Wendroff scheme for advection, which is known to be more dissipative than a donor-cell scheme, but more stable. Since our composite particle transport scheme and our advective scheme are completely different from the one in Schoonenberg & Ormel (2017) we think that reproducing these results validates our transport scheme. Finally, in order to test the other physical modules we tried to reproduce the Drążkowska & Alibert (2017) work. It should be noted that a direct comparison with Drążkowska & Alibert (2017) was not possible as the way they compute the disk evolution is completely different from ours and is independent from the dust evolution scheme (unlike ours). They use three disk models. The first is without gas and temperature evolution (their power-law model); the second is a non-irradiated disk model that does not consider stellar heating; and the third is the irradiated disk model, which evolves according to a fit to a precalculated disk that does not consider evaporation of volatiles. However, as shown in the present paper, with our model we reproduce the results of Drążkowska & Alibert (2017), namely that the planetesimals start to form at the snow line in a ring of a few astronomomical units that increases with decreasing alpha, and that a few tens to a few hundreds of Earth mass of planetesimals are formed. Finally we tested for the code convergence with respect to radial resolution. Owing to our logarithmic grid we obtain rapid convergence for a number of radial grid cells >100. In the present paper we present simulations with 300 radial grid cells.

Fig. A.1 Surface density versus distance to star. Blue solid line: solid water surface density, dashed blue line: water vapor surface density, solid red line: silicate surface density.

Appendix B Calculation of the diffusive flux between two neighboring cells

We detail how we compute the diffusive flux leaving each cell’s edge. We insist that the method we report below is notnew in any way (it is mostly a donor-cell scheme), we just report it for the sake of completeness. This method is at the heart of our codes as it compute explicitly and exactly the amount of material exchanged between two neighboring cells, in both directions, at every time-step. The net flux exchanged between two neighboring cells is the difference between the two outgoing diffusive fluxes. Let a radial cell have index j, its left cell index j − 1, and its right cell index j + 1. The left edge of cell j has index j − 1∕2 and the right edge of the central cell has index j + 1∕2. By integrating Eq. (7) over the width of a cell we can compute the variation of mass of cell j. The diffusive term can be separated into four mass flux contributions. At the left edge, we have the outward flux coming from cell j − 1 (called F_L+) and the inward flux coming fromcell j (called F_L−). At the right edge, we have the outward flux coming from cell j (called F_R−) and the inward flux coming from cell j + 1 (called F_R+). They read as follows: $$F_{\text{L+}}\simeq 2\pi r(j-1/2)D_{\text{d}}{\sigma }_{\text{g}}(j)\frac{{\sigma }_{\text{d}}(j-1)}{{\sigma }_{\text{g}}(j-1){\Delta }_r},$$(B.1) $$F_{\text{L-}}\simeq 2\pi r(j-1/2)D_{\text{d}}{\sigma }_{\text{g}}(j)\frac{{\sigma }_{\text{d}}(j)}{{\sigma }_{\text{g}}(j){\Delta }_r},$$(B.2) $$F_{\text{R+}}\simeq 2\pi r(j+1/2)D_{\text{d}}{\sigma }_{\text{g}}(j)\frac{{\sigma }_{\text{d}}(j+1)}{{\sigma }_{\text{g}}(j+1){\Delta }_r},$$(B.3) $$F_{\text{R-}}\simeq 2\pi r(j+1/2)D_{\text{d}}{\sigma }_{\text{g}}(j)\frac{{\sigma }_{\text{d}}(j)}{{\sigma }_{\text{g}}(j){\Delta }_r}.$$(B.4)

Here Δ_r is the cell radial width and the diffusion coefficient D_d depends on the dust size. Following Birnstiel et al. (2012) we assume that locally most of the mass is contained in the larger dust grains so that D_d is given by the diffusion coefficient of the largest local grains and is the local gas viscosity divided by the Schmidt number (we assume that the Schmidt number is 1+Stokes number, following Youdin & Lithwick 2007; see also Charnoz & Taillifet 2012). The Schmidt number is the ratio of the dust diffusivity to the gas diffusivity. Now we can compute the entering mass diffusive flux for a given species only, knowing that the mass fraction of species i in cell j is σ_d,i (j)∕σ_d(j), in cell j + 1 it is σ_d,i (j + 1)∕σ_d(j + 1), and in cell j − 1 it is σ_d,i(j − 1)∕σ_d(j − 1); we thus obtain $$F_{\text{L+},i}=\frac{{\sigma }_{\text{d},i}(j-1)}{{\sigma }_d(j-1)}{F}_{\text{L+}},$$(B.5) $$F_{\text{L-},i}=\frac{{\sigma }_{\text{d},i}(j)}{{\sigma }_d(j)}{F}_{\text{L-}},$$(B.6) $$F_{\text{R+},i}=\frac{{\sigma }_{\text{d},i}(j+1)}{{\sigma }_d(j+1)}{F}_{\text{R+}},$$(B.7) $$F_{\text{R-},i}=\frac{{\sigma }_{\text{d},i}(j)}{{\sigma }_d(j)}{F}_{\text{L-}}.$$(B.8)

The net variation of surface density of species i in cell j due to dust diffusion is obtained by simple mass balance: $$\frac{\partial {\sigma }_{\text{d},i}}{\partial t}=\frac{F_{\text{L+},i}+F_{\text{R+},i}-F_{\text{L-},i}-F_{\text{R-},i}}{2\pi r\Delta r}.$$(B.9)

Using a similar argument, the variation of surface density of species i in cell j due to dust advection (corresponding to the first term of Eq. (7)) is $$\begin{array}{lll}\frac{\partial {\sigma }_{\text{d},i}}{\partial t}\hfill & =\hfill & \frac{2\pi r(j-1/2){\sigma }_{\text{d},i}(j-1/2)V_{\text{d}}(j-1/2)}{2\pi r\Delta r}\hfill \\ \hfill & \hfill & -\frac{2\pi r(j+1/2){\sigma }_{\text{d},i}(j+1/2)V_{\text{d}}(j+1/2)}{2\pi r\Delta r}.\hfill \end{array}$$(B.10)

So the total surface density variation of species i is obtained by summing the two above contributions (diffusive and advective flux). Since all contributions are explicitly computed (incoming flux from adjacent cells, and outgoing flux from the central cell j), it is then very easy to transport any physical quantity between all cells (e.g., particule radius, Stokes number, composition) taking properly into account diffusion, self diffusion, and advection.

Appendix C Analytical results for a simple static and power-law disk

We estimate in this section what are the necessary conditions, in a static power-law disk, to trigger planetesimal formation. In particular we are especially interested in the role of α. We use in this section a fiducial disk model, with total mass about 0.02 Sun mass, with a gas surface density: $${\sigma }_{\text{g}}(r)=5000\text{\hspace{0.17em}kg\hspace{0.17em}}{\text{m}}^{-2}{\left(\frac{r}{\text{1au}}\right)}^{-1}.$$(C.1)

The vertically integrated dust-to-gas ratio is 1% for temperatures >150 K, and is set to 3% for temperatures <150 K in order to take into account ice condensation. We use a fixed temperature profile representing an infinitely thin irradiated disk: $$T(r)=\text{260K}{\left(\frac{r}{\text{1au}}\right)}^{-0.5}.$$(C.2)

This represent a somewhat standard disk, with a snow line at about 3 au.

We plot on the left side of Fig. C.1 the maximum Stokes number as a function of distance and α and the resulting dust-to-gas ratio in the disk midplane for two different models of the fragmentation velocity. In the top row U_f = 10 m s⁻¹ everywhere; instead, in the bottom row U_f = 1 m s⁻¹ for silicate rich dust (i.e., for T > 150 K) a value U_f = 10 m s⁻¹ for water ice dust (i.e., for T < 150 K). In each case most of the growth is limited by the fragmentation limit. Only at large distances (typically >100 au) the growth is limited by differential drift. The jump at about 3 au is due to the snow line, which implies a sudden jump in surface density and in fragmentation velocity (bottom row). In a disk with U_f = 10 m s⁻¹ everywhere, planetesimals can form for α < 5 × 10⁻³ and outward of the snow line. For a disk with U_f = 1 m s⁻¹ for silicates and 10 m s⁻¹ for ice, planetesimals can also form outward of the snow line, but cannot form inward; instead, for α ≃10⁻⁵ some pebblecan form inward of the snow line, but do not concentrate enough in the midplane to grow planetesimals. This simple calculation shows that it seems difficult in these conditions to grow particles inward of the snow line. Increasing the surface density may not change the global picture as the fragmentation limit does not depend on the surface density (see Sect. 2.7). Planetesimal formation is favored in the outer regions of the disk, because of the lower relative velocities at large distance (owing to lower local orbital velocity). However, we note that Fig. C.1 only plots the maximum achievable value of the Stokes number and does not consider the time to grow a particle to reach it. So, whereas the highest values of the Stokes number may be achieved outward, it takes more time to reach them as the growth rate of dust is proportional to the local dust concentration (Eq. (18)). In addition, solid particles in the outer region may also start drifting during their growth phase, thus enhancing their concentration at smaller distances, and thus accelerating dust growth closer to the Sun. So planetesimals may start forming at the snow line and a formation front may propagate outward as dust progressively grows. The local dust-to-gas ratio may be drastically modified due to large-scale dust transport, and static disk models may give only a partial vision of the situation.

Fig. C.1 Leftcolumn: maximum Stokes number of particles in the disk. Right column: maximum dust-to-gas ratio in the midplane. Upper row: Uf = 10 m s−1 everywhere in the disk. Lower row: Uf = 10 m s−1 for ice and Uf = 1 m s−1 for silicates.

References

All Tables

All Figures

	Fig. 1 Rosseland mean opacity as a function of the local temperature. The dust composition is given in Table 1.
In the text

Fig. 2 Initialdisk structure for Z = Z⊙. Shown is the surface density of the different species as a function of distance to the star. The thin black line stands for H-He gas surface density as a function of distance to star, the thick black line stands for planetesimal surface density, the red solid line is for solid silicates, the red dashed line is for silicate vapor, the blue dashed line is for water vapor, the blue solid line is for water ice, the cyan solid line is for CO ice, and the cyan dashed line is for CO vapor. The black dashed line plots the disk temperature at the start of the simulation (right-hand scale).

In the text

Fig. 3 Evolution of the disk with Z = 1 Z⊙ and Uf = 1 m s−1. The thin black line stands for H-He gas surface density as a function of distance to star, the thick black line stands for planetesimal surface density, the red solid line is for solid silicates, the red dashed line is for silicate vapor, the blue dashed line is for water vapor, the blue solid line is for water ice, the cyan solid line is for CO ice, and the cyan dashed line is for CO vapor. Silicate is often superposed with iron because they have the same starting concentrations. The black dashed line plots the Stokes number (right-hand scale).

In the text

	Fig. 4 Vertically averaged value of α in a disk with a dead zone as a function of distance at different times. The value of α is low in thedead zone, and decreases with increasing gas surface density.
In the text

	Fig. 5 Highest Stokes number for a particle as a function of time at different locations in the disk. The snow line crosses 1 au at 2 × 105 yr, resultingin a sudden drop of the Stokes number at this moment, due to pebble destruction into small silicate grains during evaporation.
In the text

	Fig. 6 Dust-to-gas ratio in the disk as a function of distance at T = 10 Kyr. The black line stands for the vertically integrated ratio, and the red line is the midplane ratio.
In the text

	Fig. 7 Zoom on the region close to the snow line. Top: same legend as in Fig. 3. Bottom: pressure in the disk midpane as a function of distance. The two pressure bumps are visible at 0.86 and 0.93 au. The snow line is located at 0.86 au.
In the text

	Fig. 8 (a) Gas velocity as a function of distance. (b) Vertically averaged α. (c) Local viscosity. (d) σνr1∕2. The colors designates different times: black = 0.1 Kyrs, blue = 1 Kyr, pink = 10 kyrs, orange = 10 kyrs.
In the text

	Fig. 9 Evolution of the surface density of gas (H,He) and water vapor (H2O) and vertically averaged value of α just inward of the snow line, as a function of time. Gas surface density (black line, left scale), water vapor surface density (blue line, left scale) and vertically averaged value of α (red, right scale).
In the text

Fig. 10 Evolution of the mass of different species in the disk as a function of time. Solid red line: mass of planetesimals inward of the snow line; solid blue line: mass of planetesimals outward of the snow line; dashed red line: mass of pebbles inward of the snow line; dashed blue line: mass of pebbles outward of the snow line; yellow solid line: total mass of solid material; dashed green line: total mass of water vapor.

In the text

	Fig. 11 Location of snow line as a function of time in the disk.
In the text

	Fig. 12 Edges of the dead zone, snow-line location, and accretion rate. Top: inner and outer radii of the dead zone (red line) and the snow line (dashed black line). Bottom: accretion rate (solar mass/year) vs. time. Spikes >1 Myrs are FU Orionis-like instabilities and triggers in the innermost region of the dead zone.
In the text

	Fig. 13 Composition of planetesimals at the end of the simulation for the case Z = Z⊙ and Uf = 1 m s−1. The mass fraction of silicate and iron are often the same and are superimposed in the graph.
In the text

Fig. 14 Case of a disk containing a dead zone and with Uf = 10 m s−1 everywhere and Z = Z⊙. Each panel displays the surface density as a function of distance for different types of material. The thick black line stands for H-He gas surface density. Solid and dashed lines respectively stand for the solid and vapor form of the different species. The red line is for silicates, deep blue for H2O, light blue for CO, yellow for refractory minerals, and green for iron (see legend below figure). The black dashed line plots the Stokes number (right-hand scale).

In the text

	Fig. 15 Mass of planetesimals and pebbles as a function of time in a disk with fragmentation velocity set to Uf = 10 m s−1 for silicate-rich dust.
In the text

	Fig. 16 Composition of planetesimals in a disk with Uf = 10 m s−1 as a function of distance.
In the text

Fig. 17 Evolution of the disk with Z = 2Z⊙ and Uf = 1 m s−1. The thin black line stands for H-He gas surface density as a function of distance to star, thick black for planetesimal surface density, red for solid silicates, red dashed for silicate vapor, blue dashed for water vapor, blue for water ice, cyan for CO ice, cyan dashed for CO vapor. The black dashed line plots the Stokes number (right-hand scale).

In the text

	Fig. 18 Mass of planetesimal and pebbles as a function of time in a disk with fragmentation velocity set to Uf = 1 m s−1 for silicate rich dust and Z = 2 Z⊙.
In the text

	Fig. 19 Dust-to-gas ratio in the disk as a function of distance at T = 10 Kyr. The black line stands for the vertically integrated ratio, and the red line is the midplane ratio.
In the text

	Fig. 20 Ratio of the local pebble surface density to planetesimal surface density at different distances from the star.
In the text

	Fig. A.1 Surface density versus distance to star. Blue solid line: solid water surface density, dashed blue line: water vapor surface density, solid red line: silicate surface density.
In the text

	Fig. C.1 Leftcolumn: maximum Stokes number of particles in the disk. Right column: maximum dust-to-gas ratio in the midplane. Upper row: Uf = 10 m s−1 everywhere in the disk. Lower row: Uf = 10 m s−1 for ice and Uf = 1 m s−1 for silicates.
In the text