© The Authors 2023

1 Introduction

The fact that planets with water-rich atmospheres may be common in our galaxy is a remarkable possibility, as such planets are not known to occur in our vicinity. With varying degrees of uncertainty, this idea is supported by the following points: i) planetary formation and evolution theories predicting that significant water is accreted during the embryo formation if it takes place beyond the snowline, or at a later time and possibly closer to the host star while the planet is being bombarded by icy planetesimals, or from chemical reaction between a thick primordial envelope and an oxidizing magma ocean (Venturini et al. 2020; Bitsch et al. 2021; Kimura & Ikoma 2022); ii) the constraints inferred from interior structure models on the water mass fraction for individual planets with measured mass, M_p, and radius, R_p (Valencia et al. 2010; Delrez et al. 2021); iii) a characteristic dispersion in the M_p–R_p data, suggesting a transition between rocky planets with essentially no envelopes for M_p/M_⊕≲5 and planets with a broad range of compositions for M_p/M_⊕≳5 (Otegi et al. 2020); iv) a peak in the occurrence of planets orbiting M dwarfs for bulk densities consistent with high water abundances and possibly inconsistent with planets that lack an atmosphere or that have one dominated by H₂/He (Luque & Pallé 2022).

The state of the water in the planets also stands as an open debate. Steam has been identified with transmission and emission spectroscopy in a number of exoplanets (Evans et al. 2017; Benneke et al. 2019; Tsiaras et al. 2019). In addition, water might occur in the solid, liquid, and plasma phases, or in a supercritical state (Nettelmann et al. 2010; Madhusudhan et al. 2020; Mousis et al. 2020). Significant attention will be devoted to water as a key tracer in exoplanet atmospheres in the coming years, notably because the JWST and ARIEL (Tinetti et al. 2018) telescopes have the capacity to probe the molecule in the infrared. The upcoming observations might firmly establish the existence and general properties of water-rich atmospheres.

Further insights into the water content of a strongly irradiated atmosphere can be gained by investigating its long-term stability. The mass in the atmosphere of a small planet of known Mp and Rp may be orders of magnitude larger if the gas has a high molecular weight than if it is composed of H₂/He. This means that a high-molecular weight atmosphere is, in principle, more stable even when a vigorous outflow driven by stellar irradiation is established. For example, García Muñoz et al. (2020, 2021) predicted that the sub-Neptune π Men c in its current orbit might lose its atmosphere in a few Myrs if this is composed of H₂/He, which is much less than the estimated ~5 Gyr age of the system, but that such a catastrophic fate would take longer by ×10³ if the atmospheric composition is dominated by H₂O. These arguments, along with supplementary arguments based on the detection of C II and the non-detection of H I Lyman-α in the upper atmosphere of the planet, led the authors to propose that π Men c has a high-molecular weight and (possibly) a water-rich atmosphere.

As such analyses rely on model-predicted mass loss rates, it is important to assess the completeness of the models, particularly in their treatment of the collisional-radiative processes that occur in the low-pressure layers where the outflow is established (García Muñoz et al. 2020, 2021; Shaikhislamov et al. 2020). It is worth noting that the young Earth and Venus may have had water-rich atmospheres, a possibility that has been explored to better understand their histories (Zahnle et al. 1988; Chassefière 1996; Lichtenegger et al. 2016). It is thus no surprise that a new generation of models continue to explore the outflow of water-rich atmospheres (Kurosaki et al. 2014; Luger & Barnes 2015; Guo 2019; Johnstone 2020; García Muñoz et al. 2020, 2021; Yoshida et al. 2022) and the relevance of related ideas to exoplanets.

We are currently revisiting the formulation of the collisional-radiative processes that occur in the upper atmospheres of exoplanets (García Muñoz & Schneider 2019). We believe that past modeling may not have captured all the complexity implicit in a high-temperature, multi-component plasma that transitions from optically thick to optically thin conditions. The ultimate goals of this effort are to obtain accurate mass loss rates for a range of atmospheric compositions and stellar hosts, and set up a physics-based framework for the interpretation of the spectro-scopic searches that have been conducted at some exoplanets. In particular, the current work addresses the fundamental issue of tracking the formation and slowing down of photoelectrons in water-rich atmospheres, as this partly dictates the fraction of the stellar irradiation that turns into heating and drives the outflow at those planets. That information is often quantified in the form of so-called heating efficiencies and production yields.

Calculations of heating efficiencies and production yields for general astrophysical application are available, for example, for H₂–H mixtures and for pure O and O₂ (e.g. Green et al. 1977; Victor et al. 1994; Dalgarno et al. 1999; García Muñoz 2023, and references therein). Some calculations have also been reported for exoplanets with hydrogen-dominated atmospheres (Cecchi-Pestellini et al. 2006, 2009; Shematovich et al. 2014; Guo & Ben-Jaffel 2016). We have however not found similar calculations for H₂O nor O-H mixtures. There is a long history of simulations of the slowing down of photoelectrons in the H₂O-dominated comae of comets, but they emphasize the excitation of neutrals and ions that can be remotely sensed rather than the deposition of energy; furthermore, they are mostly concerned with conditions of low fractional ionization (Ashihara 1978; Bhardwaj 2003). For completeness, it should be mentioned that some recent works related to Earth-like (Tian et al. 2008) and water-rich atmospheres (Johnstone 2020; Nakayama et al. 2022) have treated the problem of heating and ionization by photoelectrons; however, it is difficult to judge the accuracy of their approximate treatments and the significance of the phenomenon in the conditions of their simulations, as compared to more general ones.

