<figure><img src="image_1.png" /></figure>

# A Multispecies Atmospheric Escape Model with Excited Hydrogen and Helium: Application to HD209458b

Anna Ruth Taylor<fnref n="1" />, Tommi T. Koskinen<fnref n="1" />, Luca Argenti<fnref n="2" />, Nicholas Lewis<fnref n="2" />, Chenliang Huang (黃辰亮)<fnref n="3" />, Anthony Arfaux<fnref n="4" />, and Panayotis Lavvas<fnref n="4" />

## Abstract

Atmospheric escape shapes exoplanet evolution and star-planet interactions, with He I 10830 Å absorption serving as a key tracer of mass loss in hot gas giants. However, transit depths vary significantly across observed systems for reasons that remain poorly understood. HD209458b, the archetypal hot-Jupiter, exhibits relatively weak He I 10830 Å and Hα absorption, which has been interpreted as evidence for a high H/He ratio (98/2), possibly due to diffusive separation. To investigate this possibility and other processes that control these transit depths, we reassess excitation and de-excitation rates for metastable helium and explore the impact of diffusion processes, stellar activity, and tidal forces on the upper atmosphere and transit depths using a model framework spanning the whole atmosphere. Our model reproduces the observed He I transit depth and Hα upper limit, showing strong diffusive separation. We match the observations assuming a photoelectron efficiency of 20%-40%, depending on the composition of the atmosphere, corresponding to mass-loss rates of $1.9-3 \times 10^{10} \text{ g s}^{-1}$. We find that the He I 10830 Å transit depth is sensitive to both stellar activity and diffusion processes, while Hα is largely unaffected due to its strong dependence on Lyα excitation. These differences may help explain the system-to-system scatter seen in population-level studies of the He I line. While He I data alone may not tightly constrain mass-loss rates or temperatures, they do confirm atmospheric escape and help narrow the viable parameter space when interpreted with physically motivated models. Simultaneous observations of He I, Hα, and stellar activity indicators provide powerful constraints on upper atmosphere dynamics and composition, even in the absence of full transmission spectra.

*Unified Astronomy Thesaurus concepts: Upper atmosphere (1748); Exoplanet atmospheric evolution (2308); Exoplanet atmospheric structure (2310)*

*Materials only available in the online version of record: data behind figure*

## 1. Introduction

Exoplanets orbiting close to their host stars receive high-energy stellar radiation that causes their atmospheres to escape. Understanding atmospheric escape allows us to gain insights into the chemical processes occurring within exoplanetary atmospheres, and escape significantly influences the distribution and characteristics of exoplanets in terms of size, composition, and atmospheric retention. The population of known exoplanets detected by the Kepler Space Telescope and the Transiting Exoplanet Survey Satellite includes features such as the gap in planet occurrence near 1.8 $R_⊕$, known as the radius valley (B. J. Fulton et al. 2017), and the absence of Neptunes and mini-Neptunes at close-in orbits, known as the hot Neptune desert (T. Mazeh et al. 2016). Models of atmospheric escape and radius evolution show that these features can be produced by photoevaporation driven by stellar X-ray and extreme ultraviolet (XUV and EUV, respectively) radiation and envelope loss due to rapid Roche lobe overflow (e.g., J. E. Owen & Y. Wu 2017; T. T. Koskinen et al. 2022). To understand the evolution of the exoplanet population, we

must understand the escape of the upper atmospheres of exoplanets because atmospheric loss shapes their characteristics over time.

Among the exoplanets discovered to date, hot Jupiters and warm Neptunes are especially useful for investigating atmospheric escape because of their relatively large size and close proximity to the host stars. On these exoplanets, models of escaping exoplanet upper atmospheres can be validated by comparisons with observations. In the past, observations of exoplanet upper atmospheres were achieved from space with the Hubble Space Telescope (HST) and were limited to far-ultraviolet (FUV) wavelengths where the signal-to-noise ratio is relatively low, and the interpretation of the observations has proven challenging (e.g., A. Vidal-Madjar et al. 2003, 2004; L. Ben-Jaffel & S. Sona Hosseini 2010; J. L. Linsky et al. 2010; G. E. Ballester & L. Ben-Jaffel 2015). Recently, however, new observations have been accumulating, including ground- and space-based detections of excited H and He, as well as various heavy elements and metals across UV to near-IR wavelengths (e.g., A. G. Jensen et al. 2012; J. J. Spake et al. 2018; N. Casasayas-Barris et al. 2020). Simultaneous observations of H and He are especially valuable because they probe the dominant neutral species of giant planet atmospheres and constrain the thermal structure and non-local-thermodynamic-equilibrium (non-LTE) level populations in the upper

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="2">² The University of Central Florida, Department of Physics, Orlando, FL 32816, USA</footnote>

<footnote n="3">³ Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, People's Republic of China</footnote>

<footnote n="4">⁴ Groupe de Spectrométrie Moléculaire et Atmosphérique, Université de Reims Champagne Ardenne, Reims, France
Received 2025 April 14; revised 2025 May 30; accepted 2025 June 9; published 2025 August 5</footnote>
atmosphere. These observations also probe altitudes from the
middle atmosphere to the exosphere, providing powerful
constraints on the energy balance of the upper atmosphere
and mass-loss rates.

Absorption by He I at 10830 Å and the H I Balmer lines (6563, 4861, 4340 Å) have been detected on many hot-Jupiter-type exoplanets (e.g., J. J. Spake et al. 2018; V. Krishnamurthy & N. B. Cowan 2024; A. Masson et al. 2024). The He I 10830 Å line corresponds to the transition between the *n* = 2 and *n* = 3 electronic levels of the helium atom. The lower level of this transition is the metastable 2³S state, which has a half-life of ~7800 s (G. W. F. Drake 1971), a lifetime that allows He I (2³S) to persist long enough for substantial density to build up and absorption to occur. The H I Balmer lines result from transitions of electrons in the hydrogen atom from the *n* = 2 energy level to higher energy levels. None of these lines are strongly absorbed by the interstellar medium (ISM), and they can both be observed from the ground (N. Indriolo et al. 2009). This means that they are more easily observable than FUV absorption lines such as the H Lyα line that is strongly absorbed by the ISM and can only be observed from space (e.g., A. Vidal-Madjar et al. 2003; J. L. Linsky et al. 2010; T. T. Koskinen et al. 2013a).

Studies that have investigated trends in the He I 10830 Å transit depths across many systems have concluded that they correlate with stellar XUV flux (V. Krishnamurthy & N. B. Cowan 2024; J. Sanz-Forcada et al. 2025) or mass-loss rates (G. Ballabio & J. E. Owen 2025), but there is still a significant scatter of the observations deviating from these proposed trends (D. Linssen et al. 2025; A. P. Allan & A. A. Vidotto 2025). This suggests that additional factors—such as atmospheric composition, heating structure, and planetary dynamics—affect He I 10830 Å absorption in ways that are not captured by simple scaling relations. In some cases, the observed absorption probes relatively low altitudes in the thermosphere (e.g., F. J. Alonso-Floriano et al. 2019), while in others, excited helium is observed far from the planet in leading and trailing tails (e.g., J. J. Spake et al. 2018; L. Wang & F. Dai 2021a, 2021b; M. Zhang et al. 2022), and there is no straightforward explanation for the differences. Even for the archetypal hot Jupiter HD209458b, which orbits a relatively quiet Sun-like G-type star at 0.047 au (J. Southworth 2010), models have had to impose a relatively low H/He ratio in the upper atmosphere (around 98/2) to match the relatively low He I 10830 Å transit depth. This has led to the suggestion that the helium abundance in many exoplanetary atmospheres is substantially lower than solar (F. Biassoni et al. 2024), or that the helium is experiencing diffusive separation (M. Lampón et al. 2020).

Given these challenges, we revisit HD209458b to system-
atically explore the key factors controlling He I 10830 Å transit
depths. HD209458b serves as an ideal benchmark for studying
atmospheric escape due to its well-characterized system
parameters, extensive observational data, and prior transit
depth studies probing the upper atmosphere. Past studies have,
however, predicted a wide range of mass-loss rates
1–10 × 10¹⁰ g s⁻¹(A. Oklopčić & C. M. Hirata 2018;
M. Lampón et al. 2020; A. Masson et al. 2024) based on
He I 10830 Å absorption, and previous models of this absorp-
tion have often relied on ad hoc assumptions about helium
abundance to match the observations. Additionally,
HD209458b’s debated FUV transit signals highlight the need

to better understand how He I 10830 Å absorption relates to
other probes of atmospheric escape. By reanalyzing this
system with an updated treatment of some of the key physics,
we aim to clarify why the He I 10830 Å absorption by the
planet’s upper atmosphere is lower than might be expected
based on other observations and what this implies for broader
trends in observations of exoplanet atmospheres.

The models that are often used to fit He I 10830 Å transit
observations make several assumptions about thermal struc-
ture, escape mechanisms, and atmospheric composition that
affect the predicted transit depths. For example, He I (2³S)
observations are often interpreted using 1D Parker wind
models, such as the publicly available p-winds model
(L. A. Dos Santos et al. 2022) to describe the structure of the
escaping atmosphere. These models often assume an iso-
thermal outflow and do not consider much of the detailed
physics of atmospheric escape. They result in degenerate
solutions for temperature, mass-loss rate, density, and velocity,
meaning that multiple combinations of these parameters can fit
a given transit depth. Additionally, they do not properly
account for the coupling between the lower/middle atmos-
phere and the thermosphere, which can play an important role
in setting the boundary conditions for atmospheric escape
(J. M. Chadney et al. 2017; C. Huang et al. 2023). Many
current models also do not explicitly include diffusion, leading
to ad hoc assumptions about mixing, particularly in the case of
helium-to-hydrogen separation. Finally, many models do not
have the ability to simultaneously model hydrogen and helium
absorption, despite the fact that analyzing these two features
together potentially enables cross-validation of derived
thermal structures and non-LTE level populations when
observations of both lines are available.

While diffusive separation of helium has been proposed as a
possible explanation for the relatively weak He I
10830 Å absorption observed on HD209458b, only a handful
of models have explicitly included this process (e.g., L. Xing
et al. 2023; M. Schulik & J. Owen 2025). These models solve
separate momentum equations while we show in this work that
a simpler approach can reproduce similar results. Additionally,
previous models did not account for the effects of eddy
diffusion, which is expected to impact helium mixing and
transit depths. Finally, He I 10830 Å transit depths depend on
the choice of excitation and de-excitation rates, including
photoionization cross sections. Past studies have relied on
lower-resolution photoionization cross sections for the 2³S
state of He that neglect strong resonances in the EUV that
overlap with prominent stellar emission lines (L. S. Waldrop
et al. 2005; N. N. Shefov et al. 2009; A. Oklopčić &
C. M. Hirata 2018; M. Lampón et al. 2020; L. A. Dos Santos
et al. 2022). We provide updated cross sections based on
L. Argenti & R. Moccia (2006, 2008) and L. Argenti &
E. Lindroth (2023) here. We also update the Penning
ionization rate to a temperature-dependent value using
available cross sections (J. S. Cohen & N. F. Lane 1971;
H. Morgner & A. Niehaus 1979). Finally, motivated by studies
of He I (2³S) in Earth’s atmosphere (L. S. Waldrop et al. 2005;
N. N. Shefov et al. 2009), we explore the impact of
photoelectron collisions on helium excitation, but find that
their effect is negligible for HD209458b.

To address the open questions about how differences in
model atmospheres and reaction rates impact He I 10830 Å and
Hα transit depths, we use a self-consistent, 1D model of the
<table><caption>**Table 1**
HD209458b System Parameters (J. Southworth 2010)</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Stellar Mass</td><td>M<sub>*</sub></td><td>1.148 M<sub>☉</sub></td></tr><tr><td>Stellar Radius</td><td>R<sub>*</sub></td><td>1.162 R<sub>☉</sub></td></tr><tr><td>Planetary Mass</td><td>M<sub>p</sub></td><td>0.714 M<sub>J</sub></td></tr><tr><td>Planetary Radius</td><td>R<sub>p</sub></td><td>1.380 R<sub>J</sub></td></tr><tr><td>Orbital Period</td><td>P<sub>orb</sub></td><td>3.524749 days</td></tr><tr><td>Semimajor axis</td><td>a</td><td>0.04747 au</td></tr></tbody></table>

thermosphere-ionosphere and hydrodynamic escape for close-in exoplanets to compute the composition and temperature in the upper atmosphere (T. T. Koskinen et al. 2013a, 2022). We combine the results for the upper atmosphere with a radiative transfer-photochemistry model of the lower/middle atmosphere (P. Lavvas & A. Arfaux 2021) to simulate the full atmospheric structure of HD209458b. We systematically test and validate our model in the following steps: (1) We first reproduce the results of an isothermal Parker wind model for direct comparison with past studies based on this approach. (2) We then introduce our full thermosphere-ionosphere model, which accounts for nonisothermal temperature structures and diffusive separation. (3) We incorporate the lower/middle atmosphere model to fully couple the deeper layers with the escaping atmosphere, and present our best-fit model to the He I 10830 Å and Hα transit depth observations. (4) We analyze the impact of stellar activity variations, molecular and eddy diffusion, tidal effects, and the incorporation of metals on atmospheric structure, He I 10830 Å absorption, and Hα transit depths. This methodology enables us to separate the effects of different physical processes and determine which factors most significantly influence transit observations. Our results confirm that stellar activity variations, molecular diffusion, heating structure, and tidal effects all play a role in shaping He I 10830 Å absorption. We show that multiple competing processes—beyond just XUV flux scaling—control He I 10830 Å transit depths, likely explaining why some exoplanets appear as outliers in population studies.

## 2. Methods

Here, we summarize our methods to simulate the atmosphere and calculate transit depths. The stellar and planetary parameters we use for the HD209458b system are given in Table 1. Section 2.1 describes the lower/middle atmosphere model that we use to set the lower boundary conditions for the upper atmosphere. In Section 2.2, we describe the multispecies model of hydrodynamic escape that we use to model the upper atmosphere. Section 2.3 describes the inputs we use to model HD209458's stellar flux. In Section 2.4, we describe the updates to the He I (2³S) excitation/de-excitation rates. In Section 2.5, we describe the detailed balance model used to compute the excited state H populations. Finally, In Section 2.6, we describe the methods that we use to simulate transit depths based on the output of the lower/middle atmosphere model and the upper atmosphere model.

### 2.1. Lower and Middle Atmosphere

