© ESO 2022

1 Introduction

Several astrophysical phenomena are sustained by the conversion of gravitational potential energy of inflowing or outflowing matter into electromagnetic radiation, or vice versa. In this paper, we are concerned with the specific case where the primary energy source is in the form of stellar photo-ionizing radiation. The energy conversion-into heat-occurs within the atmosphere of a planet, where the ensuing thermal pressure gradient lowers the gas gravitational potential energy enough to drive an outflow (e.g., Lammer et al. 2003; Lecavelier des Etangs et al. 2004; Yelle 2004; Tian et al. 2005; Erkaev et al. 2007; Koskinen et al. 2007; Murray-Clay et al. 2009, amongst the seminal works). In this case, simple energy conservation balance yields $G{M}_{\text{p}}\dot{M}/R_{\text{p}}\simeq \pi R_{\text{XUV}}^2{F}_{\text{XUV}}$, where M_p and R_p are the planet mass and radius, F_XUV is the photo-ionizing¹ stellar flux at the planet orbital distance, R_XUV denotes the distance (from the planet) at which most of the stellar flux is absorbed, and Ṁ is the resulting mass outflow rate. It is interesting to note that the same equation holds in classic accretion theory, where the dissipation of gravitational potential energy of the accreting gas, which sinks in at a rate Ṁ, is partly converted into heat and generates an outward flux of radiation. Similar to the case of accretion flows, where much theoretical effort has gone into the determination of the radiative efficiency, bracketing the range of atmospheric outflow efficiencies remains a topic of active investigation (Lecavelier Des Etangs 2007; Lammer et al. 2009; Sanz-Forcada et al. 2011; Owen et al. 2010; Owen & Jackson 2012; Erkaev et al. 2013, 2016; Shematovich et al. 2014; Chadney et al. 2015; Khodachenko et al. 2015; Tripathi et al. 2015; Salz et al. 2016a,b; Owen & Alvarez 2016; Kubyshkina et al. 2018; Debrecht et al. 2019; Odert et al. 2020; Kubyshkina & Fossati 2021; Allan & Vidotto 2019; Vidotto & Cleary 2020).

The above energy conversion equation closely resembles the so-called energy-limited flux of escaping particles that was first derived by Watson et al. (1981) and applied to the specific case of Earth and Venus. The foundation of an energy-limited mass loss hinges on the role played by thermal conduction in limiting the radial extent of the thermosphere and thus the amount of stellar energy that is absorbed. Nevertheless, as noted by Owen (2019), the original connotation has been all but lost in the more recent literature. Rather, the energy-limited approximation is typically invoked to signify a theoretical maximum outflow rate which can be expected when a dominant fraction of the absorbed stellar radiation is converted into adiabatic expansion.

In practice, several unknowns hinder a straightforward estimate of this maximum Ṁ. Firstly, a non-negligible fraction of the stellar radiation that is absorbed by the atmosphere may go into exciting and ionizing its constituents. More broadly, radiative cooling competes with and may even offset adiabatic cooling. Second, the value of R_XUV cannot be readily estimated from first principles. Our ignorance of either parameter can be conveniently encapsulated in an efficiency term η. Replacing R_XUV with a (a priori unknown) multiple of the planetary radius gives the simplest form of the energy-limited approximation for atmospheric mass loss rate (see, e.g., Erkaev et al. 2007; Sanz-Forcada et al. 2011): $$\dot{M}=\eta \frac{3F_{\text{XUV}}}{4GK{\rho }_{\text{p}}},$$(1)

where G is the gravitational constant, ρ_p is the mean planetary mass density, and the potential energy reduction factor K < 1 accounts for the host star gravitational pull (Erkaev et al. 2007). Hereafter we shall refer to η as evaporation efficiency (we note that this is different from the heating efficiency defined by, e.g., Salz et al. 2016b, in that the latter does not include the unknown factorβ = R_XUV/R_p, which is embedded in our definition of η). A working definition of energy-limited thermal escape stems from adopting a fixed evaporation efficiency, typically in the range η ∈ [0.3 − 1] (see, e.g., Owen & Wu 2013; Odert et al. 2020). The limitations of this approximation are well documented in the literature; radiation hydrodynamics codes of varied degree of complexity all indicate that the energy-limited formula can overestimate the inferred mass outflow rates by several orders of magnitude (e.g., Krenn et al. 2021 and references therein; see, however, Cubillos et al. 2017). Whereas the evaporation efficiency of a single (or handful of) system can be derived by running dedicated photoionization hydrodynamics codes, a physically-motivated prescription for η becomes necessary for assessing the role of photoevaporation-driven mass-loss in planetary evolution studies (Baraffe et al. 2004; Ribas et al. 2005; Jin et al. 2014; Johnstone et al. 2015; Bisikalo et al. 2018), as well as its possible contribution to carving the observed gap (Fulton et al. 2017) in the radius distribution of small planets (Owen & Wu 2013; Lammer et al. 2013; Lopez & Fortney 2013; Rogers et al. 2021).

In a seminal paper, Murray-Clay et al. (2009; hereafter MC09) simulated the atmospheric escape from a Jovian-like planet subject to different irradiation levels, spanning over 5 dex in flux. They show that the energy-limited approximation tends to overestimate the inferred mass outflow rate above a flux threshold of about ≃10⁴ erg cm⁻² s⁻¹. This is interpreted as due to the onset of radiative processes (primarily in the form of hydrogen collisional excitation, that is, Lyα emission) as the dominant cooling mode. Owen & Alvarez (2016) further expanded on this and cast the transition off of energy-limited escape in terms timescales, namely the recombination time becoming comparable and eventually much shorter than the flow timescale.

Separately, Salz et al. (2016b) focused on the role of planetary gravity; they did so by running The Pluto-Cloudy Interface (TPCI; Salz et al. 2015) for a sample of more than 20 (known and artificial) planetary systems covering a wide range of masses and mass densities. This experiment reveals a sharp transition from energy-limited escape, which is valid up to specific gravitational potential energies as high as log(−Φ) ≲ 13.1 (in cgs units), to a regime where the evaporation efficiency declines sharply with increasing |Φ|.

