Open Access
Issue
A&A
Volume 562, February 2014
Article Number A80
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201322258
Published online 10 February 2014

© ESO, 2014

1. Introduction

Exoplanet transit photometry opened a new window to the interior and atmosphere of exoplanets. The biggest advantage of this technique would be that planetary radii are measured, while planetary masses are measured via other techniques, such as the radial velocity method and the transit timing variation method. Measured mass and radius relationships help us infer the internal structure and bulk composition of exoplanets theoretically, which give crucial constraints to formation and evolution processes of the planets. A growing number of small-sized exoplanets with radii of 1 to 2 R have been identified, which are often referred to as super-Earths (Batalha et al. 2013). Also, planet candidates detected by the Kepler space telescope include sub-Earth-sized objects, such as Kepler-20 e (Fressin et al. 2011), Kepler-42 b, c, d (Muirhead et al. 2012), and Kepler-37 b, c (Barclay et al. 2013). We can thus discuss the compositions of such small planets to gas giants by comparing theory with current observations.

Transiting super-Earths detected so far show a large variation in radius, suggesting diversity in composition. There are many theoretical studies on mass-radius relationships for planets with various compositions and masses (Valencia et al. 2007; Fortney et al. 2007; Sotin et al. 2007; Seager et al. 2007; Grasset et al. 2009; Wagner et al. 2011; Swift et al. 2012). A recent important finding, which compares theory to observation is that there are a significant number of low-density super-Earths that are larger in size than they would be if they were rocky. This implies that these transiting super-Earths possess components less dense than rock. From a viewpoint of planet formation, the possible components are hydrogen-rich gas and water, which make an outer envelope. A small fraction of H-rich gas or water is known to be enough to account for observed radii of the low-density super-Earths (Adams et al. 2008; Valencia et al. 2010).

The stability of the envelopes are, however, to be examined. Transiting planets are generally orbit close to their host stars (typically ≲0.1 AU), because detection probability of planetary transits is inversely proportional to the semi-major axis (e.g., Kane 2007). These close-in planets are highly irradiated and exposed to intense X-ray and ultraviolet radiation (hereafter XUV) that come from their host stars. This causes the planetary envelope to escape hydrodynamically from the planet (e.g., Watson et al. 1981). This process is often called the photo-evaporation of planetary envelopes. As for massive close-in planets, namely, hot Jupiters, the possibility of the photo-evaporation and its outcome have been investigated well both theoretically and observationally (e.g., Yelle et al. 2008 and references therein).

While the photo-evaporation may not significantly affect the evolution and final composition of hot Jupiters except for extremely irradiated or inflated hot Jupiters, its impact on small close-in planets in the sub/super-Earth mass range should be large, partly because their envelope masses are much smaller than those of hot Jupiters. For example, Valencia et al. (2010) investigated the structure and composition of the first transiting super-Earth CoRoT-7 b and discussed the sustainability of the possible H + He envelope with a mass of less than 0.01% of the total planetary mass. The envelope mass was consistent with its measured mass and radius. The estimated lifetime of the H + He envelope was, however, only 1 million years, which was much shorter than the host star’s age (2–3 Gyr). This suggests that CoRoT-7 b is unlikely to retain the H + He envelope at present.

Young main-sequence stars are known to be much more active and emit stronger XUV than the current Sun (e.g., Ribas et al. 2005). Therefore, even if a super-Earth had a primordial atmosphere initially, it may lose the atmosphere completely during its history. These discussions concerning the photo-evaporative loss of H + He envelopes were done for GJ 1214 b (Nettelmann et al. 2011; Valencia et al. 2013), super-Earths orbiting Kepler-11 (Lopez et al. 2012; Ikoma & Hori 2012), and CoRoT-7 b (Valencia et al. 2010). Systematic studies were also done by Rogers et al. (2011) and Lopez & Fortney (2013). Those studies demonstrated the large impact of the photo-evaporation on the stability of H + He envelopes for super-Earths. In particular, Lopez & Fortney (2013) performed simulations of coupled thermal contraction and photo-evaporative mass loss of rocky super-Earths with H + He envelopes. They found that there were threshold values of planetary masses and radii, below which H + He envelopes were completely stripped off. Owen & Wu (2013) also performed similar simulations with detailed consideration of the mass loss efficiency for an H + He envelope based on Owen & Jackson (2012). They argued that evaporation explained the correlation between the semi-major axes and planetary radii (or planet densities) of KOIs.

In this study, we focus on water-rich sub/super-Earths. Planet formation theories predict that low-mass planets migrate toward their host star, which is strongly supported by the presence of many close-in super-Earths, from cooler regions (e.g., Ward 1986) where they may have accreted a significant amount of water. This suggests that water/ice-rich sub/super-Earths may also exist close to host stars. Therefore, similar discussions should be done for water envelopes of close-in super-Earths. However, there are just a few studies, which treat specific sub/super-Earths such as CoRoT-7 b (Valencia et al. 2010) and Kepler-11 b (Lopez et al. 2012). No systematic study is yet to be done for the stability of water envelopes.

The purpose of this study is, thus, to examine the stability of primordial water envelopes of close-in sub/super-Earths against photo-evaporation. To this end, we simulate the thermal evolution of planets with significant fractions of water envelopes (i.e., water-worlds), incorporating the effect of stellar-XUV-driven photo-evaporative mass loss. The theoretical model is described in Sect. 2. As for the atmosphere model, the details are described in Appendix A. In Sect. 3, we show the evolutionary behavior of the water-rich planets. Then, we find threshold values of planetary masses and radii below which such water-rich planets are incapable of retaining primordial water envelopes for a period similar to ages of known exoplanet-host stars (i.e., 1–10 Gyr). In Sect. 4, we compare the theoretical mass-radius distribution of water-rich planets with that of known transiting planets. Furthermore, we compare the threshold radius with sizes of Kepler objects of interest (KOIs) to suggest that KOIs include a significant number of rocky planets. Finally, we summarize this study in Sect. 5.

thumbnail Fig. 1

Model of the planetary structure in this study.

2. Numerical models

In this study, we simulate the evolution of the mass and radius of a planet that consists of water and rock, including the effects of mass loss due to photo-evaporation. The structure model is depicted in Fig. 1. The planet is assumed to consist of three layers in spherical symmetry and hydrostatic equilibrium: namely, from top to bottom, it consisted of a water vapor atmosphere, a water envelope, and a rocky core. At each interface, the pressure and temperature are continuous.

The assumptions and equations that determine the planet’s interior structure and thermal evolution are described in Sects. 2.1 and 2.2, respectively. The equations of state for the materials in the three layers are summarized in Sect. 2.3. The structure of the atmosphere and the photoevaporative mass loss, both of which govern the planet’s overall evolution, are described in Sect. 2.4 (see also Appendix A) and Sect. 2.5, respectively. Since a goal of this study is to compare our theoretical prediction with results from transit observations, we also calculate the transit radius, which is defined in Sect. 2.6. Finally, we summarize our numerical procedure in Sect. 2.7.

2.1. Interior structure

The interior structure of the planet is determined by the differential equations (e.g. Kippenhahn & Weigert 1990), ∂PMr=GMr4πr4,∂rMr=14πr2ρ,∂TMr=GMrT4πr4P,\begin{eqnarray} \label{IE01}\frac{\partial P}{\partial M_r} &=& -\frac{GM_r}{4\pi r^4}, \\ \label{IE02}\frac{\partial r}{\partial M_r} &=& \frac{1}{4\pi r^2 \rho},\\ \label{IE03}\frac{\partial T}{\partial M_r} &=& -\frac{GM_r T}{4\pi r^4 P}\nabla, \end{eqnarray}and the equation of state, ρ=ρ(P,T),\begin{eqnarray} \label{IE04}\rho = \rho (P, T), \end{eqnarray}(4)where r is the planetocentric distance, Mr is the mass contained in the sphere with radius of r, P is the pressure, ρ is the density, T is the temperature, and G (=6.67 × 10-8 dyn cm2 g-2) is the gravitational constant. The symbol ∇ is the temperature gradient with respect to pressure. We assume that the water envelope and rocky core are fully convective and the convection is vigorous enough that the entropy S is constant; namely, =(lnTlnP)S·\begin{equation} \nabla = \left( \frac{\partial\ln T}{\partial\ln P} \right)_{S}\cdot \label{IE+02}\vspace*{-1mm} \end{equation}(5)Equations (1)–(3) require three boundary conditions. The inner one is r = 0 at Mr = 0. The outer boundary corresponds to the interface between the envelope and the atmosphere, which is called the tropopause. The tropopause pressure Pad and temperature Tad are determined from the atmospheric model; the details of which is described in Sect. 2.4 and Appendix A. The atmospheric mass is negligible, relative to the planet total mass Mp. In our calculation, the atmospheric mass is less than 0.1% of the planetary mass. Thus, the outer boundary conditions are given as P=PadandT=TadatMr=Mp.\begin{equation} \vspace*{-2mm} \begin{array}{lllll} P = P_{\rm{ad}} & \mbox{and} & T = T_{\rm{ad}} & \mbox{at} & M_r = M_{\rm p}. \label{Bd01} \vspace*{-2mm} \end{array} \end{equation}(6)As mentioned above, the pressure and temperature are also continuous at the interface between the water envelope and the rocky core.

2.2. Thermal evolution

The thermal evolution of the planet without internal energy generation is described by (e.g., Kippenhahn & Weigert 1990) ∂LMr=T∂S∂t,\begin{equation} \frac{\partial L}{\partial M_r} = -T \frac{\partial S}{\partial t}, \label{Ec04} \end{equation}(7)where L is the intrinsic energy flux passing through the spherical surface with radius of r, S is the specific entropy, and t is time. Since the entropy is constant in each layer, the integrated form of Eq. (7) is written as Lp=e∂tMcMpTdMr+c∂t0McTdMr,\begin{equation} -L_{\rm p} = \frac{\partial \bar{S}_{\rm e}}{\partial t} \int_{M_{\rm c}}^{M_{\rm p}} T {\rm d}M_r +\frac{\partial \bar{S}_{\rm c}}{\partial t} \int_{0}^{M_{\rm c}} T {\rm d}M_r , \label{Ec05} \vspace*{-2mm} \end{equation}(8)where Lp is the total intrinsic luminosity of the planet, Mc is the mass of the rocky core, and \hbox{$\bar{S}_{\rm e}$} and \hbox{$\bar{S}_{\rm c}$} are the specific entropies in the water envelope and the rocky core, respectively. In integrating Eq. (7), we have assumed L = 0 at Mr = 0.

In the numerical calculations of this study, we use the intrinsic temperature Tint, instead of Lp, which is defined by Tint4Lp4πRp2σ,\begin{equation} T_{\rm{int}}^{4} \equiv \frac{L_{\rm p}}{4\pi R_{\rm p}^{2}\sigma}, \label{Ec06}\vspace*{-1mm} \vspace*{-2mm} \end{equation}(9)where Rp is the planet photospheric radius (see Sect. 2.4 for the definition) and σ is the Stefan-Boltzmann constant (=5.67 × 10-5  erg  cm-2  K-4  s-1).

2.3. Equation of state (EOS)

In the vapor atmosphere, the temperature and pressure are sufficiently high and low, respectively, so that the ideal gas approximation is valid. We thus adopt the ideal equation of state, incorporating the effects of dissociation of H2O. In practice, we use the numerical code developed by Hori & Ikoma (2011), which calculate chemical equilibrium compositions among H2O, H2, O2, H, O, H+, O+ and e.

