THE SNOW LINE IN VISCOUS DISKS AROUND LOW-MASS STARS: IMPLICATIONS FOR WATER DELIVERY TO TERRESTRIAL PLANETS IN THE HABITABLE ZONE

The leading explanation for the delivery of water to Earth is that water-bearing, asteroid-like bodies from beyond the snow line were gravitationally scattered inward and accreted by the planet during its growth (e.g., Morbidelli et al. 2000). If the same mechanism also operates in extrasolar planetary systems, the water content of potentially habitable terrestrial planets that formed in situ depends on both the snow-line location and the extent to which water-bearing materials were scattered inward. For lower-mass M dwarf stars, N-body simulations indicated that terrestrial planets in the habitable zone would be dry (Lissauer 2007; Raymond et al. 2007), in contrast to planets that formed farther out in the disk and migrated to their current locations (Ogihara & Ida 2009). In a previous paper (Ciesla et al. 2015), we explored how the inclusion of planetesimals and more comet-like bodies in these simulations can enhance volatile delivery around these low-mass stars. In that work, the location of the snow line was a critical factor in determining how much water would be delivered to forming planets. Here, we revisit the location of the snow line based on recent protoplanetary disk observations and explore its impact on water delivery to habitable-zone terrestrial planets.

In the solar system, the location of the snow line at the time of planetesimal formation has been inferred from the transition between hydrous and anhydrous asteroids to be ∼2.5 AU (Abe et al. 2000). There is, however, a considerable uncertainty in this location if asteroids have been scattered into their current locations (e.g., Walsh et al. 2012; DeMeo & Carry 2014). There are no direct measurements of the location of a water snow line outside our solar system; hence, we do not know if it is typical for Sun-like stars. The solar system snow-line location of 2.5 AU has been taken as a reference point in previous studies and used to estimate the location around other stars of various masses (e.g., Raymond et al. 2007). The exact location of the snow line will be determined by the combined irradiative and viscous heating in the disk (e.g., Davis 2005). If mass accretion occurs through the vertical extent of the disk, the midplane temperatures at the relevant locations are set by the release of gravitational potential energy for mass accretion rates $\gt {10}^{-10}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}\;$, and the location of the snow line is a strong function of mass accretion rate (Garaud & Lin 2007; Oka et al. 2011).

The dependence of mass accretion rate on stellar mass has been well established, following a roughly quadratic relationship (e.g., Muzerolle et al. 2003, 2005; Calvet et al. 2004; Natta et al. 2006), but these relations have not been used to calculate the location of the snow line around low-mass stars. On top of that, there is an intrinsic dispersion in mass accretion rates of ∼0.4 dex around this relation (Alcalá et al. 2014). If this intrinsic dispersions reflects a range in initial disk masses, which is of the same order, 0.5 dex (Armitage et al. 2003), we expect an associated dispersion in snow-line locations (Figure 1), leading to enhanced or reduced water delivery among stars of similar mass.

Zoom In Zoom Out Reset image size

In addition to the mass accretion rate, the location of the snow line in a protoplanetary disk around a Sun-like star is very sensitive to the dust opacity (Min et al. 2011; Oka et al. 2011). While previous work on the snow-line location in the solar nebula had mainly used small sub-micron-sized grains (e.g., Davis 2005; Garaud & Lin 2007), there is abundant observational and theoretical evidence that grains in protoplanetary disks quickly grow to millimeter and centimeter sizes (e.g., Testi et al. 2014), depleting the amount of small grains by orders of magnitude (e.g., Furlan et al. 2006). Lower dust opacities trap viscous heat less efficiently, resulting in a cooler disk midplane with a snow line closer to the star, enhancing water delivery.

In this work, we explore the range of snow-line locations that may be expected for stars of different masses, factoring in variations in disk mass accretion rate and dust opacities that have been inferred for real disks (Section 2) and how this would impact the water content of terrestrial planets in the habitable zone (Section 3).

The location of the snow line, ${R}_{\mathrm{SL}\;}$, in a steady-state viscous protoplanetary disk, assuming that it is optically thick to its own radiation, is given by Min et al. (2011), their Equation (11):