It is worth noting that the decline in outflow efficiency with flux seen by MC09 is not nearly as sharp as that seen by Salz et al. (2016b) as a function of gravitational potential, suggesting that the latter plays a primary role in determining the nature of the outflow. Furthermore, at log(−Φ) = 12.94, the test planet considered by MC09 is extremely close to the threshold value identified by Salz et al., suggesting that flux plays a non-negligible (if secondary) role in the overall energy-balance.

Although several other investigations conclude that both gravity and stellar irradiation ought to have an effect on the overall process (e.g., Wang & Dai 2018; Kubyshkina et al. 2020; Kubyshkina & Fossati 2021; Krenn et al. 2021; Lampón et al. 2021, among the most recent ones), we still lack a coherent physical understanding of the limits of energy-limited escape. This is indeed the aim of this Paper. First, we introduce our sample and data (Sect. 2). Next, we focus on the effect of planetary gravity, and demonstrate analytically the existence of a threshold value above which energy-limited escape cannot take place-irrespective of the stellar irradiation level (Sect. 3).

Last, we carry out a systematic investigation of the dependence of the mass outflow rate on stellar flux over a broad range of fluxes and elucidate how the different trends exhibited by low and high-gravity planets can be broadly understood in terms of the different fractional contribution of the radiative cooling mechanisms at play (Sect. 4). We end with a summary of our main results and conclusions (Sect. 5). An analytical approximation of the best-fitting evaporation efficiency is given in the Appendix.

2 Planet sample and simulated data

We draw from a sample of 25 nearby ( ≲ 100 pc), moderately to highly irradiated gaseous planets, from sub-Neptunes² through hot Jupiters, examined by Spinelli et al. (in prep.) as part of Paper III of this series. Compared to this parent sample, we limit our investigation to 16 systems with (i) R_p > 0.15 R_J (where the subscript J indicates Jovian units), to justify the assumption of a large ( ≳ 1%) hydrogen-helium envelope (Lopez & Fortney 2014); (ii) X-ray detected host stars; and (iii) irradiation and specific gravitational potential values within the converge range of ATES (ATmospheric EScape; a publicly available photo-ionization hydrodynamics code that was presented in Paper I of this series; Caldiroli et al. 2021).

For each target Spinelli et al. present revised stellar and orbital parameters based on Gaia DR2 parallactic distances (Gaia Collaboration 2018), perform a re-analysis of all the archival X-ray (Chandra and XMM-Newton) data for the purpose of uniformly estimating F_XUV (using the scaling relations of King et al. 2018), and make use of ATES to compute the outflow temperature, ionization and density profiles, and steady-state mass loss rate. The reader is referred to Table 1 for a list of the parameters that are most relevant for this work (see Tables 1 and 3 in Spinelli et al. for a complete list of parameters).

Our targets span 1.4 dex in specific gravitational potential energy ϕ_p ≡ GM_p/R_p (hereafter defined as positive) vs. only 0.8 dex in F_XUV. To fully parse the combined effects of planetary gravity and stellar irradiation we simulated a broader range of F_XUV for a sub-sample of 8 target planets, namely GJ 3470 b, HAT-P-11 b, WASP-69 b, HD 149026 b, WASP-80 b, HD 189733 b, WASP-77A b and WASP-43 b (listed in order of increasing ϕ_p). For each of those we aimed to run ATES adopting F_XUV values spanning between 10² and 10⁶ erg cm⁻² s⁻¹, in intervals of 0.5 dex, for a total of 9 simulated fluxes per planet. In practice, whereas ATES can be run successfully across the whole flux range for GJ 3470 b, HAT-P-11 b, WASP-69 b, HD 149026 b and WASP-80 b the code fails to converge for the lowest-flux values in the case of the remaining (higher-gravity) planets, since the order of magnitude of the velocity acquired by the gas near the planet surface is comparable to the order of approximation of the numerical scheme. This limitation is fully consistent with the code convergence criterion laid out in Paper I of this series, namely for systems with log ϕ_p ≳ 12.9 + 0.17 log F_XUV (where the planet gravitational potential and stellar flux are expressed in cgs units). As a result, the lowest simulated F_XUV value is 10⁴, 10³ and 10²⁵ erg s⁻¹, respectively for WASP-43 b, WASP-77A b and HD 189733 b, for a total of 81 simulations. Figures 1a and b illustrate the simulation results for the case of a prototypical Neptune-like planet (GJ 3470 b, for which the simulations span over 4 dex in flux) and hot-Jupiter (WASP-77A b, 3 dex in flux), respectively. Whereas the latter type exhibits a sharp ionization front and small range of temperatures across the flux range, the ionization and temperature profiles of the former are highly sensitive to the irradiation level. At the highest irradiation level, both reach the typical temperature of a fully-ionized gas in photo-ionization equilibrium (≃10⁴ K).

ATES’ results - both those arising from the actual sample listed in Table 1, as well as those ensuing from the simulated FXUV values - constitute the primary data set for our investigation below. Throughout, we illustrate the former using black symbols, and the latter using colors, with bluer (resp. redder) colors corresponding to lower (higher) irradiation, following the same color scheme as laid out in Figs. 1a and b.

Fig. 1 Simulated outflow properties for the Neptune-like planet GJ 3470 b (Panel a) and for the hot Jupiter WASP-77A b (panel b). FXUV varies between 102 (deepest blue) and 106 erg cm−2 s−1 (deepest red), in intervals of 0.5 dex. The diamond symbols mark the location of the τ = 1 surface at 13.6 eV.

Fig. 2 Outflow effective efficiency is plotted against the planet specific gravitational potential for the 16 planets considered in this work (see Table 1 for a list); ηeff is defined as the ratio between the ATES-estimated mass outflow rate and the energy-limited value given by Eq. (1) for η = 1. To explore the effects of varied irradiation levels on Ṁ, (81) simulations were carried out for a broader range of FXUV for a sub-sample of eight targets. The results are represented by the colored circles, where the adopted color scheme is the same as in Figs. 1a and b.