The above arguments motivate us to review the problem of photoelectrons in water-rich atmospheres in a way that is reliable and easy to implement into atmospheric models. We focus here on O–H mixtures, namely, it is assumed that the H₂O has dissociated, but will address the more general problem of atmospheres that contain this and other molecules in undissociated form in future work. In general, this simplification is not critical, as water dissociates rapidly at close-in exoplanets when they orbit stars that are not too cool (Guo 2019; Johnstone 2020; García Muñoz et al. 2020, 2021). Even when the star is cool (and its far-ultraviolet emission weak), a calculation based on an O-H mixture provides valuable insight into what happens in the region where the XUV photons (wavelengths shorter than the Lyman continuum threshold) are deposited. We intentionally avoid approximate methods (Peterson et al. 1978) to estimate the heating efficiencies and production yields of mixtures based on the properties of their individual constituents to keep a better control of the relevant physics and because such approximations become inaccurate for the conditions of small-to-moderate photoelectron energies that are relevant here.

The paper is organized as follows. Section 2 offers a discussion of the heating efficiency and production yields in O–H mixtures. Section 3 explores the heating and ionization of the postulated water-rich atmospheres of the exoplanets π Men c and TRAPPIST-1 b, and estimates the lifetime of O(¹ D) in them. Lastly, Sect. 4 summarizes the main findings and lists some possible next steps for future studies.

2 Slowing down of photoelectrons in the O–H mixture

Photoelectrons are produced in the interaction of ionizing stellar photons with the atmospheric gas, in the process acquiring kinetic energies of up to hundreds of eV. The specifics of photoelectron formation are dictated by the properties of the local atmosphere and the stellar radiation that reaches it. Once released into the atmosphere, the nascent photoelectrons collide repeatedly with the thermal electrons and heavy particles (atoms in this work), exchanging at each collision some energy with them until the electrons are ultimately thermalized. The interactions can be elastic: $$e\prime +e_{\text{th}},X\to e\prime \prime +e_{\text{th}},X,$$

if the exchange with the thermal electrons, e_th, or the heavy particles, X, involves only kinetic energies – or inelastic: $$e\prime +X\to e\prime \prime +X^{*},e\prime \prime +e\prime \prime \prime +X^{+},$$

if the heavy particle becomes excited or ionized (by ejecting one electron in the example). Superelastic collisions, the reverse of inelastic collisions, also occur, especially for high X^* or X⁺ densities, but their significance is often minor and easier to assess on a case-by-case basis; e' refers to the incident photoelectron and e'' and e''' refer to the post-collision photoelectrons. We use the term photoelectrons as equivalent to non-thermal electrons to describe the population of electrons that cannot be ascribed to the dominant single-temperature Maxwellian distribution of velocities. We note that other works use the term photoelectrons to designate the primary electrons produced directly by pho-toionization and keep the term secondary electrons to refer to the electrons that form as the primary electrons slow down and eject additional electrons in collisions with the heavy particles. In our use of these terms, e^', e^'', and e^''' qualify as photoelectrons as long as their energies are sufficiently larger than the mean kinetic energy of the thermal electrons.

In the O–H gas that we consider, X may be any of the H, O, or O⁺ atoms. Because the electron-to-proton mass is ≪1, elastic collisions with H⁺ have a negligible effect compared to elastic electron-electron collisions and are omitted, and we also omit collisions of photoelectrons with O²⁺ and higher ions as their densities are expected to be low. In elastic collisions, the energy lost by the photoelectron is transferred to the translational energy of the thermal electrons and heavy particles, thereby heating the gas. In inelastic collisions, some of the photoelectron energy is expended to promote or eject a bound electron initially in the heavy particle. This expense represents in principle a sink for the thermal budget of the gas because it does not contribute to its translational energy. Some of that excitation or ionization energy may, however, eventually convert into heat if the heavy particle undergoes a non-radiative collision.

To simulate the slowing down of the photoelectrons in an OH mixture, we used the Monte Carlo (MC) model developed in García Muñoz (2023). The model has been expanded (the original version included only H and He atoms; see Appendix A) to take into account collisions of the photoelectrons with O and O⁺. Our numerical experiments revealed that the role of collisions with O⁺ is minor for the conditions of moderate energy and fractional ionization of interest here, in agreement with past findings for a pure O gas (Kozma & Fransson 1992; Victor et al. 1994), yet the corresponding collisional channels are kept for completeness. For reference, Fig. 1 shows the excitation and ionization cross sections for collisions with H and O. The former have been scaled by ×2 to clearly reveal the relative contributions in a O–H gas resulting from water dissociation.

Many of the MC model outcomes can be traced to the col-lisional properties of the atoms in the mixture, so it is worth noting some of them. Figure 1 shows that the inelastic (excitation + ionization) cross sections of H and O are comparable in magnitude, especially at energies of ≳100 eV. Therefore, in an O–H gas resulting from water dissociation, the photoelectrons will interact with either type of atom with comparable probabilities. For energies of ≳30 eV, the cross section for H excitation is larger by about ×2 than the cross section for H ionization. In contrast, over the same range of energies, the cross sections for O ionization are nearly an order of magnitude larger than those for O excitation. As a consequence ionization, rather than excitation, should dominate the loss of energy for the photoelectrons interacting with the O atoms.

