Issue 
A&A
Volume 624, April 2019



Article Number  A28  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834556  
Published online  03 April 2019 
Water delivery by pebble accretion to rocky planets in habitable zones in evolving disks
^{1}
EarthLife Science Institute, Tokyo Institute of Technology,
Meguroku,
Tokyo 1528550, Japan
email: ida@elsi.jp
^{2}
Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
Meguroku,
Tokyo 1528551, Japan
Received:
1
November
2018
Accepted:
14
January
2019
Context. The ocean mass of the Earth is only 2.3 × 10^{−4} of the whole planet mass. Even including water in the interior, the water fraction would be at most 10^{−3}−10^{−2}. Ancient Mars may have had a similar or slightly smaller water fraction. What controlled the amount of water in these planets has not been clear, although several models have been proposed. It is important to clarify the control mechanism to discuss water delivery to rocky planets in habitable zones in exoplanetary systems, as well as that to Earth and Mars in our solar system.
Aims. We consider water delivery to planets by icy pebbles after the snowline inwardly passes planetary orbits. We derive the water mass fraction (f_{water}) of the final planet as a function of disk parameters and discuss the parameters that reproduce a small value of f_{water} comparable to that inferred for the Earth and ancient Mars.
Methods. We calculated the growth of icy dust grains to pebbles and the pebble radial drift with a 1D model, by simultaneously solving the snowline migration and dissipation of a gas disk. With the obtained pebble mass flux, we calculated accretion of icy pebbles onto planets after the snowline passage to evaluate f_{water} of the planets.
Results. We find that f_{water} is regulated by the total mass (M_{res}) of icy dust materials preserved in the outer disk regions at the timing (t = t_{snow}) of the snowline passage of the planetary orbit. Because M_{res} decays rapidly after the pebble formation front reaches the disk outer edge (at t = t_{pff}), f_{water} is sensitive to the ratio t_{snow}∕t_{pff}, which is determined by the disk parameters. We find t_{snow}∕t_{pff} < 10 or > 10 is important. By evaluating M_{res} analytically, we derive an analytical formula of f_{water} that reproduces the numerical results.
Conclusions. Using the analytical formula, we find that f_{water} of a rocky planet near 1 au is similar to the Earth, i.e., ~10^{−4}−10^{−2}, in disks with an initial disk size of 30–50 au and an initial disk mass accretion rate of ~(10^{−8}−10^{−7}) M_{⊙} yr^{−1} for disk depletion timescale of approximately a few M yr. Because these disks may be median or slightly compact/massive disks, our results suggest that the water fraction of rocky planets in habitable zones may often be similar to that of the Earth if icy pebble accretion is responsible for water delivery.
Key words: planets and satellites: formation / planets and satellites: terrestrial planets / protoplanetary disks
© ESO 2019
1 Introduction
Earthsize planets are being discovered in habitable zones (HZs) in exoplanetary systems. These HZs are defined as a range of orbital radius, in which liquid water can exist on the planetary surface if H_{2}O exists there. However, as long as equilibrium temperature is concerned, H_{2}O ice grains condense only well beyond the HZs, because the gas pressure of protoplanetary disks is many orders of magnitude lower than that of planetary atmosphere and the condensation temperature is considerably lower than that at 1 atm. Hereafter, we simply call H_{2}O in a solid/liquid phase “water”. For the Earth, a volatile supply by the gas capture from the disk is ruled out because observed values of rareEarth elements are too low in the Earth to be consistent with the disk gas capture (e.g., Brown 1949). Therefore, the water in the Earth would have been delivered from the outer regions of the disk during planet formation.
One possible water delivery mechanism to Earth is the inward scattering of waterbearing asteroids by Jupiter (e.g., Raymond et al. 2004). If this is a dominant mechanism of water delivery, the amount of delivered water is rather stochastic and depends on configurations of giant planets in the planetary systems. If water is not delivered, a rocky planet in a HZ may not be able to be an actual habitat. On the other hand, too much water creates a planet without continents, where the supply of nutrients may not be as effective as in the Earth. The oceans of the Earth comprise only 0.023% by weight and this amount enables oceans and continents to coexist. The mantle may preserve water in the transition zone with a comparable amount of ocean (e.g., Bercovici & Karato 2003; Hirschmann 2006; Fei et al. 2017), while the core could have H equivalent to 2% by weight of H_{2} O of the Earth (Nomura et al. 2014). However, the original water fraction of the Earth would still be very small (~ 10^{−4}−10^{−2}), even with the possible water reservoir in the interior because neither stellar irradiation at ~ 1 au (Machida & Abe 2010) nor giant impacts (Genda & Abe 2005) can vaporize the majority of the water from the gravitational potential of Earth. We note that it is inferred that Mars may have subsurface water of 10^{−4} −10^{−3} of Mars mass (e.g., di Achille & Hynek 2010; Clifford et al. 2010; Kurokawa et al. 2014). From the high D/H ratio observed in the Venus atmosphere, early Venus may also have had oceans of the fraction 10^{−5} −10^{−3} and lost the H_{2}O vapor through runaway greenhouse effect (e.g., Donahue et al. 1982; Greenwood et al. 2018). The order of water fraction looks similar at least between the Earth and ancient Mars. Although the estimated total mass fraction of water in the Earth and Mars has relatively large uncertainty ranging from 10^{−4} to 10^{−2}, the range is stillmuch smaller than the dispersion predicted by the waterbearing asteroid collision model; this model ranges from 10^{−5} to 10^{−1}, depending on the formation timing, history of orbital migration/eccentricity evolution of gas giant planets, and the original surface density of planetesimals (e.g., Morbidelli et al. 2000; Lunine et al. 2003; Raymond et al. 2004; O’Brien et al. 2014; Matsumura et al. 2016). It is not clear if the similar orders of water fraction between the Earth and Mars is just a coincidence.
Sato et al. (2016) investigated water delivery by icy pebble accretion. Pebble accretion has been proposed as a new mode of planet accretion (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012). Radiative transfer calculations for viscous accretion disk models show that the water snowline at ~ 170 K may migrate to inside 1 au with the grain opacity ≳mm size, when the disk accretion rate is Ṁ_{g} ≲ 10^{−8} M_{⊙} yr^{−1} (e.g., Garaud & Lin 2007; Min et al. 2011; Oka et al. 2011), which is a typical value of Ṁ_{g} of classical T Tauri stars (Hartmann et al. 1998). After the snowline inwardly passes a planetary orbit, icy pebbles can be accreted by the planet. In situ ice condensation near the planet orbit is unlikely because the disk gas there was once in the outer region before it migrated to the inner region and icy components have been already condensed to icy grains and subtracted in the outer region (Morbidelli et al. 2016; Sato et al. 2016). Sato et al. (2016) calculated the time evolution of icy pebble mass flux, the solid surface density in the disk, and the growth of a hypothetical planet at 1 au by icy pebble accretion with a 1D model. They found that pebble accretion is so efficient that the water fraction of the planet rapidly increases after the snowline passage. They assumed a static disk and artificially set the timings of snowline passage and removal of the disk that truncates pebble accretion. They found that the water mass fraction of the final planet is very sensitive to these timings and it is zero or more than 0.1 in many cases. The modest water mass fraction of 10^{−4} −10^{−2} is possible only if the disk is compact (<100 au) and the snowline passage at 1 au later than 2–4 M yr after icy dust growth starts, which could be a narrow window of the parameters.
The sensitive dependence of the final water fraction requires that the snowline migration and the decay of the icy dust surface density must be consistently calculated in an evolving disk. We use the disk evolution model based on the selfsimilar solution for accretion disks with constant viscosity parameter α (LyndenBell & Pringle 1974)^{1}. The snowline migration, disk gas decay, and growth/drift of pebbles and the associated evolution of the icy dust surface density are simultaneously calculated by a 1D disk evolution model. The growth and drift of pebbles are tracked using the singlesize approximation formulated by Ormel (2014) and Sato et al. (2016), which enables us to perform fast calculations and survey broad ranges of parameters. Using the numerical results, we also derive an analytical formula for the final water mass fraction of the planets determined by the disk model parameters.
In Sect. 2, we describe the calculation model that we used. In Sect. 3, the results of numerical simulation are shown. We derive the semianalytical formula that successfully reproduces the numerical results in Sect. 4. In Sect. 5, using the analytical formula, we study the dependence of the planetary water fraction on the disk and pebble accretion parameters, and discuss the disk parameters to realize the water fraction of 10^{−4} −10^{−2}, which corresponds to the present Earth and ancient Mars. We show that the disk parameters are not in a narrow window and are rather realized in modest disks. Sections. 6 and 7 provide the discussion and summary.
2 Method
2.1 Gas disk model
In general, an accretion disk consists of an inner region in which the viscous heating is dominated and an outer region in which irradiation from the host star is dominated. According to Ida et al. (2016), we set the disk midplane temperature for the viscousheating dominated region (T_{vis}) and the irradiation dominated region (T_{irr}) to be
where r is the distance from the host star, Ṁ_{g} is the disk gas accretion rate, which is almost independent of r except in outermost region, L_{*} and M_{*} are the luminosity and mass of the host star, respectively, and L_{⊙} and M_{⊙} are their values of the Sun. We adopt the alpha prescription for the disk gas turbulent viscosity (Shakura & Sunyaev 1973), , where h_{g} is the disk gas scale height, defined by h_{g} = c_{s}∕Ω, c_{s} and Ω are the sound velocity and Kepler frequency, respectively, and α (< 1) is a parameter that represents the strength of the turbulence. We use slightly lower T_{irr} than that in Ida et al. (2016), assuming lower opacity with millimetersized dust grains (Oka et al. 2011), because we consider relatively inner disk regions near the snowline and those pebbles which have grown in outer regions and drifted inward. If micronsized grains are assumed, the same temperature is realized with about ten times smaller Ṁ_{g}.
One fundamental assumption behind Eq. (1) is that the rate of viscous heating per unit volume scales linearly with the gasdensity, and is therefore highest at the midplane. This assumption is questioned by magnetohydrodynamic (MHD) models of protoplanetary disks, which show that accretion heating dominantly takes place on the disk surface (Hirose & Turner 2011). Recently, Mori et al. (2019) investigated this issue using a series of MHD simulations including all nonideal MHD effects, finding that the midplane temperature derived from the simulations is generally lower than that from Eq. (1) because the heat generated near the disk surface can easily be lost through radiation. Mori et al. (2019) have also found that MHD disk winds, which are not included in our disk model, take away ≃ 30% of the magnetic energy that would be available for disk heating if the winds were absent. Although there are disk evolution models accounting for the mass and angular momentum loss due to MHD disk winds (e.g., Armitage et al. 2013; Suzuki et al. 2016; Bai et al. 2016; Hasegawa et al. 2017), none of these take into account the two effects mentioned above. For this reason, we opted to adopt the more classical viscous disk model in this study. We note that the viscous accretion model serves as a good approximation of real protoplanetary disks if some hydrodynamical instabilities drive turbulence near the midplane (see Lyra & Umurhan 2018; Klahr et al. 2018, for recent reviews on hydrodynamic instabilities of protoplanetary disks).
The gas disk aspect ratios corresponding to Eqs. (1) and (2) are
where h_{g,vis} and h_{g,irr} are the gas scale height in the viscousheating dominated and irradiation dominated regions, respectively. Hereafter, we perform simulations with L_{*} = L_{⊙} and M_{*} = M_{⊙}, while we retain their dependences in the formulas. The disk region is viscousheating dominated if T_{vis} > T_{irr}. Otherwise, the irradiation dominates. The transition radius between the viscousheating and irradiation dominated regions is given by (5)
Defining the snowline by the location at ~170 K, r_{snow} ≃max(r_{snow,vis}, r_{snow,irr}), where (6) (7)
As Ṁ_{g} decreases with time, the snowline migrates inward in the viscousheating dominated region until r_{snow,vis} becomes equal to r_{snow,irr} (Ṁ_{g} ≳ 5 × 10^{−9}M_{⊙} yr^{−1} for α = 10^{−2}). For calculating the evolution of the pebble flux, it may be enough to set up a static disk distribution as in Sato et al. (2016). However, to describe the snowline migration and disk gas depletion (which determines timings of start and termination of the icy pebble supply), we need an evolving gas disk model.
Specific orbital angular momentum is proportional to a square root of orbital radius r and most of disk mass exists in the outer irradiation dominated region. Angular momentum transfer in the entire disk that determines the snowline migration and the entire disk gas depletion is regulated by the evolution of the outer disk region. The region near the snowline is not a hot region where the viscous heating is significantly higher than the irradiation heating (Eqs. (5) and (6)). So, for our purpose, the entire disk evolution model can be approximated by a irradiation dominated disk. In the irradiation dominated disk, the viscosity ν ∝ αTr^{3∕2} ∝ αr^{15∕14} (Eq. (2)). Because it is similar to the disk with the viscosity ν ∝ r with a constant α, we adopted the selfsimilar solution with ν ∝ r by LyndenBell & Pringle (1974) for the dynamical evolution of the entire disk, while we took into account the snowline evolution regulated by time evolution of the viscous heating (Eq. (6)).
In the selfsimilar solution, well inside the initial characteristic disk size (r_{d,0}), beyond which the surface density decays exponentially, the disk accretion rate is given as a function of time by (8)
where Σ_{g} is the disk gas surface density, ν is the effective viscosity at r, , (where ν_{0} is the viscosity at r_{d,0}), and Ṁ_{g,0} is the initial disk accretion rate, respectively. Inversely, the time evolution of the surface gas density Σ_{g} is given by (9)
where we included the timedependent exponential taper in the full form of the selfsimilar solution, because we need to evaluate the total disk mass. Because ν ∝ r, Σ_{g} ∝ 1∕r for r ≪ r_{d,0}.
We add the effect of the photoevaporation with the rate Ṁ_{pe}, although the standard selfsimilar solution does not have such a term (also see Sect. 3.1). We are concerned with the region near the snowline. We assume that r ≪ r_{d,0} and the photoevaporation occurs mainly in the outer region with the constant rate of Ṁ_{pe}. Accordingly, we set
Integrating Eq. (11), the disk mass at t is given by (12)
where we use r∕3ν = t_{diff}∕r_{d,0}. The above equations show that the disk accretion rate Ṁ_{g} and the disk gas surface density Σ_{g} quickly decay when Ṁ_{g} decreases to the level of Ṁ_{pe}. This photoevaporation effect avoids longtail existence of disk gas significantly longer than a few M yr.
Here, we describe the disks with the parameters, t_{diff}, Ṁ_{g,0}, Ṁ_{pe} and r_{d,0}. The disk depletion timescale (t_{dep}) and the disk gas accretion rate onto the host star (Ṁ_{g}) are better constrained by observations than the other parameters, i.e., t_{dep}~10^{6}−10^{7} yr and Ṁ_{g}~(10^{−9}−10^{−7}) M_{⊙} yr^{−1} for solartype stars (e.g., Haisch et al. 2001; Hartmann et al. 1998, 2016; Williams & Cieza 2011). We focus on the systems around solartype stars. We assume that angular momentum transfer by turbulent viscous diffusion is a major mechanism for disk depletion rather than photoevaporation. We identify t_{dep} as t_{diff}. We use Ṁ_{g,0}, which may be slightly higher than the observed averaged values of Ṁ_{g} because we want to set the snowline beyond the orbits of planetary embryos. We perform simulations with t_{diff} = 10^{6} and 3 × 10^{6} yr, Ṁ_{g,0} = 3 × 10^{−8}M_{⊙} yr^{−1}, and 10^{−7}M_{⊙} yr^{−1}. Although r_{d,0} is not observationally well constrained, Sato et al. (2016) showed that this parameter is most important for the time evolution of the pebble mass flux. For r_{d,0}, we use r_{d,0} = 30, 100 and 300 au, which correspond to the range of the observationally inferred disk size (Williams & Cieza 2011). We adopt the photoevaporation rate Ṁ_{pe} = 10^{−9} and 10^{−8}M_{⊙} yr^{−1}.
The disk diffusion timescale is (13)
where subscript “1” expresses values at 1 au, ν_{d.0} = ν(r_{d.0}) is viscous coefficient at the outer edge of the disk, and ν = ν_{1} ⋅ (r∕r_{1}). From this equation, α is calculated as (14)
We note that the value of α that we use in our simulations depends on the choice of the parameters t_{diff} and r_{d,0}, but it is within a reasonable range, α = 10^{−3}−3 × 10^{−2}. In our formulation for disks, we prefer the setting of the observable valuables, t_{diff} and r_{d,0}, within the observationally inferred ranges to a simple assumption of a constant α.
In Fig. 1, we show time evolution of the transition radius r_{vis/irr}, the snowline r_{snow}, the disk characteristic radius , and the pebble formation front r_{pff}, which are derived in Sect. 3.1 as Eq. (36), for six typical parameters of disk evolution. In general, r_{pff} in magenta lines evolves faster than in the green lines and the pebble flux quickly decays after r_{pff} exceeds , as we discuss in Sect. 3.1. The snowline r_{snow} in the blue lines shrinks with time as long as it is in the viscousheating region (inside of r_{vis/irr} in the red lines).
Fig. 1 Time evolution of the transition radius r_{vis/irr} (Eq. (5); red lines), snowline r_{snow} = max(r_{snow,vis}, r_{snow,irr}) (Eqs. (6) and (7); blue lines), disk characteristic radius (green lines), and pebble formation front r_{pff} (Eq. (36); magenta lines) for various disk evolution parameters. The parameters, r_{d,0} in units of auand Ṁ_{g,0} in units of M_{⊙} yr^{−1} are labeled in the individual panels. The other parameters are the same for the four panels, i.e., Ṁ_{pe} =10^{−9} M_{⊙} yr^{−1} and t_{diff} = 10^{6} yr. 