Consistency with the properties of the lower and middle atmosphere can be important to correctly simulating the upper atmosphere and escape (e.g., J. M. Chadney et al. 2017;

C. Huang et al. 2023). It is also required to connect the transit depths probing the upper atmosphere to the transit depths at adjacent wavelengths. To achieve this, we employ the lower and middle atmosphere model (P. Lavvas & A. Arfaux 2021) that computes the composition and temperature at pressures ranging from 1000 bar to $10^{-6}$ bar where the latter pressure is the lower boundary of our upper atmosphere model. This model calculates photochemistry self-consistently with the temperature profile, encompassing approximately 150 molecules, atoms, and ions with more than 1600 reactions. Based on the predicted atmospheric structure and composition, the model also calculates the theoretical transit spectrum of the planet, accounting for Rayleigh scattering, gaseous absorption, and Mie scattering by aerosols when applicable. Key contributors include H₂, H₂O, CH₄, NH₃, HCN, H₂S, CO₂, Na, K, collision-induced absorptions by H₂ – H₂ and He – H₂, and many metal atomic lines. Opacity sources follow P. Lavvas & A. Arfaux (2021). For a more detailed model description, see A. Arfaux & P. Lavvas (2022).

The lower and middle atmosphere model sets the lower boundary temperature and altitude for our upper atmosphere escape model at the pressure of 1 µbar. The transit “continuum” around the He I 10830 Å and Hα lines generally probes higher pressures at altitudes below the thermosphere and tends to be set by water absorption as well as extinction by clouds and hazes. Our model includes aerosol extinction due to photochemical hazes, based on a production rate of $1 \times 10^{-14} \text{ g cm}^{-2} \text{ s}^{-1}$ in the upper atmosphere that improves the fit to the observed transmission spectrum (see Section 3.3). It is also necessary to extend the level population calculations below the lower boundary of the upper atmosphere model. To calculate the density of metastable helium in the lower and middle atmosphere, we assume a balance between production by recombination and loss to radiative decay and Penning ionization, guided by the results of our detailed balance model (see Section 2.4). We investigated other excitation mechanisms that might be important for populating excited helium in the lower/middle atmosphere, including Penning ionization with H₂, but find that these rates are not high enough to make a significant difference in the metastable helium population in the lower/middle atmosphere. For excited H, we assume LTE in the lower/middle atmosphere, multiplying the neutral hydrogen density by a Boltzmann distribution based on the temperature of the atmosphere. Both of these methods show smooth density transitions from the lower/middle atmosphere to the upper atmosphere.

### 2.2. Upper Atmosphere Escape Model

We use a 1D model of the thermosphere-ionosphere and hydrodynamic escape for close-in exoplanets (T. T. Koskinen et al. 2013a, 2013b). It solves the time-dependent coupled equations of continuity, momentum, and energy in the radial direction. The model includes photoionization by stellar XUV radiation, thermal ionization, recombination, charge exchange, advection and thermal escape of ions and neutrals, diffusion, viscous drag, heating by photoionization, adiabatic heating and cooling, heat conduction, and viscous dissipation. As we explain below, we only include H and He, their ions, and electrons in the reference simulations for this work, simplifying the model to highlight the impact of atmospheric structure and revised excitation and de-excitation rates on He I 10830 Å transit depths. For radiative cooling, we include recombination
<table><caption>**Table 2**

He(2³S) Chemistry Reaction Rate Coefficients</caption><thead><tr><th>Reaction</th><th>Rate<br>(cm<sup>3</sup> s<sup>-1</sup>)</th><th>References</th></tr></thead><tbody><tr><td>H + hν → H<sup>+</sup> + e<sup>-</sup></td><td>SC</td><td>D. G. Hummer (1962)</td></tr><tr><td>He + hν → He<sup>+</sup> + e<sup>-</sup></td><td>SC</td><td>M. Yan et al. (1998)</td></tr><tr><td>H<sup>+</sup> + e<sup>-</sup> → H + hν</td><td>4.0 × 10<sup>-12</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>P. J. Storey & D. G. Hummer (1995)</td></tr><tr><td>He<sup>+</sup> + e<sup>-</sup> → He + hν</td><td>1.48 × 10<sup>-12</sup>(T/300)<sup>-0.54</sup></td><td>R. A. Benjamin et al. (1999)</td></tr><tr><td>H + e<sup>-</sup> → H<sup>+</sup> + e<sup>-</sup> + e<sup>-</sup></td><td>2.91 × 10<sup>-8</sup>U<sup>0.39</sup>exp(-U)<sup>0.232 + U</sup>, U = 13.6/E<sub>e</sub>(eV)</td><td>G. S. Voronov (1997)</td></tr><tr><td>He + e<sup>-</sup> → He<sup>+</sup> + e<sup>-</sup> + e<sup>-</sup></td><td>1.75 × 10<sup>-8</sup>U<sup>0.35</sup>exp(-U)<sup>0.180 + U</sup>, U = 24.6/E<sub>e</sub>(eV)</td><td>G. S. Voronov (1997)</td></tr><tr><td>H + He<sup>+</sup> → H<sup>+</sup> + He</td><td>1.30 × 10<sup>-15</sup>(300/T<sub>e</sub>)<sup>-0.25</sup></td><td>P. C. Stancil et al. (1998)</td></tr><tr><td>He + H<sup>+</sup> → He<sup>+</sup> + H</td><td>1.75 × 10<sup>-11</sup>(300/T<sub>e</sub>)<sup>0.75</sup>exp(-127500/T)</td><td>S. C. O. Glover & A.-K. Jappsen (2007)</td></tr><tr><td>He<sup>+</sup> + e<sup>-</sup> → He<sup>+</sup> + e<sup>-</sup> + e<sup>-</sup></td><td>2.05 × 10<sup>-7</sup>U<sup>0.25</sup>exp(-U)<sup>1 + U<sup>0.5</sup>0.265 + U</sup>, U = 52.4/E<sub>e</sub>(eV)</td><td>G. S. Voronov (1997)</td></tr><tr><td>He<sup>2+</sup> + e<sup>-</sup> → He<sup>+</sup></td><td>1.92 × 10<sup>-11</sup>(300/T<sub>e</sub>)<sup>0.61</sup></td><td>D. A. Verner & G. J. Ferland (1996)</td></tr><tr><td>He<sup>2+</sup> + H → He<sup>+</sup> + H<sup>+</sup></td><td>1 × 10<sup>-14</sup></td><td>J. B. Kingdon & G. J. Ferland (1996)</td></tr><tr><td>He<sup>+</sup> + hν → He<sup>+</sup> + e<sup>-</sup></td><td>SC</td><td>D. A. Verner et al. (1996)</td></tr><tr><td>He(1'S) + e<sup>-</sup> → He(2³S) + e<sup>-</sup></td><td>SC</td><td>This work</td></tr><tr><td>He<sup>+</sup> + e<sup>-</sup> → He(2³S) + hν</td><td>2.10 × 10<sup>-13</sup>(T<sub>e</sub>/10<sup>4</sup>)<sup>-0.778</sup></td><td>R. A. Benjamin et al. (1999)</td></tr><tr><td>He(1'S) + e<sup>-</sup> → He(2³S) + e<sup>-</sup></td><td>2.1 × 10<sup>-8</sup>√(13.6/kT)exp(-19.81/kT)γ<sub>13</sub></td><td>I. Bray et al. (2000)</td></tr><tr><td>He(1'S) → He(2³S)</td><td>1.272 × 10<sup>-4</sup></td><td>G. W. F. Drake (1971)</td></tr><tr><td>He(2³S) + hν → He<sup>+</sup> + e<sup>-</sup></td><td>SC</td><td>This Work</td></tr><tr><td>He(2³S) + H → He(1'S) + H</td><td>T ≤ 4000 K: 1.9 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.07</sup></td><td>This Work</td></tr><tr><td>...</td><td>T > 4000 K: 9.1 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.50</sup></td><td>...</td></tr><tr><td>He(2³S) + e<sup>-</sup> → He(2¹'S) + e<sup>-</sup></td><td>2.1 × 10<sup>-8</sup>√(13.6/kT)exp(-0.80/k<sub>B</sub>T)<sup>γ<sub>13</sub>α/3</sup></td><td>I. Bray et al. (2000)</td></tr><tr><td>He(2³S) + e<sup>-</sup> → He(2¹'P) + e<sup>-</sup></td><td>2.1 × 10<sup>-8</sup>√(13.6/kT)exp(-1.40/k<sub>B</sub>T)<sup>γ<sub>13</sub>b/3</sup></td><td>I. Bray et al. (2000)</td></tr><tr><td>He(2³S) + e<sup>-</sup> → He + e<sup>-</sup></td><td>2.60 × 10<sup>-7</sup>(T<sub>e</sub>)<sup>-0.531</sup>exp(-5.43 × 10<sup>2</sup>/T<sub>e</sub>)</td><td>M. Schulik & J. Owen (2025)</td></tr><tr><td>He(2³S) + H<sup>+</sup> → He + H<sup>+</sup></td><td>5.55 × 10<sup>-24</sup>T<sub>e</sub><sup>-0.5</sup></td><td>L. Augustovičová et al. (2014)</td></tr></tbody></table>

Note. SC: Calculated self-consistently based on the model's ionization cross sections, stellar flux, and density profiles.

and H I line cooling, as described by C. Huang et al. (2023). Photochemistry and helium excitation/de-excitation are treated by the flexible open-source Kinetic Preprocessor KPP-3.0.0 chemical model that is fully coupled to the escape model (A. Sandu & R. Sander 2006). For completeness, further details of the model are given in the Appendix of T. T. Koskinen et al. (2022).

We consider two classes of upper atmosphere models in this work. Most of our simulations focus on hydrogen and helium only, including the species H($n=2$), He, H$^+$, He$^+$, He$^{2+}$, He(2³S), and electrons. These models allow for direct comparisons with many previous studies of HD209458b and isolate the effects of updated model physics and parameter variations on the He I 10830 Å transit depths. The reactions for these species are listed in Table 2, with the exception of $H(n=2)$, which is handled separately (see below). In Section 3.6, we expand on this framework by including heavy elements and ions such as O, C$^+$, Si$^{2+}$, and Fe$^+$, for which observational evidence in HD209458b's upper atmosphere has been reported (A. Vidal-Madjar et al. 2004, 2013; J. L. Linsky et al. 2010; P. E. Cubillos et al. 2020). These species can influence the energy balance and composition of the atmosphere, and we assess their impact on both He I 10830 Å and H$\alpha$ transit depths. The relevant reaction rates for metal chemistry follow the implementation described by C. Huang et al. (2023).

## 2.3. Stellar Flux

HD209458 is a G0-type inactive star with a coronal temperature comparable to that of the inactive Sun (S. Czesla

et al. 2017). The Sun's quiescent X-ray luminosity ranges from $10^{27}-10^{28}$ erg s⁻¹ (P. G. Judge et al. 2003), and the X-ray luminosity of HD209458 sits at the lower end of that range at around $1.6 \times 10^{27}$ erg s⁻¹ (S. Czesla et al. 2017). Therefore, to match observations, we use a solar minimum spectrum from 2008, accessed from LISIRD<fnref n="5" /> and scaled to the radius of HD209458 as our default input stellar spectrum. To investigate how stellar activity level impacts the He I (2³S) and H$\alpha$ Balmer line transit depths, we also employ a solar average and solar maximum spectrum, which are based on appropriate TIMED SEE spectra scaled to the radius of HD209458. Note that for each activity level, we use the LISIRD/TIMED SEE spectra at wavelengths shorter than 185 nm and a PHOENIX model ($T_{eff} = 6075$ K, log(g) = 4.38; T. O. Husser et al. 2013) at longer wavelengths. To calculate the H I level populations (see Section 2.5), the H Ly$\alpha$ flux is needed for resonant scattering calculations. For each stellar activity level, we use the solar Ly$\alpha$ profiles given by P. Lemaire et al. (2005) matched to the appropriate stellar activity.

## 2.4. Non-LTE Rate Coefficients for He I (2³S)

Building on the work presented by A. Taylor et al. (2024), we refine some of the rate coefficients that govern the population of metastable helium. The He I 10830 Å transit depth is particularly sensitive to the excitation and de-excitation processes that control the abundance of the metastable 2³S state. Among these processes, recombination of He$^+$ plays a dominant role in producing metastable helium in hot exoplanet atmospheres

<footnote n="5">⁵ https://lasp.colorado.edu/lisird/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Recombination rates to the metastable He I (2<sup>3</sup>S) triplet state vs. electron temperature from different sources.</figcaption></figure>

(S. Seager & D. D. Sasselov 2000). As shown in Figure 1, recombination rate coefficients used in some previous studies of metastable helium in different environments, including the Earth's upper atmosphere, vary by a factor of 2–3. These total recombination rates are based on the assumption that recombination from all higher-lying triplet states contributes to the population of the metastable 2<sup>3</sup>S state, given the rapid radiative decay of these states. For example, in studies of metastable helium emission in the Earth's atmosphere, L. S. Waldrop et al. (2005) used tables of partial recombination rates from D. A. Verner & G. J. Ferland (1996), while N. N. Shefov et al. (2009) appears simply to have halved the total recombination rate to estimate the recombination rate to the triplet state. Exoplanet models, on the other hand, typically use the total case B triplet recombination rate from R. A. Benjamin et al. (1999). The discrepancies seen in Figure 1 largely result from incomplete summation over the recombination pathways in the Earth models. In particular, the tables from D. A. Verner & G. J. Ferland (1996) exclude recombination to many higher-lying triplet states in the triplet branch, leading to a significantly underestimated rate. In line with previous exoplanet work, we therefore use the recombination rates from R. A. Benjamin et al. (1999). In order to calculate the recombination rate to the ground state, we assume that recombination to all singlet states rapidly decays to the ground state and thus contributes to the total ground-state recombination rate. To avoid double counting, we compute the effective ground-state recombination rate by subtracting the triplet case B rate from the total recombination rate given by R. A. Benjamin et al. (1999). The resulting rate coefficients are listed in Table 2.

Photoionization is an important loss mechanism for the 2³S
state of helium, especially at high altitudes in close-in
exoplanet atmospheres (see Section 3.4). Previous studies
have typically used a low-resolution cross section from
D. W. Norcross (1971) that does not capture strong resonances
in the EUV. We have generated a new high-resolution
photoionization cross section for the 2³S state of helium using
the B-spline K-matrix method (L. Argenti &
R. Moccia 2006, 2008; L. Argenti & E. Lindroth 2023). The
calculations were performed within a spherical quantization
box with a radius of R<sub>box</sub> = 800 Bohr radii. The radial
components of the orbitals were expanded in B-splines of