At high pressures in the water envelope, the ideal gas approximation is no longer valid, because pressure due to molecular interaction is not negligible. In this study, we mainly use the water EOS H2O-REOS (Nettelmann et al. 2008), which contains the ab initio water EOS data at high pressures of French et al. (2009). H2O-REOS covers a density range from 1.0 × 10-6  g  cm-3 to 15  g  cm-3 and a temperature range from 1.0 × 103 K to 2.4 × 104 K. For T and ρ outside the ranges that H2O-REOS covers, we use SESAME 7150 (Lyon & Johnson 1992).

The rocky core is assumed to be mineralogically the same in composition as the silicate Earth. We adopt a widely-used EOS and the Vinet EOS, and calculate thermodynamic quantities following Valencia et al. (2007).

2.4. Atmospheric model

As described above, we consider an irradiated, radiative-equilibrium atmosphere on top of the water envelope. The thermal properties of the atmosphere govern the internal structure and evolution of the planet. To integrate the atmospheric structure, we follow the prescription developed by Guillot (2010) except for the treatment of the opacity. Namely, we consider a semi-grey, plane-parallel atmosphere in local thermal equilibrium. The wavelength domains of the incoming (stellar) and outgoing (planetary) radiations are assumed to be completely separated; the former is visible, while the latter is near or mid infrared.

We solve the equation of radiative transfer by integrating the two sets (for incoming and outgoing radiations) of the zeroth and first-order moment equations for radiation with the Eddington’s closure relation: the incoming and outgoing radiations are linked through the equation of radiative equilibrium (see Eqs. (10), (11) and (17)–(19) of Guillot 2010). Guillot (2010) derived an analytical, approximate solution, which reproduced the atmospheric structure from detailed numerical simulations of hot Jupiters (see also Hansen 2008) well. The solution depends on opacities in the visible and thermal domains. Guillot (2010) also presented empirical formulae for the mean opacities of solar-composition (i.e., hydrogen-dominated) gas.

However, no empirical formula is available for opacities of water vapor of interest in this study. We take into account the dependence of the water-vapor opacity on temperature and pressure and integrate the momentum equations numerically. The details about the mean opacities and momentum equations are described in Appendix A.

The bottom of the atmosphere is assumed to be the interface between the radiative and convective zones. We use the Schwarzschild criterion (e.g., see Kippenhahn & Weigert 1990) to determine the interface. The pressure and temperature at the interface (Pad, Tad) are used as the outer boundary conditions for the structure of the convective water envelope.

The photospheric radius Rp used in Eq. (9) is the radius at which the thermal optical depth measured from infinity, τ, is 2/3; namely, τ=Rpκthrρdr=23,\begin{equation} \tau=\int_{R_{\rm p}}^{\infty}\kappa_{\mathrm{th}}^{\mathrm{r}}\rho {\rm d}r = \frac{2}{3}, \label{AT05} \end{equation}(10)where κthr\hbox{$\kappa_\mathrm{th}^{r}$} is the Rosseland mean opacity for the outgoing radiation (see Appendix A for the definition). This level is above the tropopause, the radius of which is written by Rconv (see Fig. 1). We evaluate the atmospheric thickness z (=Rp − Rconv) by integrating the hydrostatic equation from P = Pad to P = Pph using z=PadPphdP=PadPphμgTPdP,\begin{equation} z =- \int_{P_{\mathrm{ad}}}^{P_{\rm{ph}}} \frac{ {\rm d}P}{g\rho} = -\int_{P_{\mathrm{ad}}}^{P_{\rm{ph}}} \frac{\mathcal R}{\mu g} \frac{T}{P} {\rm d}P, \label{AT03} \end{equation}(11)where g is the constant gravity, ℛ (=8.31 × 107 erg K-1 g-1) is the gas constant, and μ is the mean molecular weight. Pph is the photospheric pressure that we calculate by integrating dPdτ=gκthr\begin{equation} \frac{ {\rm d}P}{ {\rm d}\tau} = \frac{g}{\kappa_{\mathrm{th}}^{\rm{r}}} \label{AT06a} \end{equation}(12)from τ = 0 to 2/3.

2.5. Mass loss

The mass loss is assumed to occur in an energy-limited fashion. Its rate, including the effect of the Roche lobe, is given by (Erkaev et al. 2007) =εFXUVRpπRXUV2GMpKtide,\begin{equation} \dot{M} =-\frac{\varepsilon F_{\mathrm{XUV}}R_{\rm p}\pi R_{\mathrm{XUV}}^2 }{GM_{\rm p} K_{\mathrm{tide}}}, \label{ML01} \end{equation}(13)where ε is the heating efficiency, which is defined as the ratio of the rate of heating that results in hydrodynamic escape to that of stellar energy absorption; FXUV is the incident flux of X-ray and UV radiation from the host star, Ktide is the potential energy reduction factor due to stellar tide; and RXUV is the effective radius at which the planet receives the incident XUV flux. In Eq. (13), we have assumed RXUV = Rp, which is a good approximation for close-in planets of interest (Lammer et al. 2013). It is noted that Lammer et al. (2013) focused on the hydrogen-helium atmosphere. Since the scale height of the vapor atmosphere is smaller than that of a hydrogen-helium atmosphere with the same temperature, RXUV ≃ Rp is a good approximation also for the vapor atmosphere.

