The evaporation valley in the Kepler planets

TL;DR

The paper explains the observed bimodal radius distribution of Kepler exoplanets through a minimal model showing photoevaporation causes a distinct valley, with implications for planet formation and evolution theories.

Contribution

It introduces a minimal model demonstrating how photoevaporation creates the observed radius valley in Kepler planets, linking it to envelope erosion timescales and planet formation conditions.

Findings

The radius valley results from envelope erosion timescales peaking at specific envelope masses.

Planets cluster around 3 Mearth in mass with initial H/He envelopes before gas disk dispersal.

Photoevaporation explains the bimodal distribution, but not the abundance of bare planets beyond 30-60 days.

Abstract

A new piece of evidence supporting the photoevaporation-driven evolution model for low-mass, close-in exoplanets was recently presented by the California-Kepler-Survey. The radius distribution of the Kepler planets is shown to be bimodal, with a ``valley' separating two peaks at 1.3 and 2.6 Rearth. Such an ``evaporation-valley' had been predicted by numerical models previously. Here, we develop a minimal model to demonstrate that this valley results from the following fact: the timescale for envelope erosion is the longest for those planets with hydrogen/helium-rich envelopes that, while only a few percent in weight, double its radius. The timescale falls for envelopes lighter than this because the planet's radius remains largely constant for tenuous envelopes. The timescale also drops for heavier envelopes because the planet swells up faster than the addition of envelope mass. Photoevaporation, therefore, herds planets into either bare cores ~1.3 Rearth, or those with double the core's radius (~2.6 Rearth). This process mostly occurs during the first 100 Myrs when the stars' high energy flux are high and nearly constant. The observed radius distribution further requires that the Kepler planets are clustered around 3 Mearth in mass, are born with H/He envelopes more than a few percent in mass, and that their cores are similar to the Earth in composition. Such envelopes must have been accreted before the dispersal of the gas disks, while the core composition indicates formation inside the ice-line. Lastly, the photoevaporation model fails to account for bare planets beyond ~30-60 days, if these planets are abundant, they may point to a significant second channel for planet formation, resembling the Solar-System terrestrial planets.

Equations

$$\frac{H}{R_c}$$

$$\rho ={\rho }_{\mathrm{rcb}}{\left[1+{\nabla }_{\mathrm{ab}}\left(\frac{G{M}_c}{c_s^2{R}_p}\right)\left(\frac{R_p}{r}-1\right)\right]}^{1/(\gamma -1)},$$

$$\rho \simeq {\rho }_{\mathrm{rcb}}{\left[{\nabla }_{\mathrm{ab}}\left(\frac{G{M}_c}{c_s^2{R}_p}\right)\left(\frac{R_p}{r}-1\right)\right]}^{1/(\gamma -1)}.$$

$$M_{\mathrm{env}}$$

$$I_2(R_c/R_p,\gamma )={\int }_{R_c/R_p}^1{x}^2{\left(x^{-1}-1\right)}^{1/(\gamma -1)}\mathrm{d}x.$$

$$I_2\approx {\nabla }_{\mathrm{ab}}{\left(\frac{\Delta R}{R_p}\right)}^{\gamma /(\gamma -1)}\approx {\nabla }_{\mathrm{ab}}{\left(\frac{\Delta R}{R_c}\right)}^{\gamma /(\gamma -1)}.$$

$$X$$

$$\frac{d\log T}{dr}={\nabla }_{\mathrm{ad}}\frac{d\log P}{dr}=-{\nabla }_{\mathrm{ad}}\frac{G{M}_c}{R_p^2{c}_s^2}.$$

$$\frac{d\log T}{dr}=-\frac{L}{4\pi R_p^2}\frac{3\kappa {\rho }_{\mathrm{rcb}}}{16\sigma T_{\mathrm{eq}}^4},$$

$$L=\frac{\mathrm{d}U}{\mathrm{d}t}\approx \frac{U}{{\tau }_{\mathrm{KH}}}\approx \frac{1}{{\tau }_{\mathrm{KH}}}\left({\int }_{R_c}^{R_p}\frac{G{M}_c\rho }{r}4\pi r^{2}dr\right).$$

$$L\approx \frac{1}{{\tau }_{\mathrm{KH}}}\frac{G{M}_c{M}_{\mathrm{env}}}{R_p}\frac{I_1(R_c/R_p)}{I_2(R_c/R_p)},$$

$$I_1(R_c/R_p)={\int }_{R_c/R_p}^{1}x(x^{-1}-1{)}^{1/(\gamma -1)}\mathrm{d}x,$$

$${\rho }_{\mathrm{rcb}}\approx \left(\frac{\mu }{k_b}\right){\left[\left(\frac{I_2}{I_1}\right)\frac{64\pi \sigma T_{\mathrm{eq}}^{3-\alpha -\beta }{R}_p{\tau }_{\mathrm{KH}}}{3{\kappa }_0{M}_{c}X}\right]}^{1/(1+\alpha )}$$

$$X$$

$${\rho }_c={\rho }_{M_{\oplus }}{\left(\frac{M_c}{1M_{\oplus }}\right)}^{1/4},$$

$$n_I$$

$$n_{\mu }$$

$$n_{\kappa }$$

$$n_T$$

$$n_{\rho }$$

$$n_M$$

$$n_a$$

$$n_b$$

$$X_2$$

$$t_{\dot{\mathrm{X}}}\equiv \frac{X}{\dot{X}}=\frac{M_{\mathrm{env}}}{{\dot{M}}_{\mathrm{env}}}.$$

$${\dot{M}}_{\mathrm{env}}=\eta \frac{\pi R_p^3{L}_{\mathrm{HE}}}{4\pi a^{2}G{M}_p}.$$

$$t_{\dot{X}}$$

$$\frac{\mathrm{d}X}{\mathrm{d}t}=-\frac{X}{t_{\dot{X}}}.$$

