© ESO 2021

1 Introduction

It is believed that planets accrete an extensive hydrogen atmosphere when they form inside the protoplanetary nebula (e.g. Stökl et al. 2016; Lammer et al. 2018; Owen et al. 2020). Following the dispersion of the nebula, the hydrogen atmosphere can be subject to significant mass loss. In the beginning, the escape is driven by a combination of high temperatures and low gravity, leading to a phase of exhaustive hydrodynamic escape (boil-off or core-powered mass loss; Stoekl et al. 2015; Owen & Wu 2016; Ginzburg et al. 2016; Fossati et al. 2017; Kubyshkina et al. 2018a; Gupta & Schlichting 2019, 2020). The high temperatures in this phase of evolution are mostly due to energy released from the contracting planetary interior and, for the most close-in planets, to the high stellar UV, optical, and infrared irradiation. The low gravity, however, is due to the large planetary radius that, at the time of the dispersal of the protoplanetary nebula, reaches as far as the Hill radius (e.g. Stoekl et al. 2015; Owen & Wu 2016; Ginzburg et al. 2016).

Following this initial phase of extreme mass loss, atmospheric escape is no longer dominated by the thermal energy and low gravity of the planet, but rather by upper atmospheric heating caused by absorption of the stellar high-energy radiation, that is, extreme X-ray and ultraviolet (XUV) radiation (Hunten & Watson 1982; Lammer et al. 2003; Yelle 2004). As the XUV luminosity of the host star declines with time following the evolution of (i.e. decrease in) the stellar rotation rate (Ribas et al. 2005; Mamajek & Hillenbrand 2008; Wright et al. 2011; Tu et al. 2015; Johnstone et al. 2021), the XUV-driven hydrodynamic escape weakens until other escape processes, such as Jeans escape or non-thermal escape processes, start to dominate mass loss.

The discovery of a large variety of planets, particularly those with radii in between those of Earth and Neptune (e.g. Mullally et al. 2015), and the detection of structures in the observed exoplanet population, such as the radius gap separating super-Earths and mini-Neptunes at roughly 1.6 Earth radii (Fulton et al. 2017), demonstrate the importance of studying and understanding planetary atmospheric evolution. Such studies gain importance considering that the Solar System planets, including the Earth, have most likely gone through substantial atmospheric evolution as well (e.g. Lammer et al. 2018, 2020). Within this context, it has been shown that planetary atmospheric escape is not only a fundamental process shaping the observed exoplanet population, but also plays a key role in determining whether a planet develops a habitable environment (e.g. Jin et al. 2014; Owen & Wu 2017; Jin & Mordasini 2018; Owen & Lai 2018; Lammer et al. 2018; Ginzburg et al. 2018; Gupta & Schlichting 2019, 2020).

Different hydrodynamic models, such as those of Yelle (2004), García Muñoz (2007), Murray-Clay et al. (2009), Koskinen et al. (2013), Salz et al. (2015), Erkaev et al. (2016), Guo & Ben-Jaffel (2016), and Kubyshkina et al. (2018a), have been used to study upper atmospheres, thermal mass loss, and how escape varies as a function of system parameters for planets hosting hydrogen-dominated atmospheres. However, these hydrodynamic models are rather complex, which makes them not particularly suitable for atmospheric evolution simulations. To overcome this, analytic approximations, in particular the energy-limited escape approximation (Watson et al. 1981; Erkaev et al. 2007), are often used.

Recently, Kubyshkina et al. (2018a,b) presented analytic functions for estimating mass-loss rates based on the outcome of hydrodynamic simulations that had already been successfully used for planetary evolution calculations (e.g. Kubyshkina et al. 2019a,b, 2020; Modirrousta-Galian et al. 2020a,b). Nevertheless, the energy-limited approximation is still widely used by the community for estimating mass-loss rates for planets over a very broad parameter space (e.g. Locci et al. 2019; Poppenhaeger et al. 2021). Therefore, we investigate here in detail the applicability of the energy-limited escape approach across a wide parameter space.

Other authors have previously studied this applicability by means of comparisons with hydrodynamic simulations. Kubyshkina et al. (2018b) compared the results obtained from a large sample of hydrodynamic models (the same used in this study) with the energy-limited mass-loss rates. They reported that the energy-limited approximation significantly underestimates mass-loss rates for planets that are very close to boil-off conditions and moderately overestimates mass-loss rates for planets with (almost) hydrostatic atmospheres. However, in contrast to our present study, Kubyshkina et al. (2018b) employed an XUV absorption radius inferred from the simulations instead of that derived from the energy-limited formalism. Murray-Clay et al. (2009) followed an approach similar to that of Kubyshkina et al. (2018b) but focused on hot Jupiters and employed a less extensive set of hydrodynamic models. They observed that the energy-limited approximation leads to a significant overestimation of the mass-loss rates, by up to two orders of magnitude for highly irradiated hot Jupiters. They concluded that this difference is due to the fact that the energy-limited approximation wrongly estimates the level of unit optical depth to XUV radiation and that conductive cooling is irrelevant in the upper atmospheres of hot gas giants. Erkaev et al. (2013) found that the energy-limited mass-loss rates significantly overestimate the escape, particularly for planets with low irradiation levels, as they fail to meet the conditions of a hydrodynamically expanding atmosphere. Finally, Tian et al. (2005) suggested that the discrepancies between the energy-limited and the hydrodynamically computed mass-loss rates lie in the single-layer absorption assumption of the energy-limited approximation, which in many cases is not sufficiently accurate.

We aim to provide an assessment of the applicability of the energy-limited approximation over a wide parameter space, highlighting the regions characterised by large discrepancies between energy-limited and hydrodynamically computed mass-loss rates. We additionally attempt to provide reasons for the failure of the energy-limited approximation. Furthermore, in contrast to all previous studies, we only use XUV absorption radii predicted within the energy-limited framework.

This paper is organised as follows. We review the energy-limited formalism and the involved assumptions in Sect. 2. In Sect. 3 we present a comparison of the results obtained employing the simplest form of the energy-limited approximation with those obtained by solving the energy-limited equations. Section 4 presents a comparison of the energy-limited mass-loss rates with those obtained from hydrodynamic simulations. In Sect. 5 we review the applicability of the energy-limited approximation across the considered parameter space and draw our conclusions in Sect. 6.

2 Review of the energy-limited atmospheric escape approach

Watson et al. (1981) introduced the concept of energy-limited atmospheric escape in the context of hydrogen atmospheres. The aim was to provide an upper limit for the mass-loss rate by equating the escape flux and the incident XUV energy. Erkaev et al. (2007) then extended the formulation of Watson et al. (1981) by considering tidal effects, obtaining $$\dot{M}=\pi \frac{\nu {\Phi }_{\text{XUV}}{r}_{\text{XUV}}^2{r}_0}{KG{M}_{\text{pl}}},$$(1)