In this study, we suppose that the host star is a G-star and adopt the empirical formula derived by Ribas et al. (2005) for FXUV: FXUV={\begin{eqnarray} F_{\rm{XUV}} = \left\{ \begin{array}{lc} 504 \left( \displaystyle \frac{a}{1\mathrm{~AU}} \right)^{-2} \mathrm{erg}~\mathrm{s}^{-1}~\mathrm{cm}^{-2} &(t<0.1\rm{~Gyr}) \\ \label{ML04} 29.7\left(\displaystyle \frac{t}{1\mathrm{~Gyr}} \right)^{-1.23} \left( \displaystyle \frac{a}{1\mathrm{~AU}} \right)^{-2} \mathrm{erg}~\mathrm{s}^{-1}~\mathrm{cm}^{-2} & (t\ge0.1\rm{~Gyr}). \label{ML02} \end{array} \right. \end{eqnarray}(14)\normalsizeWe use the formula for Ktide derived by Erkaev et al. (2007), Ktide=(η1)2(2η+1)2η3,\begin{equation} K_{\mathrm{tide}} = \frac{(\eta-1)^2(2\eta+1)}{2\eta^3}, \label{ML05} \end{equation}(15)where η is the ratio of the Roche-lobe (or Hill) radius to the planetary radius, Rp.

The value of the heating efficiency is uncertain, because minor gases such as CO2 contribute to it via radiative cooling. For photo-evaporation of hot-Jupiters, ε is estimated to be on the order of 0.1 (Yelle et al. (2008) and reference therein). Thus, we adopt ε = 0.1 as a fiducial value and investigate the sensitivity of our results to ε.

Finally, we assume that the rocky core never evaporates. That is simply because we are interested in the stability of water envelopes in this study. Whether rocky cores evaporate or not is beyond the scope of this study.

2.6. Transit radius

The planetary radius measured via transit photometry is different from the photospheric radius defined in the preceding subsection. The former is the radius of the disk that blocks the stellar light ray that grazes the planetary atmosphere in the line of sight. This radius is called the transit radius hereafter in this paper. Below we derive the transit radius, basically following Guillot (2010). Note that Guillot (2010) assumed the plane-parallel atmosphere, while we consider a spherically symmetric structure, because the atmospheric thickness is not negligibly small relative to the planetary radius in some cases in this study.

thumbnail Fig. 2

Concept of the chord optical depth.

We first introduce an optical depth that is called the chord optical depth, τch (e.g. Guillot 2010). The chord optical depth is defined as τch(r,ν)=+ρκνds,\begin{equation} \tau_{\mathrm{ch}}(r,\nu) = \int_{-\infty}^{+\infty} \rho\kappa_{\nu} {\rm d}s, \label{TR01} \end{equation}(16)where r is the planetocentric distance of the ray of interest (see Fig. 2), s is the distance along the line of sight measured from the point where the line is tangent to the sphere, and κν is the monochromatic opacity at the frequency ν. Using τch, we define the transit radius, Rtr, as τch(Rtr)=23·\begin{equation} \tau_{\mathrm{ch}}(R_{\mathrm{tr}}) = \frac{2}{3}\cdot \label{TR+00} \end{equation}(17)Let the altitude from the sphere of radius r be ztr. Then s2 = (r + ztr)2 − r2 (Fig. 2). Eq. (16) is written as τch(r,ν)=20ρκνztr+rztr2+2rztrdztr.\begin{equation} \tau_{\mathrm{ch}}(r,\nu) = 2\int_{0}^{\infty} \rho\kappa_{\nu}\frac{z_{\mathrm{tr}}+r}{\sqrt{z_{\mathrm{tr}}^2+2rz_{\mathrm{tr}}}} {\rm d}z_{\mathrm{tr}}. \label{TR02} \end{equation}(18)Furthermore for convenience, we choose pressure P as the independent variable, instead of ztr. Using the equation of hydrostatic equilibrium, dPdztr=GMpρ(r+ztr)2,\begin{equation} \frac{ {\rm d}P}{ {\rm d}z_{\mathrm{tr}}} = - \frac{GM_{\rm p}\rho}{(r+z_{\mathrm{tr}})^2}, \label{TR-01} \end{equation}(19)one obtains τch(ν,r)=2grPr0κν(1+ztr/r)3(1+ztr/r)21dP,\begin{equation} \tau_{\mathrm{ch}}(\nu,r) =-\frac{2}{g_r} \int_{P_r}^{0} \kappa_{\nu} \frac{(1+z_{\rm{tr}}/r)^3}{\sqrt{(1+z_{\rm{tr}}/r)^{2}-1}} {\rm d}P, \label{TR03} \end{equation}(20)where gr=GMr2\begin{equation} g_r = \frac{GM}{r^2} \label{TR-04b} \end{equation}(21)and Pr is the pressure at r. To integrate Eq. (20), we write ztr as a function of P. To do so, we integrate Eq. (19) and obtain 0ztrdz(r+z)2=PrPzdPGMpρ,\begin{equation} \int_{0}^{z_{\mathrm{tr}}}\frac{ {\rm d}z'}{(r+z')^2} = - \int_{P_r}^{P_z} \frac{ {\rm d}P}{GM_{\rm p}\rho}, \label{TR-02} \end{equation}(22)where Pz is the pressure at ztr. Eq. (22) is integrated as 1r+ztr=1r1r2grPzPrdPρ=1rzp(Pr,Pz)r2,\begin{eqnarray} \frac{1}{r+z_{\mathrm{tr}}} &=& \frac{1}{r} - \frac{1}{r^2g_r} \int_{P_z}^{P_r} \frac{ {\rm d}P}{\rho} \nonumber \\ \label{TR-03}&=& \frac{1}{r} - \frac{z_{\rm p}(P_r, P_z)}{r^2}, \end{eqnarray}(23)where zp(Pr,Pz)PzPrPρgrdlnP.\begin{equation} z_{\rm p}(P_r, P_z) \equiv \int_{P_z}^{P_r} \frac{P}{\rho g_r} {\rm d}\ln P. \label{TR-04a} \end{equation}(24)Thus, z is written as ztr=zp(1zpr)-1.\begin{equation} z_{\mathrm{tr}} = z_{\rm p} \left( 1-\frac{z_{\rm p}}{r} \right)^{-1}. \label{TR-05} \end{equation}(25)Note that zp corresponds to the altitude in the case of a plane-parallel atmosphere and (1 − zp/r)-1 is the correction for spherical symmetry.

2.7. Numerical procedure

To simulate the mass and radius evolution simultaneously, we integrate Eqs. (8) and (13) by the following procedure.

First, we simulate two adiabatic interior models that are separated in time by a time interval Δt for the known Mp(t) and an assumed Mp(t + Δt). To be exact, the two structures are integrated for two different values of Tint. In doing so, we integrate Eqs. (1)–(4) inward from the tropopause to the planetary center, using the fourth-order Runge-Kutta method. The inward integration is started with the outer boundary condition given by Eq. (6); Pad and Tad are calculated according to the atmospheric model described in Sect. 2.4. We then look for the solution that fulfills the inner boundary condition (i.e., r = 0 at Mr = 0) in an iterative fashion. Note that determining Pad and Tad requires the gravity in the atmosphere (or Rconv), which is obtained after the interior structure is determined. Thus, we have to find the solution in which the interior and atmospheric structures are consistent with each other also in an iterative fashion.

Then we calculate Δt from the second-order difference equation for Eq. (8), which is written as Δt=([e(t+Δt)e(t)][Θe(t+Δt)+Θe(t)]+[c(t+Δt)c(t)][Θc(t+Δt)+Θc(t)])×(Lp(t+Δt)+Lp(t))-1,\begin{eqnarray} \label{NM01} \Delta t &=& - \left( \left[ \bar{S}_{\rm e} (t+\Delta t) - \bar{S}_{\rm e} (t) \right] \left[ \Theta_{\rm e} (t+\Delta t) + \Theta_{\rm e} (t) \right] \right.\nonumber\\ &&\left. +\left[ \bar{S}_{\rm c} (t\!+\!\Delta t)\! - \!\bar{S}_{\rm c} (t) \right] \left[ \Theta_{\rm c} (t\!+\!\Delta t) \!+\! \Theta_{\rm c} (t) \right] \right)\times \left(L_{\rm p} (t\!+\!\Delta t) \!+\! L_{\rm p} (t)\right)^{-1}\!, \end{eqnarray}(26)where Θe(t)McMp(t)T(t)dMr,Θc(t)0McT(t)dMr.\begin{equation} \Theta_{\rm e} (t) \equiv \int^{M_{\rm p} (t)}_{M_{\rm c}} T (t) {\rm d}M_r, \hspace{3ex} \Theta_{\rm c} (t) \equiv \int^{M_{\rm c}}_{0} T (t) {\rm d}M_r. \label{NM02} \end{equation}(27)Using this Δt, we integrate Eq. (13) to calculate Mp(t + Δt) as Mp(t+Δt)=Mp(t)+Δt.\begin{equation} M_{\rm p} (t + \Delta t) = M_{\rm p} (t) + \dot{M} \Delta t. \label{NM03} \end{equation}(28)The assumed Mp(t + Δt) is not always equal to that obtained here. Therefore the entire procedure must be repeated until the Mp(t + Δt) in Eq. (28) coincides with that assumed for calculating Eq. (26) with satisfactory accuracy, which is ≲0.1% in our simulations.

Once we obtain the interior and atmospheric structure, we calculate the transit radius by the procedure described in Sect. 2.6. Finally, we have confirmed that our numerical code reproduces the mass and radius relationship for super-Earths well which is presented by Valencia et al. (2010).

3. Mass evolution

In this section, we show our numerical results of the mass evolution of a close-in water-rich planet. The evolution is controlled by the following five parameters: the initial total mass of the planet (Mp,0), the initial luminosity (L0), the initial water mass fraction (Xwt,0), the semi-major axis (a), and the heating efficiency (ε). Below, we adopt L0 = 1 × 1024 erg s-1, Xwt,0 = 75%, a = 0.1 AU, and ε = 0.1 as fiducial values unless otherwise noted. We also show how the five parameters affect the fate of a close-in water-rich planet.

thumbnail Fig. 3

Mass evolution of close-in water-rich planets. The solid blue lines represent planets that retain their water envelopes for 10 Gyr. In contrast, the planet shown by the dashed red line loses its water envelope completely in 10 Gyr. We set Lp,0 = 1 × 1024 erg s-1, Xwt,0 = 75%, a = 0.1 AU, and ε = 0.1 for all the planets. In this model, we assume that the rocky core never evaporates.

3.1. Examples of mass evolution

Figure 3 shows examples of the mass evolution for water-rich planets with six different initial masses; L0 = 1 × 1024 erg s-1, Xwt,0 = 75%, a = 0.1 AU, and ε = 0.1 in these simulations, as stated above. The smallest planet loses its water envelope completely in 1 Gyr (the dashed line), while more massive planets retain their water envelopes for 10 Gyr (solid lines). This means that a water-rich planet below a threshold mass ends up as a naked rocky planet.

The presence of such a threshold mass is understood in the following way. Using Eq. (13), we define a characteristic timescale of the mass loss (τM) as τM=|XwtMpp|=4GKtideXwtMpρpl3εFXUV,\begin{equation} \tau_{M} = \left| \frac{X_{\rm{wt}}M_{\rm p}}{\dot{M}_{\rm p}} \right| = \frac{4GK_{\rm{tide}}X_{\rm{wt}}M_{\rm p}\rho_{\rm{pl}}}{3\varepsilon F_{\rm{XUV}}}, \label{result_tm_01} \end{equation}(29)where ρpl is the mean density of the planet. As the planetary mass decreases, the mass-loss timescale becomes shorter. This trend is enhanced by the M − ρ relationship that the mean density decreases as Mp decreases, according to our numerical results for water-rich planets.

In addition, the time-dependence of the stellar XUV flux (see Eq. (14)) is a crucial factor to cause a striking difference in behavior between the low-mass and high-mass planets. Using Eq. (14), we obtain the following relation for τM: τM{\begin{eqnarray} \tau_{M} \simeq \left\{ \begin{array}{lcc} 3\times 10^8 f~~\rm{yr}, &\mathrm{for}&\mathit{t} < 0.1~\mathrm{Gyr}, \label{result_tm_02} \\ 3\times 10^8 \displaystyle \left( \frac{t}{0.1~\rm{Gyr}}\right)^{1.23} f~~\rm{yr}, &\mathrm{for}&\mathit{t} \ge 0.1~\mathrm{Gyr}, \label{result_tm_03} \end{array} \right. \end{eqnarray}(30)where f=1(a0.1AU)2(XwtMpM)(ρpl0.1gcm-3)(Ktide0.9)(ε0.1)-1.\begin{equation} f = 1 \left(\frac{a}{0.1\rm{~AU}} \right)^{2} \left(\frac{X_{\rm{wt}}M_{\rm p}}{M_{\oplus}} \right) \left(\frac{\rho_{\rm{pl}}}{0.1~ \rm{g}~\rm{cm}^{-3}} \right) \left( \frac{K_{\rm{tide}}}{0.9} \right) \left(\frac{\varepsilon}{0.1}\right)^{-1}. \label{result_tm_03a} \end{equation}(31)Note that 0.1 g cm-3 is a typical value of ρpl in the case of sub-Earth-mass planets with the age of 108 years, according to our calculations. As seen in Eq. (30), τM becomes longer rapidly with time. This implies that small planets that satisfy τM < 0.1 Gyr experience a significant mass loss. In other words, massive planets that avoid significant mass loss in the early phase hardly lose their mass for 10 Gyr. Thus, there exists a threshold mass below which a planet never retains its water envelope for a long period. Our numerical calculations found that the threshold mass (hereafter Mthrs) is 0.56  M for the fiducial parameter set, which is in good agreement with Mp < 0.4  M as derived from Eq. (30).

A similar threshold mass was found by Lopez & Fortney (2013) for H + He atmospheres of rocky planets. Hydrogen-rich planets are more vulnerable to the photo-evaporative mass loss than water-rich planets. According to their study, the threshold mass of the hydrogen-rich planet at 0.1 AU is ~5  M. That is, Mthrs for water-rich planets is smaller by a factor of ~10 than that of hydrogen-rich planets.

3.2. Dependence on the initial planet’s luminosity

The evolution during the first 0.1 Gyr determines the fate of a water-rich planet, as shown above. Such a trend is also shown by Lopez & Fortney (2013) for H + He atmospheres of rocky planets. This suggests that the sensitivity of the planet’s fate to the initial conditions must be checked. In particular, the initial intrinsic luminosity may affect the early evolution of the planet significantly, because the planetary radius, which has a great impact on the mass loss rate, is sensitive to the intrinsic luminosity; qualitatively, a large L0 enhances mass loss because of a large planetary radius. On the other hand, L0 is uncertain, because it depends on how the planet forms (e.g. accretion processes of planetesimals, migration processes and giant impacts). However, as shown below, the fate of the planet is insensitive to choice of L0, Fig. 4 shows Mthrs as a function of L0 for a = 0.02,0.03,0.05 and 0.1 AU. We have found that Mthrs is almost independent of L0. This is because an initially-luminous planet cools down rapidly, so that the integrated amount of water loss during the high-luminosity phase is negligible. This is confirmed by the following argument. The mass loss, ΔM, at the early stage can be estimated by ΔM~τKH,\begin{equation} \Delta M \sim \dot{M} \tau_{\rm{KH}}, \label{result_lm_01} \end{equation}(32)where τKH is the typical timescale of Kelvin-Helmholtz contraction, τKHGMp22RpLp·\begin{equation} \tau_{\rm{KH}} \simeq \frac{GM_{\rm p}^2}{2R_{\rm p}L_{\rm p}}\cdot \label{evolv_01} \end{equation}(33)With Eqs. (29) and (33) given, Eq. (32) can be written as ΔM~MpτKHτM=Mpε2Ktide·πRp2FXUVLp~3×10-2(FXUV504ergcm-2s-1)(ε0.1)(Ktide0.9)-1×(a0.1AU)-2(Rp3R)2(Lp1024ergs-1)-1Mp.\begin{eqnarray} \label{result_lm_02} \Delta M &\sim& M_{\rm p} \frac{\tau_\mathrm{KH}}{\tau_M} = M_{\rm p}\frac{\varepsilon}{2K_{\rm{tide}}} \cdot \frac{\pi R_{\rm p}^2 F_{\rm{XUV}}}{L_{\rm p}} \\ &\sim& 3\times 10^{-2} \left( \frac{F_{\rm{XUV}}}{504~\rm{erg}~\rm{cm}^{-2}~\rm{s}^{-1}} \right) \left(\frac{\varepsilon}{0.1} \right) \left( \frac{K_\mathrm{tide}}{0.9}\right)^{-1} \nonumber \\ \label{result_lm_03}&&\times \left( \frac{a}{0.1~\rm{AU}} \right)^{-2} \left( \frac{R_{\rm p}}{3~R_{\oplus}} \right)^2 \left( \frac{L_{\rm p}}{10^{24}~\rm{erg}~\rm{s}^{-1}} \right)^{-1} M_{\rm p}. \end{eqnarray}Because FXUV is constant in the early phase, ΔM decreases as Lp increases; that is, the Kelvin-Helmholtz contraction proceeds more rapidly. Therefore, the choice of the value of L0 has little effect on the total amount of water loss, as far as L0 is larger than 1024  erg  s-1. For smaller L0, Rp is insensitive to L0. Thus, Mthrs is insensitive to L0.

thumbnail Fig. 4

Threshold mass in M as a function of the initial planet’s luminosity in erg s-1 for four choices of semi-major axes. The solid (red), dashed (green), dotted (blue), and dot-dashed (purple) lines represent a = 0.02,0.03,0.05, and 0.1 AU, respectively. We have assumed Xwt = 75% and ε = 0.1.

3.3. Dependence on the initial water mass fraction

The fate of a water-rich planet also depends on the initial water mass fraction, Xwt,0. Figure 5 shows Xwt(t) at t = 10 Gyr as a function of the initial planet’s mass, Mp,0, for four different values of Xwt,0(=25%, 50%, 75%, and 100%). As Mp,0 decreases, Xwt(10 Gyr) decreases. The pure water planet (solid line) with Mp,0 < 0.82  M is completely evaporated in 10 Gyr; namely, Xwt(10  Gyr) = 0%. Otherwise, Xwt(10  Gyr) = 100%. In other cases, we find that the threshold mass, Mthrs, below which Xwt(10  Gyr) = 0%, is 0.56  M for Xwt,0 = 75%, 0.44  M for Xwt,0 = 50%, and 0.44  M for Xwt,0 = 25%.

thumbnail Fig. 5

Relationship between the initial planetary mass and the fraction of the water envelope at 10 Gyr for four initial water mass fractions of Xwt,0 = 100% (solid, red), 75% (dashed, green), 50% (dotted, blue), and 25% (dot-dashed, purple). We have assumed L0 = 1 × 1024  erg  s-1, a = 0.1 AU, and ε = 0.1.

Figure 6 shows the relationship between Xwt,0 and Mthrs for four different semi-major axes. Mthrs is found not to be a monotonous function of Xwt,0. For Xwt,0 < 25%, Mthrs decreases, as Xwt,0 increases. This is explained as follows. According to Eq. (29), the mass loss timescale , τM, depends on the absolute amount of water, XwtMp, and the planetary bulk density, ρpl. When Xwt is sufficiently small, ρpl is equal to the rocky density and is therefore constant. Thus, τM is determined only by the absolute amount of water (i.e., XwtMp). This means that, Mp must be larger for τM to be the same if Xwt,0 is small. As a consequence, Mthrs decreases with increasing Xwt,0. More exactly, Mthrs changes with Xwt,0 in such a way that Xwt,0Mthrs is constant. In contrast, when Xwt,0 is large, Xwt, Mp, and ρpl affect the mass loss timescale. For a given Mp, an increase in Xwt,0 leads to a decrease in ρpl (or, an increase in radius), which enhances mass loss. As a result, Mthrs increases with Xwt,0 for Xwt,0 > 25%. Therefore, there is a minimum value of Mthrs, which is hereafter denoted by Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}.

Similar trends can be seen in Figs. 3 and 4 of Lopez & Fortney (2013). To compare our results for water-rich planets to those for hydrogen-rich rocky planets from Lopez & Fortney (2013) in a more straightforward way, we show the relationship between the initial total mass and the fraction of the initial water envelope that is lost via subsequent photo-evaporation in 5 Gyr in Fig. 7 (see Fig. 3c of Lopez & Fortney 2013). We set L0 = 1 × 1024  erg  s-1, a = 0.1 AU, ε = 0.1, and six initial water mass fractions of Xwt,0 = 1% (solid, red), 3% (long-dashed, green), 10% (dotted, blue), 30% (dash-dotted, purple), 50% (dot-dashed, light blue), and 60% (dashed, black), which are similar to those adopted by Lopez & Fortney (2013). As mentioned above, the initial total mass needed in the H + He case is larger by a factor of ~10 than that in the water case for the same fraction of the initial envelope to survive photo-evaporation. In addition, the required initial total mass for Xwt,0 < 10% becomes significantly large in the water case. This behavior is also found in the case of the hydrogen-rich planets for Xwt,0 = 1–3%. However, the trend is less noticeable in the H + He case. This is because the density effect described above is effective even for small H + He fractions.

thumbnail Fig. 6

Relationship between the initial water mass fraction Xwt,0 in % and the threshold mass Mthrs in M for four choices of semi-major axes of 0.02 AU (solid, red), 0.03 AU (dashed, green), 0.05 AU (dotted, blue), and 0.1 AU (dot-dashed, purple). We have assumed L0 = 1 × 1024  erg  s-1 and ε = 0.1.

thumbnail Fig. 7

Relationship between the initial planetary mass and the fraction of the initial water envelope that is lost via photo-evaporation in 5 Gyr for six initial water mass fractions of Xwt,0 = 1% (solid, red), 3% (long-dashed, green), 10% (dotted, blue), 30% (dash-dotted, purple), 50% (dot-dashed, light blue), and 60% (dashed, black). We have assumed L0 = 1 × 1024  erg  s-1, a = 0.1 AU, and ε = 0.1.

3.4. Dependence on the semi-major axis

At small a, the incident stellar XUV flux becomes large. Thus, Mthrs increases, as a decreases. Certainly, the distance to the host star affects the equilibrium temperature Teq, which has an influence on ρpl: the higher Teq is, the smaller ρpl is. However, its impact on Mthrs is small, relative to that of FXUV. According to the planet’s mass and mean density relationship, ρpl differs only by a factor of ≲1.5 between 880 K and 2000 K. Therefore, increasing FXUV has a much greater impact on the mass loss than decreasing ρpl. In Fig. 6, we find Mthrs=5.2M\hbox{$M_{\rm{thrs}}^{\ast}=5.2~M_{\oplus}$} for a = 0.02 AU, Mthrs=2.5M\hbox{$M_{\rm{thrs}}^{\ast}=2.5~M_{\oplus}$} for a = 0.03 AU, Mthrs=1.2M\hbox{$M_{\rm{thrs}}^{\ast}=1.2~M_{\oplus}$} for a = 0.05 AU, and Mthrs=0.44M\hbox{$M_{\rm{thrs}}^{\ast}=0.44~M_{\oplus}$} for a = 0.1 AU.

3.5. Expected populations

Figure 8 shows the relationship between Mthrs (not Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}) and the radius that the planet with Mthrs would have at 10 Gyr without mass loss (solid line). We call this radius the threshold radius, Rthrs. We have calculated Rthrs for Xwt,0 = 100%, 75%, 50%, 25%, 10%, 5%, and 1%. In addition, the mass-radius relationships for rocky planets (dashed line) and pure-water planets (dotted line) at 0.1 AU are also drawn in Fig. 8. There are four characteristic regions in Fig. 8:

  • I

    Planets must contain components less dense than water, such ashydrogen/helium.

  • II

    Planets with water envelopes and without H/He can exist. The water envelopes survive photo-evaporative mass loss.

  • III

    Primordial water envelopes experience significant photo-evaporative mass loss in 10 Gyr.

  • IV

    Planets retain no water envelopes and are composed of rock and iron.

Only in the region II, the planet retains its primordial water envelope for 10 Gyr without significant loss. There are minimum values not only of Mthrs but also of Rthrs; the latter is denoted by Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} hereafter. Note that Rthrs\hbox{$R_{\mathrm{thrs}}^{\ast}$} is not an initial radius.