$$L_{\mathrm{HE}}=\{\begin{array}{ll}{L}_{\mathrm{sat}}& \mathrm{for}t<t_{\mathrm{sat}}\\ L_{\mathrm{sat}}{\left(\frac{t}{t_{\mathrm{sat}}}\right)}^{-1-a_0}& \mathrm{for}t\ge t_{\mathrm{sat}}.\end{array}$$

$$\frac{dN}{d\log P}\propto \{\begin{array}{ll}\mathrm{constant}& \text{for}P>7.6\text{days}\\ P^{1.9}& \text{for}P\le 7.6\text{days}.\end{array}$$

$$\frac{dN}{d{M}_c}\propto M_c{\exp }^{-M_c^2/2{\sigma }_M^2}.$$

$${\tau }_{\mathrm{KH}}=\{\begin{array}{ll}10^8\mathrm{yrs}& \mathrm{for}t<10^8\mathrm{yrs}\\ t& \mathrm{for}t\ge 10^8\mathrm{yrs}.\end{array}$$

$$t_{\dot{X}}(X=X_2,t=t_{\mathrm{sat}})\sim t_{\mathrm{sat}},$$

$$\frac{G{M}_p^2{X}_2}{8\pi R_c^3}\sim \eta t_{\mathrm{sat}}\frac{L_{\mathrm{HE}}}{a^2}\approx \eta {\mathcal{X}}_{\mathrm{HE}}$$

Full text

1. Introduction

Recent exoplanet discovery missions and targeted follow-up campaigns have fundamentally changed our understanding of what constitutes a “typical” planet (e.g. Borucki et al., 2011; Marcy et al., 2014). Specifically, the most common type of exoplanets are smaller than Neptune ($\lesssim$4 R⊕, e.g., Youdin 2011; Howard et al. 2012; Batalha et al. 2013; Petigura et al. 2013; Burke et al. 2014; Morton et al. 2016) and have masses of a few to tens of Earth masses. The fraction of sun-like stars that host at least one of these “Kepler” planets with an orbital period of less than 100 days is around 60–90% (e.g. Fressin et al., 2013; Silburt et al., 2015; Mulders et al., 2016).

Combining transit measurements of a planet’s radius with a measurement of its mass from transit-timing variations (TTVs, Carter et al. 2012; Wu & Lithwick 2013; Hadden & Lithwick 2014; Jontof-Hutter et al. 2016; Hadden & Lithwick 2016) or radial-velocity (RV) follow-up (e.g. Marcy et al., 2014; Weiss & Marcy, 2014) quickly told us that many of these planets had compositions unlike the small terrestrial planets in our own solar system (e.g. Wolfgang et al., 2016). Rather than being completely solid, such planets are likely to be composed of a dense solid core surrounded by a voluminous volatile rich envelope. However, due to degeneracies present in the mass-radius plane at low-masses (e.g. Adams et al., 2008; Rogers & Seager, 2010), we cannot infer the compositions just based on the current measured mass and radius for the majority of the observed planets. These degeneracies can be broken by considering how the local environment effects the evolution of the planet. For example, Wu & Lithwick (2013) showed that closer-in planets tend to be denser, while Ciardi et al. (2013) demonstrated that for pairs of planets in multi-planet systems, the inner planet tends to be smaller. These studies suggest that the envelopes are rich in hydrogen/helium.

Orbiting close to their parent stars, the Kepler planets can receive, over a lifetime, an integrated high-energy irradiation (high-energy exposure) that is an appreciable fraction of their gravitational binding energy (e.g. Lammer et al., 2003; Lecavelier Des Etangs, 2007; Davis & Wheatley, 2009), where this high-energy “exposure” is dominated by the first $\sim 100\,$Myr of the planet’s lifetime (e.g. Jackson et al., 2012). Planets with H/He-rich envelopes can be strongly evaporated by this irradiation (e.g. Yelle, 2004; Tian et al., 2005; Murray-Clay et al., 2009; Owen & Jackson, 2012; Johnstone et al., 2015; Erkaev et al., 2016), and indeed H/He evaporation has been observed from the low-mass planet GJ 436b (Kulow et al., 2014; Ehrenreich et al., 2015). Evaporation naturally results in planets that are smaller and denser than those at large separations (e.g. Lopez et al., 2012; Owen & Wu, 2013; Lopez & Fortney, 2013; Jin et al., 2014; Howe & Burrows, 2015). For a planet with a low enough mass and a close enough orbit, its initial low-mass H/He envelope can even be entirely stripped, leaving behind a naked solid core. This explains the dearth of planets with any envelopes at short periods (Lundkvist et al., 2016). The core’s mass and density play a primary role in controlling a planet’s evolution by setting the escape velocity (Owen & Wu, 2013; Lopez & Fortney, 2013; Owen & Morton, 2016; Zahnle & Catling, 2017), allowing one to break the compositional degeneracies in individual systems by statistical modelling (Owen & Morton, 2016). While steam atmospheres are also expected to lose mass (as in the case of early Venus – Kasting & Pollack 1983) the evaporative histories of such planets are significantly different to those that contain H/He envelopes and therefore distinguishable (Lopez, 2016).

While the evaporation theory naturally explains why Kepler planets are larger and less dense further out, it also makes another major prediction: the existence of an “evaporation valley”, a low-residence region in the radius-period plane between planets that have been completely stripped and those that are able to retain an envelope with roughly $\sim 1$% in mass. The evaporation valley was first predicted by Owen & Wu (2013) using numerical evolutionary studies for low-mass planets with pure rock (silicate) cores, and shortly after, by Lopez & Fortney (2013) for different core compositions using a different evaporation model. The evaporation valley was further reproduced by Jin et al. (2014) and Chen & Rogers (2016), again using different evaporation prescriptions and initial populations. The evaporation valley is thus a robust prediction of evaporative driven evolution of close-in H/He rich planets. This feature is largely independent of the assumed H/He evaporation model (energy-limited, recombination limited, UV driven, X-ray driven etc.), and its creation also appears insensitive to the choice of initial conditions (for at least a reasonable range of starting conditions); however, the details of these choices do control it’s properties (width, location with orbital period etc.).

The predicted occurrence valley, between stripped cores and those that retain a residual H/He envelope, is not particularly wide, with a radius width of $\sim 0.5\,$R⊕ (see Fig. 8 of Owen & Wu 2013 and Fig. 9 of Lopez & Fortney 2013). Large planetary radius errors, stemming from uncertainties in the stellar radius, have previously hampered efforts to observationally solidify its presence in the observed exoplanet population (see Owen & Wu, 2013, for a preliminary analysis that suggested its presence). Recently, the California-Kepler Survey (CKS), using spectroscopic follow-up of a large (1305) sample of planet hosting Kepler stars (Petigura et al., 2017) refined the planet parameters for 2025 Kepler planets, and reduced the typical planetary radius error to $\sim 10$% (Johnson et al., 2017). The CKS sample allowed Fulton et al. (2017) to definitely reveal a valley in the planet occurrence rate in the planet radius – period plane: close-in planets predominantly have a radius of either $\sim 1.3\,$R*⊕* or $2.6\,$R⊕, while planets with a radius of $1.8\,$R⊕ are considerably rarer, in spectacular agreement with the predicted evaporation driven evolution scenario for close-in exoplanets (Owen & Wu, 2013; Lopez & Fortney, 2013; Jin et al., 2014; Chen & Rogers, 2016).

While the presence of the evaporation valley in numerical models is robust to changes in model assumptions, a clear physical, first-principle description of its origin is missing, along with an understanding of how its properties change with model assumptions such as the core composition and evaporation model, and even when or whether it could be made to disappear. This paper serves two purposes: firstly, we clearly explain the physics behind the origin of the evaporation-valley; secondly, we will perform a preliminary investigation into how the observed evaporation valley can break many of the composition degeneracies and make inferences about both the composition of the Kepler planets and how/where they formed.

2. A minimal analytical model

Here, we build a minimal model of an evolving planet under the influence of evaporation and cooling. The planet is assumed to consist of a solid core of mass $M_{c}$, and radius $R_{c}$, surrounded by a gaseous envelope whose equation of state can be described by the ideal gas law. Such a model allows us to determine the planet outer radius as a function of time and envelope mass. We then perform numerical checks to confirm these results. While a fully numerical study is straight-forward (and has been done), such an approach obscurates the underlying properties of low-mass exoplanet atmospheres that results in the evaporation valley. Furthermore, the analytical method sheds light on many of the parameter dependencies.

2.1. Radius versus envelope mass for a low-mass planet

Our goal here is to understand how the radius of the planet ($R_{p}$) changes with its envelope mass fraction ($X=M_{\rm env}/M_{c}$) with a view to calculating the mass-loss timescale. We consider low-mass envelopes ($X<1$), so the planet’s mass is still dominated by the core. Therefore, we neglect the self-gravity of the planet’s envelope.111The role of self-gravity is to compress the envelope giving it a smaller radius than our following analysis predicts. Self-gravity becomes important when $X\sim 1$ and larger, and its impact is demonstrated in Fig 1.

Due to stellar insolation, the envelope is adiabatic (convective) in the deep interior and roughly isothermal (and radiative) near the surface. The roughly isothermal radiative cap satisfies $T\sim T_{\rm eq}$ (e.g. Rafikov, 2006; Lee & Chiang, 2015; Ginzburg et al., 2016), where the photospheric equilibrium temperature is set by stellar insolation. The radiative-convective boundary is assumed to occur at a density $\rho_{\rm rcb}$, and a radius $R_{\rm rcb}$. The planet’s radius is set by the photospheric radius and is typically $\sim 6$ pressure scale heights above $R_{\rm rcb}$ (see also Lopez & Fortney, 2014), where the scale height $H$ in the isothermal layer is,

[TABLE]

where we have estimated its value for a Sun-like host (mass 1 M⊙, radius 1 R⊙ and effective temperature 5780 K) and for an envelope with a mean-molecular weight of $\mu=2.35$ (solar composition), surrounding a core with an earth-like composition (see eq. 15). So the isothermal radiative cap is geometrically thin. In the following derivation, we take the planet radius to be approximately $R_{p}\approx R_{\rm rcb}$, but correct for the isothermal layer thickness in all relevant places (including all figures).

2.1.1 Convective interior

Adopting an equation of state for the adiabatic part of $\mathcal{P}=K\rho^{\gamma}$ with $\gamma$ and $K$ being constants,222In reality, for our ideal gas envelope $\gamma$ transitions from $\gamma=7/5$ in the upper envelope where molecules dominate to $\gamma=5/3$ in the lower envelope where molecules are dissociated. But the bulk of the planet can be considered as atomic. hydrostatic equilibrium gives a density profile of:

[TABLE]

where $\nabla_{\rm ab}\equiv(\gamma-1)/\gamma$ is the adiabatic gradient, and the isothermal sound speed $c_{s}^{2}\equiv\partial\mathcal{P}/\partial\rho|_{T}=\mathcal{P}/\rho$ is evaluated at $R_{\rm rcb}$. For planetary atmospheres that are strongly bound ($v_{\rm esc}^{2}\sim GM_{c}/R_{p}\gg c_{s}^{2}$)333Even if the initial planets may not be so, they rapidly evolve to such a state via the “boil-off” process investigated in Owen & Wu (2016)., the unity term inside the braket can be safely ignored for much of the planetary interior that contributes significantly to the envelope mass, and we can simplify the above expression into

[TABLE]

Solving for the mass enclosed in the atmosphere yields:

[TABLE]

where the dimensionless integral $I_{2}$ is:

[TABLE]

The properties of this dimensionless integral are discussed in the appendix. For $\gamma=5/3$, as applies for the bulk of the atomic interior, the integrand for $I_{2}$ peaks at $x=1/4$. So in the limit of a puffy envelope, $R_{c}/R_{p}\leq 1/4$, $I_{2}$ is fairly independent of $R_{c}/R_{p}$ and is of order unity. While in the opposite limit of a thin envelope, quantified as $\Delta R=R_{p}-R_{c}\ll R_{p}$, the integral can be approximated as

[TABLE]

With these expressions, the envelope mass fraction can be expressed as a function of the envelope thickness, $\Delta R/R_{c}$, in the two regions $\Delta R/R_{c}<1$ and $\Delta R/R_{c}>1$:

[TABLE]

where we have dropped all order unity constants for clarity.

While the above expression relates $X$ to $\Delta R/R_{c}$, it still contains a variable $\rho_{\rm rcb}$ that depends on planet radius and envelope mass, among other things. To eliminate these dependencies, we appeal to the characteristics of the radiative-convective boundary. At this location, the temperature gradient remains adiabatic by definition, or ${{d\log T}/{d\log P}}=\nabla_{\rm ad}$, giving rise to,

[TABLE]

The mode of energy transport changes at this point from advection by convective eddies to radiative diffusion, allowing us to write

[TABLE]

where we have approximated values for the temperature and radius at the radiative-convective boundary as those at the surface, the impact of this approximation is discussed in § 3.4. The internal luminosity arises from gravitational contraction of the atmosphere, i.e., changes in the gravitational binding energy $U$,444We ignore other internal heat sources, e.g., the heat capacity of the core and energy from radioactive decay. See discussion in § 4.3.

[TABLE]

Using eqs. (2) & (4), we can re-write the luminosity as

[TABLE]

where $\tau_{\rm KH}$ is the Kelvin-Helmholtz timescale (or cooling timescale). It is of order the planet age except in two cases: pre-cooling (the case considered in Owen & Wu, 2016) and high mass-loss rate (see §3.4). $I_{1}$ is another dimensionless integral given by

[TABLE]

whose relation to $I_{2}$ is discussed in the appendix. It is also shown there that the ratio $I_{1}/I_{2}$ smoothly varies from $1$ to $\sim 3$ over the parameter range of interest. Finally, we adopt an opacity law of $\kappa=\kappa_{0}{\mathcal{P}}^{\alpha}T^{\beta}$ to obtain the following expression for the density at the radiative convective boundary, by substituting eqs. (8) & (11) into eq. (9),

[TABLE]

Substituting this into eq. (7), we obtain the final monotonic dependence of envelope mass fraction on planet radius of

[TABLE]

where we have dropped all constants of physics and assumed that, due to compression, all solid cores (naked planets) have mass-radius relation as $M_{c}\propto R_{c}^{4}$ (Lopez & Fortney 2014, following, Fortney et al. 2007, Valencia et al. 2010), or

[TABLE]

where $\rho_{M\oplus}$ is the density of a 1 M⊕ core and depends only on the core composition. For terrestrial composition, $\rho_{M\oplus}=5.5~{}{\rm g\,cm^{-3}}$, while it is $11$, $4$ and $1.4$ g cm*-3* for pure iron, silicate and water/ice cores respectively (Fortney et al., 2007). The power indexes are, respectively,

[TABLE]

where we have also evaluated these expressions with $\gamma=5/3$, $\alpha=0.68$, $\beta=0.45$, the latter two identified by Rogers & Seager (2010) as the opacity law appropriate for a solar metalicity H/He envelope of a low-mass, highly irradiated planet (also see Freedman et al., 2008).

So, at a given planet size ($\Delta R/R_{c}$), the exponents in eq. (16) tell us the envelope mass is higher for a denser core composition (e.g., iron oxide vs. water ice), for an older planet, and for planets further away from their stars, all as expected. Furthermore, the nature of the opacity law makes the planet size a direct measure of the H/He envelope mass fraction, with weak sensitivities to all other factors. This has been noticed numerically in earlier works (e.g., Lopez & Fortney, 2014; Chen & Rogers, 2016). Lopez & Fortney (2014) provide power-law fits to their numerical models over a wide range of parameter space, and the power-law indices listed in eq. (16) are similar to their results (eq. 4 in that paper, though theirs are fits for an “enhanced opacity” model, not solar-metallicity). A different opacity law may lead to very different results.

Now, let us define an $X$ value at which the planet radius doubles, $\Delta R=R_{c}$, as $X_{2}$. We shall argue later that $X_{2}$ is a crucial parameter for producing the evaporation valley. In eq. (14), setting $\Delta R=R_{c}$ and re-inserting all constants that we have previously suppressed (with $\kappa=1.29\times 10^{-2}$ cm2 g*-1* at a pressure of 1 bar and a temperature of 1000 K, Rogers & Seager 2010), and for solar metallicity gas, we find:

[TABLE]

where $P$ is the orbital period. We have also set $T_{\rm eq}=(L_{*}/16\pi a^{2})^{1/4}$ and adopted an empirical mass-luminosity relation for low-mass dwarfs: $L_{*}/L_{\odot}=(M_{*}/M_{\odot})^{3.2}$ (see, e.g. Cox, 2000). In principle, the value of $I_{1}/I_{2}$ also depends on $\Delta R/R_{c}$. But since it varies by only a factor of $3$ within our range of interest, together with the weak index $n_{I}$, this can be safely ignored. Moreover, when the thickness of the isothermal layer is accounted for, the true $X_{2}$ is slightly smaller. We actually obtain that $X_{2}$ is of order a few percent at an age of a Gyrs.

In summary, the envelope mass fraction, at which the planet’s radius doubles, is of order a few percent, and is a relatively weak function of planet period, stellar mass, core mass and planet age.

2.2. Timescale for atmospheric erosion

We define a mass-loss timescale for envelope evaporation as

[TABLE]

In the following, we show that this timescale naturally peaks for envelopes with masses $X\approx X_{2}$, or those which double the core’s radius.

The mass-loss rate is given by the ratio of photo-evaporative power and the binding energy of the planet (Lecavelier Des Etangs, 2007; Erkaev et al., 2007). Let $L_{\rm HE}$ be the luminosity of high energy photons from the star, and the dimensionless factor $\eta$ be the efficiency of these photons for mass-removal,555This efficiency is defined as if the high-energy photons are intersected by a cross-section $\pi R_{p}^{2}$. The definition of this efficiency is arbitrary and many other definitions do exist in the literature. the mass-loss rate is

[TABLE]

The efficiency factor is not necessarily a constant, as has been demonstrated in multiple works (e.g. Murray-Clay et al., 2009; Owen & Jackson, 2012; Shematovich et al., 2014; Owen & Alvarez, 2016; Salz et al., 2016). However, when compared to full radiation-hydrodynamic models it takes a value of order $0.1$ for low-mass planets (Owen & Jackson, 2012; Owen & Alvarez, 2016). For simplicity, we adopt a constant $\eta=0.1$ (the so-called “energy-limited” approach, see e.g., Lopez & Fortney, 2013). However, we discuss its variation and impact on the planet population in §4.5. We do not consider an “effective absorption radius” or “expansion radius” (e.g. Baraffe et al., 2004) that accounts for the higher radius at which the high-energy flux is absorbed compared to the planet’s radius, as this effect can always be folded into the efficiency factor. Such a radius is also difficult to define for X-ray driven evaporation (which is important for low-mass planets) as different wavelength photons can be absorbed at very different radii. We also neglect the effect of stellar tides (e.g. Erkaev et al., 2007) since most observed planets are far from their Roche radii. Specifically $<13\%$ of the planets (taking them to have a mass of 3 M⊕) used in the Fulton et al. (2017) sample would have their mass-loss rates increased by $>20\%$ and $<2\%$ would have their mass-loss rates increased by a factor of $>2$ using the (Erkaev et al., 2007) prescription.

Substituting eq. (14) into the above expressions, we obtain:

[TABLE]

where we evaluate for solar-metallicity gas, and have adopted a parameter $f$ to account for the difference between $R_{\rm rcb}$ and the photospheric radius, $R_{p}=fR_{\rm rcb}$. In our model, $f$ is computed self-consistently, by locating the photosphere at where ${\mathcal{P}}=(2/3)g/\kappa$. The photosphere radius is then at $nH$ above the radiative-convective boundary, with $H$ being the local scale height, and $n=\ln(\rho_{\rm rcb}/\rho_{\rm ph})$, where $\rho_{\rm ph}$ is the density at the photosphere.

According to eq. (20), the mass-loss timescale reaches a maximum at the point at which the radius is doubled ($\Delta R\approx R_{c}$, $X=X_{2}$). Such a behaviour arises because, below $X_{2}$, the planet’s radius is dominated by that of the core and is independent of the envelope mass. So the photoevaporation timescale decreases for lower envelopes. This is an unstable situation that can cause complete stripping. While above $X_{2}$, increasing $X$ leads to a radius expansion. The expansion is so significant that, for our adopted opacity laws and a constant $\eta$, the evaporation timescale shortens. This behaviour can be shown to occur for all well known evaporation models: energy-limited evaporation, when the radius increases faster than $X^{1/3}$; photon-limited evaporation (Owen & Alvarez, 2016), where the radius increases faster than $X^{1/2}$; and recombination limited evaporation (Murray-Clay et al., 2009), where the mass-loss rate is exponentially sensitive to radius (Owen & Alvarez, 2016) such that the erosion time-scale always decreases with increasing envelope mass fraction. The fact that the envelope-mass loss time scale peaks for low-mass planets with envelope mass fractions of order a few percent was also noticed numerically by Chen & Rogers (2016) using mesa models. With our minimal model, it becomes obvious that the mass-loss timescale always peaks around $X=X_{2}$. For this feature not to happen, one requires $n_{b}<3$, or in terms of the opacity scaling with pressure, $-7/3<\alpha<-2$, which is unrealistic.

The heuristically derived eq. (20) has a sharp discontinuity at $\Delta R/R_{c}=1$. Numerically, we smooth the transition in the preceding calculations by using the fact, in reality $R_{p}=R_{c}+\Delta R$, rather than the approximation $R_{p}=R_{c}$ for $\Delta R<R_{c}$ and $R_{p}=\Delta R$ for $\Delta R>R_{c}$, presented above.

In the following, we will construct planetary structure models, based on eqs. (4), (11) and (13). In Fig. 1, we compare results from our minimal models against those produced using the mesa code (Paxton et al., 2011, 2013, 2015), suitably modified for highly irradiated low-mass planets (see Owen & Wu, 2013; Owen & Morton, 2016). The agreement between the analytic model and the mesa results is good. The mass-loss timescale peaks at $X\approx 0.02-0.03$, roughly where the planet doubles in size, increases slowly with cooling time; and the mass-loss timescale falls off as predicted on both sides of the peak. The agreement only fails for cases where $X\gtrsim 1$, where the self-gravity of the envelope (ignored in our analytical models) becomes dominant. At large envelope masses, self-gravity keeps the planet’s radius roughly constant to of order a Jupiter radius, and the mass-loss time begins to increase again as expected (it is well known hot jupiters are stable to evaporation, e.g. Yelle 2004).

3. Comparing to the observations – the “evaporation valley”

Based on refined stellar parameters from the CKS survey, Fulton et al. (2017) reported a “gap” in the planetary radius distribution (also noticed by Owen & Wu, 2013, based on the cruder KIC data): a deficit of planets at radii $1.8\pm 0.2R_{\oplus}$ (“gap”) that appears to extend from 10 days to 100 days in orbital period (“valley”). We now use our analytical model to calculate the exact valley location for different parameters. In the following, we detail our model choices and present results.

We make no effort to “fit” the observations at this stage, especially with our minimal model. That is best left to comparisons with numerical models. Nevertheless, as we will see, even our model can be used to make strong inferences about the properties and possible formation channels of the Kepler planets. Usefully, we can analytically understand the origin of these constraints.

3.1. The Integration

Starting from an initial envelope mass fraction at time zero, which we set to be 1 Myr (roughly the disk dispersal timescale, e.g. Hernández et al., 2007), we evolve the envelope mass fraction according to:

[TABLE]

according to eq. (20), and a prescription for the stellar high energy luminosity (§3.2), and the cooling of the planet (§3.3). Some example evolutionary tracks for a planet with different initial envelope mass fractions is shown in Fig. 2.

We have also applied our analytical model to the envelope evolution of the Kepler-36b/c system and compare them against the more detailed mesa calculations presented in Owen & Morton (2016). We find that while we can reproduce their general results, namely, the lower mass Kepler-36b is completely stripped off, while the higher mass Kepler-36c retains a bulky envelope, our models suffer moderately more mass loss for the latter planet during the early stages when its mass-loss rate is high, due to our assumption of constant mass-loss efficiency, an assumption not adopted by Owen & Morton (2016).

3.2. Stellar Parameters

We adopt host star masses similar to those in the CKS sample (Fulton et al., 2017), a Gaussian distribution in mass, centred at $1.3M_{\odot}$ with a variance of $0.3M_{\odot}$.

For the magnitude and the evolution of the high-energy flux ($L_{\rm HE}$, including UV through X-ray radiation), we adopt the empirical relation for main-sequence dwarfs, as summarized by Jackson et al. (2012)

[TABLE]

We further choose $a_{0}=0.5$ and $t_{\rm sat}=100$ Myrs, and $L_{\rm sat}\approx 10^{-3.5}L_{\odot}(M_{*}/M_{\odot})$, motivated by the body of observational and modelling works (e.g. Güdel et al., 1997; Ribas et al., 2005; Jackson et al., 2012; Tu et al., 2015). Since $a_{0}>0$, the time-integrated high-energy “exposure” (Lecavelier Des Etangs, 2007; Jackson et al., 2012; Owen & Wu, 2013) is dominated by that in the first $100$ Myrs. So the exact choice of $a_{0}$ has little bearing on the final planet properties.

3.3. The Planet Ensemble

We now must decide what the primordial Kepler planets look like. In this work, we consider only one population of planets – we discuss the evidence for a second population in §4.7.

We adopt the following orbital period distribution for planets around all stars,

[TABLE]

This distribution is obtained by fitting the Kepler planet sample and correcting for transit probabilities. It is similar to those obtained in earlier work (e.g. Fressin et al., 2013).

For the planetary cores, we assume them to be Earth-like ($\rho_{M_{\oplus}}=5.5\,{\rm g~{}cm^{-3}}$). In §4.2, we will vary this parameter and show that the observed CKS sample actually demands this choice. Furthermore, we take a Rayleigh distribution for the core mass

[TABLE]

The actual mass distribution of Kepler planets has not been reliably established. Early RV studies (e.g. Howard et al., 2010; Mayor et al., 2011) showed that the planet occurrence rate increased towards low-masses: planets in the mass range 3-10 M⊕ being most common, with the occurrence rate falling rapidly towards higher masses. One more recent attempt is provided by Marcy et al. (2014) where they selected 22 KOIs and measured masses (or their upper-limits) for 42 planets. The mass function for this group can roughly be fit as a Rayleigh distribution with a mode $\sigma_{M}\sim 5\pm 1\,$M⊕. However, this sample is arguably biased towards higher masses – compared to Kepler planets within the same period range (say, $<20$ days), this sample contains planets that are, on average, larger and therefore likely more massive. We therefore choose a slightly smaller mass scale of $\sigma_{M}=3\,$M⊕. We are also motivated by the position of the small-size peak. At a radius of $1.3\,$R⊕, this corresponds to a terrestrial planet with $M_{c}=3\,$M⊕.

We also need to prescribe the initial envelopes. This includes both their initial thermal time (which determines the initial entropy) and initial masses. Owen & Wu (2016) argued that after low-mass planets are born, they undergo a rapid phase of cooling and mass-loss (the “boil-off” phase) and age prematurely to $\tau_{\rm KH}\sim 100$ Myrs within a relatively short time. After about 100Myrs, they continue to cool off normally. So we take

[TABLE]

Due to the fact $a_{0}>0$ and the above form of cooling contraction, most of the envelope erosion occurs in the first 100 Myrs, and there is little change to the planet atmosphere after that.

To guide our choice for the initial envelope mass, we first ask what kind of planets would have occupied the observed gap, at an age of a few billion years. These are the grey points in Fig. 3, for cores of terrestrial composition ($\rho_{M\oplus}=5.5~{}{\rm g~{}cm^{-3}}$), with a range of core mass, initial envelope mass, and orbital period. Since these models would appear in the CKS gap, they are disfavoured by the observations. Now we make the important assumption that planets are born with the same envelope mass and core mass distributions across the range of period of interest (from a few to $100$ days). This assumption mostly reflects our ignorance on the proper initial condition, and is almost certainly not true in reality. But under this restrictive assumption, all models that are marked with grey in the right-hand panel of Fig. 3 are not permitted. In other words, there are two groups of planets that are likely progenitors for the observed Kepler planets: one group are planets that are slightly more massive ($M_{c}>3$ M⊕) and are born with at least a few percent of H/He; the other group, planets that are less massive ($M_{c}<3$ M⊕) and that are born, for all purposes, bare. The second group cannot account for planets with radii larger than $1.8$ R⊕. Furthermore, Fulton et al. (2017) indicate that the completeness of their exoplanet radius distribution becomes uncertain at small planet radii $<1.14$ R⊕. Inspection of Fig 3 indicates that planets with masses $<2$ M⊕ ($R_{p}=1.18$ R⊕) are disfavoured if they contain H/He envelopes with masses of a few percent or more. Therefore, while the current CKS sample cannot determine if planet occurrence drops significantly below 1.14 R⊕ we predict it should, unless there is a second population of planets (see §4.7). Moreover, we emphasize that the choice of the mass scale, $\sigma_{M}=3\,$M⊕, is motivated by both the radii of the bare planets, and the size distribution of planets to the right of the valley.

Guided by this insight, we adopt a logarithmically flat distribution for $X_{0}$, the initial envelope fraction, with $X_{0}\in[X_{\rm min},X_{\rm max}]$ and $X_{\rm min}=0.01$ and $X_{\rm max}=0.3$. The value of $X_{\rm max}$ matters little as long as it lies well above $1\%$, while the value of $X_{\rm min}$ is suggested by Fig. 3, for our choice of core masses.

3.4. A few details

The radius measured for a planet by the transit method is not the photospheric radius (which accounts for the radial light-path for the local black-body photons), but should be determined by the tangent light-path for transmitting stellar photons. The chord is longer for the latter group by a factor of $\sim\sqrt{8R/H}$, increasing the transit radius by $\sim\ln(\sqrt{8R/H})H\sim 3H$ (e.g. Lopez & Fortney, 2014). However, opacity for optical photons is smaller than that for infra-red photons, leading to a “deeper photosphere”. These two effects cancel each other to some degree and we decide to ignore them here.

Another complication arises from how we treat the thermal evolution of the planet in presence of mass-loss. Here, we simply assume the two are independent (eq. 25). This is inappropriate when the envelope is being evaporated faster than cooling contraction. In this case, the lifting of pressure at the top allows the envelope to expand adiabatically. The irradiated atmosphere then actually transports heat inward, maintaining the same internal entropy. The radiative-convective boundary remains fixed at the same density, and the internal entropy can be transported out with the same cooling luminosity as before. This accelerates the cooling compared to the case of no-cooling (this is the physical basis for the “boil-off” discussed in Owen & Wu 2016, see also the discussion in Ginzburg et al. 2016). We do not correct for this effect in our minimal model, but note that this may lead to enhanced mass loss in some cases; however, models that typically enter this region are on their way to be completely stripped anyway.

3.5. Planets at a few Gyrs

We evolve our initial population for three billion years under the influence of cooling and evaporation. The resulting radius distribution for planets with orbital periods less than 100 days is presented in Fig. 4. Compared to the observed 1-D radius distribution from the CKS sample, our model reproduces the observed features: the positions of the radius gap and the radius peaks, the widths of the peaks. More significantly, the observed valley in the radius-period 2-D plane is reproduced nicely by our model. This is shown in Fig. 5, alongside with the observed results by Fulton et al. (2017). Both these successes strongly support the photoevaporation theory for the evolution of Kepler planets.

Since we have not attempted to “fit” the observed distribution, but rather to match its generic features, the agreements are not perfect. Short of theoretical predictions for the initial planet properties (how they depend on, e.g., separation, stellar mass, planet core mass, etc.), and only using an approximate photoevaporation model (the energy-limited case), this not surprising. In future work, using numerical models, one may actually be able to use the observed planet population to infer these initial properties. In the following, we discuss what parameters change may impact on the agreements.

4. Discussion

We have used an analytic model to demonstrate that the evaporation valley is a robust outcome of the evolution of close-in, low mass planets with volatile envelopes. The mass-loss timescale always peaks for a H/He envelope of a few percent in mass, where the atmosphere roughly doubles the planet’s radius (see also Chen & Rogers, 2016). This simple combination of planetary structure and evaporation is the origin of the bimodal radius distribution reported in Fulton et al. (2017). Here, we explain these results in more depth, and investigate impacts on our model from various model uncertainties.

4.1. The origin of the evaporation valley

The origin for the evaporation valley is schematically shown in Fig. 6. Consider a group of identical, low-mass planets differing only in their initial envelope mass fractions ($X_{0}$) and high-energy exposure. Since after the first 100 Myrs, both the stellar flux decays and the planet cools down to a smaller size, the evaporation is dominated by that in the early stages. For a group of planets with the same high-energy exposure, if the peak mass-loss timescale – $t({\dot{X}})$ for $X=X_{2}$ – is shorter than 100 Myrs (model d in Fig. 6), all envelopes are stripped bare and we expect to see only naked cores, as is the case for the observed Kepler planets at small separations (e.g. Dressing et al., 2015). If, on the other hand, the peak mass-loss timescale is longer than $100$ Myrs (group c in Fig. 6), there is a bifurcation of final envelope masses – planets with initial envelope masses $X<X_{2}$ will be completely stripped and present as naked cores; while planets with initial $X>X_{2}$ will be herded towards $X=X_{2}$. This manifests as a bifurcation in planet radius, with peaks at both the core size and its double. This bifurcation is the origin of the evaporation valley.

Lastly, for planets in group a and b which experience too little evaporation, there are little modifications to their atmospheres unless they start with extreme envelope masses (very high or very low). There is still a shepherding toward the above mentioned two peaks, but the widths of the second peak (at twice the core radius) are broader. The observed valley shape can then be used to exclude some of these models, as is done in Fig. 3, where we show that most of the initial envelopes should have more than a few percent in mass.

This understanding allows us to derive some useful scalings. These are detailed here and below.

Since less massive planets can be stripped out to larger distances, the bare population should have a decreasing size (smaller mass and hence smaller core radius) going away from the star. We can derive this by requiring that, for a given core mass, the longest evaporation time (when $X=X_{2}$, or $\Delta R\sim R_{c}$) is of order $t_{\rm sat}$ or 100 Myr,

[TABLE]

or,

[TABLE]

where $\mathcal{X}_{\rm HE}$ is the high-energy “exposure” for a given planet. Then combining eq. (15) & (17), we can find the radius of the most massive planet that can be stripped ($R^{\rm bot}_{\rm valley}$) at a given exposure as:

[TABLE]

Using the high-energy flux dependence on stellar mass (§3.2), ${\mathcal{X}}_{\rm HE}\propto M_{*}/a^{2}$, as well as eq. (14), we convert the above relation to find $R^{\rm bot}_{\rm valley}\propto P^{-0.25}$, shown as the dotted black line in the left hand panel of Fig. 5. This explains the topology of the evaporation valley in the radius-period plane, in models with a constant evaporation efficiency $\eta$. As we discuss later, this scaling can change when the evaporation model is different, while the core composition can shift the black line vertically.

4.2. Core Composition

When the evaporation valley was first predicted by Owen & Wu (2013); Lopez & Fortney (2013) and subsequently by Jin et al. (2014), its location was conjectured to be a discriminant of core composition which could provide clues as to their formation (Lopez & Rice, 2016). And indeed our model shows clear dependency on the core composition.

Following our discussion above, we posit that the valley, independent of the core composition, always lies in-between $R_{c}$ and $2R_{c}$. Let us say, $\sqrt{2}R_{c}$. Or

[TABLE]

This is demonstrated in Fig. 7, where cores of the same mass but different composition exhibit different valleys. An order of magnitude in change in $\rho_{M_{\oplus}}$ leads to an easily detectable shift in the valley position of a factor 2, much larger than the $10\%$ radius errors achieved in the CKS sample (Johnson et al., 2017), and even large compared to the $30\%$ error in the general Kepler sample. So if the planet mean masses are known, the valley position yields the core composition.

Interestingly, even when the planet mean masses are not known, if we assume the bimodal radius distribution arises from a single population of planets, we can break the degeneracy and determine the composition. In our model, we successfully reproduce the observed valley using a single population of planets with $M_{c}\sim 3~{}{\rm M}_{\oplus}$ (see also §4.7) and terrestrial compositions. This does not seem possible for some other compositions. For instance, if the cores are made up largely of ice/water with $\rho_{M\oplus}=1.3\,{\rm g~{}cm^{-3}}$, the bare cores at $R=1.3~{}{\rm R}_{\oplus}$ will correspond to $M_{c}\sim 0.5~{}{\rm M}_{\oplus}$, and none of these cores can retain enough hydrogen to occupy the second radius peak. Similarly, for pure iron composition ($\rho_{M\oplus}=11\,{\rm g~{}cm^{-3}}$), the bare cores should correspond to $M_{c}\sim 6~{}{\rm M}_{\oplus}$. However, few of these planets could be evaporated down to naked cores, at the distances that we observe them. So the current data exclude cores that are mostly icy, and favour compositions that are terrestrial-like, i.e., silicate-iron composite. This is consistent with RV results from known bare planets (e.g. Dressing et al., 2015).

While these cores are Earth-like, we find that we cannot constrain their iron fractions to a narrow range using current data: in Fig. 8, we present the expected size distributions for different iron fractions, and all appear largely consistent with data. This could reflect either a genuine composition spread in real planets, or an intrinsically narrow composition spread that is smeared out by errors in the observed data.

In conclusion, the observed gap suggests that the most common core composition is Earth-like. But the iron fractions in these silicate-iron composite could span a wide range; with a detailed numerical fit to the observed population one maybe able to constrain the Iron fraction to a narrower range.

4.3. Core Luminosity

For simplicity, in the main analysis, we have ignored the luminosity contribution from the cooling core. This is not justified in general, since the thermal content in the more massive core can be more substantial than that in the envelope. We address this deficiency here.

We assume that the entire planet, core and atmosphere, cools with a single temperature. This is reasonably accurate for the solid interior, since even the adiabatic temperature gradient is relatively shallow. For an envelope with $\Delta R\lesssim 4R_{c}$, ($X\lesssim 0.1$) the thermal inertia is dominated by gas in the bottom few scale heights close to the core. Gas here likely share the same temperature as the solid part. We further assume that the core contributes only through its primordial heat, and ignore any radiogenic source. This is valid, especially at early times. For the heat content of the core, we adopt a heat capacity of $10^{7}$ erg g*-1* K*-1* (Valencia et al., 2010, similar to rock), and we assume that the core is made up of SiO2 molecules with a mean molecular weight of $76$ (or per particle specific heat of $\sim 7k_{B}$). The envelope is assumed to be composed of H-atoms with a per particle specific heat of $3/2k_{B}$. So including the core contribution to the internal luminosity, we have

[TABLE]

This expression explains the results shown in Fig. 3 of Lopez & Fortney (2014). We insert this factor into eq. (11) and re-calculate the photo-evaporation process. Core luminosity allows the adiabatic envelope to remain large for a longer time, enhancing mass-loss. However, the overall effect is minor. This is because the evaporation bottleneck is when the envelope mass fraction is of order a few percent ($X=X_{2}$). At this point, the core contribution to the luminosity is an order unity effect.

4.4. Atmosphere Metallicity

Some studies of transmission spectroscopy have suggested that the atmospheres of hot Neptunes (e.g., GJ1214b) are highly enriched in metals, with [Z/H] perhaps as high as 100 (e.g. Charnay et al., 2015), while others (e.g. HAT-P-26b) are consistent with solar metallicity (Wakeford et al., 2017) . Here we explore how atmosphere metallicity may affect our results.

Metallicity enters our model in two ways (eq. 14). First, higher metallicity increases the mean molecular weight ($\mu$), which reduces the scale height at a given temperature. This shrinks the atmosphere, and increases the value of $X_{2}$. Second, higher metallicity increases opacity ($\kappa_{0}$), which shifts the convective-radiative boundary. Thus, the atmosphere looks more inflated for a given age. The outcome of our model is determined by the competition between these two effects (also see Howe & Burrows, 2015). Reality may be more complicated. The constant energy-efficiency approach we adopt here may be severely invalid if metals affect the driving of the photo-evaporative flow.666For X-ray driven flow, metals (like C, O) dominate both the absorption of X-ray photons and the line cooling of the heated region (Owen & Jackson, 2012). The two effects may cancel out to a large degree. For EUV driven flow, hydrogen photo-ionization and recombination dominate the heating and cooling, though metals, if present at a large quantity, may reduce the energy efficiency as additional sources of radiative cooling. Metals may have additional thermal effects by, for instance, producing high altitude hazes/clouds.

Assuming that $\kappa_{0}\propto Z$ (e.g. Lee & Chiang, 2015), we find that, overall, raising the metallicity leads to a slightly larger radius and hence a slightly shorter mass-loss timescale. The effect maximizes (making a factor of 2 difference in timescale) when the metallicity is $10\times$ solar, but nearly vanishes when it is $100\times$ solar. Such an effect is so small it is eclipsed by other uncertainties in our current model. As a result, we can not make conclusive statement regarding atmosphere metallicity, aside from the fact that the initial atmosphere must be dominated by H/He in number (and so in molecular weight). For example, the observed position in the radius-period plane rules out the presence of “water-worlds” whose volatile envelopes are primarily water/steam. Water/steam atmospheres have a much larger mean-molecular weight ($\mu\sim 18$) compared to H/He envelopes, this significantly increases the value of $X_{2}$ to values more like $\sim 0.5$ (Lopez, 2016). The increased high-energy exposure required to evaporate an envelope with $X=0.5$ (even with an optimistic efficiency of 0.1 for steam atmospheres) results in the evaporation valley appearing at much shorter separations, around a period of 2 days (Lopez, 2016), rather than the 10 days it is observed. Therefore, water-worlds making up any reasonable fraction of the Kepler planets is clearly ruled out by the CKS observations.

4.5. Evaporation Efficiency

There are two reasons why the evaporation efficiency may not be constant for all planets. First, geometry. We define the efficiency $\eta$ to be based on the light received by the planet disk ($R_{p}$). However, the true cross-section is determined the UV/X-ray photosphere. This lies well above the planet photosphere, especially so for planets with large surface scale heights. This raises the efficiency for those planets. Second, physics. Most hydrodynamic models predict that the efficiency should drop as one moves to more massive and/or denser planets. The deeper gravitational potential in this case means it takes the flow a longer time to escape, allowing it to lose more energy by radiative cooling (e.g. Owen & Jackson, 2012). The “energy-limited” evaporation is only applicable to UV evaporation of weakly irradiated planets (where recombination-equilibrium cannot be reached, and radiative cooling is inefficient on flow-timescales) with low escape velocities (Owen & Alvarez, 2016), but not the full range of planets observed to make up the evaporation valley, where recombinations (e.g. Murray-Clay et al., 2009) and X-ray evaporation (e.g. Owen & Jackson, 2012) complicate the picture.

We use a simple model to illustrate these effects. We assume that the efficiency of mass-loss scales as777We caution the below expression is only suitable in a narrow range of parameter space. We use it here for illustrative purposes and it should not be used to describe the results from radiation-hydrodynamic evaporation models.

[TABLE]

where $v_{\rm esc}$ is measured at the planet photosphere and the normalization of $15{\rm~{}km~{}s^{-1}}$ is chosen to reproduce the observed distribution around 10 days. The power $2$ roughly coincides with the result shown in Fig. 13 of Owen & Jackson (2012), where $\eta\propto M_{c}^{-3/4}$ for earth-like cores in the range 1-10 M⊕. We adopt this efficiency scaling over other published scalings (e.g. Salz et al., 2016) as it is more appropriate for young, active stars which drive most of the mass-loss, rather than older stars where evaporation does not affect a planet’s evolution.

With such an ad-hoc efficiency law, low-mass planets can be stripped out to larger distances, while high-mass planets are more resistent to stripping. Re-working through the derivation leading to eq. (28), we obtain $R^{\rm bot}_{\rm valley}\propto P^{-0.16}$, a shallower dependence than that in the constant efficiency model. This is illustrated in Fig. 9, where one observes an extended tail of small-planets towards long periods.

To further illustrate this effect, in Fig. 10, we plot the values of $R^{\rm bot}_{\rm valley}$ for Earth-like cores and those with a composition of 1/3 Ice, 2/3 silicates, using the Fortney et al. (2007) mass-radius relationship. This figure demonstrates that, with a large enough sample of planets, and a $\sim 10$% radius errors, one could in principle distinguish between different evaporation models and different core compositions. The bottom panel translates the radius into planet core mass, demonstrating how the mass of the most massive stripped core depends on orbital period, for any of the models. Follow-up mass measurements can test these predictions.

Finally, we have not considered non-thermal escape processes such as stellar wind stripping. This is because hydrodynamic evaporation dominates for the separations we are interested in here (e.g. Murray-Clay et al., 2009). We have also not included the effects of a planetary magnetic field. Presence of a magnetic field strong enough to retain its dipolar structure (despite the evaporation) close to the planet may suppress the mass-loss rates significantly (e.g. Adams, 2011; Owen & Adams, 2014; Khodachenko et al., 2015; Owen & Adams, 2016). Unfortunately, it is unknown if the Kepler planets possess sufficiently strong fields.

4.6. Stellar Mass

The deficit of planets between $R_{c}$ and $2R_{c}$ results from interaction between planet internal structure and photo-evaporation. It is not sensitive to stellar mass. Indeed Fulton et al. (2017) report that the valley remains the same for all stars in their sample.

On the other hand, the actual morphology of the valley, i.e., the relative height of the two peaks, the widths of the peaks, the shape of the valley in period space, etc., depends on stellar mass. When plotting planet radii for a population that have a spread in host masses, it is easiest to see the evaporation valley when the horizontal axis is chosen to be the period, not bolometric insolation as commonly assumed. As eq. (20) shows, if we assume that $L_{\rm sat}\propto M_{*}$, the evaporation timescale scales with orbital period and stellar mass as,

[TABLE]

i.e., a weak dependence on stellar mass. Alternatively, if one takes the horizontal axis as the stellar insolation ($I\propto L_{*}/a^{2}$, or equivalently, planet equilibrium temperature),

[TABLE]

This latter has a stronger mass dependence; as a result, the evaporation valley is less distinct (Fig. 10 of Fulton et al., 2017). Similarly, but to a lesser degree, is plotting by semi-major axis, with $t_{\dot{X}}\propto a^{2.12}M_{*}^{-1.19}$. Ultimately, the stellar mass independent variable to use is high energy exposure (the total received high-energy flux over the planet’s lifetime); however, the down side is high energy exposure is a model dependent quantity and not observationally accessible for individual planets. The current high energy flux received by a planet at billions of years evolution is not representative of its high energy exposure (e.g. Tu et al., 2015).

4.7. Two populations of planets?

Envelope stripping becomes decreasingly important as one moves away from the star. For our assumed population of gas-rich planets, this means, going outward, there should be fewer bare cores, and these bare cores should be lower in mass (smaller in size). For example, Fig. 5 (the dotted black line) indicates that there should be no bare planets of mass $3~{}{\rm M}_{\oplus}$ ($R_{c}=1.3~{}{\rm R}_{\oplus}$) beyond orbital periods of $30$ days and no bare planets of mass $1\,$M⊕ beyond orbital periods of $60$ days (see also Lopez & Rice, 2016).888The variation of $\eta$ (§4.5) may extend this value out to $50$ days (Fig. 9) for bare 3 M⊕ planets and out to $100$ days for bare 1 M⊕ planets. Furthermore, Owen & Jackson (2012) showed that the hydrodynamical outflow we posit here do not happen for planets much outside $30$ days. Rather, the evaporative outflows is so rarefied the gas stops being collisional before the sonic point is reached and the mass-loss rate is instead determined by Jeans’ escape and falls much below the hydrodynamic value.

Interestingly, despite the selection effects against detecting small planets at large distances, they appear to be present in the Kepler sample.999The CKS sample does not contain enough detection in this range, and it is unclear what the stellar masses of individual planets are. The radius error bars remain large. And the abundance of such planets are not yet solidly established. However, if future observations confirm that such planets are abundant, this will require a separate population of planets than that posited here. This population is born bare.

5. Conclusions

We have developed a minimal analytical model which allows us to efficiently follow the evolution of low-mass planets, under the combined effects of cooling contraction and mass loss by evaporation, aside from illuminating the controlling physics. This model show that the mass-loss timescale peaks at around where planet sizes are doubled by their H/He envelopes, also where the envelope mass is of order a few percent. The timescale drops below this value because, while the envelope becomes more tenuous, the planet radii remain largely constant and so do the photo-evaporating fluxes they receive. The timescale also drops above this value because the planet swells up faster than the addition of envelope mass.

As a result, photoevaporation naturally gives rise to a final planet distribution that is bimodal in radius, peaking at the naked core size and twice its value. This then explains the observed radius “valley” in the California Kepler Survey sample (Fulton et al., 2017), as well as the steep fall-off of planets beyond $\sim 3R_{\oplus}$ in the general Kepler catalogue, for a single population of planets that are born with at least a few percent of H/He envelopes (“water-worlds” with water/steam envelopes are ruled out), this amount of gas lies above the most optimistic estimate for outgassing (Elkins-Tanton, 2008; Elkins-Tanton & Seager, 2008) and suggests that the envelopes were accreted from the protoplanetary discs. In the latter scenario, a few percent or more envelope mass may be natural (Rafikov, 2006). It is roughly the mass of an adiabatic envelope maintained by a planet of a few Earth masses at $0.1$ AU (but see Lee et al., 2014; Ginzburg et al., 2016, for modifications).

The positions of the peaks and valley also lead us to the following conclusions: the planet masses can be described by a Rayleigh distribution with a mode at $3~{}{\rm M}_{\oplus}$, the cores have a composition that is ice-poor. Ice-rich cores move the gap too far above the observed value. The absence of icy cores indicates local assembly of the Kepler planets, as opposed to formation beyond the ice-line (followed by large-scale migration). Furthermore, as most of the erosion occurs before 100 Myrs, the Kepler planets must have reached their current orbital locations well before that time. This rules out the possibility that they were migrated in at late times by dynamical processes (e.g. high-eccentricity migration) and also rules out the suggestion that the dearth of intermediate sized planets on short orbits (the “sub-jovian pampas”) is entirely due to tidal stripping during dynamical migration (Matsakos & Königl, 2016). Lastly, current data allow the silicate-iron cores to have a wide spread in iron fraction.

Planets that can be stripped bare have vanishingly small mass at larger orbital periods. Our energy-limited model does not produce planets of size $\sim 1.3~{}{\rm R}_{\oplus}$ (or $3~{}{\rm M}_{\oplus}$) much beyond $30$ days. We investigate the impacts on this (and other) prediction by factors such as atmosphere metallicity, core luminosity, stellar mass, and evaporation efficiency. None made much difference except perhaps the last one. In our ad-hoc model where the efficiency is higher for less bound atmospheres, we are able to extend the above period limit to $\sim 50$ days. This gives hope that a large enough observational sample could discriminate various evaporation models.

Observationally, there is some evidence small planets appear to exist at orbital periods longer than $50$ days in the overall Kepler sample. They are possibly numerous given their low detection probability. If this is confirmed by future observations, then our one-population model, though appearing to explain many of the observed features inward of $100$ days, may be insufficient. There may need to be another population of low-mass planets that are born with essentially no envelopes. This separate population will need to have masses $M\leq 3~{}{\rm M}_{\oplus}$, so as not to fill in the evaporation valley. Orbital instabilities have been suggested to be common among Kepler systems (e.g., Pu & Wu, 2015). Giant impacts among short period planets will happen at speeds well above their surface escape velocities (Agnor & Asphaug, 2004; Marcus et al., 2009), and will therefore disperse much of the original planetary envelopes. Since these impacts will occur after the natal disk has dispersed, they may give rise to this new population.

Future improvements are needed to further solidify results here. Theoretically, our results are obtained using a minimal analytical model. This was a deliberate choice in order to provide a basic understanding for the origin of the evaporation valley. We have calibrated our model against mesa calculations, but it still contains a number of assumptions (e.g., constant evaporation efficiency). Work should be carried out with more accurate numerical planetary structure models and more physically motivated evaporation model to determine the exact nature of the birth properties of short-period, low-mass exoplanets. Observationally, planet properties will be refined further by future precision stellar data, allowing more detailed comparison. We have adopted a mass function for the Kepler planets that is Rayleigh with a mode of $3M_{\oplus}$. This prediction should be tested by future mass measurements.

Notes in proof. While this paper was in review, Jin & Mordasini (2017) submitted a paper which used numerical modelling to reach much of the same conclusions about the composition of the Kepler planets as we do here. An interesting development comes from Dong et al. (2017) where they reported a population of “Hop-tunes”, Neptune-sized planets at close distances from their host stars, where we predict most planets should have been evaporated to bare cores. Intriguingly, this population only exists around metal-rich stars. This presents currently an unsolved puzzle.

The authors are grateful to the referee for comments which improved the manuscript. We thank Nikhil Mahajan for performing some of the early explorations. We are grateful to BJ Fulton for permission to reproduce his figure and useful discussions. We acknowledge Tim Morton, Lauren Weiss and Josh Winn for helpful insights. JEO acknowledges support by NASA through Hubble Fellowship grant HST-HF2-51346.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. YW acknowledges support from NSERC.

Appendix A The dimensionless integrals $I_{1}$ and $I_{2}$

In § 2.1.1 we have introduced dimensionless integrals of the form:

[TABLE]

In the limit $\Delta R\ll 1$, or $x\sim 1$ then $I_{n}$ is independent of $n$ and can be approximated as:

[TABLE]

This result implies that the ratio $I_{1}/I_{2}\approx 1$ for the case of a thin envelope. For the case of a puffy envelope, we need to consider how the $I_{n}$ varies with $R_{c}/R_{p}$. For $\gamma=5/3$, the choice appropriate in our model, the integrand of eq. (A1) is, to first order, $\propto x^{n-3/2}$ at small $x$. Therefore, for $n>1/2$ the integral is not dominated close to $x\sim R_{c}/R_{p}$ and eq. (A1) becomes independent of $R_{c}/R_{p}$ as the envelope becomes large and $R_{c}/R_{p}$ becomes small. Thus in the limit of large atmospheres ($\Delta R>1$), both $I_{1}$ and $I_{2}$ are constant and as such their ratio $I_{1}/I_{2}$ is also constant. In our model we require $I_{2}$ and the ratio $I_{1}/I_{2}$ explicitly. The full numerical solution of $I_{2}$ and the approximations discussed above are shown in the left-hand panel of Figure 11 and the full numerical solution of the ratio $I_{1}/I_{2}$ is shown in right-hand panel Figure 11.