where Ṁ is the mass-loss rate, ν the efficiency of XUV heating (see e.g. Shematovich et al. 2014; Salz et al. 2016), Φ_XUV the incident stellar XUV flux, r_XUV the radius at which the XUV flux is absorbed in the atmosphere, r₀ an arbitrary lower boundary radius usually defined as the photospheric radius of the planet (R_pl), G the gravitational constant, and M_pl the planetary mass. Here, K is the correction factor, ranging between zero and one, which was introduced by Erkaev et al. (2007) to account for the tidal effects that were first discussed by Gu et al. (2003). The existence of the Roche lobe means that to escape the gravitational potential of a planetary body, atmospheric particles must only reach the Roche lobe, which in turn means that the energy each individual particle needs to escape is lower than if no host star were present. The correction factor K, which is particularly important for close-in low-mass planets, is defined as (Erkaev et al. 2007) $$K=1-\frac{3R_{\text{pl}}}{R_{\text{Rl}}}+\frac{R_{\text{pl}}^3}{2R_{\text{Rl}}^3},$$(2)

where the Roche radius (R_Rl) is $$R_{\text{Rl}}\approx d_0{\left(\frac{M_{\text{pl}}}{3M_{\text{star}}}\right)}^{\frac{1}{3}}.$$(3)

In Eq. (3), d₀ is the planetary orbital separation and M_star is the stellar mass.

The term r_XUV is usually unknown and, particularly in extended hydrogen envelopes, is significantly larger than r₀ (Watson et al. 1981; Kubyshkina et al. 2018a). Employing a variety of assumptions and simplifications, detailed below, Watson et al. (1981) derived a system of equations that allows the simultaneous calculation of Ṁ and r_XUV with numerical methods. Below we summarise the framework of Watson et al. (1981).

We remark that throughout this work we consider the atmospheric mass-loss rate (Ṁ) in kg s⁻¹, while Watson et al. (1981) primarily considered the particle escape rate (F) in particles sr⁻¹ s⁻¹. The relationbetween F and Ṁ is $$\dot{M}=4\pi mF.$$(4)

2.1 Model description and assumptions

Watson et al. (1981) defined dimensionless variables for atmospheric height (λ), temperature (τ), velocity (Ψ), particle escape (ζ), and the energy parameter (β), which, accounting for tidal effects, are defined as $$\lambda =\frac{KG{M}_{\text{pl}}m}{k{T}_{0}r},$$(5) $$\tau =\frac{T}{T_0},$$(6) $$\Psi =\frac{m{u}^2}{k{T}_0},$$(7) $$\zeta =\dot{M}\frac{k^2{T}_0}{4\pi {\kappa }_{0}KG{M}_{\text{pl}}{m}^2},$$(8)

and $$\beta =\nu {\Phi }_{\text{XUV}}\frac{KG{M}_{\text{pl}}m}{k{T}_0^2{\kappa }_0},$$(9)

respectively. In Eqs. (5)–(9), k is the Boltzmann constant, m the mass of an escaping particle, T₀ the temperature at r₀, κ₀ the thermal conductivity at r₀, r an arbitrary radius at or above r₀, T the temperature at r, and u the bulk velocity at r. We note that λ is the Jeans escape parameter and that, when using the photospheric radius R_pl as r and the equilibrium temperature T_eq as T₀, Eq. (5) reduces to the ‘restricted Jeans escape parameter’ (Fossati et al. 2017): $$\Lambda =\frac{KG{M}_{\text{pl}}m}{k{T}_{\text{eq}}{R}_{\text{pl}}}.$$(10)

Watson et al. (1981) parametrised the thermal conductivity of a neutral gas as $$\kappa ={\kappa }_0{\tau }^s,$$(11)

where s = 0.7 is usually chosen for a hydrogen-dominated gas. Using Eqs. (5)–(7), the steady state equations for mass, momentum, and energy conservation (see Eqs. (3)–(5) in Watson et al. 1981) can be reduced to the following two equations $$\left(1-\frac{\tau }{\Psi }\right)\frac{\text{d}\Psi }{\text{d}\lambda }=2\left(1-\frac{2\tau }{\lambda }-\frac{\text{d}\tau }{\text{d}\lambda }\right)$$(12)

and $$\frac{{\tau }^s}{\zeta }\frac{\text{d}\tau }{\text{d}\lambda }=ϵ+\lambda -\frac{5\tau }{2}-\frac{\Psi }{2},$$(13)

where ϵ is the energy each escaping particle carries away from the system and is defined as $$ϵ={ϵ}_{\infty }-\frac{k}{\zeta {\kappa }_{0}KG{M}_{\text{pl}}m}{\displaystyle {\int }_r^{\infty }q}{r}^2\text{d}r.$$(14)

In Eq. (14), q is the volume heating rate and ϵ_∞ denotes the energy flow at infinity. In the special case where ϵ_∞ = 0, the energy flowing outwards is just enough to lift the gas from the gravitational field but does not leave any excess energy at very large distances. Watson et al. (1981) assumed this to be the case for every escaping particle since any additional energy carried away would only decrease the maximum particle escape flux. Therefore, ϵ_∞ is set equal to zero.

Three assumptions need to be made to define ϵ: (1) The incident XUV energy is absorbed in a single narrow atmospheric layer; (2) the fraction of the absorbed radiation spent on escape is equal to the heating efficiency (ν); and (3) the incident XUV energy is averaged over a sphere by a factor of 1/4 as this is the ratio of the cross-section and surface area of a sphere (Erkaev et al. 2007). With these assumptions, Eq. (14) becomes $$ϵ=-\frac{\pi m\nu {\Phi }_{\text{XUV}}{r}_{\text{XUV}}^2}{\dot{M}k{T}_0}.$$(15)

Employing Eqs. (8) and (9), the hydrodynamic Eqs. (12) and (13), valid for the atmospheric region lying below r_XUV, can be further reduced to (Watson et al. 1981) $$\frac{{\tau }^s}{\zeta }\frac{\text{d}\tau }{\text{d}\lambda }=\lambda -\frac{\beta }{\zeta {\lambda }_{\text{XUV}}^2}-\frac{5\tau }{2}-\frac{\Psi }{2}.$$(16)

By assuming r_XUV to be well below the sonic level, namely where the velocity of the gas reaches the sound speed, the velocity term can be set equal to zero. Furthermore, close to the lower boundary, λ ≫ τ. With these additional considerations, Eq. (16) can be further simplified to (Watson et al. 1981) $$\frac{{\tau }^s}{\zeta }\frac{\text{d}\tau }{\text{d}\lambda }=\lambda -\frac{\beta }{\zeta {\lambda }_{\text{XUV}}^2}.$$(17)

Finally, Watson et al. (1981) made assumptions regarding the shape of the temperature profile in the atmosphere. They argued that it has to be such that the escaping energy, the downward conduction of energy, and adiabatic cooling are in balance. A maximum particle escape rate can therefore be obtained by allowing the temperature in the thermosphere to drop to a 0 K minimum somewhere between the lower boundary and the XUV absorption region. This is because in this case the least possible amount of the absorbed energy is lost to downward conduction and the largest possible amount of absorbed energy is used to fuel the escape. Then, the fact that at the 0 K minimum both τ and $\frac{\text{d}\tau }{\text{d}\lambda }$ must be zero leads to the derivation of a system of equations that allows one to simultaneously calculate the escape parameter (ζ) and the Jeans escape parameter at the absorption height (λ_XUV) for this special temperature profile as (Watson et al. 1981) $$\zeta =\frac{2}{s+1}{\left[\frac{{\left(\frac{{\lambda }_{\text{XUV}}}{2}\right)}^{\frac{s+1}{2}}+1}{\Lambda -{\lambda }_{\text{XUV}}}\right]}^2$$(18)