3 The role of planetary gravity: ϕ_p-limited thermal escape

Following Salz et al. (2016b), in Fig. 2 we plot the outflow “effective efficiency”, η_eff, defined as the ratio between the mass outflow rate estimated by ATES and the energy-limited expression given by Eq. (1) for η = 1, as a function of ϕ_p. This quantity enables a swift estimate of how far off from the energy-limited expression a given system is. To the extent that ATES correctly models all of the relevant physics, it also yields a quantitative estimate of the actual energy-limited evaporation efficiency.

This exercise confirms the existence of a limiting ϕ_p value above which the effective efficiency begins to dive well below theoretical energy-limited value. Our results are qualitatively and quantitatively consistent with the conclusions drawn by Salz et al. (2016b), who interpret this behavior in terms of a transition from thermally-driven outflows to a regime where radiative cooling becomes comparable to (and eventually dominates over) adiabatic cooling, thus making the outflow progressively more and more inefficient compared to the energy-limited approximation. For log ϕ_p ≲ 13 the bulk of the absorbed stellar flux is mostly converted into expansion work. Working against a weak gravitational pull, this yields efficient atmospheric escape. In contrast, higher gravitational potential planets host more tightly bound atmospheres and a dominant fraction of the absorbed stellar flux is re-emitted locally, primarily through Lyα and bremsstrahlung emission.

This widely accepted scenario, while phenomenologically correct, falls short of a quantitative explanation for the existence of such a sharp threshold, as opposed to a more gentle decline in efficiency. The heretofore unexplained numerical value of the threshold is also of significant interest. In what follows we show analytically that the existence of a limiting ϕ_p at ≃ 10¹³ erg g⁻¹ can be understood in terms of a quantitative balance between the atmospheric particles’ binding energy, the photo-ionization energy budget, and the remaining energy reservoir that is available to initiate and sustain the expansion.

Let us consider the effects of an external photo-ionizing source (that is, the stellar radiation field) on an thin layer of atmospheric gas of a planet with specific gravitational potential energy ϕ_p. Borrowing from cosmological re-ionization theory (see, e.g., Puchwein et al. 2019), the mean excess energy ∆E_xs within this layer can be defined as the ratio between the photo-heating rate (erg s⁻¹) and the photo-ionization rate (s⁻¹) per ion: $$\text{Δ}{E}_{\text{XS}}=\frac{{\displaystyle {\int }_{v_t}^{\infty }\text{d}v\frac{F_v}{hv}{\sigma }_{v}h\left(v-v_t\right)}}{{\displaystyle {\int }_{v_t}^{\infty }\text{d}v\frac{F_v}{hv}{\sigma }_v}},$$(2)

where F_v = F_v,0 exp(−τ_v) is the local, absorbed stellar flux at frequency v, σ_v is the photo-ionization cross section, τ_v is the optical depth, v_t is the photo-ionization frequency threshold and h is the Planck constant. In practical terms, this quantity represents the amount of energy that is still available to heat up the gas after a single photo-ionization event (see, e.g., Osterbrock & Ferland 2006). We postulate that, for an atmospheric particle of mass m to be able to escape the system under the effect of photo-heating, the volume-averaged mean excess energy must exceed the particle binding energy: 〈∆E_xs〉 ≳ ϕ_pm.

The crux of this approach is that ∆E_xs is independent of the normalization of the photo-ionizing flux, as it only depends on the shape of the radiation spectrum close to the ionization threshold of the element(s) under consideration. For simplicity, we start by considering the case of a pure hydrogen atmosphere. To first order the shape of the photo-ionization cross-section close 13.6 eV scales with the frequency as v^−2.75, such that, given an external, power-law shaped photo-ionizing radiation spectrum ∝ v^−ß, the mean excess energy can be approximated as ∆E_xs ≃ 13.6/(ß + 1.75) eV. In the optically-thin regime, such as for a typical unabsorbed stellar spectrum with ß ≃ 1, ∆E_xs ≃ 5 eV. However, owing to the effect of atmospheric absorption, the spectrum becomes progressively harder deeper into the atmosphere. In fact, F_v becomes nearly flat above 13.6 eV just across the ionization front (see, e.g., Fig. 5 in Haardt & Madau 1996 for a cosmological application), implying that ∆E_xs ≃ 8 eV in the layer where the bulk of stellar radiation is absorbed.

For a more robust, quantitative estimate we used ATES to simulate the atmospheres of eight case studies, spanning 1.4 dex in ϕ_p. As illustrated in Fig. 3a, ∆E_xs varies within a remarkably narrow range (between 5 and 7 eV) within the optically thin portion of the outflow (defined by τ ≤ 1 at 13.6 eV). Only the most energetic (and thus rare) photons penetrate deeper, into the optically thick portion, producing (rare) ionization events with extremely high values of ∆E_xs. Using the atmospheric density and ionization profiles obtained by ATES for the same eight planets, also under different irradiation regimes (for a total of 81 simulations, spanning 4 dex in F_XUV) we calculated the volume-averaged mean excess energy by integrating ∆E_xs over the free-electron number density. This yields values in the range 〈∆E_xs〉∈ [8−15] eV.

The minimal condition for the onset of a thermally-driven outflow can thus be recast in terms of an upper limit to the specific gravitational potential. By setting m equal to the proton rest mass m_p we obtain: ${\varphi }_{\text{p}}\lesssim {\varphi }_{\text{p}}^{\text{thr}}\approx [7.7-14.4]\times {10}^{12}\text{erg\hspace{0.17em}}{\text{g}}^{-1}$, or log ${\varphi }_{\text{p}}\lesssim \log {\varphi }_{\text{p}}^{\text{thr}}\approx [12.9-13.2]$ (in cgs units).

