© The Authors 2023

1 Introduction

The analysis of Kepler satellite data has revealed a dearth of 1.9 R_⊕ planets, often referred to as a valley or gap between the two populations of sub-Neptunes and super-Earths (Fulton et al. 2017; Van Eylen et al. 2018). Atmospheric escape models had predicted this dearth as an evaporation valley prior to the observational discovery (Lopez & Fortney 2013; Owen & Wu 2013; Jin et al. 2014). While the properties of the valley are now observationally quite well known, its origin is still debated. The leading hypotheses are X-ray and extreme-ultraviolet (XUV)-driven atmospheric photoevaporation (Owen & Wu 2013; Jin et al. 2014) and atmospheric loss driven by the cooling of the core (Gupta & Schlichting 2019). Alternatively, it might also be a direct imprint of formation, separating dry planets that have formed inside the ice line from volatile-rich ones that have migrated in from beyond the ice line (Venturini et al. 2020; Izidoro et al. 2022), or a consequence of impact-driven atmospheric erosion (Wyatt et al. 2020). It might also be caused by primordial gas accretion alone (Lee et al. 2022).

However, even within the context of just the atmospheric escape models, the details of the atmospheric mass loss driving the evolution, such as the different escape regimes (boil-off, blow-off, and Jeans escape) and limiting physical processes such as energy- or radiation-recombination-limited escape, are still poorly understood (for recent reviews on atmospheric escape, see Johnson et al. 2008; Zahnle & Catling 2017; Owen 2019).

This complexity is not surprising as escaping atmospheres in the Solar System (e.g. Lichtenegger et al. 2010) are also not yet fully understood in spite of in situ observations due to the complexity of molecular kinetic interactions, which include hydrodynamically escaping atmospheres. It is important to note that the main escape mechanisms in the Solar System are of a non-thermal nature; whereas, for close-in exoplanets, thermal escape mechanisms dominate.

Examples of atmospheric escape in the Solar System include hydrodynamic escape in the Earth’s H exosphere (Blamont et al. 1975), on early Earth and Venus (Watson et al. 1981), on Titan (Johnson et al. 2016) and a blend of Jeans and energy-limited escape on the Kuiper belt objects (Johnson et al. 2015), as well as plasma-driven escape on Mars (Jakosky et al. 2018; Leblanc et al. 2019) and Mercury (Gamborino & Wurz 2018; Thomas et al. 2004).

While the latter is also conceivable at a close-in irradiated exoplanet system (Oza et al. 2019; Gebek & Oza 2020), especially given the correlation in X-ray luminosity (McDonald et al. 2019), here we focus on the former mechanism, thermally driven hydrodynamic escape. We note that there is a growing interest in assessing the atmospheric escape of young bodies as well, such as protoplanets where the atmosphere is thought to be sourced directly from a magma ocean (e.g. Charnoz et al. 2021).

Energy-limited escape has long been used to approximate hydrodynamic escape to first order (Watson et al. 1981). It has the advantage of simplicity, hiding the complex physics in the evaporation efficiency factor η_XUV. Especially in planet evolution calculations, η_XUV is often assumed to be a constant, in contrast to the results of direct hydrodynamic simulations (e.g. Salz et al. 2016; Owen & Wu 2017). Nevertheless, the energy-limited approximation can fare rather well compared to some extensive kinetic simulations (Johnson et al. 2013) up to a critical threshold in the reduced heating rate. Above this threshold, when the escape transitions to a more Jeans-like regime (Volkov et al. 2011), the energy-limited escape approximation overestimates the escape rate by orders of magnitude (Salz et al. 2015). A second limitation is that at high EUV fluxes, the escape becomes radiation-recombination-limited rather than energy-limited (Murray-Clay et al. 2009). Third, even at intermediate XUV fluxes, the energy-limited approximation is not applicable for planets with a particularly high or low gravitational potential (Krenn et al. 2021). Finally, in the initial evolutionary phase of planets immediately after the dissipation of the natal protoplanetary gas disk, escape of primordial H and He envelopes is driven by a combination of low gravity and high atmospheric temperatures. This leads to very vigorous boil-off (Stökl et al. 2015; Owen & Wu 2016; Fossati et al. 2017), which is also neglected in a purely XUV-driven, energy-limited approach.

Therefore, in light of recent observations characterising the radius valley (Van Eylen et al. 2018; Martinez et al. 2019; Petigura et al. 2022; Chen et al. 2022; Ho & Van Eylen 2023), we seek to statistically test both the energy- and recombination-limited model as well as a direct numerical treatment of hydro-dynamic escape (Kubyshkina et al. 2018) which overcomes the assumptions and limitations of the energy-limited formula, against these observational constraints. For this study, we thus worked under the assumption that the valley is a consequence of atmospheric escape.

Regarding the observational data, we make use of the analysis of the California–Kepler Survey in tandem with parallaxes from the Gaia mission (Gaia Collaboration 2018) by Martinez et al. (2019) and Petigura et al. (2022). The latest Gaia release is essential since the new parallaxes provide a more accurate determination of planetary radii on a population-wide scale. We also compared our simulations to the observations from Van Eylen et al. (2018) which predates the second Gaia data release, but used astroseismology to determine accurate stellar parameters. Finally, we also used the valley locus as determined by Ho & Van Eylen (2023) based on short cadence Kepler data.

Our paper is organised as follows: in Sect. 2 we describe our theoretical model and describe the simulation setup in Sect. 3. In Sect. 4, using a rectangular grid of initial conditions, we compare the locus and slope of the evaporation valley in a radius–period diagram and demonstrate on a population-wide level how the hydrodynamic escape model, but not the energy-and recombination-limited model, leads to excellent agreement with the observed slope. The physical reason for this eventually become clear in Sects. 4.2 and 4.3 where we study selected individual evolutionary tracks, which highlight the differences between the two models, and we compare the models on the entire grid, respectively.

In Sect. 5, we determine the slope of the valley as found with the two evaporation models, though using initial conditions derived from the observed Kepler population instead of the rectangular grid. We end our paper with a summary and the conclusions (Sect. 6). In Appendix A we address the impact of the post-formation entropy on the valley locus.

2 Model

Our approach to modelling planetary evolution under the effect of atmospheric escape is two-fold: on the one hand, we use a simpler semi-analytical energy and radiation-recombination limited escape model of XUV-driven atmospheric photoevap-oration. This escape model was already used in the first paper of the series, Mordasini (2020). Here, we slightly update it, as described later in this section. The model is itself based on Jin et al. (2014) and Jin & Mordasini (2018). On the other hand, we now also use the tabulated escape rates obtained with a sophisticated numerical hydrodynamic escape model (Kubyshkina et al. 2018). In both cases, we couple these escape models to our model of temporal planetary interior evolution (cooling and contraction).

We then use this model to evolve a population of close-in planets. In both approaches, the interior evolution component of the calculations are performed with the COMPLETO21 planet evolution model that was described in details in the first paper (Mordasini 2020, hereafter Paper I).

This evolution model simulates the temporal thermodynamical and compositional evolution of the planet by solving the classical 1D spherically symmetric interior structure equations. The planets consist of an iron and silicate core described with the EOS of Seager et al. (2007) and a H and He envelope described with the EOS of Saumon et al. (1995). The atmosphere is described with an improved version of the double grey model of Guillot (2010), as described in Jin et al. (2014). This interior structure model yields, together with the mass loss, in particular the radius of the planet as a function of time as well as the remaining H and He envelope mass. In this work, as stated, we implement a new coupling to the hydrody-namic escape model described in Kubyshkina et al. (2018) that we shall next summarise along with the standard energy- and recombination-limited model that we used previously.

2.1 Energy- and radiation-recombination-limited escape model

Energy-limited (EL) escape assumes that the energy is lost most efficiently by gas expansion to space rather than conduction (downwards) or radiation (upwards to space). An in-depth analysis of the assumptions underlying the EL formalism can be found in Krenn et al. (2021), while a detailed description of our energyand radiation-recombination-limited escape model is given in Paper I.

In its simplest form, EL escape is written as $${\dot{M}}_{\text{EL}}~\frac{Q\left(R_{\text{a}}\right)}{U\left(R_{\text{p}}\right)},$$(1)

where U = GM_p/R_p is the specific binding energy of matter in the potential well of a planet of mass M_p and radius R_p with G the gravitational constant. R_a is the effective radius at which incoming radiation is absorbed on the planet. In our model, the radius where the EUV radiation is absorbed is calculated as described in Murray-Clay et al. (2009).

The complexity then arises in the heating rate Q, which is assumed to occur in the upper atmosphere. For XUV-driven escape, we can thus approximate the energy-limited escape due to upper atmospheric heating as $${\dot{M}}_{\text{U}}={\eta }_{\text{XUV}}\frac{\pi R_{\text{p}}{R}_{\text{a}}^2{F}_{\text{XUV}}}{G{M}_{\text{p}}{K}_{\text{tide}}},$$(2)

where we assume the only energy absorbed by the planetary envelope cross-section $\pi R_{\text{a}}^2$ driving escape is stellar X-ray and EUV radiation (collectively: XUV) which is written as the flux at the planet’s position F_XUV. It is further assumed that only a fraction of the total flux of the star drives mass loss, given by the evaporative efficiency factor η_XUV. This efficiency factor is, as discussed above, problematic as it oversimplifies the heating and cooling specific to each planet. Finally, the K_tide factor corrects for the gravity due to the stellar tide described in Erkaev et al. (2007).