and $${\lambda }_{\text{XUV}}=\sqrt{\frac{\beta }{\zeta }{\left[\Lambda -\sqrt{\frac{2}{(s+1)\zeta }}\right]}^{-1}}.$$(19)

Below we summarise the major assumptions and simplifications employed by Watson et al. (1981) to derive Eqs. (18) and (19).

2.1.1 Hydrogen-dominated atmosphere

The atmosphere is assumed to be dominated by hydrogen, meaning that above r₀ the mixing ratios of species heavier than hydrogen must be significantly smaller than one. Watson et al. (1981) used an arbitrarily defined lower boundary height (r₀) that lies well above the homopause, but in an atmospheric region that is still dense enough to have an optical depth to XUV much larger than one. However, as long as hydrogen is the dominant atmospheric gas, the extent of the heterosphere is negligible in comparison to the extent of the whole atmosphere, allowing one to use R_pl as the lower boundary height (r₀). In planetarysystems with extensive heterospheres, the energy-limited equations may only be used in regions well above the homopause.

2.1.2 Hydrodynamic escape

The XUV-driven hydrodynamic escape is assumed to be the dominating escape mechanism in the atmosphere. Other escape processes, such as Jeans escape or non-thermal escape processes (e.g. sputtering, ion pickup, and charge exchange), must be negligible. The hydrodynamic escape must also be fuelled entirely by the incident XUV energy of the host star.

2.1.3 Narrow absorption region

All the incident XUV energy is assumed to be absorbed in a narrow region near r_XUV. Furthermore, the optical depth to XUV radiation at r_XUV is equal to one.

2.1.4 Gas properties

The gas is assumed to be non-viscous and to have a constant mean molecular weight. The gas pressure is assumed to be isotropic and to decline towards zero at infinity. Watson et al. (1981) argued that allowing for the pressure at infinity to remain above zero would decrease the maximum particle escape rate because it applies an additional force working against the outward expanding gas.

2.1.5 No excess energy

The escaping particles are assumed to carry no excess energy away from the planet. This implies that each individual particle gets just enough energy to reach escape velocity, but no more. Watson et al. (1981) argued that additional excess energy would decrease the maximum particle escape rate because it would decrease the efficiency of XUV-driven hydrodynamic escape.

2.1.6 Tightly bound condition

The atmosphere at and below the lower boundary (r₀) is assumed to be ‘tightly bound’. This terminology, which was first introduced by Parker (1964), means that the atmosphere is quasi-static and expands subsonically upwards. Conversely, if the atmosphere is not tightly bound, the thermal velocity of the gas is comparable to or higher than the escape velocity, leading to supersonic catastrophic escape. Watson et al. (1981) defined the tightly bound condition in accordance with Parker (1964) to be met in all systems with Λ ≳ 10. Volkov et al. (2011) showed that catastrophic escape sets in at values of Λ < 2.4, but that significant deviations from tightly bound models can already be observed at Λ < 6.

2.1.7 Subsonic assumption

The XUV absorption height (r_XUV) is assumed to be below the sonic level. This assumption enables one to neglect the velocity term in Eqs. (12) and (13) and to further place limits to the temperature in the thermosphere. Watson et al. (1981) argued that this assumption may break down for strongly heated planets with high mass-loss rates.

2.2 Thermospheric temperature profile

The atmospheric temperature profile obtained through the energy-limited escape approach is one of the keys needed to understand the framework’s validity and limitations. As mentioned above, the energy-limited approach assumes that the temperature in the thermosphere drops to a 0 K minimum before the atmosphere is reheated by the absorbed incident XUV radiation. The energy acquired by XUV absorption is then redistributed within the gas as that heat is then used to drive the escape. We summarise here the mathematical description of the temperature profile in the region between the lower boundary and the XUV absorption height (r_XUV).

The starting point is Eq. (17), which is a differential equation that can be solved via integration. The conditions λ = Λ and τ = 1 are used for the lower boundary, while for the upper boundary the arbitrary values λ′ and τ′ are used. The solution is obtained through the following steps $${\tau }^s\text{d}\tau =(\lambda \zeta -\frac{\beta }{{\lambda }_{\text{XUV}}^2})\text{d}\lambda ,$$(20) $${\displaystyle {\int }_{\text{1}}^{{\tau }^{\prime }}{\tau }^s}\text{d}\tau ={\displaystyle {\int }_{\Lambda }^{{\lambda }^{\prime }}(}\lambda \zeta -\frac{\beta }{{\lambda }_{\text{XUV}}^2})\text{d}\lambda ,$$(21) $$\frac{{\tau }^{s+1}}{s+1}{|}_1^{{\tau }^{\prime }}=(\frac{{\lambda }^2}{2}\zeta -\frac{\beta }{{\lambda }_{\text{XUV}}^2}\lambda ){|}_{\Lambda }^{{\lambda }^{\prime }},$$(22) $$\frac{{{\tau }^{\prime }}^{s+1}}{s+1}-\frac{1}{s+1}=(\frac{{{\lambda }^{\prime }}^2}{2}\zeta -\frac{\beta }{{\lambda }_{\text{XUV}}^2}{\lambda }^{\prime })-(\frac{{\Lambda }^2}{2}\zeta -\frac{\beta }{{\lambda }_{\text{XUV}}^2}\Lambda ),$$(23)

and $$\frac{{{\tau }^{\prime }}^{s+1}-1}{s+1}=\frac{\zeta }{2}({{\lambda }^{\prime }}^2-{\Lambda }^2)-\frac{\beta }{{\lambda }_{\text{XUV}}^2}({\lambda }^{\prime }-\Lambda ).$$(24)

Solving Eq. (24) for τ′ yields the atmospheric temperature as a function of height in the atmosphere (λ′). For the sake of simplicity, all primed values can be replaced with unprimed values, resulting in $$\tau ={\{1+(s+1)[\frac{\zeta }{2}({\lambda }^2-{\Lambda }^2)-\frac{\beta }{{\lambda }_{\text{XUV}}^2}(\lambda -\Lambda )]\}}^{\frac{1}{s+1}}.$$(25)

The temperature profile we derived above is slightly different than that given by Watson et al. (1981). The difference lies in the choice of the lower boundary used for the integration of Eq. (21) that we consider to be the lower atmospheric boundary (i.e. the planetary radius), while Watson et al. (1981) considered it to be the location of the 0 K minimum. We modified the computation of the temperatureprofile for the following reasons: (1) the temperature profile described by Eq. (25) allows for the calculation of the atmospheric temperature at the XUV absorption radius (T_XUV), whether or not there is a 0 K minimum within the thermosphere; (2) the T_XUV values computed by employing Eq. (25) are physically more reasonable compared to those computed with the equation derived by Watson et al. (1981), which in several cases predicts T_XUV values reaching as high as 10⁵ K; (3) finally, Eq. (25) takes the physical conditions of the atmosphere at the lower boundary into account, in contrast to the equation derived by Watson et al. (1981). We remark that Eq. (25) meets all the assumptionsof the energy-limited formalism.