order $k = 8$ (H. Bachau et al. 2001), while spherical harmonics represented the angular components. The configuration space comprises two sectors: (1) a close-coupling sector, including all channels of the form $n\ell\epsilon_{\ell'}$ with $n \le 3$ and angular momenta $\ell, \ell'$ compatible with the overall $^3S^e$ and $^3P^o$ symmetries; and (2) a dynamic-correlation sector, composed of configuration state functions of the form $n\ell n'\ell'$ with $n, n' \le 28$ and $\ell, \ell' \le 8$, constructed from orbitals localized within 48 Bohr radii of the origin. Cross sections were computed in the velocity gauge (B. H. Bransden & C. J. Joachain 2003). We make the cross sections available in a machine-readable format with this paper and on The University of Arizona's Research Data Repository.<fnref n="6" /> Figure 2 shows the calculated cross section compared with the D. W. Norcross (1971) cross section. In the right panel of Figure 2, our high-resolution cross sections at EUV wavelengths show sharp resonance features that are orders of magnitude stronger than those in the previously used low-resolution data, areas where they overlap with strong stellar coronal emission lines.

Penning ionization with hydrogen is another important loss
mechanism for HeI (2<sup>3</sup>S), especially in the lower thermo-
sphere and below it (see Section 3.4). Previous studies have
used a temperature-independent value for the rate coefficient
(A. Oklopčić & C. M. Hirata 2018; M. Lampón et al. 2020)
based on a study on the formation and destruction of HeH<sup>+</sup> in
astrophysical plasmas (W. Roberge & A. Dalgarno 1982). In
the latter study, they find that the rate coefficient for the sum of
associative ionization reaction [He(2<sup>3</sup>S) + H → HeH<sup>+</sup> + e<sup>-</sup>]
and the Penning ionization reaction [He(2<sup>3</sup>S) + H → He(1<sup>1</sup>S)
+ H<sup>+</sup> + e<sup>-</sup>] is ~5 × 10<sup>-10</sup> cm<sup>3</sup>s<sup>-1</sup>. Cross sections are
available for the sum of the associative ionization and the
Penning ionization reaction, which can be used to calculate a
temperature-dependent reaction rate by using the Maxwell-
Boltzmann distribution. The cross sections we use for this
calculation are from H. Morgner & A. Niehaus (1979) and
J. S. Cohen & N. F. Lane (1971). We fit the resulting
temperature-dependent rate to obtain the rate coefficient in
Table 2. This revised rate coefficient leads to an updated
assessment of the relative role of Penning ionization in the
depletion of metastable helium (see the Appendix for a more
detailed discussion of these calculations).

Free electron collisions are interactions between free
electrons (those not bound to atoms or molecules) and other
particles that can transfer energy between the electron and the
particles, resulting in an excited or de-excited target particle.
For collisions with thermal electrons, we use the rate
expressions that have been referenced in other works (see
Table 2; A. Oklopčić & C. M. Hirata 2018; M. Lampón et al.
2020) with temperature-dependent oscillator strengths from
I. Bray et al. (2000). Note that the oscillator strengths (denoted
by Υ in Table 2) have a limited temperature range, with the
lowest available temperature tabulations at 6000 K. As our
temperature profiles range from 3000–9000 K, we use only the
oscillator strengths at 6000 K. The 2<sup>3</sup>S → 2<sup>1</sup>S transition due to
thermal electron collisions is an important loss mechanism for
metastable helium in the thermosphere. The value of the
oscillator strength for this transition is 2.389. We also include
the thermal electron collisions that result in the 1<sup>1</sup>S → 2<sup>3</sup>S
transition and the 2<sup>3</sup>S → 2<sup>1</sup>P transition, although these are not
significant production or loss mechanisms for metastable

<footnote n="6">⁶ https://redata.arizona.edu/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The photoionization cross section for He I (2³S) vs. wavelength. The black line indicates what has been used in previous literature (D. W. Norcross 1971), while the pink lines show high-resolution cross sections used in our model. The data behind this figure is available in machine-readable format.
(The data used to create this figure are available in the online article.)</figcaption></figure>

<figure><img src="image_5.png" /></figure>

helium. The oscillator strengths at 6000 K are 6.198 × 10⁻² for the ¹¹S → ²³S transition and 7.965 × 10⁻¹ for the ²³S → ²¹P transition. Note that since the He I (2¹S) and (2¹P) states are so short-lived, we assume that they decay to ground-state helium instantaneously.

Photoelectrons are released by photoionization events, and their typical energies far exceed those of thermal electrons. Photoelectron impact on helium is an important production mechanism of metastable helium in the Earth's atmosphere (L. S. Waldrop et al. 2005), but it has been mostly ignored in exoplanet studies thus far. Here, we estimate a rough upper limit on the contribution of photoelectron excitation to the metastable helium population, based mostly on the production of primary photoelectrons. This involves using known photoionization cross sections and stellar fluxes to calculate the energy-dependent production of primary photoelectrons throughout the upper atmosphere. Cross sections for electron impact excitation to the ²³S state are taken from laboratory measurements and fitted following the methods of L. S. Waldrop et al. (2005) and F. J. de Heer & R. H. J. Jansen (1977; see the Appendix for more details). To more realistically account for the contribution of high-energy photoelectrons, we adopt a simplified scheme for estimating the contribution by secondary electrons. Specifically, for each primary photoelectron with energy $E_p$ exceeding twice the characteristic excitation energy $E_{cs}$ of metastable helium (roughly 36 eV, where the cross section for metastable excitation peaks), we assume that $[E_p/E_{cs}] - 1$ secondary electrons are produced with energy $E_{cs}$. These secondary electrons are then added to the appropriate energy bin before computing the total excitation rate. This approach ignores the interaction of the primary photoelectrons with other particles and assumes that secondary electron energies coincide with the peak interaction probability for exciting He. We followed this simplified approach because we found that photoelectron excitation does not significantly affect the transit depth for HD209458b, even with this upper limit.

In addition to the above revisions, we have included proton de-excitation to the ground state of helium. For this process, we use cross sections from L. Augustovičová et al. (2014), and calculate a temperature-dependent reaction rate by using the Maxwell-Boltzmann distribution. We fit the resulting rate coefficients with analytical equations to derive the rate expressions shown in Table 2. This means that collisional rates

to these states are combined into one rate for collisions to the ground state. We note that the proton de-excitation rates do not significantly impact the He I (2³S) transit depth on HD209458b and are included here for completeness. The Appendix provides a more detailed discussion of the relevant calculations.

## 2.5. H I Level Populations

In addition to the He I 10830 Å line, we model the Hα absorption line of hydrogen using a previously developed non-LTE model to simulate the level populations of hydrogen (C. Huang et al. 2017, 2023). As demonstrated by previous studies (C. Huang et al. 2017; A. García Muñoz & P. C. Schneider 2019), the excitation of the $H(n=2)$ state is primarily driven by radiative excitation, with the population of $H(2p)$ largely excited by the $Lyα$ intensity in the atmosphere. In this study, we employ the $Lyα$ Monte Carlo resonant scattering radiative transfer model developed by C. Huang et al. (2017) to determine the $Lyα$ mean intensity within a plane-parallel atmosphere.

Due to the computational intensity of this calculation, it is not feasible to directly integrate the non-LTE model for H with our model of the upper atmosphere. Instead, we calculate the $Lyα$ radiation field separately based on the upper atmosphere model output and then use the results to update the $H(n=2)$ density in the upper atmosphere. Thus, the hydrodynamic model and radiation transfer calculations are performed iteratively until the electron densities and temperature profiles in the upper atmosphere model converge, typically after three iterations. This iterative process ensures that the level populations are consistent with the atmospheric conditions. We note that ionization of $H(n=2)$ is a significant source of free electrons that contributes to the heating rate of the upper atmosphere.

## 2.6. Simulated Transit Depths

To compare our atmosphere models to observations of the H I and He I 10830 Å transit depths, we calculate the transit depths based on the modeled pressure, temperature, and density profiles from the upper atmosphere and lower/middle atmosphere models. To compute the relative diminution of the flux from the star caused by the presence of the planet, we use
Equation (11) from T. M. Brown (2001):

$$ \frac{\partial F(k)}{F(k)} = \frac{1}{\pi r_*^2} \int_0^{z_{\text{max}}} 2\pi(r_p + z^*) \times [1 - \exp(-\tau(z^*, k))] dz^* \quad (1) $$

where $r_*$ is the radius of the star, $z^*$ is the ray depth, $r_p$ is the radius of the planet, $k$ is wavenumber, and $z_{\text{max}}$ is the maximum ray depth. The optical depth $\tau$ is given by

$$ \tau(z', k) = 2 \int_{0}^{x_{\text{max}}} \kappa(z[z', x]) dx \quad (2) $$

where $z'$ is the minimum height of a ray above the 1 bar pressure level, $\kappa$ is the opacity, and $x_{\text{max}}$ is the maximum line-of-sight distance in the atmosphere. We use line properties for atomic lines from the NIST database.<fnref n="7" /> Here, we only focus on the He I 10830 Å and H$\alpha$ transit depths, while the details for the lower/middle atmosphere transit depth calculation and the relevant molecular line lists are given in Section 2.1.

For this study, we calculate center-of-transit transit depths by accounting for thermal and natural broadening using a Voigt profile. We include velocity broadening by incorporating the line-of-sight outflow velocity as an additive term in the Gaussian component of the Voigt profile and add rotational broadening of 2 km s⁻¹ based on the spin of the planet. We follow this relatively simple approach to include velocity and rotational broadening to facilitate an easier comparison with previous studies and to provide a clear constraint on the required broadening to match observations, which our parameters do (see Section 3). To model the He I 10830 Å and H$\alpha$ features, we use density and temperature profiles from both the lower/middle and upper atmosphere models. We then combine the He I 10830 Å and H$\alpha$ line profiles and the lower/middle atmosphere transit spectrum, so that the excited hydrogen and helium absorption features appear as enhancements on top of the continuum transit spectrum (see Section 3.3). This approach ensures that the atomic lines are treated consistently with the atmospheric structure and that the resulting transit depths are normalized to the expected continuum level from the lower atmosphere.

We note that transit depths are also influenced by factors such as 3D dynamics, possible planetary magnetic fields, and interactions with stellar winds that are not included in our 1D model (J. E. Owen et al. 2020). One-dimensional models can, however, be used to explore the physics and processes in detailed balance models that help to identify key inputs to more complex models. One-dimensional models are also much faster and facilitate population-level surveys of broad theoretical trends. Here, we also note recent studies showing that a hotter dayside means that transmission spectra are more reflective of the average conditions on the dayside than along the terminator (A. Y. Jaziri et al. 2024). This is especially true of hot, extended upper atmospheres because their large extent means that radiation on close-in exoplanets can penetrate past the terminator to the night side, suggesting that a significant fraction of the upper atmosphere exhibits dayside-like behavior (T. T. Koskinen et al. 2007).

### 3. Results

Here, we present results from our atmosphere models for HD209458b, focusing on critical parameters such as temperature and velocity distributions, atmospheric composition, metastable helium production and loss processes, and transport mechanisms. We begin by examining a baseline model configured to replicate the framework of the isothermal Parker wind model, characterized by an isothermal temperature profile and the absence of diffusion processes, using this setup to fit existing observational data. We then introduce our upper atmosphere model where the temperature profile is self-consistently calculated from the energy equation rather than assuming an isothermal structure. Building upon this, we integrate the lower and middle atmosphere model to produce a fully coupled atmosphere model, enabling direct comparisons with observational constraints from HST and JWST. Then, we present our best-fit model for the transit depths of He I 10830 Å and H$\alpha$, which constrain the photoelectron heating efficiency and mass-loss rate. We also examine how key parameter variations—including the stellar activity level, molecular and eddy diffusion, and tidal forcing—shape the thermal structure, composition, and transit depths predicted by our models. Finally, we incorporate heavy elements into the upper atmosphere to investigate their impact on the atmospheric structure and observable transit signatures.

#### 3.1. Isothermal Upper Atmosphere Model

To establish a foundation for our analysis, we first use our model to reproduce the approach of previous studies, such as those utilizing the isothermal p-winds framework and the modeling techniques presented by L. A. Dos Santos et al. (2022), M. Lampón et al. (2020), and A. Oklopčić & C. M. Hirata (2018). This baseline model, which we call Model A (see Table 3 for a summary of model configurations), assumes a constant temperature of 8500 K with altitude, bulk transport with no diffusion processes, and the updated excitation and de-excitation rates for He I (2$^3$S) as described in Section 2.4. We model radiative transfer in the atmosphere by using a solar zenith angle of 60° for a dayside average and further dividing the incident stellar flux by a factor of 2 to approximate global redistribution (F. L. Smith & C. Smith 1972) that may be more suitable for the terminator that is probed in transit than the overhead illumination adopted in several other 1D studies. Unlike L. Xing et al. (2023), we do not include the substellar tide by default because it does not apply at the terminator (see Section 3.5 for further discussion on this).

We set the temperature to 8500 K because this is at the lower end of temperatures used to match the FUV transit observations (T. T. Koskinen et al. 2010) and because it is consistent with our best-fit full atmosphere model. Using this temperature to match the He I 10830 Å observations (F. J. Alonso-Floriano et al. 2019), we determine a mass-loss rate of $1.86 \times 10^{10}$ g s⁻¹ that is consistent with best-fit models from M. Lampón et al. (2020). We note that temperature and mass-loss rate are degenerate in isothermal fits like this (see below). Figure 3 presents the temperature and velocity profiles (left panel), and the composition of the atmosphere (right panel). Figure 4 shows the resulting transit depth compared with observations for the He I 10830 Å line (left panel) and the H$\alpha$ line (right panel). We explored the effect of diffusion

<footnote n="7">⁷ https://physics.nist.gov/PhysRefData/ASD/lines_form.html</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Isothermal model (Model A) results for HD209458b. The left panel shows the temperature (pink) and velocity (orange) profiles as a function of planetary radii, with a mass-loss rate of $1.86 \times 10^{10}$ g s$^{-1}$. The right panel presents the number densities of key atmospheric species, including neutral hydrogen (H), ionized hydrogen (H$^+$), neutral helium (He), ionized helium (He$^+$), doubly ionized helium (He$^{2+}$), metastable helium He I (2$^3$S), and free electrons.</figcaption></figure>