Open with DEXTER 
2.2 Dust/pebble evolution model
We adopt the method with the single size approximation formulated by Ormel (2014) and Sato et al. (2016) to calculate the pebble mass flux. The size distribution of icy particles obtained by a fullsize simulation is generally peaky at some size and the peak size depends on r. In the single size approximation, icy dust particles have a single size that depends on r, corresponding to the peak size.
We briefly summarize the dust evolution model (for more details, see Sato et al. 2016). We set the initial particle surface density as Σ_{p} = Z_{0}Σ_{g}. We adopt Z_{0} = 0.01. The evolution of the dust/pebble surface density (Σ_{p}) and the peaked massof the particles (m_{p}) are calculated by the equation,
where is the particle radius, ρ_{int} is the internal density of icy particles, v_{r,d} and Δv_{pp} are the radial and relative velocities of the particles at the midplane, respectively, and h_{p} is the scale height of the particles given by (e.g. Youdin & Lithwick 2007) (17)
The Stokes number defined by St = t_{s}Ω is an important dimensionless number in pebble growth and radial drift, where t_{s} is the stopping time that represents the timescale of particles’s momentum relaxation due to gas drag, given by (18)
where is the thermal velocity, λ_{mfp} is the mean free path of gas particles, and m_{g} is the gas molecule mass. The mean free path is expressed by λ_{mfp} = m_{g}∕(σ_{mol}ρ_{g}), where σ_{mol} = 2.0 × 10^{−15}cm^{2} is the molecular collision cross section.
Relative velocity of dust particles Δv_{pp} is given by (19)
where Δv_{B}, Δv_{r}, Δv_{ϕ}, Δv_{z}, and Δv_{t} are the relative velocities induced by Brownian motion, radial drift, azimuthal drift, vertical settling, and turbulence, respectively (for detailed expressions, see Sato et al. 2016). We assume perfect sticking for Δv_{pp} < 30m s^{−1}; otherwise, we set the righthand side of Eq. (16) to be zero.
The radial drift velocity of dust particles is (20)
where v_{K} = rΩ_{K} is the Kepler velocity and η is a deviation of gas rotation velocity from Kepler velocity, which we set (21)
If the effect of the disk gas accretion is taken into account, 2St in the numerator of Eq. (20) is replaced by (2St + α) (e.g., Ida et al. 2016). However, it is neglected here, because the numerical simulation shows St≳0.1 for migrating pebbles (Sect. 3.1). The radial drift timescale is (22)
The pebble mass flux through the disk is given by (23)
where Σ_{p} is the pebble surface density in the pebble migrating region and we assume St ^{2} ≪ 1.
2.3 Pebble accretion onto planets
We use the same formulas as Sato et al. (2016) for the pebble accretion rate onto planets, assuming that the planets are already large enough for pebble accretion in the “settling regime” (e.g., Ormel & Klahr 2010; Guillot et al. 2014). Ormel & Klahr (2010) derived the cross section of pebble accretion as (24)
where M_{pl} is the planetary embryo mass and Δv is the relative velocity between the embryo and pebbles. The 3D pebble accretion rate Ṁ_{pl} onto the planetary embryo is given by (25)
where ρ_{p} is the spatial mass density of the pebbles and we used Eqs. (23) and (24). The parameter f_{flt} expresses the mass fraction of the accretion flow onto the planet in the pebble flux through the disk before the accretion, which is called a “filtering factor.” We note that f_{flt} does not directly depend on Δv, b_{set}, and St.
When b_{set} > h_{p}, the accretion is 2D and in Eq. (25) is replaced by 2b_{set}Σ_{p}. Accordingly, a complete filtering factor is (26)
While f_{flt} depends on b_{set} in 2D mode, it still does not depend directly on St and Δv. In our simulations, we set the planetary embryos with the masses and semimajor axes identical to Venus, Earth, and Mars (eccentricities are set to be zero). In these cases, the relative velocity is given by Δv = (3∕2)b_{set}Ω (“Hill regime”). Substituting this into Eq. (24), (27)
In the numerical simulations, we use a more general formula including “Bondi regime” and the cutoff parameter () for large St cases (Ormel & Kobayashi 2012). The dependence on Δv cancels in the accretion rate both in 2D and 3D cases and Ṁ_{pl} does not directly depend on b_{set} and St in the 3D case. For a small planet, the accretion is in 3D mode. The accretion becomes 2D mode when (28)
The filtering factor is given by (29)
where we used Eq. (21).
We take into account the decrease in dust surface density and pebble flux due to pebble accretion onto the planets. After the snowline passage, Eq. (15) is rewritten in the grid, where the planets exist, as (30)
where Δr is the radial grid size. After the icy pebble accretion onto the planet starts, the pebble mass flux decreases discontinuously at the planet orbits, according to the accretion onto the planets.
2.4 Simulation parameters
The parameters for our simulations are the initial disk characteristic size r_{d,0}, the diffusion timescale of the disk t_{diff}, the initial disk gas accretion rate Ṁ_{g,0}, and the photoevaporation rate Ṁ_{pe}. As we explained in Sect. 2.1, the other disk parameters are calculated by these parameters. We perform simulations with
 (1)