If the height of the XUV absorption (λ_XUV) is known or fixed, Eq. (25) can be used to investigate the dependence of the escape rate on the temperature profile. As long as the escape flux is below a critical value, the temperature will monotonically rise between Λ and λ_XUV. If the escape flux increases further, the remaining XUV energy, which is conducted down to the lower atmosphere, is not sufficient to counteract adiabatic cooling, and a temperature minimum will appear (Watson et al. 1981).

Equation (25) analytically describes the temperature profile only if λ_XUV is known, which therefore must be given as input. We remark that Eq. (25) is strongly dependent on the assumption that the height of the XUV absorption lies well below the sonic level, and thus the derived temperature profile is valid only below that point.

3 Consequences of approximating the location of the XUV absorption radius

Having described the details and assumptions involved in the energy-limited escape approach, we analyse now the resulting planetary properties (escape rate, absorption height, pressure, and temperature) obtained from applying it as a function of planetary mass, radius, and incident stellar XUV flux. To this end, we employed about 250 000 synthetic planets extracted from a planetary grid, whose characteristics are listed in Table 1. From this grid, we then considered only planets complying with the following four conditions: (1) 10 ≤ Λ ≤ 200, (2) average density between 30 and 10⁴ kg m⁻³, (3) R_Rl > 1.5 R_pl, and (4) $\log (-{\Phi }_G)=\log \left(\frac{G{M}_{\text{pl}}}{R_{\text{pl}}}\right)\le 9.11$, in accordance with Salz et al. (2016), who found that planets with higher gravitational potential become hydrodynamically stable.

For all calculations, we adopted a heating efficiency of ν = 0.15 (Shematovich et al. 2014; Salz et al. 2016; Kubyshkina et al. 2018a). We remark that the question of what ν value to adopt is still strongly debated. Attempts to constrain it for a range of planets have been performed for example by Shematovich et al. (2014) and Salz et al. (2016). These studies have found that, although slight variations are to be expected, ν remains around a value of roughly 20% for the planets considered in this work. However, the aim of this paper is not to provide an algorithm for calculating mass-loss rates, which would require an accurate estimate of ν for each individual planet, but rather to compare the results of the energy-limited approximation with results of hydrodynamic simulations, both of which require ν as input. Therefore, we chose a ν to match the one employed to compute the hydrodynamic simulations.

Often, the energy-limited approach is used to approximate the XUV absorption radius (r_XUV) with the planetary radius (R_pl; e.g. Ehrenreich & Désert 2011; Owen & Wu 2017; Dorn & Heng 2018; Jin & Mordasini 2018). With this additional approximation, the resulting energy-limited mass-loss rate is $${\dot{M}}_{{\text{R}}_{\text{pl}}}=\pi \frac{\nu {\Phi }_{\text{XUV}}{R}_{\text{pl}}^3}{KG{M}_{\text{pl}}}.$$(26)

The top panel of Fig. 1 shows the ratio of the mass-loss rate obtained by numerically solving Eqs. (18) and (19) (i.e. Ṁ_ζ) to the mass-loss rate obtained from Eq. (26) (i.e. ${\dot{M}}_{{\text{R}}_{\text{pl}}}$) and thus assuming r_XUV = R_pl as a function of Ṁ_ζ. For large values of Ṁ_ζ the two approaches lead to comparable mass-loss rates. However, ${\dot{M}}_{{\text{R}}_{\text{pl}}}$ underestimates the energy-limited mass-loss rate for values of Ṁ_ζ smaller than 10⁷ kg s⁻¹ by up to more than an order of magnitude. The middle panel of Fig. 1 shows how this behaviour appears in terms of the ratio of the XUV absorption radius (r_XUV) obtained with the energy-limited approach to the planetary radius (R_pl).

The bottom panel of Fig. 1 illustrates how the r_XUV/R_pl radius ratio behaves as a function of incident stellar XUV flux. The two radii get closer as the stellar XUVflux increases. This means that the increase in the mass-loss rate with increasing Φ_XUV is slowed down by the decreasing r_XUV.

Figure 2 shows the atmospheric temperature (T_XUV) and pressure (P_XUV) at the XUV absorption radius as a function of both stellar XUV flux and Λ. The T_XUV values range between 10³ and 50 000 K, and T_XUV increases with increasing stellar XUV flux and Λ.

The P_XUV values are not directly obtained from the energy-limited formalism, but they can be derived as follows. In a static and isothermal atmosphere, the condition defining the level of unit optical depth to XUV radiation is (Watson et al. 1981) $$n_{\text{XUV}}{H}_{\text{XUV}}{\sigma }_{\text{XUV}}\approx 1,$$(27)

where $$H_{\text{XUV}}=\frac{k{T}_{\text{XUV}}{r}_{\text{XUV}}^2}{KGMm}$$(28)

is the scale height at the XUV absorption radius, n_XUV is the particle density at the XUV absorption radius, and σ_XUV ≈ 2 × 10⁻²² m² is the absorption cross-section of atomic hydrogen to XUV radiation at an XUV photon energy of 20 eV (Murray-Clay et al. 2009). We remark that, because the atmosphere is not static but accelerating upwards and because the temperature is declining from r_XUV upwards, the true scale height is smaller than that given by Eq. (28), and thus Eq. (27) becomes (Watson et al. 1981) $$n_{\text{XUV}}{H}_{\text{XUV}}{\sigma }_{\text{XUV}}>1.$$(29)

However, Eq. (27) can still be used to estimate the pressure at the XUV absorption height that Fig. 2 shows lying in the 10⁻¹¹ to 10⁻⁹ bar range. The energy-limited approach indicates that P_XUV increases with increasing stellar XUV flux and Λ.

Fig. 1 Ratios of calculated mass-loss rates and XUV absorption radii using different versions of the energy-limited approximation. Top: ratio between the mass-loss rate obtained by numerically solving Eqs. (18) and (19) (Ṁζ) and the mass-loss rate obtained through Eq. (26) (${\dot{M}}_{{\text{R}}_{\text{pl}}}$) as a functionof Ṁζ. Middle: same as the top panel, but for the ratio between rXUV, obtained through solving Eqs. (18) and (19), and Rpl. Bottom: same as the middle panel, but as a function of the incident stellar XUV flux ΦXUV. For reference, the vertical line shows the XUV flux incident on HD 97658 b (Kubyshkina et al. 2018b). In each panel, the red line marks a ratio of one for reference.

Fig. 2 Temperature (top) and pressure (bottom) at the XUV absorption layer as a function of Λ and of the incident XUV flux. For reference, the purple X marks the position of HD97658 b (Kubyshkina et al. 2018b).

4 Comparison with hydrodynamic simulations

We comparehere the results obtained through the energy-limited approach with those derived from hydrodynamic atmosphericmodelling. Kubyshkina et al. (2018a) presented the results of one-dimensional hydrodynamic simulations of hydrogen-dominated planetary atmospheres for about 7000 synthetic non-magnetised planets. By construction, the hydrodynamic simulations account for Jeans escape, XUV hydrodynamic escape, and boil-off (i.e. core-powered mass loss). They also developed an interpolation routine that allows one to derive mass-loss rates (Ṁ_hc) for planets lying within the grid boundaries.