We carried out a systematic investigation of the slowing down of photoelectrons in O–H mixtures. In the calculations, each mixture contains the atoms H, H⁺, O, and O⁺, plus thermal electrons, e_th. The number densities of hydrogen and oxygen nuclei are [H_n]=[H]+[H+] and [O_n]=[O]+[O⁺], respectively, and for the mixture [X_n]=[H_n]+[O_n]. We assume at first that H, O, and O⁺ are all in their ground levels, but we relax this assumption in the steps detailed in Sect. 3.2 to explore the lifetime of the metastable O(¹ D). We assume charge neutrality, [H⁺]+[O⁺]=[e_th], and that [H⁺]/[H]=[O⁺]/[O]. The latter is expected to be a good approximation because charge exchange H⁺+O↔H+O⁺ in the atmosphere will tend to rapidly establish the statistical balance ([H⁺]/[H])/([O⁺]/[O])=9/8. This simplification reduces the free parameters in the systematic investigation. With the above, the relative proportions of the mixture components are fully specified by the ratio ƒ_On/Hn=[O_n]/[H_n] and the fractional ionization x_e=[e_th][X_n]. We note that ƒ_On/Hn =0 for a pure H gas, =1/2 for the O–H gas resulting from water dissociation, and ∞ for a pure O gas. Also, x_e=0 for a neutral gas, and =1, when [H]=[O]=0. After some algebra, the fraction of each atom in the mixture is [H]/[X_n]=(1 – x_e)/(1 + ƒ_On/Hn), [H⁺]/[[X_n]=x_e/(l + ƒ_On/Hn), [O]/[X_n]=ƒ_On/Hn(1 – x_e)/(1 + ƒ_On/Hn), and [O⁺]/[X_n]=ƒOn/Hnx_e/(1 + ƒ_On/Hn). We adopt [X_n]=10¹⁰ cm⁻³ in the MC model, although the actual choice has a minor impact on the calculations because, except for the electron-electron cross sections, all the other factors in the formulation are independent of the absolute densities. Calculations were performed for ƒ_On/Hn =0, 1/6, 1/2, 3/2, ∞, and for log₁₀(x_e)=−6, −5, −4, −3, −2, and −1, for incident photoelectron energies below 10 000 eV.

As it slows down, a photoelectron of energy, E₀, loses a fraction, E_elas, in elastic collisions and another fraction, Σj E_inel,j, in inelastic collisions, where the summation extends over all the availa_able inelastic channels. From energy conservation, E₀=E_elas+∑_j E_inel,j,¹. It is usually accepted that all of E_elas contributes to heating the gas and that some of ∑_jE_inel,j may also contribute locally. Using E_h to denote the part of E₀ that goes into heating, it follows that E₀≥E_h≥E_elas. The ratio E_h/E₀ is a heating efficiency, although this denomination may be confusing as it is used in the literature to describe a variety of related but non-identical concepts.

Figure 2 summarizes some of the calculations for an O–H mixture with ƒ_On/Hn=1/2, as well as for pure H and pure O. The top panel shows E_elas/E₀ and the bottom one shows the normalized yield Φ_i/E₀. The ionization yield, Φ_i, represents the number of ions produced per injected photoelectron of energy, E₀, regardless whether they originated from H, O, or O⁺ and of their charge. Because the MC model considers double ionization for O(³P), which leads to O²⁺(³P) plus two electrons, Φ_i is slightly larger than the number of secondary electrons produced.

We do not discuss the specifics of our calculations for pure H (see e.g., García Muñoz 2023, for that) or for pure O. For the latter, however, we note that they are consistent with what has been reported in past works (see e.g. Dalgarno & Lejeune 1971; Green et al. 1977; Slinker et al. 1988; Kozma & Fransson 1992; Victor et al. 1994). For example, we estimate from Fig. 14 of Dalgarno & Lejeune (1971) that the amount of energy transferred elastically to the thermal electrons for E₀=50 eV is about 25, 13, 8, 4, and 2.5 eV for log₁₀(x_e)=−2, −3, −4, −5, and −6, respectively. For the same conditions, we obtained the following values: 25.5, 14.9, 8.4, 4.6, and 3.4 eV.

It is seen in Fig. 2 (top) that E_elas/E₀ depends on the gas composition and also on the photoelectron initial energy, E₀, and the fractional ionization of the mixture, x_e. For E₀< 10.2 eV, E_elas/E₀≡1 for pure H, but <1 for mixtures that contain O, especially for low x_e. The reason is that the O atom has multiple excitation channels open between 10.2 and 1.96 eV (the threshold for O(¹ D) excitation) that remain competitive against electron-electron collisions for small-to-moderate x_e. For E₀>100 eV and x_e~10⁻³–10⁻², E_elas/E₀ is typically higher for an O-H mixture with ƒ_On/Hn=1/2 than for pure H. The trend does, however, reverse for both lower and higher fractional ionizations. For pure H and E₀≳50eV, the heating efficiency E_elas/E₀ becomes nearly independent of the fractional ionization for x_e≲10⁻⁴. In contrast, E_elas/E₀ drops to ~0.07 (not shown in Fig. 2) in the O–H mixture with ƒ_On/Hn = 1/2 and x_e=10⁻⁶ because for such small fractional ionizations excitation into O(¹ D) remains efficient at extracting energy from the slow photoelectrons. It is also apparent in Fig. 2 (bottom) that adding O enhances Φ_i and therefore the number of newly created ions with respect to the case of pure H. The reason for this is that at energies ≳30 eV, the collisions of photoelectrons with O tend to favor its ionization rather than its excitation.

Figure 3 is specific to an O–H mixture with ƒ_On/Hn = 1/2, and shows the ratio E_inel,j/E₀ for the various excitation and ionization channels – except ionization from either O(³P) or O+(⁴S^o) to O²⁺, which remain minor. We note that excitation is typically more significant than ionization for H, whereas the reverse occurs for O. Amongst the O excitation channels, excitation into the metastable O(¹ D) plus a smaller contribution to O(¹S) stand out as important, especially at low energies. Significant amounts of energy go into the simultaneous ionization and excitation of O. O⁺ (⁴P) is likely to radiate rapidly to the ground level O⁺ (⁴S^o). In contrast, the metastables O+(²D^o) and O+(²P^o) may return some of their excitation energy to the ground level ion in collisions with the other gases in the atmosphere. Interestingly, most of the newly created ions are produced through O ioniza-tion rather than through H ionization, even though the ratio of number densities in these calculations is [O]/[H]=1/2.

Fig. 1 Cross sections for excitation (H*, O*), ionization (single, H+, O+, and double, O2+) and total inelastic (=excitation + ionization) collisions of photoelectrons with ground-level H and O atoms (LG90 set, see Appendix A). The H cross sections have been multiplied by ×2 to better reveal the significance of collisions with both types of atoms in an O–H mixture with ƒOn/Hn = 1/2. We also show the adopted electron-electron pseudo-cross sections multiplied by ×10−3.

Fig. 2 Some calculated properties from the slowing down of photoelectrons. Top: fraction of the photoelectron initial energy E0 that is transferred to the gas in elastic collisions. Bottom: number of secondary ions created by a photoelectron of initial energy E0, normalized by E0. The three color-style combinations refer to pure H (black, long dashed), pure O (black, dotted), and an O–H mixture with ƒOn/Hn = 1/2 (blue, solid). For each mixture, the four curves that are displayed (distinguishable by their thickness) refer to fractional ionizations log10(xe) = −1 (thickest), −2, −3, and −4 (thinnest). Although not shown for the sake of clarity, the Eelas/E0 ratio for E0≳100eV for pure H becomes essentially independent of fractional ionization at xe<10−4. In contrast, Eelas/E0 continues to drop with xe for the other O–H mixtures. For example, for E0>200eV and xe~10−6, Eelas/E0~0.05, and 0.07 for pure O and the O–H mixture with ƒOn/Hn=1/2, respectively. Calculations were done with the LG90 set of cross sections.

3 Application to π Men c and TRAPPIST-1 b

The heating and ionization of an atmosphere are partly dictated by the interaction of the photoelectrons with the gas while they slow down to thermal energies. That information is encoded in E_elas/E₀ and Φ_i. A more complete description requires, in addition, tracking the fate of the excited and ionized levels produced by the photoelectrons. For example, the energy E_h that goes into heat will exceed E_elas if the population of excited levels produced by the photoelectrons is rapidly quenched (collisionally deexcited) to a lower-energy level. The problem is complex and non-local, and involves solving self-consistently the equations of mass, momentum and energy conservation, and radiative transfer in the atmosphere.

Here, we focus on the simpler problem of estimating the heating and ionization that occurs locally. These quantities depend on the local density, composition and fractional ioniza-tion of the gas in ways that are described below. We devote the bulk of our study to the exoplanets π Men c (Gandolfi et al. 2018; Huang et al. 2018) and TRAPPIST-1 b (Gillon et al. 2017), whose atmospheres might prove to be water-rich, according to some analyses (García Muñoz et al. 2021; Agol et al. 2021).

For π Men c (M_p/M_⊕=4.52+0.81, R_p/R_⊕=2.06+0.03; orbital distance a = 0.067 AU; Gandolfi et al. 2018), we take the number densities in the upper atmosphere from a model published by García Muñoz et al. (2021) that assumes ƒ_On/Hn≈1/2 (their case 6). TRAPPIST-1 b (M_p/M_⊕=1.374+0.069, $R_{\text{p}}/R_{\otimes }={1.116}_{-0.012}^{+0.014}$; a=(1.154±0.010)×10⁻² AU; Agol et al. 2021) is part of a system of seven known planets, some of which orbit in the habitable zone of their ultracool dwarf hots star. For TRAPPIST-1 b, we used the same photochemical-hydrodynamical model as in García Muñoz et al. (2021) to predict the number densities. We prescribe at the bottom boundary of the model (pressure of 1 dyn cm⁻²) that the atmosphere is pure H₂O and let the model evolve its composition at higher altitudes. We adopted the TRAPPIST-1 spectrum published by Wilson et al. (2021), which shows that the far-ultraviolet emission of the star (a wavelength range usually key in the photodissociation of molecules) is rather faint compared to the XUV emission.

Figure 4 (top) shows the number densities of H, H⁺, O, O⁺ and e_th for both planets, which we adopt throughout our analysis as shown there. At pressures of ≳10⁻³ dyn cm⁻², the chemistry of TRAPPIST-1 b becomes particularly interesting because the weak photodissociation at far-ultraviolet wavelengths facilitates the survival of molecules. Indeed, H₂O and O₂ are abundant in that region amongst the neutrals, with H₃O⁺ and ${\text{O}}_2^{+}$ amongst the ions. At lower pressures, the H₂O molecule is no longer shielded by the gas column overhead and its abundance drops rapidly, and so does the abundance of the other molecules. For TRAPPIST-1 b, we focus our analysis on pressures ≲ 10⁻³ dyn cm⁻², for which the chemistry remains mainly atomic. We recall that the model in García Muñoz et al. (2021) does not take into account the non-thermal effects described here or, equivalently, that it assumes E_h/E₀≡1 and Φ_i≡0. Our analysis below is therefore not self-consistent because it does not allow for feedbacks that might modify the atmosphere. Even with these caveats, which will be fixed in follow-up work, the analysis clearly demonstrates and partially quantifies the importance of non-thermal effects in water-rich atmospheres.

We calculated the production rates for primary photoelec-trons from H + hc/λ → H⁺ + e⁻ and O + hc/λ → O⁺ + e⁻. We note that hc is the product of Planck's constant and the speed of light, and λ is the wavelength of the ionizing photon. They are shown in Fig. 4 (middle) as curves I and II, while curves III show the total production rates of primary photoelectrons (=I+II). The O atom contributes most of the primary photoelectrons, even though its abundance is about half that of H. The reason for this is that the cross section for O photoionization is generally larger than the cross section for H photoionization. The disparity in the production rates between the O and H atoms is exacerbated in the case of TRAPPIST-1 b by the large fraction of XUV photons emitted by its host star at the shortest wavelengths.

Each primary photoelectron of energy E₀ (≈hc/λ–13.6eV for both H and O, where for consistency with the photochemical-hydrodynamical model it is assumed that O photoionization results in ground-level O⁺) produces a number of secondary electrons that depends on E₀ and the local composition and fractional ionization of the atmosphere. We used the number densities of H, H⁺, O, O⁺, and eth from Fig. 4, plus the energy spectrum of the primary photoelectrons, to calculate the yield Φ_i and, ultimately, the net ionization rate at each altitude following conventional methods (García Muñoz 2023). Figure 5 shows the energy spectrum of the primary photoelectrons in the atmosphere of π Men c over a range of pressures. At low pressures, photoionization of the gas by photons of wavelengths near the Lyman continuum threshold results mainly in low-energy photoelectrons. Deeper in the atmosphere, where only very energetic photons penetrate, the nascent photoelectrons become increasingly energized.

3.1 Ionization and heating rates

The total production rates, that is, primary plus secondary electrons, are shown in Fig. 4 (middle) as curves IV. The main finding for π Men c is that the photoelectrons enhance the ionization rate by a factor of a few at p≳10⁻² dyn cm⁻². This region of the atmosphere is key because it defines where the stellar energy is deposited and the gas becomes accelerated to escape the planet. A change in the net ionization rate will surely have implications on these important aspects. For TRAPPIST-1 b, at p~2× 10⁻³ dyn cm⁻² the contribution of photoelectrons to the overall ionization rate is about 40%. It is justifiable to assume that as for π Men c, their relative contribution will increase towards the deeper atmosphere. The high abundance of molecules there (in particular H₂0) may, however, introduce specific features in the ionization balance that are worth investigating in the future.

For the conditions encountered in the hydrogen-dominated atmospheres of close-in exoplanets at pressures of ≲l dyn cm⁻², the excited levels of hydrogen will likely radiate away their energy, therefore not contributing to the local heating of the atmosphere (see e.g., García Muñoz 2023, for the specific conditions of the ultrahot Jupiter HAT-P-32 b)². This occurs even if the lines in the Lyman series, which connect with the ground level H(l s), are optically thick. Also, the excited level H(2s), which does not connect with H(l s) through an optically allowed transition, will typically radiate by conversion first into H(2p) via collisions with protons. Similar arguments should be valid for most of the excited levels of O and O⁺, except (possibly) for the metastables O(¹D), O(¹S), O⁺(²D), and O⁺(²P), which have long radiative lifetimes and may become quenched before they radiate.

When considering the contribution of the metastables to heating, two limits can be identified. In the no-quenching limit, which should be valid when collisions are rare and quenching is inefficient, only the kinetic energy transferred by the photoelectrons to the thermal electrons and heavy particles in elastic collisions contributes to heating, $$E_{\text{h}}^{\text{nq}}=E_{\text{elas}}.$$

Alternatively, in the full-quenching limit, which should apply when many deexcitation collisions occur before the atom radiates, $$\begin{array}{c}{E}_{\text{h}}^{\text{fq}}=E_{\text{elas}}+E_{\text{inel,O}\left({}^3\text{P}\right)\to \text{O}\left({}^1\text{D}\right)}+E_{\text{inel,O}\left({}^3\text{P}\right)\to \text{O}\left({}^1\text{S}\right)}\\ +E_{\text{inel,O}\left({}^3\text{P}\right)\to {\text{O}}^{+}\left({}^2\text{D}\right)}\frac{\left(16.93-13.61\right)}{16.93}\\ +E_{\text{inel,O}\left({}^3\text{P}\right)\to {\text{O}}^{+}\left({}^{2\text{P}}\right)}\frac{\left(18.63-13.61\right)}{18.63}\\ +E_{{\text{inel,O}}^{+}\left({}^4\text{S}\right)\to {\text{O}}^{+}\left({}^2\text{D}\right)}+E_{{\text{inel,O}}^{+}\left({}^4\text{S}\right)\to {\text{O}}^{+}\left({}^2\text{P}\right)}.\end{array}$$

The latter assumes that when simultaneous ionization-excitation occurs, quenching converts into heat only the excitation energy above the ground level of the ion.

We calculated the heating rate in the atmosphere as follows: $$\mathcal{Q}={\displaystyle \int \left(\left[\text{H}\right]{\sigma }_{\text{H}}^{\text{pi}}\left(\lambda \right)+\left[\text{O}\right]{\sigma }_{\text{O}}^{\text{pi}}\left(\lambda \right)\right){ℱ}_{\star }\left(\lambda \right)\left(1-\lambda /{\lambda }_{\text{Ly,c}}\right)E_{\text{h}}/E_0\text{d}\lambda },$$

where σ^pi stands for the photoionization cross section, λ_Ly,c is the Lyman continuum threshold (to practical effects, the H and O atoms have the same ionization potential), and the stellar radiation flux ℱ_★(λ) is expressed in energy units. The results for $E_{\text{h}}=E_{\text{h}}^{\text{nq}}$ and $E_{\text{h}}=E_{\text{h}}^{\text{fq}}$ are shown in Fig. 4 (bottom) as two different curves. They should bracket the actual heating rate in the atmosphere because in general the metastables lose a part of their energy to radiation and another part to collisional deexcitation. For reference, the heating rate calculated for E_h/E₀≡l, as in the photochemical-hydrodynamical model of García Muñoz et al. (2021), is also shown.

One conclusion from this figure is that metastable quenching may enhance E_h over E_elas by minor amounts in the cases under study here. More importantly perhaps, the figure shows that neglecting photoelectrons may overestimate the heating rate by a factor of a few. A more accurate determination of this requires considering possible feedbacks in the atmosphere.

Fig. 3 For an O–H mixture, ƒOn/Hn=1/2 is the fraction of the photoelectron energy, E0, that goes into various inelastic channels specified by their initial and end levels. Then, H* and O* refer to the summation over all excited levels within the corresponding neutral atoms. For each panel, the four curves that are displayed (distinguishable by their thickness) refer to fractional ionizations log10(xe)=−1 (thickest), −2, −3, and −4 (thinnest). Calculations done with the LG90 set of cross sections.

Fig. 4 Some properties of the atmospheres of π Men c and TRAPPIST-1 b. Top: adopted number densities of H, H+, O, O+ and eth,. See text for details. For TRAPPIST-1 b, the densities of additional molecules (H2O, O2, H3O+, ${\text{O}}_2^{+}$) that remain relatively abundant up to 10−3 dyn cm−2 are also shown. Middle: production rates of primary and secondary ions. Bottom: heating rates in the no- and full-quenching limits (see text), and under the assumption that photoelectrons transfer all their kinetic energy into heat.

Fig. 5 Energy spectra of the nascent photoelectrons produced through photoionization of O and H atoms at different levels in the atmosphere of π Men c for ƒOn/Hn=l/2. The quoted numbers in the plot represent the pressure levels.

3.2 Metastables O(¹ D) and O(¹ S)

The two lowermost excited levels of atomic oxygen, O(¹D) and O(¹S), lie at only 1.96 and 4.18 eV above the ground-level O(³P). They are metastable, with total radiative lifetimes of 134 and 0.75 s, respectively. They form through various processes, including collisional excitation by thermal and non-thermal electrons, photodissociation of H₂O or O₂, and dissociative recombination of ${\text{O}}_2^{+}$. As both O(¹D)) and O(¹S) may act as efficient coolants by emitting in the optical and ultraviolet, it is relevant to ask whether collisions with photoelectrons might shorten their effective lifetimes, thereby potentially modifying the energy budget of the atmosphere.