$$\begin{eqnarray}&&{R}_{\mathrm{SL}\;}={\left(\displaystyle \frac{3\mu {m}_{p}{({{GM}}_{\star })}^{3/2}{\dot{M}}^{2}{\kappa }_{R}}{128{\pi }^{2}{k}_{b}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{ice}}f\alpha }\right)}^{2/9},\end{eqnarray} \tag{ 1 }$$

where μ is the mean molecular weight, m_p the proton mass, G the gravitational constant, $\dot{M}$ the mass accretion rate, ${M}_{\star }$ the stellar mass, ${\kappa }_{R}$ the Rosseland mean opacity, k_b the Boltzmann constant, ${\sigma }_{\mathrm{SB}}$ the Stefan–Boltzmann constant, ${T}_{\mathrm{ice}}$ the temperature where water-ice condenses, f the gas-to-dust ratio, and α the turbulent mixing strength. This equation is based on the estimate of the midplane temperature of a viscous disk by Hubeny (1990). Min et al. (2011) showed that this equation is in good agreement with the snow line computed through detailed radiative transfer modeling for a Sun-like star. Even though the role of irradiation becomes larger when using the pre-main-sequence luminosity and when considering lower-mass disks around lower-mass stars, we have verified in the Appendix that viscous heating indeed dominates the thermal budget of the disk midplane for the range of stellar masses and disk mass accretion rates explored in this paper.

The ice sublimation temperature is thought to lie between 150 and 170 K (e.g., Podolak & Zucker 2004), depending on the water vapor pressure in the disk. The radiative transfer models in Min et al. (2011) include a detailed treatment of solid sublimation as described in Kama et al. (2009). Min et al. (2011) show that, for a large range of mass accretion rates—and hence surface densities and partial vapor pressures—the location of the snow line is well predicted by assuming a single sublimation temperature of ${T}_{\mathrm{ice}}=160\;{\rm{K}}$. We therefore do not solve the radial diffusion equation (e.g., Ciesla & Cuzzi 2006), as we expect the variations of the partial water pressure by inward drift of evaporating ices to be within the tested range and not lead to large variations in the location of the snow line.

Grouping all constants together, taking μ = 2.3, we obtain the more manageable expression

$$\begin{eqnarray}{R}_{\mathrm{SL}\;} & = & 2.1\;\mathrm{AU}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{\dot{M}}{{10}^{-8}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}\;}\right)}^{4/9}\\ & & \times {\left(\displaystyle \frac{{\kappa }_{R}}{770\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}}\right)}^{2/9}{\left(\displaystyle \frac{f}{100}\right)}^{-2/9}\\ & & \times {\left(\displaystyle \frac{\alpha }{0.01}\right)}^{-2/9}{\left(\displaystyle \frac{{T}_{\mathrm{ice}}}{160\;{\rm{K}}}\right)}^{-10/9}.\end{eqnarray} \tag{ 2 }$$

In the remainder of the section we explain the origin and applicability of the parameters in Equation (2) that lead to the range of snow lines shown in Figure 1.

2.1. Disk Mass Accretion Rate

The mass accretion rate scales roughly with stellar mass squared (e.g., Muzerolle et al. 2003, 2005; Calvet et al. 2004; Natta et al. 2006), with an observed dispersion of up to an order of magnitude around this relation. We will use the observed mass accretion rates from Alcalá et al. (2014) since these have the smallest observed dispersion, and we will explain in Section 2.2 how we treat this dispersion.

The median mass accretion rate for the ∼3 ± 1 Myr old (Comerón 2008) Lupus star-forming region, over a range of stellar masses from $0.03\;{M}_{\odot }$ to $1.0\;{M}_{\odot }$, is given by Alcalá et al. (2014):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Lupus}}({M}_{\star })=5.6\cdot {10}^{-9}{\left[\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right]}^{1.81}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}\;.\end{eqnarray} \tag{ 3 }$$

We derive a generic expression for the median mass accretion rate as a function of time by appending this equation with a factor ${t}^{-3/2}$ to take into account the age of the cluster:

$$\begin{eqnarray}&&{\dot{M}}_{m}\;({M}_{\star },t)=3\cdot {10}^{-8}{\left[\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right]}^{1.81}{\left[\displaystyle \frac{t}{\mathrm{Myr}\;}\right]}^{-3/2}\;{M}_{\odot }\;{\mathrm{yr}}^{-1}\;.\end{eqnarray} \tag{ 4 }$$

This accretion rate is consistent with the derived accretion rate for the 1 Myr old Taurus region of $\sim {10}^{-8}{M}_{\odot }\;{\mathrm{yr}}^{-1}\;$ for $\sim 0.5\;{M}_{\odot }$ stars (Hartmann et al. 1998) and matches well with the observed decay of the mass accretion rate between 0.3 and 30 Myr (Sicilia-Aguilar et al. 2010, their Figure 2). The exponent of $\eta =1.5$ comes from a viscously evolving disk with fixed temperature (Hartmann et al. 2006) and matches well with the analytical model of Chambers (2009) that includes a changing disk temperature (η = 20/13 = 1.54), which is more similar to the steady-state accretion disk models employed in deriving Equation (1) by Min et al. (2011). Equation (4) is only valid at timescales much longer than the viscous timescale and starts significantly overpredicting the mass accretion rate at $t\lt 0.3$ Myr for α = 0.01 compared to the Chambers (2009) disk model.

2.2. Range of Snow-line Locations

A dispersion in disk mass accretion rate at a given age will lead to a dispersion in snow-line locations. Observed mass accretion rates show a large dispersion at any given stellar mass. Part of this dispersion may be attributed to variations in the accretion flow onto the star (Costigan et al. 2012) and use of secondary tracers of accretion such as emission lines (Rigliaco et al. 2012). An intrinsic dispersion may reflect an age spread and variations in the initial disk mass (e.g., Hartmann et al. 1998). The tightest constraint on the intrinsic dispersion currently comes from Alcalá et al. (2014). Using the X-shooter spectrograph, the authors compute the mass accretion rates using the UV excess emission, a direct tracer of accretion, in the Lupus star-forming region over a range of stellar masses from $0.03\;{M}_{\odot }$ to $1.0\;{M}_{\odot }$. The dispersion in mass accretion rates at a given stellar mass is characterized by a lognormal distribution with a standard deviation ${\sigma }_{\dot{M}}=0.4$ dex, a factor of two lower than in previous works. We assume that this dispersion reflects a range of initial disk masses as proposed by Hartmann et al. (1998).

This assumption is supported by two results. First, the diversity in protoplanetary disk properties, in particular the outer radius-disk mass relation (Andrews et al. 2010), and the exponential decay of the disk fraction versus time (e.g., Mamajek 2009) argue for an intrinsic spread in initial disk masses. Since the inferred spread in disk masses from the time-dependent disk fraction is of the same order, ∼0.5 dex (Armitage et al. 2003), we take the observed dispersion in mass accretion rate to be the intrinsic one. Second, although accretion rates measured from spectral lines are known to be variable on timescales of days to years, this variation is smaller than the observed dispersion and driven by the rotational modulation of the accretion flow onto the star (Costigan et al. 2014) and does not reflect an intrinsic variation in the disk mass accretion rate. It is possible that disks show larger variations at timescales of decades or longer but smaller in magnitude than FUor outbursts, or that it reflects a large age spread in the Lupus cloud, in which case the intrinsic scatter could be smaller.

With the assumption that the observed spread in mass accretion rates is the intrinsic one, the probability distribution of the snow line R around the median, ${R}_{\mathrm{SL}\;}$ (Equation (2)), is given by a lognormal distribution:

$$\begin{eqnarray}&&{f}_{\mathrm{SL}}(R,{R}_{\mathrm{SL}\;})=\displaystyle \frac{1}{\sqrt{2\pi }R{\sigma }_{R}}\mathrm{exp}\left(-\displaystyle \frac{{(\mathrm{log}R-\mathrm{log}{R}_{\mathrm{SL}\;})}^{2}}{2{\sigma }_{R}^{2}}\right)\end{eqnarray} \tag{ 5 }$$

with a standard deviation ${\sigma }_{R}=4/9*{\sigma }_{\dot{M}}\approx 0.18$ dex. This distribution is shown as the shaded regions in Figure 1.