We employed the energy-limited formalism to derive Ṁ_ζ for roughly 97 000 planets contained within a grid we constructed by employing boundaries equal to those of the grid built by Kubyshkina et al. (2018a). The grid parameters are given in Table 2. Furthermore, this grid also complies with the conditions listed in Sect. 3, but we applied a stricter upper limit on Λ (i.e. Λ ≤80.0) to better match the grid of Kubyshkina et al. (2018a). We remark that Kubyshkina et al. (2018a) performed a comparison similar to that shown here, but with the key difference that Kubyshkina et al. (2018a) computed Ṁ_ζ by employing r_XUV values extracted from the hydrodynamic simulations, while we employed the r_XUV values derived from the energy-limited formalism (i.e. solving Eqs. (18) and (19)).

The top panel of Fig. 3 shows the Ṁ_ζ -to-Ṁ_hc mass-loss ratio as a function of the mass-loss rate obtained from interpolating across the grid of hydrodynamic models (Ṁ_hc). In general, the mass-loss ratio ranges over several orders of magnitude and decreases with increasing Ṁ_hc. The energy-limited formalism only reproduces the hydrodynamic results well for planets with moderate mass-loss rates, while it overestimates and underestimates escape by several orders of magnitude for weakly and strongly mass-losing planets, respectively. This is because for planets with low mass-loss rates the energy-limited formalism fails as a result of the hydrodynamic escape assumption being invalid. On the contrary, for planets with high mass-loss rates, the energy-limited formalism fails mostly as a consequenceof boil-off taking over as the main escape-driving mechanism (see also Fossati et al. 2017; Kubyshkina et al. 2018a). This is confirmed by the bottom panel of Fig. 3, which shows that for planets with low Λ values, which are prone to boil-off escape, the energy-limited formalism severely underestimates mass-loss rates.

Furthermore, Fig. 3 shows that the energy-limited formalism also fails to reproduce the results of the hydrodynamic simulations for some planets with moderate mass-loss rates and intermediate Λ values. We further explore this result in Fig. 4, which shows how Ṁ_ζ /Ṁ_hc varies as a function of both Λ and stellar XUV flux irradiation.

For planets with Λ smaller than ≈30, the energy-limited approximation on average significantly underestimates the mass-loss rate. As mentioned above, this is mostly due to the fact that, for planets with low Λ values, boil-off becomes the main escape-driving mechanism. In practice, the energy-limited approach assumes that the temperature profile presents a 0 K minimum. This implies that any thermal energy already available in the atmosphere as a result of, for example, gravitational contraction and/or optical and infrared stellar irradiation is not taken into account. For planets characterised by low Λ values, this thermalenergy is significant with respect to the planetary gravitational energy, which drives boil-off and falsifies the assumption set by the energy-limited approach on the temperature profile. In contrast, the energy-limited approximation appears to perform better at low Λ values and for high stellar XUV irradiation conditions. However, this stems from the fact that for such planets the energy-limited approach ignores boil-off and overestimates the upper atmospheric temperature, which have opposite effects on the mass-loss rate, finally returning mass-loss rates comparable to those given by the hydrodynamic simulations. We come back to this point in Sect. 5.

Figure 4 shows that for planets with higher Λ values, on average the energy-limited approximation reproduces the mass-loss rates given by hydrodynamic simulations within roughly an order of magnitude. This is the regime in which most of the assumptions involved in the energy-limited approach are valid. However, for a rather large number of planets characterised by low or high XUV flux irradiation levels, the energy-limited approximation appears to significantly overestimate mass loss, independently of Λ. This can be interpreted as follows. The main parameter setting the magnitude of the mass-loss rate of a given planet is the gravitational potential. Therefore, planets with different gravitational potentials, all else being equal, have significantly different mass-loss rates. However, Ṁ_ζ and Ṁ_hc behave differently as a function of the planetary gravitational potential (Φ_Grav). This is shown in Fig. 5, which presents, as an example, Ṁ_ζ and Ṁ_hc as a function of Φ_Grav for two sets of planets of the same mass, stellar XUV irradiation, and equilibrium temperature. In both cases, Ṁ_ζ and Ṁ_hc behave differently as a function of Φ_Grav. At low Φ_Grav values, Ṁ_hc is consistently higher than Ṁ_ζ as a result of the fact that Ṁ_ζ does not account for boil-off. Instead, at high Φ_Grav values, Ṁ_ζ is comparable to Ṁ_hc in the case of lower mass planets and significantly larger than Ṁ_hc in the case of higher mass planets. Furthermore, the relative behaviour of Ṁ_ζ and Ṁ_hc varies as a function of the other system parameters. This explains why planets lying in the same place in the Λ-Φ_XUV plane can present widely different Ṁ_ζ /Ṁ_hc values.

Therefore, the energy-limited escape approach is not suitable for describing planets with Λ values smaller than about 30 because the assumptions involved in the approach lead to an inappropriate description of the actual thermospheric temperature profile. Instead, for planets with higher Λ values, the energy-limited escape approach can be used to provide an upper limit on the mass-loss rates, but significant overestimation may be expected for some planets as a result of the improper dependence of the mass-loss rates on planetary gravitational potential. This indicates that the energy-limited formalism may not even be safe to use in cases where the main assumptionslisted in Sect. 2.1 appear to be met; we analyse this fact in more detail in the next section.

Fig. 3 Ratio of the mass-loss rates obtained following the energy-limited formalism (Ṁζ) and from hydrodynamic simulations (Ṁhc). Top: plotted as a function of Ṁhc. Bottom: plotted as a function of Λ. The dark green and bright blue points represent planets that do and do not comply with the subsonic condition, respectively(see Sect. 5.1). In both panels, the red and blue lines mark the regions within which Ṁζ and Ṁhc differ by less than a factor of two and ten, respectively.

Fig. 4 Ratio of the mass-loss rates obtained following the energy-limited formalism (Ṁζ) and from hydrodynamic simulations (Ṁhc) as a functionof Λ and incidentstellar XUV flux. For reference, the purple X shows the position of HD 97658 b (Kubyshkina et al. 2018b).

Fig. 5 Energy-limited (green dots) and hydrodynamic (red crosses) mass-loss rates as a function of gravitational potential (GMpl /Rpl) for two samples of planets with the same mass, stellar XUV irradiation, and equilibrium temperature (the parameters are listed at the top of each panel). The green and red lines show quadratic fits to the energy-limited and hydrodynamic mass-loss rates, respectively, in the form log Ṁ = $p_0+p_1{\Phi }_{\text{Grav}}+p_2{\Phi }_{\text{Grav}}^2$. The coefficients p0, p1, and p2 of the fits are listed in Table 3.

5 The potential use of the energy-limited approximation in estimating mass-loss rates

5.1 Inconsistencies within the energy-limited escape approach