Those minimum values are helpful to discuss whether planets can possess water components or not, because the uncertainty in water mass fractions can be removed. Since Mthrs and Rthrs depend on semi-major axis, we also compare those threshold values with observed M − a and R − a relationships in the next section.

thumbnail Fig. 8

Relationship between the threshold mass and the threshold radius. The latter is defined by the radius that the planet with Mthrs would have at 10 Gyr without ever experiencing mass loss ( denoted by Rthrs). The squares, which are connected with a solid line, are Mthrs and Rthrs for 0.1 AU and seven different initial water mass fractions Xwt,0 = 100%, 75%, 50%, 25%, 10%, 5%, 1%, and 0.5%. The dashed and dotted lines represent mass-radius relationships, respectively, for rocky planets and pure-water planets at 0.1 AU. Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} and Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} represent the minimum values of Mthrs and Rthrs, respectively.

4. Implications for distributions of observed exoplanets

Figure 9 compares the relationship between the threshold mass, Mthrs, and threshold radius, Rthrs with measured masses and radii of super-Earths around G-type stars identified so far. Here we show three theoretical relationships for a = 0.02, 0.05, and 0.1 AU. As discussed above, only planets on the right side of the theoretical line (i.e., in region II) for a given a are able to retain their water envelopes without significant loss for 10 Gyr.

For future characterizations, planets in region III would be of special interest, because our results suggest that planets should be rare in region III. Three out of the 14 planets, 55 Cnc e, Kepler-20 b, and CoRoT-7 b might be in region III, although errors and the uncertainty in ε (see also the lower panel of Fig. 10 for the sensitivity of Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} to ε) are too large to conclude so. There are at least three possible scenarios for the origin of planets in region III. One is that those planets are halfway to complete evaporation of their water envelopes. Namely, some initial conditions happen to make planets in region III, although such conditions are rare. The second possible scenario is that those planets had formed far from and migrated toward their host stars recently. The third is that those planets are in balance between degassing from the rocky core and the atmospheric escape. Thus, deeper understanding of the properties of those super-Earths via future characterization will provide important constraints on their origins.

thumbnail Fig. 9

Relationship between the threshold mass Mthrs and radius Rthrs (lines; see text for definitions) compared to masses and radii of observed transiting super-Earths around G-type stars (points with error bars; exoplanets.org (Wright et al. 2011), as of June 29, 2013, ). The dotted (blue), dashed (green), and solid (red) lines represent the Mthrs and Rthrs relationships for orbital periods of 11 days (=0.1 AU), 4 days (=0.05 AU), and 1 day (=0.02 AU), respectively. The dash-dotted (brown) line represents the planet composed of rocks. Note that black points represent planets whose orbital periods are longer than 11days. In those calculations, we have assumed the heating efficiency ε = 0.1 and the initial luminosity L0 = 1 × 1024  erg  s-1. “CoR” are short for CoRoT and “Kep” are short for Kepler.

thumbnail Fig. 10

