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/00046361/201322258  
Published online  10 February 2014 
Impact of photoevaporative mass loss on masses and radii of waterrich sub/superEarths ^{⋆}
^{1}
Department of Earth and Planetary ScienceThe University of
Tokyo,
731 Hongo, Bunkyoku, 1130033
Tokyo,
Japan
email:
kkurosaki@eps.s.utokyo.ac.jp
^{2}
Division of Theoretical Astronomy, National Astronomical
Observatory of Japan, 2211 Osawa,
Mitaka, 1818588
Tokyo,
Japan
Received:
11
July
2013
Accepted:
3
December
2013
Context. Recent progress in transit photometry opened a new window to the interior of superEarths. From measured radii and masses, we can infer constraints on planetary internal compositions. It has been recently revealed that superEarths orbiting close to host stars (i.e., hot superEarths) are diverse in composition. This diversity is thought to arise from diversity in volatile content.
Aims. The stability of the volatile components, which we call the envelopes, is to be examined, because hot superEarths, which are exposed to strong irradiation, undergo photoevaporative mass loss. While several studies investigated the impact of photoevaporative mass loss on hydrogenhelium envelopes, there are few studies as to the impact on watervapor envelopes, which we investigate in this study. To obtain theoretical prediction to future observations, we also investigate the relationships among masses, radii, and semimajor axes of waterrich superEarths and also subEarths that have undergone photoevaporative mass loss.
Methods. We simulate the interior structure and evolution of highlyirradiated sub/superEarths that consist of a rocky core surrounded by a water envelope, which include mass loss due to the stellar XUVdriven energylimited hydrodynamic escape.
Results. We find that the photoevaporative mass loss has a significant impact on the evolution of hot sub/superEarths. With a widelyused empirical formula for XUV flux from typical Gstars and the heating efficiency of 0.1 for example, the planets of less than 3 Earth masses orbiting 0.03 AU have their water envelopes completely stripped off. We then derive the threshold planetary mass and radius below which the planet loses its water envelope completely as a function of the initial water content and find that there are minimums of the threshold mass and radius.
Conclusions. We constrain the domain in the parameter space of planetary mass, radius, and the semimajor axis in which sub/superEarths never retain water envelopes in 1–10 Gyr. This would provide an essential piece of information for understanding the origin of closein, lowmass planets. The current uncertainties in stellar XUV flux and its heating efficiency, however, prevent us from deriving robust conclusions. Nevertheless, it seems to be a robust conclusion that Kepler planet candidates contain a significant number of rocky sub/superEarths.
Key words: planets and satellites: composition / planets and satellites: interiors
Appendix A is available in electronic form at http://www.aanda.org
© 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 smallsized exoplanets with radii of 1 to 2 R_{⊕} have been identified, which are often referred to as superEarths (Batalha et al. 2013). Also, planet candidates detected by the Kepler space telescope include subEarthsized objects, such as Kepler20 e (Fressin et al. 2011), Kepler42 b, c, d (Muirhead et al. 2012), and Kepler37 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 superEarths detected so far show a large variation in radius, suggesting diversity in composition. There are many theoretical studies on massradius 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 lowdensity superEarths that are larger in size than they would be if they were rocky. This implies that these transiting superEarths possess components less dense than rock. From a viewpoint of planet formation, the possible components are hydrogenrich gas and water, which make an outer envelope. A small fraction of Hrich gas or water is known to be enough to account for observed radii of the lowdensity superEarths (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 semimajor axis (e.g., Kane 2007). These closein planets are highly irradiated and exposed to intense Xray 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 photoevaporation of planetary envelopes. As for massive closein planets, namely, hot Jupiters, the possibility of the photoevaporation and its outcome have been investigated well both theoretically and observationally (e.g., Yelle et al. 2008 and references therein).
While the photoevaporation may not significantly affect the evolution and final composition of hot Jupiters except for extremely irradiated or inflated hot Jupiters, its impact on small closein planets in the sub/superEarth 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 superEarth CoRoT7 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 CoRoT7 b is unlikely to retain the H + He envelope at present.
Young mainsequence 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 superEarth had a primordial atmosphere initially, it may lose the atmosphere completely during its history. These discussions concerning the photoevaporative loss of H + He envelopes were done for GJ 1214 b (Nettelmann et al. 2011; Valencia et al. 2013), superEarths orbiting Kepler11 (Lopez et al. 2012; Ikoma & Hori 2012), and CoRoT7 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 photoevaporation on the stability of H + He envelopes for superEarths. In particular, Lopez & Fortney (2013) performed simulations of coupled thermal contraction and photoevaporative mass loss of rocky superEarths 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 semimajor axes and planetary radii (or planet densities) of KOIs.
In this study, we focus on waterrich sub/superEarths. Planet formation theories predict that lowmass planets migrate toward their host star, which is strongly supported by the presence of many closein superEarths, from cooler regions (e.g., Ward 1986) where they may have accreted a significant amount of water. This suggests that water/icerich sub/superEarths may also exist close to host stars. Therefore, similar discussions should be done for water envelopes of closein superEarths. However, there are just a few studies, which treat specific sub/superEarths such as CoRoT7 b (Valencia et al. 2010) and Kepler11 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 closein sub/superEarths against photoevaporation. To this end, we simulate the thermal evolution of planets with significant fractions of water envelopes (i.e., waterworlds), incorporating the effect of stellarXUVdriven photoevaporative 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 waterrich planets. Then, we find threshold values of planetary masses and radii below which such waterrich planets are incapable of retaining primordial water envelopes for a period similar to ages of known exoplanethost stars (i.e., 1–10 Gyr). In Sect. 4, we compare the theoretical massradius distribution of waterrich 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.
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 photoevaporation. 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), $\begin{array}{ccc}\frac{\mathit{\partial P}}{\mathit{\partial}{\mathit{M}}_{\mathit{r}}}& \mathrm{=}& \mathrm{}\frac{\mathit{G}{\mathit{M}}_{\mathit{r}}}{\mathrm{4}\mathit{\pi}{\mathit{r}}^{\mathrm{4}}}\mathit{,}\\ \frac{\mathit{\partial r}}{\mathit{\partial}{\mathit{M}}_{\mathit{r}}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}{\mathit{r}}^{\mathrm{2}}\mathit{\rho}}\mathit{,}\\ \frac{\mathit{\partial T}}{\mathit{\partial}{\mathit{M}}_{\mathit{r}}}& \mathrm{=}& \mathrm{}\frac{\mathit{G}{\mathit{M}}_{\mathit{r}}\mathit{T}}{\mathrm{4}\mathit{\pi}{\mathit{r}}^{\mathrm{4}}\mathit{P}}\mathrm{\nabla}\mathit{,}\end{array}$and the equation of state, $\begin{array}{ccc}\mathit{\rho}\mathrm{=}\mathit{\rho}\mathrm{\left(}\mathit{P,T}\mathrm{\right)}\mathit{,}& & \end{array}$(4)where r is the planetocentric distance, M_{r} 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 cm^{2} 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, $\mathrm{\nabla}\mathrm{=}{\left(\frac{\mathit{\partial}\mathrm{ln}\mathit{T}}{\mathit{\partial}\mathrm{ln}\mathit{P}}\right)}_{\mathit{S}}\mathrm{\xb7}$(5)Equations (1)–(3) require three boundary conditions. The inner one is r = 0 at M_{r} = 0. The outer boundary corresponds to the interface between the envelope and the atmosphere, which is called the tropopause. The tropopause pressure P_{ad} and temperature T_{ad} 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 M_{p}. In our calculation, the atmospheric mass is less than 0.1% of the planetary mass. Thus, the outer boundary conditions are given as $\begin{array}{c}\\ \mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{ad}}& \mathrm{and}& \mathit{T}\mathrm{=}{\mathit{T}}_{\mathrm{ad}}& \mathrm{at}& {\mathit{M}}_{\mathit{r}}\mathrm{=}{\mathit{M}}_{\mathrm{p}}\mathit{.}\end{array}$(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) $\frac{\mathit{\partial L}}{\mathit{\partial}{\mathit{M}}_{\mathit{r}}}\mathrm{=}\mathrm{}\mathit{T}\frac{\mathit{\partial S}}{\mathit{\partial t}}\mathit{,}$(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 $\mathrm{}{\mathit{L}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{\partial}\mathit{S\u0305}\mathrm{e}}{\mathit{\partial t}}{\mathrm{\int}}_{{\mathit{M}}_{\mathrm{c}}}^{{\mathit{M}}_{\mathrm{p}}}\mathit{T}\mathrm{d}{\mathit{M}}_{\mathit{r}}\mathrm{+}\frac{\mathit{\partial}\mathit{S\u0305}\mathrm{c}}{\mathit{\partial t}}{\mathrm{\int}}_{\mathrm{0}}^{{\mathit{M}}_{\mathrm{c}}}\mathit{T}\mathrm{d}{\mathit{M}}_{\mathit{r}}\mathit{,}$(8)where L_{p} is the total intrinsic luminosity of the planet, M_{c} is the mass of the rocky core, and and are the specific entropies in the water envelope and the rocky core, respectively. In integrating Eq. (7), we have assumed L = 0 at M_{r} = 0.
In the numerical calculations of this study, we use the intrinsic temperature T_{int}, instead of L_{p}, which is defined by ${\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{\equiv}\frac{{\mathit{L}}_{\mathrm{p}}}{\mathrm{4}\mathit{\pi}{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}\mathit{\sigma}}\mathit{,}$(9)where R_{p} is the planet photospheric radius (see Sect. 2.4 for the definition) and σ is the StefanBoltzmann 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 H_{2}O. In practice, we use the numerical code developed by Hori & Ikoma (2011), which calculate chemical equilibrium compositions among H_{2}O, H_{2}, O_{2}, 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 H_{2}OREOS (Nettelmann et al. 2008), which contains the ab initio water EOS data at high pressures of French et al. (2009). H_{2}OREOS covers a density range from 1.0 × 10^{6} g cm^{3} to 15 g cm^{3} and a temperature range from 1.0 × 10^{3} K to 2.4 × 10^{4} K. For T and ρ outside the ranges that H_{2}OREOS 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 widelyused 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, radiativeequilibrium 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 semigrey, planeparallel 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 firstorder 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 solarcomposition (i.e., hydrogendominated) 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 watervapor 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 (P_{ad}, T_{ad}) are used as the outer boundary conditions for the structure of the convective water envelope.
The photospheric radius R_{p} used in Eq. (9) is the radius at which the thermal optical depth measured from infinity, τ, is 2/3; namely, $\mathit{\tau}\mathrm{=}{\mathrm{\int}}_{{\mathit{R}}_{\mathrm{p}}}^{\mathrm{\infty}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}\mathit{\rho}\mathrm{d}\mathit{r}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathit{,}$(10)where ${\mathit{\kappa}}_{\mathrm{th}}^{\mathit{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 R_{conv} (see Fig. 1). We evaluate the atmospheric thickness z (=R_{p} − R_{conv}) by integrating the hydrostatic equation from P = P_{ad} to P = P_{ph} using $\mathit{z}\mathrm{=}\mathrm{}{\mathrm{\int}}_{{\mathit{P}}_{\mathrm{ad}}}^{{\mathit{P}}_{\mathrm{ph}}}\frac{\mathrm{d}\mathit{P}}{\mathit{g\rho}}\mathrm{=}\mathrm{}{\mathrm{\int}}_{{\mathit{P}}_{\mathrm{ad}}}^{{\mathit{P}}_{\mathrm{ph}}}\frac{\mathrm{\mathcal{R}}}{\mathit{\mu g}}\frac{\mathit{T}}{\mathit{P}}\mathrm{d}\mathit{P,}$(11)where g is the constant gravity, ℛ (=8.31 × 10^{7} erg K^{1} g^{1}) is the gas constant, and μ is the mean molecular weight. P_{ph} is the photospheric pressure that we calculate by integrating $\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{\tau}}\mathrm{=}\frac{\mathit{g}}{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}}$(12)from τ = 0 to 2/3.
2.5. Mass loss
The mass loss is assumed to occur in an energylimited fashion. Its rate, including the effect of the Roche lobe, is given by (Erkaev et al. 2007) $\mathit{M\u0307}\mathrm{=}\mathrm{}\frac{\mathit{\epsilon}{\mathit{F}}_{\mathrm{XUV}}{\mathit{R}}_{\mathrm{p}}\mathit{\pi}{\mathit{R}}_{\mathrm{XUV}}^{\mathrm{2}}}{\mathit{G}{\mathit{M}}_{\mathrm{p}}{\mathit{K}}_{\mathrm{tide}}}\mathit{,}$(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; F_{XUV} is the incident flux of Xray and UV radiation from the host star, K_{tide} is the potential energy reduction factor due to stellar tide; and R_{XUV} is the effective radius at which the planet receives the incident XUV flux. In Eq. (13), we have assumed R_{XUV} = R_{p}, which is a good approximation for closein planets of interest (Lammer et al. 2013). It is noted that Lammer et al. (2013) focused on the hydrogenhelium atmosphere. Since the scale height of the vapor atmosphere is smaller than that of a hydrogenhelium atmosphere with the same temperature, R_{XUV} ≃ R_{p} is a good approximation also for the vapor atmosphere.
In this study, we suppose that the host star is a Gstar and adopt the empirical formula derived by Ribas et al. (2005) for F_{XUV}: $\begin{array}{ccc}{\mathit{F}}_{\mathrm{XUV}}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(14)\normalsizeWe use the formula for K_{tide} derived by Erkaev et al. (2007), ${\mathit{K}}_{\mathrm{tide}}\mathrm{=}\frac{\mathrm{(}\mathit{\eta}\mathrm{}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{2}\mathit{\eta}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{2}{\mathit{\eta}}^{\mathrm{3}}}\mathit{,}$(15)where η is the ratio of the Rochelobe (or Hill) radius to the planetary radius, R_{p}.
The value of the heating efficiency is uncertain, because minor gases such as CO_{2} contribute to it via radiative cooling. For photoevaporation of hotJupiters, ε 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 planeparallel 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.
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 ${\mathit{\tau}}_{\mathrm{ch}}\mathrm{\left(}\mathit{r,\nu}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{}\mathrm{\infty}}^{\mathrm{+}\mathrm{\infty}}\mathit{\rho}{\mathit{\kappa}}_{\mathit{\nu}}\mathrm{d}\mathit{s,}$(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, R_{tr}, as ${\mathit{\tau}}_{\mathrm{ch}}\mathrm{\left(}{\mathit{R}}_{\mathrm{tr}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{\xb7}$(17)Let the altitude from the sphere of radius r be z_{tr}. Then s^{2} = (r + z_{tr})^{2} − r^{2} (Fig. 2). Eq. (16) is written as ${\mathit{\tau}}_{\mathrm{ch}}\mathrm{\left(}\mathit{r,\nu}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{\rho}{\mathit{\kappa}}_{\mathit{\nu}}\frac{{\mathit{z}}_{\mathrm{tr}}\mathrm{+}\mathit{r}}{\sqrt{{\mathit{z}}_{\mathrm{tr}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{r}{\mathit{z}}_{\mathrm{tr}}}}\mathrm{d}{\mathit{z}}_{\mathrm{tr}}\mathit{.}$(18)Furthermore for convenience, we choose pressure P as the independent variable, instead of z_{tr}. Using the equation of hydrostatic equilibrium, $\frac{\mathrm{d}\mathit{P}}{\mathrm{d}{\mathit{z}}_{\mathrm{tr}}}\mathrm{=}\mathrm{}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{\rho}}{\mathrm{(}\mathit{r}\mathrm{+}{\mathit{z}}_{\mathrm{tr}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(19)one obtains ${\mathit{\tau}}_{\mathrm{ch}}\mathrm{\left(}\mathit{\nu ,r}\mathrm{\right)}\mathrm{=}\mathrm{}\frac{\mathrm{2}}{{\mathit{g}}_{\mathit{r}}}{\mathrm{\int}}_{{\mathit{P}}_{\mathit{r}}}^{\mathrm{0}}{\mathit{\kappa}}_{\mathit{\nu}}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathrm{tr}}\mathit{/}\mathit{r}{\mathrm{)}}^{\mathrm{3}}}{\sqrt{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathrm{tr}}\mathit{/}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\mathrm{}\mathrm{1}}}\mathrm{d}\mathit{P,}$(20)where ${\mathit{g}}_{\mathit{r}}\mathrm{=}\frac{\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}$(21)and P_{r} is the pressure at r. To integrate Eq. (20), we write z_{tr} as a function of P. To do so, we integrate Eq. (19) and obtain ${\mathrm{\int}}_{\mathrm{0}}^{{\mathit{z}}_{\mathrm{tr}}}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime}}}{\mathrm{(}\mathit{r}\mathrm{+}{\mathit{z}}^{\mathrm{\prime}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{=}\mathrm{}{\mathrm{\int}}_{{\mathit{P}}_{\mathit{r}}}^{{\mathit{P}}_{\mathit{z}}}\frac{\mathrm{d}\mathit{P}}{\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{\rho}}\mathit{,}$(22)where P_{z} is the pressure at z_{tr}. Eq. (22) is integrated as $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{r}\mathrm{+}{\mathit{z}}_{\mathrm{tr}}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}\mathrm{}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}{\mathit{g}}_{\mathit{r}}}{\mathrm{\int}}_{{\mathit{P}}_{\mathit{z}}}^{{\mathit{P}}_{\mathit{r}}}\frac{\mathrm{d}\mathit{P}}{\mathit{\rho}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}\mathrm{}\frac{{\mathit{z}}_{\mathrm{p}}\mathrm{\left(}{\mathit{P}}_{\mathit{r}}\mathit{,}{\mathit{P}}_{\mathit{z}}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}\end{array}$(23)where ${\mathit{z}}_{\mathrm{p}}\mathrm{\left(}{\mathit{P}}_{\mathit{r}}\mathit{,}{\mathit{P}}_{\mathit{z}}\mathrm{\right)}\mathrm{\equiv}{\mathrm{\int}}_{{\mathit{P}}_{\mathit{z}}}^{{\mathit{P}}_{\mathit{r}}}\frac{\mathit{P}}{\mathit{\rho}{\mathit{g}}_{\mathit{r}}}\mathrm{d}\mathrm{ln}\mathit{P}\mathit{.}$(24)Thus, z is written as ${\mathit{z}}_{\mathrm{tr}}\mathrm{=}{\mathit{z}}_{\mathrm{p}}{\left(\mathrm{1}\mathrm{}\frac{{\mathit{z}}_{\mathrm{p}}}{\mathit{r}}\right)}^{1}\mathit{.}$(25)Note that z_{p} corresponds to the altitude in the case of a planeparallel atmosphere and (1 − z_{p}/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 M_{p}(t) and an assumed M_{p}(t + Δt). To be exact, the two structures are integrated for two different values of T_{int}. In doing so, we integrate Eqs. (1)–(4) inward from the tropopause to the planetary center, using the fourthorder RungeKutta method. The inward integration is started with the outer boundary condition given by Eq. (6); P_{ad} and T_{ad} 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 M_{r} = 0) in an iterative fashion. Note that determining P_{ad} and T_{ad} requires the gravity in the atmosphere (or R_{conv}), 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 secondorder difference equation for Eq. (8), which is written as $\begin{array}{ccc}\mathrm{\Delta}\mathit{t}& \mathrm{=}& \mathrm{}(\left[\mathit{S\u0305}\mathrm{e}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{}\mathit{S\u0305}\mathrm{e}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right]\left[{\mathrm{\Theta}}_{\mathrm{e}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{+}{\mathrm{\Theta}}_{\mathrm{e}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right]\\ & & \mathrm{+}\left[\mathit{S\u0305}\mathrm{c}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{}\mathit{S\u0305}\mathrm{c}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right]\left[{\mathrm{\Theta}}_{\mathrm{c}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{+}{\mathrm{\Theta}}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right])\mathrm{\times}{\mathrm{(}{\mathit{L}}_{\mathrm{p}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{+}{\mathit{L}}_{\mathrm{p}}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{)}}}^{1}\mathit{,}\end{array}$(26)where ${\mathrm{\Theta}}_{\mathrm{e}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\equiv}{\mathrm{\int}}_{{\mathit{M}}_{\mathrm{c}}}^{{\mathit{M}}_{\mathrm{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}{\mathit{M}}_{\mathit{r}}\mathit{,}{\mathrm{\Theta}}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\equiv}{\mathrm{\int}}_{\mathrm{0}}^{{\mathit{M}}_{\mathrm{c}}}\mathit{T}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}{\mathit{M}}_{\mathit{r}}\mathit{.}$(27)Using this Δt, we integrate Eq. (13) to calculate M_{p}(t + Δt) as ${\mathit{M}}_{\mathrm{p}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta}\mathit{t}\mathrm{)}\mathrm{=}{\mathit{M}}_{\mathrm{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}\mathit{M\u0307}\mathrm{\Delta}\mathit{t}\mathit{.}$(28)The assumed M_{p}(t + Δt) is not always equal to that obtained here. Therefore the entire procedure must be repeated until the M_{p}(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 superEarths 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 closein waterrich planet. The evolution is controlled by the following five parameters: the initial total mass of the planet (M_{p,0}), the initial luminosity (L_{0}), the initial water mass fraction (X_{wt,0}), the semimajor axis (a), and the heating efficiency (ε). Below, we adopt L_{0} = 1 × 10^{24} erg s^{1}, X_{wt,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 closein waterrich planet.
Fig. 3 Mass evolution of closein waterrich 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 L_{p,0} = 1 × 10^{24} erg s^{1}, X_{wt,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 waterrich planets with six different initial masses; L_{0} = 1 × 10^{24} erg s^{1}, X_{wt,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 waterrich 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 ${\mathit{\tau}}_{\mathit{M}}\mathrm{=}\left\frac{{\mathit{X}}_{\mathrm{wt}}{\mathit{M}}_{\mathrm{p}}}{\mathit{M\u0307}\mathrm{p}}\right\mathrm{=}\frac{\mathrm{4}\mathit{G}{\mathit{K}}_{\mathrm{tide}}{\mathit{X}}_{\mathrm{wt}}{\mathit{M}}_{\mathrm{p}}{\mathit{\rho}}_{\mathrm{pl}}}{\mathrm{3}\mathit{\epsilon}{\mathit{F}}_{\mathrm{XUV}}}\mathit{,}$(29)where ρ_{pl} is the mean density of the planet. As the planetary mass decreases, the massloss timescale becomes shorter. This trend is enhanced by the M − ρ relationship that the mean density decreases as M_{p} decreases, according to our numerical results for waterrich planets.
In addition, the timedependence of the stellar XUV flux (see Eq. (14)) is a crucial factor to cause a striking difference in behavior between the lowmass and highmass planets. Using Eq. (14), we obtain the following relation for τ_{M}: $\begin{array}{ccc}{\mathit{\tau}}_{\mathit{M}}\mathrm{\simeq}\{\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$(30)where $\mathit{f}\mathrm{=}\mathrm{1}{\left(\frac{\mathit{a}}{\mathrm{0.1}\mathrm{AU}}\right)}^{\mathrm{2}}\left(\frac{{\mathit{X}}_{\mathrm{wt}}{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{\oplus}}}\right)\left(\frac{{\mathit{\rho}}_{\mathrm{pl}}}{\mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{3}}\right)\left(\frac{{\mathit{K}}_{\mathrm{tide}}}{\mathrm{0.9}}\right){\left(\frac{\mathit{\epsilon}}{\mathrm{0.1}}\right)}^{1}\mathit{.}$(31)Note that 0.1 g cm^{3} is a typical value of ρ_{pl} in the case of subEarthmass planets with the age of 10^{8} 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 M_{thrs}) is 0.56 M_{⊕} for the fiducial parameter set, which is in good agreement with M_{p} < 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. Hydrogenrich planets are more vulnerable to the photoevaporative mass loss than waterrich planets. According to their study, the threshold mass of the hydrogenrich planet at 0.1 AU is ~5 M_{⊕}. That is, M_{thrs} for waterrich planets is smaller by a factor of ~10 than that of hydrogenrich planets.
3.2. Dependence on the initial planet’s luminosity
The evolution during the first 0.1 Gyr determines the fate of a waterrich 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 L_{0} enhances mass loss because of a large planetary radius. On the other hand, L_{0} 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 L_{0}, Fig. 4 shows M_{thrs} as a function of L_{0} for a = 0.02,0.03,0.05 and 0.1 AU. We have found that M_{thrs} is almost independent of L_{0}. This is because an initiallyluminous planet cools down rapidly, so that the integrated amount of water loss during the highluminosity phase is negligible. This is confirmed by the following argument. The mass loss, ΔM, at the early stage can be estimated by $\mathrm{\Delta}\mathit{M}\mathrm{~}\mathit{M\u0307}{\mathit{\tau}}_{\mathrm{KH}}\mathit{,}$(32)where τ_{KH} is the typical timescale of KelvinHelmholtz contraction, ${\mathit{\tau}}_{\mathrm{KH}}\mathrm{\simeq}\frac{\mathit{G}{\mathit{M}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}_{\mathrm{p}}{\mathit{L}}_{\mathrm{p}}}\mathrm{\xb7}$(33)With Eqs. (29) and (33) given, Eq. (32) can be written as $\begin{array}{ccc}\mathrm{\Delta}\mathit{M}& \mathrm{~}& {\mathit{M}}_{\mathrm{p}}\frac{{\mathit{\tau}}_{\mathrm{KH}}}{{\mathit{\tau}}_{\mathit{M}}}\mathrm{=}{\mathit{M}}_{\mathrm{p}}\frac{\mathit{\epsilon}}{\mathrm{2}{\mathit{K}}_{\mathrm{tide}}}\mathrm{\xb7}\frac{\mathit{\pi}{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{XUV}}}{{\mathit{L}}_{\mathrm{p}}}\\ & \mathrm{~}& \mathrm{3}\mathrm{\times}{\mathrm{10}}^{2}\left(\frac{{\mathit{F}}_{\mathrm{XUV}}}{\mathrm{504}\mathrm{erg}{\mathrm{cm}}^{2}{\mathrm{s}}^{1}}\right)\left(\frac{\mathit{\epsilon}}{\mathrm{0.1}}\right){\left(\frac{{\mathit{K}}_{\mathrm{tide}}}{\mathrm{0.9}}\right)}^{1}\\ & & \mathrm{\times}{\left(\frac{\mathit{a}}{\mathrm{0.1}\mathrm{AU}}\right)}^{2}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{\mathrm{3}{\mathit{R}}_{\mathrm{\oplus}}}\right)}^{\mathrm{2}}{\left(\frac{{\mathit{L}}_{\mathrm{p}}}{{\mathrm{10}}^{\mathrm{24}}\mathrm{erg}{\mathrm{s}}^{1}}\right)}^{1}{\mathit{M}}_{\mathrm{p}}\mathit{.}\end{array}$Because F_{XUV} is constant in the early phase, ΔM decreases as L_{p} increases; that is, the KelvinHelmholtz contraction proceeds more rapidly. Therefore, the choice of the value of L_{0} has little effect on the total amount of water loss, as far as L_{0} is larger than 10^{24} erg s^{1}. For smaller L_{0}, R_{p} is insensitive to L_{0}. Thus, M_{thrs} is insensitive to L_{0}.
Fig. 4 Threshold mass in M_{⊕} as a function of the initial planet’s luminosity in erg s^{1} for four choices of semimajor axes. The solid (red), dashed (green), dotted (blue), and dotdashed (purple) lines represent a = 0.02,0.03,0.05, and 0.1 AU, respectively. We have assumed X_{wt} = 75% and ε = 0.1. 
3.3. Dependence on the initial water mass fraction
The fate of a waterrich planet also depends on the initial water mass fraction, X_{wt,0}. Figure 5 shows X_{wt}(t) at t = 10 Gyr as a function of the initial planet’s mass, M_{p,0}, for four different values of X_{wt,0}(=25%, 50%, 75%, and 100%). As M_{p,0} decreases, X_{wt}(10 Gyr) decreases. The pure water planet (solid line) with M_{p,0} < 0.82 M_{⊕} is completely evaporated in 10 Gyr; namely, X_{wt}(10 Gyr) = 0%. Otherwise, X_{wt}(10 Gyr) = 100%. In other cases, we find that the threshold mass, M_{thrs}, below which X_{wt}(10 Gyr) = 0%, is 0.56 M_{⊕} for X_{wt,0} = 75%, 0.44 M_{⊕} for X_{wt,0} = 50%, and 0.44 M_{⊕} for X_{wt,0} = 25%.
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 X_{wt,0} = 100% (solid, red), 75% (dashed, green), 50% (dotted, blue), and 25% (dotdashed, purple). We have assumed L_{0} = 1 × 10^{24} erg s^{1}, a = 0.1 AU, and ε = 0.1. 
Figure 6 shows the relationship between X_{wt,0} and M_{thrs} for four different semimajor axes. M_{thrs} is found not to be a monotonous function of X_{wt,0}. For X_{wt,0} < 25%, M_{thrs} decreases, as X_{wt,0} increases. This is explained as follows. According to Eq. (29), the mass loss timescale , τ_{M}, depends on the absolute amount of water, X_{wt}M_{p}, and the planetary bulk density, ρ_{pl}. When X_{wt} 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., X_{wt}M_{p}). This means that, M_{p} must be larger for τ_{M} to be the same if X_{wt,0} is small. As a consequence, M_{thrs} decreases with increasing X_{wt,0}. More exactly, M_{thrs} changes with X_{wt,0} in such a way that X_{wt,0}M_{thrs} is constant. In contrast, when X_{wt,0} is large, X_{wt}, M_{p}, and ρ_{pl} affect the mass loss timescale. For a given M_{p}, an increase in X_{wt,0} leads to a decrease in ρ_{pl} (or, an increase in radius), which enhances mass loss. As a result, M_{thrs} increases with X_{wt,0} for X_{wt,0} > 25%. Therefore, there is a minimum value of M_{thrs}, which is hereafter denoted by ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$.
Similar trends can be seen in Figs. 3 and 4 of Lopez & Fortney (2013). To compare our results for waterrich planets to those for hydrogenrich 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 photoevaporation in 5 Gyr in Fig. 7 (see Fig. 3c of Lopez & Fortney 2013). We set L_{0} = 1 × 10^{24} erg s^{1}, a = 0.1 AU, ε = 0.1, and six initial water mass fractions of X_{wt,0} = 1% (solid, red), 3% (longdashed, green), 10% (dotted, blue), 30% (dashdotted, purple), 50% (dotdashed, 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 photoevaporation. In addition, the required initial total mass for X_{wt,0} < 10% becomes significantly large in the water case. This behavior is also found in the case of the hydrogenrich planets for X_{wt,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.
Fig. 6 Relationship between the initial water mass fraction X_{wt,0} in % and the threshold mass M_{thrs} in M_{⊕} for four choices of semimajor axes of 0.02 AU (solid, red), 0.03 AU (dashed, green), 0.05 AU (dotted, blue), and 0.1 AU (dotdashed, purple). We have assumed L_{0} = 1 × 10^{24} erg s^{1} and ε = 0.1. 
Fig. 7 Relationship between the initial planetary mass and the fraction of the initial water envelope that is lost via photoevaporation in 5 Gyr for six initial water mass fractions of X_{wt,0} = 1% (solid, red), 3% (longdashed, green), 10% (dotted, blue), 30% (dashdotted, purple), 50% (dotdashed, light blue), and 60% (dashed, black). We have assumed L_{0} = 1 × 10^{24} erg s^{1}, a = 0.1 AU, and ε = 0.1. 
3.4. Dependence on the semimajor axis
At small a, the incident stellar XUV flux becomes large. Thus, M_{thrs} increases, as a decreases. Certainly, the distance to the host star affects the equilibrium temperature T_{eq}, which has an influence on ρ_{pl}: the higher T_{eq} is, the smaller ρ_{pl} is. However, its impact on M_{thrs} is small, relative to that of F_{XUV}. 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 F_{XUV} has a much greater impact on the mass loss than decreasing ρ_{pl}. In Fig. 6, we find ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{=}\mathrm{5.2}{\mathit{M}}_{\mathrm{\oplus}}$ for a = 0.02 AU, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{=}\mathrm{2.5}{\mathit{M}}_{\mathrm{\oplus}}$ for a = 0.03 AU, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{=}\mathrm{1.2}{\mathit{M}}_{\mathrm{\oplus}}$ for a = 0.05 AU, and ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{=}\mathrm{0.44}{\mathit{M}}_{\mathrm{\oplus}}$ for a = 0.1 AU.
3.5. Expected populations
Figure 8 shows the relationship between M_{thrs} (not ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$) and the radius that the planet with M_{thrs} would have at 10 Gyr without mass loss (solid line). We call this radius the threshold radius, R_{thrs}. We have calculated R_{thrs} for X_{wt,0} = 100%, 75%, 50%, 25%, 10%, 5%, and 1%. In addition, the massradius relationships for rocky planets (dashed line) and purewater 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 photoevaporative mass loss.

III
Primordial water envelopes experience significant photoevaporative 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 M_{thrs} but also of R_{thrs}; the latter is denoted by ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ hereafter. Note that ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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 M_{thrs} and R_{thrs} depend on semimajor axis, we also compare those threshold values with observed M − a and R − a relationships in the next section.
Fig. 8 Relationship between the threshold mass and the threshold radius. The latter is defined by the radius that the planet with M_{thrs} would have at 10 Gyr without ever experiencing mass loss ( denoted by R_{thrs}). The squares, which are connected with a solid line, are M_{thrs} and R_{thrs} for 0.1 AU and seven different initial water mass fractions X_{wt,0} = 100%, 75%, 50%, 25%, 10%, 5%, 1%, and 0.5%. The dashed and dotted lines represent massradius relationships, respectively, for rocky planets and purewater planets at 0.1 AU. ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ and ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ represent the minimum values of M_{thrs} and R_{thrs}, respectively. 
4. Implications for distributions of observed exoplanets
Figure 9 compares the relationship between the threshold mass, M_{thrs}, and threshold radius, R_{thrs} with measured masses and radii of superEarths around Gtype 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, Kepler20 b, and CoRoT7 b might be in region III, although errors and the uncertainty in ε (see also the lower panel of Fig. 10 for the sensitivity of ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\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 superEarths via future characterization will provide important constraints on their origins.
Fig. 9 Relationship between the threshold mass M_{thrs} and radius R_{thrs} (lines; see text for definitions) compared to masses and radii of observed transiting superEarths around Gtype 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 M_{thrs} and R_{thrs} relationships for orbital periods of 11 days (=0.1 AU), 4 days (=0.05 AU), and 1 day (=0.02 AU), respectively. The dashdotted (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 L_{0} = 1 × 10^{24} erg s^{1}. “CoR” are short for CoRoT and “Kep” are short for Kepler. 
Fig. 10 Upper panel: theoretical distribution of masses and semimajor 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 photoevaporation. The green line is the minimum threshold masses, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$. Here, we have adopted ε = 0.1. Lower panel: distribution of masses and semimajor axes (or incident fluxes) of detected exoplanets compared to the minimum threshold mass, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$, derived in this study (see Sect. 3.3 for definition). We have shown three ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{}\mathit{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 Ktype stars and Gtype stars). Planets are colored according to their zeroalbedo equilibrium temperatures in K. In the planet names, “CoR” and “Kep” stand for CoRoT and Kepler, respectively. 
Fig. 11 Upper panel: theoretical distribution of radii and semimajor 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 photoevaporation 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, ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$. Here, we have adopted ε = 0.1. Lower panel: distribution of radii and semimajor axes (or incident fluxes) of Kepler planetary candidates, compared to the threshold radius, ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ (see Sect. 3.3 for definition). We have shown three ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{}\mathit{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 (Gtype stars). 
In this study, lowmass exoplanets, whose masses are ≤20 M_{⊕} and radii ≤4 R_{⊕}, are of special interest. (We call them superEarths below.) While there are only 14 superEarths whose masses and radii were both measured (see Fig. 9), the minimum masses (M_{p}sini) and the orbital periods were measured for about 22 superEarths around Gtype stars (see Fig. 10). Also, over 1000 sub/superEarthsized planet candidates have been identified by the Kepler space telescope (Batalha et al. 2013). The size and semimajor axis distribution of those objects is known. It is, thus, interesting to compare our theoretical prediction with the observed M_{p}–a and R_{p}–a distributions.
Before doing so, we demonstrate that ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ and ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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., X_{wt,0} = 25,50,75 and 100% and log (M_{p,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 waterrich planets are clearly separated by the ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ and ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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 M_{p}sini and a of lowmass exoplanets detected around Gtype and Ktype stars so far, as compared with ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ for three choices of ε. Among them, α Cen B b, Kepler10 b and CoRoT7 b are well below the ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ line for ε = 0.1. Thus, the three planets are likely to be rocky, provided ε = 0.1. However, the uncertainty in ε (and F_{XUV}) prevents us from deriving a robust conclusion. An orderofmagnitude difference in ε is found to change ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ by a factor of three. The aforementioned three planets are between the two ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ lines for ε = 0.01 and 0.1. This demonstrates quantitatively how important determining ε and F_{XUV} more accurately is for understanding the composition of superEarths 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 radialvelocity method, the apparent gap would be unlikely to be due to observational bias. Thus, the gap might suggest that the actual ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ line lies between those two ones.
In Fig. 11b, we show the distribution of R_{p} and a of KOIs, which is compared with ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ for three choices of ε. Many planets are found to be below the ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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 ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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 superEarths (see Figs. 8 and 9 of Lopez & Fortney 2013). For the horizontal axis, they adopted the incident stellar flux, instead of semimajor axis. In Figs. 10 and 11, we have also indicated another scale of the incident flux calculated from the relationship between the semimajor axis a and the incident flux F, $\mathit{F}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{star}}}{\mathrm{4}\mathit{\pi}{\mathit{a}}^{\mathrm{2}}}\mathrm{=}{\mathit{F}}_{\mathrm{Earth}}\left(\frac{{\mathit{L}}_{\mathrm{star}}}{{\mathit{L}}_{\mathrm{\odot}}}\right){\left(\frac{\mathit{a}}{\mathrm{1}\mathrm{AU}}\right)}^{2}\mathit{,}$(36)where L_{star} is the luminosity of the host star and F_{Earth} 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 H_{2}O 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 = 10^{3}F_{⊕}, while it is for H_{2}O 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 F_{p} ≥ 100 F_{⊕} (corresponding to a ≤ 0.1 AU) should be low as a consequence of photoevaporative 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 waterrich planets with radii of 1.5–3.0 R_{⊕} are relatively common, because they are able to sustain their water envelopes against photoevaporation. 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 R_{p}a diagram shown in Fig. 11a. Thus, those KOIs may be waterrich 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 pickup induced by stellar winds and coronal mass ejections may be effective in stripping off atmospheres of closein planets, as discussed for closein planets with hydrogenrich atmospheres (e.g. Lammer et al. 2013). Such nonthermal effects lead to increase in ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$. This implies that the ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ obtained in this study is a lower limit on survival of waterrich planets.
5. Summary
In this study, we have investigated the impact of photoevaporative mass loss on masses and radii of waterrich sub/superEarths with short orbital periods around Gtype stars. We simulated the interior structure and the evolution of highlyirradiated sub/superEarths that consist of a rocky core surrounded by a water envelope, including the effect of mass loss due to the stellar XUVdriven energylimited hydrodynamic escape (see Sect. 2).
The findings from this study are summarized as follows. In Sect. 3, we have investigated the mass evolution of waterrich sub/superEarths, and then found a threshold planet mass M_{thrs}, 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 M_{thrs} (Sect. 3.2). We have found that there is a minimum value, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$, for given a and ε (Sect. 3.4). Waterrich planets with initial masses smaller than ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ lose their water envelopes completely in 10 Gyr, independently of initial water mass fraction. The threshold radius, R_{thrs}, is defined as the radius that the planet of mass M_{thrs} 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, ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\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 F_{XUV}, respectively is needed for deriving robust prediction for planetary composition. Nevertheless, the comparison between ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ and radii of KOIs in the R_{p} − a plane suggests that KOIs contain a significant number of rocky planets.
In this study, we have demonstrated that photoevaporative mass loss has a significant impact on the evolution of water envelopes of sub/superEarths, especially with short orbital periods, and that of H + He envelopes of superEarths. Since the M_{thrs} 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 M_{thrs} and R_{thrs} will provide valuable information for future searches of rocky Earthlike planets.
Online material
Appendix A: Atmospheric model
First, we describe opacity models for the water vapor atmosphere. We define the Plancktype (κ^{P}) and the Rosselandtype mean opacities (κ^{r}) as $\begin{array}{ccc}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}& \mathrm{=}& \int \mathrm{visible}{\mathit{\kappa}}_{\mathit{\nu}}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathit{\star}}\mathrm{\right)}\mathrm{d}\mathit{\nu}/{\mathrm{\int}}_{\mathrm{visible}}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathit{\star}}\mathrm{\right)}\mathrm{d}\mathit{\nu ,}\\ \frac{\mathrm{1}}{{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}}& \mathrm{=}& \\ {\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}& \mathrm{=}& \int \mathrm{thermal}{\mathit{\kappa}}_{\mathit{\nu}}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right)}\mathrm{d}\mathit{\nu}/{\mathrm{\int}}_{\mathrm{thermal}}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right)}\mathrm{d}\mathit{\nu ,}\\ \frac{\mathrm{1}}{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}}& \mathrm{=}& {\mathrm{\int}}_{\mathrm{thermal}}\frac{\mathrm{1}}{{\mathit{\kappa}}_{\mathit{\nu}}}\frac{\mathrm{d}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right)}}{\mathrm{d}\mathit{T}}\mathrm{d}\mathit{\nu}/{\mathrm{\int}}_{\mathrm{thermal}}\frac{\mathrm{d}{\mathit{B}}_{\mathit{\nu}}\mathrm{\left(}{\mathit{T}}_{\mathrm{atm}}\mathrm{\right)}}{\mathrm{d}\mathit{T}}\mathrm{d}\mathit{\nu ,}\end{array}$where ν is the frequency; κ_{ν} the monochromatic opacity at a given ν; T_{⋆} the stellar effective temperature; T_{atm} 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 powerlaw functions of P and T, using the least squares method; $\begin{array}{ccc}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}& \mathrm{=}& \mathrm{1.94}\mathrm{\times}{\mathrm{10}}^{\mathrm{4}}{\left(\frac{\mathit{P}}{\mathrm{1}\mathrm{bar}}\right)}^{\mathrm{0.01}}{\left(\frac{\mathit{T}}{\mathrm{1000}\mathrm{K}}\right)}^{\mathrm{1.0}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{g}}^{1}\mathit{,}\\ {\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}& \mathrm{=}& \mathrm{2.20}{\left(\frac{\mathit{P}}{\mathrm{1}\mathrm{bar}}\right)}^{\mathrm{1.0}}{\left(\frac{\mathit{T}}{\mathrm{1000}\mathrm{K}}\right)}^{0.4}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{g}}^{1}\mathit{,}\\ {\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}& \mathrm{=}& \mathrm{4.15}\mathrm{\times}{\mathrm{10}}^{\mathrm{5}}{\left(\frac{\mathit{P}}{\mathrm{1}\mathrm{bar}}\right)}^{\mathrm{0.01}}{\left(\frac{\mathit{T}}{\mathrm{1000}\mathrm{K}}\right)}^{1.1}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{g}}^{1}\mathit{,}\\ {\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}& \mathrm{=}& \mathrm{3.07}\mathrm{\times}{\mathrm{10}}^{\mathrm{2}}{\left(\frac{\mathit{P}}{\mathrm{1}\mathrm{bar}}\right)}^{\mathrm{0.9}}{\left(\frac{\mathit{T}}{\mathrm{1000}\mathrm{K}}\right)}^{4.0}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{g}}^{1}\mathit{,}\end{array}$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, planeparallel 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 $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{H}}_{\mathrm{v}}}{\mathrm{d}\mathit{m}}& \mathrm{=}& {\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}{\mathit{J}}_{\mathrm{v}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{K}}_{\mathrm{v}}}{\mathrm{d}\mathit{m}}& \mathrm{=}& {\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}{\mathit{H}}_{\mathrm{v}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{H}}_{\mathrm{th}}}{\mathrm{d}\mathit{m}}& \mathrm{=}& {\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}\mathrm{(}{\mathit{J}}_{\mathrm{th}}\mathrm{}{\mathit{B}}^{\mathrm{)}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{K}}_{\mathrm{th}}}{\mathrm{d}\mathit{m}}& \mathrm{=}& {\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}{\mathit{H}}_{\mathrm{th}}\mathit{,}\end{array}$and the atmosphere in radiative equilibrium satisfies $\begin{array}{ccc}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}{\mathit{J}}_{\mathrm{v}}\mathrm{+}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}\mathrm{(}{\mathit{J}}_{\mathrm{th}}\mathrm{}{\mathit{B}}^{\mathrm{)}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(A.13)where J_{v} (J_{th}), H_{v} (H_{th}), and K_{v} (K_{th}) are, respectively, the zeroth, first, and secondorder 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 frequencyintegrated Planck function, $\mathit{B}\mathrm{\equiv}{\mathrm{\int}}_{\mathrm{thermal}}{\mathit{B}}_{\mathit{\nu}}\mathrm{d}\mathit{\nu}\mathrm{~}\frac{\mathit{\sigma}}{\mathit{\pi}}{\mathit{T}}^{\mathrm{4}}\mathit{,}$(A.14)where σ is the StefanBoltzmann 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 $\begin{array}{ccc}\mathrm{\left(}{\mathit{J}}_{\mathrm{v}}\mathit{,}{\mathit{H}}_{\mathrm{v}}\mathit{,}{\mathit{K}}_{\mathrm{v}}\mathrm{\right)}& \mathrm{\equiv}& \int \mathrm{visible}\mathrm{\left(}{\mathit{J}}_{\mathit{\nu}}\mathit{,}{\mathit{H}}_{\mathit{\nu}}\mathit{,}{\mathit{K}}_{\mathit{\nu}}\mathrm{\right)}\mathrm{d}\mathit{\nu ,}\\ \mathrm{\left(}{\mathit{J}}_{\mathrm{th}}\mathit{,}{\mathit{H}}_{\mathrm{th}}\mathit{,}{\mathit{K}}_{\mathrm{th}}\mathrm{\right)}& \mathrm{\equiv}& {\mathrm{\int}}_{\mathrm{thermal}}\mathrm{\left(}{\mathit{J}}_{\mathit{\nu}}\mathit{,}{\mathit{H}}_{\mathit{\nu}}\mathit{,}{\mathit{K}}_{\mathit{\nu}}\mathrm{\right)}\mathrm{d}\mathit{\nu ,}\end{array}$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: $\begin{array}{ccc}\mathit{J}& \mathrm{\equiv}& \int \begin{array}{c}\mathrm{\infty}\\ \mathrm{0}\end{array}{\mathit{J}}_{\mathit{\nu}}\mathrm{d}\mathit{\nu}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathrm{d}\mathit{\nu}{\mathrm{\int}}_{1}^{\mathrm{1}}\mathrm{d}\mathit{\mu}{\mathit{I}}_{\mathit{\nu ,\mu}}\mathrm{=}{\mathit{J}}_{\mathrm{v}}\mathrm{+}{\mathit{J}}_{\mathrm{th}}\mathit{,}\\ \mathit{H}& \mathrm{\equiv}& \int \begin{array}{c}\mathrm{\infty}\\ \mathrm{0}\end{array}{\mathit{H}}_{\mathit{\nu}}\mathrm{d}\mathit{\nu}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathrm{d}\mathit{\nu}{\mathrm{\int}}_{1}^{\mathrm{1}}\mathrm{d}\mathit{\mu}{\mathit{I}}_{\mathit{\nu ,\mu}}\mathit{\mu}\mathrm{=}{\mathit{H}}_{\mathrm{v}}\mathrm{+}{\mathit{H}}_{\mathrm{th}}\mathit{,}\\ \mathit{K}& \mathrm{\equiv}& {\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}{\mathit{K}}_{\mathit{\nu}}\mathrm{d}\mathit{\nu}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathrm{d}\mathit{\nu}{\mathrm{\int}}_{1}^{\mathrm{1}}\mathrm{d}\mathit{\mu}{\mathit{I}}_{\mathit{\nu ,\mu}}{\mathit{\mu}}^{\mathrm{2}}\mathrm{=}{\mathit{K}}_{\mathrm{v}}\mathrm{+}{\mathit{K}}_{\mathrm{th}}\mathit{,}\end{array}$where I_{ν,μ} is the specific intensity and θ the angle of a intensity with respect to the zaxis, μ = cosθ. The energy conservation of the total flux implies $\mathit{H}\mathrm{=}{\mathit{H}}_{\mathrm{v}}\mathrm{+}{\mathit{H}}_{\mathrm{th}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}\mathit{\sigma}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{,}$(A.20)where T_{irr} is the irradiation temperature given by ${\mathit{T}}_{\mathrm{irr}}\mathrm{=}{\mathit{T}}_{\mathit{\star}}\sqrt{\frac{{\mathit{R}}_{\mathit{\star}}}{\mathit{a}}}\mathit{,}$(A.21)where R_{⋆} is the radius of the host star and a the semimajor axis.
For the closure relations, we use the Eddington approximation (e.g. Chandrasekhar 1960), namely, $\begin{array}{ccc}{\mathit{K}}_{\mathrm{v}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}}{\mathit{J}}_{\mathrm{v}}\mathit{,}\\ {\mathit{K}}_{\mathrm{th}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}}{\mathit{J}}_{\mathrm{th}}\mathit{.}\end{array}$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): $\begin{array}{ccc}{\mathit{H}}_{\mathrm{v}}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}\mathit{\sigma}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\mathit{,}\\ {\mathit{H}}_{\mathrm{v}}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \mathrm{}\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}{\mathit{J}}_{\mathrm{v}}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{,}\\ {\mathit{H}}_{\mathrm{th}}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{J}}_{\mathrm{th}}\mathrm{(}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{.}\end{array}$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 T–P profile of the water vapor atmosphere. We assume that the boundary is at P_{0} = 1 × 10^{5} bar. The choice of P_{0} (≤1 × 10^{5}bar) has little effect on the atmospheric temperaturepressure structure. T_{0} is determined in an iterative fashion until abs(T_{0} − [πB(m = 0,P_{0},T_{0})/σ] ^{1/4}) ≤ 0.01 is fulfilled. Then we integrate Eqs. (A.9)–(A.13) over m by the 4thorder RungeKutta method, until we find the point where dlnT/dlnP ≥ ∇_{ad}. The pressure and temperature, P_{ad} and T_{ad}, are the boundary conditions for the convectiveinterior structure (see Sect. 2.1).
In Fig. A.1, we show the P–T profile for the solarcomposition atmosphere with g = 980 cm s^{2}, T_{int} = 300 K, and T_{irr} = 1500 K (dotted line). In this calculation, we take ${\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}$ and ${\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}$ as functions of P and T from Freedman et al. (2008) and calculate ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}$ and ${\mathit{\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 H_{2}, He, H_{2}O, CO, CH_{4}, Na, and K for the solar abundance respectively as $\begin{array}{ccc}{\mathit{\kappa}}_{\mathrm{v}}\mathrm{=}\{\begin{array}{c}\\ & & & & \\ & & & & \\ & & & & \\ & & & & \end{array}& & \end{array}$(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 P–T 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 P–T 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 P–T 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 P–T 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.
Fig. A.1 Temperature–pressure profiles for a solarcomposition 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}, T_{int} = 300 K, and T_{irr} = 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 ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}$ and ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}$ to obtain analytical solutions for J_{v} and H_{v}. Here we assume ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}$ and ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}$ are constant throughout the atmosphere. We differentiate (A.9) and (A.10) by m and obtain $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{J}}_{\mathrm{v}}}{\mathrm{d}{\mathit{m}}^{\mathrm{2}}}& \mathrm{=}& \frac{{\mathit{H}}_{\mathrm{v}}}{{\mathit{\mu}}^{\mathrm{2}}}\frac{\mathrm{d}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}}{\mathrm{d}\mathit{m}}\mathrm{+}\frac{{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}}{{\mathit{\mu}}^{\mathrm{2}}}{\mathit{J}}_{\mathrm{v}}\mathit{,}\\ \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{v}}}{\mathrm{d}{\mathit{m}}^{\mathrm{2}}}& \mathrm{=}& {\mathit{J}}_{\mathrm{v}}\frac{\mathrm{d}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}}{\mathrm{d}\mathit{m}}\mathrm{+}\frac{{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}}{{\mathit{\mu}}^{\mathrm{2}}}{\mathit{H}}_{\mathrm{v}}\mathit{,}\end{array}$where μ^{2} = K_{v}/J_{v}. Assuming J_{v} = H_{v} = 0 as m → ∞, we obtain $\begin{array}{ccc}\mathrm{(}{\mathit{J}}_{\mathrm{v}}\mathit{,}{\mathit{H}}_{\mathrm{v}}\mathrm{)}\mathrm{=}\mathrm{(}{\mathit{J}}_{\mathrm{v}\mathit{,}\mathrm{0}}\mathit{,}{\mathit{H}}_{\mathrm{v}\mathit{,}\mathrm{0}}\mathrm{)}\mathrm{exp}\left(\mathrm{}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\mathit{\mu}}\mathit{m}\right)\mathit{,}& & \end{array}$(A.30)where $\overline{{\mathit{\kappa}}_{\mathrm{v}}}\mathrm{=}\sqrt{{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{r}}}$ and J_{v,0} and H_{v,0} are the values of J_{v} and H_{v} 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 ${\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}\mathrm{(}{\mathit{J}}_{\mathrm{th}}\mathrm{}{\mathit{B}}^{\mathrm{)}}\mathrm{+}{\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}{\mathit{J}}_{\mathrm{v}}\mathrm{=}\mathit{Q}\mathit{.}$(A.31)We integrate Eq. (A.31) and obtain $\mathit{H}\mathrm{=}{\mathit{H}}_{\mathrm{\infty}}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{\left(}\mathit{m,}\mathrm{\infty}\mathrm{\right)}\mathit{,}$(A.32)where H_{∞} is the value of H evaluated at m → ∞ and $\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{\left(}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{m}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{{\mathit{m}}_{\mathrm{1}}}^{{\mathit{m}}_{\mathrm{2}}}\mathit{Q}\mathrm{\left(}{\mathit{m}}^{\mathrm{\prime}}\mathit{,\mu ,\phi}\mathrm{\right)}\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathit{.}$(A.33)To obtain H_{th} and J_{th}, we substitute Eq. (A.31) in Eqs. (A.11) and (A.12) and integrate by m. Then we obtain $\begin{array}{ccc}{\mathit{H}}_{\mathrm{th}}& \mathrm{=}& {\mathit{H}}_{\mathrm{\infty}}\mathrm{}{\mathit{H}}_{\mathrm{v}\mathit{,}\mathrm{0}}\mathrm{exp}\left(\mathrm{}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\mathit{\mu}}\mathit{m}\right)\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{\left(}\mathit{m,}\mathrm{\infty}\mathrm{\right)}\\ {\mathit{J}}_{\mathrm{th}}& \mathrm{=}& {\mathit{J}}_{\mathrm{th}\mathit{,}\mathrm{0}}\mathrm{}\frac{{\mathit{H}}_{\mathrm{v}\mathit{,}\mathrm{0}}}{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{\mathrm{\int}}_{\mathrm{0}}^{\mathit{m}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}\mathrm{exp}\left(\mathrm{}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\mathit{\mu}}{\mathit{m}}^{\mathrm{\prime}}\right)\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\\ & & \mathrm{+}\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{\mathrm{\int}}_{\mathrm{0}}^{\mathit{m}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}\mathrm{\{}{\mathit{H}}_{\mathrm{\infty}}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{(}{\mathit{m}}^{\mathrm{\prime}}\mathit{,}\mathrm{\infty}{\mathrm{)}}^{\mathrm{\}}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathit{,}\end{array}$where f_{Kth} = K_{th}/J_{th}, f_{Hth} = H_{th}/J_{th}, and ${\mathit{J}}_{\mathrm{th}\mathit{,}\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{H}\mathrm{th}}}\mathrm{\{}{\mathit{H}}_{\mathrm{\infty}}\mathrm{}{\mathit{H}}_{\mathrm{v}\mathit{,}\mathrm{0}}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{(}\mathrm{0}\mathit{,}\mathrm{\infty}{\mathrm{\left)}}^{\mathrm{\right\}}}\mathit{.}$(A.36)That is, we obtain $\begin{array}{ccc}\mathit{B}& \mathrm{=}& {\mathit{H}}_{\mathrm{\infty}}\left[\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{H}\mathrm{th}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{\mathit{\tau}}_{\mathrm{th}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right]\\ & & \mathrm{}{\mathit{H}}_{\mathrm{v}\mathit{,}\mathrm{0}}\left[\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{H}\mathrm{th}}}\mathrm{+}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\mathit{\mu}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{\mathit{\tau}}_{\mathrm{ext}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right]\mathrm{+}\mathit{E}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathit{,}\end{array}$(A.37)where $\begin{array}{ccc}{\mathit{\tau}}_{\mathrm{th}}\mathrm{\left(}\mathit{m}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{m}\\ \mathrm{0}\end{array}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathit{,}\\ {\mathit{\tau}}_{\mathrm{ext}}\mathrm{\left(}\mathit{m}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{m}\\ \mathrm{0}\end{array}\left(\overline{{\mathit{\kappa}}_{\mathrm{th}}}\mathrm{2}\mathrm{}\frac{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{{\mathit{\mu}}^{\mathrm{2}}}\overline{{\mathit{\kappa}}_{\mathrm{v}}}\mathrm{2}\right)\frac{\mathrm{1}}{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}}\mathrm{exp}\left(\mathrm{}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\mathit{\mu}}{\mathit{m}}^{\mathrm{\prime}}\right)\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathit{,}\\ \mathit{E}\mathrm{\left(}\mathit{m}\mathrm{\right)}& \mathrm{=}& \mathrm{}\left[\frac{\mathit{Q}}{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{f}}_{\mathit{K}\mathrm{th}}}{\mathrm{\int}}_{\mathrm{0}}^{\mathit{m}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{\left(}{\mathit{m}}^{\mathrm{\prime}}\mathit{,}\mathrm{\infty}\mathrm{\right)}\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathrm{+}\frac{\begin{array}{c}\mathrm{\u02dc}\\ \mathit{Q}\end{array}\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{\infty}\mathrm{\right)}}{{\mathit{f}}_{\mathit{H}\mathrm{th}}}\right]\mathit{,}\u2001\u2001\u2001\end{array}$and $\overline{{\mathit{\kappa}}_{\mathrm{th}}}\mathrm{=}\sqrt{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}}$. In our conditions, we assume $\mathit{\mu}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$, f_{Kth} = 1/3, f_{Hth} = 1/2 and Q = 0. Consequently, we obtain
the temperature profile as ${\mathit{T}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\left[\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}{\mathit{\tau}}_{\mathrm{th}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right]\mathrm{+}\frac{\sqrt{\mathrm{3}}}{\mathrm{4}}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}\left[\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\frac{\overline{{\mathit{\kappa}}_{\mathrm{v}}}}{\sqrt{\mathrm{3}}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}}\mathrm{+}{\mathit{\tau}}_{\mathrm{ext}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\right]$(A.41)where$\begin{array}{ccc}{\mathit{\tau}}_{\mathrm{ext}}\mathrm{\left(}\mathit{m}\mathrm{\right)}& \mathrm{=}& {\mathrm{\int}}_{\mathrm{0}}^{\mathit{m}}\frac{\overline{{\mathit{\kappa}}_{\mathrm{th}}}\mathrm{2}\mathrm{}\overline{{\mathit{\kappa}}_{\mathrm{v}}}\mathrm{2}}{{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}}\mathrm{exp}\left(\mathrm{}\sqrt{\mathrm{3}}\overline{{\mathit{\kappa}}_{\mathrm{v}}}{\mathit{m}}^{\mathrm{\prime}}\right)\mathrm{d}{\mathit{m}}^{\mathrm{\prime}}\mathit{.}\end{array}$(A.42)If we assume ${\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{p}}\mathrm{=}{\mathit{\kappa}}_{\mathrm{th}}^{\mathrm{r}}$ and ${\mathit{\kappa}}_{\mathrm{v}}^{\mathrm{p}}\mathrm{=}{\mathit{\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 (H_{2}OEOS) 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 GrantsinAid 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 GrantinAid for JSPS Fellows (No. 23003491) from MEXT, Japan.
References
 Adams, E. R., Seager, S., & ElkinsTanton, L. 2008, ApJ, 673, 1160 [NASA ADS] [CrossRef] [Google Scholar]
 Barclay, T., Rowe, J. F., Lissauer, J. J., et al. 2013, Nature, 494, 452 [NASA ADS] [CrossRef] [Google Scholar]
 Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative transfer (New York: Dover) [Google Scholar]
 Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 French, M., Mattsson, T. R., Nettelmann, N., & Redmer, R. 2009, Phys. Rev. B, 79, 054107 [NASA ADS] [CrossRef] [Google Scholar]
 Fressin, F., Torres, G., Désert, J.M., et al. 2011, ApJS, 197, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Grasset, O., Schneider, J., & Sotin, C. 2009, ApJ, 693, 722 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T. 2010, A&A, 520, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hansen, B. M. S. 2008, ApJS, 179, 484 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, MNRAS, 420, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
 Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Iro, N., Bézard, B., & Guillot, T. 2005, A&A, 436, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kane, S. R. 2007, MNRAS, 380, 1488 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution, Astron. Astrophys. Lib. (Berlin, Heidelberg, New York: SpringerVerlag) [Google Scholar]
 Lammer, H., Erkaev, N. V., Odert, P., et al. 2013, MNRAS, 430, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Lyon, S., & Johnson, J. D. 1992, LANL Tech. Rep. LAUR923407 (Los Alamos: LANL) [Google Scholar]
 Muirhead, P. S., Hamren, K., Schlawin, E., et al. 2012, ApJ, 750, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Nettelmann, N., Holst, B., Kietzmann, A., et al. 2008, ApJ, 683, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Nettelmann, N., Fortney, J. J., Kramm, U., & Redmer, R. 2011, ApJ, 733, 2 [Google Scholar]
 Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
 Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rothman, L. S., Gordon, I. E., Barbe, A., et al. 2009, J. Quant. Spectr. Rad. Transf., 110, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., Kuchner, M., HierMajumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Swift, D. C., Eggert, J. H., Hicks, D. G., et al. 2012, ApJ, 744, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Sasselov, D. D., & O’Connell, R. J. 2007, ApJ, 656, 545 [Google Scholar]
 Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, F. W., Sohl, F., Hussmann, H., Grott, M., & Rauer, H. 2011, Icarus, 214, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1986, Icarus, 67, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412 [NASA ADS] [CrossRef] [Google Scholar]
 Yelle, R., Lammer, H., & Ip, W.H. 2008, Space Sci. Rev., 139, 437 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Model of the planetary structure in this study. 

In the text 
Fig. 2 Concept of the chord optical depth. 

In the text 
Fig. 3 Mass evolution of closein waterrich 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 L_{p,0} = 1 × 10^{24} erg s^{1}, X_{wt,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 
Fig. 4 Threshold mass in M_{⊕} as a function of the initial planet’s luminosity in erg s^{1} for four choices of semimajor axes. The solid (red), dashed (green), dotted (blue), and dotdashed (purple) lines represent a = 0.02,0.03,0.05, and 0.1 AU, respectively. We have assumed X_{wt} = 75% and ε = 0.1. 

In the text 
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 X_{wt,0} = 100% (solid, red), 75% (dashed, green), 50% (dotted, blue), and 25% (dotdashed, purple). We have assumed L_{0} = 1 × 10^{24} erg s^{1}, a = 0.1 AU, and ε = 0.1. 

In the text 
Fig. 6 Relationship between the initial water mass fraction X_{wt,0} in % and the threshold mass M_{thrs} in M_{⊕} for four choices of semimajor axes of 0.02 AU (solid, red), 0.03 AU (dashed, green), 0.05 AU (dotted, blue), and 0.1 AU (dotdashed, purple). We have assumed L_{0} = 1 × 10^{24} erg s^{1} and ε = 0.1. 

In the text 
Fig. 7 Relationship between the initial planetary mass and the fraction of the initial water envelope that is lost via photoevaporation in 5 Gyr for six initial water mass fractions of X_{wt,0} = 1% (solid, red), 3% (longdashed, green), 10% (dotted, blue), 30% (dashdotted, purple), 50% (dotdashed, light blue), and 60% (dashed, black). We have assumed L_{0} = 1 × 10^{24} erg s^{1}, a = 0.1 AU, and ε = 0.1. 

In the text 
Fig. 8 Relationship between the threshold mass and the threshold radius. The latter is defined by the radius that the planet with M_{thrs} would have at 10 Gyr without ever experiencing mass loss ( denoted by R_{thrs}). The squares, which are connected with a solid line, are M_{thrs} and R_{thrs} for 0.1 AU and seven different initial water mass fractions X_{wt,0} = 100%, 75%, 50%, 25%, 10%, 5%, 1%, and 0.5%. The dashed and dotted lines represent massradius relationships, respectively, for rocky planets and purewater planets at 0.1 AU. ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ and ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ represent the minimum values of M_{thrs} and R_{thrs}, respectively. 

In the text 
Fig. 9 Relationship between the threshold mass M_{thrs} and radius R_{thrs} (lines; see text for definitions) compared to masses and radii of observed transiting superEarths around Gtype 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 M_{thrs} and R_{thrs} relationships for orbital periods of 11 days (=0.1 AU), 4 days (=0.05 AU), and 1 day (=0.02 AU), respectively. The dashdotted (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 L_{0} = 1 × 10^{24} erg s^{1}. “CoR” are short for CoRoT and “Kep” are short for Kepler. 

In the text 
Fig. 10 Upper panel: theoretical distribution of masses and semimajor 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 photoevaporation. The green line is the minimum threshold masses, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$. Here, we have adopted ε = 0.1. Lower panel: distribution of masses and semimajor axes (or incident fluxes) of detected exoplanets compared to the minimum threshold mass, ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}$, derived in this study (see Sect. 3.3 for definition). We have shown three ${\mathit{M}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{}\mathit{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 Ktype stars and Gtype stars). Planets are colored according to their zeroalbedo equilibrium temperatures in K. In the planet names, “CoR” and “Kep” stand for CoRoT and Kepler, respectively. 

In the text 
Fig. 11 Upper panel: theoretical distribution of radii and semimajor 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 photoevaporation 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, ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$. Here, we have adopted ε = 0.1. Lower panel: distribution of radii and semimajor axes (or incident fluxes) of Kepler planetary candidates, compared to the threshold radius, ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}$ (see Sect. 3.3 for definition). We have shown three ${\mathit{R}}_{\mathrm{thrs}}^{\mathrm{\ast}}\mathrm{}\mathit{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 (Gtype stars). 

In the text 
Fig. A.1 Temperature–pressure profiles for a solarcomposition 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}, T_{int} = 300 K, and T_{irr} = 1500 K. 

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