Ṁ_{g,0} = 3 × 10^{−8}, 10^{−7}M_{⊙} yr^{−1},
 (2)
Ṁ_{pe} = 10^{−9}, 10^{−8}M_{⊙} yr^{−1},
 (3)
r_{d,0} = 30, 100, 300 au,
 (4)
t_{diff} = 10^{6}, 3 × 10^{6} yr.
The initial total gas disk mass is given by Eq. (9) as (31)
which ranges from 0.06 M_{⊙} to 0.6 M_{⊙} with the above parameters. Because we start with relatively large Ṁ_{g,0}, the initial disk mass is relatively large. As we show in Sect. 4, initially massive disks tend to produce dry planets, while small mass disks tend to produce waterrich planets.
The initial dust surface density distribution is simply given by 0.01 Σ_{g,0}. With the gas disk evolution, we simultaneously calculate the pebble growth from the dust disk and its migration. The dust disk is depleted by the formation and radial drift of pebbles.
We set Venus, Earth, and Mars analogs in circular orbits with the same masses and semimajor axes as the current Venus, Earth, and Mars in solar system. For the results presented in this work, we assume M_{*} = M_{⊙} and L_{*} = L_{⊙}, while we retain the dependence on M_{*} and L_{*} in the equations. For othermass stars, the dependence of Ṁ_{g,0}, Ṁ_{pe}, r_{d,0}, and t_{diff} on the stellar mass also have to be considered.
3 Numerical results
3.1 Icy grain/pebble evolution
We calculate the growth and radial drift of icy grains and pebbles following the method by Sato et al. (2016), which is described in Sect. 2.2. The evolution of icy particles (grains and pebbles) found by simulations of Sato et al. (2016) is summarized as follows:
 1.
The growth timescale of a particle with mass m_{p} is well approximated by the simple formula with Δv_{pp} ≃ Δv_{t} in the Epstein regime (also see Takeuchi & Lin 2005; Brauer et al. 2008) as (32)
where Z_{0} is the initial particletogas ratio in the disk. The timescale of growth from μmsized grains to pebbles by several orders of magnitude in radius is (33)
 2.
The drift timescale becomes shorter as m_{p} (equivalently, Stokes number St) increases (Eq. (22)), while the growth timescale is independent of m_{p} (Eq. (32)). When t_{drift} becomes smaller than t_{grow}, the particle drift effectively starts and Σ_{p} starts being sculpted. From Eqs. (22) and (32) with Z_{0} replaced by the particletogas ratio of drifting pebbles (Z), the equilibrium Stokes number of drifting pebbles is (34)
where we used a typical value of the solidtogas ratio of migrating pebble, Z ~ 0.1 Z_{0} (Ida et al. 2016). We adopt Z_{0} as the initial particletogas ratio and adopt Z_{0} = 0.01 as a nominal value.
 3.