In contrast to Paper I, where a constant η_XUV was assumed, we now use an η_XUV that depends on the escape speed v_esc, as suggested by approximate fits to the mass-loss simulations of Owen & Jackson (2012). The same functional form was also used in Owen & Wu (2017) and Rogers & Owen (2021). For the specific values of the parameters in Eq. (3), we use the ones which were found in Wu (2019) to lead to the best reproduction of the observed Kepler planet population: $${\eta }_{\text{XUV}}=0.17{\left(\frac{{\upsilon }_{\text{esc}}}{23\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^{-0.42}.$$(3)

These values are consistent with the ones found in Rogers & Owen (2021).

In practice, it is however found that η_XUV remains in the range of 0.1−0.3 because of the small exponent, and because the escape speed does not change by orders of magnitudes. Consequently, the differences to a simulation with a constant η_XUV are very limited. In particular, the slope of the valley is not affected by this modification and remains virtually identical to the value found in Paper I with a constant η_XUV. In this paper, it is also found that fixed globally higher or lower values of η_XUV do not affect the slope, but rather shift the valley up and down as a whole. This behaviour is in turn in perfect agreement with the predictions of the analytical model derived in Paper I (see Eq. (36) in that work).

As described in Jin et al. (2014) and Jin & Mordasini (2018), heating by UV and X-rays are treated separately in the model, using the criterion of Owen & Jackson (2012) to identify the dominant process. In the radiation-recombination-limited (RR) regime (Murray-Clay et al. 2009) that occurs at high EUV fluxes, the escape rate is given by the equilibrium of photoionisation with radiative recombination. In this regime, we closely follow Murray-Clay et al. (2009) to calculate the escape rate. The final escape rate is taken to be the minimum of the energy-limited and the radiation-recombination-limited escape rates (Lopez 2017).

The fact that the numerical results obtained with this evaporation model can be very well understood with an analytical model based on the energy-limited formula only (Paper I), shows that the importance of the recombination-limited regime is small for the planets studied here.

2.2 Hydrodynamic escape model

To estimate atmospheric escape within a more sophisticated direct hydrodynamic approach, we employ the grid of planetary upper atmosphere models presented by Kubyshkina et al. (2018). The grid consists of roughly 7000 models, each corresponding to a planet, and covers the following parameter space: planetary mass (Mp) between 1 and 39 Earth masses; planetary radius (R_p) between 1 and 10 Earth radii; planetary equilibrium temperature (T_eq) between 300 K and 2000 K; stellar mass between 0.4 and 1.3 solar masses; and XUV flux between 0.4 and 10⁴ the one experienced by present Earth because of solar irradiation, with values scaled for the specific stellar masses.

The range of orbital separations covered by the grid was set on the basis of the stellar mass and planetary equilibrium temperature, thus stellar radius (R_*) and effective temperature (T_eff). R_* and T_eff were derived considering the range of radii and effective temperatures covered by a star of each considered mass along the main-sequence on the basis of stellar evolutionary tracks (Yi et al. 2001). Considering all stellar masses, the orbital separation ranges between 0.002 and 1.3 AU.

The basic hydrodynamic model used to construct the grid is an updated version of the model developed in Erkaev et al. (2016). It considers a pure hydrogen atmosphere subject to heating and cooling processes, including radiative Ly-α cooling following Yelle (2004) and ${\text{H}}_3^{+}$ cooling following Miller et al. (2013) as well as adiabatic cooling (see details in Kubyshkina et al. 2018). These cooling processes are not explicitly included in the energy-limited approximation but are (incorrectly) supposed to be captured in the (constant) efficiency factor η_XUV introduced in Sect. 2.1.

The model numerically solves a full set of hydrodynamic equations, including energy and momentum conservation laws and continuity equations accounting for the full atmospheric hydrogen chemistry comprising dissociation, recombination, and ionisation. The complete list of chemical reactions is given in Kubyshkina et al. (2018). The model does not account for the presence of metals, which could induce additional heating and/or cooling that have been shown to be effective for ultra-hot primary (Fossati et al. 2021), secondary (e.g. Earth-like; Johnstone et al. 2018) and rock vapour (Ito & Ikoma 2021) atmospheres. However, conditions at planets in the grid are such that condensation may occur in the lower atmosphere, limiting the penetration of heavy elements in the upper atmosphere (see, e.g. Charnay et al. 2021). At the early evolution stages, when the extreme atmospheric escape takes place, the small amount of metals in a hydrogen-dominated upper atmosphere are dragged away by the hydrogen outflow and ultimately have little impact on the atmospheric mass-loss rates (Zahnle & Kasting 1986; Hunten et al. 1987; Odert et al. 2018; Lammer et al. 2020), while at the later stages, metal abundances can be fractionated by more moderate hydrodynamic escape or non-thermal escape processes (e.g. Gronoff et al. 2020). Even more importantly, one has also to consider that metal abundances may vary significantly between individual planets, already from the formation stage. Therefore, our current knowledge remains limited regarding the precise composition of sub-Neptunes, and the uncertainty in metal abundances does not enable one to place reasonable assumptions valid for planets spanning over a wide parameter space. Furthermore, the possibility to include metals into consideration is limited from the practical numerical side; the computational costs of the hydrodynamic models allowing for a proper metal treatment (i.e. including a detailed chemical framework and a photoionisation treatment including the explicit calculation of the energy levels populations) are at the moment still too high for computing a large and dense grid of mass-loss rates similar to that used in the present study.

The boundaries of the model are the photospheric radius of the planet (lower boundary) and its Roche lobe (upper boundary) $$R_{\text{roche}}=a{\left[\frac{M_{\text{p}}}{3\left(M_{\text{p}}+M_{\ast }\right)}\right]}^{1/3},$$(4)

where a is the planet’s orbital distance and M_*, the stellar mass. The model accounts for stellar heating in two wavelength intervals: EUV and X-ray ranges, assuming that the integrated flux of each range is emitted at a single wavelength (60 and 5 nm, respectively). The heating is included into the energy conservation equation as an external source given by $$Q_{\text{m}}=\frac{1}{2}{\eta }_{\text{XUV}}^{*}{\sigma }_{\text{m}}{n}_{\text{H}}{\displaystyle {\int }_0^{\frac{\pi }{2}+\arccos \left(\frac{1}{r}\right)}{J}_{\text{m}}\left(r,\theta \right)\cdot \sin \left(\theta \right)\text{d}\theta ,}$$(5)

where m stands for either X-ray or EUV radiation, σ_m is an absorption cross-section of hydrogen for the specific wavelength, n_H is the hydrogen (H + H₂) density, r is the distance to the planetary centre, and J_m(r, θ) is a function with spherical coordinates describing the spacial variations of the XUV flux due to atmospheric absorption $$J_{\text{m}}\left(r,\theta \right)=\exp \left(-{\displaystyle {\int }_r^{R_{\text{roche}}}{\sigma }_{\text{m}}{n}_{\text{H}}\left(\xi \right)\sqrt{{\xi }^2-r^2\sin \left(\theta \right)\xi \text{\hspace{0.17em}d}\xi }}\right),$$(6)

which is approximately equivalent to the optical depth at θ = 0.

We note that the ${\eta }_{\text{XUV}}^{*}$ in Eq. (5) for the hydrodynamic model is not the same as η_XUV given by Eq. (3) for the energy-limited model. This is because ${\eta }_{\text{XUV}}^{*}$ does not account for any additional cooling processes or other physical mechanisms supposed to be captured by (or hidden in) η_XUV, as they are included self-consistently in the hydrodynamic model (Kubyshkina et al. 2018). Instead, ${\eta }_{\text{XUV}}^{*}$ accounts solely for the efficiency of the photoionisation heating, and is not an overall evaporation efficiency as η_XUV in the energy-limited model.

Given that a self-consistent calculation of ${\eta }_{\text{XUV}}^{*}$ is currently too time-consuming for computing a large grid, it was set to be equal to a constant value of 15%, which is a reasonable assumption for the considered M_p range (e.g. Shematovich et al. 2014; Salz et al. 2016). We note that despite this compelled simplification, the hydrodynamic code remains a superior model relative to the energy-limited approximation, as the latter approach relies on many more assumptions than just a constant heating efficiency and the absence of the explicitly modelled radiative cooling processes. In particular, it omits the contribution from the thermal energy of the planet atmosphere and the stellar VIS/IR irradiation (as discussed below) and makes crude assumptions on the atmospheric structure, which is calculated self-consistently by the hydrodynamic model (for a more thorough discussion, see Krenn et al. 2021).

Typically, for hydrodynamic planetary and stellar wind models, the initially subsonic outflow (we set the bulk velocity V_bulk equal to zero at the lower boundary) is accelerated to supersonic velocities before the flow reaches the Roche lobe. Within our grid of models, it happens typically at a distance of a few planetary radii. To ensure that the atmospheres of planets in the grid remain collisional throughout the simulation, we calculate the Knudsen number for each point of the atmospheric profiles a posteriori. The atmospheric mass-loss rate is finally defined as the flow through the sphere of radius r in a unit of time (m_Hn_H(r)V_bulk(r)) multiplied by the surface of this sphere. As the outflow is continuous, for the computation of the mass-loss rate the specific distance r is not relevant (except for the small region at the lower boundary), but for convenience it is taken at the Roche radius.

The predictions of our model are comparable to those made by other hydrodynamic models, including the more sophisticated ones (such as those calculating self-consistently the heating efficiency in various approaches as Murray-Clay et al. 2009 and Salz et al. 2016, and models accounting for the detailed spectral energy distribution as Guo & Ben-Jaffel 2016, or 3D geometry as Carolan et al. 2021). Further details about the physical model and the grid, including the comparison to observations and to the results of other literature models, can be found in Kubyshkina et al. (2018).

By construction, the hydrodynamic model accounts for Jeans escape, XUV hydrodynamic escape, and boil-off escape regimes. This is, as we shall see below, of central importance for the results for the valley found here. The model also transitions smoothly from one escape regime to the other depending on the system parameters. To ease distinguishing between the latter two regimes, it is convenient to employ the restricted Jeans parameter (Fossati et al. 2017), which is a combination of the physical planetary parameters and is defined as $$\text{Λ}=\frac{G{M}_{\text{p}}{m}_{\text{H}}}{k_{\text{B}}{T}_{\text{eq}}{R}_{\text{p}}},$$(7)

with m_H the mass of the hydrogen atom and k_B the Boltzmann constant.

Planets with a Λ smaller than 15–35 are in the boil-off regime, where the escape is driven by the atmospheric thermal energy and low planetary gravity (Owen & Wu 2016; Fossati et al. 2017; Ginzburg et al. 2016). The specific critical value depends on the stellar mass and orbital separation. Such planets are typically just released from the protoplanetary gas disk. A Λ of 20 is for an isothermal gas identical to the condition derived by Owen & Wu (2016) for the occurrence of boil-off, namely that the planet radius is larger than about 0.1 times the Bondi radius.

For the boil-off regime, it is crucial that the hydrodynamic model also accounts for the stellar continuum (dominated by VIS/IR) heating that can drive escape, in contrast to the energy-and recombination-limited model that is driven by XUV heating only. The continuum heating is implicitly included by fixing the temperature at the lower boundary equal to T_eq. We have verified that the photospheric temperatures of our model planets as predicted by the interior structure model is always very close to T_eq. The largest difference (a temperature that is about 4% higher) occurs for the most massive planets we model at the beginning of the simulations, which is due to the contribution of the intrinsic luminosity. Overall, the difference is, however, much smaller and generally less than 1%.

Studies comparing the hydrodynamical model used here with the energy-limited escape have found the following (Kubyshkina et al. 2018; Krenn et al. 2021): For planets with Λ less than about 20, the energy-limited formalism on average severely underestimates mass-loss rates, because it lacks the continuum (VIS/IR) heating. For higher Λ, the energy-limited rate provides an upper limit on the mass-loss rate, with significant overestimations possible depending on a planet’s gravitational potential.

Model outputs for each planet in the grid are profiles of the main atmospheric parameters, which allow deriving the effective radii of the stellar XUV absorption, and atmospheric escape rates. To finally obtain the mass-loss rates for any planet during its evolution, we linearly interpolate among the grid points. Our interpolation scheme is simpler than the one in Kubyshkina et al. (2018), but allows to fully exploit all grid data including the borders of the tabulated regions.

Fig. 1 Temporal evolution of the bolometric (blue), X-ray (green), and EUV luminosity (brown) of a 1 M⊙ star as assumed in our model. The bolometric luminosity is divided by a factor 1000 to bring it on a similar scale as LX and LEUV. The two black bars near 4.5 Gyr show the range of our Sun’s LX over the course of a solar cycle. The grey dashed lines show for comparison the LX of Tu et al. (2015) for the 10th, 50th, and 90th percentiles of the rotational distribution.

2.3 Stellar XUV luminosity as a function of time

A modification of our theoretical model relative to Paper I is the usage of a more recent description of the stellar XUV luminosity as a function of time. In Paper I, we used the data of Ribas et al. (2005). In the updated model, we use instead McDonald et al. (2019). These authors compiled observationally derived relations extracted from several studies (Jackson et al. 2012; Shkolnik & Barman 2014) of the X-ray luminosity of stars as a function of time and stellar type. We use their mean X-ray luminosity as function of time L_X(t) and convert it into the extreme UV-luminosity L_EUV(t) with the relation of Sanz-Forcada et al. (2011).

Figure 1 shows LX, LEUV, and the bolometric luminosity Lbol of a 1 M star in our model. At young ages, the XUV luminosity is of the order of 10⁻³ the bolometric luminosity, as expected (e.g. Güdel 2020), and approximately constant in time, except for a certain drop at around the time (40 Myr) when the star reaches the main sequence. Afterwards, it decreases approximately following a power law. At 4.6 Gyr, the predicted L_X is compatible with the Sun’s measured L_X at its activity maximum.

We compare our model with the one of Tu et al. (2015). They calculated the L_X predicted for stars on the 10th, 50th, and 90th percentiles of the stellar rotational distribution. We see that our relation is similar to their 50th percentile case. At the earliest epochs, our L_X is about a factor 2 lower than theirs, while at high ages, the fall off is a bit faster in Tu et al. (2015). The impact of varying the stellar XUV on the locus of the evaporation valley was studied in Paper I (see also Ketzer & Poppenhaeger 2023).

3 Procedure

We used the model to evolve a large number of close-in low-mass planets. To set their initial (post-formation) properties, we followed two approaches: a rectangular grid, and initial conditions derived from the planetary population detected of the Kepler satellite. The rectangular grid allows us to see clearly (but also under idealised assumptions) the population-wide imprints of the two evaporation models. The initial conditions derived from the Kepler population give us an understanding if these imprints remain visible also when the initial conditions are more complex, in particular when there is a spread in the post-formation envelope mass at a fixed core mass.

3.1 Rectangular grid of models and initial conditions

Our first approach was the same as in Paper I: for the two escape models, we simulated a rectangular grid of 6000 planets each, equally spaced in semi-major axis a and core mass M_core ranging from 0.01 to 0.6 AU in 0.01 AU increments and from 1 to 20 M_⊕ in 0.2 M_⊕ increments. We let these planets evolve from 3 Myr (a typical lifetime of protoplanetary disks, Mamajek 2009) to 10 Gyr around a 1M_⊕ star. With this data, we can analyse the slope and temporal evolution of the evaporation valley.

The initial (post-formation) H and He envelope mass M_env,0 was estimated as $$M_{\text{env,0}}=0.024M_{\oplus }{\left(\frac{M_{\text{core}}}{1M_{\oplus }}\right)}^{2.23}{\left(\frac{a}{1\text{\hspace{0.17em}AU}}\right)}^{0.72}.$$(8)

As described in Paper I, this relation was found as a typical mean value from planet formation simulations by Mordasini et al. (2014) based on the core accretion paradigm which find the envelope mass similarly as in Pollack et al. (1996), but include many additional effects such as orbital migration and disk evolution. The results of Paper I employing the XUV-driven energy-and recombination-limited escape model indicate, however, that using a different initial envelope mass within plausible ranges should not strongly influence the location of the valley. An additional argument for this weak dependency in the context of boil-off is that if the initial envelope mass is larger, then the radius is larger and thus the boil-off escape is larger as well. Therefore, at the end of the boil-off phase, an initially larger and an initially smaller planet end up with a similar radius because the larger planet had a stronger escape compared to the smaller one (Kubyshkina et al. 2020). We nevertheless investigate the impact of the initial envelope mass further in Sect. 5. The initial intrinsic luminosity of the planets was also estimated as a function of core and envelope mass based on the same formation simulations. This is the same approach as in Paper I. In Appendix A, we study the impact of different post-formation luminosities or entropies, finding an only weak influence on the valley locus.

Regarding their bulk composition, all cores have an Earthlike composition with a 2:1 silicate-to-iron mass fraction. Such a composition is in agreement with the composition of planets below the valley (Owen & Wu 2017; Jin & Mordasini 2018).

Fig. 2 Initial conditions derived from Kepler survey. Left and middle panel: histogram of the distribution of the core masses and initial (post-formation) envelope mass fractions Menv,0/Mcore for the comparison with the Kepler planet population. Both these distributions are taken from Rogers & Owen (2021). Right panel: average detection probability of the Kepler satellite pdet as function of orbital period and planet radius (Petigura et al. 2018). The final probability to find a planet is given by pdet times its geometric transit probability ptr.

3.2 Initial conditions derived from Kepler survey

The initial conditions on the rectangular grid are not tuned to reproduce the observed Kepler planet population, which makes the comparison with observations less straightforward. The grid also assumes in an idealised way that there is a unique value of the initial envelope mass as a function of core mass and orbital distance. Formation models (e.g. Mordasini et al. 2014), but also inference analyses of the post-formation properties of the Kepler planets (Rogers & Owen 2021), indicate in contrast a spread in post-formation envelope masses.

Our second approach for the initial conditions was thus to adopt distributions for the orbital period, core mass, and initial envelope mass that have been derived from fitting the observed properties of the close-in low-mass population found by the Kepler satellite through inference analyses (Gupta & Schlichting 2020; Rogers & Owen 2021). In these works, the core and initial envelope mass distributions that were derived lead – after evolution under the effects of core-powered mass loss and photo-evaporation, respectively – to a synthetic population that agrees with the observed period-radius distribution (the CKS data, Fulton & Petigura 2018). The period distribution is also derived from the Kepler observations and given as (Rogers et al. 2021) $$\frac{\text{d}N}{\text{d\hspace{0.17em}log\hspace{0.17em}}P}\propto \{\begin{array}{ll}{P}^2,\hfill & P<8\text{\hspace{0.17em}days,}\hfill \\ \text{constant,}\hfill & P>8\text{\hspace{0.17em}days.}\hfill \end{array}$$(9)

The core mass distribution we adopted is the one inferred in Rogers & Owen (2021) in their preferred Model III. It peaks at a core mass of about 4 M_⊕, with a tail extending to about 100 M_⊕. The post-formation envelope mass fraction was also taken from this source. It is a distribution peaking at an envelope mass fraction of about 4%, but covering a significant range. In contrast to the theoretical relation (Eq. (8)), the envelope mass fraction is here an independent quantity. Both these distributions are shown in the left and middle panel of Fig. 2. With these initial conditions, we calculated the evolution of 37242 and 37416 planets from 3 Myr to 10 Gyr for the energy-limited and hydrodynamic evaporation model, respectively.

To understand if the imprints of the different evaporation models remain observable, we applied a simple synthetic detection bias of the Kepler satellite to the model output. In this way, we got the detectable synthetic population. For each synthetic planet, we computed the detectability as a function of planet size and orbital period. It has two components. The first component is the geometric transit probability p_tr. For it, following Petigura et al. (2018), we used that a randomly inclined planet on a circular orbit transits with an impact parameter b < 0.9 with a probability p_tr = 0.9R_⊙/a, where R_⊙ is the radius of the Sun (we only consider 1 M_⊙ stars in this paper). The second component is the detection probability p_det, which depends mainly on the S/N of the observations. Here we took the average p_det also from Petigura et al. (2018), which is based on the transit injection and recovery study of Christiansen et al. (2015). This is shown in the right panel of Fig. 2. The total probability is then the product of the two probabilities, p_tr × p_det. By comparing p_tr × p_det with a random number drawn from the standard uniform deviate, we obtained the detectable synthetic planets. To have enough detectable synthetic planets despite the low detection probability of the transit method, we oversampled 100 times, that is to say we ran through the list of synthetic planets 100 times, obtaining each time different detectable planets. This means that the same planet can end up several times in the final list of detectable planets. However, for the statistical analysis at hand, this is not an issue. In this way, we ended up with 114 634 and 113 465 detectable synthetic planets for the energy-limited and hydrodynamic model, respectively.

The initial condition distributions derived in Gupta & Schlichting (2020) and Rogers & Owen (2021) lead with their relative theoretical forward models to synthetic populations agreeing with the Kepler data. However, the distributions that the two works infer differ from each other. This reflects that these ‘fitting’ initial conditions are also a function of the forward model. Here, we use again another forward model (or even two, counting the two different evaporation models). Thus, we cannot expect that we find with our forward model in the end a detectable subpopulation agreeing equally well with the actual Kepler population. As shown in Sect. 5 our synthetic detectable populations do, however, still share key properties with the observed population, such as a bimodal radius distribution. We could, in principle, conduct a similar hierarchical inference process as Gupta & Schlichting (2020) and Rogers & Owen (2021) to derive our own fitting initial conditions. However, practically this would be difficult because of the much higher computational cost of our forward model compared to theirs. It would also be beyond the scope of this paper, which addresses the comparison of two evaporation models.

3.3 Quantifying the locus of the valley

Following the approach of several previous papers (e.g. Van Eylen et al. 2018; Lopez & Rice 2018; Jin & Mordasini 2018; Martinez et al. 2019; Rogers et al. 2021; Petigura et al. 2022; Ho & Van Eylen 2023), we quantify the valley locus with a power law. Normalising at an orbital period of 10 days, we express for the rectangular grid simulations the planetary radius Rp of the largest bare core (i.e. most massive planet which has completely lost its H and He envelope) at a given orbital period p as $$R_{\text{b}}\left(p\right)={\tilde{R}}_{\text{b}}{\left(\frac{p}{10\text{\hspace{0.17em}days}}\right)}^{\alpha }$$(10)

where ${\tilde{R}}_{\text{b}}$ is the value at 10 days and the slope is $$\alpha =\frac{\text{d\hspace{0.17em}log\hspace{0.17em}}{R}_{\text{p}}}{\text{d\hspace{0.17em}}\log p}.$$(11)

Using this definition, we are consistent with previous works on the same topic.

A power law dependency has also been analytically found by several theoretical works, for example by Owen & Wu (2017) and Paper I for photoevaporation or by Gupta & Schlichting (2019) for core-driven escape. These works show that the slope of the valley is a good indicator for the dependence of the evaporation rate on the planets’ distance from the host star and therefore the underlying evaporation mechanism.

It is important to specify that for the results obtained with the rectangular grid of initial conditions, R_b(p) is the lower boundary of the observed valley. In contrast, observational studies, but also our theoretical results obtained for initial conditions derived from the observed Kepler planets, report the middle of the valley.

For the initial conditions derived from the Kepler survey, the absence of a completely empty gap (or a well-defined largest bare core as a function of period as in the rectangular grid) means that we cannot as simply quantify the valley position as for the rectangular grid. On the other hand, compared to the observed population, where various statistical methods must be used such as support vector machines (e.g. Van Eylen et al. 2018) to determine the valley position and its slope, we are in the advantageous position to have a very large data set. We have therefore proceeded in the following simple way to determine the middle of the valley: We have binned the planets according to orbital period, with a bin width of 0.2 dex in log P, with partially overlapping bins at log P = 0.6, 0.7,… 1.7, 1.8 (or 1.9 for the unbiased case). For each bin, we represent the radius distribution with a kernel density estimate, and get the position of the gap centre (local minimum in the radius distribution) from the zero point of the derivative. This procedure is similar to the one employed by Petigura et al. (2022).

Finally, to determine ${\tilde{R}}_{\text{b}}$ and α from the simulations, we simply make, as in Paper I, a least-square power law fit to the largest bare core radius (for the rectangular grid) respectively gap centre (for the Kepler initial conditions) as a function of orbital period at a given age, typically at 5 Gyr.

4 Results for the rectangular grid

We present our planet evolution simulations by first examining the locus of the valley in the period-radius diagrams as predicted by the two evaporation models for the rectangular grid of initial conditions in Sect. 4.1. By performing case studies on individual planets in Sect. 4.2, we are then able to identify the cause of the distinct valley slopes. We then extend this analysis to the entire grid (Sect. 4.3).

4.1 Period-radius diagrams

Figure 3 shows the simulated grid in the orbital period-transit radius plane for the two escape models at an age of5 Gyr. Clearly visible in both is the evaporation valley running diagonally downwards, that is the gap in radius between the super-Earth planets whose envelopes have fully evaporated and the sub-Neptunes which still have an envelope. For the planets still possessing H and He, the colour code shows the fraction of the initial H and He envelope that has evaporated. The closer to the valley, the higher this fraction, as expected.

With the divide being very distinct, a good fit of the slope can be achieved. We note that there are some quasi-regular patterns emerging in the dots above the valley, and in the hydro-dynamic model, some linear patterns are visible. These patterns are simply the result from the regular grid, and, in the case of the hydrodynamic simulations, also consequences of the interpolation in the grid of tabulated evaporation rates. For the hydrodynamic model, this also translates in the largest bare core as a function of period not being a completely smooth power law function, as it is the case for the energy- and recombination-limited model. However, also for the hydrodynamic model, the simulations of the individual planets presented below exhibiting clear, physically understandable outcomes that are not dominated by interpolation artefacts, which, together with the small scales of the non-smoothness of R_b, mean that the interpolation does not significantly affect the quantities we are interested in (${\tilde{R}}_{\text{b}}$ and α).

In both models, there is a region with an absence of (sub-) Neptunian planets (i.e. planets with H and He) at periods smaller than 2 or 3 days. This corresponds to the evaporation (also called the sub-Neptunian) desert (Lundkvist et al. 2016; Mazeh et al. 2016; Bourrier et al. 2018). It should be noted that the specific period marking the onset of the desert in our simulations shown here depends also on the fact that the most massive core we simulate has a mass of 20 M_⊙. The minimum period is thus model dependent. To quantify the valley locus, we only considered the planets having the smallest period for which there are still planets with an envelope.

The comparison of the two panels shows that while both models lead to a very similar position of the valley at 10 days ${\tilde{R}}_{\text{b}}$ of about 1.8 to 1.9 R_⊕, they differ in the slope, that is in α as is directly visually apparent. In the energy- and recombination-limited model, the slope is clearly steeper than in the hydrody-namic model. This is one of the key outcomes of this study. In both panels, the magenta lines show the aforementioned power law fit to the largest, numerically found bare core as a function of period. It makes the shallower slope in the hydrodynamic case even more apparent.

In Table 1, we report the parameters of these fits, ${\tilde{R}}_{\text{b}}$ and α, together with the results of the observational studies of Van Eylen et al. (2018), Martinez et al. (2019), and Ho & Van Eylen (2023). The data confirm the visual impression: the difference in ${\tilde{R}}_{\text{b}}$ (1.84 and 1.88 R_⊕ in the energy- and recombination-limited and hydrodynamic model, respectively) is very small. This difference is comparable to, or smaller than, the observational error bars, and thus not significant.

Our theoretical values for ${\tilde{R}}_{\text{b}}$ are for the bottom of the valley, while the two observational studies are for the middle. Correcting for this difference would shift the theoretical ${\tilde{R}}_{\text{b}}$ to larger values by about 0.2−0.3 R_⊕ (Mordasini 2020; David et al. 2021), that is to a radius between 2.0 and 2.1 R_⊕. This is larger than the nominal observational values of about 1.9 R_⊕. This could be an indication that both models overestimate escape. A possible explanation for this could be that the interior structure and the escape models do not account for metals. Because the planets considered here orbit late-type stars, metals do not cause much heating, but may lead to significant cooling, which would lower the escape (e.g. García Muñoz 2007). A metal-enriched instead of a pure H and He envelope would also affect the interior structure model via the equation of state and the opacities. The impact of enriched envelopes was studied in Paper I. It was found that a gas with a metal mass fraction of 10% and 30% would lead to a downward shift of the valley of about 0.1 and 0.2 R_⊕. This would bring the theoretical predictions into better agreement with the observations. Such enrichments could be a consequence of the formation process (Fortney et al. 2013; Brouwers & Orme 2020) or result from evolutionary magma-hydrogen interactions at the core-envelope interface (Misener et al. 2023).

In contrast to ${\tilde{R}}_{\text{b}}$, the values of the slope α differ clearly between the two theoretical models: in the energy- and recombination-limited model, a slope of 0.18 is found, while for the hydrodynamic model, the slope is 0.11. This corresponds to a fractional difference of about 50%, a difference that is clearly observationally relevant, given the error bars reported in the observational studies. The α found in the updated numerical energy- and recombination-limited model used here is identical to the one derived analytically for a purely energy-limited model with a constant η_XUV derived in Paper I. This shows that the introduction of an escape velocity dependent η_XUV does not significantly affect the results.

Comparing with the observationally inferred values of α which are $-{0.09}_{-0.04}^{+0.02}$ in Van Eylen et al. (2018) and −0.11 ± 0.02 in both Martinez et al. (2019) and Petigura et al. (2022), we see that the energy- and recombination-limited model yields with 0.18 a slope that is clearly too steep by more than two or three sigma. The hydrodynamic model in contrast predicts a slope that is in excellent agreement with these observations.

The finding that a simple energy-limited model yields a steeper slope than a hydrodynamic model is not new: Owen & Wu (2017) had already found an α = −0.25 for their simple energy-limited model with a constant efficiency factor, while their full hydrodynamic evaporation model (Owen & Jackson 2012) without these assumptions yielded α = −0.16. Both these values are, however, steeper than the observed values, in contrast to our new findings with the Kubyshkina et al. (2018) hydrody-namic model. Another theoretical energy- and recombination-limited XUV photoevaporation model (Lopez & Rice 2018) found with α = −0.15 also a steeper slope than observed. This makes it worthwhile to further investigate the result found here.

The Owen & Jackson (2012) model predates studies such as Stökl et al. (2015), Owen & Wu (2016), and Fossati et al. (2017), which found the boil-off phase as the first escape regime occurring after the dissipation of the disk. Instead, their initial evaporation regime is X-ray driven, in contrast to our hydrodynamic model where boil-off is included.

We report here also the masses of the largest bare core as a function of period in the two models, found in the same way as the radii with least-square power law fits. This is an information of interest for radial velocity studies. One finds for the energy-and recombination-limited model $$M_{\text{b,enRR}}\left(p\right)=9.6{\left(\frac{p}{10\text{\hspace{0.17em}days}}\right)}^{-0.64}{M}_{\oplus }$$(12)

and for the hydrodynamic model $$M_{\text{b,hyd}}\left(p\right)=10.6{\left(\frac{p}{10\text{\hspace{0.17em}days}}\right)}^{-0.41}{M}_{\oplus }.$$(13)

These values essentially reflect the mass-radius relation of silicate (MgSiO₃)-iron planets. The analytical model of Paper I for an energy-limited model predicts for comparison that the exponent should be −0.66. This is very close to the value obtained numerically in the present paper (−0.64). In the hydrodynamic model, the exponent is as expected significantly lower (−0.41).

Fig. 3 Results of planet evolution simulations in the orbital period-transit radius plane at 5 Gyr (left panel: energy- and recombination-limited model; right panel: hydrodynamic model). Each dot represents a planet which is coloured based on the fraction of its initial H and He envelope mass that was evaporated. Grey symbols indicate a complete loss of the envelope, corresponding to sub-Neptunes that have evolved into super-Earths. The magenta lines show the fit to determine the slope α of the valley as represented by the largest super-Earth at a given period (i.e. the bottom of the valley). We note the shallower slope in the hydrodynamic model. The squares show two individual cases (one far planet in lilac, one close in black). Both start with identical initial conditions in the two models and are discussed in Sect. 4.2. They end up on different sides of the valley in the two models in a way to make the slope shallower in the hydrodynamic model.

4.2 Evolution of specific cases

The fact that ${\tilde{R}}_{\text{b}}$ is similar in the two models while the valley is clearly shallower in the hydrodynamic model means that the hydrodynamic model does not generate as large bare cores inside the 10-day period line as the energy- and recombination-limited model does. The opposite is true outside the 10-day period. This is quickly verified in Fig. 3.

Therefore, to understand the reasons for the different slopes, it is helpful to study two cases of individual planets. The first one is a far planet at a period of about 133 days, which remains a subNeptune in the energy- and recombination-limited model, but becomes a super-Earth in the hydrodynamic model (Sect. 4.2.1). The second one is a close planet at 3 days, which becomes a super-Earth in the energy- and recombination-limited model, but stays a sub-Neptune in the hydrodynamic model (Sect. 4.2.2). In other words, we study cases that end up on opposites sides of the valley in the two models in such a way that the slope is shallower in the hydrodynamic model. We select cases that are close to or at the upper boundary of the super-Earths for the two distances. In Fig. 3, these individual planets are shown with squares (black for the close case, lilac for the far case).

4.2.1 Distant planet

Figure 4 shows the temporal evolution of the evaporation rate, restricted Jeans parameter, transit radius, and total mass of a planet at an orbital distance of 0.51 AU (period of 133 days) with an initial total mass of approximately 3.86 M_⊕. The planet consists of a 3.6 M_⊕ silicate-iron core and a 0.26 M_⊕ H and He envelope.

The important result here is that in the hydrodynamic model, the planet completely loses its H and He envelope at about 2.5 Gyr, transforming the planet into a super-Earth, whereas in the energy- and recombination-limited model, the planet keeps about 60% of the initial envelope till the end of the simulations and thus remains a sub-Neptune. The final radii (at 5 Gyr) in the two cases are about 1.4 R_⊕ and 2.8 R_⊕ (the former equal to the core radius), typical for planets below and above the evaporation valley.

The reason for this different evolution can be seen in the evaporation rate (top left panel in Fig. 4). We see that the hydro-dynamic model initially predicts a much higher evaporation rate. At the very beginning, the evaporation rate is more than 2 orders of magnitude higher. It remains larger than the one in the energy-and recombination-limited model to about 40 Myr. The reason for the higher evaporation rate can be seen in the Jeans-escape parameter (top right panel): we see that initially, Λ is about 12. This puts the planet firmly into the boil-off regime (Krenn et al. 2021), which leads to the very high escape rates in the hydrody-namic model. During this time, the stellar continuum irradiation (mainly in VIS/IR), rather than the XUV irradiation, catalyses the escape. A tell-tale sign of this is the local maximum in the escape rate in the hydrodynamic model at about 30 Myr, which is caused by the local maximum of the star’s bolometric luminosity (see Fig. 1). During boil-off, the planet loses about half its H and He envelope in the first 3 Myr. The escape gradually transitions into XUV-driven escape at about 30–50 Myr. By this time, the radius of the planet has shrunk to a size that is comparable to one tenth of the Bondi radius (which is about 40 R_⊕), as predicted by Owen & Wu (2016). By this time, Λ has increased to about 30. At later times, the hydrodynamic model predicts a comparable or somewhat smaller escape rate than the energy- and recombination-limited model, as can be seen in the top left panel. However, by that time, the properties of the planets (namely the radii) have already diverged significantly between the two models, thus it is difficult to compare them. We do this later when we compare the escape rates at fixed planet properties.

To summarise, we see that for these distant low-mass planets, the hydrodynamic model predicts that comparatively more massive planets become super-Earth than in the energy- and recombination-limited, shifting the valley to larger radii because of boil-off. This evaporation regime is included in the hydrody-namic model, but not in the energy- and recombination-limited one, which makes the difference. This planet is actually a typical example of the first category of planets where the energy-limited approximation consistently fails by under-predicting the escape rate (Krenn et al. 2021): planets with low-to-intermediate XUV irradiation and low gravitational potential.

4.2.2 Close-in planet

The second case, a close-in massive planet, is shown in Fig. 5. This is a planet at 0.04 AU (orbital period of about 3 days), with an initial total mass of about 20.9 M_⊕. The initial envelope mass is 1.72 M_⊕. This may seem a small envelope mass for the significant core mass, but it is a consequence of the also very small orbital distance, that reduces in formation simulations the ability to accrete gas (see Eq. (8)).

Here, the key result is that the outcome is opposite to the distant case. As a matter of fact, in the energy- and recombination-limited model, the envelope is completely lost by 800 Myr. Instead, in the hydrodynamic model, the planet keeps an envelope to the end of the simulation at 5 Gyr; see Fig. 3. The remaining envelope mass at this time is actually tiny, but sufficient to lead to a radius of about 2.2 versus 2.5 R_⊕ for the two cases. In the case we study here, the difference of the two models is by construction small, since we have chosen a case that is just above respectively below the valley in the two models.

The reason why the models lead to such different outcomes can be seen in the top left panel of Fig. 6. We see that as in the case of the distant planet, the hydrodynamic model initially predicts a stronger mass loss, which is again due to boil-off. However, compared to the far case, the boil-off is here less extreme, leading to a difference in the evaporation rates of initially a bit more than one order of magnitude. This is due to the fact that as the planet already starts with a higher Λ of about 17. Already after about 10 Myr, the evaporation rates become comparable in the two models, and after about 30 Myr, the escape rate in the energy- and recombination-limited model is consistently higher than in the hydrodynamic model. Neither the radii nor the masses (dominated by the core mass anyway) of the planets differ strongly at this point. The similar M_p, R_p, and the identical T_eq in the two simulations imply that Λ is also similar¹. Therefore, one can directly compare the evaporation rates in the two models. The lower rate in the hydrodynamic model is thus not merely a consequence of different planet properties, but a genuine consequence of the more complex physics included in the hydrodynamic model, and more specifically in the different temperature structure compared to that assumed in the energy-limited approximation, as we shall further discuss in the following section. The difference in the predicted escape rates is not very large (factor 2 to 3), but integrated over time, this is sufficient for complete evaporation in the energy- and recombination-limited, but not the hydrodynamic model.

The underlying reason is that this planet is a typical example of the second regime where the energy-limited approach consistently fails by over-predicting the escape rate (Krenn et al. 2021), namely planets characterised by high XUV irradiation and high gravitational potential.

Fig. 4 Temporal evolution of the evaporation rate (top left), restricted escape parameter Λ (top right), transit radius (bottom left), and total mass (bottom right) for a planet at 0.51 AU and initial mass of approximately 3.86 M⊕. The hydrodynamic and the energy- and recombination-limited model are shown. We note the very high mass loss in the hydrodynamic model in the initial boil-off phase. It is sufficient to eventually cause the total loss of the envelope. In the energy- and recombination-limited model which lacks this phase, the planet can in contrast keep a significant fraction of its envelope. The planet therefore resides in the end on different sides of the valley in the two models (super-Earth for the hydrodynamic model, sub-Neptune for the energy- and recombination-limited model). This is shown by the lilac squares in Fig. 3.

Fig. 5 Temporal evolution of the evaporation rate (top left), restricted escape parameter Λ (top right), transit radius (bottom left), and total mass (bottom right) for a planet at 0.04 AU and initial mass of approximately 20.9 M⊕. The hydrodynamic and the energy- and recombination-limited model are shown. Despite the lack of the initial boil off phase in the energy- and recombination-limited model, the higher mass-loss rate in this model after about 30 Myr results in complete envelope loss, in contrast to the hydrodynamic model. The planets therefore reside in the end on different sides of the valley (super-Earth for the energy- and recombination-limited model, sub-Neptune for the hydrodynamic model). The underlying reason is the negligence of conduction in the energy- and recombination-limited model, leading to too high escape rates. This planet is shown by the black squares in Fig. 3.

4.3 Comparison of the escape models on the entire grid

The previous section demonstrated the different outcomes of two selected planets using the two escape models, and which effects played a significant role. We now generalise this comparison to the entire grid of planets we simulated (for another systematic comparison, see also Krenn et al. 2021).

Figure 6 compares the two models at an age of 50 Myr. At this time, boil-off in the hydrodynamic model has already ceased, but strong XUV-driven evaporation (because we are still at early times) is ongoing. Four panels are shown in the periodradius plane, colour-coding different quantities. In panels a, b, and c, the results of the energy- and recombination-limited model are shown with coloured dots, in panel d they show the results of the hydrodynamic model. In all cases we see the valley, which is, however, not yet at the same position as in Fig. 3, which shows the situation at 5 Gyr.

Panel a shows colour-coded the ratio of the escape rate predicted by the hydrodynamic model over the escape rate predicted by the energy- and recombination-limited model, Ṁ_hyd/Ṁ_enRR· The latter is the rate that is actually used to model the evolution of the planets shown in the panel. The former is merely calculated given the properties of the planets at this moment, but is not used for the evolution. This allows to compare the two models at fixed planet properties which was not easily possible in the analysis of the two individual cases.

The plot reveals the two aforementioned shortcomings of the energy- and recombination-limited approximation. The first regime is shown by the blue points: for close-in compact and massive planets with high gravitational potential exposed to high XUV irradiation, the energy- and recombination-limited model overestimates the escape rate relative to the hydrody-namic model. This is shown in panel b, which colour codes the gravitational potential for the energy- and recombination-limited case. The reason for the incorrect results of the energy-limited approximation in this particular regime has been described in Krenn et al. (2021, their Sect. 5.2): the assumed thermospheric temperatures underlying the energy-limited approximation are much higher than the ones found when directly solving the governing equations in the hydrodynamic model. The discrepancy is a consequence of the lack of downward heat conduction underlying the energy-limited approximation, leading to excessively high temperatures and thus loss rates. In reality, in the dense atmospheres at the XUV absorption height of planets with high gravitational potential, thermal conduction is a significant process, leading to lower temperatures and thus escape rates.

The second regime is shown by red and yellow points: for planets with low-to-intermediate XUV irradiation and low gravitational potential, the energy- and recombination-limited model underestimates the escape rate relative to the hydrodynamic one. Here, the energy-limited approximation fails again because of the incorrect assumed temperature structure (Krenn et al. 2021): for such planets, boil-off is the dominant escape mechanism. However, the energy-limited approximation implicitly neglects thermal energy already available in the atmosphere resulting from optical and infrared stellar irradiation. When Λ is low (less than about 30), this thermal energy is comparable to the binding energy, leading to boil-off. It is more intense for planets with lower masses, while planets more massive than approximately 10 M_⊕ are less affected (Owen & Wu 2016). The impact of boil-off is illustrated by panels c and d. In panel c, we see that the regime where the hydrodynamic model predicts significantly higher escape rates corresponds to the planets with the lowest Λ values. For them, boil-off and rapid mass loss would occur in the hydrodynamic model, but this is neglected in the energy- and recombination-limited approximation. This strong mass loss rapidly reduces the planet radius, increasing A until it approaches about 30, when boil-off stops. This explains what is seen in panel d: here, in the hydrodynamic model, there are no planets with Λ less than about 30, because the excess envelope mass has already boiled-off by 50 Myr (or even much faster, as we saw for the two individual cases studied above).

Apart from these two regimes of discrepancies, there is also a significant part where the two models yield similar escape rate (light blue - cyan - green colours in panel a). The discrepant regimes, are, however, the ones setting the valley slope.

Fig. 6 Comparison of the two escape models at 50 Myr in the orbital period-transit radius plane. In panels a, b, and c, the energy- and recombination-limited model is displayed while panel d displays the hydrodynamic model. In panel a, the colour code shows the ratio of the instantaneous escape rates predicted by the two models, Ṁhyd/ṀenRR One notes how the hydrodynamic model predicts higher escape rates for the distant small planets, but lower ones for the close-in planets. Panel b shows the gravitational potential of the planets. Panels c and d colour code the restricted Jeans parameter Λ. On the right, no Λ ≲ 30 occur, as the excess envelope has already boiled-off. Grey points are planets that have lost the entire H and He envelope.

5 Results for initial conditions derived from Kepler observations

The goal of this part is to understand if the main result obtained with the rectangular grid of initial conditions – the distinct slopes of the valley – also persist if we use more complex (and more realistic) distributions of the initial conditions, and apply a synthetic detection bias of the Kepler satellite survey to the model output. For this, we analyse the synthetic planet population obtained with the initial conditions described in Sect. 3.2, considering both the biased and the unbiased data.

The top row of Fig. 7 shows the scatter plot of transit radius as a function of orbital period for the two evaporation models at 5 Gyr. No detection bias was applied. Compared to the equivalent plot for the rectangular grid of initial conditions (Fig. 3), we see a number of similarities, but also differences. Similarities are the scarcity of close-in planets with large radii in the top left corner of the plots (the evaporation desert), and the presence of the evaporation valley running diagonally downwards. An important similarity regarding the main subject of the paper is that the slope of the valley is steeper for the energy- and recombination-limited model than for the hydrodynamic model. We quantify the slopes in the following. We also see the following differences: the period distribution is by construction different, with an increase in planet frequency from 1 to about 8 days of period, followed by a distribution constant in log P. This simply reflects the initial conditions (Eq. (9)). A more important difference concerns the presence of an overdensity of planets in the region immediately above the valley. This populates the sub-Neptunian peak in the radius histogram. At even larger radii (≳3.5 R_⊕) the frequency of planets drops strongly (the cliff, Kite et al. 2019). Both these features are important aspects of the observed planet distribution (e.g. Petigura 2020), but were absent in the rectangular grid. We see that with the inferred core and envelope mass distribution of Rogers & Owen (2021), we find these features also with our different forward (escape and interior structure) model. As a last difference we see that the valley is not fully empty, but contains some planets. There are two types of planets in the valley: First, massive bare cores (≳20 – 30 M_⊕) that started with very small post-formation envelopes (0.01 M_⊕ or less), such that they were fully evaporated despite the large core mass. These planets populate the gap from ‘below’ and are dominant in the lower half of the depleted gap area. Second, lower-mass planets that are in the process of losing the final part of their envelope. They only contain less than ~1% of their initial envelope mass. These planets populate the gap from ‘above’ and dominate in the upper half of the gap.

Our procedure to obtain the gap locus (centre) was described in Sect. 3.3. The approach employing a running mean is illustrated with Fig. 8. It shows the Kernel Density Estimate of the radius distribution for the hydrodynamic model, including all detectable planets (grey dashed line) as well as planets within three different period intervals. We see how the valley position systematically moves to smaller radii with increasing orbital distance.

From the minima in the distributions, we obtain 14 (13) positions for the middle of the valley, for the unbiased and biased population, respectively. These positions are shown with dots in the middle and bottom panels of Fig. 7. These panels show 2D Gaussian Kernel Density Estimations of the unbiased (middle) and biased (bottom) population. The impact of the detection bias which removes distant small planets is clearly visible. The two types of planets (super-Earths and sub-Neptunes) and the cliff are also clearly visible in this representation.

Finally, we have made least-square fits to these points. They are shown with white lines in the figure. For the unbiased case, we find slopes of α = −0.18 and −0.11 for the energy-and recombination-limited and the hydrodynamic model, respectively. These values are identical to those derived for the rectangular grid (Table 1). We thus find that using these very different and more realistic initial conditions does not affect the main result found with the idealised rectangular grid of initial conditions. This indicates that the imprint of the different evaporation models is quite solid, and not strongly dependent on the initial conditions, such as the assumed post-formation envelope mass.

In the biased case, which is the one most directly comparable with observations, we find slopes of −0.16 and −0.10 for the energy- and recombination-limited and the hydrodynamic model, respectively. Applying the detection bias thus makes the slopes slightly less steep for both cases, an effect that should be kept in mind when comparing (unbiased) model predictions and observations, although the difference is tiny (see also Petigura et al. 2022). More importantly, however, these values still compare in the same way to the observed values as was already found with the rectangular grid: the slope found with the hydrody-namic model is in very good agreement with the observed slope (covering a 1−σ range of about −0.13 to −0.07 depending on reference), whereas the energy- and recombination-limited model yields a too steep slope. Thus, applying a detection bias does also not alter the main conclusion of the study based on the idealised rectangular grid.

Regarding the absolute position of the valley at 10 days period, we find that the middle of the gap is predicted to be at about 2.3 R_⊕ for both evaporation models (see Table 1). As for the rectangular grid, these are larger radii than observed (1.9 ± 0.2 R_⊕ according to Van Eylen et al. 2018). Thus, our theoretical model seems to overestimate in a general way the strength of evaporation. As already discussed in Sect. 4.1, the presence of a lot of metals as coolants in the atmospheres might explain the difference. In Paper I it was found that a metal mass fraction of about Z = 0.5 would shift the valley downwards by approximately 0.4 R_⊕. This calculation did, however, employ a highly uncertain scaling of the escape rate with Z derived from photoevaporation models of protoplanetary disks (Ercolano & Clarke 2010). Systematic tabulations of atmospheric escape rates found with hydrodynamic models as the one used here but now as a function of Z (represented, e.g. as scaled solar composition) including very high values instead of pure hydrogen would help to further clarify this point. Observationally, future measurements of the atmospheric composition of sub-Neptunes with JWST will also be useful for a better understanding.

In the radius histogram including all detectable planets obtained from the model (dashed line in Fig. 8), the super-Earth peak is almost three times as high as the sub-Neptune peak. Observationally, they are in contrast of similar height (Fulton & Petigura 2018; Zhu & Dong 2021). It is not surprising that we get such a discrepancy, because the initial condition distributions we use were derived from an inference analysis utilising another evaporation (forward) model. Modifying the initial conditions would, however, allow to change this ratio: by shifting the core mass distribution to more massive values, a higher fraction of planets would be massive enough to keep a H and He envelope and populate the sub-Neptunian peak. The minimum core mass necessary to keep some H and He at a given orbital distance was analytically derived in Paper I (their Eq. (29)).

Fig. 7 Transit radius as a function of orbital period at 5 Gyr for distributions of the initial envelope mass, core mass, and orbital period derived from Kepler observations (Rogers & Owen 2021). Left column: energy- and radiation-recombination-limited escape model. Right column: hydrodynamic escape model. Top row: raw scatter plot of the unbiased synthetic populations. Middle row: 2D Gaussian Kernel Density Estimation of the unbiased synthetic populations. Bottom row: as in the middle, but after applying a detection bias representative of the Kepler survey, which disfavours small distant planets. White dots and lines indicate the valley position. One notes in all cases the shallower slope in the hydrodynamic model.

6 Summary and conclusions

For this work, we tested both a simpler XUV-driven energy-and recombination-limited escape model and a complex hydro-dynamic escape model (Kubyshkina et al. 2018) against a key observational constraint, the valley slope. The latter model includes the boil-off, blow-off, and Jeans escape regimes. The comparison was done by coupling the two escape models to a model for the temporal evolution of planetary interiors. This interior model solves the classical spherically symmetric interior structure equations. With these models, we simulated the evolution of 6000 planets on an idealised rectangular grid in orbital period and mass, and for about 37 000 planets with initial conditions (period, core, and envelope mass) derived from an inference analysis of the Kepler survey planet population (Rogers & Owen 2021). We studied the valley locus predicted by the two escape models at 5 Gyr.

We find that the hydrodynamic model leads to a valley slope d log R/d log p = −0.11 both for the rectangular grid and the unbiased synthetic Kepler planet population. For the hydrody-namic model, applying a simple detection bias of the Kepler survey (Petigura et al. 2018) leads to a slope of −0.10. These slopes agree closely with the observational result derived by Van Eylen et al. (2018, $-{0.09}_{-0.04}^{+0.02}$) and Martinez et al. (2019) and Petigura et al. (2022, −0.11 ± 0.02). As past photoevap-oration models, the simple energy- and recombination-limited escape model in contrast predicts a slope that is too steep of −0.18 for the rectangular grid and the unbiased synthetic Kepler population, and of −0.16 for the biased synthetic population. Regarding the radius of the lower boundary of the valley at a fixed 10-day orbital period, both models yield similar values for the rectangular grid, namely about 1.8–1.9 R_⊕.

The slope that is too steep in the energy- and recombination-limited escape model is caused by two limitations of this approximation (Krenn et al. 2021), as is found by comparing the escape rates for both individual planets and the entire grid: In particular, it underestimates escape rates for distant, fluffy low-mass planets while simultaneously overestimating it for close-in, compact massive planets. The former is caused by the omission of the boil-off escape regime in the purely XUV-driven energy- and recombination-limited model, while the latter can be explained by its negligence of heat conduction in the atmosphere.

Boil-off (Stökl et al. 2015; Owen & Wu 2016; Fossati et al. 2017) causes a rapid mass decrease in the first few megayears for fluffy planets with a low restricted Jeans-escape parameter. These are planets with considerable thermal energy stored in their atmosphere relative to their gravitational potential. This initial mass loss is significant enough to alter the slope of the valley by evaporating the atmosphere of more massive planets at larger distances. It is interesting to note that the escape rate in the boil-off regime depends on the stellar continuum irradiation (VIS/IR) via the planetary equilibrium temperature, and not on the XUV irradiation. This is a property it shares with core-driven escape, contrasting the purely XUV-driven energy- and recombination-limited model. Our results suggest that a combination of aspects of both models (namely heating both in VIS/IR and XUV) yield a valley slope in agreement with observations.

The second limitation, the negligence of heat conduction in the energy- and recombination-limited approximation produces higher temperatures in the atmosphere than when conduction is cooling the upper atmosphere, as it is the case in the hydro-dynamic model. The energy- and recombination-limited model therefore overestimates the temperature, leading to an excessive mass-loss rate. This effect is prevalent for massive close-in planets, which are highly XUV irradiated and feature compact atmospheres (Krenn et al. 2021). By including conduction-cooling, the hydrodynamic model predicts lower mass-loss rates over time for such planets, leaving lower-mass planets still with a H and He envelope, which also alters the valley slope. In combination, the two shortcomings act together in the same direction: the evaporation that is too weak at larger distances (resulting in smaller evaporated cores) and evaporation that is too strong at smaller distances (resulting in larger evaporated cores) give the valley a slope that is too steep in the energy- and recombination-limited model.

Our results indicate that the more realistic description of the thermospheric temperature structure in the hydrodynamic model relative to the energy- and recombination-limited model is critical. It allows one of the most important observational constrains for escape models to be reproduced, the valley slope.

Future work will address the evaporation valley’s dependency on host star mass (see also Gupta et al. 2022) and the effect of including metals, which may act as coolants. When compared with observational studies exploring the temporal (Berger et al. 2020a; Sandoval et al. 2021; David et al. 2021; Petigura et al. 2022; Ho & Van Eylen 2023) and stellar mass dependency (Fulton & Petigura 2018; Wu 2019; Berger et al. 2020b; Petigura et al. 2022; Chen et al. 2022), this should allow one to develop an even better understanding of the origin of the radius valley.

Fig. 8 Kernel density estimate of the distribution of the radii in the biased synthetic population obtained with the hydrodynamic escape model. The grey dashed line includes all detectable planet (at all orbital periods). The green, orange, and blue lines include planets with log (period/day) of 0.8 ± 0.1, 1.3 ± 0.1, and 1.8 ± 0.1. One sees how the centre of the valley shifts to smaller radii with increasing orbital period.

Acknowledgements

L.A., A.V.O., and C.M. acknowledge support from the Swiss National Science Foundation under grant 200021_204847 ‘PlanetsInTime’. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. Calculations were performed on the Horus cluster at the University of Bern. The research described in this paper was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics Space Administration.

Appendix A Impact of the post-formation luminosity

In this appendix, we evaluate the sensitivity of our results for the valley locus on the post-formation luminosity L₀ or a closely related quantity, the specific entropy s₀ in the inner convective zone.

A.1 Parameterisation and expected range of s₀

As in Paper I, our nominal approach is to assume an initial luminosity L₀ that is given by a fit to results of planet formation population syntheses (Mordasini et al. 2014),

In this equation, is the intrinsic luminosity of Jupiter today (about 8.7 × 10⁻¹⁰L_⊙). It is clear that in reality, depending on the particular formation history of a planet, the post-formation luminosity may vary (e.g. Marley et al. 2007; Mordasini et al. 2017; Cumming et al. 2018; Marleau et al. 2019). Mordasini et al. (2017) for example found a spread in post-formation entropy s₀ in the planet mass range of interest here of about 1 to 1.5 k_B /baryon at fixed envelope mass.

These earlier works investigated the post-formation entropy mainly in the context of giant planets and their detectability with direct imaging. More recently, the impact of the post-formation thermodynamic state was also addressed for evaporating planets: Owen (2020) showed how the post-formation entropies of young evaporating planets might be constrained by observations. Kubyshkina & Vidotto (2021) found that the initial entropies of planets have a minor or even absent effect on most of the population of evolved planets with ages of ~ 1 Gyr. Only for planets suffering extremely strong atmospheric mass loss, s₀ was found to be of importance. A low importance of the entropy is also expected from the rather weak dependency of the thickness of the convective zone of H and He envelopes on the age and thus entropy (Lopez & Fortney 2014).

Our approach here is to re-run the evolutionary simulations on the rectangular grid of initial conditions with the two evaporation models, but instead of using Eq. A.1, we cover a wide range of s₀, including the one suggested by more modern formation models. We can then systematically study how s₀ affects the valley location. For this, we generalise the parameterisation of s₀ of Malsky & Rogers (2020), $$s_0=s_{0,\text{n}}+\frac{M_{\text{p}}}{25M_{\oplus }}{k}_{\text{B}}/\text{baryon.}$$(A.2)

Malsky & Rogers (2020) fixed s_0,n to 6 k_B /baryon. Here, we generalise this and use s_0,n = 6, 7, 8, and 9 k_B/baryon.

Before discussing the results of these grid simulations with different s_0,n, we compare the s₀ obtained in this way with the ones predicted by the recent comprehensive planet population synthesis simulation NGPPS (Emsenhuber et al. 2021b). These simulations represent a much improved update to the ones used to derive the original fitting equation (Mordasini et al. 2014). These NGPPS results are generated with the Generation III Bern Model (Emsenhuber et al. 2021a) which is a complex end-to-end formation and evolution model based on the core accretion paradigm. The model solves the 1D internal structure equations during the formation (both attached and detached state), and evolutionary phases. In the luminosity calculation it takes into account the accretion of planetesimals and gas, the cooling and contraction of the envelope, radiogenic heating, as well as giant impacts. The model also takes the concurrent formation of several planets in one disk into account, in contrast to the older Mordasini et al. (2014) syntheses, which used the one-embryo-per-disk approximation. This leads to more diverse formation pathways (Emsenhuber et al. 2021b).

The left panel of Fig. A.1 shows as black dots the entropy at the core-envelope boundary of the planets the NGPPS simulation. The nominal synthetic population NG76 is shown at the moment when the gas disk dissipates, which corresponds to ages between 1 and 10 Myr. The host star mass is 1 M_⊕. We see that generally speaking, the (mean) entropy is an increasing function of the planet mass. Especially at smaller masses, there is significant spread in s₀. Given the high density of points in this mass range, it is however difficult to get a quantitative picture of the distribution of s₀ from the scatter plot. Thus, in the right panel we additionally show the kernel density estimation of the distribution of s₀ for three mass ranges of interest for our study. We see that the mode indeed increases with mass, lying at about 7.3, 8.0, and 8.3 k_B/baryon for the low-, mid-, and high-mass range. The FWHM is about 1 to 1.5 k_B/baryon in the three cases. This spread is thus similar to the one in the older syntheses.

The grey points show the s₀ obtained with the nominal approach (Eq. A.1). Since we here specify L₀ and not directly the entropy, and given that the atmospheric boundary conditions (in particular the temperature) also affect the relation between L₀ and s₀ (see Marleau & Cumming 2014 and Kubyshkina & Vidotto 2021), a range of s₀ occur. The points fall on lines of fixed orbital distance (or equilibrium temperature), with the high s₀ values corresponding to the closest distances. We also see that the majority of the grey points falls into a similar range as also the majority of the black points do. Thus, the simple fit derived from the older population synthesis still seem to capture — at least in a rough way — the new NGPPS results for s₀. The four coloured lines finally show Eq. A.2 with the four values of s_0,n. One sees that the two extreme values (s_0,n = 6 and 9 k_B /baryon) are not representative of the predictions of the formation simulations, but are clearly too low respectively too high in comparison (one should here keep in mind the logarithmic nature of the entropy. It means that a numerically small difference in s_0,n actually corresponds to a very significant change of the gravothermal heat content). As visible from Fig. A.1, in the formation simulations there are, in particular, virtually no low- and intermediate-mass planets with s₀ as low (high) as 6 (9) k_B /baryon.

A.2 An example case

Figure A.2 shows the temporal evolution of a specific planet from the rectangular grid for the four s_0,n. The energy- and recombination-limited escape model is used, but qualitatively equivalent effects are also occurring for the hydrodynamic model. This planet has an orbital distance of 0.1 AU, M_core = 20M_⊕, and M_env,0 = 3.64M_⊕. The L₀ is 0.34, 12.1, 602.4, and for the four s₀ studied. The latter value is certainly extremely high for a planet of this mass (Mordasini et al. 2017). Equation A.1 yield for comparison .

This planet was chosen because it illustrates with a specific example the two main findings of the grid analysis of the valley location as a function of s_0,n in the next section: namely a weak impact of s_0,n for the three lower entropy values, and envelope overflow for some planets for s_0,n = 9 k_B/baryon.

In the left panel we see the radius as a function of time. The initial radius is as expected the larger the higher s_0,n is. This has the well-known consequence (e.g. Owen 2020) that stronger XUV-driven atmospheric escape occurs for a higher s₀ at young ages, such that at high ages, when the initial s if forgotten, the planet has a smaller radius because it has a lower-mass envelope. For the highest s₀ case, there is, however, an additional effect: the huge initial radius is here larger than R_roche, meaning that some envelope gas is unbound. In the model, we then remove at each time-step one third of the mass outside of R_roche. This factor smaller than unity (which would in principle be the value to use) was chosen for numerical stability. The exact value is, however, inconsequential: in any case, extremely rapid and strong mass loss occurs until the outer radius becomes smaller than R_roche, and only the time duration until this happens varies somewhat with the specific fraction chosen. On the other hand, in a situation of such rapid mass loss as here, quantitative results of our 1D strictly hydrostatic model with a radially constant luminosity at a given time should be regarded with caution.

Fig. A.1 Post-formation entropy distributions. Left panel: Post-formation entropy s0 as a function of planet mass. The black dots show s0 predicted by core accretion planet formation simulations in the New Generation Planet Population Synthesis NGPPS (Emsenhuber et al. 2021b). Grey dots are the s0 in the nominal case (Eq. A.1). The coloured lines finally represent the four parameterisations used in this appendix, which are a generalisation of Malsky & Rogers (2020). Right panel: Kernel density estimation of the distributions of s0 in the NGPPS of three intervals of the total planet mass. One sees how the mode of the distribution shifts to higher values with increasing planet mass.

Fig. A.2 Temporal evolution of the outer radius (left), total mass (middle), and Kelvin-Helmholtz timescale (right panel) for a 23.64 M⊕ planet (Mcore = 20M⊕) at 0.1 AU for the four s0,n indicated in the plots. Time is measured relative to the moment when the simulation starts. In the left panel, the grey line shows the Hill sphere radius. The planet with the (unrealistically) high s0,n = 9 kB/baryon initially overflows the Hill sphere, leading to a strong reduction of the envelope mass. At late times, this leads to a smaller radius. The other three cases which lack this overflow phase given in contrast similar values, with a slight anti-correlation of the radius at late time and the initial entropy.

As is visible in the middle panel, this overflow phase removes about one third of M_env,0 on an extremely short timescale which is of the order of just 100 years. At late times, this Roche lobe overflow has the consequence that the planet has a clearly smaller radius (4.05 R_⊕) and mass than the other three cases. For them, the radius varies only between 4.39 and 4.47 R_⊕. The largest radius corresponds to the lowest s₀ because this planet suffered from less ‘normal’ XUV-driven escape because of its higher mean density at young ages, as mentioned.

The occurrence of such an overflow phase is indicative of an unrealistic combination of initial conditions for the evolutionary phase in terms of core mass, envelope mass, and luminosity. In reality, during the precedent formation phase, while embedded in the nebula, a core of such a high luminosity (caused for example by a burst of solid accretion) would not posses an envelope of this mass. Instead, potential excess gas would get expelled out of the Hill sphere back into the surrounding disk, and M_env would be lower than assumed here. This effect is by construction not taken into account when s₀ is assumed to be only a function of the total mass, as it is the case both for Eq. A.1 and A.2. In the formation simulations solving the internal structure equations, this is in contrast automatically taken into account. Thus, whenever possible, s₀ should be estimated in evolutionary models not only based on the total planet mass, but the core and envelope mass separately. Such a prescription for L₀(M_core, M_env) can be found in Paper I.

Fig. A.3 Transit radius as a function of orbital period at 5 Gyr for the four initial entropies s0,n = 6, 7, 8, and 9 kB /baryon. The lower three values lead to virtually identical valley slopes. For the highest initial entropy, unstable initial conditions with Roche lobe overflow result, which strongly removes mass especially at the larger orbital distances. This affects the valley slope. Such a high entropy is, however, not in agreement with the predictions of formation models.

The right panel of Fig. A.2 shows the Kelvin-Helmholtz timescale. It is calculated with the actual numerically obtained total energy of the planets and not the approximation $G{M}_p^2/R_p$. This approximation can yield very different incorrect values for planets with a very large, extremely tenuous outer envelope, as it is the case here at early times. We see that s_0,n = 9k_B/baryon corresponds to an extremely small T_KH of less than 10⁴ yr. A planet would thus extremely quickly evolve away from such conditions, making it an unlikely state to exist exactly at the moment of disk dissipation. The lowest s_0,n = 6 k_B/baryon on the other hand yields an extremely long initial T_KH ~ 10¹⁰ yr. The radius hardly changes for about 1 Gigayear. Such an extremely cold start seems also unlikely given the energy liberated when accreting a solid core of 20 M_⊕.

A.3 Valley locus as function of s₀

The example of this individual planet suggests that the impact of s₀ should be rather small, except if an unexpected high entropy is used. Figure A.3 showing the rectangular grid of simulations indeed reflects a similar pattern. The figure, which can be compared with Fig. 3 using the nominal L₀, shows for the four s_0,n the radius at 5 Gyr as a function of orbital period, colour-coding the total mass. The hydrodynamic escape model is used. As for the nominal rectangular grid (Sect. 4), we have made a least square fit to determine the valley slope a and the normalisation radius at 10 days period, ${\tilde{R}}_b$. The values for the hydrodynamic model are given in Table A.1. The result for the energy- and recombination-limited model are similar.

The panel in the bottom left corner of Fig. A.3 differs clearly from the other three, which are in contrast similar to each other. We see an absence of planets in the upper right corner. The iso-mass curves visible through the colour-code are significantly shifted, especially at larger orbital distances. The valley slope is also less steep than in other three cases. These differences are the consequence of intense mass loss resulting from Roche lobe overflow right at the beginning of the simulations, as seen in the example case. It affects both planets far from the valley (as in the example), but also planets close to it. An interesting point is here that the maximum radii are limited to about 3.5 to 4 R_⊕ which corresponds to the radius above which observationally the frequency of planets drop strongly (the cliff, Kite et al. 2019). This suggest that in the very high entropy case, it is not necessary to fine-tune the initial (i.e. post-formation) H and He masses to reproduce the cliff. This echoes the suggestion of Owen & Wu (2016) that the ‘boil-off’ process could be partially responsible for the lack of larger planets.

Another small feature visible in Fig. A.3 is the absence of points in the bottom left corner. This is a consequence of the following: for these very close-in, low-mass planets, no structure was found for the requested s0. Because of Eq. 8, these planets have tenuous atmospheres approaching for lower s₀ an isothermal structure. The given equilibrium temperature then excludes certain combinations of core mass, envelope mass, and s₀. This is in contrast to colder and more massive planets with an inner convective zone.

As discussed in the previous section, our quantitative results for planets undergoing Roche lobe overflow should be taken with caution, given our model’s capabilities. However, this process only affects planets with s_0,n = 9 k_B/baryon, which is not a likely initial condition for low-mass planets. The important conclusion from examining the role of s₀ is rather the following: for the relevant, very wide range of s_0,n from 6 to 8 k_B /baryon, the impact of the post-formation entropy on the final valley slope is only very small, as can be seen in Table A.1. The slope a hardly changes with values between −0.113 and −0.117. This is also the same as found for the nominal L₀. We do see that ${\tilde{R}}_b$ shifts to higher values with increasing s₀ as expected, but the shift from 1.81 to 1.90 R_⊕ is rather small. This is especially the case when one considers that formation models predict a spread of about only 1, and not 3 k_B /baryon at fixed planet mass.

To summarise, we find that varying s_0,n over a wide range of 6 to 8 kB/baryon has virtually no impact on the valley slope, and shifts ${\tilde{R}}_b$ only by a rather small increment. This range of initial entropies includes those suggested by formation models and leads to stable initial conditions where the initial planet radius is smaller than the Roche lobe. Only when using a for this mass range unrealistically high s_0,n = 9 k_B/baryon, the impact becomes significant, because mass loss via Roche lobe overflow occurs immediately at the beginning of the simulations. Such unstable initial conditions are, however, not predicted by formation models.

References

All Tables

All Figures

Fig. 1 Temporal evolution of the bolometric (blue), X-ray (green), and EUV luminosity (brown) of a 1 M⊙ star as assumed in our model. The bolometric luminosity is divided by a factor 1000 to bring it on a similar scale as LX and LEUV. The two black bars near 4.5 Gyr show the range of our Sun’s LX over the course of a solar cycle. The grey dashed lines show for comparison the LX of Tu et al. (2015) for the 10th, 50th, and 90th percentiles of the rotational distribution.

In the text

Fig. 2 Initial conditions derived from Kepler survey. Left and middle panel: histogram of the distribution of the core masses and initial (post-formation) envelope mass fractions Menv,0/Mcore for the comparison with the Kepler planet population. Both these distributions are taken from Rogers & Owen (2021). Right panel: average detection probability of the Kepler satellite pdet as function of orbital period and planet radius (Petigura et al. 2018). The final probability to find a planet is given by pdet times its geometric transit probability ptr.

In the text

Fig. 3 Results of planet evolution simulations in the orbital period-transit radius plane at 5 Gyr (left panel: energy- and recombination-limited model; right panel: hydrodynamic model). Each dot represents a planet which is coloured based on the fraction of its initial H and He envelope mass that was evaporated. Grey symbols indicate a complete loss of the envelope, corresponding to sub-Neptunes that have evolved into super-Earths. The magenta lines show the fit to determine the slope α of the valley as represented by the largest super-Earth at a given period (i.e. the bottom of the valley). We note the shallower slope in the hydrodynamic model. The squares show two individual cases (one far planet in lilac, one close in black). Both start with identical initial conditions in the two models and are discussed in Sect. 4.2. They end up on different sides of the valley in the two models in a way to make the slope shallower in the hydrodynamic model.

In the text

Fig. 4 Temporal evolution of the evaporation rate (top left), restricted escape parameter Λ (top right), transit radius (bottom left), and total mass (bottom right) for a planet at 0.51 AU and initial mass of approximately 3.86 M⊕. The hydrodynamic and the energy- and recombination-limited model are shown. We note the very high mass loss in the hydrodynamic model in the initial boil-off phase. It is sufficient to eventually cause the total loss of the envelope. In the energy- and recombination-limited model which lacks this phase, the planet can in contrast keep a significant fraction of its envelope. The planet therefore resides in the end on different sides of the valley in the two models (super-Earth for the hydrodynamic model, sub-Neptune for the energy- and recombination-limited model). This is shown by the lilac squares in Fig. 3.

In the text

Fig. 5 Temporal evolution of the evaporation rate (top left), restricted escape parameter Λ (top right), transit radius (bottom left), and total mass (bottom right) for a planet at 0.04 AU and initial mass of approximately 20.9 M⊕. The hydrodynamic and the energy- and recombination-limited model are shown. Despite the lack of the initial boil off phase in the energy- and recombination-limited model, the higher mass-loss rate in this model after about 30 Myr results in complete envelope loss, in contrast to the hydrodynamic model. The planets therefore reside in the end on different sides of the valley (super-Earth for the energy- and recombination-limited model, sub-Neptune for the hydrodynamic model). The underlying reason is the negligence of conduction in the energy- and recombination-limited model, leading to too high escape rates. This planet is shown by the black squares in Fig. 3.

In the text

Fig. 6 Comparison of the two escape models at 50 Myr in the orbital period-transit radius plane. In panels a, b, and c, the energy- and recombination-limited model is displayed while panel d displays the hydrodynamic model. In panel a, the colour code shows the ratio of the instantaneous escape rates predicted by the two models, Ṁhyd/ṀenRR One notes how the hydrodynamic model predicts higher escape rates for the distant small planets, but lower ones for the close-in planets. Panel b shows the gravitational potential of the planets. Panels c and d colour code the restricted Jeans parameter Λ. On the right, no Λ ≲ 30 occur, as the excess envelope has already boiled-off. Grey points are planets that have lost the entire H and He envelope.

In the text

Fig. 7 Transit radius as a function of orbital period at 5 Gyr for distributions of the initial envelope mass, core mass, and orbital period derived from Kepler observations (Rogers & Owen 2021). Left column: energy- and radiation-recombination-limited escape model. Right column: hydrodynamic escape model. Top row: raw scatter plot of the unbiased synthetic populations. Middle row: 2D Gaussian Kernel Density Estimation of the unbiased synthetic populations. Bottom row: as in the middle, but after applying a detection bias representative of the Kepler survey, which disfavours small distant planets. White dots and lines indicate the valley position. One notes in all cases the shallower slope in the hydrodynamic model.

In the text

Fig. 8 Kernel density estimate of the distribution of the radii in the biased synthetic population obtained with the hydrodynamic escape model. The grey dashed line includes all detectable planet (at all orbital periods). The green, orange, and blue lines include planets with log (period/day) of 0.8 ± 0.1, 1.3 ± 0.1, and 1.8 ± 0.1. One sees how the centre of the valley shifts to smaller radii with increasing orbital period.

In the text

Fig. A.1 Post-formation entropy distributions. Left panel: Post-formation entropy s0 as a function of planet mass. The black dots show s0 predicted by core accretion planet formation simulations in the New Generation Planet Population Synthesis NGPPS (Emsenhuber et al. 2021b). Grey dots are the s0 in the nominal case (Eq. A.1). The coloured lines finally represent the four parameterisations used in this appendix, which are a generalisation of Malsky & Rogers (2020). Right panel: Kernel density estimation of the distributions of s0 in the NGPPS of three intervals of the total planet mass. One sees how the mode of the distribution shifts to higher values with increasing planet mass.

In the text

Fig. A.2 Temporal evolution of the outer radius (left), total mass (middle), and Kelvin-Helmholtz timescale (right panel) for a 23.64 M⊕ planet (Mcore = 20M⊕) at 0.1 AU for the four s0,n indicated in the plots. Time is measured relative to the moment when the simulation starts. In the left panel, the grey line shows the Hill sphere radius. The planet with the (unrealistically) high s0,n = 9 kB/baryon initially overflows the Hill sphere, leading to a strong reduction of the envelope mass. At late times, this leads to a smaller radius. The other three cases which lack this overflow phase given in contrast similar values, with a slight anti-correlation of the radius at late time and the initial entropy.

In the text

Fig. A.3 Transit radius as a function of orbital period at 5 Gyr for the four initial entropies s0,n = 6, 7, 8, and 9 kB /baryon. The lower three values lead to virtually identical valley slopes. For the highest initial entropy, unstable initial conditions with Roche lobe overflow result, which strongly removes mass especially at the larger orbital distances. This affects the valley slope. Such a high entropy is, however, not in agreement with the predictions of formation models.

In the text