Upper panel: theoretical distribution of masses and semi-major axes (or incident fluxes) of planets at 10 Gyr with various initial masses and water mass fractions. Crosses (red) represent planets that lost their water envelopes completely in 10 Gyr, while open squares (blue) represent planets that survive significant loss of the water envelopes via photo-evaporation. The green line is the minimum threshold masses, Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}. Here, we have adopted ε = 0.1. Lower panel: distribution of masses and semi-major axes (or incident fluxes) of detected exoplanets compared to the minimum threshold mass, Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}, derived in this study (see Sect. 3.3 for definition). We have shown three Mthrsa\hbox{$M_{\rm{thrs}}^{\ast}{-}a$} relationships for different heating efficiencies: ε = 1 (solid line), ε = 0.1 (dashed line), and ε = 0.01 (dotted line). Filled circles with error bars represent observational data (from http://exoplanet.org (Wright et al. 2011)) for planets orbiting host stars with effective temperature of 5000–6000 K (relatively early K-type stars and G-type stars). Planets are colored according to their zero-albedo equilibrium temperatures in K. In the planet names, “CoR” and “Kep” stand for CoRoT and Kepler, respectively.

thumbnail Fig. 11

Upper panel: theoretical distribution of radii and semi-major axes (or incident fluxes) of planets at 10 Gyr with various initial masses and water mass fractions. Crosses (red) represent planets that lost their water envelopes completely due to the photo-evaporation in 10 Gyr, while open squares (blue) represent planets that survive significant loss of the water envelopes. The green line is the minimum threshold radii, Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$}. Here, we have adopted ε = 0.1. Lower panel: distribution of radii and semi-major axes (or incident fluxes) of Kepler planetary candidates, compared to the threshold radius, Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} (see Sect. 3.3 for definition). We have shown three Rthrsa\hbox{$R_{\rm{thrs}}^{\ast}{-}a$} relationships for different heating efficiencies: ε = 1 (solid red line), ε = 0.1 (dashed green line), and ε = 0.01 (dotted blue line). Filled squares represent observational data (http://kepler.nasa.gov, as of June 29, 2013) for planets orbiting host stars with effective temperature of 5300–6000 K (G-type stars).

In this study, low-mass exoplanets, whose masses are ≤20  M and radii ≤4  R, are of special interest. (We call them super-Earths below.) While there are only 14 super-Earths whose masses and radii were both measured (see Fig. 9), the minimum masses (Mpsini) and the orbital periods were measured for about 22 super-Earths around G-type stars (see Fig. 10). Also, over 1000 sub/super-Earth-sized planet candidates have been identified by the Kepler space telescope (Batalha et al. 2013). The size and semi-major axis distribution of those objects is known. It is, thus, interesting to compare our theoretical prediction with the observed Mpa and Rpa distributions.

Before doing so, we demonstrate that Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} and Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} are good indicators for constraining the limits below which evolved planets retain no water envelopes. Figure 10a and 11a show the theoretical distributions of masses and radii of planets that evolved for 10 Gyr, starting with various initial water mass fractions and planetary masses (i.e., Xwt,0 = 25,50,75 and 100% and log (Mp,0/M) = −1 + 0.1j with j = 0,1,···   ,21). The crosses (red) and open squares (blue) represent the planets that lost their water envelopes completely (i.e., rocky planets) and those which survive significant loss of their water envelopes, respectively. As seen in these figures, two populations of rocky planets and water-rich planets are clearly separated by the Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} and Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} lines. Note that there are some planets that retain their water envelope below the threshold line. These planets just retain ≲1% water mass fraction at 10 Gyr. However, such planets are found to be obviously rare.

In Fig. 10b, we show the distribution of Mpsini and a of low-mass exoplanets detected around G-type and K-type stars so far, as compared with Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} for three choices of ε. Among them, α Cen B b, Kepler-10 b and CoRoT-7 b are well below the Mthrs\hbox{$M_\mathrm{thrs}^\ast$} line for ε = 0.1. Thus, the three planets are likely to be rocky, provided ε = 0.1. However, the uncertainty in ε (and FXUV) prevents us from deriving a robust conclusion. An order-of-magnitude difference in ε is found to change Mthrs\hbox{$M_\mathrm{thrs}^\ast$} by a factor of three. The aforementioned three planets are between the two Mthrs\hbox{$M_\mathrm{thrs}^{\ast}$} lines for ε = 0.01 and 0.1. This demonstrates quantitatively how important determining ε and FXUV more accurately is for understanding the composition of super-Earths only with measured masses. It would be worth mentioning that few planets are found between the lines for ε = 0.1 and ε = 1. Since all the planets in Fig. 10b were found by the radial-velocity method, the apparent gap would be unlikely to be due to observational bias. Thus, the gap might suggest that the actual Mthrs\hbox{$M_\mathrm{thrs}^\ast$} line lies between those two ones.

In Fig. 11b, we show the distribution of Rp and a of KOIs, which is compared with Rthrs\hbox{$R_\mathrm{thrs}^{\ast}$} for three choices of ε. Many planets are found to be below the Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} lines. We are unable to constrain the fraction of rocky planets quantitatively, because of the uncertainty in ε. However, since there are many points below the Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} line for ε of as small as 0.01, it seems to be a robust conclusion that KOIs contain a significant number of rocky planets. Note that the distribution must include rocky planets that were formed rocky without ever experiencing mass loss. This means that there are more rocky planets in reality than we have predicted in this study.

As mentioned in the Introduction, Lopez & Fortney (2013) performed a similar investigation of the threshold mass and radius concerning H + He atmospheres on rocky super-Earths (see Figs. 8 and 9 of Lopez & Fortney 2013). For the horizontal axis, they adopted the incident stellar flux, instead of semi-major axis. In Figs. 10 and 11, we have also indicated another scale of the incident flux calculated from the relationship between the semi-major axis a and the incident flux F, F=Lstar4πa2=FEarth(LstarL)(a1AU)-2,\begin{equation} F = \frac{L_{\mathrm{star}}}{4\pi a^2} = F_{\mathrm{Earth}}\left( \frac{L_{\mathrm{star}}}{L_{\odot}} \right) \left( \frac{a}{1\mathrm{~AU}}\right)^{-2}, \end{equation}(36)where Lstar is the luminosity of the host star and FEarth is the current bolometric flux that the Earth receives from the Sun. Comparing their results for the H + He envelope, we find that the threshold value of the initial mass (or incident flux) for H2O is smaller by a factor of about 10 than that for H + He although a similar linear dependence is found. For example, the threshold mass for H + He is ~30   M (derived by Eq. (6) of Lopez & Fortney 2013) in the case of F = 103F, while it is for H2O is ~2   M.

In Fig. 9 of Lopez & Fortney (2013), it has also been suggested that the frequency of planets with radii of 1.8–4.0 R for Fp ≥ 100  F (corresponding to a ≤ 0.1  AU) should be low as a consequence of photo-evaporative mass loss. Owen & Wu (2013) also found a deficit of planets around 2  R in their planet distribution (see Fig. 8 of Owen & Wu 2013). In contrast, our results suggest that water-rich planets with radii of 1.5–3.0 R are relatively common, because they are able to sustain their water envelopes against photo-evaporation. This seeming disagreement on the predicted distribution demonstrates the influence of the envelope composition on the predicted distribution. Indeed, there are many KOIs found in such a domain in the Rp-a diagram shown in Fig. 11a. Thus, those KOIs may be water-rich planets, although it is also possible that they are rocky planets without ever experiencing mass loss.

Finally, we focused in this study on the thermal escape of the upper atmosphere due to stellar XUV irradiation. In addition, ion pick-up induced by stellar winds and coronal mass ejections may be effective in stripping off atmospheres of close-in planets, as discussed for close-in planets with hydrogen-rich atmospheres (e.g. Lammer et al. 2013). Such non-thermal effects lead to increase in Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}. This implies that the Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} obtained in this study is a lower limit on survival of water-rich planets.

5. Summary

In this study, we have investigated the impact of photo-evaporative mass loss on masses and radii of water-rich sub/super-Earths with short orbital periods around G-type stars. We simulated the interior structure and the evolution of highly-irradiated sub/super-Earths that consist of a rocky core surrounded by a water envelope, including the effect of mass loss due to the stellar XUV-driven energy-limited hydrodynamic escape (see Sect. 2).

The findings from this study are summarized as follows. In Sect. 3, we have investigated the mass evolution of water-rich sub/super-Earths, and then found a threshold planet mass Mthrs, below which the planet has its water envelope stripped off in 1–10 Gyr (Sect. 3.1). The initial planet’s luminosity has little impact on Mthrs (Sect. 3.2). We have found that there is a minimum value, Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}, for given a and ε (Sect. 3.4). Water-rich planets with initial masses smaller than Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} lose their water envelopes completely in 10 Gyr, independently of initial water mass fraction. The threshold radius, Rthrs, is defined as the radius that the planet of mass Mthrs would have at 10 Gyr if it evolved without undergoing mass loss. We have also found that there is a minimum value of the threshold radius, Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} (Sect. 3.5). Finally, we have discussed the composition of observed exoplanets in Sect. 4 by comparing the threshold values to measured masses and radii of the exoplanets. Then, we have confirmed quantitatively that more accurate determination of planet masses and radii, ϵ and FXUV, respectively is needed for deriving robust prediction for planetary composition. Nevertheless, the comparison between Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} and radii of KOIs in the Rp − a plane suggests that KOIs contain a significant number of rocky planets.

In this study, we have demonstrated that photo-evaporative mass loss has a significant impact on the evolution of water envelopes of sub/super-Earths, especially with short orbital periods, and that of H + He envelopes of super-Earths. Since the Mthrs for water envelope models is larger by a factor of 10, relative to that for H + He envelope models by Lopez & Fortney (2013), the stability limit for water envelopes gives more robust constraints on the detectability of rocky planets. Thus, the Mthrs and Rthrs will provide valuable information for future searches of rocky Earth-like planets.

Online material

Appendix A: Atmospheric model

First, we describe opacity models for the water vapor atmosphere. We define the Planck-type (κP) and the Rosseland-type mean opacities (κr) as κvp=visibleκνBν(T)dν/visibleBν(T)dν,1κvr=κthp=thermalκνBν(Tatm)dν/thermalBν(Tatm)dν,1κthr=thermal1κνdBν(Tatm)dTdν/thermaldBν(Tatm)dTdν,\appendix \setcounter{section}{1} \begin{eqnarray} \kappa_{\mathrm{v}}^{\mathrm{p}} &=& \int_\mathrm{visible} \kappa_{\nu} B_{\nu}(T_{\star}) {\rm d}\nu~\bigg/ \int_\mathrm{visible} B_{\nu}(T_{\star}) {\rm d}\nu, \\[1mm] \frac{1}{ \kappa_{\mathrm{v}}^{\mathrm{r}}} &=& \int_\mathrm{visible} \frac{1}{\kappa_{\nu}}\frac{{\rm d}B_{\nu}(T_{\star})}{{\rm d}T} {\rm d}\nu ~\bigg/ \int_\mathrm{visible} \frac{{\rm d}B_{\nu}(T_{\star})}{{\rm d}T} {\rm d}\nu \label{kapvr}, \\[1mm] \kappa_{\mathrm{th}}^{\mathrm{p}} &=& \int_\mathrm{thermal} \kappa_{\nu} B_{\nu}(T_{\mathrm{atm}}) {\rm d}\nu~\bigg/ \int_\mathrm{thermal} B_{\nu}(T_{\mathrm{atm}}) {\rm d}\nu, \\[1mm] \frac{1}{\kappa_{\mathrm{th}}^{\mathrm{r}}} &=& \int_\mathrm{thermal} \frac{1}{\kappa_{\nu}}\frac{{\rm d}B_{\nu}(T_{\mathrm{atm}})}{{\rm d}T} {\rm d}\nu~\bigg/ \int_\mathrm{thermal} \frac{{\rm d}B_{\nu}(T_{\mathrm{atm}})}{{\rm d}T} {\rm d}\nu, \end{eqnarray}where ν is the frequency; κν the monochromatic opacity at a given ν; T the stellar effective temperature; Tatm the atmospheric temperature of the planet; and Bν the Planck function. The subscripts, “th” and “v”, mean opacities in the thermal and visible wavelengths, respectively. In this study, we assume T = 5780 K. We adopt HITRAN opacity data for water (Rothman et al. 2009) and calculate mean opacities for 1000 K, 2000 K, and 3000 K at 1, 10, 100 bar. Mean opacities are fitted to power-law functions of P and T, using the least squares method; κvp=1.94×104(P1bar)0.01(T1000K)1.0cm2g-1,κvr=2.20(P1bar)1.0(T1000K)-0.4cm2g-1,κthp=4.15×105(P1bar)0.01(T1000K)-1.1cm2g-1,κthr=3.07×102(P1bar)0.9(T1000K)-4.0cm2g-1,\appendix \setcounter{section}{1} \begin{eqnarray} \label{kap1}\kappa_{\mathrm{v}}^{\mathrm{p}} &=& 1.94\times10^4 \left(\frac{P}{1\mathrm{~bar}} \right)^{0.01} \left( \frac{T}{1000\mathrm{~K}} \right)^{1.0} \mathrm{cm}^2~\mathrm{g}^{-1}, \\[1mm] \label{kap2}\kappa_{\mathrm{v}}^{\mathrm{r}} &=& 2.20 \left(\frac{P}{1\mathrm{~bar}} \right)^{1.0} \left( \frac{T}{1000\mathrm{~K}} \right)^{-0.4} \mathrm{cm}^2~\mathrm{g}^{-1}, \\[1mm] \label{kap3}\kappa_{\mathrm{th}}^{\mathrm{p}} &=& 4.15\times10^{5} \left(\frac{P}{1\mathrm{~bar}} \right)^{0.01} \left(\frac{T}{1000\mathrm{~K}} \right)^{-1.1} \mathrm{cm}^2~\mathrm{g}^{-1}, \\[1mm] \label{kap4}\kappa_{\mathrm{th}}^{\mathrm{r}} &=&3.07 \times 10^{2} \left(\frac{P}{1\mathrm{~bar}} \right)^{0.9} \left(\frac{T}{1000\mathrm{~K}} \right)^{-4.0} \mathrm{cm}^2~\mathrm{g}^{-1}, \end{eqnarray}where P is the pressure and T the temperature.