To answer that question, we have gathered the basic set of cross sections shown in Fig. A.1 for collisions of O(¹D). They include deexcitation into O(³P), excitation into O(¹S) and ionization into O⁺(⁴S^o). Following the small abundance approximation described in García Muñoz (2023), we determined the lifetimes for O(¹D) in deexcitation, excitation, and ionization collisions with photoelectrons. We know that O(¹D) photoionizes at wavelengths ≲83OÅ, a spectral region effectively shielded by the dominating neutral gases. Our calculations confirm that O(¹D) photoionization is slow compared to spontaneous emission or quenching by thermal electrons or heavy particles and it is thus omitted from the discussion.

The results of the calculations are shown in Fig. 6 for π Men c and TRAPPIST-1 b. For comparison, we also show the O(¹D) radiative lifetime, the quenching lifetimes in collisions with thermal electrons (rate coefficients based on Barklem 2007), O atoms (Jamieson et al. 1992), and H atoms (Krems et al. 2006). For TRAPPIST-1 b, we also give the lifetime for chemical reaction with H₂O molecules (rate coefficient from Burkholder et al. 2019). Clearly, the lifetime of O(¹D) is negligibly affected by photoelectron collisions. Rather, deexcitation collisions with O and H atoms and chemical reactions with H₂O control the O(¹D) lifetime in the mostly neutral atmosphere; while deexcitation collisions with thermal electrons take the control when the fractional ionization reaches a few times 10⁻³. The situation is reminiscent of the modeling of stellar atmospheres, where it is thought that the population of oxygen levels is affected by inelastic collisions with both hydrogen atoms and thermal electrons, with the importance of each type of collision depending on the stellar effective temperature (Allende Prieto et al. 2003; Barklem 2018).