Several effects complicate the matter, such as the presence of heavier elements than hydrogen (primarily helium), the possible presence of hydrogen molecules in the lower atmosphere, and the fact that the outflow is launched past R_p, implying a slightly reduced binding energy. Each is expected to carry a correction factor of order a few to the above criterion. Overall, we find a remarkably good agreement with the empirical value of ~ 10¹³ erg g⁻¹ found by Salz et al. (2016b) and subsequent work (including the estimates presented in Fig. 2), all of which make use of radiation hydrodynamics codes that properly model the wind launching mechanism and also account for the presence of helium. Hereafter, we shall loosely refer to planets with log ϕ_p < 12.9 and log ϕ_p > 13.2 as low-and high-gravity planets, respectively.

Introducing the concept of mean excess energy enables a novel definition of the outflow launching point, as follows. At each radius, ∆E_xs represents the mean energy carried by a photo-electron; this energy is thermalized locally through collisions, such that the mean kinetic energy gained by a proton in the process can be written as: $$m_{\text{p}}\frac{{\upsilon }_{\text{XS}}^2}{2}\equiv \frac{n_e}{n_{\text{H}}}\text{Δ}{E}_{\text{XS}}.$$(3)

Here n_H and n_e are the hydrogen and free electron number density, respectively, and the “excess velocity” υ_xs is the mean velocity that is acquired locally by a proton. The outflow launching point is identified by the radius (as in distance from the planet center) at which the excess velocity exceeds the escape velocity: (υ_xs/υ_esc)² ≳ 1, where υ_esc is calculated through the following expression: $$\frac{{\upsilon }_{\text{esc}}^2\left(r\right)}{2}=\left[\text{Φ}\left(R_{\text{RL}}\right)-\text{Φ}\left(r\right)\right].$$(4)

Φ(r) is the gravitational potential of the planet-star system, defined in equation 2 of Caldiroli et al. (2021), R_RL refers to the Roche Lobe radius, while r is the distance from the planet center.

The above arguments can be easily generalized for a mixed hydrogen-helium atmosphere. In this case, Eq. (3) becomes: $$\mu \frac{{\upsilon }_{\text{XS}}^2}{2}\equiv f_{\text{HII}}\text{Δ}{E}_{\text{HI}}+Y\left(f_{\text{HeII}}\text{Δ}{E}_{\text{HeI}}+2f_{\text{HeIII}}\text{Δ}{E}_{\text{HeII}}\right),$$(5)

where H I and H II denote neutral and ionized hydrogen, while He I, He II and He III denote neutral, single-ionized and double-ionized helium; f_i and ∆E_i are the ionization fraction and mean excess energy for the relevant i species, Y ≡ n_He/n_H is the helium to hydrogen number abundance ratio and µ ≡ (1 + 4Y)m_p.

Figure 3b compares the location of the outflow launching point as identified by the (υ_xs/υ_esc)² ≳ 1 criterion with the location of the τ = 1 surface (calculated at 13.6 eV throughout), for the same eight planets shown in Fig. 3a. Interestingly, our newly defined outflow launching point approximately coincides with the τ = 1 surface only for the intermediate-gravity planet WASP-80 b, whereas it is located within (outside) the τ = 1 surface for low- (high-) gravity planets. This has profound consequences for the outflow evaporation efficiency, since the τ = 1 surface sets the location where the stellar energy absorption starts to become significant (in fact, the τ = 1 surface is often adopted as a working definition for the outflow launching point in the literature). For low-gravity planets, the shallow ionization profiles (see, e.g., Fig. 1) guarantee that a sizable fraction of the absorption occurs within an extended layer located between the τ ≳ 1 and the υ_xs ≃ υ_esc radius. This means that most of the energy is disposed of where (on average) the ions have sufficiently high kinetic energy to escape the system. As illustrated in the next Section, this regime corresponds to evaporation efficiencies in the range ∈[0.1−1], depending on the irradiation level.

In contrast, the sharp ionization fronts of high-gravity planets (see, e.g., Fig. 1) imply that only a small fraction of the stellar flux is absorbed above the (υ_xs/υ_esc)² = 1 line, where the ion-ization fractions are high and the densities are low. Most of the energy absorption occurs within a very narrow layer where (after thermalization) the average ion kinetic energy is insufficient to warrant escape. Combined, these effects are responsible for the sharp decline in evaporation efficiency above ϕ_p, down to values as low as 10⁻⁴.

Fig. 3 H I mean excess energy defined in Eq. (2) (panel a) as a function of planetary distance for different systems and Second power of the ratio between the excess and escape velocity defined in Eqs. (3) and (4) (panel b); r and RRL are expressed in units of the planet radius, such that the planetary distance on the x axis is given as a fraction of the Roche lobe radius RRL. The dashed line represents the threshold for wind launching, i.e., where (υxs/υesc)2 = 1. In both panels, the diamond symbols indicate the location of the τ = 1 surface at 13.6 eV.

Fig. 4 Simulated mass loss rate Ṁ for an artificial planet with Mp = 1.3 × 1030 g and Rp = 1010 cm, subject to different levels of stellar XUV irradiation (panel a) and simulated mass outflow rates for eight planets subject to different levels of stellar XUV irradiation, for a total of 81 simulations (panel b). In panel a, the dashed red line represents the results of the simulations carried out by Murray-Clay et al. (2009), under the assumption of a pure hydrogen atmosphere and mono-energetic irradiation field at 20 eV. The blue line represents the results obtained by ATES under the same assumptions. The green line represents the ATES results for a pure hydrogen atmosphere subject to a power-law irradiation spectrum with slope v−1. The black line refers to a primordial hydrogen-helium atmosphere with a power-law spectrum, that is, the standard setup for ATES. The adopted color scheme in panel b is the same as shown in Figs. 1a, b and 2. In both panels, the dotted lines represent linear proportionality between mass loss rate and flux.

4 Flux-dependent evaporation efficiency

We turn our attention to the functional dependence of the outflow rate on irradiation. MC09 showed that for an artificial planet with M_p = 1.3 × 10³⁰ g and R_p = 10¹⁰ cm (or log ϕ_p = 12.94) the transition from energy-limited thermal escape to a regime where radiative cooling starts to dominate occurs close to ≃10⁴ erg cm⁻² s⁻¹. To start with, we verify that ATES reproduces the same results as found by MC09. To do so, we modify the code to simulate a monochromatic radiation field and pure hydrogen atmosphere. The results of this comparison are illustrated in Fig. 4a. While relaxing those approximations results in higher Ṁ across the flux range, we recover the same behavior as described MC09, namely that Ṁ scales linearly with F_XUV below ≲ 10⁴ erg cm⁻² s⁻¹. The slope flattens above it, signaling a departure off of energy-limited thermal escape.