In this study, we basically follow the prescription developed by Guillot (2010) except for the treatment of the opacity. We consider a static, plane-parallel atmosphere in local thermodynamic equilibrium. We assume that the atmosphere is in radiative equilibrium between an incoming visible flux from the star and an outgoing infrared flux from the planet. Thus, the radiation energy equation and radiation momentum equation are written as dHvdm=κvpJv,dKvdm=κvrHv,dHthdm=κthp(JthB),dKthdm=κthrHth,\appendix \setcounter{section}{1} \begin{eqnarray} \label{RMT-01}\frac{{\rm d}H_{\mathrm{v}}}{{\rm d}m} &=& \kappa_{\mathrm{v}}^{\mathrm{p}}J_{\mathrm{v}}, \\[1mm] \label{RMT-02}\frac{{\rm d}K_{\mathrm{v}}}{{\rm d}m} &=& \kappa_{\mathrm{v}}^{\mathrm{r}}H_{\mathrm{v}}, \\[1mm] \label{RMT-03}\frac{{\rm d}H_{\mathrm{th}}}{{\rm d}m} &=& \kappa_{\mathrm{th}}^{\mathrm{p}}\left( J_{\mathrm{th}} - B \right), \\[1mm] \label{RMT-04}\frac{{\rm d}K_{\mathrm{th}}}{{\rm d}m} &=& \kappa_{\mathrm{th}}^{\mathrm{r}}H_{\mathrm{th}}, \end{eqnarray}and the atmosphere in radiative equilibrium satisfies κvpJv+κthp(JthB)=0,\appendix \setcounter{section}{1} \begin{eqnarray} \kappa_{\mathrm{v}}^{\mathrm{p}}J_{\mathrm{v}} + \kappa_{\mathrm{th}}^{\mathrm{p}}\left( J_{\mathrm{th}} -B\right) = 0, \label{RMT-05} \end{eqnarray}(A.13)where Jv (Jth), Hv (Hth), and Kv (Kth) are, respectively, the zeroth-, first-, and second-order moments of radiation intensity in the visible (thermal) wavelengths, m the atmospheric mass coordinate, dm = ρdz, where z is the altitude from the bottom of the atmosphere, ρ the density, and B the frequency-integrated Planck function, BthermalBνdν~σπT4,\appendix \setcounter{section}{1} \begin{equation} B \equiv \int_\mathrm{thermal} B_\nu {\rm d}\nu \sim \frac{\sigma}{\pi}T^4, \end{equation}(A.14)where σ is the Stefan-Boltzmann constant. We assume here that thermal emission from the atmosphere at visible wavelengths are negligible, so that Bν ~ 0 in the visible region. The six moments of the radiation field are defined as (Jv,Hv,Kv)visible(Jν,Hν,Kν)dν,(Jth,Hth,Kth)thermal(Jν,Hν,Kν)dν,\appendix \setcounter{section}{1} \begin{eqnarray} (J_\mathrm{v}, H_\mathrm{v}, K_\mathrm{v}) &\equiv& \int_\mathrm{visible} (J_\nu, H_\nu, K_\nu) {\rm d}\nu, \\ (J_\mathrm{th}, H_\mathrm{th}, K_\mathrm{th}) &\equiv& \int_\mathrm{thermal} (J_\nu, H_\nu, K_\nu) {\rm d}\nu, \end{eqnarray}where Jν is the mean intensity, 4πHν the radiation flux, and 4πKν/c the radiation pressure (c is the speed of light).

We integrate three moments of specific intensity, Jν,Hν and Kν, over all the frequencies: J0Jνdν=120dν-11dμIν,μ=Jv+Jth,H0Hνdν=120dν-11dμIν,μμ=Hv+Hth,K0Kνdν=120dν-11dμIν,μμ2=Kv+Kth,\appendix \setcounter{section}{1} \begin{eqnarray} J &\equiv& \int_0^{\infty} J_\nu {\rm d}\nu = \frac{1}{2} \int_0^{\infty} {\rm d}\nu \int_{-1}^{1} {\rm d}\mu I_{\nu,\mu} = J_{\mathrm{v}} + J_{\mathrm{th}}, \\ H &\equiv& \int_0^{\infty} H_\nu {\rm d}\nu = \frac{1}{2} \int_0^{\infty} {\rm d}\nu \int_{-1}^{1} {\rm d}\mu I_{\nu,\mu}\mu = H_{\mathrm{v}} + H_{\mathrm{th}}, \\ K &\equiv& \int_0^{\infty} K_\nu {\rm d}\nu = \frac{1}{2} \int_0^{\infty} {\rm d}\nu \int_{-1}^{1} {\rm d}\mu I_{\nu,\mu}\mu^2 = K_{\mathrm{v}} + K_{\mathrm{th}}, \end{eqnarray}where Iν,μ is the specific intensity and θ the angle of a intensity with respect to the z-axis, μ = cosθ. The energy conservation of the total flux implies H=Hv+Hth=14πσTint4,\appendix \setcounter{section}{1} \begin{equation} H = H_{\mathrm{v}} + H_{\mathrm{th}} = \frac{1}{4\pi} \sigma T_{\mathrm{int}}^4, \label{AT01a} \end{equation}(A.20)where Tirr is the irradiation temperature given by Tirr=TRa,\appendix \setcounter{section}{1} \begin{equation} T_{\mathrm{irr}} = T_{\star}\sqrt{\frac{R_{\star}}{a}}, \end{equation}(A.21)where R is the radius of the host star and a the semi-major axis.