<table><caption>**Table 3**
Summary of Model Configurations</caption><thead><tr><th>Feature</th><th>Model A</th><th>Model B</th><th>Model C</th><th>Best-fit Model</th><th>Metals Model</th></tr></thead><tbody><tr><td>Name</td><td>Isothermal Model<br>...</td><td>Self-consistent<br>Upper Atmosphere</td><td>Full Atmosphere Model<br>...</td><td>Best-Fit Model<br>...</td><td>Metal Run<br>...</td></tr><tr><td>Temperature Structure</td><td>Fixed (8500 K)<br>...</br>...</td><td>Energy equation<br>...</br>...</td><td>Energy equation<br>and lower/middle<br>atmosphere</td><td>Energy equation<br>and lower/middle<br>atmosphere</td><td>Energy equation<br>and lower/middle<br>atmosphere</td></tr><tr><td>Photoelectron Heating<br>Efficiency</td><td>N/A<br>...</td><td>100%<br>...</td><td>100%<br>...</td><td>40%<br>...</td><td>23%<br>...</td></tr><tr><td>Lower Boundary Pressure<br>Upper Atmosphere Model</td><td>7 nbar<br>...</td><td>10<sup>-6</sup> bar<br>...</td><td>10<sup>-6</sup> bar<br>...</td><td>10<sup>-6</sup> bar<br>...</td><td>10<sup>-6</sup> bar<br>...</td></tr><tr><td>Lower Boundary Altitude<br>Upper Atmosphere Model</td><td>0 km<br>...</td><td>0 km<br>...</td><td>5152 km<br>...</td><td>0 km<br>...</td><td>0 km<br>...</td></tr><tr><td>Diffusion Processes</td><td>None</td><td>Multispecies</td><td>Multispecies</td><td>Multispecies</td><td>Multispecies</td></tr><tr><td>Revised Rates</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>He/H Ratio</td><td>Fixed (92/8)</td><td>Evolving</td><td>Evolving</td><td>Evolving</td><td>Evolving</td></tr><tr><td>Lower/Middle Atmosphere</td><td>Not included</td><td>Not included</td><td>Included</td><td>Included</td><td>Included</td></tr><tr><td>Metals Included</td><td>No</td><td>No</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>Mass-loss Rate (g s<sup>-1</sup>)</td><td>1.9 &times; 10<sup>10</sup></td><td>3.5 &times; 10<sup>10</sup></td><td>5.5 &times; 10<sup>10</sup></td><td>1.9 &times; 10<sup>10</sup></td><td>3.0 &times; 10<sup>10</sup></td></tr></tbody></table>

processes on Model A, but since the lower boundary pressure of 7 nbar in this model is relatively low, and the temperature is high everywhere, the outflow velocities are substantial even from the bottom boundary, and we see no effect from diffusive separation as a result.

While Model A provides a reasonable fit to the data, real atmospheres are not isothermal. The assumption of a constant temperature neglects key physical processes such as the dependency of the radiative heating and cooling rates on altitude, adiabatic expansion, and heat conduction that shape the actual thermal structure of the atmosphere. In addition, the choice of the mass-loss rate and temperature in the isothermal Parker wind model dictates the lower boundary pressure, which, in our best-fit case, is 7 nbar. This dependence arises because the Parker wind model enforces a specific relationship between temperature, velocity, and pressure, which constrains

the lower boundary conditions in a way that does not necessarily reflect actual atmospheric structure. To understand this constraint, we briefly outline the fundamental equations governing the solution. The basic Parker wind model uses a set of equations to describe the hydrodynamic escape of an atmosphere, focusing on mass and momentum conservation. For an isothermal wind, the momentum equation and mass conservation can be expressed as:

$$ \frac{d\psi}{d\xi}\left(1 - \frac{\tau}{\psi}\right) = -2\xi^2 \frac{d}{d\xi}\left(\frac{\tau}{\xi^2}\right) - \frac{2\lambda}{\xi^2} \quad (3) $$

where $\psi = \frac{1}{2} \frac{m_i v^2}{kT_0}$ is a dimensionless measure of kinetic energy, $\xi = \frac{r}{r_0}$ is the dimensionless radial distance, $\tau = \frac{T(r)}{T_0}$ represents the temperature profile, and $\lambda = \frac{GM_s m_i}{2r_0 k T_0}$ is the thermal escape
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Transit depth results from our baseline model (Model A), comparing the modeled He I 10830 Å (left) and H I Balmer α (right) transit depths from our baseline model (orange solid lines) with the best-fit model to the observed transit depth for the He I 10830 Å line from F. J. Alonso-Floriano et al. (2019; blue line), respectively, demonstrating that the baseline model reasonably captures the overall absorption features.</figcaption></figure>

parameter. Parker's model assumes that the temperature remains constant (isothermal condition) up to a certain radius $r_b$ beyond which the temperature drops to negligible levels, although this implicitly assumed decrease in temperature is not included in the models. The transonic Parker wind solution, assumed by the above models, includes a critical point ($\xi_c, \psi_c$) where $\xi_c = \lambda/2$ and $\frac{d\psi}{d\xi} = 0 \rightarrow \psi_c = 1$. As was pointed out by J. W. Chamberlain (1960), for a given escape parameter $\lambda$ (i.e., temperature), there is only one unique transonic solution (i.e., only one value of the lower boundary kinetic energy parameter $\psi_0$ that will provide the transonic solution). This means that the velocity at the lower boundary depends on the temperature and cannot vary. The lower boundary density and therefore pressure depend on the assumed mass-loss rate, given by

$$ \dot{M} = 4\pi r_0^2 m_i n_0 v_0. \quad (4) $$

This naturally leads to a degeneracy between mass-loss rate and temperature when Parker wind models are fitted to transit observations. In other words, setting the mass-loss rate $\dot{M}$ in Parker wind models only changes the density/pressure profile rather than the velocity profile if the temperature is held fixed. This allows for arbitrary changes in density/pressure that may not be consistent with the energy balance of the upper atmosphere. This is because the parameter $\lambda$ dictates the flow's velocity profile, but this parameter has no dependence on density. Thus, changing the mass-loss rate does not alter the velocity profile that is governed by $\psi(\xi)$ and the critical point condition. In models such as the isothermal p-winds (L. A. Dos Santos et al. 2022), the user therefore selects the temperature and mass-loss rate to match observations, and this selection fixes the lower boundary pressure. We emphasize that the lower boundary density should be set by the conditions at the actual wind launching radius, which depend on the detailed properties of the atmosphere.

The lower boundary pressure of 7 nbar in our best-fit isothermal model is located in the middle of the photoionization region for He (where the helium photoionization rate peaks) that is ionized deeper than H due to its higher ionization

potential and the majority of He ionization occurs at altitudes near and above the lower boundary rather than higher up in the atmosphere. This artificially enhances the ionized helium, and therefore the metastable helium density at lower radii (see Section 3.2) while suppressing it in the extended upper atmosphere, leading to a reduced He I transit depth compared to a self-consistent thermal model with a more appropriate lower boundary pressure below the ionization peak. The Hα differential transit depth is also lower, partly because the transit continuum is set to the 7 nbar level while the true continuum is much deeper.

## 3.2. Self-consistent Upper Atmosphere Model

Motivated by these considerations, we move beyond the isothermal approximation by solving the full energy equation. We also include multispecies transport with molecular and eddy diffusion. This Model B (see Table 3 for a summary of model configurations) is initialized with a lower boundary pressure of $10^{-6}$ bar, a more physically motivated choice that coincides with the transition between the middle and upper atmosphere. In addition, previous studies fitting the H Lyα transit of HD209458b required a lower boundary pressure near $10^{-6}$ bar, above which temperature increases with altitude, to achieve a match with observations (T. T. Koskinen et al. 2010, 2013a). Model B does not yet include the lower and middle atmosphere. Instead, the lower boundary altitude is 0 km with a temperature of 1270 K. The model assumes a photoelectron heating efficiency of 100%, based on some previous studies (e.g., R. A. Murray-Clay et al. 2009), and an eddy diffusion coefficient of $10^5 \text{ m}^2\text{s}^{-1}$ (A. Arfaux & P. Lavvas 2022). The left panel of Figure 5 shows the temperature and bulk outflow profiles for Model B. The temperature peaks at ~9500 K where $r = 1.65 R_p$ ($p = 1.45 \times 10^{-10}$ bar). The energy deposition rate peaks at a radius of $1.18 R_p$ and pressure of $3.94 \times 10^{-9}$ bar. The outflow velocity is $3.6 \text{ km s}^{-1}$ at $4 R_p$ (the approximate Roche lobe radius). We note that the altitude of the sonic point is strongly dependent on the details of the temperature profile (T. T. Koskinen et al. 2013a). It also depends on Roche lobe effects that we do not include in our simulations here. If we include tidal forcing and run our model at the substellar point,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left panel: temperature and bulk outflow velocity predicted by our full upper atmosphere model for HD209458b (Model B). Right panel: composition predicted for HD209458b by Model B. Note that the H+ profile is practically identical to the electron (e⁻) density profile.</figcaption></figure>

we find a sonic point of 4.4 $R_p$, consistent with past models (L. Xing et al. 2023). We find a steady-state mass-loss rate ($\dot{M} = 4\pi v\rho r^2$) of $3.5 \times 10^{10} \text{ g s}^{-1}$. We calculate a pressure-averaged temperature of the hydrodynamic model above the 7 nbar to the top of the model at $3 \times 10^{-4} \text{ nbar}$ to compare with the isothermal model:

$$ \bar{T}_p = \frac{\int_{P_1}^{P_2} T(P)d(\ln p)}{\ln(P_2/P_1)}, \qquad (5) $$

where $P_1$ and $P_2$ are the pressure bounds for the integration. Using this method, we obtain a temperature of 8860 K. The right panel of Figure 5 shows the composition of major species in Model B. Neutral helium and hydrogen dominate toward the base of the atmosphere, but as photoionization becomes stronger at higher altitudes, hydrogen and helium ions dominate over their neutral counterparts. The photoionization timescale becomes shorter than the transport timescale at a radius of around 2 $R_p$, allowing the ions to begin to dominate over the neutral atoms. The He I (2³S) density is significantly lower than that of the dominant species and follows the general shape of the He⁺ profile due to recombination being the dominant production process.

M. Lampón et al. (2020) argued for an He/H ratio lower than 90/10 in the upper atmosphere of HD209458b, possibly due to diffusive separation, and we confirm that diffusive separation is possible. Our model produces similar results to those of L. Xing et al. (2023) and M. Schulik & J. Owen (2025), who solved separate momentum equations for different species, showing that the diffusion approximation to multi-species transport that we use works reasonably well. Figure 6 shows the hydrogen-helium ratio in Model B decreasing from 92/8 at the base of the thermosphere to a ratio of about 96/4 at high altitudes. The strongest separation happens around 1.05 $R_p$, which is also where the helium diffusion timescale is shorter than or comparable to the advection timescale and where the helium molecular diffusion coefficient becomes larger than the eddy diffusion coefficient (see the middle and right panels of Figure 6). We note that this result depends on the stellar XUV flux and the assumed eddy mixing coefficient,

where the latter was not included by L. Xing et al. (2023) or M. Schulik & J. Owen (2025). A higher stellar flux would lead to a higher mass-loss rate, while a higher eddy diffusion coefficient would lead to more efficient mixing of He to altitudes where it can more easily escape. In both cases, the effect of diffusive separation would become less pronounced.

The differences between the isothermal (Model A) and self-consistent models (Model B) become evident when comparing the He I (2³S) densities and transit depths (Figure 7). The peak in excited helium density is lower in altitude but higher in magnitude in Model A compared to Model B. Since atmospheric extent has changed due to the change in lower boundary pressure, the peak in the excited helium density moved to higher altitudes producing a higher transit depth. While the mass-loss rate in Model B is less than a factor of 2 higher than that of the Model A, the two models produce different transit depths, as shown in the middle and right panels of Figure 7. Additionally, the pressure-averaged temperature above the 7 nbar lower boundary to the top of the model at $3 \times 10^{-4} \text{ nbar}$ in Model B (8860 K) is comparable to the constant temperature used in the isothermal model (8500 K) reinforcing that differences in atmospheric structure—rather than temperature or mass-loss rate alone—can also drive differences in He I 10830 Å transit depths. Furthermore, the observed He I and $H\alpha$ transit depths are inconsistent with a photoelectron heating efficiency of 100%, suggesting that lower efficiencies (20%–40%) are required to reproduce observations.

### 3.3. Full Atmosphere Model

Building on the upper atmosphere model, we now incorporate the lower and middle atmosphere to construct a coupled model of HD209458b's vertical atmospheric structure, which we call Model C (see Table 3 for a summary of model configurations). By including the lower/middle atmosphere, we aim to assess its influence on the He I 10830 Å and $H\alpha$ transit depths and compare the resulting visible to infrared spectrum predicted by the full model with observations from HST (D. K. Sing et al. 2016) and JWST (Q. Xue et al. 2024), in addition to exploring the He I 10830 Å and $H\alpha$ line observations probing the upper atmosphere. By combining
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The elemental H/He ratio (left panel) and transport and diffusion timescales for helium (middle panel) and the eddy vs. molecular diffusion coefficients (right panel) for our upper atmosphere model (Model B).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left panel: excited helium number density comparison between our isothermal model (Model A) with no multispecies processes in orange and our upper atmosphere model (Model B) in teal. Middle panel: predicted transit depth for He I (2^3S) for our upper atmosphere model (Model B) in teal, compared to the best-fit fit line to the observations (F. J. Alonso-Floriano et al. 2019) in blue and the isothermal model in orange (Model A). Right panel: predicted transit depth for Hα for our upper atmosphere model (Model B) in teal, compared to the isothermal model (Model A) in orange.</figcaption></figure>

the upper atmosphere transit depths with the full spectrum, we obtain a self-consistent normalization for model differential transit depths that are compared with the observations (e.g., Figure 7), instead of simply normalizing the model transit depths with the assumption that the atmosphere is fully opaque below the lower boundary of the upper atmosphere model with negligible extent. Despite their critical role in regulating energy deposition and escape processes (J. M. Chadney et al. 2017; C. Huang et al. 2023), lower/middle atmosphere models have typically not been included in previous studies of He I 10830 Å absorption, potentially leading to incomplete or biased interpretations of exoplanet outflows.

The inclusion of the lower and middle atmosphere alters the overall atmospheric structure, as shown in Figure 8. The left panel presents the temperature-pressure profile of the full atmosphere model, illustrating a well-defined lower atmosphere that smoothly transitions into the thermosphere. In Model C, the total mass-loss rate is found to be 5.5 × 10^10 g s⁻¹, compared to 3.5 × 10^10 g s⁻¹ in Model B. Additionally, we compute the pressure-averaged temperature above the 7 nbar level and the 1 µbar level (up to the top of the model at 3 × 10^⁻⁴ nbar) for both models. In Model B, the pressure-averaged temperature above 7 nbar is found to be T = 8860 K, while in Model C, it is slightly lower at T = 8640 K. Similarly, above 1 µbar, the pressure-averaged temperature is T = 6040 K for Model B and T = 5610 K for Model C. Figure 8 compares the He I 10830 Å transit depth (middle panel) and H I Balmer α transit depth (right panel) for Models