MC09 attributed the change is slope in the Ṁ:F_XUV relation seen around F_XUV ≃ 10⁴ erg cm⁻² s⁻¹ to the transition from an energy-limited regime to one where radiative cooling, specifically in the form of Lyα, dominates the flow. This argument was expanded upon by Owen & Alvarez (2016), who differentiate between an energy-limited regime and a recombination-limited regime. In the former the recombination time is much longer than the flow timescale, the outflow exhibits a very extended ionization front, and the mass outflow rate can be expected to scale linearly with the ionizing flux. This regime is typical of low-gravity planets at low irradiation. At high irradiation levels the ionization front becomes sharp, the flow is in a radiative-recombination equilibrium state, and the mass-loss rate scales approximately as the square root of the ionizing flux. Owen & Alvarez (2016) also show that the transition is rather slow, in the sense that it occurs over a large range of fluxes, with higher gravity planets transitioning to the recombination-limited regime at lower fluxes (we note that the highest-gravity planet simulated by Owen & Alvarez 2016 has log ϕ_p = 12.6, and is thus in the low-gravity regime as per Sect. 3).

In order to further investigate the relation between Ṁ and F_XUV across the full range of ϕ_p and F_XUV we rely on the simulations described in Sect. 2. As a reminder, these include four low-gravity planets (GJ 3470 b, WASP-69 b, HAT-P-11 b and HD 149026 b), three high-gravity planets (HD 189733 b, WASP-77A b and WASP-43 b), plus WASP-80 b, which, at log ϕ_p = 13.02 can be considered as an intermediate case. All of those were simulated over a wide range of F_XUV spanning from 10² to 10⁶ erg cm⁻² s⁻¹, in intervals of 0.5 dex, for a total of 81 simulations. The results are illustrated in Fig. 4b. Overall, it is interesting to note how different planets exhibit different slopes as well as normalizations – both compared to the Ṁ: F_XUV relation seen for the test-planets considered by MC09 (Fig. 4a) and Owen & Alvarez (2016), as well as different from one another. This is best seen in Fig. 5, where we plot the K-reduced mass outflow rate³ as a function of the stellar irradiation to planetary density ratio (F_XUV/ρ_p) for the sample under consideration. Since, for an energy-limited outflow, the product KṀ is expected to scale linearly with this quantity along the η = 1 line (Eq. (1)) downward deviations readily pinpoint any departures from the energy-limited formalism. Our simulations confirm that energy-limited escape is attained by low-gravity planets only up to F_XUV/ρ_p ≲ 10^3–3.5g cm s⁻¹ g⁻¹. The curves flatten at higher fluxes, and the outflow efficiency decreases (by a factor up to 10) with increasing irradiation.

In contrast, the high-gravity planets all exhibit much lower effective efficiencies (η ≲ 0.05). However, their outflows tend to become more efficient for progressively higher irradiation values, that is the slope of the KM vs F_XUV/ρ_p curves for high-gravity planet is steeper than linear (the intermediate gravity planet WASP-80 b –solid circles in Fig. 4b – exhibits somewhat of a hybrid behavior). The same effect can be visualized by looking at the colored points in Fig. 2, where different colors correspond to different irradiation levels (the adopted color scheme is the same as in Figs. 4b and 5). Whereas low-gravity planets show a decrease in effective efficiency with increased fluxes the opposite is true for high-gravity planets. We refer to this behavior as ‘evaporation efficiency inversion’.

Next, we examine this phenomenon in terms of fractional contributions from the competing cooling processes to the overall energy budget, and specifically adiabatic, advection and radiative cooling. This is shown in Fig. 6, where the volumetric contributions of each are plotted as a function of flux for four case studies.

Fig. 5 Simulated, K-reduced mass outflow rate, shown as a function of the irradiation to planet mass density ratio, for the same planets shown in Fig. 4b. Energy-limited escape yields KṀ = FXUV/ρp (i.e., η = 1 in Eq. (1)). The adopted color scheme is the same as shown in Figs. 1a, b and 2.

4.1 Low-gravity planets: flux-limited thermal escape

We start by discussing the behavior of GJ 3470 b as our low-gravity case study (log ϕ_p = 12.3). Planets such as GJ 3470 b exhibit a wide range of atmospheric temperatures and ioniza-tion fractions as the irradiation level changes. The atmospheric temperature at the τ = 1 surface varies by about one order of magnitude as the flux increases from 10² to 10⁶ erg s⁻¹ cm⁻² (see Fig. 1a). At the lowest flux, the ionization front is very broad and virtually 100% of the absorbed stellar radiation is converted into adiabatic expansion (top left panel of Fig. 6). Adiabatic cooling (dashed green line) starts to dive (albeit gently) below ≃100% already at ≃10³ erg s⁻¹ cm⁻². Above this flux level, the flow ceases to be formally energy-limited and its efficiency decreases with increasing flux owing to the progressive increase of (i) advective cooling (solid green line) and (ii) radiative cooling (solid blue), primarily in the form of collisional excitation (Lyα). The fractional contribution of advective cooling (which scales a ṀT/µ) peaks at ≃10⁴ erg s⁻¹ cm⁻² (where it accounts for 26% of the total cooling), and decreases at higher fluxes. This happens as a result of the fact that the flow temperature approaches the thermostat value of 10⁴ K, while the flow approaches full ionization.

However, Lyα starts to skyrocket⁴ above recombinations and bremsstrahlung losses right around the same flux where advective cooling reaches its peak fractional contribution. At F_XUV ≳ 10⁵ Lyα dominates over advective cooling.