The sensitive rdependence of t_{grow,peb} results in “an insideout formation” of pebbles; formation is earlier in the inner region and the formation front migrates outward (also see Brauer et al. 2008; Okuzumi et al. 2012; Birnstiel et al. 2012; Lambrechts & Johansen 2014). After thepebble formation front reaches the characteristic disk radius (r_{d,0}), Σ_{p} rapidly decays uniformly in the disk because the supply from further outer regions to the region at r ~ r_{d,0} is limited.
As we show below, the time (t_{pff}) at which the pebble formation front reaches the characteristic disk radius is a very important parameter. The timescale for the pebble formation front to reach the radius of r_{d,0} is given with Eq. (33) by (35)
In other words, the pebble formation front radius is given as a function of time (t) by (36)
The timescale t_{pff} is determined only by r_{d,0} in the disk parameters. It is independent of the other disk parameters such as t_{diff} and Ṁ_{g,0}.
Once the formation front reaches r_{d,0}, the supply of solid materials from further outer regions is limited. Accordingly, the pebble formation and drift from there results in the decay of Σ_{p} near r_{d,0}. The solid surface density at r < r_{d,0} also decays because it is contributed from drifting pebbles formed near the formation front. Thus, the decay rate of Σ_{p} for t > t_{pff} is regulated by the pebble growth near r_{d,0}, and is approximated as . From this relation, Z = Σ_{p}∕Σ_{g} ∝ t^{−1} at t > t_{pff}, which is approximated as . However, the numerical results show a faster decrease due to the effect of finite r_{d,0} (Sato et al. 2016). The numerically obtained time evolution of Z at r= r_{d,0} with r_{d,0} = 30, 100 and 300 au is shown in Fig. 2. For r_{d,0} = 30, 100, and 300 au, t_{pff} ≃ 3.3 × 10^{4}, 2.0 × 10^{5}, and 1 × 10^{6} yr, respectively (Eq. (36)). We fit the numerical results as (37)
Because t_{pff} ≳ t_{diff} and the effect of a finite r_{d,0} is not important for r_{d,0}≳300 au (Eq. (36)), the r_{d,0}dependence as a factor of (300 au∕r_{d,0}) is reasonable. Although the value of γ_{2} may includean uncertainty due to the disk model and the singlesize approximation, we show in Sect. 5 that the predicted function of the water fraction depends only weakly on γ_{2}.
We note that we assume that collisions between icy pebbles always result in sticking and the pebble size is determined by the drift limit. If the pebble growth is limited by bouncing collisions or collisional fragmentation, St is determined by the threshold velocity for bouncing or fragmentation, which is lower than the value in Eq. (34), and Eq. (36) depends on the threshold St. It is a widely accepted idea that pebbles made of H_{2} O ice grains have a high sticking threshold velocity (Wada et al. 2009; Gundlach & Blum 2015), but this is not the case if the grains are mantled by poorly sticky CO_{2} (Musiolik et al. 2016). Recent studies have found that this CO_{2}induced fragmentation can have important implications for the growth of pebbleaccreting protoplanets (Johansen et al. 2019)and for the observational appearance of protoplanetary disks (Okuzumi & Tazaki 2019). In a future work, we will study how the bouncing and fragmentation barriers affect the water delivery to rocky planets.
Fig. 2 Time evolution of Z = Σ_{p}∕Σ_{g} obtained by the numerical simulation with Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 10^{6} yr. The magenta, blue, and green solid lines show Z obtained at r_{d,0} by the numerical results of with r_{d,0} = 30, 100 and 300 au. The dashed lines indicate gradients for individual cases given by Eqs. (37) and (38), where the absolute values of Z are arbitrary. 

Open with DEXTER 
3.2 Water fraction
With the simulated pebble mass flux Ṁ_{peb}, we calculate the growth rate of a planet due to icy pebble accretion (Ṁ_{pl}). The filtering factor is defined by f_{flt} = Ṁ_{pl}∕Ṁ_{peb}. When the snowline passes each planetary orbit, we switch on the icy pebble accretion onto the planet. We assume that the 1:1 ratio of rocky to icy fraction of icy pebbles. We set the rocky planets with Venus, Earth, and Mars masses when the snowline passes 0.72, 1.0, and 1.52 au, respectively.
After the snowline passes a planetary orbit, icy pebble accretion starts. When the cumulative accreted mass by icy pebbles is ΔM_{pl}, the total ice mass in the planet is (1∕2)ΔM_{pl}. The water fraction of the planet at the mass M_{pl} is given by (39)
where M_{pl,0} is the planetary mass at the snowline passage, M_{pl} = M_{pl,0} + ΔM_{pl}. If ΔM_{pl} becomes much larger than M_{pl,0}, the water fraction saturates to f_{water} ≃ 1∕2.
We repeat the simulations with different disk parameters t_{diff}, Ṁ_{g,0}, r_{d,0}, and Ṁ_{pe} to investigate how the water fraction of the final planets depends on these parameters and which values of the parameters produce the water fraction consistent with the terrestrial planets in the solar system. As we show in Sect. 4, the water fraction of final planets is insensitive to M_{pl,0}.
Figure 3 shows the time evolution of water fraction for the models with t_{diff} = 10^{6} yr, Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. The left, middle, and right panels show the results of r_{d,0} = 30, 100, and 300 au, respectively. We first explain general evolution pattern of the water fraction. In all cases, the water fraction rapidly increases once the snowline passes the planetary orbit and the icy pebble accretion starts. It is saturated to its asymptotic value, even sufficiently before completion of disk gas depletion. From Eqs. (6) and (14), (40)
At t = 0 with Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, the snowline is located outside the orbit of Mars. As shown in Fig. 1, the snowline migrates inward and passes through the planetary orbits one after another as the disk accretion rate Ṁ_{g} decreases with time as Eq. (10). The snowline passage always occurs in the order of Mars at 1.52 au, Earth at 1.0 au, and Venus at 0.72 au. Because the initial r_{snow,vis} is closer to the planetary orbits for larger r_{d,0}, the snowline passage is earlier for larger r_{d,0}, as shown in Fig. 3. The water fraction is saturated when Σ_{p} and Ṁ_{peb} significantly decay. The rapid decay starts at t ~ t_{pff} and Eq. (33) shows that t_{pff} is as short as approximately a few × 10^{5} yr for r_{d,0} = 100 au. Even for r_{d,0} = 300 au, the decay starts at t_{pff} ≃ 10^{6} yr before disk depletion. Therefore, the evolution of the water fraction evolution is mainly reduced by the consumption of icy dust grain reservoir rather than by the disk gas depletion.
Figure 3 shows a clear trend that the final water fraction is lower for smaller r_{d,0}. This trend is explained by a comparison between the snowline passage time t_{snow} and t_{pff}, as follows. The water fraction due to pebble accretion rapidly increases until Σ_{p} decays by more than an order of magnitude, which corresponds to, say, t > 10 t_{pff} (Eq. (37)). If t_{snow} > 10 t_{pff}, f_{water} can be ≪ 1. While t_{pff} is smaller for a smaller r_{d,0} (Eq. (33)), t_{snow} is larger (Eq. (40)). The latter implies that the disk is warmer for a smaller r_{d,0}. The viscousheating increases as the blanketing effect by optical depth (∝ Σ_{g}) increases. In Fig. 3, Ṁ_{g,0} and t_{diff} are fixed. Smaller r_{d,0} means smaller α (Eq. (13)) and larger Σ_{g} (Eq. (9)), resulting in a warmer disk. Thereby, the final water fraction is lower for smaller r_{d,0}.
We examine the condition of t_{snow}∕t_{pff} > 10 or <10 in more detail. In the case of r_{d,0} = 30 au, t_{pff} ≃ 4 × 10^{4} yr (Eq. (33)) and t_{snow} is identified by the timing at which f_{water} starts rapid increase, which is ≃1 × 10^{6}, 2.5 × 10^{6}, and 5 × 10^{6} yr for the Mars, Earth and Venus analogs, respectively (the left panel of Fig. 3). Because t_{snow} > 10 t_{pff} in this case, f_{water} ~ 10^{−2} even for the outermost Mars analog. For the Earth and Venus analogs, f_{water} is further smaller. In the case of r_{d,0} = 300 au, t_{pff} ≃ 1 × 10^{6} yr (Eq. (33)) and t_{snow} ≃ 1 × 10^{5}, 1 × 10^{6} and 2 × 10^{6} yr for the Mars, Earth, and Venus analogs, respectively (right panel of Fig. 3). Because t_{snow} < 10 t_{pff}, f_{water} ~ 1∕2 for all of Mars, Earth, and Venus. In the middle panel of Fig. 3, r_{d,0} = 100 au and t_{pff} ≃ 2 × 10^{5} yr (Eq. (33)). In this case, t_{snow} < 10 t_{pff} and f_{water} ≪ 1 for the Venus analog, while t_{snow} > 10 t_{pff} and f_{water} ~ 1∕2 for the Mars analog.
The condition of t_{snow}∕t_{pff} > 10 or < 10 also explains other results of the final water fraction of planets because this condition represents how much icy materials remain at the snowline passage. (In Sect. 5, we revisit this condition.) Figure 4 is the time evolution of water fraction of t_{diff} = 10^{6} yr (left panel) and t_{diff} = 3 × 10^{6} yr (right panel), respectively. The other disk parameters are the same. For both cases, 10 t_{pff} ≃ 2 × 10^{6} yr (Eq. (33)). The snowline is already inside Mars’ orbit from the beginning (t = 0) of the calculations, that is, t_{snow} = 0, which results in f_{water} ≃ 1∕2 for the Mars analog. The snowline passage time t_{snow} is smaller than10 t_{pff} only for the Venus analog in the case of t_{diff} = 10^{6}yr (left panel), and for both the Earth and Venus analogs in the case of t_{diff} = 3 × 10^{6} yr (right panel). This explains the results in Fig. 4.
Figure 5 shows the results of Ṁ_{g,0} = 3 × 10^{−8} M_{⊙} yr^{−1} (left) and 10^{−7} M_{⊙} yr^{−1} (right) with r_{d,0} = 100 au, t_{diff} = 10^{6} yr, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. Again, 10 t_{pff} ≃ 2 × 10^{6} yr for both cases. Smaller Ṁ_{g,0} with the fixed t_{diff} = 10^{6} yr^{−1} means earlier passage of the snowline, resulting in a lower water fraction. Figure 6 shows the results of Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} (left) and Ṁ_{pe} = 3 × 10^{−9}M_{⊙} yr^{−1} (middle) and Ṁ_{pe} = 10^{−8}M_{⊙} yr^{−1}. With larger Ṁ_{pe}, the outer disk region is truncated at shorter radius. The disk truncation reduces the reservoir of icy materials, resulting in a lower water fraction, for the planets (the Earth and Venus analogs) with t_{snow} ≳ 10 t_{pff} ≃ 2 × 10^{6} yr.
In Sect. 4, we derive a semianalytical formula to predict the water fraction of planets after disk depletion. We show that the total mass (M_{res}) of the icy dust materials preserved in the disk determines the final water fraction because they eventually drift to the inner regions and pass the planetary orbits^{2}. Since our gas disk model is analytical, we can analytically evaluate t_{snow}. We already know the analytical expression of t_{pff} given by Eq. (36) is a good approximation. We also semianalytically derive how the icy grain surface density evolves as Eq. (37). A synthesis of these results enables us to derive the semianalytical formula for f_{water}.
Finally we point out that f_{water} is higher in the order of Mars, Earth, and Venus analogs according to the snowline passage timing, except for the fully saturated cases where f_{water} ~ 1∕2 for all the planets. We show that the final water fraction is insensitive to the embryo mass. The timing is the most important factor for f_{water}.
In the next section, we show that the derived semianalytical expression of f_{water} reproduces the numerical results. Using the analytical expression, we clarify the dependence of the final water fraction on the disk parameters and pebble accretion parameters such as the initial planetary mass and Stokes number of accreting pebbles in Sect. 5. We also survey the disk parameter range that may reproduce f_{water} comparable to that of the terrestrial planets in our solar system.
Fig. 3 Time evolution of water fraction for the models with t_{diff} = 10^{6} yr, Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. Left, center, and right panel: results of r_{d.0} = 30 100, and 300 au, respectively. The blue, red, and green lines represent the Earth, Mars and Venus analogs, respectively. 

