Free Access
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/0004-6361/201834556
Published online 03 April 2019

© ESO 2019

1 Introduction

Earth-size 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 H2O exists there. However, as long as equilibrium temperature is concerned, H2O 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 H2O 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 rare-Earth 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 water-bearing 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 H2 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 H2O 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 water-bearing 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 self-similar solution for accretion disks with constant viscosity parameter α (Lynden-Bell & 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 single-size 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 semi-analytical 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 viscous-heating dominated region (Tvis) and the irradiation dominated region (Tirr) 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 hg is the disk gas scale height, defined by hg = cs∕Ω, cs 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 Tirr than that in Ida et al. (2016), assuming lower opacity with millimeter-sized 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 micron-sized 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 hg,vis and hg,irr are the gas scale height in the viscous-heating 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 viscous-heating dominated if Tvis > Tirr. Otherwise, the irradiation dominates. The transition radius between the viscous-heating and irradiation dominated regions is given by (5)

Defining the snowline by the location at ~170 K, rsnow ≃max(rsnow,vis, rsnow,irr), where (6) (7)

As g decreases with time, the snowline migrates inward in the viscous-heating dominated region until rsnow,vis becomes equal to rsnow,irr (g ≳ 5 × 10−9M 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 ναTr3∕2αr15∕14 (Eq. (2)). Because it is similar to the disk with the viscosity νr with a constant α, we adopted the self-similar solution with νr by Lynden-Bell & 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 self-similar solution, well inside the initial characteristic disk size (rd,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 rd,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 time-dependent exponential taper in the full form of the self-similar solution, because we need to evaluate the total disk mass. Because νr, Σg ∝ 1∕r for rrd,0.

We add the effect of the photoevaporation with the rate pe, although the standard self-similar solution does not have such a term (also see Sect. 3.1). We are concerned with the region near the snowline. We assume that rrd,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ν = tdiffrd,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 long-tail existence of disk gas significantly longer than a few M yr.

Here, we describe the disks with the parameters, tdiff, g,0, pe and rd,0. The disk depletion timescale (tdep) and the disk gas accretion rate onto the host star (g) are better constrained by observations than the other parameters, i.e., tdep~106−107 yr and g~(10−9−10−7) M yr−1 for solar-type stars (e.g., Haisch et al. 2001; Hartmann et al. 1998, 2016; Williams & Cieza 2011). We focus on the systems around solar-type stars. We assume that angular momentum transfer by turbulent viscous diffusion is a major mechanism for disk depletion rather than photoevaporation. We identify tdep as tdiff. 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 tdiff = 106 and 3 × 106 yr, g,0 = 3 × 10−8M yr−1, and 10−7M yr−1. Although rd,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 rd,0, we use rd,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−8M yr−1.

The disk diffusion timescale is (13)

where subscript “1” expresses values at 1 au, νd.0 = ν(rd.0) is viscous coefficient at the outer edge of the disk, and ν = ν1 ⋅ (rr1). 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 tdiff and rd,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, tdiff and rd,0, within the observationally inferred ranges to a simple assumption of a constant α.

In Fig. 1, we show time evolution of the transition radius rvis/irr, the snowline rsnow, the disk characteristic radius , and the pebble formation front rpff, which are derived in Sect. 3.1 as Eq. (36), for six typical parameters of disk evolution. In general, rpff in magenta lines evolves faster than in the green lines and the pebble flux quickly decays after rpff exceeds , as we discuss in Sect. 3.1. The snowline rsnow in the blue lines shrinks with time as long as it is in the viscous-heating region (inside of rvis/irr in the red lines).

thumbnail Fig. 1

Time evolution of the transition radius rvis/irr (Eq. (5); red lines), snowline rsnow = max(rsnow,vis, rsnow,irr) (Eqs. (6) and (7); blue lines), disk characteristic radius (green lines), and pebble formation front rpff (Eq. (36); magenta lines) for various disk evolution parameters. The parameters, rd,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 tdiff = 106 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 full-size 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 = Z0Σg. We adopt Z0 = 0.01. The evolution of the dust/pebble surface density (Σp) and the peaked massof the particles (mp) are calculated by the equation,

where is the particle radius, ρint is the internal density of icy particles, vr,d and Δvpp are the radial and relative velocities of the particles at the midplane, respectively, and hp is the scale height of the particles given by (e.g. Youdin & Lithwick 2007) (17)

The Stokes number defined by St = tsΩ is an important dimensionless number in pebble growth and radial drift, where ts 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 mg is the gas molecule mass. The mean free path is expressed by λmfp = mg∕(σmolρg), where σmol = 2.0 × 10−15cm2 is the molecular collision cross section.

Relative velocity of dust particles Δvpp is given by (19)

where ΔvB, Δvr, Δvϕ, Δvz, and Δvt 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 Δvpp < 30m s−1; otherwise, we set the right-hand side of Eq. (16) to be zero.

The radial drift velocity of dust particles is (20)

where vK = 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 Mpl 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 fflt 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 fflt does not directly depend on Δv, bset, and St.

When bset > hp, the accretion is 2D and in Eq. (25) is replaced by 2bsetΣp. Accordingly, a complete filtering factor is (26)

While fflt depends on bset 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)bsetΩ (“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 bset 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 rd,0, the diffusion timescale of the disk tdiff, the initial disk gas accretion rate g,0, and the photo-evaporation 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−7M yr−1,

  • (2)

    pe = 10−9, 10−8M yr−1,

  • (3)

    rd,0 = 30, 100, 300 au,

  • (4)

    tdiff = 106, 3 × 106 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 water-rich 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 other-mass stars, the dependence of g,0, pe, rd,0, and tdiff 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 mp is well approximated by the simple formula with Δvpp ≃ Δvt in the Epstein regime (also see Takeuchi & Lin 2005; Brauer et al. 2008) as (32)

    where Z0 is the initial particle-to-gas ratio in the disk. The timescale of growth from μm-sized grains to pebbles by several orders of magnitude in radius is (33)

  • 2.

    The drift timescale becomes shorter as mp (equivalently, Stokes number St) increases (Eq. (22)), while the growth timescale is independent of mp (Eq. (32)). When tdrift becomes smaller than tgrow, the particle drift effectively starts and Σp starts being sculpted. From Eqs. (22) and (32) with Z0 replaced by the particle-to-gas ratio of drifting pebbles (Z), the equilibrium Stokes number of drifting pebbles is (34)

    where we used a typical value of the solid-to-gas ratio of migrating pebble, Z ~ 0.1 Z0 (Ida et al. 2016). We adopt Z0 as the initial particle-to-gas ratio and adopt Z0 = 0.01 as a nominal value.

  • 3.

    The sensitive r-dependence of tgrow,peb results in “an inside-out 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 (rd,0), Σp rapidly decays uniformly in the disk because the supply from further outer regions to the region at r ~ rd,0 is limited.

As we show below, the time (tpff) 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 rd,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 tpff is determined only by rd,0 in the disk parameters. It is independent of the other disk parameters such as tdiff and g,0.

Once the formation front reaches rd,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 rd,0. The solid surface density at r < rd,0 also decays because it is contributed from drifting pebbles formed near the formation front. Thus, the decay rate of Σp for t > tpff is regulated by the pebble growth near rd,0, and is approximated as . From this relation, Z = Σp∕Σgt−1 at t > tpff, which is approximated as . However, the numerical results show a faster decrease due to the effect of finite rd,0 (Sato et al. 2016). The numerically obtained time evolution of Z at r= rd,0 with rd,0 = 30, 100 and 300 au is shown in Fig. 2. For rd,0 = 30, 100, and 300 au, tpff ≃ 3.3 × 104, 2.0 × 105, and 1 × 106 yr, respectively (Eq. (36)). We fit the numerical results as (37)

where (38)

Because tpfftdiff and the effect of a finite rd,0 is not important for rd,0≳300 au (Eq. (36)), the rd,0-dependence as a factor of (300 au∕rd,0) is reasonable. Although the value of γ2 may includean uncertainty due to the disk model and the single-size 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 H2 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 CO2 (Musiolik et al. 2016). Recent studies have found that this CO2-induced fragmentation can have important implications for the growth of pebble-accreting 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.

thumbnail Fig. 2

Time evolution of Z = Σp∕Σg obtained by the numerical simulation with g,0 = 10−7M yr−1, pe =10−9M yr−1, and tdiff = 106 yr. The magenta, blue, and green solid lines show Z obtained at rd,0 by the numerical results of with rd,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 fflt = plpeb. 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 ΔMpl, the total ice mass in the planet is (1∕2)ΔMpl. The water fraction of the planet at the mass Mpl is given by (39)

where Mpl,0 is the planetary mass at the snowline passage, Mpl = Mpl,0 + ΔMpl. If ΔMpl becomes much larger than Mpl,0, the water fraction saturates to fwater ≃ 1∕2.

We repeat the simulations with different disk parameters tdiff, g,0, rd,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 Mpl,0.

Figure 3 shows the time evolution of water fraction for the models with tdiff = 106 yr, g,0 = 10−7M yr−1, and pe = 10−9M yr−1. The left, middle, and right panels show the results of rd,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−7M 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 rsnow,vis is closer to the planetary orbits for larger rd,0, the snowline passage is earlier for larger rd,0, as shown in Fig. 3. The water fraction is saturated when Σp and peb significantly decay. The rapid decay starts at t ~ tpff and Eq. (33) shows that tpff is as short as approximately a few × 105 yr for rd,0 = 100 au. Even for rd,0 = 300 au, the decay starts at tpff ≃ 106 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 rd,0. This trend is explained by a comparison between the snowline passage time tsnow and tpff, 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 tpff (Eq. (37)). If tsnow > 10 tpff, fwater can be ≪ 1. While tpff is smaller for a smaller rd,0 (Eq. (33)), tsnow is larger (Eq. (40)). The latter implies that the disk is warmer for a smaller rd,0. The viscousheating increases as the blanketing effect by optical depth (∝ Σg) increases. In Fig. 3, g,0 and tdiff are fixed. Smaller rd,0 means smaller α (Eq. (13)) and larger Σg (Eq. (9)), resulting in a warmer disk. Thereby, the final water fraction is lower for smaller rd,0.

We examine the condition of tsnowtpff > 10 or <10 in more detail. In the case of rd,0 = 30 au, tpff ≃ 4 × 104 yr (Eq. (33)) and tsnow is identified by the timing at which fwater starts rapid increase, which is ≃1 × 106, 2.5 × 106, and 5 × 106 yr for the Mars, Earth and Venus analogs, respectively (the left panel of Fig. 3). Because tsnow > 10 tpff in this case, fwater ~ 10−2 even for the outermost Mars analog. For the Earth and Venus analogs, fwater is further smaller. In the case of rd,0 = 300 au, tpff ≃ 1 × 106 yr (Eq. (33)) and tsnow ≃ 1 × 105, 1 × 106 and 2 × 106 yr for the Mars, Earth, and Venus analogs, respectively (right panel of Fig. 3). Because tsnow < 10 tpff, fwater ~ 1∕2 for all of Mars, Earth, and Venus. In the middle panel of Fig. 3, rd,0 = 100 au and tpff ≃ 2 × 105 yr (Eq. (33)). In this case, tsnow < 10 tpff and fwater ≪ 1 for the Venus analog, while tsnow > 10 tpff and fwater ~ 1∕2 for the Mars analog.

The condition of tsnowtpff > 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 tdiff = 106 yr (left panel) and tdiff = 3 × 106 yr (right panel), respectively. The other disk parameters are the same. For both cases, 10 tpff ≃ 2 × 106 yr (Eq. (33)). The snowline is already inside Mars’ orbit from the beginning (t = 0) of the calculations, that is, tsnow = 0, which results in fwater ≃ 1∕2 for the Mars analog. The snowline passage time tsnow is smaller than10 tpff only for the Venus analog in the case of tdiff = 106yr (left panel), and for both the Earth and Venus analogs in the case of tdiff = 3 × 106 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 rd,0 = 100 au, tdiff = 106 yr, and pe = 10−9M yr−1. Again, 10 tpff ≃ 2 × 106 yr for both cases. Smaller g,0 with the fixed tdiff = 106 yr−1 means earlier passage of the snowline, resulting in a lower water fraction. Figure 6 shows the results of pe = 10−9M yr−1 (left) and pe = 3 × 10−9M yr−1 (middle) and pe = 10−8M 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 tsnow ≳ 10 tpff ≃ 2 × 106 yr.

In Sect. 4, we derive a semi-analytical formula to predict the water fraction of planets after disk depletion. We show that the total mass (Mres) 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 orbits2. Since our gas disk model is analytical, we can analytically evaluate tsnow. We already know the analytical expression of tpff given by Eq. (36) is a good approximation. We also semi-analytically derive how the icy grain surface density evolves as Eq. (37). A synthesis of these results enables us to derive the semi-analytical formula for fwater.

Finally we point out that fwater is higher in the order of Mars, Earth, and Venus analogs according to the snowline passage timing, except for the fully saturated cases where fwater ~ 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 fwater.

In the next section, we show that the derived semi-analytical expression of fwater 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 fwater comparable to that of the terrestrial planets in our solar system.

thumbnail Fig. 3

Time evolution of water fraction for the models with tdiff = 106 yr, g,0 = 10−7M yr−1, and pe = 10−9M yr−1. Left, center, and right panel: results of rd.0 = 30 100, and 300 au, respectively. The blue, red, and green lines represent the Earth, Mars and Venus analogs, respectively.

Open with DEXTER
thumbnail Fig. 4

Time evolution of water fraction for the models with rd.0 = 100 au, g,0 = 3 × 10−8M yr−1 and pe = 10−9M yr−1. Left and right panels: results of tdiff = 106 and 3 × 106yr, 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 ΔMpl (Eq. (39)). We can simply estimate it as ΔMplfflt(Mpl) × Mres, where Mres 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 tsnow to infinity (effectively, to tdiff) must be similar to Mres.

Equation (39) is approximated to be (41)

where we use fflt(Mpl) ≃ fflt(Mpl,0) for ΔMplMpl,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 ffltMpl. In that case, Eq. (41) shows that fwater is independent of Mpl both in the cases of ΔMpl < Mpl,0 and ΔMpl > Mpl,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 fwater is still weak: it is proportional to for ΔMpl < Mpl,0 and independentof Mpl otherwise. We can combine the two limits in Eq. (41) by a simple formula as

thumbnail Fig. 5

Time evolution of water fraction for the models with rd.0 = 100 au, tdiff = 106 yr and pe = 10−9 M yr−1. Left and right panels: results of g,0 = 3 × 10−8 and 10−7M yr−1, respectively.

Open with DEXTER
thumbnail Fig. 6

Time evolution of water fraction for the models with rd.0 = 100 au, tdiff = 106 yr and ini = 10−7M yr−1. Left, center, and right panels: results of pe = 10−9, 3 × 10−9, and 10−8M yr−1, respectively.

Open with DEXTER
(42)

In Fig. 7, fwater,sim obtained byour simulations is compared with Eq. (42) with the simulated Mres. This figure showsthat Eq. (42) reproduces the numerical results very well. Both the numerical results and Eq. (42) show fwaterMres for fwater ≪ 1∕2. The icy dust mass Mres depends sensitively on the disk parameters, g,0, pe, tdiff, and rd,0. On the other hand, fflt in 2D is independent of these disk parameters, while that in 3D depends only weakly (). In the case of ΔMplMpl,0, fwater ≃ (ffltMpl,0)Mres. As we discussed, ffltMpl,0 is almost constant, so that fwater should be almost proportional to Mres for ΔMplMpl,0.

Because the analytical solution given by Eq. (42) reproduces the numerical results, we next derive an analytical formula for the icy dust mass Mres at the snowline passage as follows. From Eq. (6), the snowline passes the planet orbit at apl 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 Mres ~ Z Mg,snow, we obtain (47)

where we approximated Z as from Eq. (37), (Eq. (36)), and γ = 1 + 0.15(300 au∕rd,0) (Eq. (38)). Substituting the filtering factor fflt given by Eq. (29) with St = 0.1 and Mres estimated above into Eq. (42), we can analytically estimate the water fraction from the disk parameters, g,0, pe, tdiff, and rd,0 as (Eqs. (29) and (39))

thumbnail Fig. 8

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

Open with DEXTER
(48)

where Mres is given by Eq. (47), fflt = min(fflt,3D, fflt,2D), and

Except for the value of γ, which was fitted by the numerical results, the other derivations are analytical. We note that fwaterMresZ0. Around metal-rich stars, fwater would be larger.

In Fig. 8, we compare the analytically estimated water fraction fwater,anly with the numerically simulated fwater,sim. In addition to the Earth analogs (Mpl,0 = 1M) at 1 au, we also plot the results for the Mars analogs (Mpl,0 = 0.11M) at 1.52 au and the Venus analogs (Mpl,0 = 0.82M) at 0.72 au. For the Earth and Mars analogs, fwater,anly reproduces fwater,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, fwater,anly is larger by a factor of a few to a few tens than fwater,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.

thumbnail Fig. 7

Water fraction fwater,sim and icy dust mass preserved in the outer disk at the snowline passage Mres obtained by our simulations. The blue squares show the simulation results for the Earth analogs with different values of g,0, pe, tdiff, and rd,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 (fwater) is determined by the disk parameters. Figure 9 shows fwater for a planet at apl = 1.0 au as a function of the disk parameters, g,0 and rd,0. The other parameters are pe = 10−9M yr−1 and tdiff = 3 × 106 yr−1. The panels a show the dependence on Mpl,0 for St = 0.1. The panels b show the dependence on St for Mpl,0 = 1M. The planetsformed with the parameters in the red region are very dry, fwater ≲ 10−4. Those in the green region have a modest amount of water, fwater ~ 0.1, and those in the blue region are icy planets, fwater ≃ 1∕2. The yellow and orange regions represent fwater ~ 10−4−10−2, which corresponds to the water fraction of the current Earth and that estimated for ancient Mars.

We find that fwater is the most sensitive to g,0 and rd,0. The common features in the contours in Fig. 9 are that (i) fwater is lower for smaller rd,0 and larger g,0, and (ii) fwater ~ 10−4−10−2 is realized at rd,0 ~ 30−50 au and g,0 ≳ 10−8M yr−1. This is consistent with the conclusion by Sato et al. (2016) that such fwater 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 fwater ~ 10−4−10−2 is not in a very narrow window.

As we showed in Sect. 4, fwater is regulated by Mres and Mres is sensitive to g,0 and rd,0, especially through the parameter tsnowtpff (Eq. (47)). Because pebble accretion is fast to realize fwater ~ 10−4−10−2, Mres has to have significantly decayed until t ~ tsnow. For small rd,0, tpff is small (Eq. (36)), while tsnow is large (Eq. (45)) due to small g,snow (Eq. (44)). For large g,0, tsnow is small (Eq. (44)) while tpff is the same (Eq. (36)). As a result, Σp decays more quickly (Eq. (37)) for small rd,0 and large g,0.

Figure 9 shows the dependences of fwater on the pebble accretion parameters: the initial planetary mass (Mpl,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 > tpff. 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 fwater is almost independent of Mpl,0 and St. The weak dependence on Mpl,0 is explained by the following arguments. The water fraction fwater is ~ (ffltMpl,0)Mres in the case of ΔMplMpl,0. The dust mass Mres is independent of Mpl. As already explained in Sect. 4, the factor ffltMpl,0 is independent of Mpl,0 in 3D accretion and only weakly depends on Mpl,0 in 2D accretion ().

The weak dependence of fwater on Stokes number shown in Fig. 9 is the result of the assumption that Mres is independent of St; it depends on St only weakly (Eq. (29)) through fflt. The parameter γ could be dependent on St because the sculpture rate of Σp at r = rpff may depend on St. Figure 10 shows the dependence on γ2, assuming that the functional form of γ = 1 + γ2(300 au∕rd,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 fwater 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, tdiff and pe. For larger tdiff, 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, Mres is smaller (Eq. (47)) because tsnowtdiff and γ ≥ 1. For smaller pe, the disk is hotter and the snowline passage is later (tsnow is larger), resulting in smaller Mres. If Mres is smaller, fwater 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 rd,0 and lower g,0 for smaller tdiff and larger pe, because Mres 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 fwater is plotted for apl = 0.72, 1.00 and 1.52 au (the Venus, Earth, and Mars analogs, respectively). In the outer region, fwater is generally larger due to early snowline passage (smaller tsnow). In other words, g,snow is larger for larger apl. As mentioned in Sect. 4, the Venus analogs may have further lower fwater, 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 rd,0 ~ 30−50 au and g,0 ≳ 2 × 10−8 M yr−1, fwater ~ 10−4−10−2 for both the Earth and Mars analogs.

In Sect. 3, we pointed out that the simple condition of tsnowtpff < 10 or > 10 discriminates between the water-rich case (fwater ~ 1∕2) and the water-poor case. The pebble formation timescale tpff is given by Eq. (36) and the snowline passage time tsnow is given by Eq. (45). Both tpff and tsnow are independent of Mpl and St, which is consistent with Fig. 9. The condition of tsnow = 10 tpff is shown on the rd,0g,0 plane in Fig. 13. We show the dependences on tdiff, pe, and apl. A comparison of this figure with Figs. 11 and 12 show that the simple condition approximately reproduces the more detailed evaluation except for rd,0 > 100 au. As shown in Fig. 1, for rd,0 > 100 au, the disk radius expands with the pebble formation front radius rpff, so that the sculpture of icy dust reservoir is delayed, compared to the estimation of tpff at rd,0. The calculation of Mres does not have this problem.

As we have shown, fwater is the most sensitive to g,0 and rd,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, fwater ~ 10−4−10−2, is realized at rd,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.

thumbnail Fig. 9

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,0) and the initial disk accretion rate (g,0). Panels a: dependence on the planet mass (Mpl,0). Panels b: Stokes number (St) of pebbles. The other parameters are pe = 10−9M yr−1 and tdiff = 3 × 106 yr. In the panelsa, we use St = 0.1. In the panelsb, Mpl,0 = 1M. The middle panels in a and b are identical. The color scales are log10fwater.

Open with DEXTER
thumbnail Fig. 10

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,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, Mpl,0 = 1M, , and tdiff = 3 × 106 yr. The color scales are log10fwater.

Open with DEXTER
thumbnail Fig. 11

Analytical estimate of fwater at 1 au as a function of rd,0 and g,0. Panels a: The dependence on the disk diffusion timescale tdiff. Panels b: dependence on the photoevaporation rate (pe). The other parameters are St = 0.1 and Mpl,0 = 1M. In the panels a, we used pe = 10−9M yr−1, and tdiff = 3 × 106 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 log10 fwater.

Open with DEXTER
thumbnail Fig. 12

Analytical estimate of fwater as a functionof rd,0 and g,0, at 0.72, 1.00, and 1.52 au. The other parameters are St = 0.1, Mpl,0 = 1M, pe =10−9M yr−1, and tdiff = 3 × 106 yr. The color scales are log10fwater.

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 (Mpl) 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” (Mpeb,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 Mpl > Mpeb,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, Mpeb,iso ≃ 4.1M for tdiff = 3 × 106 yr, g,0 = 10−8M yr−1, rd,0 au, hr ≃ 0.0246 (at 1 au), and α ≃ 3 × 10−3. If Mpl,0 ≳several M, the effect of pebble isolation mass is not negligible. Figure 14 shows the water fraction for Mpl,0 = 3, 5, and 10 M. We stop pebble accretion when Mpl reaches Mpeb,iso. Because α is smaller for smaller rd,0 (Eq. (14)) and hr is lower for smaller g,0 and higher α (Eq. (3)), Mpeb,iso is smaller for smaller rd,0 and g,0. The dry regime (red-colored regime) in the left bottom part of the plots represent the cases of Mpl,0 > Mpeb,iso. Because tsnow is small in the low g,0 regions, fwater rapidly increases because of a high pebble flux until Mpl increases to Mpeb,iso. Thereby, the edge in fwater at Mpl,0Mpeb,iso is sharp.

For the same tdiff, g,0, and rd,0 as the above, . The filtering rate for Mpeb,iso is in 2D case (Eq. (29)). If ice giants or cores of gas giants are formed, even before their mass exceeds Mpeb,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 Mres ≳ 10 M. From Eq. (47) with Md,0 ~ 10−8 M yr−1 and tdiff ~ 3 × 106 yr, , which is for rd,0 ~ 100 au. Hence, the Jupiter’s core formation time tjup must be ≲ tdiff. 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 fwater rapidly increases after t = tsnow and becomes saturated before t ~ tdiff. If tsnow > tjup, fwater = 0 for terrestrial planets. Otherwise, it is likely that fwater is already close to the saturated value because the increase of fwater 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 tsnow ~ 10 tpff 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 cosmo-chemical data.

thumbnail Fig. 13

Condition of tsnow = 10 tpff as a functionof rd,0 and g,0. Panels a: dependence on the disk diffusion timescale tdiff. Panels b: photoevaporation rate (pe). Panels c: planetary orbital radius (apl). The right regions from the curves represent water-rich regions (fwater ~ 1∕2). In the panelsa, we used pe = 10−9M yr−1 and apl= 1 au, tdiff = 3 × 106 yr and apl= 1 au in the panelsb, and pe = 10−9M yr−1 and tdiff = 3 × 106 yr in the panels c.

Open with DEXTER
thumbnail Fig. 14

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,0) and the initial disk accretion rate (g,0), including the effect of the pebble isolation mass. The other parameters are pe =10−9M yr−1, tdiff = 3 × 106 yr, and St = 0.1. Pebble accretion is halted when Mpl reaches Mpeb,iso. The dry regime (red) in the left bottom part of the plots represents the cases of Mpl,0 > Mpeb,iso.

Open with DEXTER

7 Summary

If water is not delivered to rocky planets in HZs, the planets cannot be actual habitats because H2 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 solar-type 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 = tsnow) 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 self-similar solution with constant viscous α (Lynden-Bell & 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 (fwater), the snowline passage time (tsnow) relative to the time (tpff) 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 (Mres) preserved in the disk outer region at t = tsnow determines the water fraction of the final planets. The accreted ice mass to the planet is estimated by ~ (1∕2)ffltMres, where the filtering factor fflt is the fraction of the pebble mass flux that is accreted onto the planet. Because Mres rapidly decreases after t ~ tpff, tsnowtpff > 10 or < 10 is crucial for final value of fwater. If tsnowtpff > 10, Mres should have significantly decayed when icy pebble accretion starts at t = tsnow.

Using these numerical results, we derived an analytical formula for fwater by icy pebble accretion. In the formula, fwater is explicitly given as a function of the ratio tsnowtpff and the disk parameters. The parameter tsnowtpff is also determined by the disk parameters. As a result, fwater is predictedby the disk parameters, especially the initial disk mass accretion rate g,0 and initial disk size rd,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 fwater ~ 10−4−10−2, which may correspond to the value of the current Earth, in disks with an initial disk size of rd,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−8M yr−1, both the Earth and Mars analogs have fwater ~ 10−4−10−2, while fwater 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 solar-type 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

  1. Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bercovici, D., & Karato, S.-i. 2003, Nature, 425, 39 [NASA ADS] [CrossRef] [Google Scholar]
  5. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Brown, H. 1949, in The Atmospheres of the Earth and Planets, ed. G. P. Kuiper (Chicago: University of Chicago), 258 [Google Scholar]
  9. Clifford, S. M., Lasue, J., Heggy, E., et al. 2010, J. Geophys. Res. Planets, 115, E07001 [NASA ADS] [CrossRef] [Google Scholar]
  10. di Achille, G., & Hynek, B. M. 2010, Lunar Planet. Sci. Conf., 41, 2366 [NASA ADS] [Google Scholar]
  11. Donahue, T. M., Hoffman, J. H., Hodges, R. R., & Watson, A. J. 1982, Science, 216, 630 [NASA ADS] [CrossRef] [Google Scholar]
  12. Fei, H., Yamazaki, D., Sakurai, M., et al. 2017, Sci. Adv., 3, e1603024 [NASA ADS] [CrossRef] [Google Scholar]
  13. Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [NASA ADS] [CrossRef] [Google Scholar]
  14. Genda, H., & Abe, Y. 2005, Nature, 433, 842 [NASA ADS] [CrossRef] [Google Scholar]
  15. 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]
  16. Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [NASA ADS] [CrossRef] [Google Scholar]
  18. Haisch, Jr. K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
  19. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
  20. Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31 [NASA ADS] [CrossRef] [Google Scholar]
  22. Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hirschmann, M. M. 2006, Ann. Rev. Earth Planet. Sci., 34, 629 [NASA ADS] [CrossRef] [Google Scholar]
  24. Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Johansen, A., Ida, S., & Brasser, R. 2019, A&A, 622, A202 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Klahr, H., Pfeil, T., & Schreiber, A. 2018, Handbook of Exoplanets (London: Springer), 138 [Google Scholar]
  27. Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Nat. Acad. Sci., 114, 6712 [NASA ADS] [Google Scholar]
  28. Kurokawa, H., Sato, M., Ushioda, M., et al. 2014, Earth Planet. Sci. Lett., 394, 179 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Lunine, J. I., Chambers, J., Morbidelli, A., & Leshin, L. A. 2003, Icarus, 165, 1 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lyra, W., & Umurhan, O. 2018, PASP, submitted [Google Scholar]
  35. Machida, R., & Abe, Y. 2010, ApJ, 716, 1252 [NASA ADS] [CrossRef] [Google Scholar]
  36. Marty, B. 2012, Earth Planet. Sci. Lett., 313, 56 [NASA ADS] [CrossRef] [Google Scholar]
  37. Matsumura, S., Brasser, R., & Ida, S. 2016, ApJ, 818, 15 [NASA ADS] [CrossRef] [Google Scholar]
  38. Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [NASA ADS] [CrossRef] [Google Scholar]
  40. Morbidelli, A., Chambers, J., Lunine, J. I., et al. 2000, Meteorit. Planet. Sci., 35, 1309 [NASA ADS] [CrossRef] [Google Scholar]
  41. Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mori, S., Bai, X.-N., & Okuzumi, S. 2019, ApJ, 872, 98 [NASA ADS] [CrossRef] [Google Scholar]
  43. Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 818, 16 [NASA ADS] [CrossRef] [Google Scholar]
  44. Nomura, R., Hirose, K., Uesugi, K., et al. 2014, Science, 343, 522 [NASA ADS] [CrossRef] [Google Scholar]
  45. O’Brien, D. P., Walsh, K. J., Morbidelli, A., Raymond, S. N., & Mandell, A. M. 2014, Icarus, 239, 74 [NASA ADS] [CrossRef] [Google Scholar]
  46. Oka, A., Nakamoto, T., & Ida, S. 2011, ApJ, 738, 141 [NASA ADS] [CrossRef] [Google Scholar]
  47. Okuzumi, S., & Tazaki, R. 2019, ApJ, submitted [Google Scholar]
  48. Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
  49. Ormel, C. W. 2014, ApJ, 789, L18 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Ormel, C. W., & Kobayashi, H. 2012, ApJ, 747, 115 [NASA ADS] [CrossRef] [Google Scholar]
  52. Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icarus, 168, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  53. Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  55. Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Takeuchi, T., & Lin, D. N. C. 2005, ApJ, 623, 482 [NASA ADS] [CrossRef] [Google Scholar]
  57. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
  58. Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
  59. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]

1

The wind-driven disk accretion that has recently been proposed (e.g., Suzuki et al. 2016; Bai et al. 2016) is commented on in Sect. 2.1.

2

The grains at r > 300−500 au would not undergo radial drift sufficiently, because their growth timescale is longer than tdiff. We consider disks with rd,0 ≤ 300 au.

All Figures

thumbnail Fig. 1

Time evolution of the transition radius rvis/irr (Eq. (5); red lines), snowline rsnow = max(rsnow,vis, rsnow,irr) (Eqs. (6) and (7); blue lines), disk characteristic radius (green lines), and pebble formation front rpff (Eq. (36); magenta lines) for various disk evolution parameters. The parameters, rd,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 tdiff = 106 yr.

Open with DEXTER
In the text
thumbnail Fig. 2

Time evolution of Z = Σp∕Σg obtained by the numerical simulation with g,0 = 10−7M yr−1, pe =10−9M yr−1, and tdiff = 106 yr. The magenta, blue, and green solid lines show Z obtained at rd,0 by the numerical results of with rd,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
thumbnail Fig. 3

Time evolution of water fraction for the models with tdiff = 106 yr, g,0 = 10−7M yr−1, and pe = 10−9M yr−1. Left, center, and right panel: results of rd.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
thumbnail Fig. 4

Time evolution of water fraction for the models with rd.0 = 100 au, g,0 = 3 × 10−8M yr−1 and pe = 10−9M yr−1. Left and right panels: results of tdiff = 106 and 3 × 106yr, respectively.

Open with DEXTER
In the text
thumbnail Fig. 5

Time evolution of water fraction for the models with rd.0 = 100 au, tdiff = 106 yr and pe = 10−9 M yr−1. Left and right panels: results of g,0 = 3 × 10−8 and 10−7M yr−1, respectively.

Open with DEXTER
In the text
thumbnail Fig. 6

Time evolution of water fraction for the models with rd.0 = 100 au, tdiff = 106 yr and ini = 10−7M yr−1. Left, center, and right panels: results of pe = 10−9, 3 × 10−9, and 10−8M yr−1, respectively.

Open with DEXTER
In the text
thumbnail Fig. 7

Water fraction fwater,sim and icy dust mass preserved in the outer disk at the snowline passage Mres obtained by our simulations. The blue squares show the simulation results for the Earth analogs with different values of g,0, pe, tdiff, and rd,0 (Sect. 2.4). The dotted curve represents the analytical solution given by Eq. (42).

Open with DEXTER
In the text
thumbnail Fig. 8

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

Open with DEXTER
In the text
thumbnail Fig. 9

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,0) and the initial disk accretion rate (g,0). Panels a: dependence on the planet mass (Mpl,0). Panels b: Stokes number (St) of pebbles. The other parameters are pe = 10−9M yr−1 and tdiff = 3 × 106 yr. In the panelsa, we use St = 0.1. In the panelsb, Mpl,0 = 1M. The middle panels in a and b are identical. The color scales are log10fwater.

Open with DEXTER
In the text
thumbnail Fig. 10

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,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, Mpl,0 = 1M, , and tdiff = 3 × 106 yr. The color scales are log10fwater.

Open with DEXTER
In the text
thumbnail Fig. 11

Analytical estimate of fwater at 1 au as a function of rd,0 and g,0. Panels a: The dependence on the disk diffusion timescale tdiff. Panels b: dependence on the photoevaporation rate (pe). The other parameters are St = 0.1 and Mpl,0 = 1M. In the panels a, we used pe = 10−9M yr−1, and tdiff = 3 × 106 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 log10 fwater.

Open with DEXTER
In the text
thumbnail Fig. 12

Analytical estimate of fwater as a functionof rd,0 and g,0, at 0.72, 1.00, and 1.52 au. The other parameters are St = 0.1, Mpl,0 = 1M, pe =10−9M yr−1, and tdiff = 3 × 106 yr. The color scales are log10fwater.

Open with DEXTER
In the text
thumbnail Fig. 13

Condition of tsnow = 10 tpff as a functionof rd,0 and g,0. Panels a: dependence on the disk diffusion timescale tdiff. Panels b: photoevaporation rate (pe). Panels c: planetary orbital radius (apl). The right regions from the curves represent water-rich regions (fwater ~ 1∕2). In the panelsa, we used pe = 10−9M yr−1 and apl= 1 au, tdiff = 3 × 106 yr and apl= 1 au in the panelsb, and pe = 10−9M yr−1 and tdiff = 3 × 106 yr in the panels c.

Open with DEXTER
In the text
thumbnail Fig. 14

Analytical estimate of fwater at 1 au as a function of the initial disk radius (rd,0) and the initial disk accretion rate (g,0), including the effect of the pebble isolation mass. The other parameters are pe =10−9M yr−1, tdiff = 3 × 106 yr, and St = 0.1. Pebble accretion is halted when Mpl reaches Mpeb,iso. The dry regime (red) in the left bottom part of the plots represents the cases of Mpl,0 > Mpeb,iso.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.