B and C. Model B (teal) predicts a slightly stronger He I 10830 Å and Hα absorption signal compared to Model C (pink). The decreased He I 10830 Å and Hα transit depths result from the change in lower boundary altitude and modified velocity structure, which reduce the excited hydrogen and helium abundances at higher altitudes. While the inclusion of the lower and middle atmosphere modifies the velocity and temperature structure and increases the total mass-loss rate by about 60%, the predicted He I 10830 Å and Hα transit depths remain relatively similar. This suggests that coupling to the lower/middle atmosphere provides improvements in consistency and mass-loss estimates but does not dramatically change the modeled transit signatures of HD209458b. We also emphasize that the conclusions drawn in the previous subsection regarding the Parker wind models are largely independent of assumptions made about the lower and middle atmosphere and apply as long as the upper atmosphere is modeled self-consistently above the 1 µbar level.

Since our atmosphere model includes both the upper and lower/middle atmosphere, we can generate a complete transmission spectrum for HD209458b. Figure 9 presents the full transmission spectrum derived from the lower/middle atmosphere model, overlaid with observational data from HST and JWST. The lower/middle atmosphere model establishes the transit “continuum” around key atomic absorption lines, providing a reference for differential transit depth calculations. This spectrum captures several prominent spectral features, including the Na D and K I lines originating in the lower/middle
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left panel: P-T profile of our full atmosphere model (Model C). Middle panel: transmission spectrum around the He I triplet. In teal is the self-consistent upper atmosphere model (Model B), and in pink is the full atmosphere model (Model C). Right panel: transmission spectrum around the H I Balmer α line. In teal is the self-consistent upper atmosphere model (Model B), and in pink is the full atmosphere model (Model C).</figcaption></figure>

atmosphere, as well as the He I 10830 Å triplet and H I Balmer α line, which probe the extended escaping atmosphere. Note that the Na D and K I lines are a result of only the lower/middle atmosphere model, and we do not calculate transit depths for these lines based on the upper atmosphere model. The model successfully reproduces the continuum shape across a broad wavelength range, with the HST and JWST data aligning well with the modeled continuum, particularly in the infrared, further validating the atmospheric structure used in this study. To accurately compare our model predictions with observations from HST and JWST, we accounted for a systematic offset between the transit depths reported in the literature. Previous analyses have noted discrepancies between transit depth measurements, but the origin of these offsets remains unclear. We derived the optimal shift using the Bayesian information criterion (BIC; R. E. Kass & A. E. Raftery 1995) to determine the most likely alignment between the data sets. We define the optimal shift Δd as the value that minimizes the BIC when applied to the JWST transit depths. The BIC is given by

$$ \text{BIC} = k \log(n) - 2 \log(\mathcal{L}), \quad (6) $$

where k = 1 is the number of free parameters (the shift), n is the number of observed data points, and L is the likelihood function given by

$$ \log(\mathcal{L}) = -\frac{1}{2} \sum_i \left[ \left( \frac{T_{o,i} - T_{m,i}}{\Delta T_{o,i}} \right) + \log(2\pi\Delta T_{o,i}^2) \right] \quad (7) $$

where T_m is the model transit depth, T_o is the observed transit depth, and its uncertainty ΔT_o. We performed a minimization of the BIC as a function of the shift parameter, separately for the JWST Eureka! and Sparta data sets. We found optimal shift values of Δd = 0.000169 for the JWST Eureka! data and Δd = 0.000179 for the JWST Sparta data. These values were applied to the JWST transit depths in our analysis to align them with the HST observations and our model predictions.

### 3.4. Best-fit Model

After exploring the impact of different assumptions about atmospheric structure, we now present our best-fit H/He model. The free parameters of this model include the photoelectron heating efficiency, the eddy diffusion coefficient, and the activity level of the star, none of which are well constrained for this system. We find that achieving a fit to the

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Transit depth vs. wavelength for our full atmosphere model of HD209458b (Model C). We zoom in past the extent of the Hα and He I 10830 Å line to show the comparison with HST observations (D. K. Sing et al. 2016) and JWST observations (with different data reduction pipelines in different colors; Q. Xue et al. 2024). The JWST observations analyzed can be accessed via DOI: 10.17909/h04g-kt93. Note that the Na D and K I lines are a result of only the lower/middle atmosphere model, and we do not calculate transit depths for these lines based on the upper atmosphere model.</figcaption></figure>

He I 10830 Å transit depth that does not violate the upper limit on the Hα transit depth set by A. G. Jensen et al. (2012) requires a photoelectron heating efficiency of 40% or lower, leading to a reduced mass-loss rate compared to our previous estimates. The left panel of Figure 10 shows the temperature and velocity profiles for the best-fit model. The temperature peaks at approximately 8200 K near 1.7 R_p, lower than in Model C of Figure 5 due to the reduced heating efficiency. The outflow velocity reaches 1.1 km s⁻¹ at 4 R_p, indicating a slower, yet still significant, atmospheric escape compared to the run with 100% photoelectron heating efficiency. The corresponding mass-loss rate is 1.9 × 10¹⁰ g s⁻¹, which is lower than in Model C due to the reduced energy input but still within a physically plausible range for HD209458b.

The pressure-averaged temperature between 1 µbar and the top of the model at a pressure of 3 × 10⁻⁴ nbar is 5336 K, and above 7 nbar it is 7970 K. Somewhat by design, the temperature and mass-loss rate in our best-fit model are similar to those derived from the isothermal Parker wind solution (Model A). This indicates that it is possible for the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Left panel: temperature and bulk outflow velocity predicted by our best-fit atmosphere model for HD209458b. Middle panel: composition predicted for HD209458b. Right panel: the elemental H/He ratio for our best-fit model.</figcaption></figure>

isothermal model to produce a fit that is close to a more self-consistent solution. The isothermal Parker wind framework, however, yields a family of solutions in which mass-loss rate, temperature, and helium abundance can be freely adjusted to match observations, without enforcing a physically consistent atmospheric structure. In contrast, our self-consistent model determines a solution by coupling energy balance, multi-species transport, and radiative processes, ensuring that the temperature, velocity, and density profiles emerge from the underlying physics rather than being arbitrarily prescribed. Additionally, our approach allows us to estimate the photo-electron heating efficiency by matching the observed transit depth, which is particularly valuable given the challenges and uncertainties associated with calculating such heating efficiencies from first principles in exoplanetary atmospheres. The composition of the escaping atmosphere is shown in the middle panel of Figure 10, highlighting the number densities of key species. The distribution of neutral and ionized species follows trends consistent with our previous models, with stronger diffusive separation of helium occurring in the right panel of Figure 10. The elemental He/H ratio decreases from 8% at the base of the thermosphere to approximately 2.5% at high altitudes. This enhanced separation results from the lower mass-loss rate allowing molecular diffusion to dominate over a larger region of the atmosphere.

The loss and production rates relevant to He I (2³S) are presented in Figure 11, respectively. Recombination remains the dominant production mechanism at all altitudes. Although photoelectron collisions can be important near the lower boundary of our model, they do not significantly affect the He I 10830 Å transit depth on HD209458b. Penning ionization dominates loss at the lowest altitudes ($r < 1.5 R_p$), with collisional de-excitation briefly dominating from ($1.5 R_p < r < 2 R_p$), and finally photoionization dominating past $2 R_p$. While our best-fit model uses updated reaction rates (outlined in Section 2.4), uncertainties in these rates can significantly impact the predicted He I 10830 Å transit depth. To test the sensitivity of our results to reaction rates, we run our best-fit model with reaction rates outlined in M. Lampón et al. (2020). As shown in the left panel of Figure 12, differences in these rates lead to nontrivial changes in the He I 10830 Å transit depth. The changes in the He I 10830 Å transit depth are mainly due to the new photoionization cross section and the ground-state helium recombination rate. M. Lampón et al. (2020) used a lower ground-state recombination rate than our model, allowing there to be more ionized helium and electrons in the atmosphere than in our model, leading to more

recombination to the metastable state. This comparison is limited to the rates listed by M. Lampón et al. (2020), and comparisons with other works are expected to differ.

The best-fit model successfully reproduces the observed He I 10830 Å absorption feature (F. J. Alonso-Floriano et al. 2019). Our results are also consistent with the upper limit of about -1% on $H\alpha$ absorption at line center set by the observations of A. G. Jensen et al. (2012), with our model predicting an $H\alpha$ excess absorption of -0.9% at line center. Our $H\alpha$ transit depth is also consistent with the broad upper limit from J. N. Winn et al. (2004). Using the same observations as J. N. Winn et al. (2004), however, N. Astudi-llo-Defru & P. Rojo (2013) inferred a transit depth of -0.123 ± 0.012% detection with a 1.125 Å bandpass, while the transit predicted by our model for this bandpass is higher at -0.4%. Based on the observations of N. Casasayas-Barris et al. (2021), we may infer a rough upper limit of -0.3% for a 0.5 Å bandpass, while our model predicts a transit depth of -0.73% for this bandpass. We note that the $H\alpha$ observations are not simultaneous with the He I 10830 Å transit observations. Future high-resolution observations could help further constrain the $H\alpha$ absorption and its variability to provide additional insight into the hydrogen ionization structure in escaping atmospheres.

## 3.5. Sensitivity of Transit Depths to Stellar Activity and Other Factors

To assess the key physical mechanisms driving the transit depths, we investigate the effects of stellar activity, diffusion processes, and tidal forcing. Stellar activity variations influence the high-energy flux incident on the atmosphere, altering the heating rate, mass-loss rate, and photoionization rates. Diffusion processes, including molecular and eddy diffusion, regulate the vertical transport of species and determine the degree of diffusive separation between hydrogen and helium. Tidal forcing modifies the gravitational potential and velocity structure, affecting the efficiency of mass loss. By systematically varying these processes, we aim to identify their individual and combined effects on atmospheric structure and the observed transit depths.

### 3.5.1. Stellar Activity Level

To investigate the effect of changing stellar activity levels on the upper atmosphere and the He I 10830 Å and $H\alpha$ transit depths, we run our best-fit model using stellar spectra corresponding to solar minimum, solar average, and solar
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Left panel: production mechanisms relevant for the formation of He I (2³S). Right panel: loss mechanisms relevant for the formation of He I (2³S).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Transit depth results from our best-fit model, comparing the modeled He I 10830 Å (left) and H I Balmer α (right) transit depths from our best-fit model (green solid lines) with the best fit to the observed transit depth for the He I 10830 Å line from F. J. Alonso-Floriano et al. (2019; blue line). He I 10830 Å transit depth results using previously used reaction rates shown in magenta (M. Lampón et al. 2020).</figcaption></figure>

maximum conditions (as described in Section 2.3). These spectra differ in both their overall XUV flux and overall spectral energy distribution and are all scaled to the orbital distance and stellar surface radius of HD209458b. We find that increasing stellar activity leads to higher mass-loss rates. Specifically, the mass-loss rate rises from $1.9 \times 10^{10}$ g s⁻¹ in the solar minimum case to $3.2 \times 10^{10}$ g s⁻¹ in the solar average case, and to $4.2 \times 10^{10}$ g s⁻¹ in the solar maximum case. This enhancement is primarily due to increases in both the outflow velocity and density (for example, taken here at 4 $R_p$, the approximate Roche Lobe radius). In the solar minimum case, the velocity and density at 4 $R_p$ are 1.1 km s⁻¹ and $8.78 \times 10^{-15}$ kg m⁻³, respectively, whereas in the solar average case, they increase to $1.5 \times 10^{-1}$ and $1.10 \times 10^{-14}$ kg m⁻³, respectively, and in the solar maximum case to $1.7 \times 10^{-1}$ and $1.25 \times 10^{-14}$ kg m⁻³, respectively. The pressure-averaged temperature between 1 µbar and the top of the model at $3 \times 10^{-4}$ nbar is 5720 K in the average case and 5906 K in the maximum case, compared to 5336 K in the minimum case. The pressure-averaged temperature between 7 nbar and the top

of the model at $3 \times 10^{-4}$ nbar is 8690 K in the average activity case and 9021 K in the maximum case, compared to 7970 K in the minimum case.

The degree of diffusive separation decreases with higher stellar activity, with the elemental He/H ratio at the top of the atmosphere increasing from 0.025 in the solar minimum case to 0.04 and 0.055 in the solar average and maximum cases, respectively. As predicted by previous studies (A. Oklopčić 2019; V. Krishnamurthy & N. B. Cowan 2024), the He I 10830 Å transit depth is highly sensitive to the stellar XUV flux. Interestingly, the He I 10830 Å transit depth is much more sensitive to the stellar XUV flux than the Hα transit depth, as seen in Figure 13. These results demonstrate that stellar activity levels significantly influence the He I 10830 Å transit depth, while the effect on the Hα transit depth is more limited. If HD209458b experiences an activity cycle similar to that of the Sun, variations in XUV flux over time could lead to periods with significantly larger transit depths than those observed. This suggests that future observations of HD209458b, taken at different phases of its host star's activity cycle, may reveal substantial variability in its
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Transmission spectrum around the He I triplet (left) and the Hα Balmer line (right). The transit depth produced using a maximum, average, and minimum solar spectrum is in pink, orange, and yellow, respectively.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Transmission spectrum around the He I triplet (left) and the Hα Balmer line (right). The model in purple, blue, and green have eddy diffusion coefficients of $10^7 \text{ m}^2 \text{ s}^{-1}$, $10^6 \text{ m}^2 \text{ s}^{-1}$, and $10^5 \text{ m}^2 \text{ s}^{-1}$, respectively. The model with no vertical diffusion is shown in pink.</figcaption></figure>

atmospheric escape signature. Additionally, the difference in sensitivity to stellar activity between the He I 10830 Å line and the Hα line presents an opportunity to better understand the interplay between stellar radiation and planetary escape mechanisms. By simultaneously observing both He I 10830 Å and Hα over multiple epochs, we can isolate the role of high-energy stellar radiation in driving escape.

### 3.5.2. Molecular and Eddy Diffusion

Studies have suggested that diffusive separation plays a key role in shaping He I 10830 Å transit depths (M. Lampón et al. 2020). To assess the importance of molecular and eddy diffusion, we investigate how they affect the structure of the upper atmosphere and the resulting He I 10830 Å and Hα transit depths. Running our model without any vertical diffusion predicts an He I 10830 Å transit depth nearly 5 times larger than observed, as seen in the left panel of Figure 14. Eddy diffusion ($K_{zz}$) is an essential component of upper atmosphere models, representing the aggregate effects of global circulation, turbulence, and breaking gravitational