Open with DEXTER 
Fig. 4 Time evolution of water fraction for the models with r_{d.0} = 100 au, Ṁ_{g,0} = 3 × 10^{−8}M_{⊙} yr^{−1} and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. Left and right panels: results of t_{diff} = 10^{6} and 3 × 10^{6}yr, respectively. 

Open with DEXTER 
4 Analytical formula for planetary water fraction
The planetary water fraction is calculated by estimating the cumulative mass of accreted icy pebbles ΔM_{pl} (Eq. (39)). We can simply estimate it as ΔM_{pl} ≃ f_{flt}(M_{pl}) × M_{res}, where M_{res} is the total icy dust mass preserved in outer disk regions at the snowline passage because the pebble flux integrated from the snowline passage time t_{snow} to infinity (effectively, to t_{diff}) must be similar to M_{res}.
Equation (39) is approximated to be (41)
where we use f_{flt}(M_{pl}) ≃ f_{flt}(M_{pl,0}) for ΔM_{pl} ≪ M_{pl,0} in the upper line. For the range of the disk parameters we use in numerical simulations, the accretion is mostly in the 3D regime (Eq. (28)) and f_{flt} ∝ M_{pl}. In that case, Eq. (41) shows that f_{water} is independent of M_{pl} both in the cases of ΔM_{pl} < M_{pl,0} and ΔM_{pl} > M_{pl,0}, while it depends on the disk parameters. Even if the transition to 2D accretion occurs at a smaller planet mass (which happens for smaller values of α), the mass dependence of f_{water} is still weak: it is proportional to for ΔM_{pl} < M_{pl,0} and independentof M_{pl} otherwise. We can combine the two limits in Eq. (41) by a simple formula as
Fig. 5 Time evolution of water fraction for the models with r_{d.0} = 100 au, t_{diff} = 10^{6} yr and Ṁ_{pe} = 10^{−9} M_{⊙} yr^{−1}. Left and right panels: results of Ṁ_{g,0} = 3 × 10^{−8} and 10^{−7}M_{⊙} yr^{−1}, respectively. 

Open with DEXTER 
Fig. 6 Time evolution of water fraction for the models with r_{d.0} = 100 au, t_{diff} = 10^{6} yr and Ṁ_{ini} = 10^{−7}M_{⊙} yr^{−1}. Left, center, and right panels: results of Ṁ_{pe} = 10^{−9}, 3 × 10^{−9}, and 10^{−8}M_{⊙} yr^{−1}, respectively. 

Open with DEXTER 
In Fig. 7, f_{water,sim} obtained byour simulations is compared with Eq. (42) with the simulated M_{res}. This figure showsthat Eq. (42) reproduces the numerical results very well. Both the numerical results and Eq. (42) show f_{water} ∝ M_{res} for f_{water} ≪ 1∕2. The icy dust mass M_{res} depends sensitively on the disk parameters, Ṁ_{g,0}, Ṁ_{pe}, t_{diff}, and r_{d,0}. On the other hand, f_{flt} in 2D is independent of these disk parameters, while that in 3D depends only weakly (). In the case of ΔM_{pl} ≪ M_{pl,0}, f_{water} ≃ (f_{flt}∕M_{pl,0})M_{res}. As we discussed, f_{flt} ∕M_{pl,0} is almost constant, so that f_{water} should be almost proportional to M_{res} for ΔM_{pl} ≪ M_{pl,0}.
Because the analytical solution given by Eq. (42) reproduces the numerical results, we next derive an analytical formula for the icy dust mass M_{res} at the snowline passage as follows. From Eq. (6), the snowline passes the planet orbit at a_{pl} when (43)
Using Eq. (14), (44)
From Eq. (10), the snowline passage time is (45)
From Eq. (12), the remaining gas disk mass at the snowline passage is given by (46)
Because M_{res} ~ Z M_{g,snow}, we obtain (47)
where we approximated Z as from Eq. (37), (Eq. (36)), and γ = 1 + 0.15(300 au∕r_{d,0}) (Eq. (38)). Substituting the filtering factor f_{flt} given by Eq. (29) with St = 0.1 and M_{res} estimated above into Eq. (42), we can analytically estimate the water fraction from the disk parameters, Ṁ_{g,0}, Ṁ_{pe}, t_{diff}, and r_{d,0} as (Eqs. (29) and (39))
Fig. 8 Comparison of the water fraction obtained by the numerical simulations (f_{water,sim}) with the analytical estimate (f_{water,anly}). The blue squares, red circles, and green cross represent the results for the Earth, Mars, and Venus analogs, respectively. If f_{water,sim} obtained by the numerical simulation is smaller than 10^{−5}, we put f_{water,sim} = 10^{−5} because the numerical results would include a numerical uncertainty for such small values of f_{water,sim}. 

Open with DEXTER 
where M_{res} is given by Eq. (47), f_{flt} = min(f_{flt,3D}, f_{flt,2D}), and
Except for the value of γ, which was fitted by the numerical results, the other derivations are analytical. We note that f_{water} ∝ M_{res} ∝ Z_{0}. Around metalrich stars, f_{water} would be larger.
In Fig. 8, we compare the analytically estimated water fraction f_{water,anly} with the numerically simulated f_{water,sim}. In addition to the Earth analogs (M_{pl,0} = 1M_{⊕}) at 1 au, we also plot the results for the Mars analogs (M_{pl,0} = 0.11M_{⊕}) at 1.52 au and the Venus analogs (M_{pl,0} = 0.82M_{⊕}) at 0.72 au. For the Earth and Mars analogs, f_{water,anly} reproduces f_{water,sim} within a factor of several, while the water fraction varies by several orders of magnitude, except for some runs for which the water fraction is overestimated by the analytical formula. The agreement both for the Earth and Mars analogs strongly suggests that the mass and semimajor axis dependences are also reproduced. In the case of Venus analogs, f_{water,anly} is larger by a factor of a few to a few tens than f_{water,sim}. When the Earth analog increases its mass by capture of pebbles, it captures a significant fraction of the pebble flux and the accretion of the Venus analogs that reside in inner region of the Earth analogs can be significantly decreased. This effect is included in the numerical simulation, while it is not taken into account in the analytical formula.
Fig. 7 Water fraction f_{water,sim} and icy dust mass preserved in the outer disk at the snowline passage M_{res} obtained by our simulations. The blue squares show the simulation results for the Earth analogs with different values of Ṁ_{g,0}, Ṁ_{pe}, t_{diff}, and r_{d,0} (Sect. 2.4). The dotted curve represents the analytical solution given by Eq. (42). 