For the closure relations, we use the Eddington approximation (e.g. Chandrasekhar 1960), namely, Kv=13Jv,Kth=13Jth.\appendix \setcounter{section}{1} \begin{eqnarray} K_{\mathrm{v}} &=& \frac{1}{3}J_{\mathrm{v}}, \\ K_{\mathrm{th}} &=& \frac{1}{3}J_{\mathrm{th}}. \end{eqnarray}For an isotropic case of both the incoming and outgoing radiation fields, we find boundary conditions of the moment equations as follows (see also Guillot 2010, for details): Hv(m=0)=1314πσTirr4,Hv(m=0)=13Jv(m=0),Hth(m=0)=12Jth(m=0).\appendix \setcounter{section}{1} \begin{eqnarray} \label{AT01b}H_{\mathrm{v}}(m=0) &=& -\frac{1}{\sqrt{3}}\frac{1}{4\pi} \sigma T_{\mathrm{irr}}^4, \\ \label{AT01c}H_{\mathrm{v}}(m=0) &=& -\frac{1}{\sqrt{3}} J_{\mathrm{v}}(m=0), \\ \label{AT01d}H_{\mathrm{th}}(m=0) &=& \frac{1}{2} J_{\mathrm{th}}(m=0). \end{eqnarray}Thus, we integrate Eqs. (A.9)–(A.13) over m numerically, using mean opacities of (A.5)–(A.8) and boundary conditions of (A.24)–(A.26), and then determine a TP profile of the water vapor atmosphere. We assume that the boundary is at P0 = 1 × 10-5 bar. The choice of P0 (≤1 × 10-5bar) has little effect on the atmospheric temperature-pressure structure. T0 is determined in an iterative fashion until abs(T0 −  [πB(m = 0,P0,T0)/σ1/4) ≤ 0.01 is fulfilled. Then we integrate Eqs. (A.9)–(A.13) over m by the 4th-order Runge-Kutta method, until we find the point where dlnT/dlnP ≥ ∇ad. The pressure and temperature, Pad and Tad, are the boundary conditions for the convective-interior structure (see Sect. 2.1).

In Fig. A.1, we show the PT profile for the solar-composition atmosphere with g = 980 cm s-2, Tint = 300 K, and Tirr = 1500 K (dotted line). In this calculation, we take κthr\hbox{$\kappa_{\mathrm{th}}^{\mathrm{r}}$} and κthp\hbox{$\kappa_{\mathrm{th}}^{\mathrm{p}}$} as functions of P and T from Freedman et al. (2008) and calculate κvp\hbox{$\kappa_\mathrm{v}^\mathrm{p}$} and κvr\hbox{$\kappa_\mathrm{v}^\mathrm{r}$}, for P = 1 × 10-3,0.1,1,10 bar and T = 1500 K from HITRAN and HITEMP data that include H2, He, H2O, CO, CH4, Na, and K for the solar abundance respectively as κv={\appendix \setcounter{section}{1} \begin{eqnarray} \kappa_{\mathrm{v}} = \left\{ \begin{array}{llclc} 1.51\times 10^{-5} & \mathrm{cm}^2 \, \mathrm{g}^{-1} & (10^{-3} &\le P[\mathrm{bar}]),& \\[2mm] 3.88\times 10^{-4} & \mathrm{cm}^2 \, \mathrm{g}^{-1} & (10^{-3} &< P[\mathrm{bar}] \le& 10^{-1}), \\[2mm] 3.05\times 10^{-3} & \mathrm{cm}^2 \, \mathrm{g}^{-1} & (10^{-1} &< P[\mathrm{bar}] \le& 1), \\[2mm] 2.65\times 10^{-2} & \mathrm{cm}^2 \, \mathrm{g}^{-1} & (1 &> P[\mathrm{bar}]), & \\ \end{array} \right. \end{eqnarray}(A.27)by use of (A.2). The thin and thick parts of the dotted line represent the radiative and convective zones, respectively.

In addition, we test our atmosphere model by comparing it with the PT profile derived by Guillot (2010) with γ = κv/κth = 0.4 (solid line), which reproduces more detailed atmosphere models by Fortney et al. (2005) and Iro et al. (2005, see Fig. 6 of Guillot 2010). As seen in Fig. A.1, our atmospheric model yields a PT profile similar to that from Guillot (2010). In our model, temperatures are relatively low compared with the Guillot (2010) model at P ≲ 40 bar, which is due to difference in opacity. In our model, deep regions of P ≳ 40 bar are convective, while there is no convective region in the Guillot (2010) model because of constant opacity. We have compared our PT profile with the Fortney et al. (2005)’s and Iro et al. (2005)’s profiles, which are shown in Fig. 6 of Guillot (2010) and confirmed that our PT profile in the convective region is almost equal to their profiles. Of special interest in this study is the entropy at the radiative/convective boundary, because it governs the thermal evolution of the planet. In this sense, it is fair to say that our atmospheric model yields appropriate boundary conditions for the structure of the convective interior.

thumbnail Fig. A.1

Temperature–pressure profiles for a solar-composition atmosphere (see the details in text). The solid (red) and dotted (green) lines represent both Guillot (2010)’s (γ = 0.4) and our model’s, respectively. The thin and thick parts of the dotted line represent the radiative and convective regions, respectively. We have assumed g = 980 cm s-2, Tint = 300 K, and Tirr = 1500 K.

Finally, we describe an analytical expression for our atmospheric model. We basically follow the prescription developed by Heng et al. (2012), except for the treatment of the opacity. As Heng et al. (2012) mentioned, it would be a challenging task without assumption of constant κvp\hbox{$\kappa_{\mathrm{v}}^{\mathrm{p}}$} and κvr\hbox{$\kappa_{\mathrm{v}}^{\mathrm{r}}$} to obtain analytical solutions for Jv and Hv. Here we assume κvp\hbox{$\kappa_{\mathrm{v}}^{\mathrm{p}}$} and κvr\hbox{$\kappa_{\mathrm{v}}^{\mathrm{r}}$} are constant throughout the atmosphere. We differentiate (A.9) and (A.10) by m and obtain d2Jvdm2=Hvμ2dκvrdm+κvrκvpμ2Jv,d2Hvdm2=Jvdκvpdm+κvrκvpμ2Hv,\appendix \setcounter{section}{1} \begin{eqnarray} \label{RMT-01-1}\frac{{\rm d}^2 J_{\mathrm{v}}}{{\rm d}m^2} &=& \frac{H_{\mathrm{v}}}{\mu^2} \frac{ {\rm d}\kappa_{\mathrm{v}}^{\mathrm{r}}}{{\rm d}m} + \frac{\kappa_{\mathrm{v}}^{\mathrm{r}}\kappa_{\mathrm{v}}^{\mathrm{p}}}{\mu^2} J_{\mathrm{v}}, \\[2mm] \label{RMT-02-1}\frac{{\rm d}^2 H_{\mathrm{v}}}{{\rm d}m^2} &=&J_{\mathrm{v}} \frac{ {\rm d}\kappa_{\mathrm{v}}^{\mathrm{p}}}{{\rm d}m} + \frac{\kappa_{\mathrm{v}}^{\mathrm{r}}\kappa_{\mathrm{v}}^{\mathrm{p}}}{\mu^2} H_{\mathrm{v}}, \end{eqnarray}where μ2 = Kv/Jv. Assuming Jv = Hv = 0 as m → ∞, we obtain (Jv,Hv)=(Jv,0,Hv,0)exp(κv¯μm),\appendix \setcounter{section}{1} \begin{eqnarray} (J_{\mathrm{v}},~H_{\mathrm{v}}) = (J_{\mathrm{v},0},~H_{\mathrm{v},0})\exp\left(-\frac{\bar{\kappa_{\mathrm{v}}}}{\mu}m \right), \end{eqnarray}(A.30)where κv¯=κvpκvr\hbox{$\bar{\kappa_{\mathrm{v}}} = \sqrt{ \kappa_{\mathrm{v}}^{\mathrm{p}}\kappa_{\mathrm{v}}^{\mathrm{r}} }$} and Jv,0 and Hv,0 are the values of Jv and Hv evaluated at m = 0, respectively. In general, the heat transportation, such as circulation, produces a specific luminosity of heat. Heng et al. (2012) introduced the specific luminosity as Q, which has units of erg s-1 g-1. Q can be related to the moments of the specific intensity and we obtain κthp(JthB)+κvpJv=Q.\appendix \setcounter{section}{1} \begin{equation} \kappa_{\mathrm{th}}^{\mathrm{p}}\left(J_{\mathrm{th}} - B \right) + \kappa_{\mathrm{v}}^{\mathrm{p}} J_{\mathrm{v}} = Q. \label{RMT-05-1} \end{equation}(A.31)We integrate Eq. (A.31) and obtain H=H˜Q(m,),\appendix \setcounter{section}{1} \begin{equation} H = H_{\infty} - \tilde{Q}(m,\infty) \label{RMT-05-2}, \end{equation}(A.32)where H is the value of H evaluated at m → ∞ and ˜Q(m1,m2)=m1m2Q(m,μ,φ)dm.\appendix \setcounter{section}{1} \begin{equation} \tilde{Q}(m_1,m_2) = \int_{m_1}^{m_2} Q(m',\mu,\phi) {\rm d}m'. \end{equation}(A.33)To obtain Hth and Jth, we substitute Eq. (A.31) in Eqs. (A.11) and (A.12) and integrate by m. Then we obtain Hth=HHv,0exp(κv¯μm)˜Q(m,)Jth=Jth,0Hv,0fKth0mκthrexp(κv¯μm)dm+1fKth0mκthr{H˜Q(m,)}dm,\appendix \setcounter{section}{1} \begin{eqnarray} H_{\mathrm{th}}&= & H_{\mathrm{\infty}} -H_{\mathrm{v},0} \exp\left( -\frac{\bar{\kappa_{\mathrm{v}}} }{\mu} m \right) - \tilde{Q}(m,\infty) \\ J_{\mathrm{th}} &=& J_{\mathrm{th},0} - \frac{H_{\mathrm{v},0}}{ f_{K\mathrm{th}}} \int_{0}^{m} \kappa_{\mathrm{th}}^{\mathrm{r}} \exp\left(-\frac{\bar{\kappa_{\mathrm{v}}}}{\mu} m' \right) {\rm d}m' \nonumber \\ &&+ \frac{1}{f_{K\mathrm{th}}} \int_{0}^{m} \kappa_{\mathrm{th}}^{\mathrm{r}} \left\{ H_{\mathrm{\infty}} - \tilde{Q} (m', \infty) \right\} {\rm d}m', \end{eqnarray}where fKth = Kth/Jth, fHth = Hth/Jth, and Jth,0=1fHth{HHv,0˜Q(0,)}.\appendix \setcounter{section}{1} \begin{equation} J_{\mathrm{th},0} = \frac{1}{f_{H\mathrm{th}}}\left\{ H_{\infty}-H_{\mathrm{v},0}-\tilde{Q}(0,\infty) \right\}. \end{equation}(A.36)That is, we obtain B=H[1fHth+1fKthτth(m)]Hv,0[1fHth+κv¯μκthp+1fKthτext(m)]+E(m),\appendix \setcounter{section}{1} \begin{eqnarray} B &=& H_{\infty}\left[ \frac{1}{f_{H\mathrm{th}}} + \frac{1}{f_{K\mathrm{th}}} \tau_{\mathrm{th}}(m) \right] \nonumber \\ &&-H_{\mathrm{v},0} \left[ \frac{1}{f_{H\mathrm{th}}} +\frac{\bar{\kappa_{\mathrm{v}}}}{\mu \kappa_{\mathrm{th}}^{\mathrm{p}}} +\frac{1}{f_{K\mathrm{th}}} \tau_{\mathrm{ext}}(m) \right] + E(m), \end{eqnarray}(A.37)where τth(m)=m0κthrdm,τext(m)=m0(κth¯2fKthμ2κv¯2)1κthpexp(κv¯μm)dm,E(m)=[Qκthp+1fKth0mκthr˜Q(m,)dm+˜Q(0,)fHth],\appendix \setcounter{section}{1} \begin{eqnarray} \tau_{\mathrm{th}}(m) &=& \int_{0}^{m}\kappa_{\mathrm{th}}^{\mathrm{r}} {\rm d}m', \\ \tau_{\mathrm{ext}}(m) &=& \int_{0}^{m} \left( \bar{\kappa_{\mathrm{th}}}^2 - \frac{ f_{K\mathrm{th}} }{\mu^2} \bar{\kappa_{\mathrm{v}}}^2 \right) \frac{1}{\kappa_{\mathrm{th}}^{\mathrm{p}}}\exp\left( -\frac{\bar{\kappa_{\mathrm{v}}}}{\mu}m' \right) {\rm d}m' , \\ E(m) & = & -\left[ \frac{Q}{\kappa_{\mathrm{th}}^{\mathrm{p}}} + \frac{1}{f_{K\mathrm{th}}} \int_{0}^{m} \kappa_{\mathrm{th}}^{\mathrm{r}} \tilde{Q}(m',\infty){\rm d}m' + \frac{\tilde{Q}(0,\infty)}{f_{H\mathrm{th}}} \right],\quad\quad\quad \end{eqnarray}and κth¯=κthpκthr\hbox{$\bar{\kappa_{\mathrm{th}}} = \sqrt{\kappa_{\mathrm{th}}^{\mathrm{p}}\kappa_{\mathrm{th}}^{\mathrm{r}}}$}. In our conditions, we assume μ=1/3\hbox{$\mu = 1/\sqrt{3}$}, fKth = 1/3, fHth = 1/2 and Q = 0. Consequently, we obtain

the temperature profile as T4=34Tint4[23+τth(m)]+34Tirr4[23+κv¯3κthp+τext(m)]\appendix \setcounter{section}{1} \begin{equation} T^4 = \frac{3}{4}T_{\mathrm{int}}^4 \left[ \frac{2}{3} + \tau_{\mathrm{th}}(m) \right] + \frac{\sqrt{3}}{4} T_{\mathrm{irr}}^4 \left[ \frac{2}{3} + \frac{\bar{\kappa_{\mathrm{v}}}}{\sqrt{3}\kappa_{\mathrm{th}}^{\mathrm{p}}} + \tau_{\mathrm{ext}}(m) \right] \label{TP_ana} \end{equation}(A.41)whereτext(m)=0mκth¯2κv¯2κthpexp(3κv¯m)dm.\appendix \setcounter{section}{1} \begin{eqnarray} \tau_{\mathrm{ext}}(m) &=& \int_{0}^{m} \frac{\bar{\kappa_{\mathrm{th}}}^2 - \bar{\kappa_{\mathrm{v}}}^2}{\kappa_{\mathrm{th}}^{\mathrm{p}}} \exp\left(-\sqrt{3} \bar{\kappa_{\mathrm{v}}} m' \right) {\rm d}m'. \end{eqnarray}(A.42)If we assume κthp=κthr\hbox{$\kappa_{\mathrm{th}}^{\mathrm{p}} = \kappa_{\mathrm{th}}^{\mathrm{r}} $} and κvp=κvr\hbox{$\kappa_{\mathrm{v}}^{\mathrm{p}} = \kappa_{\mathrm{v}}^{\mathrm{r}}$}, Eq. (A.41) agrees with Eq. (27) of Heng et al. (2012).

Acknowledgments

We thank N. Nettelmann for providing us with tabulated data for the equation of state of water (H2O-EOS) and S. Ida and T. Guillot for fruitful advice and discussions. We also thank the anonymous referee for his/her careful reading and constructive comments that helped us improve this paper greatly. We also thank Y. Ito and Y. Kawashima for providing us with the opacity data and fruitful suggestions about the atmospheric structure. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. This study is supported by Grants-in-Aid for Scientific Research on Innovative Areas (No. 23103005) and Scientific Research (C) (No. 25400224) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. K. K. is supported by a grant for the Global COE Program, “From the Earth to “Earths””, of MEXT, Japan. Y. H. is supported by the Grant-in-Aid for JSPS Fellows (No. 23003491) from MEXT, Japan.

References

  1. Adams, E. R., Seager, S., & Elkins-Tanton, L. 2008, ApJ, 673, 1160 [NASA ADS] [CrossRef] [Google Scholar]
  2. Barclay, T., Rowe, J. F., Lissauer, J. J., et al. 2013, Nature, 494, 452 [NASA ADS] [CrossRef] [Google Scholar]
  3. Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24 [NASA ADS] [CrossRef] [Google Scholar]
  4. Chandrasekhar, S. 1960, Radiative transfer (New York: Dover) [Google Scholar]
  5. Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69 [NASA ADS] [CrossRef] [Google Scholar]
  7. Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
  8. Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
  9. French, M., Mattsson, T. R., Nettelmann, N., & Redmer, R. 2009, Phys. Rev. B, 79, 054107 [NASA ADS] [CrossRef] [Google Scholar]
  10. Fressin, F., Torres, G., Désert, J.-M., et al. 2011, ApJS, 197, 5 [NASA ADS] [CrossRef] [Google Scholar]
  11. Grasset, O., Schneider, J., & Sotin, C. 2009, ApJ, 693, 722 [NASA ADS] [CrossRef] [Google Scholar]
  12. Guillot, T. 2010, A&A, 520, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Hansen, B. M. S. 2008, ApJS, 179, 484 [NASA ADS] [CrossRef] [Google Scholar]
  14. Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, MNRAS, 420, 20 [NASA ADS] [CrossRef] [Google Scholar]
  15. Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
  16. Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66 [NASA ADS] [CrossRef] [Google Scholar]
  17. Iro, N., Bézard, B., & Guillot, T. 2005, A&A, 436, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Kane, S. R. 2007, MNRAS, 380, 1488 [NASA ADS] [CrossRef] [Google Scholar]
  19. Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution, Astron. Astrophys. Lib. (Berlin, Heidelberg, New York: Springer-Verlag) [Google Scholar]
  20. Lammer, H., Erkaev, N. V., Odert, P., et al. 2013, MNRAS, 430, 1247 [NASA ADS] [CrossRef] [Google Scholar]
  21. Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lyon, S., & Johnson, J. D. 1992, LANL Tech. Rep. LA-UR-92-3407 (Los Alamos: LANL) [Google Scholar]
  24. Muirhead, P. S., Hamren, K., Schlawin, E., et al. 2012, ApJ, 750, L37 [NASA ADS] [CrossRef] [Google Scholar]
  25. Nettelmann, N., Holst, B., Kietzmann, A., et al. 2008, ApJ, 683, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  26. Nettelmann, N., Fortney, J. J., Kramm, U., & Redmer, R. 2011, ApJ, 733, 2 [Google Scholar]
  27. Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [NASA ADS] [CrossRef] [Google Scholar]
  28. Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
  29. Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680 [NASA ADS] [CrossRef] [Google Scholar]
  30. Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
  31. Rothman, L. S., Gordon, I. E., Barbe, A., et al. 2009, J. Quant. Spectr. Rad. Transf., 110, 533 [NASA ADS] [CrossRef] [Google Scholar]
  32. Seager, S., Kuchner, M., Hier-Majumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279 [NASA ADS] [CrossRef] [Google Scholar]
  33. Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [NASA ADS] [CrossRef] [Google Scholar]
  34. Swift, D. C., Eggert, J. H., Hicks, D. G., et al. 2012, ApJ, 744, 59 [NASA ADS] [CrossRef] [Google Scholar]
  35. Valencia, D., Sasselov, D. D., & O’Connell, R. J. 2007, ApJ, 656, 545 [Google Scholar]
  36. Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10 [NASA ADS] [CrossRef] [Google Scholar]
  38. Wagner, F. W., Sohl, F., Hussmann, H., Grott, M., & Rauer, H. 2011, Icarus, 214, 366 [NASA ADS] [CrossRef] [Google Scholar]
  39. Ward, W. R. 1986, Icarus, 67, 164 [NASA ADS] [CrossRef] [Google Scholar]
  40. Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150 [NASA ADS] [CrossRef] [Google Scholar]
  41. Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412 [NASA ADS] [CrossRef] [Google Scholar]
  42. Yelle, R., Lammer, H., & Ip, W.-H. 2008, Space Sci. Rev., 139, 437 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Model of the planetary structure in this study.

In the text
thumbnail Fig. 2

Concept of the chord optical depth.

In the text
thumbnail Fig. 3

Mass evolution of close-in water-rich planets. The solid blue lines represent planets that retain their water envelopes for 10 Gyr. In contrast, the planet shown by the dashed red line loses its water envelope completely in 10 Gyr. We set Lp,0 = 1 × 1024 erg s-1, Xwt,0 = 75%, a = 0.1 AU, and ε = 0.1 for all the planets. In this model, we assume that the rocky core never evaporates.

In the text
thumbnail Fig. 4

Threshold mass in M as a function of the initial planet’s luminosity in erg s-1 for four choices of semi-major axes. The solid (red), dashed (green), dotted (blue), and dot-dashed (purple) lines represent a = 0.02,0.03,0.05, and 0.1 AU, respectively. We have assumed Xwt = 75% and ε = 0.1.

In the text
thumbnail Fig. 5

Relationship between the initial planetary mass and the fraction of the water envelope at 10 Gyr for four initial water mass fractions of Xwt,0 = 100% (solid, red), 75% (dashed, green), 50% (dotted, blue), and 25% (dot-dashed, purple). We have assumed L0 = 1 × 1024  erg  s-1, a = 0.1 AU, and ε = 0.1.

In the text
thumbnail Fig. 6

Relationship between the initial water mass fraction Xwt,0 in % and the threshold mass Mthrs in M for four choices of semi-major axes of 0.02 AU (solid, red), 0.03 AU (dashed, green), 0.05 AU (dotted, blue), and 0.1 AU (dot-dashed, purple). We have assumed L0 = 1 × 1024  erg  s-1 and ε = 0.1.

In the text
thumbnail Fig. 7

Relationship between the initial planetary mass and the fraction of the initial water envelope that is lost via photo-evaporation in 5 Gyr for six initial water mass fractions of Xwt,0 = 1% (solid, red), 3% (long-dashed, green), 10% (dotted, blue), 30% (dash-dotted, purple), 50% (dot-dashed, light blue), and 60% (dashed, black). We have assumed L0 = 1 × 1024  erg  s-1, a = 0.1 AU, and ε = 0.1.

In the text
thumbnail Fig. 8

Relationship between the threshold mass and the threshold radius. The latter is defined by the radius that the planet with Mthrs would have at 10 Gyr without ever experiencing mass loss ( denoted by Rthrs). The squares, which are connected with a solid line, are Mthrs and Rthrs for 0.1 AU and seven different initial water mass fractions Xwt,0 = 100%, 75%, 50%, 25%, 10%, 5%, 1%, and 0.5%. The dashed and dotted lines represent mass-radius relationships, respectively, for rocky planets and pure-water planets at 0.1 AU. Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$} and Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} represent the minimum values of Mthrs and Rthrs, respectively.

In the text
thumbnail Fig. 9

Relationship between the threshold mass Mthrs and radius Rthrs (lines; see text for definitions) compared to masses and radii of observed transiting super-Earths around G-type stars (points with error bars; exoplanets.org (Wright et al. 2011), as of June 29, 2013, ). The dotted (blue), dashed (green), and solid (red) lines represent the Mthrs and Rthrs relationships for orbital periods of 11 days (=0.1 AU), 4 days (=0.05 AU), and 1 day (=0.02 AU), respectively. The dash-dotted (brown) line represents the planet composed of rocks. Note that black points represent planets whose orbital periods are longer than 11days. In those calculations, we have assumed the heating efficiency ε = 0.1 and the initial luminosity L0 = 1 × 1024  erg  s-1. “CoR” are short for CoRoT and “Kep” are short for Kepler.

In the text
thumbnail Fig. 10

Upper panel: theoretical distribution of masses and semi-major axes (or incident fluxes) of planets at 10 Gyr with various initial masses and water mass fractions. Crosses (red) represent planets that lost their water envelopes completely in 10 Gyr, while open squares (blue) represent planets that survive significant loss of the water envelopes via photo-evaporation. The green line is the minimum threshold masses, Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}. Here, we have adopted ε = 0.1. Lower panel: distribution of masses and semi-major axes (or incident fluxes) of detected exoplanets compared to the minimum threshold mass, Mthrs\hbox{$M_{\rm{thrs}}^{\ast}$}, derived in this study (see Sect. 3.3 for definition). We have shown three Mthrsa\hbox{$M_{\rm{thrs}}^{\ast}{-}a$} relationships for different heating efficiencies: ε = 1 (solid line), ε = 0.1 (dashed line), and ε = 0.01 (dotted line). Filled circles with error bars represent observational data (from http://exoplanet.org (Wright et al. 2011)) for planets orbiting host stars with effective temperature of 5000–6000 K (relatively early K-type stars and G-type stars). Planets are colored according to their zero-albedo equilibrium temperatures in K. In the planet names, “CoR” and “Kep” stand for CoRoT and Kepler, respectively.

In the text
thumbnail Fig. 11

Upper panel: theoretical distribution of radii and semi-major axes (or incident fluxes) of planets at 10 Gyr with various initial masses and water mass fractions. Crosses (red) represent planets that lost their water envelopes completely due to the photo-evaporation in 10 Gyr, while open squares (blue) represent planets that survive significant loss of the water envelopes. The green line is the minimum threshold radii, Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$}. Here, we have adopted ε = 0.1. Lower panel: distribution of radii and semi-major axes (or incident fluxes) of Kepler planetary candidates, compared to the threshold radius, Rthrs\hbox{$R_{\rm{thrs}}^{\ast}$} (see Sect. 3.3 for definition). We have shown three Rthrsa\hbox{$R_{\rm{thrs}}^{\ast}{-}a$} relationships for different heating efficiencies: ε = 1 (solid red line), ε = 0.1 (dashed green line), and ε = 0.01 (dotted blue line). Filled squares represent observational data (http://kepler.nasa.gov, as of June 29, 2013) for planets orbiting host stars with effective temperature of 5300–6000 K (G-type stars).

In the text
thumbnail Fig. A.1

Temperature–pressure profiles for a solar-composition atmosphere (see the details in text). The solid (red) and dotted (green) lines represent both Guillot (2010)’s (γ = 0.4) and our model’s, respectively. The thin and thick parts of the dotted line represent the radiative and convective regions, respectively. We have assumed g = 980 cm s-2, Tint = 300 K, and Tirr = 1500 K.

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.