waves. Without eddy diffusion, species mixing is driven only by bulk outflow and molecular diffusion, which may not accurately represent real atmospheric dynamics. To assess the impact of eddy diffusion, we run our best-fit model with three different values of $K_{zz}$—$10^5 \text{ m}^2 \text{ s}^{-1}$ (best-fit model), $10^6 \text{ m}^2 \text{ s}^{-1}$, and $10^7 \text{ m}^2 \text{ s}^{-1}$—and compare these to the case with no diffusion. The choice of these values is motivated by previous studies that have estimated a wide range of eddy diffusion coefficients for exoplanet atmospheres. J. I. Moses et al. (2011) suggested values between $10^8$ and $10^9 \text{ m}^2 \text{ s}^{-1}$ for hot Jupiter atmospheres, while V. Parmentier et al. (2013) proposed lower values around $10^6-10^7 \text{ m}^2 \text{ s}^{-1}$. For our best-fit model, we adopt a value of $10^5 \text{ m}^2 \text{ s}^{-1}$, consistent with the lower and middle atmosphere model used in this study and based on a recent assessment of eddy diffusion in exoplanet atmospheres (A. Arfaux & P. Lavvas 2023). Increasing the eddy diffusion coefficient by a factor of 10 reduces the mass-loss rate by a factor of 1.1 and decreases the degree of diffusive separation. The elemental He/H ratio at the upper boundary increases from 0.025 for $K_{zz} = 10^5 \text{ m}^2 \text{ s}^{-1}$ to 0.035 for $K_{zz} = 10^6 \text{ m}^2 \text{ s}^{-1}$ and 0.04 for $K_{zz} = 10^7 \text{ m}^2 \text{ s}^{-1}$.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Transmission spectrum around the He I triplet (left panel) and the Hα Balmer line (right panel). The model in green is our best-fit model, in blue is the zenith angle of 60° with tidal forcing, and in purple is the model run at the substellar line with tidal forcing.</figcaption></figure>

Consequently, higher eddy diffusion enhances the He I 10830 Å transit depth, as shown in the left panel of Figure 14. In contrast, the effect of diffusion is much weaker for the Hα transit depth (right panel of Figure 14). Since hydrogen constitutes the bulk of the atmosphere, the Hα transit depth remains largely unaffected even when diffusion is absent or with changing eddy diffusion. Comparing both lines—alongside independent constraints on stellar activity—can provide valuable insights into the role of molecular and eddy diffusion in shaping exoplanetary upper atmospheres. This analysis also points to a degeneracy between the eddy diffusion coefficient and the photoelectron heating efficiency: a lower $K_{zz}$ could, in principle, be compensated for by a higher heating efficiency, and vice versa. The adopted value of $10^5 \text{ m}^2\text{s}^{-1}$ represents a conservative choice that balances diffusive separation and observational constraints.

### 3.5.3. Tidal Forcing

Tidal forces (i.e., Roche lobe effects) play an important role in shaping the atmospheric dynamics of close-in exoplanets by modifying the pressure structure and accelerating mass loss. Some previous studies have incorporated tidal effects in models that are appropriate for the substellar line but not necessarily a good assumption for a global average (e.g., L. Xing et al. 2023). To investigate this, we implement tidal forcing in our upper atmosphere model by including the Roche potential in the momentum equation and running simulations at two different locations:

1. Substellar Case with Tide—At the substellar point, where tidal effects are strongest. Here, we divide the incident stellar flux by 4 to approximate global heat redistribution.

2. 60° Stellar Zenith Angle Case with Substellar Tide—At a stellar zenith angle of 60°, dividing stellar radiation by 2 (matching the best-fit model) to directly compare to the best-fit model and isolate tidal effects.

We find that introducing tidal forcing with a stellar zenith angle of 60° increases the mass-loss rate from $1.9 \times 10^{10}$ in the best-fit model to $3.7 \times 10^{10} \text{ g s}^{-1}$, then moving to the substellar point lowers this mass-loss rate to $2.5 \times 10^{10} \text{ g s}^{-1}$.

This occurs due to a combination of differences in outflow velocity, density structure, and the altitude of peak heating. The outflow velocity at $4 R_p$ increases significantly in the 60° stellar zenith angle case, reaching $6 \text{ km s}^{-1}$ compared to $1.1 \text{ km s}^{-1}$ in the best-fit model and $3.8 \text{ km s}^{-1}$ in the substellar case. However, the mass-loss rate depends on both velocity and density. In the 60° stellar zenith angle case, the expansion of the atmosphere leads to a lower density at $4 R_p$, while in the substellar case, the density remains slightly higher. The key factor driving these differences is the location of the heating peak. In the substellar case, the heating peak shifts to lower altitudes, decreasing the efficiency of the atmospheric expansion and dampening escape. Additionally, the sonic point is at $4.4 R_p$ in the substellar case and $4.2 R_p$ in the 60° stellar zenith angle case.

Interestingly, we show that tidal effects increase the He I 10830 Å transit depth, as shown in Figure 15. This arises from an increase in the metastable helium population at radii of $<2.5 R_p$ in the model with Roche lobe overflow that is driven by an increase in the $He^+$ population at these altitudes, producing a higher recombination rate to the metastable state. In the global average model with tidal forcing and a solar zenith angle of 60°, the He I 10830 Å transit depth is enhanced compared to the best-fit case without Roche lobe overflow, showing a slightly larger transit depth than the substellar model with Roche lobe overflow. This suggests that the increased outflow velocity and density changes from tidal effects enhance metastable helium absorption, but the exact magnitude depends on the incident stellar flux distribution. In contrast, the Hα transit depth substantially decreases in both tidal runs. The reduction in Hα transit depth in the tidal forcing models is primarily driven by a decrease in the Lyα mean intensity in the atmosphere, which reduces radiative excitation of hydrogen to the $n=2$ state. This occurs because tidal expansion lowers the overall density of the atmosphere, including H and electrons, decreasing collisional excitation rates and the scattering optical depth to Lyα photons. As a result, fewer Lyα photons are available to populate the H($n=2$) state, leading to weaker Hα absorption. The effect is most pronounced in the substellar case, where the density
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Left panel: temperature and bulk outflow velocity predicted by our atmosphere model for HD209458b including heavy elements (blue) compared to our best-fit model (green). Middle and right panels: transmission spectrum around the He I triplet (middle) and the Hα Balmer line (right). The transit depth produced by our best-fit model is in green compared to the model with heavy elements in blue.</figcaption></figure>

drops the most, causing the strongest suppression of neutral H density and thus the lowest Hα transit depth.

While tidal effects enhance the outflow velocity and alter the density structure, it is important to recognize that transit observations primarily probe the terminator region rather than the substellar line. Previous 1D models applying tidal effects along the substellar line do not necessarily reflect the geometry of transit observations. Our results demonstrate that models using substellar conditions produce transit depths that differ from those using globally averaged conditions. Since globally averaged conditions may better approximate the terminator region, while substellar conditions clearly do not, using substellar conditions in fits to transit observations can be problematic. In the end, 3D models will be required to address these differences (e.g., M. L. Khodachenko et al. 2021).

### 3.6. Metals

Heavy elements and metal species can play a key role in regulating the thermal structure and dynamics of the escaping atmosphere (C. Huang et al. 2023), yet they are often omitted from models of giant planet atmospheric escape. Observations of HD209458b suggest the presence of metals such as Fe, O, C, Mg, and Si in the upper atmosphere (A. Vidal-Madjar et al. 2004, 2013; J. L. Linsky et al. 2010; P. E. Cubillos et al. 2020). While the analysis and interpretation of these observations have been subject to debate and could be significantly affected by changes in stellar activity, it is reasonable to explore the impact of metals on atmospheric structure, atmospheric escape processes, and the He I 10830 Å and Hα transit depths. Therefore, we run our best-fit model incorporating Mg, Fe, Si, O, C, N, Na, K, S, Ca, and relevant ions. This model assumes that none of these species are removed from the upper atmosphere by condensation processes in the lower and middle atmosphere. The relevant details of this version of the model are given in full by C. Huang et al. (2023).

As can be seen in the temperature profiles in the left panel of Figure 16, the inclusion of metals results in a higher peak temperature in the upper atmosphere, as additional heating and changes in radiative cooling processes alter the thermal balance. The temperature peaks at approximately 9000 K near 1.7 R<sub>p</sub>, higher than in the best-fit model, which peaks at 8200 K at 1.7 R<sub>p</sub>. The corresponding mass-loss rate is 3 × 10<sup>10</sup> g s<sup>-1</sup>, slightly higher than in the best-fit model, which is 1.9 × 10<sup>10</sup> g s<sup>-1</sup>, but still within a physically plausible range for HD209458b. The pressure-averaged

temperature above 1 µbar to the top of the model at 3 × 10<sup>-4</sup> nbar is 5724 K, and above 7 nbar to the top of the model is 8614 K. The elemental He/H ratio decreases from 8% at the base of the thermosphere to approximately 5.5% at high altitudes, which is a less enhanced diffusive separation than in the best-fit model. As expected based on the diffusive separation of helium in the model, we also see diffusive separation between hydrogen and all of the heavier elements.

The most direct observational impact of metals in the atmosphere is reflected in the He I 10830 Å transit depths (middle panel Figure 16). The He I 10830 Å absorption depth in the metal-rich model is notably enhanced compared to the best-fit model. This increase is predominantly driven by increased ionization caused by the presence of metals that elevate the electron density and recombination to the metastable He I (2<sup>3</sup>S) state. In order to match the He I 10830 Å observations with this model, we now require a photoelectron heating efficiency of only 23%. On the other hand, the Hα transit depth (right panel Figure 16) remains essentially unchanged between models. This is because Hα line formation primarily depends on Lyα radiative excitation rate, and the neutral hydrogen density profile, both of which remain relatively stable even in the presence of additional metals.

### 4. Discussion

Our results highlight the importance of comprehensive atmosphere modeling in interpreting the escaping atmosphere of HD209458b. Assumptions about temperature structure, composition, and transport processes can significantly influence the inferred mass-loss rates and predicted transit depths. Parker wind models can provide a useful analytical framework for constraining the parameter space of the solution while self-consistent models help to narrow the solution and offer more detailed insights into upper atmospheric structure and processes. Including the lower and middle atmosphere increases the mass-loss rate by roughly 60%, while the He I 10830 Å and Hα transit depths remain largely unchanged. Although this suggests limited sensitivity to the properties of the lower and middle atmosphere of HD209458b, the structure and composition of the lower and middle atmosphere may play a more important role in other planets—particularly those with more extended lower and middle atmospheres or those with less extended lower and middle atmospheres but enhanced metallicities that feed higher abundances of heavy elements to
the upper atmosphere. Finally, we also explore how different excitation and de-excitation rates, the inclusion of metals, and stellar variability affect the predicted line profiles.

Parker wind models have been widely used to account for hydrodynamic atmospheric escape in forward model fits to observations of escaping atmospheres. This approach has led to mass-loss estimates for HD209458b of $8 \times 10^{10}$ g s$^{-1}$ (A. Oklopčić & C. M. Hirata 2018; L. A. Dos Santos et al. 2022), $0.42-1 \times 10^{11}$ g s$^{-1}$ (M. Lampón et al. 2020), and a high upper limit of $<2.4 \times 10^{11}$ g s$^{-1}$ (A. Masson et al. 2024). Our best-fit models require a photoelectron heat efficiency of 20%–40% and predict a mass-loss rate of $1.9-3 \times 10^{10}$ g s$^{-1}$, depending on whether metals are included or not in the upper atmosphere. We note that photoelectron heating efficiency here refers specifically to the fraction of photoelectron energy deposited as heat, rather than the net heating efficiency of the atmosphere. The overall heating efficiency, which accounts for all energy deposition processes relative to the incoming XUV flux and determines the energy-limited mass-loss rate, is 15% for our best-fit model. Our best-fit models also show diffusive separation between hydrogen and helium. Our results underscore the importance of self-consistently modeling the upper atmosphere at least above the 1 µbar level. This is because temperature gradients and dynamics shape the structure and composition of the upper atmosphere, influencing He I 10830 Å and Hα absorption. The processes that govern them must be treated carefully. While the lower and middle atmosphere help establish the continuum level for differential transit depth comparisons, the structure of the upper atmosphere remains the dominant factor in shaping the He I 10830 Å and Hα line cores on HD209458b. Our best-fit models are also consistent with the upper limit on Hα transit depth set by A. G. Jensen et al. (2012). However, as noted in Section 3, our model slightly exceeds other published constraints on Hα absorption in this system (N. Astudillo-Defru & P. Rojo 2013; N. Casasayas-Barris et al. 2021). This modest discrepancy may arise from geometric effects—combined with tidal forcing—which we show can strongly influence the Hα line, or from variations in stellar activity, to which the measured Hα transit depth is somewhat sensitive. That said, we refrain from drawing strong conclusions until longer-baseline observations of this line become available.

As shown in Sections 3.3 and 3.4, significant diffusive separation was observed in our models, and the degree of diffusive separation increased in our best-fit model. This demonstrates that reducing the heating efficiency enhances diffusive separation due to the decrease in mass-loss rate. We believe, however, that the degree of diffusive separation is surprising, given theoretical expectations based on the crossover mass equation (D. M. Hunten et al. 1987). The crossover mass equation provides a rough estimate of the limiting mass-loss rate, $\dot{M}_{\text{lim}}$, required for a species of mass $M_c$ (in units of $m_{\text{H}}$) to escape along with hydrogen via neutral-neutral collisions (T. T. Koskinen et al. 2014):

$$ \dot{M}_{\text{lim}} = 4\pi m_{\text{H}}^2 G M_p (M_c - 1) \frac{n D_c}{kT}. \quad (8) $$

Here, $D_c$ is the mutual diffusion coefficient for species c with H, given by T. R. Marrero & E. A. Mason (1972):

$$ \log(D_c p) = \log(A) + s \log(T) \quad (9) $$

where A and s are constants given by $14.2 \times 10^{-5} \text{ atm} \cdot \text{cm}^2$ and $1.721$, respectively, for He with H. At a temperature of 5336 K (the pressure-averaged temperature of our best-fit model), the mass-loss rate required to mix He with H is approximately $6.4 \times 10^9$ g s$^{-1}$. Given that our best-fit model estimates a mass-loss rate of $1.9 \times 10^{10}$ g s$^{-1}$, which is well above the limit, diffusive separation should be suppressed, yet we still observe significant separation between H and He in our model. This is despite the fact that the above equation is based on neutral-neutral collisions, while the upper atmosphere of HD209458b is strongly ionized, and this should further reduce the importance of diffusive separation.