Open with DEXTER 
5 Dependence of water fraction on disk parameters
Using the analytical formula, we investigate how the water fraction (f_{water}) is determined by the disk parameters. Figure 9 shows f_{water} for a planet at a_{pl} = 1.0 au as a function of the disk parameters, Ṁ_{g,0} and r_{d,0}. The other parameters are Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and t_{diff} = 3 × 10^{6} yr^{−1}. The panels a show the dependence on M_{pl,0} for St = 0.1. The panels b show the dependence on St for M_{pl,0} = 1M_{⊕}. The planetsformed with the parameters in the red region are very dry, f_{water} ≲ 10^{−4}. Those in the green region have a modest amount of water, f_{water} ~ 0.1, and those in the blue region are icy planets, f_{water} ≃ 1∕2. The yellow and orange regions represent f_{water} ~ 10^{−4}−10^{−2}, which corresponds to the water fraction of the current Earth and that estimated for ancient Mars.
We find that f_{water} is the most sensitive to Ṁ_{g,0} and r_{d,0}. The common features in the contours in Fig. 9 are that (i) f_{water} is lower for smaller r_{d,0} and larger Ṁ_{g,0}, and (ii) f_{water} ~ 10^{−4}−10^{−2} is realized at r_{d,0} ~ 30−50 au and Ṁ_{g,0} ≳ 10^{−8}M_{⊙} yr^{−1}. This is consistent with the conclusion by Sato et al. (2016) that such f_{water} can be realized for compact disks with sizes <100 au and late snowline passage (>2−4 M yr). We show that the disk parameters for f_{water} ~ 10^{−4}−10^{−2} is not in a very narrow window.
As we showed in Sect. 4, f_{water} is regulated by M_{res} and M_{res} is sensitive to Ṁ_{g,0} and r_{d,0}, especially through the parameter t_{snow}∕t_{pff} (Eq. (47)). Because pebble accretion is fast to realize f_{water} ~ 10^{−4}−10^{−2}, M_{res} has to have significantly decayed until t ~ t_{snow}. For small r_{d,0}, t_{pff} is small (Eq. (36)), while t_{snow} is large (Eq. (45)) due to small Ṁ_{g,snow} (Eq. (44)). For large Ṁ_{g,0}, t_{snow} is small (Eq. (44)) while t_{pff} is the same (Eq. (36)). As a result, Σ_{p} decays more quickly (Eq. (37)) for small r_{d,0} and large Ṁ_{g,0}.
Figure 9 shows the dependences of f_{water} on the pebble accretion parameters: the initial planetary mass (M_{pl,0}) and the Stokes number (St) of pebbles that accrete onto the planet.
In the numerical simulations, the Stokes number of pebbles that accrete onto the planet is calculated by growth and radial drift of pebbles in an evolving disk. According to the simulations, the Stokes number of radially drifting particles is ~ 0.1 at early times (Eq. (34)), but decreases with time as the dust and gas disks evolve. As Eq. (37) shows, Σ_{p} decays more rapidly thanΣ_{g} for t > t_{pff}. Accordingly,Z decreases and the equilibrium Stokes number of migrating pebbles become smaller (Eq. (34); also see Sato et al. 2016). As we pointed out in Sect. 3.1, if a fragmentation/rebound barrier limits the icy pebble growth, St also becomes small. With this reasoning, we also showed plots with St = 0.01 and 0.001 inFig. 9.
Figure 9 shows that f_{water} is almost independent of M_{pl,0} and St. The weak dependence on M_{pl,0} is explained by the following arguments. The water fraction f_{water} is ~ (f_{flt}∕M_{pl,0})M_{res} in the case of ΔM_{pl} ≪ M_{pl,0}. The dust mass M_{res} is independent of M_{pl}. As already explained in Sect. 4, the factor f_{flt}∕M_{pl,0} is independent of M_{pl,0} in 3D accretion and only weakly depends on M_{pl,0} in 2D accretion ().
The weak dependence of f_{water} on Stokes number shown in Fig. 9 is the result of the assumption that M_{res} is independent of St; it depends on St only weakly (Eq. (29)) through f_{flt}. The parameter γ could be dependent on St because the sculpture rate of Σ_{p} at r = r_{pff} may depend on St. Figure 10 shows the dependence on γ_{2}, assuming that the functional form of γ = 1 + γ_{2}(300 au∕r_{d,0}) still holds. Because the sculpture rate may be lower for smaller values of St, we tested the cases of γ_{2} = 0.05 and 0.1, in addition to the nominal case of γ_{2} = 0.15. The result is insensitive to γ_{2} for γ_{2} ≳ 0.1, while the sculpture rate is lower and accordingly the f_{water} is relatively higher for γ_{2} ~ 0.05. The detailed functional form of γ is left for future works.
Figure 11 shows the dependence on the other disk parameters, t_{diff} and Ṁ_{pe}. For larger t_{diff}, the evolution of the snowline is slower and the snowline passage is later (Eq. (45)). Before the passage, the icy dust disk has been more sculpted. As a result, M_{res} is smaller (Eq. (47)) because t_{snow} ∝ t_{diff} and γ ≥ 1. For smaller Ṁ_{pe}, the disk is hotter and the snowline passage is later (t_{snow} is larger), resulting in smaller M_{res}. If M_{res} is smaller, f_{water} is smaller. The water fraction corresponding to the current Earth is 10^{−4}−10^{−2} (the yellow and orange colored regions) are only slightly shifted to larger r_{d,0} and lower Ṁ_{g,0} for smaller t_{diff} and larger Ṁ_{pe}, because M_{res} is smaller for these parameters (Eq. (47)). We note that the parameter region with Ṁ_{g,0} < Ṁ_{pe} is empty in the right panel of Fig. 11b because the disks do not exist under that condition.
In Fig. 12, the analytically estimated f_{water} is plotted for a_{pl} = 0.72, 1.00 and 1.52 au (the Venus, Earth, and Mars analogs, respectively). In the outer region, f_{water} is generally larger due to early snowline passage (smaller t_{snow}). In other words, Ṁ_{g,snow} is larger for larger a_{pl}. As mentioned in Sect. 4, the Venus analogs may have further lower f_{water}, if we take into account the decrease in the pebble flux due to accretion by the Earth analog. As shown in Eq. (7), the snowline cannot reach the region inside 0.53 au in our disk model. The planets there are completely dry. Although the dependence on the orbital radius exists, in the case of r_{d,0} ~ 30−50 au and Ṁ_{g,0} ≳ 2 × 10^{−8} M_{⊙} yr^{−1}, f_{water} ~ 10^{−4}−10^{−2} for both the Earth and Mars analogs.
In Sect. 3, we pointed out that the simple condition of t_{snow}∕t_{pff} < 10 or > 10 discriminates between the waterrich case (f_{water} ~ 1∕2) and the waterpoor case. The pebble formation timescale t_{pff} is given by Eq. (36) and the snowline passage time t_{snow} is given by Eq. (45). Both t_{pff} and t_{snow} are independent of M_{pl} and St, which is consistent with Fig. 9. The condition of t_{snow} = 10 t_{pff} is shown on the r_{d,0} –Ṁ_{g,0} plane in Fig. 13. We show the dependences on t_{diff}, Ṁ_{pe}, and a_{pl}. A comparison of this figure with Figs. 11 and 12 show that the simple condition approximately reproduces the more detailed evaluation except for r_{d,0} > 100 au. As shown in Fig. 1, for r_{d,0} > 100 au, the disk radius expands with the pebble formation front radius r_{pff}, so that the sculpture of icy dust reservoir is delayed, compared to the estimation of t_{pff} at r_{d,0}. The calculation of M_{res} does not have this problem.
As we have shown, f_{water} is the most sensitive to Ṁ_{g,0} and r_{d,0} and almost independentof the other parameters of disks and pebble accretion. A robust result would be that the water fraction inferred for the present Earth and ancient Mars, f_{water} ~ 10^{−4}−10^{−2}, is realized at r_{d,0} ~ 30−50 au and Ṁ_{g,0} ≳ 2 × 10^{−8} M_{⊙} yr^{−1}, which may correspond to median disks of classical T Tauri stars or slightly compact and massive disks.
Fig. 9 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}). Panels a: dependence on the planet mass (M_{pl,0}). Panels b: Stokes number (St) of pebbles. The other parameters are Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and t_{diff} = 3 × 10^{6} yr. In the panelsa, we use St = 0.1. In the panelsb, M_{pl,0} = 1M_{⊕}. The middle panels in a and b are identical. The color scales are log_{10}f_{water}. 

Open with DEXTER 
Fig. 10 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}), where γ_{2} = 0.05, 0.10, and 0.15. The other parameters are St = 0.1, M_{pl,0} = 1M_{⊕}, , and t_{diff} = 3 × 10^{6} yr. The color scales are log_{10}f_{water}. 

Open with DEXTER 
Fig. 11 Analytical estimate of f_{water} at 1 au as a function of r_{d,0} and Ṁ_{g,0}. Panels a: The dependence on the disk diffusion timescale t_{diff}. Panels b: dependence on the photoevaporation rate (Ṁ_{pe}). The other parameters are St = 0.1 and M_{pl,0} = 1M_{⊕}. In the panels a, we used Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 3 × 10^{6} yr in the panels b. The middle panels of panels a and the right panel of panels b are identical. In the right panel of panel b, the region with Ṁ_{g,0} < Ṁ_{pe} is empty because the disks do not exist under that condition. The color scales are log_{10} f_{water}. 

Open with DEXTER 
Fig. 12 Analytical estimate of f_{water} as a functionof r_{d,0} and Ṁ_{g,0}, at 0.72, 1.00, and 1.52 au. The other parameters are St = 0.1, M_{pl,0} = 1M_{⊕}, Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 3 × 10^{6} yr. The color scales are log_{10}f_{water}. 