2.3. Dust Opacity

Since a major source of uncertainty in the calculated location of the snow line is the dust opacity (Min et al. 2011), we explore two limiting cases: the first is where the grains are ISM-like in size as in most previous work, with a Rosseland mean opacity at temperatures relevant for water-ice condensation of ${\kappa }_{R}=570\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}$, as in Min et al. (2011); the second is where dust growth is assumed to occur, resulting in a size distribution from micron- to centimeter-sized grains, yielding an opacity of ${\kappa }_{R}=20\;{\mathrm{cm}}^{2}\;{{\rm{g}}}^{-1}$. Such a lower opacity will also increase the amount of solids available for planetesimal formation interior to the snow line; see Min et al. (2011) for details. Because the dependence of the snow-line location on this opacity is relatively weak (${R}_{\mathrm{SL}\;}\propto {\kappa }_{R}^{2/9}$, Equation (2)), the snow line moves in by about a factor of two between the two cases considered here (top and bottom panels of Figure 1). Table 1 lists the mass accretion and corresponding snow-line locations for the range of stellar masses for the two opacities considered here. Figure 1 shows the median snow lines for the same stellar mass range, compared to the recent habitable-zone estimate from Kopparapu et al. (2013). Throughout this paper, we use a gas-to-dust ratio f = 66 based on a condensation sequence at solar metallicity, following Min et al. (2011).

Table 1.  Mass Accretion Rates and Location of the Viscous Snow Line as a Function of Stellar Mass Used in This Work, Including 1σ Ranges

M*	$\mathrm{log}\dot{M}$	RSL (1σ)		RSL
		Small Grains	Large Grains	MS	Pre-MS
(${M}_{\odot }$)	($\mathrm{log}{M}_{\odot }\;{\mathrm{yr}}^{-1}$)	(AU)	(AU)	(AU)	(AU)
1.0	−7.5 ± 0.4	3.3 (2.2...5.0)	1.6 (1.1...2.4)	2.5	1.35
0.8	−7.7 ± 0.4	2.6 (1.7...3.9)	1.2 (0.82...1.9)	1.23	1.15
0.6	−7.9 ± 0.4	1.9 (1.2...2.8)	0.9 (0.59...1.3)	0.61	0.95
0.4	−8.3 ± 0.4	1.2 (0.78...1.8)	0.6 (0.37...0.84)	0.32	0.8
0.2	−8.8 ± 0.4	0.5 (0.36...0.81)	0.2 (0.17...0.38)	0.16	0.5

Note. For comparison, the last two columns show the snow line locations calculated from the stellar luminosity used in previous work: main sequence (MS; Raymond et al. 2007) and pre-MS (Ciesla et al. 2015).

Download table as:  ASCIITypeset image

2.4. Snow-line Location

2.5. Disk Viscosity and Time of Planetesimal Formation

Throughout this paper, we assume $t=1\;\mathrm{Myr}\;$ and $\alpha =0.01$ for setting the location of the snow line. Both quantities are not well constrained observationally. The strength of viscosity in protoplanetary disks is inferred to be of order α ∼ 0.01 based on disk masses, sizes, and lifetimes (Hartmann et al. 1998). The time that defines the location of the snow line depends on both the planet formation timescale and the extent to which (icy) grains continue to be incorporated into larger bodies as the snow line moves in. The theoretical planetesimal formation timescale at ∼1 AU is of order 0.01–0.1 Myr (e.g., Wetherill & Stewart 1993), and the earliest differentiated bodies formed within $\lt 1\;\mathrm{Myr}$ (Kleine et al. 2009), suggesting an early epoch of planetesimal formation. On the other hand, the age of the chondrules in chondritic parent bodies indicates that planetesimal formation continued to take place over a 3 Myr timescale (e.g., Johansen et al. 2014, p. 547). If water vapor remains present in the disk, an inward-moving snow line may give rise to rapid but late icy planetesimal formation (Kennedy et al. 2006). The age of $t=1\;\mathrm{Myr}\;$ is a compromise between two extremes.