The cross sections for the deexcitation, excitation, and ionization of O(¹S) in collisions with photoelectrons are comparable or smaller than those for O(¹D) (Barklem 2007; Tayal & Zatsarinny 2016). Given that the radiative lifetime of O(¹S) is two orders of magnitude shorter than that of O(¹D), it turns out that the lifetime of O(¹S) is also negligibly affected by photo-electron collisions. The corresponding calculations of lifetimes are not presented.

Fig. 6 Some lifetimes for O(1D). I: Radiative emission, O(1 D)→O(3P)+hv; II, III, and IV: Quenching e′+O(1 D)→e"+O(3P), excitation e'+O(1D)→e"+O(1S), and ionization e′+O(1D)→e"+e"'+O+ by photoelectrons, respectively; V, VI and VII: quenching by thermal electrons eth+O(1D)→eth+O(3P), quenching by oxygen atoms O(3P)+O(1D)→O(3P)+O(3P), quenching by hydrogen atoms H(1)+O(1D)→H(1)+O(3P), respectively; VIII (for TRAPPIST-1 b only): chemical transformation O(1D)+H2O→2OH.

4 Summary and perspectives

We simulated the slowing down of photoelectrons in O-H mixtures, determining the yields for each elastic and inelastic collisional channel and the heating efficiencies. Our work extends the calculations done for hydrogen-dominated atmospheres by Cecchi-Pestellini et al. (2006, 2009) and Shematovich et al. (2014) to atmospheric compositions that may be more relevant in the characterization of small exoplanets.