Open with DEXTER 
6 Discussion
In this section, we comment on the effect of pebble isolation mass, which we did not include in our simulations. A planet with relatively large mass (M_{pl}) can make a dip in the gas disk along the planetary orbit to prevent the pebbles from passing the orbit. The threshold mass is called “pebble isolation mass” (M_{peb,iso}; Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018). When the planetary mass reaches the isolation mass, pebble accretion onto the planet with M_{pl} > M_{peb,iso} and other planets inside the planetary orbit is truncated and the increases in their water fraction are stalled.
Bitsch et al. (2018) derived a detailed expression of the pebble isolation mass, (51)
For example, M_{peb,iso} ≃ 4.1M_{⊕} for t_{diff} = 3 × 10^{6} yr, Ṁ_{g,0} = 10^{−8}M_{⊙} yr^{−1}, r_{d,0} au, h∕r ≃ 0.0246 (at 1 au), and α ≃ 3 × 10^{−3}. If M_{pl,0} ≳several M_{⊕}, the effect of pebble isolation mass is not negligible. Figure 14 shows the water fraction for M_{pl,0} = 3, 5, and 10 M_{⊕}. We stop pebble accretion when M_{pl} reaches M_{peb,iso}. Because α is smaller for smaller r_{d,0} (Eq. (14)) and h∕r is lower for smaller Ṁ_{g,0} and higher α (Eq. (3)), M_{peb,iso} is smaller for smaller r_{d,0} and Ṁ_{g,0}. The dry regime (redcolored regime) in the left bottom part of the plots represent the cases of M_{pl,0} > M_{peb,iso}. Because t_{snow} is small in the low Ṁ_{g,0} regions, f_{water} rapidly increases because of a high pebble flux until M_{pl} increases to M_{peb,iso}. Thereby, the edge in f_{water} at M_{pl,0} ≃ M_{peb,iso} is sharp.
For the same t_{diff}, Ṁ_{g,0}, and r_{d,0} as the above, . The filtering rate for M_{peb,iso} is in 2D case (Eq. (29)). If ice giants or cores of gas giants are formed, even before their mass exceeds M_{peb,iso}, the pebble flux is reduced by ~50% by individual ice giants.
If Jupiter’s core is formed in the outer region, it shuts down the pebble mass flux into the terrestrial planet region. Morbidelli et al. (2016) proposed that the pebble flux truncation by the formation of Jupiter accounts forthe dichotomy of our solar system – the total solid mass contained in Jupiter, Saturn, Uranus, and Neptune is about 50 times larger than the total mass of terrestrial planets. Jupiter’s core must be formed at M_{res} ≳ 10 M_{⊕}. From Eq. (47) with M_{d,0} ~ 10^{−8} M_{⊙} yr^{−1} and t_{diff} ~ 3 × 10^{6} yr, , which is for r_{d,0} ~ 100 au. Hence, the Jupiter’s core formation time t_{jup} must be ≲ t_{diff}. If this is the case, the dichotomy of our solar system could be created. However, the water fraction of the terrestrial planets is a problem in this case, in addition to a problem of too fast type I migration of the core (Matsumura et al. 2017). Our results show that f_{water} rapidly increases after t = t_{snow} and becomes saturated before t ~ t_{diff}. If t_{snow} > t_{jup}, f_{water} = 0 for terrestrial planets. Otherwise, it is likely that f_{water} is already close to the saturated value because the increase of f_{water} is very rapid. It is difficult for giant planet formation to directly produce modestly low values of water mass fraction (~ 10^{−4}−10^{−2}) corresponding to the Earth and ancient Mars. The modestly low values are attained by disk parameters with t_{snow} ~ 10 t_{pff} as we showed.
We also point out that D/H ratio is expected to be radially uniform among the Earth, asteroids, and comets if water is delivered by icy pebbles that form in disk outer regions and drift all the way to the host star. However, observations show that in our solar system, D/H ratios of Oort cloud comets are clearly higher than the Earth (e.g., Marty 2012). One possibility to reconcile the discrepancy is the shutdown of the pebble flux by Jupiter formation (Kruijer et al. 2017). Because pebble formation front migrates outward, the isotope ratios of drifting pebbles before and after Jupiter formation, which are the building materials for terrestrial planets and those for comets, respectively, should be different. A more careful comparison should be necessary between planet formation model and cosmochemical data.
Fig. 13 Condition of t_{snow} = 10 t_{pff} as a functionof r_{d,0} and Ṁ_{g,0}. Panels a: dependence on the disk diffusion timescale t_{diff}. Panels b: photoevaporation rate (Ṁ_{pe}). Panels c: planetary orbital radius (a_{pl}). The right regions from the curves represent waterrich regions (f_{water} ~ 1∕2). In the panelsa, we used Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and a_{pl}= 1 au, t_{diff} = 3 × 10^{6} yr and a_{pl}= 1 au in the panelsb, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and t_{diff} = 3 × 10^{6} yr in the panels c. 

Open with DEXTER 
Fig. 14 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}), including the effect of the pebble isolation mass. The other parameters are Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, t_{diff} = 3 × 10^{6} yr, and St = 0.1. Pebble accretion is halted when M_{pl} reaches M_{peb,iso}. The dry regime (red) in the left bottom part of the plots represents the cases of M_{pl,0} > M_{peb,iso}. 