Watson et al. (1981) employed a variety of assumptions to justify the simplifications that are necessary to derive the energy-limited approximation from the hydrodynamic equations. Therefore, it would appear that the energy-limited approximation can only be considered to reliably estimate mass-loss rates when the assumptions are met. However, this is not exactly the case. As a matter of fact, the energy-limited approximation sometimes appears to also reproduce the mass-loss rates provided by hydrodynamic simulations fairly well for planets that do not appear to comply with the energy-limited assumptions. For these cases, the match of the two mass-loss rates (i.e. Ṁ_ζ and Ṁ_hc) seems to occur by chance. In reality, this happens either because the true planetary conditions comply with the energy-limited assumptions, though the results of the energy-limited equations indicate the opposite, or because wrong assumptions within the energy-limited approach cancel one another out. The problem lies in the fact that without hydrodynamic simulations it is not possible to tell whether the energy-limited mass-loss rates, computed for planets for which the energy-limited approach indicates that the assumptions are not met, are correct or not. This nullifies the need to use the energy-limited approach. This is why without hydrodynamic simulations it is nearly impossible to determine a priori whether certain planetary conditions would be such that the energy-limited approximation returns a mass-loss rate compatible, or not, with that given by hydrodynamic simulations or to determine why it is even more difficult to conclude whether a certain planet would comply with the assumptions of the energy-limited approach. We present below a few examples describing what is stated above.

The subsonic condition, which assumes that the entirety of the stellar XUV energy is absorbed below the sonic point, is important within the energy-limited approach to justify neglecting the velocity term in the hydrodynamic equations. Within the energy-limited framework one can assess whether a planet complies with the subsonic assumption by using the following condition on the temperature at the XUV absorption radius (T_XUV), which is derived (in SI units) following Watson et al. (1981) and Eq. (25): $$T_{\text{XUV}}<T_{\text{limit}}=\frac{T_{\text{eq}}^2}{1.2605{(\zeta {\kappa }_0)}^2}.$$(30)

We remark that since this condition has been derived within the energy-limited framework, it is valid only as long as all the assumptions of the energy-limited approximation are met. This implies that it cannot be used to assess the applicability of the energy-limited approximation, but it can be used to check whether the energy-limited framework is consistent within itself once it is applied to a specific planet.

The bottom panel of Fig. 3 highlights in green planets that comply with the subsonic assumption and in light blue those that do not. Figure 3 shows that there is a large number of planets in the grid to which the subsonic assumption does not apply, though their energy-limited mass-loss rates are in agreement with those given by hydrodynamic simulations. For these planets, the physical conditions at the XUV absorption height would imply a supersonic expansion of the gas, meaning that the velocity term in the hydrodynamic equations can no longer be neglected. As shown by Fig. 3, however, this inconsistency does not automatically imply that the energy-limited approximation fails to correctly estimate the mass-loss rates.

We show this in practice by comparing, in Table 4, two planets that have slightly different masses (19 and 13 M_⊕) but similar values of Λ, similar Ṁ_ζ -to-Ṁ_hc mass-loss ratios (close to 1), similar equilibrium temperatures (r_XUV), and similar T_XUV values. Despite the many similarities, the heavier planet complies with the subsonic assumption, while the lighter planet does not. Since the subsonic assumption appears to be violated for the lighter planet, one would wonder why the energy-limited approximation returns a mass-loss rate in agreement with that of the hydrodynamic simulations. The reason is that the subsonic assumption is only violated within the energy-limited framework. In fact, because of the higher incident XUV flux of the lighter planet, the energy-limited approach artificially moves r_XUV closer to the planet, while in reality XUV absorption happens higher up in the atmosphere, where it is easier for the particles to escape, resulting in a more expanded atmosphere and similar mass-loss rates without the need to increase T_XUV. What happens is that the combinationof a larger r_XUV and a lower T_XUV allows for the velocity of the gas at r_XUV to remain subsonic, justifying neglecting the velocity term in the hydrodynamic equations. In this case, the first problem of the energy-limited approximation is not that it produces a significantly wrong mass-loss rate, but that it is not consistent within itself as it predicts a violation of its own assumptions. The second problem is that the user of the energy-limited approximation cannot realise this unless an independent, more reliable, and computationally expensive approach is used, cancelling out the benefits of using the energy-limited approximation.

Another assumption that does not seem to be met for a significant number of planets is the presence of a 0 K minimum in the thermospheric temperature profile. Watson et al. (1981) argue that, although not realistic, such a temperature profile is the most efficient at maximising the amount of energy used to drive escape. As a matter of fact, this is the assumption that led them to the fact that the energy-limited approximation returns an upper limit on the XUV-driven mass-loss rate. A steeper temperatureprofile would result in negative temperatures, which are not physical, while a shallower profile would drain more of the absorbed XUV energy into conduction, heating the lower parts of the atmosphere. However, for example in the case of planets with very high mass-loss rates, the energy-limited approximation estimates that r_XUV lies very close to the planetary surface (middle panel of Fig. 1) and that T_XUV is very high (top panel of Fig. 2). This results in extremely steep temperature profiles very close to the lower boundary, where the atmosphere is still very dense. Both conduction and convection in the collisionally dominated thermosphere make such asteep thermospheric temperature profile unlikely.

The results of such an inadequate prediction of the temperature profile can, for example, be seen in the case of planets with large Λ values that are subject to strong XUV irradiation. For these planets, the energy-limited approximation predicts that r_XUV lies close tothe lower boundary (bottom panel of Fig. 1) and that T_XUV reaches as high as 50 000 K (top-right corner of the top panel of Fig. 2). The resulting steep temperature profile is very efficient at driving particle escape, leading to high mass-loss rates. However, when looking at the results of hydrodynamic simulations (Fig. 4), one can see that the mass-loss rate for these planets is on average overestimated by up to three orders of magnitude. In reality, what happens is that conduction and convection distribute the absorbed stellar energy across several atmospheric layers, leading to atmospheric expansion, which in turn results in a larger r_XUV and a shallower thermospheric temperature profile, and thus in lower T_XUV and more moderate mass-loss rates. For these planets, the energy-limited approximation fails to correctly predict the thermospheric temperature profile. Hydrodynamic simulations predict that the thermospheric temperatures of planets with low equilibrium temperatures remain below 1000 K, while those of hot Jupiters do not exceed 20 000 K (e.g. Tian et al. 2005; Murray-Clay et al. 2009; Kubyshkina et al. 2018a). Therefore, there are cases for which the energy-limited approximation returns thermospheric temperatures that are significantly higher than what is predicted by more sophisticated models, which is caused by the fact that the energy-limited approximation is unable to correctly estimate temperature profiles and XUV absorption heights.

Finally, we briefly examine how well the energy-limited approximation estimates the location of r_XUV with respect to the atmospheric pressure, and thus optical depth to XUV radiation. The bottom panel of Fig. 2 shows that, according to the energy-limited approximation, the position of the XUV absorption height lies at pressures ranging between 10⁻⁹ and 10⁻¹¹ bar. Murray-Clay et al. (2009) found that in a hydrogen atmosphere P_XUV can be estimated as (in SI units) $$P_{\text{XUV}}=5.3\times {10}^{16}m\frac{KG{M}_{\text{pl}}}{R_{\text{pl}}^2},$$(31)

which typically results in pressures of about 10⁻⁹ bar, but they can also reach up to two orders of magnitude lower for planets with low equilibrium temperatures (~300 K) and two orders of magnitude higher for planets with high equilibrium temperatures. Therefore, there are planets within the considered grid for which the value of P_XUV estimated by the energy-limited approach is too low to be able to absorb the majority of the incoming XUV radiation.