We find that the fraction E_elas/E₀ of the photoelectron energy that goes into elastic collisions (and thus into heat) for an O–H mixture is not very different from what is predicted for pure H. The exception to this behaviour occurs at fractional ionizations x_e≪l0⁻⁴, as in these conditions E_elas/E₀ can be significantly lower than for pure H. The difference arises from the fact that excitation in the O atom remains possible down to energies of about 2 eV, unlike the case of the H atom. We also find that for a given rate of primary photoelectrons, the total (=primary + secondary) number of electrons produced in an O–H mixture is higher than for pure H.

We have investigated the significance of non-thermal effects in the upper atmospheres of two possibly water-rich exoplanets. We predict that the atmospheres that contain abundant O will ionize faster than those that are composed of pure H. The reason for this is that the O atom has cross sections for photoionization and photoelectron-collision ionization that are notably larger than those for H. Faster ionization of O, along with subsequent ionization of H through charge exchange in the atmosphere, will surely have an impact on the detectability of H atoms escaping the atmospheres of water-rich exoplanets through Lyman-α transmission spectroscopy. This finding reinforces the idea postulated in García Muñoz et al. (2020) that the published detections of Lyman-α absorption are correlated with bulk densities of the planets below 2 and 3 g cm⁻³, which is likely indicative of low-metallicity atmospheres. At the same time, the non-detections are consistent with higher bulk densities consistent with high-metallicity atmospheres. Subsequent detections of Lyman-α absorption for HAT-P-11 b (p_p~1.6 g cm⁻³; Ben-Jaffel et al. 2022) and HD 63433 c (~2.l g cm⁻³; Zhang et al. 2022), and non-detections for HD 63433 b (~3 g cm⁻³; Zhang et al. 2022), K2-25 b (~3.3 g cm⁻³; Rockcliffe et al. 2021), GJ 9827 b (~7.5 g cm⁻³; Carleo et al. 2021) are consistent with the idea, even though the detection of Lyman-α absorption may depend on other factors as well such as the overall stellar XUV flux or the stellar wind strength (e.g., Shaikhislamov et al. 2020).