As the flux further increases the outflow efficiency decreases at a progressively greater pace (down to about 10%) due to the nonlinear growth of the Lyα cooling contribution. While the actual flux above which Lyα losses dominate over advec-tive cooling depends on the planet gravity, the same qualitative description applies to all low-gravity planets (log ϕ_p ≲ 12.9). For those, advective and radiative (Lyα) losses, combined, drive the progressive decline of the evaporation efficiency with increasing irradiation. More specifically, as also shown by Owen & Alvarez (2016), the transition toward a radiation-dominated regime occurs at lower fluxes for higher gravity planets (top right panel of Fig. 6 for the case of HD 149026 b, with log ϕ_p = 12.8).

4.2 High-gravity planets

As we argue in Sect. 3, high-gravity planets can never sustain energy-limited outflows. As illustrated in Fig. 1b for the case of WASP-77A b, the outflow density is typically so low that the ionization structure resembles an inverted Strömgren sphere, exhibiting a very sharp transition between an outer fully ionized region and an inner, neutral region (this happens even at relatively low irradiation, contrary to the case of light planets). Irrespective of the irradiation level, the atmosphere temperature close to the ionization front approaches (1−2) × 10⁴ K, that is, the thermostat temperature of a highly ionized gas in ion-ization equilibrium. As a consequence Lyα cooling approaches a saturation level, in the sense that, unlike for the low-gravity planets, increasing the irradiation level does not result in a sensible increase of the fractional contribution in Lyα cooling. This can be seen in the bottom left panel of Fig. 6. Lyα accounts for close to 90% of the overall cooling across the whole flux range; if anything, its fractional contribution decreases slightly, from 92 to 88%, in going from F_XUV = 10³ to 10⁶ erg cm⁻² s⁻¹, while the fractional contributions from recombinations and bremsstrahlung increase from 3 to 5% and from 3 to 4%, respectively, over the same range. The remaining energy (1−2%) goes into adiabatic expansion; advection losses are completely negligible in this regime.

The case of WASP-80 b (bottom right panel of Fig. 6) is intermediate between low- and high-gravity planets. The relative contributions from adiabatic and radiative cooling are similar up to F_XUV ~ 10⁴ erg cm⁻² s⁻¹, above which Lyα losses take over. Advective losses are always subdominant. This causes the efficiency to increase ever so slightly with flux up to F_XUV ~ 10⁴ erg cm⁻² s⁻¹, above which it starts to decrease with increasing flux, thus resembling low-gravity planets.

In summary, as the irradiation flux increases, the fractional contribution of adiabatic cooling to the overall process - however small - grows more rapidly than the radiative component. As a result, the evaporation efficiency - however small - tends to increase with flux in high-gravity planets.

Fig. 6 Volume integrals of the various cooling processes of our model as functions of the impinging XUV flux. From top-left, moving clockwise, the four panels represent different planets with increasing φp, from : GJ 3470 b (log φp = 12.30), HD 149026 b (log φp = 12.84), WASP-80 b (log φp = 13.02) and WASP-77A b (log φp = 13.40). The various contributions have been calculated using the simulations presented in Sect. 2. The first two objects fall in the range of low gravity planets, the third one is an intermediate gravity planet, while the last one is a genuine high gravity planet. The advection contribution is not shown for WASP-77A b, since it is of the same order of magnitude the order of the numerical approximation of the code.

5 Summary and conclusions

In this Paper, we make use of the ATES 1D photo-ionization hydrodynamics code (Caldiroli et al. 2021) to carry out a systematic assessment of the combined effects of planetary gravity and stellar irradiation upon atmospheric escape of primordial (hydrogen-helium) atmospheres, with a particular focus on the limits of validity of the so called energy-limited approximation (Eq. (1)). To investigate as broad as possible a parameter space, we simulated a set of 16 known nearby sub-Neptunes and hot-Jupiters (Table 1) over a broad range of irradiation levels. The target planets have specific gravitational potential energy between 12.15 ≲ log ϕ_p ≲ 13.56, and are selected to have radii below 0.15R_J, so as to reasonably fulfill the assumption of an extended hydrogen-helium atmosphere (Lopez & Fortney 2014). For each, we run ATES for input values of stellar XUV flux at the orbital distance (F_XUV) varying between 3 and 4 dex, in intervals of 0.5 dex. Our main results can be summarized as follows:

• Borrowing from cosmological re-ionization theory, we define the atmospheric mean excess energy as the ratio between the photo-heating rate and the photo-ionization rate per ion, and postulate that, for an atmospheric particle to be able to escape the system efficiently under the effect of photo-heating, the volume-averaged mean excess energy must exceed the particle binding energy. This simple criterion yields a (flux-independent) limiting ${\varphi }_{\text{p}}^{\text{thr}}$ below which energy-limited outflows may be attained: $\log {\varphi }_{\text{p}}\gtrsim \log {\varphi }_{\text{p}}^{\text{thr}}\approx [12.9-13.2]$ (in erg g⁻¹, see Fig. 2). This range agrees remarkably well with the phenomenological threshold that was previously identified in the literature (Salz et al. 2016b), and cements the notion that the energy-limited formula (i.e., Eq. (1) with η ≃ 0.3−1) severely overestimates the mass outflow rates for gas giants.

• Whether planets with ${\varphi }_{\text{p}}\lesssim {\varphi }_{\text{p}}^{\text{thr}}$ can be expected to drive energy-limited outflows mainly depends on the XUV irradiation level. As also shown by previous work (Owen & Alvarez 2016), the transition from energy-limited to a regime where cooling is dominated by Lyα is rather slow and occurs over a broad range of fluxes (F_XUV ≃ 10^3.5–5 erg cm⁻² s⁻¹), with higher gravity planets becoming radiation-dominated at lower fluxes.

• The introduction of an atmospheric mean excess energy translates into a novel definition of the outflow launching site, that is, the location where the mean excess velocity acquired by an ion by thermalization of a photo-electron exceeds the local escape velocity (Fig. 3b). Crucially, most of the stellar energy absorption occurs above this height in the case of low gravity planets, implying that the evaporation efficiency is largely independent of ϕ_p for those planets. Conversely, in the case of high-gravity planets only a small fraction of the energy is absorbed at heights where the ion excess velocity is sufficient high to warrant escape. As a result, the overall evaporation efficiency declines sharply as ϕ_p increases.