5.2 Ranges of applicability of the energy-limited escape approach

To provide a more practical assessment of the parameter space in which the energy-limited approximation can be safely used to estimate atmospheric mass-loss rates, we computed the probability of obtaining Ṁ_ζ values comparable to those obtained by hydrodynamic simulations for different combinations of planetary parameters within the grid defined by Table 2. Figure 6 shows the probability of Ṁ_ζ and Ṁ_hc agreeing within a factor of two (top panel) or ten (bottom panel) as a function of Λ and Φ_XUV.

The top panel of Fig. 6 gives the impression that the energy-limited approximation can never be used to accurately predict mass-loss rates as there is essentially no region of the parameter space in which one can be sure of obtaining values within a factor of two of what is given by hydrodynamic simulations. Indeed, for the vast majority of the parameter space the probability of having Ṁ_ζ and Ṁ_hc agreeing within a factor of two is often lower than 50%. However, we remark that the energy-limited approximation has not been developed to correctly estimate mass-loss rates, but rather only to provide reliable upper limits. Furthermore, Watson et al. (1981) developed the energy-limited approximation to study the case of Earth’s atmospheric escape and could not foresee its use in a completely different context, such as that of exoplanets.

The bottom panel of Fig. 6 shows that the energy-limited approximation does a much better job when it is used to provide an order of magnitude estimate of the mass-loss rates. Indeed, there is a wide area within the considered parameter space in which Ṁ_ζ and Ṁ_hc agree within a factor of ten at least 75% of the time. However, there does not seem to be a safe range of system parameters in which there is the certainty that Ṁ_ζ and Ṁ_hc match within a factor of ten.

We assess here in greater detail the ranges of applicability of the energy-limited approximation. There are two distinct regions of the parameter space in which the energy-limited approximation appears to systematically provide incorrect mass-loss rates. The first one is located at the bottom-left corner of the bottom panel of Fig. 6 and is characterised by planets with low gravitational potential and low-to-intermediate XUV irradiation. The discrepancy in this regime between the energy-limited approximation and hydrodynamic simulations can be explained by looking at the resulting temperature profiles estimated within the energy-limited framework. We computed the linear temperature gradient from the predicted 0 K minimum to the XUV absorption height (i.e. the maximum temperature) as a measure of the steepness of the corresponding temperature profile. The position of the 0 K minimum is determined as (Watson et al. 1981) $${\lambda }_{\text{0K}}=\frac{\beta }{\zeta {\lambda }_{\text{XUV}}^2}.$$(32)

The linear temperature gradient is then given by $$\Delta T_{\text{0K}}=\frac{T_{\text{XUV}}}{R_{\text{XUV}}-R_{\text{0K}}}.$$(33)

The top panel of Fig. 7 shows the probability of the energy-limited mass-loss rates being within a factor of ten of those given by hydrodynamic simulations as a function of ΔT_0K and Λ. Again, the bottom-left corner of this plot (low Λ and low-to-intermediate ΔT_0K values) is populated by planets for which the energy-limited approximation returns mass-loss rates that are significantly different from those given by the hydrodynamic simulations. Therefore, we interpret this deviation as being related to the fact that the energy-limited approximation predicts too rapid a cooling of the lowest part of the atmosphere for planets that are not heated very strongly. As a matter of fact, the assumed 0 K minimum of the energy-limited approximation leads to the amount of energy already present in the lower atmosphere not being used to drive escape, forcing the atmosphere to cool to 0 K and then reheat up to the XUV absorption height. However, in reality, the lowest part of the atmosphere cools much more slowly and behaves nearly hydrostatically, leading to a more expanded atmosphere and an absorption of XUV energy farther from the lower boundary than predicted by the energy-limited approximation. The resulting temperature profile for planets with Λ ≈ 10 does not show a local temperature minimum but rather decreases monotonously from the lower boundary to the XUV absorption height. Furthermore, because of the low gravitational potential of these planets, the atmosphere already possesses almost enough energy at the lower boundary to drive the escape, and the additional energy from XUV absorption goes to increasing the already high mass-loss rates. Once the input energy (equilibrium temperature and XUV irradiation) decreases or the gravitational potential increases, the available energy is no longer sufficient to counteract adiabatic cooling at the lower boundary, leading to the appearance of a temperature minimum. The deeper the location of the minimum, the closer the mass-loss rate gets to what is predicted by the energy-limited formalism, generating a transition region in which the energy-limited estimates get increasingly more realistic with increasing Λ and/or Φ_XUV values. Also, even for planets with low gravitational potentials, if Φ_XUV reaches sufficiently high levels, the energy-limited approximation predicts small XUV absorption radii, leading to steeper temperatureprofiles that are closer to reality and, thus, more realistic mass-loss rate estimates.

The second distinct range of parameters in which the energy-limited approximation consistently fails in reproducing the results of hydrodynamic simulations is covered by planets characterised by high gravitational potential and high XUV irradiation. The middle and bottom panels of Fig. 7 show the probability of the energy-limited mass-loss rates being within a factor of ten of those given by hydrodynamic simulations as a function of Φ_grav, T_eq, and Φ_XUV, with the last two parameters being tightly related through the orbital separation. Planets with high gravitational potential and high irradiation are located in the top-right corner of each panel. We remark that the combinationof a high gravitational potential and high equilibrium temperature still results in moderate Λ values, which is why Fig. 6 does not highlight this part of the parameter space. These planets absorb a large amount of energy at the XUV absorption height, but their atmospheres are rather dense, because of the high gravitational potential, leading to an absorption height lying close to the lower boundary. The energy-limited approximation still assumes the presence of a 0 K minimum, but a significant amount of XUV radiation is then absorbed at a short distance from the temperature minimum, and most of this energy must be used to drive escape. If only a small amount of the absorbed energy is used to drive the escape, the atmosphere will conduct more energy towards the lower boundary, preventing the formation of the 0 K minimum. Therefore, the restriction imposed by the energy-limited approximation to downward conduction leads to unphysically high temperatures at the XUV absorption height (up to 50 000 K) and thus to excessively large mass-loss rates. In reality, thermal conduction is significant, particularly for these planets, and planetary atmospheres never reach a 0 K minimum, leading to lower temperatures at the XUV absorption height, shallower temperature profiles, and less energy used to drive escape, finally resulting in mass-loss rates lower than those predicted by the energy-limited approximation.

We conclude that in the case of low-to-intermediate equilibrium temperatures and incident XUV fluxes, the energy-limited formalism does not provide reliable mass-loss rates for planets characterised by excessively high or low gravitational potential. We find that the main reason for which the energy-limited mass-loss rates deviate most strongly from those given by hydrodynamic simulations is the inadequate representation of the temperature profile, mostly driven by the assumption of the presence of a 0 K minimum. Outside of those regions of the parameter space, the energy-limited mass-loss rates reproduce those given by hydrodynamic simulations fairly well in most cases. However, we remark that the boundaries within which the results of the energy-limited formalism are comparable to those of hydrodynamic simulations are not particularly sharp, and that there are extended transition regions in which energy-limited mass-loss rates are unreliable.