The discrepancy between our results and the predictions of the crossover mass equation highlights the need to reconsider key assumptions in its derivation. The standard formulation assumes (1) no eddy diffusion, (2) a constant temperature, and (3) a volume mixing ratio for helium that remains constant with altitude. However, none of these assumptions hold in our self-consistent model. We find that the atmosphere exhibits a strong temperature gradient rather than a uniform temperature structure. Additionally, our simulations include nonzero eddy diffusion, which enhances vertical mixing and alters the competition between advection and molecular diffusion. Furthermore, there exists a region in our model where the molecular diffusion timescale is shorter than or comparable to both vertical advection and eddy diffusion, indicating that helium cannot maintain a constant mixing ratio throughout the atmosphere. These deviations from the idealized conditions of the crossover mass equation explain why we observe significant diffusive separation of helium even at mass-loss rates that exceed the theoretical limit for full mixing. Previous studies have also found that diffusive separation could be happening on HD209458b, but they did not include eddy diffusion (L. Xing et al. 2023; M. Schulik & J. Owen 2025). The choice of the eddy diffusion coefficient remains an important uncertainty that is degenerate with the mass-loss rate in controlling the degree of diffusive separation.

Excitation/de-excitation rates and photoionization cross sections obviously play an important role in determining metastable helium abundances and, consequently, He I 10830 Å transit depths. In this work, we present updates to the excitation/de-excitation rates and cross sections relevant to the formation and destruction of the metastable $2^3\text{S}$ state of He I 10830 Å in exoplanet atmospheres. The most notable revisions that we present are new photoionization cross sections and a new Penning ionization rate, described in Section 2.4. We also consider excitation/de-excitation mechanisms such as photoelectron collisions, which are important in Earth's atmosphere, and proton de-excitation. We provided an upper limit on photoelectron excitation for HD209458b and included proton de-excitation rates but did not find that they made a significant difference to the results (see Figure 11). For demonstration, we compare our best-fit model using the (de-)excitation rate scheme from M. Lampón et al. (2020) to our updated model that includes the new excitation scheme. The combined effect of these updates lowers the He I 10830 Å transit depth by a factor of ~2.5. The dominant factor in this reduction is the new high-resolution photoionization cross section, which includes the previously neglected EUV resonances, and the update to the ground-state recombination rate, which controls the population of ionized helium. These changes reduce the
overall metastable helium population and, consequently, the predicted transit depth. We note, however, that this difference is strongly model- and rate-dependent. In particular, in isothermal models where temperatures are fixed, the reduction in He I 10830 Å transit depth is expected to be smaller.

Metals are expected to influence the thermodynamics and composition of the escaping atmosphere, yet they are often omitted from models of atmospheric escape. Our model that includes heavy elements shows that the addition of metals has a relatively small effect on the temperature profile and mass-loss rate. The presence of metals that have lower ionization potentials than H and He, however, increases the electron density and boosts the recombination rate to the He I (2³S) state. This leads to deeper He I 10830 Å transit depths, requiring us to lower the photoelectron heating efficiency to 23% in order to maintain the match to the He I 10830 Å observations, while the Hα line shows little change from the best-fit model. In these models, we see significant diffusive separation of metals in the lower thermosphere, as expected based on our results for the He/H ratio in the upper atmosphere.

One-dimensional models are commonly used to interpret observations of exoplanet upper atmospheres, even though they do not account for the multidimensional nature of Roche lobe overflow or temperature and compositional differences between the day and night sides that affect global mass-loss rates. In addition to these limitations, other multidimensional phenomena—such as leading and trailing tails, asymmetries induced by stellar winds, or magnetic confinement—can introduce further complexities that are obviously not captured by 1D models (L. Wang & F. Dai 2021a, 2021b; M. L. Khodachenko et al. 2021). These effects are especially important in systems where extended or directional escape signatures are observed. Multidimensional simulations are essential for fully characterizing such cases. Many 1D models assume substellar illumination and Roche lobe overflow, which we show in Section 3.5.3, overestimate the mass-loss rate, and misrepresent the spatial distribution of helium and hydrogen at the transit terminator, impacting the modeled transit depths. Assumptions about the radiation field in 1D models are also important. We run our best-fit model using a dayside-averaged radiation field with global redistribution, achieved by assuming a solar zenith angle of 60° and dividing the incident stellar flux by a factor of 2. These are standard assumptions for globally averaged 1D simulations and may be closer to the conditions around the terminator, which is probed in transit. Even without the impact of Roche lobe overflow, the resulting transit depths differ from those in substellar geometry. Therefore, 1D forward models of the transit observations should always include an additional uncertainty in the derived parameters due to the lack of an accurate representation of multidimensional effects. In addition to these multidimensional considerations, interactions with the stellar wind could also influence mass-loss rates and the escape of different species. However, we do not expect stellar wind interactions to be a dominant factor in shaping the observed He I 10830 Å and Hα transit depths. The relatively low transit depths for both species are consistent with a transit obstacle that is significantly smaller than the size of the Roche lobe, reducing the likelihood of strong star-planet interactions (see also M. L. Khodachenko et al. 2021).

The presence of a planetary magnetic field, however, can influence atmospheric escape by modifying the structure of the outflow and dictating the regions where mass loss is permitted. The strength of this influence depends on the balance between thermal, magnetic, and dynamic pressures in the upper atmosphere. To better understand this interaction, we estimate the size of the magnetosphere to constrain the regions where the stellar wind and planetary outflow interact. While more complex models of the magnetosphere exist (M. L. Khodachenko et al. 2021), it is also useful to provide simple order-of-magnitude estimates that are more readily applicable to other planetary systems. One way to estimate the size of the magnetosphere is through the Chapman-Ferraro (C-F) boundary, which represents the distance at which the planetary magnetic pressure balances the incoming stellar wind pressure. The C-F standoff distance is given by:

$$R_{CF} \approx \zeta \left( \frac{m_B^2}{2\mu_0\rho_{sw}v_{sw}^2} \right)^{1/6} \quad (10)$$

where $\zeta$ is a scaling factor (~1.26), $m_B$ is the planet's magnetic moment, $\rho_{sw}$ and $v_{sw}$ are the density and velocity of the stellar wind for which we use $v_{sw} = 220 \text{ km s}^{-1}$, $\rho_{sw} = 4000 \text{ g cm}^{-3}$ (M. L. Khodachenko et al. 2021). We estimate the C-F stagnation distance to be $4.03 r/R_p$, with the polar standoff distance at $8.06 r/R_p$, marking the location where the planet's magnetosphere balances the incoming stellar wind pressure. However, the C-F model assumes a static balance between magnetic and dynamic pressures and does not account for the contribution of internal plasma pressure or the additional effects of planetary outflows. To address this limitation, we calculate the plasma beta ($\beta$), the ratio of thermal to magnetic pressure, which determines how much the escaping atmosphere influences the magnetosphere's structure. The plasma beta is given by:

$$\beta = \frac{2\mu_0 p}{B^2}, \quad (11)$$

where $\mu_0$ is the permeability of free space, p is thermal pressure, and B is the magnetic field strength, estimated by a dipole given by

$$B(r, \theta) = \frac{m_B}{r^3} (1 + 3 \cos^2 \theta)^{1/2}. \quad (12)$$

$m_B$ is the magnetic moment, which we estimate to be 10% of Jupiter's magnetic moment (K. G. Kislyakova et al. 2014), giving us a value of $2.83 \times 10^{19} \text{ T m}^3(\text{J. E. P. Connerney et al. 2017})$, r is radius, and $\theta$ is colatitude. The magnetic moment estimate used here is based on fitting Ly$\alpha$ transit observations with a model of the escaping atmosphere (M. L. Khodachenko et al. 2021). However, the interpretation of these observations is not unique (e.g., T. T. Koskinen et al. 2010), so the inferred magnetic moment should be considered uncertain. Based on calculations at the equator, $\beta$ significantly exceeds unity at all radii, as shown in Figure 17.

This implies that the thermal pressure of the plasma dominates over the magnetic pressure, leading to an extended and dynamic magnetosphere. This result is consistent with solar system giant planet magnetospheres that extend beyond the C-F standoff distance when significant internal plasma
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Plasma beta (β) as a function of radius (r/Rp) at the magnetic equator (pink) and pole (teal). The high values of β indicate that thermal pressure dominates over magnetic pressure at all altitudes, leading to an extended and dynamic magnetosphere influenced by planetary outflows.</figcaption></figure>

sources are present. The dynamic effects of the escaping planetary atmosphere must therefore be considered when estimating the true size of the magnetospheric boundary. To incorporate the effects of planetary outflow, we estimate the location of the magnetic “hydropause,” which is the location where the combined thermal, dynamic, and magnetic pressures of the stellar wind balance the combined thermal, dynamic, and magnetic pressures of the escaping planetary atmosphere. Unlike in solar system studies, where this boundary is often simplified to a static balance of magnetic pressures, exoplanetary systems require a full pressure balance equation that includes all relevant forces. This balance is described by:

$$ \rho_{sw} v_{sw}^2 + p_{sw} + \frac{B_{sw}^2}{2\mu_0} = \rho_{pl} v_{pl}^2 + p_{pl} + \frac{B_{pl}^2}{2\mu_0} \quad (13) $$

where the subscript *sw* refers to the stellar wind, and the subscript *pl* refers to the planetary magnetosphere. For the stellar wind magnetic field strength, we use $5 \times 10^{-9}$ T (M. L. Khodachenko et al. 2021). Balancing these terms in our model, we find that the hydropause is located at $4.7 r/R_p$, extending slightly beyond the C-F boundary. In addition, a stagnation distance of 4–5 $R_p$ suggests a polar cap boundary near a magnetic latitude of 60°–63°, indicating that approximately 27% of the planetary surface feeds open magnetic field lines along which escape is not suppressed. This suggests that the planetary outflow exerts enough pressure to expand beyond the traditional magnetospheric interaction region, reinforcing the notion of a dynamic magnetosphere shaped by both stellar wind interactions and atmospheric escape.

The upper atmosphere of HD209458b has also been probed by various FUV lines, including H Lyα, O I, C II, and Si III (A. Vidal-Madjar et al. 2003, 2004; L. Ben-Jaffel & S. Sona Hosseini 2010; J. L. Linsky et al. 2010), although the reality of the detection of Si III has been questioned (G. E. Ballester & L. Ben-Jaffel 2015). These observations provided some of the first detections of atmospheric escape on exoplanets and inspired many of the earlier modeling studies. Our best-fit model predicts a mass-loss rate of $1.9 \times 10^{10}$ g s⁻¹, which is within the mass-loss range of about $1–10 \times 10^{10}$ in other

modeling studies of HD209458b (e.g., A. Vidal-Madjar et al. 2003; R. V. Yelle 2004; F. Tian et al. 2005; A. García Muñoz 2007; T. T. Koskinen et al. 2010; J. L. Linsky et al. 2010), albeit at the lower end of this range. Our best-fit model predicts a pressure-averaged temperature of 7180 K when integrated from 1 µbar to 3 $R_p$, which is comparable to the temperatures required to match the FUV observations. The mass-loss rates and net heating efficiencies required to match the He I 10830 Å and Hα observations, however, are generally lower than the values fitted to the FUV transit observations (T. T. Koskinen et al. 2013b). Our best-fit model also predicts diffusive separation between He and H, which presents a possible challenge for matching the FUV observations of escaping heavy elements. Observations of O I and C II in the upper atmosphere suggest that these heavier species are not significantly depleted relative to hydrogen. If diffusive separation of He is occurring, we would expect these species to be less abundant at high altitudes, yet they are observed in the FUV with relatively large transit depths.

The He I 10830 Å transit depth is influenced not only by atmospheric properties but also by variations in the stellar XUV flux, as demonstrated in this (see Figure 13) and previous studies (A. Oklopčić & C. M. Hirata 2018; A. Oklopčić 2019; V. Krishnamurthy & N. B. Cowan 2024; G. Ballabio & J. E. Owen 2025; J. Sanz-Forcada et al. 2025). Our results suggest that repeated observations of HD209458b should reveal significant changes in the He I 10830 Å absorption unless the star is exceptionally quiet, with activity variations far lower than those of the Sun. This behavior contrasts with the relative stability of the Hα transit depth, which remains less sensitive to stellar activity cycles. This is because the H(n=2) population is predominantly excited by Lyα radiation, and the solar Lyα flux exhibits much smaller variability compared to the short-wavelength EUV flux over the solar cycle. We believe that the possible range of variability prevents direct comparisons between He I 10830 Å and the existing H Lyα or other FUV transit depths, as attempted in previous studies (e.g., A. Oklopčić & C. M. Hirata 2018). This is because the observations were taken at different times and, for example, H Lyα measurements are themselves sensitive to changes in stellar activity, in addition to the changes in the planet's atmosphere due to the changes in the incident XUV flux (e.g., T. T. Koskinen et al. 2010, 2013a). Instead, our model here is intended to be consistent with the He I 10830 Å transit depth and roughly consistent with the upper limit on the Hα Balmer line depth, both of which are relatively low, possibly reflecting a period of lower stellar activity. Generally, monitoring exoplanet He I 10830 Å and Hα transit depths across multiple epochs, together with host star activity indicators, provides valuable constraints on the host star's activity level and its temporal influence on atmospheric escape. Having said that, observers have struggled to detect the Hα line on HD209458b, and we have never observed He I 10830 Å transit depth consistent with moderate or higher stellar activity that we simulated based on the XUV flux variations of the Sun, a similarly quiet G-type star.

Finally, our results provide insights into why helium absorption detections do not always follow clear trends across exoplanet populations. While previous studies have suggested correlations between He I 10830 Å transit depths and planetary or stellar parameters such as irradiation levels (the likely extent of the planetary atmosphere) and mass-loss rates,
significant outliers exist (A. P. Allan & A. A. Vidotto 2025; G. Ballabio & J. E. Owen 2025; D. Linssen et al. 2025; J. Sanz-Forcada et al. 2025). The intermittency of helium detections can likely be attributed to a combination of planetary atmospheric conditions and stellar variability. One key factor is the impact of stellar activity variations. As our results also show, the He I 10830 Å transit depth is highly sensitive to changes in the stellar XUV flux. This means that repeated observations of the same planet could yield varying transit depths, depending on the phase of stellar activity at the time of observation. This variability complicates direct comparisons between different observations, as a single snapshot may not be representative of the true long-term behavior of the atmosphere. Unless HD209458 is significantly less active than the Sun, we would expect transit depths to vary across multiple epochs, reinforcing the need for continued monitoring. In addition to stellar activity, atmospheric properties such as diffusive separation and vertical mixing processes remain uncertain and vary between planets, as well as with stellar activity, even for the same planet. A study of these processes in different environments would provide broader context for understanding the role of vertical mixing, atmospheric composition, and mass-loss rates in shaping transit observations. However, improved models of stellar XUV fluxes and spectra are also necessary to make such a study more valuable, as uncertainties in the stellar input strongly affect He I 10830 Å predictions. Observationally, simultaneous He I 10830 Å and Hα transit measurements, along with possible other species, combined with stellar activity indicators would provide necessary insights into the interplay between planetary mass loss and host star variability. These measurements would help distinguish between intrinsic planetary properties and external stellar influences on transit depth variations. Additionally, repeated observations of HD209458b and other planets would be useful in building a time-averaged picture of mass loss, mitigating the effects of stellar activity variability and enabling more robust comparisons across exoplanet systems. Since we are aware of the time limitations of the observatories likely required for this purpose, we advocate for this goal to be achieved for a few well-selected systems to start with.