• To further elucidate the dependence of evaporation efficiency (defined in Eq. (1)) on irradiation we examine the relative contribution of adiabatic vs. advective and radiative cooling to the overall energy budget (Fig. 6). Energy-limited escape (≃100% efficiency) only occurs for planets below ${\varphi }_{\text{p}}^{\text{thr}}$ at the lowest irradiation, that is when advective and radiative losses are negligible. As F_XUV increases above ≳ 10^3–35 erg cm⁻² s⁻¹ the evaporation efficiency starts to decline, by up to an order of magnitude. This trend stems from a progressive increase in the relative contribution of advective cooling, followed by a sharp surge in Lyα losses (Fig. 6).

• In spite of the much lower absolute values (a few per cent or lower) compared with lower gravity planets, the evaporation efficiency of planets above ${\varphi }_{\text{p}}^{\text{thr}}$ increases with irradiation (Fig. 5). These systems are characterized by sharp ionization fronts and (nearly) fully ionized outflows, which implies a saturation of the Lyα cooling channel. As a result, the fractional contribution of adiabatic cooling (albeit low) increases with flux at a faster pace than Lyα.

In closing, these results allow for the rapid, physically motivated characterization of the expected atmospheric evaporation efficiency η for a given planet, thus enabling the community to move past the “one-size-fits-all” approach. An analytical approximation of the best-fitting effective efficiency (as a function of the system reduced gravitational potential, irradiation and planetary density) is given in the Appendix.

Notwithstanding the important limitation that our results only hold for atomic H + He atmospheres, they can be applied, for example, to perform systematic studies that investigate the role of photoevaporation-driven mass loss in shaping the observed distribution of planetary masses and radii. Additionally, they can be employed to compare the measured spectro-scopic transit signatures (or the lack thereof) against realistic mass outflow rates; and to efficiently select prime targets for future transit spectroscopy campaigns.

Acknowledgements

We thank the reviewer for their suggestions which have significantly improved the content and clarity of the manuscript.

Appendix A Planet-dependent evaporation efficiency

Based on the results presented in this work, we aim to provide an analytical approximation for the evaporation efficiency η_eff as a function of the relevant parameters at play. The ensuing best-fit efficiency is meant to replace the “one-size-fits-all” value that is routinely adopted within the framework of the energy-limited escape formalism, in the sense that it allows to employ the expression below regardless of whether the mass outflow rate is indeed energy-limited (according to the common definition): $$\dot{M}={\eta }_{\text{eff}}\frac{3F_{\text{XUV}}}{4GK{\rho }_{\text{P}}},$$(A.1)

where F_XUV is the XUV flux at the (average) orbital distance, ρ_p is the mean planetary mass density and the factor K accounts for the host star tidal forces⁵.

For a direct application, the interested reader is directed to ATES online repository⁶ where a dedicated script enables the straightforward estimate of the best-fitting η_eff (and ensuing Ṁ) for a given planetary system. For completeness, we provide a quantitative description of the best-fitting analytical approximation in § A.1.

Appendix A.1 Analytical approximation

Noting that the energy-limited formalism implies a linear scaling of the mass loss rate with the irradiation to planetary density ratio, as well as the inverse of the K-factor, we aim to derive a best-fitting approximation as a function of those two variables. Concurrently, we wish to incorporate the dependence on planetary gravitational potential energy ϕ_p = GM_p/R_p highlighted by our analysis (§ 3). Thus, for a given planetary system with known (i) planetary radius R_p (defined as the optical transit radius); (ii) planetary mass, M_p; (iii) host stellar massM_⋆; (iv) stellar XUV flux F_XUV (see § A.2 below) measured at the (v) orbital distance a, we define the following quantities:

• ф_red ≡ Κф_p [erg g⁻¹]

• ф_red,0 = 10^13.22 [erg g⁻¹]

• ${\tilde{F}}_2$ ≡ F_XUV/ρ_p in units of 10² erg cm² s⁻¹ g⁻¹.

The functional dependence of η_eff on ϕ_red and ${\tilde{F}}_2$ is chosen to reproduce the behavior of Figure 2, that is: at low-gravity ( ≲ 13.2, in log cgs units) η_eff is approximately constant ≃ η₀) for a fixed value of ${\tilde{F}}_2$. At at high-gravity ( ≳ 13.2, in log cgs units), η_eff decreases sharply with increasing ${\tilde{F}}_2$, where the actual rate of decrease is in turn a function of ${\tilde{F}}_2$. The low- and high-gravity functional forms are joined through an appropriate weighting function, that is a sigmoid curve in the log ϕ_red − log η_eff plane.

The resulting functional shape is as follows: $$\text{log}{\eta }_{\text{eff}}=A{\varphi }_{\text{red}}^{\alpha }\sigma +{\eta }_{\text{o}}\left(1-\sigma \right)$$(A.2)

with $$\sigma ={\left[1+{\left(\frac{{\varphi }_{\text{red}}}{{\varphi }_{\text{red,o}}}\right)}^{\beta }\right]}^{-1},\beta <0,$$(A.3)

where A, α, η₀, and ß are expressed as power-law combinations of ${\tilde{F}}_2$ and log $K=1-\frac{3}{2\xi }+\frac{1}{2{\xi }^3}$, as appropriate, through Equations A.4–A.7. The numerical coefficients were evaluated through a least-squares fitting Matlab routine⁷ making use of 81 η_eff values resulting from the simulations discussed in § 3 and § 4 of the Paper, yielding: $$A=1.682{\left(\log {\tilde{F}}_2\right)}^{0.2802}-5.488,$$(A.4) $$\alpha =0.02489{\tilde{F}}_2^{-0.0860}-0.01007{\tilde{F}}_2^{-0.9543},$$(A.5) $${\eta }_0=-0.03973{\left(\log {\tilde{F}}_2\right)}^{2.173}-0.01359,$$(A.6) $$\beta =-0.01799{\tilde{F}}_2^{0.1723}-3.3875{\tilde{F}}_2^{0.0140}.$$(A.7)