Open with DEXTER 
7 Summary
If water is not delivered to rocky planets in HZs, the planets cannot be actual habitats because H_{2} O ice condenses in the disk regions well beyond the HZs. In the pebble accretion model, accretion of icy pebbles after the snowline passage may be a primary mechanism to deliver water to the rocky planets.
In this paper, we have investigated the water delivery to rocky planets by pebble accretion around solartype stars through a 1D simulation of the growth of icy dust grains to pebbles and the pebble radial drift in an evolving disk. We assume that the planetary embryos did not migrate significantly and consist of pure rock, which means that accretion of ice starts when the snowline migrates inward and passes the planetary orbit due to disk evolution. Our previous paper, Sato et al. (2016), pointed out that the water fraction of the final planets are determined by the timings of the snowline passage through the planetary orbit (t = t_{snow}) and disk gas depletion because pebble accretion is fast and efficient. While Sato et al. (2016) used a simple static disk model, we used the evolving disk model due to viscous diffusion based on the selfsimilar solution with constant viscous α (LyndenBell & Pringle 1974) and simultaneously calculated pebble formation/drift/accretion and snowline migration with the disk model. Because the snowline migration is correlated with global disk diffusion in the evolving disk model, we found that for a water fraction of the final planets (f_{water}), the snowline passage time (t_{snow}) relative to the time (t_{pff}) at which pebble formation front reaches the disk outer edge is more important than that relative to the disk gas depletion timescale. Our simulation shows that the ice dust mass (M_{res}) preserved in the disk outer region at t = t_{snow} determines the water fraction of the final planets. The accreted ice mass to the planet is estimated by ~ (1∕2)f_{flt}M_{res}, where the filtering factor f_{flt} is the fraction of the pebble mass flux that is accreted onto the planet. Because M_{res} rapidly decreases after t ~ t_{pff}, t_{snow}∕t_{pff} > 10 or < 10 is crucial for final value of f_{water}. If t_{snow}∕t_{pff} > 10, M_{res} should have significantly decayed when icy pebble accretion starts at t = t_{snow}.
Using these numerical results, we derived an analytical formula for f_{water} by icy pebble accretion. In the formula, f_{water} is explicitly given as a function of the ratio t_{snow}∕t_{pff} and the disk parameters. The parameter t_{snow}∕t_{pff} is also determined by the disk parameters. As a result, f_{water} is predictedby the disk parameters, especially the initial disk mass accretion rate Ṁ_{g,0} and initial disk size r_{d,0}. It is insensitive to pebble accretion parameters such as the planet mass and Stokes number of drifting pebbles.
We found that the expected water fraction of an Earth analog near 1 au has f_{water} ~ 10^{−4}−10^{−2}, which may correspond to the value of the current Earth, in disks with an initial disk size of r_{d,0} ~ 30–50 au and the initial disk mass accretion rate Ṁ_{g,0} ~ (10^{−8}−10^{−7})M_{⊙} yr^{−1}. For Ṁ_{g,0} ≳ 2 × 10^{−8}M_{⊙} yr^{−1}, both the Earth and Mars analogs have f_{water} ~ 10^{−4}−10^{−2}, while f_{water} is generally larger for Mars than for Earth. Because these disks may be median or slightly compact/massive disks among classical T Tauri stars, our results suggest that rocky planets in HZs in exoplanatery systems around solartype stars often have a water fraction similar to the Earth if pebble accretion is responsible for water delivery, while large diversity in the water fraction due to the variations of disk conditions may also exist.
Acknowledgements
We thank Michiel Lambrechts for detailed and helpful comments. This work was supported by JSPS KAKENHI 15H02065 and 16K17661 and by MEXT KAKENHI 18H05438.
References
 Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bai, X.N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Bercovici, D., & Karato, S.i. 2003, Nature, 425, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, H. 1949, in The Atmospheres of the Earth and Planets, ed. G. P. Kuiper (Chicago: University of Chicago), 258 [Google Scholar]
 Clifford, S. M., Lasue, J., Heggy, E., et al. 2010, J. Geophys. Res. Planets, 115, E07001 [NASA ADS] [CrossRef] [Google Scholar]
 di Achille, G., & Hynek, B. M. 2010, Lunar Planet. Sci. Conf., 41, 2366 [NASA ADS] [Google Scholar]
 Donahue, T. M., Hoffman, J. H., Hodges, R. R., & Watson, A. J. 1982, Science, 216, 630 [NASA ADS] [CrossRef] [Google Scholar]
 Fei, H., Yamazaki, D., Sakurai, M., et al. 2017, Sci. Adv., 3, e1603024 [NASA ADS] [CrossRef] [Google Scholar]
 Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [NASA ADS] [CrossRef] [Google Scholar]
 Genda, H., & Abe, Y. 2005, Nature, 433, 842 [NASA ADS] [CrossRef] [Google Scholar]
 Greenwood, J. P., Karato, S.i., Vander Kaaden, K. E., Pahlevan, K., & Usui, T. 2018, Space Sci. Rev., 214, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Haisch, Jr. K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Hirschmann, M. M. 2006, Ann. Rev. Earth Planet. Sci., 34, 629 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Ida, S., & Brasser, R. 2019, A&A, 622, A202 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H., Pfeil, T., & Schreiber, A. 2018, Handbook of Exoplanets (London: Springer), 138 [Google Scholar]
 Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Nat. Acad. Sci., 114, 6712 [NASA ADS] [Google Scholar]
 Kurokawa, H., Sato, M., Ushioda, M., et al. 2014, Earth Planet. Sci. Lett., 394, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lunine, J. I., Chambers, J., Morbidelli, A., & Leshin, L. A. 2003, Icarus, 165, 1 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W., & Umurhan, O. 2018, PASP, submitted [Google Scholar]
 Machida, R., & Abe, Y. 2010, ApJ, 716, 1252 [NASA ADS] [CrossRef] [Google Scholar]
 Marty, B. 2012, Earth Planet. Sci. Lett., 313, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., Brasser, R., & Ida, S. 2016, ApJ, 818, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Chambers, J., Lunine, J. I., et al. 2000, Meteorit. Planet. Sci., 35, 1309 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Mori, S., Bai, X.N., & Okuzumi, S. 2019, ApJ, 872, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 818, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Nomura, R., Hirose, K., Uesugi, K., et al. 2014, Science, 343, 522 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, D. P., Walsh, K. J., Morbidelli, A., Raymond, S. N., & Mandell, A. M. 2014, Icarus, 239, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., & Tazaki, R. 2019, ApJ, submitted [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W. 2014, ApJ, 789, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Kobayashi, H. 2012, ApJ, 747, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icarus, 168, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Takeuchi, T., & Lin, D. N. C. 2005, ApJ, 623, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
The winddriven disk accretion that has recently been proposed (e.g., Suzuki et al. 2016; Bai et al. 2016) is commented on in Sect. 2.1.
All Figures
Fig. 1 Time evolution of the transition radius r_{vis/irr} (Eq. (5); red lines), snowline r_{snow} = max(r_{snow,vis}, r_{snow,irr}) (Eqs. (6) and (7); blue lines), disk characteristic radius (green lines), and pebble formation front r_{pff} (Eq. (36); magenta lines) for various disk evolution parameters. The parameters, r_{d,0} in units of auand Ṁ_{g,0} in units of M_{⊙} yr^{−1} are labeled in the individual panels. The other parameters are the same for the four panels, i.e., Ṁ_{pe} =10^{−9} M_{⊙} yr^{−1} and t_{diff} = 10^{6} yr. 

Open with DEXTER  
In the text 
Fig. 2 Time evolution of Z = Σ_{p}∕Σ_{g} obtained by the numerical simulation with Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 10^{6} yr. The magenta, blue, and green solid lines show Z obtained at r_{d,0} by the numerical results of with r_{d,0} = 30, 100 and 300 au. The dashed lines indicate gradients for individual cases given by Eqs. (37) and (38), where the absolute values of Z are arbitrary. 

Open with DEXTER  
In the text 
Fig. 3 Time evolution of water fraction for the models with t_{diff} = 10^{6} yr, Ṁ_{g,0} = 10^{−7}M_{⊙} yr^{−1}, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. Left, center, and right panel: results of r_{d.0} = 30 100, and 300 au, respectively. The blue, red, and green lines represent the Earth, Mars and Venus analogs, respectively. 

Open with DEXTER  
In the text 
Fig. 4 Time evolution of water fraction for the models with r_{d.0} = 100 au, Ṁ_{g,0} = 3 × 10^{−8}M_{⊙} yr^{−1} and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}. Left and right panels: results of t_{diff} = 10^{6} and 3 × 10^{6}yr, respectively. 

Open with DEXTER  
In the text 
Fig. 5 Time evolution of water fraction for the models with r_{d.0} = 100 au, t_{diff} = 10^{6} yr and Ṁ_{pe} = 10^{−9} M_{⊙} yr^{−1}. Left and right panels: results of Ṁ_{g,0} = 3 × 10^{−8} and 10^{−7}M_{⊙} yr^{−1}, respectively. 

Open with DEXTER  
In the text 
Fig. 6 Time evolution of water fraction for the models with r_{d.0} = 100 au, t_{diff} = 10^{6} yr and Ṁ_{ini} = 10^{−7}M_{⊙} yr^{−1}. Left, center, and right panels: results of Ṁ_{pe} = 10^{−9}, 3 × 10^{−9}, and 10^{−8}M_{⊙} yr^{−1}, respectively. 

Open with DEXTER  
In the text 
Fig. 7 Water fraction f_{water,sim} and icy dust mass preserved in the outer disk at the snowline passage M_{res} obtained by our simulations. The blue squares show the simulation results for the Earth analogs with different values of Ṁ_{g,0}, Ṁ_{pe}, t_{diff}, and r_{d,0} (Sect. 2.4). The dotted curve represents the analytical solution given by Eq. (42). 

Open with DEXTER  
In the text 
Fig. 8 Comparison of the water fraction obtained by the numerical simulations (f_{water,sim}) with the analytical estimate (f_{water,anly}). The blue squares, red circles, and green cross represent the results for the Earth, Mars, and Venus analogs, respectively. If f_{water,sim} obtained by the numerical simulation is smaller than 10^{−5}, we put f_{water,sim} = 10^{−5} because the numerical results would include a numerical uncertainty for such small values of f_{water,sim}. 

Open with DEXTER  
In the text 
Fig. 9 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}). Panels a: dependence on the planet mass (M_{pl,0}). Panels b: Stokes number (St) of pebbles. The other parameters are Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and t_{diff} = 3 × 10^{6} yr. In the panelsa, we use St = 0.1. In the panelsb, M_{pl,0} = 1M_{⊕}. The middle panels in a and b are identical. The color scales are log_{10}f_{water}. 

Open with DEXTER  
In the text 
Fig. 10 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}), where γ_{2} = 0.05, 0.10, and 0.15. The other parameters are St = 0.1, M_{pl,0} = 1M_{⊕}, , and t_{diff} = 3 × 10^{6} yr. The color scales are log_{10}f_{water}. 

Open with DEXTER  
In the text 
Fig. 11 Analytical estimate of f_{water} at 1 au as a function of r_{d,0} and Ṁ_{g,0}. Panels a: The dependence on the disk diffusion timescale t_{diff}. Panels b: dependence on the photoevaporation rate (Ṁ_{pe}). The other parameters are St = 0.1 and M_{pl,0} = 1M_{⊕}. In the panels a, we used Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 3 × 10^{6} yr in the panels b. The middle panels of panels a and the right panel of panels b are identical. In the right panel of panel b, the region with Ṁ_{g,0} < Ṁ_{pe} is empty because the disks do not exist under that condition. The color scales are log_{10} f_{water}. 

Open with DEXTER  
In the text 
Fig. 12 Analytical estimate of f_{water} as a functionof r_{d,0} and Ṁ_{g,0}, at 0.72, 1.00, and 1.52 au. The other parameters are St = 0.1, M_{pl,0} = 1M_{⊕}, Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, and t_{diff} = 3 × 10^{6} yr. The color scales are log_{10}f_{water}. 

Open with DEXTER  
In the text 
Fig. 13 Condition of t_{snow} = 10 t_{pff} as a functionof r_{d,0} and Ṁ_{g,0}. Panels a: dependence on the disk diffusion timescale t_{diff}. Panels b: photoevaporation rate (Ṁ_{pe}). Panels c: planetary orbital radius (a_{pl}). The right regions from the curves represent waterrich regions (f_{water} ~ 1∕2). In the panelsa, we used Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and a_{pl}= 1 au, t_{diff} = 3 × 10^{6} yr and a_{pl}= 1 au in the panelsb, and Ṁ_{pe} = 10^{−9}M_{⊙} yr^{−1} and t_{diff} = 3 × 10^{6} yr in the panels c. 

Open with DEXTER  
In the text 
Fig. 14 Analytical estimate of f_{water} at 1 au as a function of the initial disk radius (r_{d,0}) and the initial disk accretion rate (Ṁ_{g,0}), including the effect of the pebble isolation mass. The other parameters are Ṁ_{pe} =10^{−9}M_{⊙} yr^{−1}, t_{diff} = 3 × 10^{6} yr, and St = 0.1. Pebble accretion is halted when M_{pl} reaches M_{peb,iso}. The dry regime (red) in the left bottom part of the plots represents the cases of M_{pl,0} > M_{peb,iso}. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.