We argue that omitting the effect of photoelectrons can artificially boost the heating of water-rich atmospheres by a factor of a few in the layers where most of the high-energy stellar photons are deposited and the atmospheric outflow is established. The quenching of metastable levels of O and O⁺ produced by the photoelectrons does not appreciably change these findings. A weaker heating rate will surely extend the long-term stability of water-rich atmospheres.

As we strive to understand the nature of water-rich exoplan-ets, the models that we are developing ought to explore the main atmospheric processes that play a role in them. That mission has not been accomplished yet and in what follows, we list some of the caveats of the present work. Our treatment of the collisional-radiative processes remains incomplete as it does not consider the excitation of O and O⁺ in collisions with thermal electrons and other heavy particles, nor diffuse continuum radiation that might be important at the altitudes where the atmosphere transitions from optically thick to optically thin at XUV wavelengths. Additionally, a change in the heating efficiency of the atmosphere attributable to non-thermal processes will modify the thermal cooling of the gas and, inversely, changes in the temperature due thermal cooling will affect the overall structure of the atmosphere and its fractional ionization - thereby modifying the effect of non-thermal processes. Our exploration of the TRAPPIST-1 b system reveals that the conversion of H₂O into its atomic components depends strongly on the specifics of the stellar spectrum. We have not explored how efficiently the photoelectrons produced by deposition of XUV photons will dissociate and ionize the H₂O molecule in the transition region between the lower and upper atmospheres, as this calls for additional work. Lastly, we have not explored the feedbacks between thermal and non-thermal processes and the hydrodynamics of the flow. The combined effects of these processes may contribute constructively or destructively in ways that need to be quantified, and we will explore these issues in due course.

Appendix A Collisional and radiative properties of the model

The MC model for the slowing down of photoelectrons in a pure H gas is described in García Muñoz (2023). It considers: one elastic channel for electron-electron collisions; one elastic channel for electron-H collisions; six inelastic channels for electron collisions with the H atom in its ground level that lead to excitation into levels with principal quantum numbers up to five and into a pseudo-level representative of higher-energy levels; one inelastic channel for ionization. The model can optionally treat the excitation, deexcitation and ionization of various excited levels of the H atom.

In what follows, we describe the expansion of the MC model to take into account collisions of photoelectrons with O and O⁺. Typically, a calculation involves one elastic channel for electron-O(³P) collisions, 66 inelastic channels for collisions with O(³P) and four inelastic channels for collisions with O⁺(⁴ S^o).

A.1 Collisions with the ground-level atoms O(³P) and O⁺(⁴S^o). Sensitivity to different sets of cross sections

Momentum-transfer cross sections for elastic collisions of the ground-level atom O(³P) have been reported in the literature, for instance: Itikawa (1978) and Tayal & Zatsarinny (2016). Both sets are consistent and we adopted the most recent one, which covers a broader range of energies. Laher & Gilmore (1990) compiled a set (hereafter the LG90 set) of inelastic cross sections for O(³P) that remains comprehensive to date. After correction for auto-ionization, it includes 61 channels for excitation within the neutral atom, one channel for ionization into the ground level ion O⁺(⁴S^o), 3 channels for ionization-excitation into O⁺(²D^o, ²P^o, ⁴P), and one channel for double ionization into O²⁺(³P). The latter is included for completeness, but given its small cross sections, the effect of double ionization is minor in our calculations.

There have been more recent calculations of some of the inelastic cross sections in the LG90 set (e.g., Barklem 2007; Tayal & Zatsarinny 2016), which motivated us to assess the sensitivity to them of our results. Our LG90-TZ16 set of cross sections relies on the LG90 set but, when possible, we replaced the older determinations by the Tayal & Zatsarinny (2016) calculations. In particular, we adopted, from the latter work, the total cross sections for single ionization, which include auto-ionization. For the relative branching ratios into ⁴S^o, ²D^o,²P^o, and ⁴P, we implemented the branching ratios from the LG90 set. For the LG90-TZ16 set, we also updated the cross sections for the channels with end levels (excitation energy in parenthesis, in eV): 2p^{4 1}D (1.96), 2p^{4 1}S (4.18), 2p³(⁴S^o)3s ⁵S^o (9.14), 2p³(⁴S^o)3s ³S^o (9.51), 2p³(⁴S^o)3p ⁵P (10.73), 2p³(⁴S^o)3p ³P (10.98), 2p³(⁴S^o)3d ³D^o (12.08), 2p³(²D^o)3s ³D^o (12.53), 2p³(²P^o)3s ³P^o (14.11), and 2s2p^{5 3}P^o (15.65). The update includes excitation into O(¹ D), which has the largest cross section of all excitation channels. In particular, we assumed for the LG90-TZ16 set that excitation into 2p³(²P^o)3s ³P^o occurs without ionization; we also assumed that excitation into 2s2p^{5 3}P^o occurs with the excited level always auto-ionizing and therefore its excitation cross section is zero. We extracted the data from Tayal & Zatsarinny (2016) by scanning the corresponding graphs. As needed, the Laher & Gilmore (1990) data were used at higher energies.

We additionally considered inelastic collisions of the pho-toelectrons with the ground level ion O⁺ (⁴S^o), including ion-ization into O²⁺(³P) (cross sections from Bell et al. 1983) and excitation into O+(²D^o, ²P^o, ⁴P) (cross sections from Itikawa et al. 1985). The excitation cross sections were extrapolated beyond their stated range of validity, although this should have a negligible effect as ionization seems to dominate at high energies. The cross sections for collisions with O⁺ (⁴S^o) are presented in Fig. A.1 (top).

