Heating and ionization by non-thermal electrons in the upper atmospheres of water-rich exoplanets

Context. The long-term evolution of an atmosphere and the remote detectability of its chemical constituents are susceptible to how the atmospheric gas responds to stellar irradiation. The response remains poorly characterized for water and its dissociation products, however, this knowledge is relevant to our understanding of hypothetical water-rich exoplanets. Aims: Our work investigates the effect of photoelectrons, namely, the non-thermal electrons produced by photoionizing stellar radiation on the heating and ionization of extended atmospheres dominated by the dissociation products of water. Methods: We used a Monte Carlo model and up-to-date collision cross sections to simulate the slowing down of photoelectrons in O-H mixtures for a range of fractional ionizations and photoelectron energies. Results: We find that that the fraction of energy of a photoelectron that goes into heating is similar in a pure H gas and in O-H mixtures, except for very low fractional ionizations, whereby the O atom remains an efficient sink of energy. The O-H mixtures will go on to produce more electrons because the O atom is particularly susceptible to ionization. We quantified all that information and present it in a way that can be easily incorporated into photochemical-hydrodynamical models. Conclusions: Neglecting the role of photoelectrons in models of water-rich atmospheres will result in overestimations of the atmospheric heating and, foreseeably, the mass-loss rates as well. It will also underestimate the rate at which the atmospheric gas becomes ionized, which may have implications for the detection of extended atmospheres with Lyman-{\alpha} transmission spectroscopy. Our simulations for the small exoplanets {\pi} Men c and TRAPPIST-1 b reveal that they respond very differently to irradiation from their host stars, with water remaining in molecular form at lower pressures in the latter case.


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. This process would occur 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 2 /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 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 M p and R p may be orders of magnitude larger if the gas has a high molecular weight than if it is composed of H 2 /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. ( , 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 2 /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 3 if the atmospheric composition is dominated by H 2 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 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;Yoshida et al., 2022) and the relevance of related ideas to exoplanets.
We are currently revisiting the formulation of the collisionalradiative 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 spectroscopic 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 2 -H mixtures and for pure O and O 2 (e.g. Green et al., 1977;Victor et al., 1994;Dalgarno et al., 1999;García Muñoz, 2023, and refs. therein). Some calculations have also been reported for exoplanets with hydrogen-dominated atmospheres (Cecchi-Pestellini et al., 2006, 2009Shematovich et al., 2014;Guo & Ben-Jaffel, 2016). We have however not found similar calculations for H 2 O nor O-H mixtures. There is a long history of simulations of the slowing down of photoelectrons in the H 2 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 2 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. 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( 1 D) in them. Lastly, Section 4 summarizes the main findings and lists some possible next steps for future studies.

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: if the exchange with the thermal electrons, e th , or the heavy particles, X, involves only kinetic energies -or inelastic: e ′ + X − → e ′′ + X * , e ′′ + e ′′′ + 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 photoionization 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 2+ 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 We also show the adopted electron-electron pseudo-cross sections multiplied by ×10 −3 .
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 O-H 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 collisional 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 . We adopt [X n ]=10 10 cm −3 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 f O n /H n =0, 1/6, 1/2, 3/2, ∞, and for log 10 (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 0 , 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 available inelastic channels. From energy conservation, E 0 =E elas + j E inel, j . 1 It is usually accepted that all of E elas contributes to heating the gas and that some of j E inel, j may also contribute locally. Using E h to denote the part of E 0 that goes into heating, it follows that E 0 ≥E h ≥E elas . The ratio E h /E 0 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 f O n /H n =1/2, as well as for pure H and pure O. The top panel shows E elas /E 0 and the bottom one shows the normalized yield Φ i /E 0 . The ionization yield, Φ i , represents the number Fig. 2. Some calculated properties from the slowing down of photoelectrons. Top: Fraction of the photoelectron initial energy E 0 that is transferred to the gas in elastic collisions. Bottom: Number of secondary ions created by a photoelectron of initial energy E 0 , normalized by E 0 . The three color-style combinations refer to pure H (black, long dashed), pure O (black, dotted), and an O-H mixture with f O n /H n =1/2 (blue, solid). For each mixture, the four curves that are displayed (distinguishable by their thickness) refer to fractional ionizations log 10 (x e )=−1 (thickest), −2, −3, and −4 (thinnest). Although not shown for the sake of clarity, the E elas /E 0 ratio for E 0 ≳100 eV for pure H becomes essentially independent of fractional ionization at x e <10 −4 . In contrast, E elas /E 0 continues to drop with x e for the other O-H mixtures. For example, for E 0 ≳200 eV and x e ∼10 −6 , E elas /E 0 ∼0.05, and 0.07 for pure O and the O-H mixture with f O n /H n =1/2, respectively. Calculations were done with the LG90 set of cross sections. of ions produced per injected photoelectron of energy, E 0 , regardless whether they originated from H, O, or O + and of their charge. Because the MC model considers double ionization for O( 3 P), which leads to O 2+ ( 3 P) plus two electrons, Φ i is slightly larger than the number of secondary electrons produced.
It is seen in Fig. 2 (top) that E elas /E 0 depends on the gas composition and also on the photoelectron initial energy, E 0 , and the fractional ionization of the mixture, x e . For E 0 <10.2 eV, E elas /E 0 ≡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( 1 D) excitation) that remain competitive against electronelectron collisions for small-to-moderate x e . For E 0 ≳100 eV and x e ∼10 −3 -10 −2 , E elas /E 0 is typically higher for an O-H mixture with f O n /H n =1/2 than for pure H. The trend does, however, re- Fig. 3. For an O-H mixture, f O n /H n =1/2 is the fraction of the photoelectron energy, E 0 , 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 log 10 (x e )=−1 (thickest), −2, −3, and −4 (thinnest). Calculations done with the LG90 set of cross sections. verse for both lower and higher fractional ionizations. For pure H and E 0 ≳50 eV, the heating efficiency E elas /E 0 becomes nearly independent of the fractional ionization for x e ≲10 −4 . In contrast, E elas /E 0 drops to ∼0.07 (not shown in Fig. 2) in the O-H mixture with f O n /H n =1/2 and x e =10 −6 because for such small fractional ionizations excitation into O( 1 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 f O n /H n =1/2, and shows the ratio E inel,j /E 0 for the various excitation and ionization channels -except ionization from either O( 3 P) or O + ( 4 S o ) to O 2+ , which remain minor. We note that excitation is typically more significant than ionization for H, whereas the reverse oc-

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 0 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 ionization 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).  2021)) is part of a system of seven known planets, some of which orbit in the habitable zone of their ultracool dwarf host 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 −2 ) that the atmosphere is pure H 2 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 −3 dyn cm −2 , the chemistry of TRAPPIST-1 b becomes particularly interesting because the weak photodissociation at far-ultraviolet wavelengths facilitates the survival of molecules. Indeed, H 2 O and O 2 are abundant in that region amongst the neutrals, with H 3 O + and O + 2 amongst the ions. At lower pressures, the H 2 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 −3 dyn cm −2 , 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 0 ≡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 photoelectrons 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 0 (≈hc/λ−13.6 eV for both H and O, where for consistency with the photochemicalhydrodynamical model it is assumed that O photoionization results in ground-level O + ) produces a number of secondary electrons that depends on E 0 and the local composition and fractional ionization of the atmosphere. We used the number densities of H, H + , O, O + , and e th 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.

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 −2 dyn cm −2 . 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 −3 dyn cm −2 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 2 O) 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 ≲1 dyn cm −2 , 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). 2 This occurs even if the lines in the Lyman series, which connect with the ground level H(1s), are optically thick. Also, the excited level H(2s), which does not connect with H(1s) 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( 1 D), O( 1 S ), O + ( 2 D), and O + ( 2 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, Alternatively, in the full-quenching limit, which should apply when many deexcitation collisions occur before the atom radiates, 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:  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 0 ≡1, 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. . Following the small abundance approximation described in García Muñoz (2023), we determined the lifetimes for O( 1 D) in deexcitation, excitation, and ionization collisions with photoelectrons. We know that O( 1 D) photoionizes at wavelengths ≲830 Å, a spectral region effectively shielded by the dominating neutral gases. Our calculations confirm that O( 1 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( 1 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 2 O molecules (rate coefficient from Burkholder et al. (2019)). Clearly, the lifetime of O( 1 D) is negligibly affected by photoelectron collisions. Rather, deexcitation collisions with O and H atoms and chemical reactions with H 2 O control the O( 1 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 −3 . 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( 1 S ) in collisions with photoelectrons are comparable or smaller than those for O( 1 D) (Barklem, 2007;Tayal & Zatsarinny, 2016). Given that the radiative lifetime of O( 1 S ) is two orders of magnitude shorter than that of O( 1 D), it turns out that the lifetime of O( 1 S ) is also negligibly affected by photoelectron collisions. The corresponding calculations of lifetimes are not presented.

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 0 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 <<10 −4 , as in these conditions E elas /E 0 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  that the published detections of Lyman-α absorption are correlated with bulk densities of the planets below 2-3 g cm −3 , 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 ∼1.6 g cm −3 ; Ben-Jaffel et al. (2022) 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 exoplanets, 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 2 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 2 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.     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( 3 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 + ( 4 S o ), 3 channels for ionization-excitation into O + ( 2 D o , 2 P o , 4 P), and one channel for double ionization into O 2+ ( 3 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 4 S o , 2 D o , 2 P o , and 4 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 3 ( 4 S o )3s 5 S o (9.14), 2p 3 ( 4 S o )3s 3 S o (9.51), 2p 3 ( 4 S o )3p 5 P (10.73), 2p 3 ( 4 S o )3p 3 P (10.98), 2p 3 ( 4 S o )3d 3 D o (12.08), 2p 3 ( 2 D o )3s 3 D o (12.53), 2p 3 ( 2 P o )3s 3 P o (14.11), and 2s2p 5 3 P o (15.65). The update includes excitation into O( 1 D), which has the largest cross section of all excitation channels. In particular, we assumed for the LG90-TZ16 set that excitation into 2p 3 ( 2 P o )3s 3 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 photoelectrons with the ground level ion O + ( 4 S o ), including ionization into O 2+ ( 3 P) (cross sections from Bell et al., 1983) and excitation into O + ( 2 D o , 2 P o , 4 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 + ( 4 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 E(O( 3 P))=13.3 eV as an average of the equivalent factors reported in Itikawa & Ichimura (1990) for the single ionization channels and, for simplicity, we adopted the same value for the double ionization channel. In double ionization 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( 1 D))=0.6×IP(O( 1 D))=7 eV andĒ(O + ( 4 S o ))=0.6×IP(O + ( 4 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 0 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( 1 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( 1 D)
We considered the superelastic channel for deexcitation into O( 3 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( 1 S ) and single ionization into O + ( 4 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( 1 D) are comparable in magnitude to those for O( 3 P), and given the much lower abundances expected for O( 1 D), the corresponding collisional channel needs not be included. The cross sections for collisions with O( 1 D) are presented in