Fig. A.1 Analytic approximation of the function ηeff(FXUV/ρp, Kϕp) as given by Eqs. (A.2)–(A.7). The data points are color-coded to reflect the intensity of the XUV irradiation, according to the same color scheme as in Fig. 2 where FXUV varies between 102 (deepest blue) and 106 (deepest red) erg s−1 cm−2. The thick red line represents the convergence limit of the ATES code described in Sect. 6.1 of Caldiroli et al. (2021).

Figure A.1 illustrates the full functional dependence of the best-fitting η_eff on F_XUV/ρ_p and Kϕ; its projection on the F_XUV/ρ_p − η_eff plane is shown in Figure A.2. The maximum discrepancy between the Ṁ values obtained by ATES (and discussed in § 3 and § 4 of the paper) and the best-fit values obtained through Equation A.1 is ≲ 0.16 dex; the mean deviation is ≃ 0.075 dex. We remark that the above formulas are formally valid within the simulated range of parameters, that is for 10² ≲ F_xuv/ρ_p ≲ 10⁶ and 10^12.17 ≲ Kϕ_p ≲ 10^13.29 (in cgs units).

Appendix A.2 F_XUV estimate

The mass loss rates based upon which the above fit is derived are obtained by ATES under the assumption of a power-lawlike photoionizing spectrum. Specifically, ATES assumes a piece-wise power law of the form F ∝ E⁻¹, both in the EUV ([13.6 – 124] eV) and in the X-ray band ([0.124 – 12.4] keV). In each band, the spectrum is normalized to the (user-provided) input EUV and X-ray luminosity values. As discussed in Section 6 of Caldiroli et al. (2021), this choice is motivated by the fact that the steady-state mass loss rates mainly depend on the total number of photoionizing photons, rather than the chosen stellar spectral shape.

Fig. A.2 Projection of the analytical approximation illustrated in Fig. A.1 on the FXUV/ρp − ηeff plane.

References

All Tables

All Figures

	Fig. 1 Simulated outflow properties for the Neptune-like planet GJ 3470 b (Panel a) and for the hot Jupiter WASP-77A b (panel b). FXUV varies between 102 (deepest blue) and 106 erg cm−2 s−1 (deepest red), in intervals of 0.5 dex. The diamond symbols mark the location of the τ = 1 surface at 13.6 eV.
In the text

Fig. 2 Outflow effective efficiency is plotted against the planet specific gravitational potential for the 16 planets considered in this work (see Table 1 for a list); ηeff is defined as the ratio between the ATES-estimated mass outflow rate and the energy-limited value given by Eq. (1) for η = 1. To explore the effects of varied irradiation levels on Ṁ, (81) simulations were carried out for a broader range of FXUV for a sub-sample of eight targets. The results are represented by the colored circles, where the adopted color scheme is the same as in Figs. 1a and b.

In the text

Fig. 3 H I mean excess energy defined in Eq. (2) (panel a) as a function of planetary distance for different systems and Second power of the ratio between the excess and escape velocity defined in Eqs. (3) and (4) (panel b); r and RRL are expressed in units of the planet radius, such that the planetary distance on the x axis is given as a fraction of the Roche lobe radius RRL. The dashed line represents the threshold for wind launching, i.e., where (υxs/υesc)2 = 1. In both panels, the diamond symbols indicate the location of the τ = 1 surface at 13.6 eV.

In the text

Fig. 4 Simulated mass loss rate Ṁ for an artificial planet with Mp = 1.3 × 1030 g and Rp = 1010 cm, subject to different levels of stellar XUV irradiation (panel a) and simulated mass outflow rates for eight planets subject to different levels of stellar XUV irradiation, for a total of 81 simulations (panel b). In panel a, the dashed red line represents the results of the simulations carried out by Murray-Clay et al. (2009), under the assumption of a pure hydrogen atmosphere and mono-energetic irradiation field at 20 eV. The blue line represents the results obtained by ATES under the same assumptions. The green line represents the ATES results for a pure hydrogen atmosphere subject to a power-law irradiation spectrum with slope v−1. The black line refers to a primordial hydrogen-helium atmosphere with a power-law spectrum, that is, the standard setup for ATES. The adopted color scheme in panel b is the same as shown in Figs. 1a, b and 2. In both panels, the dotted lines represent linear proportionality between mass loss rate and flux.

In the text

	Fig. 5 Simulated, K-reduced mass outflow rate, shown as a function of the irradiation to planet mass density ratio, for the same planets shown in Fig. 4b. Energy-limited escape yields KṀ = FXUV/ρp (i.e., η = 1 in Eq. (1)). The adopted color scheme is the same as shown in Figs. 1a, b and 2.
In the text

Fig. 6 Volume integrals of the various cooling processes of our model as functions of the impinging XUV flux. From top-left, moving clockwise, the four panels represent different planets with increasing φp, from : GJ 3470 b (log φp = 12.30), HD 149026 b (log φp = 12.84), WASP-80 b (log φp = 13.02) and WASP-77A b (log φp = 13.40). The various contributions have been calculated using the simulations presented in Sect. 2. The first two objects fall in the range of low gravity planets, the third one is an intermediate gravity planet, while the last one is a genuine high gravity planet. The advection contribution is not shown for WASP-77A b, since it is of the same order of magnitude the order of the numerical approximation of the code.

In the text

Fig. A.1 Analytic approximation of the function ηeff(FXUV/ρp, Kϕp) as given by Eqs. (A.2)–(A.7). The data points are color-coded to reflect the intensity of the XUV irradiation, according to the same color scheme as in Fig. 2 where FXUV varies between 102 (deepest blue) and 106 (deepest red) erg s−1 cm−2. The thick red line represents the convergence limit of the ATES code described in Sect. 6.1 of Caldiroli et al. (2021).

In the text

	Fig. A.2 Projection of the analytical approximation illustrated in Fig. A.1 on the FXUV/ρp − ηeff plane.
In the text