To give the reader an idea of how these uncertainties impact the location of the snow line, Figure 2 shows how the location of the snow line varies around a Sun-like star with a median mass accretion rate following Equation (2). Earlier times and lower turbulent mixing strength correspond to higher surface densities and mass accretion rates that lead to a hotter midplane and a more distant snow line. Higher mass accretion rates and later times have the opposite effect, leading to a colder, optically thinner disk with a closer-in snow line. In this regime, irradiation may become important as a heating source in the midplane.

Zoom In Zoom Out Reset image size

To explore how the range in snow-line locations may impact the water content of terrestrial planets in the habitable zone, we use the set of N-body simulations from Ciesla et al. (2015). This approach follows that of O'Brien et al. (2006), Raymond et al. (2007), and Lissauer (2007) for forming terrestrial planets via giant impacts as in the solar system. The simulations are started from a $5\;{M}_{\oplus }$ disk with equal-mass planetary embryos (1/20 Earth mass) and planetesimals (1/20 embryo mass), spaced such that they define a 1/R power law in surface density, as in O'Brien et al. (2006), between 0.5 and 4.0 AU. We do not include giant planets, but note that they may be rare around low-mass stars (Johnson et al. 2010) and water delivery can take place in their absence (Quintana & Lissauer 2014). Four simulations are run with near-identical initial conditions to take into account stochastics. These are repeated for lower stellar masses (0.2, 0.4, 0.6, and 0.8 ${M}_{\odot }$) by scaling all masses (total, embryo, and planetesimal) to the stellar mass and adjusting the inner and outer edge to encompass both snow line and habitable zone. After 200 Myr, a number of terrestrial planets has formed in the inner regions; see Ciesla et al. (2015) for details.

The water content of these planets depends on the composition of accreted planetary embryos and planetesimals, which in turn depends on their starting locations. We assume that the division between icy and rocky building blocks is located at the snow line, i.e., that planetesimals and embryos do not move radially during their growth process, and that the transition from pebbles to planetesimals and embryos is short compared to the timescale at which the snow line moves inward (e.g., Carrera et al. 2015). Assuming that the water content of icy building blocks is low (5%, similar to carbonaceous chondrites), the extra mass outside of the snow line does not affect the dynamics of the system, allowing the snow line to be inserted a posteriori. By keeping track of the starting locations of planetary building blocks that comprise the final planets, their water content can be calculated for different snow-line locations from a single simulation. Figure 3 gives an example of this approach, which shows the water content of planets in the $0.6\;{M}_{\odot }$ simulations calculated from three different locations of the snow line.

Zoom In Zoom Out Reset image size

By calculating the water content of planets in the habitable zone for a range of snow-line locations from the simulations, and assigning to each snow line a probability based on the observed dispersion in mass accretion rates (Equation (5); see also Figure 1), we calculate the probability distribution of the water fraction of planets in the habitable zone. These are shown as cumulative probability distributions in Figure 4.

Zoom In Zoom Out Reset image size

We take a water mass fraction of 0.1% as a reference point for early Earth, between 0.03% and 0.1% currently in the mantle (Lécuyer 2013) and as much as 1% in the past according to Abe et al. (2000). In the simulations with Sun-like stars, for a snow line set by small grains, 20%–40% of planets in the habitable zone remain dry, as a fraction of stars have a mass accretion rate high enough to put the snow line too far away from the habitable zone (4–5 AU) to allow water delivery via giant impacts. The fraction of dry planets is higher in simulations with lower-mass stars, reaching 40%–60% for 0.6 and 0.8 solar-mass stars, and close to 100% for $\lt 0.4$ solar-mass stars, in line with previous results (Raymond et al. 2007; Ciesla et al. 2015). The 1σ range of water fractions spans more than an order of magnitude, significantly larger than the stochastic dispersion in the N-body simulations (a factor of ∼2 to ∼5).

Using larger grains for calculating the snow line location significantly changes the results. Because the snow line is located closer to the star, water delivery is more efficient: the fraction of dry habitable-zone planets in our simulations decreases to zero for 0.8 and 1.0 solar-mass stars, and to 10% for $0.6\;{M}_{\odot }$. A small fraction of habitable-zone planets around low-mass stars does receive Earth-like amounts of water, around 20% for $0.2\;{M}_{\odot }$ and $0.4\;{M}_{\odot }$.