For the probability distribution function of the energy of post-collision electrons, we used the formulation described in García Muñoz (2023). In particular, we implemented Ē (O(³ P)=13.3 eV as an average of the equivalent factors reported in Itikawa & Ichimura (1990) for the single ion-ization channels and, for simplicity, we adopted the same value for the double ionization channel. In double ioniza-tion collisions, the MC model determines the energy of one of the post-collision photoelectrons directly from the probability distribution function P(E''; E'); it is then assumed that the other two post-collision photoelectrons take each of them half of the remaining energy. When needed, we took (Kozma & Fransson 1992), Ē (O(¹ D))=0.6×IP(O(¹ D))=7 eV and Ē (O⁺(⁴S^o))=0.6×IP(O+(⁴S^o))=21 eV.

We carried out the calculations presented in Figs. 2 and 3 with both the LG90 and LG90-TZ16 sets of cross sections. We found that the calculated E_elas/E₀ and Φ_i differed by typically less than 10%. This amount that serves as a reasonable guess for the accuracy of our calculations, which should be dominated by the accuracy of the cross sections - and not by the solution method. The discrepancies are largely dictated by the dominant inelastic channels, namely: excitation into O(¹ D) below ~20 eV and ionization at higher energies. The rest of calculations presented in this work were done with the LG90 set.

A.2 Collisions with O(¹ D)

We considered the superelastic channel for deexcitation into O(³P), the cross sections which were obtained from the application of detailed balancing to the excitation cross sections, along with two inelastic channels for excitation into O(¹S) and single ionization into O+(⁴S^o). We took the inelastic cross sections from Tayal & Zatsarinny (2016) and confirmed them with Fig. 1 of Barklem (2007), finding that it is probably safe to neglect other excitation channels. The momentum-transfer cross sections for O(¹ D) are comparable in magnitude to those for O(³P), and given the much lower abundances expected for O(¹ D), the corresponding collisional channel needs not be included. The cross sections for collisions with O(¹ D) are presented in Fig. A.1 (bottom).

Fig. A.1 Adopted cross sections for the excitation and ionization of O+(4So) (top) and for (de)excitation and ionization of O(1 D) (bottom). Note: in the top panel, the black curve labeled O+,*· includes all the excitation channels that are shown in the graph.

References

All Figures

Fig. 1 Cross sections for excitation (H*, O*), ionization (single, H+, O+, and double, O2+) and total inelastic (=excitation + ionization) collisions of photoelectrons with ground-level H and O atoms (LG90 set, see Appendix A). The H cross sections have been multiplied by ×2 to better reveal the significance of collisions with both types of atoms in an O–H mixture with ƒOn/Hn = 1/2. We also show the adopted electron-electron pseudo-cross sections multiplied by ×10−3.

In the text

Fig. 2 Some calculated properties from the slowing down of photoelectrons. Top: fraction of the photoelectron initial energy E0 that is transferred to the gas in elastic collisions. Bottom: number of secondary ions created by a photoelectron of initial energy E0, normalized by E0. The three color-style combinations refer to pure H (black, long dashed), pure O (black, dotted), and an O–H mixture with ƒOn/Hn = 1/2 (blue, solid). For each mixture, the four curves that are displayed (distinguishable by their thickness) refer to fractional ionizations log10(xe) = −1 (thickest), −2, −3, and −4 (thinnest). Although not shown for the sake of clarity, the Eelas/E0 ratio for E0≳100eV for pure H becomes essentially independent of fractional ionization at xe<10−4. In contrast, Eelas/E0 continues to drop with xe for the other O–H mixtures. For example, for E0>200eV and xe~10−6, Eelas/E0~0.05, and 0.07 for pure O and the O–H mixture with ƒOn/Hn=1/2, respectively. Calculations were done with the LG90 set of cross sections.

In the text

Fig. 3 For an O–H mixture, ƒOn/Hn=1/2 is the fraction of the photoelectron energy, E0, that goes into various inelastic channels specified by their initial and end levels. Then, H* and O* refer to the summation over all excited levels within the corresponding neutral atoms. For each panel, the four curves that are displayed (distinguishable by their thickness) refer to fractional ionizations log10(xe)=−1 (thickest), −2, −3, and −4 (thinnest). Calculations done with the LG90 set of cross sections.

In the text

Fig. 4 Some properties of the atmospheres of π Men c and TRAPPIST-1 b. Top: adopted number densities of H, H+, O, O+ and eth,. See text for details. For TRAPPIST-1 b, the densities of additional molecules (H2O, O2, H3O+, ${\text{O}}_2^{+}$) that remain relatively abundant up to 10−3 dyn cm−2 are also shown. Middle: production rates of primary and secondary ions. Bottom: heating rates in the no- and full-quenching limits (see text), and under the assumption that photoelectrons transfer all their kinetic energy into heat.

In the text

	Fig. 5 Energy spectra of the nascent photoelectrons produced through photoionization of O and H atoms at different levels in the atmosphere of π Men c for ƒOn/Hn=l/2. The quoted numbers in the plot represent the pressure levels.
In the text

Fig. 6 Some lifetimes for O(1D). I: Radiative emission, O(1 D)→O(3P)+hv; II, III, and IV: Quenching e′+O(1 D)→e"+O(3P), excitation e'+O(1D)→e"+O(1S), and ionization e′+O(1D)→e"+e"'+O+ by photoelectrons, respectively; V, VI and VII: quenching by thermal electrons eth+O(1D)→eth+O(3P), quenching by oxygen atoms O(3P)+O(1D)→O(3P)+O(3P), quenching by hydrogen atoms H(1)+O(1D)→H(1)+O(3P), respectively; VIII (for TRAPPIST-1 b only): chemical transformation O(1D)+H2O→2OH.

In the text

	Fig. A.1 Adopted cross sections for the excitation and ionization of O+(4So) (top) and for (de)excitation and ionization of O(1 D) (bottom). Note: in the top panel, the black curve labeled O+,*· includes all the excitation channels that are shown in the graph.
In the text