Planetesimal formation in an evolving protoplanetary disk with a dead zone

Context. When and where planetesimals form in a protoplanetary disk are highly debated questions. Streaming instability is considered the most promising mechanism, but the conditions for its onset are stringent. Disk studies show that the planet forming region is not turbulent because of the lack of ionization forming possibly dead zones (DZs). Aims. We investigate planetesimal formation in an evolving disk, including the DZ and thermal evolution. Methods. We used a 1D time-evolving stratiﬁed disk model with composite chemistry grains, gas and dust transport, and dust growth. Results. Accretion of planetesimals always develops in the DZ around the snow line, due to a combination of water recondensation and creation of dust traps caused by viscosity variations close to the DZ. The width of the planetesimal forming region depends on the disk metallicity. For Z = Z (cid:12) , planetesimals form in a ring of about 1 au width, while for Z > 1 . 2 Z (cid:12) planetesimals form from the snow line up to the outer edge of the DZ (cid:39) 20 au. The efﬁciency of planetesimal formation in a disk with a DZ is due to the very low effective turbulence in the DZ and to the efﬁcient piling up of material coming from farther away; this material accumulates in region of positive pressure gradients forming a dust trap due to viscosity variations. For Z = Z (cid:12) the disk is always dominated in terms of mass by pebbles, while for Z > 1 . 2 Z (cid:12) planetesimals are always more abundant than pebbles. If it is assumed that silicate dust is sticky and grows up to impact velocities ∼ 10 ms − 1 , then planetesimals can form down to 0.1 au (close to the inner edge of the DZ). In conclusion the DZ seems to be a sweet spot for the formation of planetesimals: wide scale planetesimal formation is possible for Z > 1 . 2 Z (cid:12) . If hot silicate dust is as sticky as ice, then it is also possible to form planetesimals well inside the snow line.


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  or direct growth of very porous aggregates (Arakawa &Nakamoto 2016 andTatsuuma 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 St > 10 −2 ) 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 −1 (Blum & Wurm 2008) is commonly considered, whereas for water ice U f 10 m s −1 (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.

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 multispecies, planetesimal formation, dust growth, dust sublimation, and layering due to the DZ. We detail the different aspects below.

System
The protoplanetary disk is made of two components: gas and dust.  (25, 150, 1500, 1550, 1650. 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 2 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.

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 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) 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 the value of σ g,i is evolved as 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: with D g standing for the diffusivity of the gas, taken to be equal to the disk viscosity (Hueso & Guillot 2005or Fromang & Papaloizou 2006

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 explicitly compute 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 different materials surface density in dust form. We have by mass conservation: 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)) 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 αC 2 s /(Ω k S c) with α taken in the dead zone (as most dust resides in the dead zone rather than in the active layer), and S c 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) where V g is given by Eq.
(2), while (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. The dust composition is given in Table 1.

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 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 −3 < α < 10 −2 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 −5 -10 −3 in the midplane, and could be >10 −2 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 −2 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 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).

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

Dust sublimation and condensation
As previously noted, the sublimation temperature at which the different dust species (REF, Fe, Si, H 2 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  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.

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 −3 , ρ Fe = 7800 kg m −3 , ρ Si = 3000 kg m −3 , ρ H 2 O = 900 kg m −3 , ρ CO = 900 kg m −3 . 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: 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 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 : -the maximum radius due to drifting (particle drift before they can grow further): where γ is the local pressure gradient -the maximum particle radius achievable before fragmentation due to differential drag: 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). . 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 dustto-gas ratio in the midplane to be >1. This quantity, denoted β can be approximated as follows (assuming that sedimentation equilibrates with turbulent diffusion): 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 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 −3 as a fiducial value. Once planetesimals are formed, they no longer evolve and stay at their birth location.
For simplicity, we explore two values of the fragmentation velocity for silicates: U f = 1 m s −1 and U f = 10 m s −1 . The value of U f is held constant at 10 m s −1 for icy grains. We note that recent papers on CAI formation show that U f 10 m s −1 explains the maximum 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 renders every collision more dissipative, thus increasing the threshold velocity for fragmentation (e.g., Metzler 2012; Jacquet 2014).

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 2 O, CO, and CH 4 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 4 (g) and H 2 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 4 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).

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 −1 and U f = 10 m s −1 .

Case of low fragmentation velocity for silicate dust (1 ms) and Z = Z
We ran a simulation containing a dead zone using α = 10 −5 in the midplane and α = 10 −2 in the active layer. We also assumed that U f = 1 m s −1 for silicate rich dust and 10 m s −1 for water ice rich dust. The initial disk structure is displayed in Fig. 2 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 −5 up to 10 −3 . 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 5 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 −1 ) 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 order to 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 νσ g r 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 νσ g r 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 σ g (r) = σ g,0 (r/r 0 ) p , the dust surface density σ d (r) = σ d,0 (r/r 0 ) p , the sound velocity C(r) = C 0 (r/r 0 ) q , and that α(r) = α 0 (r/r 0 ) a , where r 0 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): 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 : 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: 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). 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) For S t α 0 (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 5 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 5 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 5 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 A50, page 8 of 16 S. Charnoz et al.: Planetesimal formation in an evolving disk  occurs here well after planetesimal formation, its description is beyond the scope of the current paper.

Case of high fragmentation velocity for silicate dust:
10/ms and Z = Z Here we explore the effect of setting U f = 10 m s −1 for silicaterich 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 in Appendix C and with Drążkowska et al. (2016); these studies show that, for a large   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 2 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 −1 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(Gillon et al. , 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).

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 −1 . 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 4 yr we see that higher metallicity favors a faster growth of dust (Fig. 17a) and a higher dust-togas 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 4 yr. This contrasts sharply with the 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.

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 5 yr. As in the previous cases, as the snow line moves outward, planetesimals that formed beyond the snow line end up inward of the snow line. At 3 × 10 5 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 close to 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 −1 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.

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 A50, page 10 of 16  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).  Drążkowska & Alibert (2017) studied planetesimal formation in a fully active disk. Our results should be compared to the A50, page 11 of 16 A&A 627, A50 (2019) Fig. 17. Evolution of the disk with Z = 2Z and U f = 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). 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 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. 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%.  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.

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 nonphysical. 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-togas 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 A50, page 13 of 16 A&A 627, A50 (2019) 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.

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 −1 ) 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 −1 ) planetesimals do not form significantly inward of the snow line. If U f = 10 m s −1 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 5 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.

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.

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 not new 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 from cell 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 R+ 2πr( j + 1/2)D d σ g ( j) σ d ( j + 1) σ g ( j + 1)∆ r , (B.3) F R− 2πr( j + 1/2)D d σ g ( j) σ d ( j) σ g ( j)∆ 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 The net variation of surface density of species i in cell j due to dust diffusion is obtained by simple mass balance: 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 ∂σ d,i ∂t = 2πr( j − 1/2)σ d,i ( j − 1/2)V d ( j − 1/2) 2πr∆r (B.10) − 2πr( j + 1/2)σ d,i ( j + 1/2)V d ( j + 1/2) 2πr∆r .
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.