We derived the fraction of water-bearing terrestrial planets in the habitable zone by assuming that they are formed in situ. There is, however, growing observational evidence that migration of planets or their building blocks plays an important role in the formation of super-Earths and mini-Neptunes at short orbital periods (Hansen & Murray 2012; Swift et al. 2013; Raymond & Cossou 2014; Mulders et al. 2015). Whether this mechanism is able to form smaller, Earth-sized planets farther out remains an open question. If this were the case, the population of terrestrial planets in the habitable zone might be a mixture of planets formed in situ and through migration. As migrating planets typically form beyond the snow line (e.g., Cossou et al. 2014; Izidoro et al. 2014), they are expected to have water mass fractions of order 50% (Ogihara & Ida 2009; Tian & Ida 2015). Hence, the fractions of wet terrestrial planets quoted in this paper can be considered as lower limits.

An additional concern for the habitability of habitable-zone planets around lower-mass stars is that water loss may be more efficient as these planets form hotter than those around Sun-like stars (Lissauer 2007). Ramirez & Kaltenegger (2014) show that duing the bright pre-main-sequence phase of low-mass stars, stellar fluxes are high enough to trigger a runaway greenhouse effect that leads to enhanced water loss. Tian & Ida (2015) come to a similar conclusion by modeling water loss through hydrodynamc escape in a planet population synthesis model. For a planet to remain habitable, water has to arrive later during the pre-main-sequence evolution of the star. The timing and delivery mechanism of water may be crucial here: in our simulations, water delivery persists over ∼100 Myr timescales. In addition, water is not delivered in the form of ices that would be directly deposited in the atmosphere. Rather, it arrives as hydrated minerals that would enter the atmosphere later via vulcanic outgassing. A detailed study of arrival times and atmospheric release of water may be neccessary to address whether water can be deliverd late enough and in sufficient quantities to avoid escape.

Despite water delivery to low-mass planets in the habitable zone being less efficient than around Sun-like stars (Lissauer 2007; Raymond et al. 2007), we show that a small fraction may still receive Earth-like amounts of water due to the dispersion in snow-line locations. When taking into account the larger number of M stars with respect to G stars (factor of 10), their higher planet-occurrence rates (factor of two; Mulders et al. 2015), and increased water delivery from more comet-like icy bodies (Ciesla et al. 2015), the majority of water-rich terrestrial planets may still be found around low-mass M stars.

This work is partially supported by the National Aeronautics and Space Administration under Agreement No. NNX15AD94G, Earths in other Solar Systems, issued through the Science Mission Directorate interdivisional initiative Nexus for Exoplanet System Science.

Throughout this paper, we have assumed that the midplane temperature at the location of the snow line is determined by viscous heating (Equation (1)), and irradiation is negligable. For Sun-like stars, Min et al. (2011) have shown that this is a good approximation for the mass accretion rates considered in this paper (see also Garaud & Lin 2007; Oka et al. 2011), assuming that accretion is not layered (e.g., Lesniak & Desch 2011). Around lower-mass stars, lower mass accretion rates and a brighter pre-main-sequence phase may lead to a larger role of irradiation. To verify the assumption that Equation (1) also holds for low-mass stars, we compute the original radiative transfer model from Min et al. (2011) for lower-mass stars. For the stellar photosphere, we use the Baraffe et al. (1998) evolutionary tracks at t = 1 Myr. We explore a ±4σ range in mass accretion rates for the stellar masses used in this paper. Figure 5 shows the location of the snow line (i.e., the location where half the water is condensed into ice) divided by the predicted location from Equation (1).

Zoom In Zoom Out Reset image size

For stellar masses larger than $0.4\;{M}_{\odot }$, the location of the snow line is predicted accurately down to 2σ below the median within the precision of the radiative transfer model (∼5%). At smaller mass accretion rates, irradiation becomes imporant and Equation (1) starts underpredicting the location of the snow line. For the $0.2\;{M}_{\odot }$ case, irradiation has a larger influence, and the predicted location of the snow line starts deviating at 2σ below the median.