## 5. Conclusion

In this study, we develop a self-consistent model of HD209458b's full atmosphere, combining results from a multispecies hydrodynamic escape code with results from a photochemical model of the lower and middle atmosphere. By updating key excitation and de-excitation rates for metastable helium, and incorporating nonisothermal structure, stellar activity variations, diffusion, and metals, we identify key processes that shape He I 10830 Å and Hα transit depths. Our model reproduces the observed transit signals with a photoelectron heating efficiency of 20%–40% corresponding to a mass-loss rate of $1.9–3 \times 10^{10} \text{ g s}^{-1}$ and confirms that diffusive separation of helium occurs in the upper atmosphere. These results highlight the fact that interpreting helium absorption trends requires accounting for atmospheric structure and mixing, in addition to stellar variability. Simultaneous or near-simultaneous observations of He I, Hα, other species, and stellar activity indicators such as X-rays or Ca II H&K emission will be essential to fully constrain mass loss and

atmospheric escape processes across different exoplanet systems. While detailed modeling is required to fully interpret He I 10830 Å transit observations and simultaneous observations of multiple absorbers help to refine models of the upper atmosphere, observations of metastable helium alone still provide valuable insights. They are sensitive to the presence of an extended, escaping atmosphere that has consequences on the temperature and escape rate of the upper atmosphere. Even without other complementary constraints, observations of the He I 10830 Å line can help identify whether escape is likely and guide the plausible range of physical conditions, making them a potentially powerful diagnostic of the state of the upper atmosphere. Simple Parker wind models may offer a basic framework for estimating the parameter space for temperature, mass-loss rate, and the H/He ratio, but their assumptions are often too limiting. Self-consistent models are required to meaningfully narrow the solutions and rule out unphysical scenarios. Finally, we emphasize that planet-specific properties —such as gravity, irradiation, and composition—can lead to diverse atmospheric outcomes. Although the inclusion of metals had a limited effect on HD209458b, metal-rich atmospheres may exhibit significantly different thermal and compositional structures, including higher mean molecular weights and altered energy balance in the upper atmosphere.

## Acknowledgments

A.R.T. and T.T.K. acknowledge support by NASA/XRP 80NSSC23K0260. N.L. and L.A. acknowledge support by NSF grant No. PHY-1912507. C.H. is sponsored by Shanghai Pujiang Program (grant NO. 23PJ1414900). P.L. acknowledges support from the Programme National de Planétologie de CNRS/INSU under the project "Temperate Exoplanets."

*Software:* Astropy (Astropy Collaboration et al. 2018, 2022), Matplotlib (J. D. Hunter 2007), NumPy and SciPy (P. Virtanen et al. 2020), Pandas (W. McKinney 2010).

## Appendix

Here, we elaborate on how we obtained some of the revised de/excitation rates for metastable helium in this work. Figure 18 illustrates a schematic overview of all of the radiative and collisional transitions considered in our calculation. The radiative transition between the metastable state and the $2^3\text{P}$ state, responsible for the 10830 Å line, is not included in the modeled reaction network. This omission is due to the fact that this process does not significantly depopulate the metastable state; the atom in the $2^3\text{P}$ state simply decays back into the $2^3\text{S}$ state, preserving the triplet configuration. The resulting emission is isotropic and, therefore, does not significantly affect the observed transit depth.

For Penning ionization and proton de-excitation, we use cross sections to calculate a temperature-dependent rate to fit an equation already implemented into our KPP routine for simplicity using the Python package scipy. curve_fit and implement it into our atmosphere code. Our cross sections are shown and cited in Figure 19. The per volume per time reaction rate, $N_{12}$, is given by
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Atomic structure of the helium atom, indicating the radiative and collisional transitions included in our analysis. Adapted from A. Oklopčić & C. M. Hirata (2018).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Cross sections vs. energy for proton de-excitation (top left; L. Augustovičová et al. 2014), and penning ionization (top right; J. S. Cohen & N. F. Lane 1971; H. Morgner & A. Niehaus 1979). Reaction rates (cm³ s⁻¹) vs. temperature (K) for proton de-excitation (bottom left), and penning ionization (bottom right). The numerical solution for each rate is plotted in black, while the blue or green curves fit the rate.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Photoelectron cross section as a function of energy for the impact reaction between photoelectrons and helium.</figcaption></figure>

(e.g., T. I. Gombosi 1994)

$$N_{12} = 4\pi n_1 n_2 \left(\frac{m_{12}}{2\pi kT}\right)^{\frac{3}{2}} \int_{0}^{\infty} v^3 \sigma(v) \exp \left( -\left(\frac{m_{12}}{2kT}\right)^2 \right) dv \quad (A1)$$

where $n_1$ and $n_2$ are the number densities of the reactants, $v$ is velocity, $m_{12}$ is reduced mass, and $\sigma$ is the cross section. Using this method, we derive a numerical solution for the temperature-dependent reaction rate and implement it into KPP by curve-fitting the data to an equation already implemented into our KPP routine for simplicity using the Python package `scipy.optimize.curve_fit`. Figure 19 shows the numerical solutions compared with the curve fits we have implemented into our upper atmosphere code.

For the photoelectron calculation described in Section 2.4, for the cross section between a photoelectron and He I (2³S), we use discrete values given by F. J. de Heer & R. H. J. Jansen (1977) between the threshold energy at 19.82 eV and 80 eV, and an interpolated exponential function is used to describe the cross section beyond 80 eV, calculated using `scipy.interpolate.interp1d`. Once the interpolation reaches 0 at 145 eV, we assume the cross section is zero at all higher energies. Figure 20 shows this cross section versus energy.

## ORCID iDs

Anna Ruth Taylor https://orcid.org/0009-0003-5882-9663

Tommi T. Koskinen https://orcid.org/0000-0003-3071-8358

Luca Argenti https://orcid.org/0000-0003-1093-4379

Nicholas Lewis https://orcid.org/0000-0001-9580-3487

Chenliang Huang

(黄辰亮) https://orcid.org/0000-0001-9446-6853

Anthony Arfaux https://orcid.org/0000-0001-5097-4784

Panayotis Lavvas https://orcid.org/0000-0002-5360-3660

## References

Allan, A. P., & Vidotto, A. A. 2025, MNRAS, 539, 2144

Alonso-Floriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019, A&A, 629, A110

Arfaux, A., & Lavvas, P. 2022, MNRAS, 515, 4753

Arfaux, A., & Lavvas, P. 2023, MNRAS, 522, 2525

Argenti, L., & Lindroth, E. 2023, PhRvR, 5, 033047

Argenti, L., & Moccia, R. 2006, JPhB, 39, 2773

Argenti, L., & Moccia, R. 2008, JPhB, 41, 035002

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astudillo-Defru, N., & Rojo, P. 2013, A&A, 557, A56

Augustovičová, L., Kraemer, W. P., & Soldán, P. 2014, JQSRT, 148, 27

Bachau, H., Cormier, E., Decleva, P., Hansen, J. E., & Martín, F. 2001, RPPh, 64, 1815

Ballabio, G., & Owen, J. E. 2025, MNRAS, 537, 1305

Ballester, G. E., & Ben-Jaffel, L. 2015, ApJ, 804, 116

Ben-Jaffel, L., & Sona Hosseini, S. 2010, ApJ, 709, 1284

Benjamin, R. A., Skillman, E. D., & Smits, D. P. 1999, ApJ, 514, 307

Biassoni, F., Caldiroli, A., Gallo, E., et al. 2024, A&A, 682, A115

Bransden, B. H., & Joachain, C. J. 2003, Physics of Atoms and Molecules, 2nd edn. (Reading, MA: Addison-Wesley)

Bray, I., Burgess, A., Fursa, D. V., & Tully, J. A. 2000, A&AS, 146, 481

Brown, T. M. 2001, ApJ, 553, 1006

Casasayas-Barris, N., Palle, E., Stangret, M., et al. 2021, A&A, 647, A26

Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2020, A&A, 640, C6

Chadney, J. M., Koskinen, T. T., Galand, M., Unruh, Y. C., & Sanz-Forcada, J. 2017, A&A, 608, A75

Chamberlain, J. W. 1960, ApJ, 131, 47

Cohen, J. S., & Lane, N. F. 1971, CPL, 10, 623

Connerney, J. E. P., Adriani, A., Allegrini, F., et al. 2017, Sci, 356, 826

Cubillos, P. E., Fossati, L., Koskinen, T., et al. 2020, AJ, 159, 111

Czesla, S., Salz, M., Schneider, P. C., Mittag, M., & Schmitt, J. H. M. M. 2017, A&A, 607, A101

de Heer, F. J., & Jansen, R. H. J. 1977, JPhB, 10, 3741

Dos Santos, L. A., Vidotto, A. A., Vissapragada, S., et al. 2022, A&A, 659, A62

Drake, G. W. F. 1971, PhRvA, 3, 908

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109

García Muñoz, A. 2007, P&SS, 55, 1426

García Muñoz, A., & Schneider, P. C. 2019, ApJL, 884, L43

Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1

Gombosi, T. I. 1994, Gaskinetic Theory, Cambridge Atmospheric and Space Science Series (Cambridge: Cambridge Univ. Press)

Huang, C., Arras, P., Christie, D., & Li, Z.-Y. 2017, ApJ, 851, 150

Huang, C., Koskinen, T., Lavvas, P., & Fossati, L. 2023, ApJ, 951, 123

Hummer, D. G. 1962, MNRAS, 125, 461

Hunten, D. M., Pepin, R. O., & Walker, J. C. 1987, Icar, 69, 532

Hunter, J. D. 2007, CSE, 9, 90

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Indriolo, N., Hobbs, L. M., Hinkle, K. H., & McCall, B. J. 2009, ApJ, 703, 2131

Jaziri, A. Y., Pluriel, W., Bocchieri, A., et al. 2024, A&A, 684, A25

Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86

Judge, P. G., Solomon, S. C., & Ayres, T. R. 2003, ApJ, 593, 534

Kass, R. E., & Raftery, A. E. 1995, JASA, 90, 773

Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2021, MNRAS, 507, 3626

Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205

Kislyakova, K. G., Holmström, M., Lammer, H., Odert, P., & Khodachenko, M. L. 2014, Sci, 346, 981

Koskinen, T. T., Aylward, A. D., & Miller, S. 2007, Natur, 450, 845

Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013a, Icar, 226, 1678

Koskinen, T. T., Lavvas, P., Harris, M. J., & Yelle, R. V. 2014, RSPTA, 372, 20130089

Koskinen, T. T., Lavvas, P., Huang, C., et al. 2022, ApJ, 929, 52

Koskinen, T. T., Yelle, R. V., Harris, M. J., & Lavvas, P. 2013b, Icar, 226, 1695

Koskinen, T. T., Yelle, R. V., Lavvas, P., & Lewis, N. K. 2010, ApJ, 723, 116

Krishnamurthy, V., & Cowan, N. B. 2024, AJ, 168, 30

Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13

Lavvas, P., & Arfaux, A. 2021, MNRAS, 502, 5643

Lemaire, P., Emerich, C., Vial, J.-C., et al. 2005, AdSpR, 35, 384

Linsky, J. L., Yang, H., France, K., et al. 2010, ApJ, 717, 1291

Linssen, D., Oklopčić, A., & MacLeod, M. 2025, A&A, 698, A112

Marrero, T. R., & Mason, E. A. 1972, JPCRD, 1, 3

Masson, A., Vinatier, S., Bézard, B., et al. 2024, A&A, 688, A179

Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75
McKinney, W. 2010, in Proc. the 9th Python in Science Conf., ed. S. van der Walt & J. Millman (Austin, Texas: SciPy), 56

Morgner, H., & Niehaus, A. 1979, JPhB, 12, 1805

Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15

Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23

Norcross, D. W. 1971, JPhB, 4, 652

Oklopčić, A. 2019, ApJ, 881, 133

Oklopčić, A., & Hirata, C. M. 2018, ApJL, 855, L11

Owen, J. E., Shaikhislamov, I. F., Lammer, H., Fossati, L., & Khodachenko, M. L. 2020, SSRv, 216, 129

Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29

Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91

Roberge, W., & Dalgarno, A. 1982, ApJ, 255, 489

Sandu, A., & Sander, R. 2006, ACP, 6, 187

Sanz-Forcada, J., López-Puertas, M., Lampón, M., et al. 2025, A&A, 693, A285

Schulik, M., & Owen, J. 2025, MNRAS, 732

Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916

Shefov, N. N., Semenov, A. I., & Yurchenko, O. T. 2009, Ge&Ae, 49, 93

Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Natur, 529, 59

Smith, F. L., III, & Smith, C. 1972, JGR, 77, 3592

Southworth, J. 2010, MNRAS, 408, 1689

Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Natur, 557, 68

Stancil, P. C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41

Taylor, A., Koskinen, T., Argenti, L., & Lewis, N. 2024, Assessment of He I Triplet Absorption at 10830 Å in Escaping Atmospheres of Hot Gaseous Exoplanets (Leiden, Netherlands: Poster Presented at the Exoplanets V Conference)

Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, ApJ, 621, 1049

Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Vidal-Madjar, A., Désert, J. M., Lecavelier des Etangs, A., et al. 2004, ApJL, 604, L69

Vidal-Madjar, A., Huitson, C. M., Bourrier, V., et al. 2013, A&A, 560, A54

Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J. M., et al. 2003, Natur, 422, 143

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Voronov, G. S. 1997, ADNDT, 65, 1

Waldrop, L. S., Kerr, R. B., GonzáLez, S. A., et al. 2005, JGRA, 110, A08304

Wang, L., & Dai, F. 2021a, ApJ, 914, 98

Wang, L., & Dai, F. 2021b, ApJ, 914, 99

Winn, J. N., Suto, Y., Turner, E. L., et al. 2004, PASJ, 56, 655

Xing, L., Yan, D., & Guo, J. 2023, ApJ, 953, 166

Xue, Q., Bean, J. L., Zhang, M., et al. 2024, ApJL, 963, L5

Yan, M., Sadeghpour, H. R., & Dalgarno, A. 1998, ApJ, 496, 1044

Yelle, R. V. 2004, Icar, 170, 167

Zhang, M., Knutson, H. A., Wang, L., Dai, F., & Barragán, O. 2022, AJ, 163, 67