Fig. 6 Probability of Ṁζ agreeing with Ṁhc within a certain factor as a function of Λ and ΦXUV. The probabilities are computed considering planets lying within 30 Λ bins and 30 ΦXUV bins. Bins comprising fewer than ten planets have been discarded. Top: agreement within a factor of two. Bottom: agreement within a factor of ten.

Fig. 7 Probability of Ṁζ being within a factor of ten of Ṁhc as a function of different planetary parameters. Top: plotted as a function of Λ and ΔT0 K, i.e. the linear temperature gradient between the 0 K minimum of the energy-limited temperature profile and the XUV absorption height. Middle: plotted as a function of the gravitational potential of the planet corrected for tidal effects (Φgrav) and Teq. Bottom: plotted as a function of Φgrav and ΦXUV. For all panels the probabilities are computed considering planets lying within 30 bins of each parameter, except for the equilibrium temperature. Here just 17 bins have been considered as only 18 distinct values of Teq are represented in the planetary grid. Bins comprising fewer than ten planets have been discarded.

6 Conclusion

We reviewed the energy-limited atmospheric escape approach (Watson et al. 1981; Erkaev et al. 2007), which is widely used for estimating mass-loss rates for a broad range of planets that host hydrogen-dominated atmospheres, as well as for performing atmospheric evolution calculations. In particular, we listed and thoroughly described the assumptions employed by Watson et al. (1981) to simplify the hydrodynamic equations in the energy-limited formalism. We also studied the consequences of further approximating the energy-limited mass-loss formula by substituting the XUV absorption radius with the planetary radius. We found that this additional approximation leads to an underestimate of the energy-limited mass-loss rates by a factor of 25 at most. This result was obtained by considering planets that cover a very wide parameter space in terms of planetary mass, radius, and equilibrium temperature, the values of which range between 0.5 M_⊕ and 1 M_⊕, 1 R_⊕ and 1 R_⊕, and 300 and 2000 K, respectively.

Finally, we considered a smaller grid of planets (planetary mass, radius, and equilibrium temperature ranging between 1 M_⊕ and 39 M_⊕, 1 R_⊕ and 10 R_⊕, and 300 and 2000 K, respectively) to compare in detail the results of the energy-limited approximation with those of hydrodynamic simulations. For this comparison we considered the grid of models computed by Kubyshkina et al. (2018a), who also provided an interpolation routine that enables the extraction of key simulation results for any planet lying within the grid boundaries.

This comparison led us to conclude that the energy-limited approximation can be used as an order of magnitude estimate of the mass-loss rate for planets with intermediate gravitational potentials and low-to-intermediate equilibrium temperatures and XUV irradiation fluxes. The energy-limited formalism fails to approximate mass-loss rates for planets outside of this parameter space because for these planets it wrongly predicts the atmospheric temperature profile as a result of forcing the presence of a 0 K temperatureminimum in the thermosphere.

This result allows the energy-limited approximation to still be used for single planets, as long as their planetary parameters are within the parameter space of reliable estimation presented here. However, it precludes the use of the energy-limited approximation for planetary evolution calculations that require computing mass-loss rates for planets covering a wide parameter space; it is very likely that the approximation is going to fail in some of those cases, leading to significantly inaccurate planetary atmospheric evolutionary tracks (e.g. Kubyshkina et al. 2020). Therefore, in order to reliably estimate mass-loss rates over a broad range of parameters, we recommend using different tools, such as hydrodynamic models, or approximations obtained on the basis of hydrodynamic simulations, such as those presented by Kubyshkina et al. (2018a,b), which, being analytical, can be incorporated into planetary evolution models, as done for example by Kubyshkina et al. (2018a, 2019a,b, 2020), and Modirrousta-Galian et al. (2020a,b).

Acknowledgements

D.K. acknowledges funding received from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 817540, ASTROFLOW).

References

All Tables

All Figures

Fig. 1 Ratios of calculated mass-loss rates and XUV absorption radii using different versions of the energy-limited approximation. Top: ratio between the mass-loss rate obtained by numerically solving Eqs. (18) and (19) (Ṁζ) and the mass-loss rate obtained through Eq. (26) (${\dot{M}}_{{\text{R}}_{\text{pl}}}$) as a functionof Ṁζ. Middle: same as the top panel, but for the ratio between rXUV, obtained through solving Eqs. (18) and (19), and Rpl. Bottom: same as the middle panel, but as a function of the incident stellar XUV flux ΦXUV. For reference, the vertical line shows the XUV flux incident on HD 97658 b (Kubyshkina et al. 2018b). In each panel, the red line marks a ratio of one for reference.

In the text

	Fig. 2 Temperature (top) and pressure (bottom) at the XUV absorption layer as a function of Λ and of the incident XUV flux. For reference, the purple X marks the position of HD97658 b (Kubyshkina et al. 2018b).
In the text

Fig. 3 Ratio of the mass-loss rates obtained following the energy-limited formalism (Ṁζ) and from hydrodynamic simulations (Ṁhc). Top: plotted as a function of Ṁhc. Bottom: plotted as a function of Λ. The dark green and bright blue points represent planets that do and do not comply with the subsonic condition, respectively(see Sect. 5.1). In both panels, the red and blue lines mark the regions within which Ṁζ and Ṁhc differ by less than a factor of two and ten, respectively.

In the text

	Fig. 4 Ratio of the mass-loss rates obtained following the energy-limited formalism (Ṁζ) and from hydrodynamic simulations (Ṁhc) as a functionof Λ and incidentstellar XUV flux. For reference, the purple X shows the position of HD 97658 b (Kubyshkina et al. 2018b).
In the text

Fig. 5 Energy-limited (green dots) and hydrodynamic (red crosses) mass-loss rates as a function of gravitational potential (GMpl /Rpl) for two samples of planets with the same mass, stellar XUV irradiation, and equilibrium temperature (the parameters are listed at the top of each panel). The green and red lines show quadratic fits to the energy-limited and hydrodynamic mass-loss rates, respectively, in the form log Ṁ = $p_0+p_1{\Phi }_{\text{Grav}}+p_2{\Phi }_{\text{Grav}}^2$. The coefficients p0, p1, and p2 of the fits are listed in Table 3.

In the text

	Fig. 6 Probability of Ṁζ agreeing with Ṁhc within a certain factor as a function of Λ and ΦXUV. The probabilities are computed considering planets lying within 30 Λ bins and 30 ΦXUV bins. Bins comprising fewer than ten planets have been discarded. Top: agreement within a factor of two. Bottom: agreement within a factor of ten.
In the text

Fig. 7 Probability of Ṁζ being within a factor of ten of Ṁhc as a function of different planetary parameters. Top: plotted as a function of Λ and ΔT0 K, i.e. the linear temperature gradient between the 0 K minimum of the energy-limited temperature profile and the XUV absorption height. Middle: plotted as a function of the gravitational potential of the planet corrected for tidal effects (Φgrav) and Teq. Bottom: plotted as a function of Φgrav and ΦXUV. For all panels the probabilities are computed considering planets lying within 30 bins of each parameter, except for the equilibrium temperature. Here just 17 bins have been considered as only 18 distinct values of Teq are represented in the planetary grid. Bins comprising fewer than ten planets have been discarded.

In the text