# A Hydrodynamic Study of the Escape of Metal Species and Excited Hydrogen from the Atmosphere of the Hot Jupiter WASP-121b

Chenliang Huang, Tommi Koskinen, Panayotis Lavvas, Luca Fossati

► To cite this version:

Chenliang Huang, Tommi Koskinen, Panayotis Lavvas, Luca Fossati. A Hydrodynamic Study of the Escape of Metal Species and Excited Hydrogen from the Atmosphere of the Hot Jupiter WASP-121b. The Astrophysical Journal, 2023, 951 (2), pp.123. 10.3847/1538-4357/accd5e. hal-04282389

HAL Id: hal-04282389

https://cnrs.hal.science/hal-04282389v1

Submitted on 13 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A hydrodynamic study of the escape of metal species and excited hydrogen from the atmosphere of the hot Jupiter WASP-121b

CHENLIANG HUANG (黄辰亮),¹,² TOMMI KOSKINEN,¹ PANAYOTIS LAVVAS,³ AND LUCA FOSSATI⁴

¹Lunar and Planetary Laboratory, University of Arizona, 1629 East University Boulevard, Tucson, AZ 85721

²Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China

³Groupe de Spectrométrie Moléculaire et Atmosphérique UMR CNRS 6089, Université Reims Champagne-Ardenne, 51687 Reims, France

⁴Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria

(Received; Revised; Accepted)

## ABSTRACT

In the near-UV and optical transmission spectrum of the hot Jupiter WASP-121b, recent observations have detected strong absorption features of Mg, Fe, Ca, and Hα, extending outside of the planet's Roche lobe. Studying these atomic signatures can directly trace the escaping atmosphere and constrain the energy balance of the upper atmosphere. To understand these features, we introduce a detailed forward model by expanding the capability of a one-dimensional model of the upper atmosphere and hydrodynamic escape to include important processes of atomic metal species. The hydrodynamic model is coupled to a Lyα Monte Carlo radiative transfer calculation to simulate the excited hydrogen population and associated heating/ionization effects. Using this model, we interpret the detected atomic features in the transmission spectrum of WASP-121b and explore the impact of metals and excited hydrogen on its upper atmosphere. We demonstrate the use of multiple absorption lines to impose stronger constraints on the properties of the upper atmosphere than the analysis of a single transmission feature can provide. In addition, the model shows that line broadening due to atmospheric outflow driven by the Roche lobe overflow is necessary to explain the observed line widths and highlights the importance of the high mass-loss rate caused by the Roche lobe overflow that requires careful consideration of the structure of the lower and middle atmosphere. We also show that metal species and excited state hydrogen can play an important role in the thermal and ionization balance of ultra-hot Jupiter thermospheres.

## 1. INTRODUCTION

Heating of the thermospheres of short-period exoplanets by stellar X-ray and extreme UV (XUV) radiation drives mass loss by atmospheric escape that shapes the population of known exoplanets and contributes to features such as the hot Neptune desert and radius gap between Mini-Neptunes with significant H/He envelopes and Super-Earths that have lost their envelopes (e.g., Owen & Wu 2017). An improved understanding of atmospheric escape and how it affects a broad range of planetary atmospheres is therefore necessary for a better understanding of the demographics of planetary systems and their evolution. In this respect, ultra-hot Jupiters such as WASP-121b in extremely close-in or-

bits around their host stars are particularly interesting because they offer a unique window to rapid hydrodynamic escape that can be observed with transmission spectroscopy. In particular, deep absorption features by atoms and atomic ions directly probe the escaping atmosphere (e.g., Ehrenreich et al. 2015; Ballester & Ben-Jaffel 2015; Spake et al. 2018; Yan & Henning 2018). Comparing these observations with models of atmospheric escape sheds light on the key physical processes and the properties of escaping atmospheres in extreme environments (e.g., Koskinen et al. 2013b; Khodachenko et al. 2019; García Muñoz & Schneider 2019; Shaikhislamov et al. 2020; Wang & Dai 2021a).

Due to its relatively low surface gravity and proximity to an F6V-type host star, the atmosphere of WASP-121b is likely to undergo rapid hydrodynamic escape enhanced by Roche lobe overflow, with a mass-loss rate that is among the highest for known exoplanets (Koski-
nen et al. 2022). The distance to the system (~270 pc) and the relatively high effective temperature of the host star do not favor the detection of the escaping atmosphere of WASP-121b at Lyman α or through the He 1083 nm feature that is often used to probe the upper atmospheres (Oklopčić 2019). However, the high spectral resolution and high signal-to-noise (S/N) detection of other atomic and ion lines in the optical and near-UV (NUV) spectrum of the planet means that WASP-121b offers some of the best observations of an escaping UHJ atmosphere. The interpretation of these observations, however, is complicated by the fact that they include several minor species and excited states that are not typically included in models of exoplanet atmospheric escape. New, more comprehensive models of exoplanet upper atmospheres are therefore required to unlock the full potential of these observations to characterize exoplanet atmospheres.

Using the Space Telescope Imaging Spectrograph (STIS) on board the Hubble Space Telescope (HST), Sing et al. (2019) detected strong absorption features of Mg II and Fe II in the NUV spectrum of WASP-121b, with transit depths corresponding to radii larger than the limb of the planetary Roche lobe (RL). Broad-band NUV observations of WASP-121b have also been obtained with the SWIFT telescope (Salz et al. 2019), where a tentative excess absorption signature was also seen. Using three archival transits acquired by the High Accuracy Radial velocity Planet Searcher (HARPS), atomic/ion absorption features of Fe I, Mg I, Ca I, Na I, Cr I, Ni I, V I, Fe II were detected using the cross-correlation method (Hoeijmakers et al. 2020; Ben-Yami et al. 2020; Bourrier et al. 2020), along with the absorption line profiles of Hα and Na doublets (Cabot et al. 2020). Observations of two transits with UV-Visual Echelle Spectrograph (UVES) at the Very Large Telescope (VLT) confirmed many of the previous detections and further detected the features of K I, Ca II, Sc II (Gibson et al. 2020; Merritt et al. 2021). Finally, Borsa et al. (2021) analyzed two transits observed with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations ESPRESSO/VLT and observed the absorption line profiles of Na doublets, Hα, Hβ, Ca II H&K, Mg I, K, Li, and Mn. The signal of the second transit was collected from all four Unit Telescopes (UTs) of the VLT simultaneously (4-UT mode), yielding a 16 m-equivalent telescope. Using the same ESPRESSO data, Azevedo Silva, T. et al. (2022) further detected Ba II, Co, and Sr II, and also pointed out that Ca II H&K absorption features are asymmetric with broader blue wings.

In this work, we present a new multi-species model of the upper atmosphere and escape for WASP-121b. Since escaping upper atmospheres do not exist independently of the lower and middle atmosphere, we use a photochemical model to simulate the composition of the atmosphere from 100 bar pressure to the lower boundary of the escape model at 1 µbar pressure. In order to compare the results with the observations, we use the combined model density and temperature profiles to simulate transit spectra. As shown by Koskinen et al. (2022), proper accounting for the extent of the whole atmosphere is critical for calculating the mass-loss rate from planets such as WASP-121b, and is also important for correctly explaining the observations. Here, we focus mostly on the NUV spectrum of the planet where the Fe II and Mg II lines clearly probe the escaping plasma. However, we also compare our model results with the absorption lines in the optical observed with the 4-UT mode of the VLT due to the unparalleled S/N of these data and the constraints that they provide on the atomic hydrogen distribution in the atmosphere.

Many of the recent theoretical studies on the escape of hot Jupiter atmospheres simulate the escape of hydrogen and protons, in some cases also including He, and focus on interpreting a single observed feature, such as the Hα line or the He 10830Å line (e.g. García Muñoz & Schneider 2019; Lampón et al. 2020; Wang & Dai 2021b). For example, Yan et al. (2021) constructed a model of hydrodynamic escape for WASP-121b that only included hydrogen and used the observed Hα line alone to constrain the properties of the upper atmosphere. In addition, they did not attempt to include the lower and middle atmosphere in their model. Our model extends the capabilities of similar previous models in two key ways. First, we include multiple minor species and compare our results with observations of several absorption lines simultaneously. Second, we account for the extent of the lower and middle atmosphere that affects the mass-loss rate and the interpretation of the observations. In summary, our focus on metal line absorption with model density profiles that span the whole atmosphere extends our understanding of WASP-121b's upper atmosphere.

The models used in this work are computationally expensive. By necessity at this point, therefore, they are one-dimensional, which can limit their ability to fully interpret the observations. Nevertheless, they provide a valuable tool for identifying the key physics necessary to model UHJ upper atmospheres. They also demonstrate both the capability and limitations of 1D models in explaining observations of planets like WASP-121b. All of this provides guidance for future developments of multi-dimensional models of exoplanet thermospheres
that include the relevant minor species with their excited states. This guidance comes in two forms. First, 1D models can include more of the relevant physics and thus the analysis of the results can be used to infer simplified, faster schemes for multi-dimensional models. Second, they point to features in the observations that require more complex models to interpret. In this way, our work provides motivation and focus for the development of 3D models of atmospheric escape that are particularly relevant to planets such as WASP-121b that undergo Roche lobe overflow.

This paper is organized as follows. First, we summarize the setup of the atmospheric model in Section 2 and discuss a reference atmosphere model for WASP-121b in Section 3. Next, we explain how we simulate the transmission spectrum that we compare with the observations in Section 4. Then we discuss our approach to simulating Roche lobe overflow in Section 5 and outline several factors that influence the results and present models that can in principle explain the observed features in Section 6. Finally, we discuss the key physical processes and properties of the best-fit model and its implications for the long-term evolution of WASP-121b's atmosphere in Section 7.

## 2. THE ATMOSPHERE MODEL

Separating the atmosphere of WASP-121b into two sections, we use a hydrodynamic model to simulate the upper atmosphere at $P < 1$ µbar, and a photochemical hydrostatic model to simulate the middle and lower atmosphere at $P > 1$ µbar. The transmission spectrum can then be calculated based on the combination of the results from the two models. In this section, we first discuss the planetary system parameters (Section 2.1) and the stellar spectra (Section 2.2) used in the hydrodynamic model and the photochemical model. An outline of the hydrodynamic atmospheric model is given in Section 2.3. Then, we explain the new physical mechanisms added to the hydrodynamic model, including the ionization balance of metal species (Section 2.4) and radiative cooling (Section 2.5). The photochemical hydrostatic model is described in Section 2.6. Finally, we describe the Lyα radiation transfer module that provides the excited state H population, as well as the associated heating effect, in Section 2.7.

### 2.1. System Parameters

In this work, we adopted the physical parameters of the planetary system given by Delrez et al. (2016) and summarized in Table 1. It is worth noting that the planetary radius given in the abstract of Delrez et al. (2016) is the planetary mean radius after asphericity correction

<table><caption>**Table 1.** Orbital and Physical Parameters of WASP-121b and WASP-121b (Delrez et al. 2016)</caption><thead><tr><th>Item</th><th>Value</th></tr></thead><tbody><tr><td>Stellar radius R<sub>*</sub></td><td>1.4572 ± 0.03 R<sub>☉</sub></td></tr><tr><td>Stellar mass M<sub>*</sub></td><td>1.3521<sup>+0.080</sup><sub>-0.079</sub> M<sub>☉</sub></td></tr><tr><td>Orbital period P</td><td>1.2749255 day</td></tr><tr><td>Planetary mass M<sub>p</sub></td><td>1.1824<sup>+0.064</sup><sub>-0.062</sub> M<sub>J</sub></td></tr><tr><td>Planetary radius R<sub>p</sub></td><td>1.766 ± 0.044 R<sub>J</sub></td></tr><tr><td>Distance D</td><td>263.2 ± 0.7 pc <sup>a</sup></td></tr><tr><td>Orbital separation a</td><td>0.02544<sup>+0.00049</sup><sub>-0.00050</sub> AU</td></tr></tbody></table>

(see Section 5.1) instead of the transit radius. In addition, Delrez et al. (2016) used the mean radius of Jupiter as the planetary radius unit, instead of the equatorial radius as suggested by the XXIXth IAU General Assembly (Mamajek et al. 2015). Using values provided in their Table 4, we list the system parameters of WASP-121b using the units suggested by the IAU in Table 1.

### 2.2. Stellar Spectrum

In order to constrain the XUV spectrum of the host star, we use the integrated XUV flux at the planet $F_{XUV} = 1.6 \times 10^6$ erg cm$^{-2}$ s$^{-1}$ from Salz et al. (2019), which is based on the measured stellar rotation rate. The spectral energy distribution (SED) at wavelengths shorter than 1700 Å is generated based on the solar SED. The SED is rescaled by the integrated XUV flux and the stellar radius considering that WASP-121 is larger than the Sun. The shape of the spectrum is consistent with observations of the star Procyon by the Extreme Ultraviolet Explorer (EUVE) but the flux is about 10 times higher than that of Procyon in the same wavelength range, which is expected given that Procyon is likely less active, while WASP-121 is a fast rotator and therefore likely more active (Fossati et al. 2018). The SED at longer wavelengths is generated by LLmodels (Shulyak et al. 2004), which has been specifically designed for simulating the spectra of F, A, and B type stars. The spectra between 1210 Å and 1220 Å have 0.5 Å wavelength bin width, while the rest parts of the spectra have 5 Å bin width.

The synthetic spectrum of WASP-121 is shown by the black line in Figure 1. The average solar spectrum from Koskinen et al. (2013a), a 6459 K blackbody, and the observed out-of-transit NUV spectrum (Sing et al. 2019) are shown for comparison. We note that the Lyman continuum flux of the spectrum is $2.69 \times 10^5$ erg cm$^{-2}$ s$^{-1}$. We assume full redistribution of energy in the planetary atmosphere model, approximated by using a zenith

<footnote><sup>a</sup> (Gaia Collaboration 2020)</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1.** A comparison between the synthetic WASP-121 stellar spectrum at the orbital distance of WASP-121b used in the simulation (black solid line), a 6459 K black body (black dashed line), and the solar spectrum (blue dashed line). The magenta line shows the spectrum with reduced high-energy intensity used in the best-fit model (see Section 6.1). The observed out-of-transit NUV stellar spectrum is shown by the red line (Sing et al. 2019).

angle of 60° for incident radiation, with a flux that is divided by a factor of two (Smith & Smith 1972).

The stellar Lyα emission line has a strong impact on the excited-state hydrogen population in the planet's atmosphere (see Section 2.7). The Lyα emission line of WASP-121 is not observable due to absorption by interstellar hydrogen along the long line of sight to the system from the Earth. Instead, we construct the Lyα line based on the Lyα profile of a similar F5V star—HD 106516—provided by Wood et al. (2005), assuming that they have the same surface Lyα intensity. We fit the line profile with a double-Voigt profile, shown as the dotted line in Figure 2. The Lyα flux inferred in this manner at the orbit of WASP-121b is $1.0 \times 10^5 \text{ erg cm}^{-2} \text{ s}^{-1}$. The Lyα emission line for moderately active Sun (Lemaire et al. 2005) is also shown for comparison.

To validate the stellar spectrum model described above, we compare the model integrated XUV fluxes in the wavelengths bands of 70 - 80 nm, 80 - 91.2 nm, and 91.2 - 117 nm with the fluxes based on the empirical XUV/Lyα flux equations suggested by Linsky et al. (2014). Based on the empirical flux ratio, we can estimate that the stellar fluxes in the bands of 70 - 80 nm, 80 - 91.2 nm, and 91.2 - 117 nm at the distance of the planet should be $1.2 \times 10^4$, $4.0 \times 10^4$, and $1.3 \times 10^4 \text{ erg cm}^{-2} \text{ s}^{-1}$ respectively. For comparison, the fluxes of the synthetic spectrum used in the calculation at these three bands are $1.4 \times 10^4$, $3.6 \times 10^4$, and $6.3 \times 10^4 \text{ erg cm}^{-2} \text{ s}^{-1}$, respectively.

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

**Figure 2.** Stellar Lyα profile at the planetary orbit. Black solid line: Reconstructed HD 106516 Lyα flux (Wood et al. 2005) scaled to the stellar radius of WASP-121. Black dotted line: A double-Voigt profile fit that is used in this calculation. Cyan dashed line: Moderate solar Lyα flux scaled by the radius of WASP-121.

The excess emission in the 91.2 - 117.0 nm band of the synthetic spectrum is dominated by emission in the Lyβ, C III 97.7 nm, and O IV 103.2 and 103.8 nm lines.

## 2.3. Hydrodynamic Atmospheric Escape Model

We use a grid-based time-dependent one-dimensional multi-species hydrodynamic atmosphere model CE-TIMB (Koskinen et al. 2013a,b; Koskinen et al. 2022) to simulate the coupled thermosphere and ionosphere of WASP-121b, including hydrodynamic escape. The model includes photoionization, thermal ionization, recombination, charge exchange, advection and thermal escape of ions and neutrals, diffusion, viscous drag, heating by photoionization, adiabatic heating and cooling, conduction, viscous dissipation, recombination cooling, radiative cooling, and free-free cooling. The general details of the model are described in Appendix B by Koskinen et al. (2022). Here, we include a description of photoionization, thermal ionization, and related heating (Section 2.4.1), charge exchange (Section 2.4.2), recombination and recombination cooling (Section 2.4.3), and radiative and free-free cooling (Section 2.5) that apply specifically to the simulations of WASP-121b.

The model has 580 stretched grid cells in the radial direction, with a spacing $\Delta r_0 = 10$ km at the lower boundary and a stretch factor of 1.014, spanning from the radius at 1 $\mu$bars up to 31.7 R$_J$ above the altitude of the bottom level. The location of the upper boundary is always outside the sonic point and the atmosphere undergoes hydrodynamic escape. Therefore, the model uses
the so-called outflow upper boundary conditions by extrapolating outflow velocity, species densities, and temperature from the model domain to the upper boundary point at a constant slope (Tian et al. 2005). Given that the upper boundary is above the sonic point, the solution at lower altitudes is not affected.

## 2.4. Ionization, Recombination and Charge Exchange

A total of 12 elements that are relatively abundant in a solar composition atmosphere or have spectral features that have been detected on WASP-121b are included in the model, including H, He, Mg, Fe, Si, O, C, N, S, Ca, Na, and K. Solar abundances (Asplund et al. 2009) are applied to these species by default. Second ionization states of Mg, Fe, Si, and Ca have been taken into account in the ionization/recombination balance. For the other elements, only neutral atoms and first ionization states are included. In the ionization level calculation, we consider photoionization, collisional ionization, and radiative and dielectronic recombination of all species listed, as well as charge exchange processes between these species whose rates are available. In the following, the reaction rates of these processes used in the model are described in detail.

### 2.4.1. Photoionization and Collisional Ionization

The photoionization cross-section tables of atoms and ions are built upon the sum of the fitting formulae for their outer shell (ground state) electron photoionization cross-sections (Verner et al. 1996), and their inner shell electron cross-sections (Verner & Yakovlev 1995). In addition, high-resolution tabulated photoionization cross-sections of H I, He I, C I, N I, O I, Na I, Mg I, Mg II, Si I, Si II, S I, Ca I, H I(2s), and H I(2p) from the Opacity Project (TOP, Cunto et al. 1993) are used to replace the fitting functions near the ionization threshold. In some cases, the ionization threshold listed by TOP is inconsistent with the value given by Verner & Yakovlev (1995). In these cases, we scale the photon energy from TOP to match the correct ionization threshold in Verner and Yakovlev (1995). The cross-sections in TOP are set to zero for photon energies that are lower than the photoionization threshold in TOP. For Fe I, we replace the ground state photoionization cross-section near the ionization threshold with the more accurate calculation by Zatsarinny et al. (2019). Tabulated cross-sections only cover the photon energy near the ionization threshold and often provide a simple power-law trend at the high-energy end. In these cases, cross-sections from fitting formulae are applied for energy higher than transitional energies that are chosen so that the two results can be smoothly connected. Although the cross-section of Ca II is also available in TOP, it is not adopted because we

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Photoionization cross-section of S I, Ca I, and Ca II. Solid lines show the cross-sections as given by Verner & Yakovlev (1995) and Verner et al. (1996). Dashed lines show the tabulated cross-sections given by TOP. The cross-section of Ca II from TOP is not adopted because it is described by a power law that does not agree with the formula from Verner & Yakovlev (1995). Averaged to the same wavelength bins as the input stellar spectrum, the dots show the cross-section value used by the model.</figcaption></figure>

doubt the reliability of the power-law trend given by TOP which is also significantly inconsistent with the fitting formula.

The high-resolution tabulated cross-sections are averaged to the same wavelength bins as the input stellar spectrum. As an example, the merged and averaged cross-sections of S I, Ca I, and Ca II used in our model are shown in Figure 3.

If photoelectrons produced by high-energy photons have sufficient energy, they can ionize and excite atoms and ions by collisions. The fraction of photoelectron energy that is converted to heat through collisions is the photoelectron heating efficiency. In this work, we assume a constant photoelectron heating efficiency of 93%, which is suggested for 50 eV photons at an electron mixing ratio of 0.1 (Cecchi-Pestellini et al. 2009). A relatively high heating efficiency overall is expected in hot, strongly ionized atmospheres in which the mean free path for Coulomb collisions is relatively short for all but the highest photoelectron energies.

Electron collisional ionization rates from Voronov (1997) are applied to all species that are included in the calculation. The energy consumed by the collisional ionization of hydrogen is accounted for in the calculation of the cooling rate in the atmosphere. The cooling rates resulting from collisional ionization of other, more minor atoms are assumed to be negligible.

### 2.4.2. Charge Exchange
A total of 65 charge exchange reactions between different atoms and ions are included in the model. The relevant rate expressions and references are listed in Table 4. Among these reactions, there are several processes for which reaction rates are only available in the exothermic direction. For instance, Kingdon & Ferland (1996) only lists the rate of the forward reaction Mg + H⁺, but not the rate of the reverse reaction Mg⁺ + H. The rates of the reverse reaction are likely to be small because they require a large kinetic energy input from the reactants, but for completeness we nevertheless estimate the reverse rates assuming microscopic balance.

According to microscopic balance, the ratio of the forward reaction to the reverse reaction rate is

$$K(T) = \frac{k_f}{k_r} = \exp\left(-\frac{\Delta G}{RT}\right), \quad (1)$$

where

$$\Delta G = G_{\text{products}} - G_{\text{reactants}} = \Delta(H - ST) \quad (2)$$

is the change in Gibbs free energy through the forward reaction, R is the gas constant, H is the enthalpy, and S is the entropy. The enthalpy and entropy of atoms and ions are calculated using fits to the results obtained by using empirical equations given by Gordon & McBride (1999). Here, the fitting function for $ln(K) = ln(a) + bln(T) - c/T$, where a, b, and c are the fitted parameters.

### 2.4.3. Recombination

We apply the case B recombination rate for H (e.g, Murray-Clay et al. 2009). Because the electron number density varies slowly in the upper atmosphere (see Figure 12), we ignore the weak dependence of the recombination rate coefficient on the electron density. Fitting the values given by Storey & Hummer (1995) at $n_e = 10^8$ cm⁻³, we have

$$\alpha_B = 1.00 \times 10^{-11} (T/300)^{-1.02} \text{ cm}^3 \text{ s}^{-1}. \quad (3)$$

For each hydrogen recombination, the average contribution to atmospheric cooling is (Draine 2011)

$$\langle E_{rr} \rangle = [0.684 - 0.0416 \ln(T_4)] k_B T, \quad (4)$$

where $T_4 = T/10^4$ K.

For He, we use the recombination rate from Yelle (2004). For Mg, Si, S, and Ca, we use the recombination rates from Shull & van Steenberg (1982). The fitting formulae from UMIST (Woodall et al. 2007) are used for O, C, and N. The recombination rates for Na and K are from Verner & Ferland (1996) and Landini & Fossi (1991) respectively.

We fit functions to the tabulated radiative and dielectric recombination rates of Fe I in Nahar et al. (1997) and rates of Fe II in Nahar (1997), and obtain

$$\begin{aligned} \alpha_{FeI} = 2.833 \times 10^{-8} T^{-1.5} \exp\left(-\frac{5.731 \times 10^4 \text{K}}{T}\right) & \\ \times \left(1 + 1.383 \times 10^4 \exp\left(-\frac{120.4 \text{K}}{T}\right)\right) & \\ + 1.248 \times 10^{-12} (T/10^4 \text{K})^{-0.485} \text{ cm}^3 \text{ s}^{-1}, & \end{aligned} \quad (5)$$

and

$$\begin{aligned} \alpha_{FeII} = 1.094 \times 10^{-5} T^{-1.5} \exp\left(-\frac{1.490 \times 10^4 \text{K}}{T}\right) & \\ \times \left(1 + 36.74 \exp\left(-\frac{1.153 \times 10^5 \text{K}}{T}\right)\right) & \\ + 1.728 \times 10^{-12} (T/10^4 \text{K})^{-0.618} \text{ cm}^3 \text{ s}^{-1} & \end{aligned} \quad (6)$$

for the recombination rates of Fe I and Fe II.

### 2.5. Radiative Cooling of Mg, Na, Ca, Fe, and H

When an atom or ion de-excites radiatively from an excited state that is collisionally excited by, say, an electron or hydrogen atom, the kinetic energy of the colliding particle is converted to radiation, and this effectively cools the atmosphere. To estimate the radiative cooling rates, we consider an approximate two-level system for each atom or ion, consisting of the two energy levels associated with the relevant transitions.

Under the two-level approximation, the lower and upper level populations satisfy,

$$\begin{aligned} n_l(n_e C_{lu}^{(e)} + n_H C_{lu}^{(H)} + B_{lu} \bar{J}_{lu}) & \\ = n_u(A_{ul} + n_e C_{ul}^{(e)} + n_H C_{ul}^{(H)}), & \end{aligned} \quad (7)$$

where $\bar{J}_{lu}$ is the line profile averaged mean intensity of the radiation field, $C_{ul}$ and $C_{lu}$ are collisional excitation and de-excitation rates, respectively, and the superscript indicates the type of collider particle, whether it is an electron or a hydrogen atom. Then, the radiative cooling rate of the transition can be written as

$$\begin{aligned} \Lambda(T, n_e, n_H) n_e n_l &= \Delta E(C_{lu}^{(e)} n_e n_l + C_{lu}^{(H)} n_H n_l \\ &\quad - C_{ul}^{(e)} n_e n_u - C_{ul}^{(H)} n_H n_u) \\ &= n_e n_l \Delta E \frac{C_{ul}^{(e)} + (n_H/n_e)C_{ul}^{(H)}}{A_{ul} + n_e C_{ul}^{(e)} + n_H C_{ul}^{(H)}} \\ &\quad \times \left( \frac{g_u}{g_l} \exp\left(-\frac{\Delta E}{k_B T}\right) A_{ul} - B_{lu} \bar{J}_{lu} \right), \end{aligned} \quad (8)$$

where $\Delta E$ is the energy difference between two levels.
Adopting the same method, Huang et al. (2017) estimated the radiative cooling rates of Mg, Na, K, O, C, Si, and S in the upper atmosphere of the hot Jupiter HD 189733b. The radiative excitation rates $B_{lu}J_{lu}$ are ignored because permitted transitions in the optical and near-UV bands, which dominate the radiative cooling, are always associated with strong absorption features in the stellar spectrum. The results show that Mg is the major coolant throughout the upper atmosphere, while cooling by Na can dominate at $P>0.5$ µbar where Na is less strongly ionized. Cooling by Na may have limited relevance, however, because on many Hot Jupiters, strong molecular coolants are likely to dominate in the middle atmosphere.

In this work, we estimate the Mg I, Mg II, and Na I cooling rates in the same way as Huang et al. (2017). The line cooling of Mg I λ2853 transition is included for Mg I, and Mg II λ2796 and λ2804 transitions are included for Mg II. For ion species, we also consider the free-free cooling rate using the formula (Osterbrock & Ferland 2006; Draine 2011)

$$ \Lambda_{ff} = 1.9095 \times 10^{-25} Z^2 T_4^{0.55} (\text{erg cm}^{-3} \text{s}^{-1}), \quad (9) $$

where Z is the charge number of the ion. To validate the calculated cooling rate, we compare our Mg I and Mg II cooling rates with the rates estimated by Gnat & Ferland (2012) using Cloudy (Ferland 1996, version 10.00), shown in Figure 4. Gnat & Ferland (2012)'s calculation, which also neglects radiative excitation, applies to a low-density gas ($n < 1 \text{ cm}^{-3}$) with a temperature above 10,000K. This shows that our simplified approach provides a good approximation of the most important metal line cooling rate in the upper atmosphere, which is due to Mg II.

Although radiative cooling by Ca and Fe is not considered by Huang et al. (2017) because of their low abundance in the upper atmosphere of HD 189733b, the detection of Ca and Fe absorption lines in the upper atmosphere of WASP-121b indicates that their radiative cooling should be evaluated in more detail. To estimate the cooling rate of Ca II, we use the electron collisional rates of λ3934 and λ3968 transitions in the CHIANTI database (Dere et al. 1997; Del Zanna et al. 2021) as well as the free-free cooling.

We also use the electron collisional rates from the CHIANTI database to estimate the Fe II radiative cooling rate. Of the 4339 Fe II transitions whose collisional rates are provided in the database, 1105 transitions that also have Einstein A coefficient or non-zero oscillator strength are included in the model.

Since the CHIANTI database does not include the collisional excitation rates of Fe I, we calculate the electron

<figure><img src="image_5.png" /><figcaption>Figure 4. Comparison of Mg I and Mg II cooling rates used in the hydrodynamic atmosphere model (blue and black solid line) with the tabulated value in Gnat & Ferland (2012).</figcaption></figure>

collisional excitation rates of 1025 permitted transitions with a minimum Einstein A of $10^6 \text{ s}^{-1}$ by using the van Regemorter formula (van Regemorter 1962; Jefferies 1968). The formula,

$$ C_{lu}^{(e)} = 2.16\alpha^{-1.68}e^{-\alpha}T^{-3/2}f \text{ cm}^3 \text{s}^{-1}, \quad (10) $$

where $\alpha = \Delta E/kT$, expresses the collisional excitation rate in terms of the oscillator strength $f$, which are acquired from the NIST atomic spectra database (Kramida et al. 2020). Although the van Regemorter formula does not apply to the collisional rates for dipole-forbidden transitions, their Einstein A are small so that the contribution of these transitions to cooling is generally negligible. In addition, the electron and neutral hydrogen collisional excitation rates of fine-structure and forbidden lines at 1.44, 1.36, 14.2, 22.3, 24.0, and 34.2 µm given by Hollenbach & McKee (1989) are included.

Besides the collisional excitation rates and spontaneous decay rates of each transition, the radiative cooling rate of a species also depends on the population of the lower levels that may of course differ from the ground state. For Fe I and Fe II, we assume that the energy levels that do not connect to any lower states with electric dipole transitions follow the Boltzmann distribution. For example, the lowest energy level that has an electric dipole transition with the ground state of Fe II is the 64th level ($3d^64p^6D^o$) at the energy of 4.8 eV. Thus, we assume that the 64th level of Fe II and all states above it are depleted, while states below follow the Boltzmann distribution. Similarly, we assume that the 54th level of Fe I ($3d^64s4p^5D^o$) at the energy of 3.2eV and all states above it are depleted. The same level population treat-
ment is applied to the transmission spectrum calculation described in Section 4.2.4.

Using a more sophisticated level population model in Cloudy, Verner et al. (1999) calculated the departure of each Fe II level from its LTE population in an iron dominant gas. In the radiation-dominated case, such as the situation discussed in this work, they show that the 64th level and levels above it are depleted, while the lowest 63 levels are only slightly under-populated compared to LTE (< 30% for the 45th level).

We likewise compare the cooling rates of Fe I and Fe II thus calculated with the results given by Gnat & Ferland (2012), shown in Figure 5 and 6. In the low number density environment ($n_e = n_{Fe} = 1 \text{ cm}^{-3}$, $n_H = 10^{-15} \text{ cm}^{-3}$) that they discuss, the spontaneous decay time is much shorter than the collisional excitation-time and thus the gas is in the coronal limit, meaning that only the ground state is significantly populated.

To compare with our Fe cooling rate at a broader temperature range, we repeat Gnat & Ferland (2012)'s calculation using Cloudy (version 17.01) and extend the temperature sampling down to 1500 K. Besides the low number densities used in Gnat & Ferland (2012), we also calculate the cooling rates at higher number density environments that are closer to the planetary upper atmosphere for comparison. Due to the restriction of the output information of Cloudy, to estimate the cooling rate, the number density of the studied ion species has to be much larger than the density of hydrogen. Therefore, we apply $n_e = n_H = 10^8 \text{ cm}^{-3}$ and $n_{FeI} = 10^{15} \text{ cm}^{-3}$ for Fe I, $n_e = n_H = 10^{10} \text{ cm}^{-3}$ and $n_{FeII} = 10^{15} \text{ cm}^{-3}$ for Fe II.

Previous hydrodynamic atmosphere models have shown that radiative cooling due to Lyα emission from the planetary atmosphere is important (e.g. Murray-Clay et al. 2009). In contrast to the cooling rate given by Black (1981) that has been used in many previous models, the Lyα cooling rate in this work uses the collisional excitation rates from the CHIANTI database, which are based on updated calculations by Anderson et al. (2002). Figures 7 compares the cooling rates obtained by the two methods, where Black (1981) gives a rate that is about 10 times higher compared to the new results here.

## 2.6. Lower and Middle Atmosphere

To reduce complexity, our hydrodynamic model does not include molecules such as CH$_4$, CO, NH$_3$, or H$_2$O, and associated photochemistry. Therefore, we set the bottom boundary of the hydrodynamic model at a pressure of 1 µbar, above which these molecules are completely dissociated. A hydrostatic photochem-

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

**Figure 5.** Comparison of Fe I cooling rate used in the hydrodynamic atmosphere model (black solid line) with the tabulated value in Gnat & Ferland (2012) (green line with circles), reproduced value using a newer version of Cloudy at the original $n_H$ and $n_e$ number densities (red solid line), and higher number densities (red dashed line). The contributions to cooling by dipole transitions (magenta) and fine-structure/forbidden lines (blue) are also shown separately.

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

**Figure 6.** Comparison of Fe II cooling rate used in the hydrodynamic atmosphere model (black solid line) with the tabulated value in Gnat & Ferland (2012) (green line with circles), reproduced value using a newer version of Cloudy at original number densities (red solid line), and higher number densities (red dashed line). The contributions to cooling by dipole transitions (blue) are shown separately. The free-free cooling is important at $T < 3000$ K.

ical/thermochemical model, which includes diffusion, condensation, more than 100 species, and more than 1500 reactions (Lavvas et al. 2014), is used to model the atmosphere at pressures of $10^{-6} < P < 100$ bar to pro-
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Radiative cooling rate due to the emission of Lyα photons. The rate provided by Black (1981) (blue) is compared to the cooling rate derived from the collisional excitation rate in the CHIANTI database (black).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Temperature profiles applied in the lower and middle atmosphere photochemical model. The lower atmospheric temperature profile presented in Mikal-Evans et al. (2019) as well as the dayside and nightside atmospheric temperature profiles presented in Mikal-Evans et al. (2022) are also shown for comparison.</figcaption></figure>

vide the radius and mixing ratios of the atoms at the lower boundary of the hydrodynamic escape model.

Shown by the solid blue line in Figure 8, the temperature profile of the photochemical model is assigned as follows. In the pressure range $10^{-4} < P < 10$ bar, we adopt the Markov chain Monte Carlo (MCMC) median temperature profile of the secondary eclipse observation retrieval analysis given by Mikal-Evans et al. (2019), shown as the magenta dashed line in Figure 8. On the high pressure end, this retrieved temperature is slightly below the Fe condensation temperature, making the model yield a very low abundance of Fe that would not be detectable in the upper atmosphere. To avoid Fe condensation, we increase the temperature at 10 bar by ~200 K, which is comparable to the $1-\sigma$ uncertainty in Mikal-Evans et al. (2019). A constant temperature is extrapolated from $P > 10$ bar up to $P = 100$ bar. Between $10^{-6} < P < 10^{-4}$ bar, the temperature is smoothly connected to the calculated upper atmosphere temperature profile at $1 \mu\text{bar}$.

The applied temperature profile is in a similar temperature range as the dayside temperature profile presented in Mikal-Evans et al. (2022) derived from observations in longer wavelengths, which are also shown in Figure 8.

We take the pressure at the visual broadband transit radius $R_p$ to be 4 mbar. The transit radius in the optical continuum shown in Section 4.1 supports this choice. The radius, temperature, and elemental abundances in the photochemical model at the bottom boundary of the hydrodynamic upper atmosphere model are consistent with each other. The lower and middle atmosphere

simulations also verify that all molecules are depleted at $1 \mu\text{bar}$ of this UHJ atmosphere as assumed.

## 2.7. Lyα Radiative Transfer, Excited State H and Associated Heating and Ionization

The strong $H\alpha$ absorption feature observed on WASP-121b calls for careful consideration of the distribution of excited state H. As have been shown in previous works (Huang et al. 2017; García Muñoz & Schneider 2019; Yan et al. 2021; Miroshnichenko et al. 2021; Yan et al. 2022), the excitation of the $H(n=2)$ state is dominated by radiative excitation and the population of $H(2p)$ is mostly determined by Lyα intensity through

$$n_{2p} = n_{1s} \frac{g_{2p} c^2}{g_{1s} 2h\nu^3} \bar{J}_{Ly\alpha}, \quad (11)$$

where $\bar{J}_{Ly\alpha} = \int_{0}^{\infty} J_{Ly\alpha}(\nu)\phi(\nu)d\nu$ is the Voigt profile averaged mean intensity.

In this work, we use the Lyα Monte Carlo radiative transfer calculation developed in Huang et al. (2017) to calculate the Lyα mean intensity inside a plane-parallel atmosphere. This calculation is computationally intensive and cannot realistically be coupled to our atmosphere model. Instead, we calculate the level populations separately based on the model output and iterate until the results are self-consistent.

Unlike the complete frequency redistribution function used by García Muñoz & Schneider (2019) to approximate the probability distribution of scattered photon frequency and direction for a given incident photon, we
apply the more accurate Hummer-IIB partial redistribution function (Hummer 1962). It accounts for the dipole angular dependence of resonant scattering when the initial and final states are the H(1s) state, and the intermediate H (2p) state has a finite lifetime and its fine-structure splitting is negligible. The Stokes matrix for Rayleigh scattering is used to track the polarization of propagating photons, which further determines the angular distribution of outgoing photons at scattering. Recoil is included in computing the new frequency of the photon after scattering.

Unpolarized Lyα photons are ejected into the scattering calculation through three physical channels. Stellar Lyα photons based on the line profile and flux given in Section 2.2 are incident vertically on the top boundary of the medium. The stellar Lyα intensity is divided by a factor of 2, which accounts for the uniform redistribution of stellar Lyα photons across the day-side of the atmosphere. Two internal Lyα photon generation channels are also included, electron-impact excitation and recombination cascade. We assume that every recombination of H (case B) produces one Lyα photon.

Besides exciting an H atom in the ground state, Lyα photons may photoionize Si, Mg, Fe, Ca, Na, and K, or photodissociate H₂, and thus end the scattering cycle. Although the absorption cross-section of H₂ is complex, Lyα spectrum within the atmosphere has a broad plateau near the line center within the velocity of ~±70 km s⁻¹ (Huang et al. 2017). Averaging the cross-section of H₂ between 1215.4 Å and 1216 Å at 2500 K, we have σH₂ = 4.0 × 10⁻¹⁹ cm² (Koskinen et al. 2021). Although Lyα can also photoionize excited-state H, we do not include this process in the calculation because of the following two reasons. First, the number density of excited-state H is much lower compared to the metal species listed above, thus it is not an effective sink for Lyα photon. In addition, because the Lyα intensity is negligible compared to the Balmer continuum, the effect of this process on the photoionization of excited hydrogen is also negligible.

Instead of resonant scattering, Lyα photons may also undergo complete redistribution if the excited H(2p) transits to H(2s) through collisional ℓ-mixing, which is followed by the reverse process back to H(2p). We fit the ℓ-mixing rate provided by Seaton (1955) for temperatures between 100 and 100,000 K using an analytical function and thus, the electron and proton collisional mixing rates are

$$C_{2s \to 2p}^{(e)} = 1.394 \times 10^{-3} T^{-0.4055} \exp(-13.35/T) + 8.723 \times 10^{-4} T^{-0.3803} \exp(-23.97/T) \text{ cm}^3 \text{ s}^{-1} \quad (12)$$

and

$$C_{2s \to 2p}^{(p)} = 3.706 \times 10^{-3} T^{-0.2509} \exp(-188.6/T) + 3.654 \times 10^{-4} T^{-0.07835} \exp(-1080/T) \text{ cm}^3 \text{ s}^{-1}, \quad (13)$$

respectively. The scattering Lyα photons are considered destroyed if H(2p) is photoionized by Balmer continuum photons or de-excited by electron collisions instead of radiative decay. If the collisional ℓ-mixing to H(2s) is followed by photoionization, collisional de-excitation, or two-photon decay, the Lyα photons are also considered destroyed.

To reduce the computational cost of the Monte Carlo simulation, we smooth the J̄L<sub>Lyα</sub> as a function of altitude within the atmosphere with B-spline. The procedure removes the high-frequency statistical noise while preserving the overall trend of the simulated function. For comparison, we note that Yan et al. (2022) undertook a 3-D Monte Carlo Lyα radiation transfer calculation with a spherically symmetric planetary atmosphere model. In their study, the Lyα radiation transfer calculation was conducted with the LaRT code (Seon & Kim 2020), which applies the Hummer-IIB partial redistribution function. The LaRT model does not account for the absorption of photons by metals or any other processes of H(2p) besides the emission of Lyα photons. Their results suggest, in agreement with our assumptions, that the distribution of excited state H is approximately spherically symmetric. Additionally, even when only considering the stellar Lyα source, where 3-D effects are most pronounced, the results from our plane-parallel simulation are consistent with the results of the 3-D radiative transfer model of Yan et al. (2022) in the substellar direction. With the same stellar Lyα flux, the peak intensity produced by our plane-parallel model differs by only ~10% from the peak intensity along the substellar direction presented in Fig. 13 of their work.

Based on the Hummer-IIB partial redistribution function, if the incident photon is a line-wing photon whose frequency is far from the line center, the scattered photon is most likely to be a line-wing photon whose frequency is within one Doppler width of the incident photon. In combination with the small scattering cross-section to line-wing photons, this allows the photons to travel large distances with a few scattering events (see Draine 2011, chap. 15.7). In comparison, the frequency of photons based on the complete redistribution function is always in the vicinity of the line core, making it much less likely that the photons escape. As a result, the Lyα intensity calculated using the complete frequency redistribution function for the resonant scattering process is typically higher than that obtained using the Hummer-IIB partial redistribution function.
Solving the rate equilibrium, we obtain the H(2s) number density

$$n_{2s} = \frac{n_{2p}(n_e C_{2p \to 2s}^{(e)} + n_p C_{2p \to 2s}^{(p)}) + n_H n_e C_{1s \to 2s} + n_p n_e \alpha_{2s}}{A_{2s} + \Gamma_{2s} + n_e C_{2s \to 2p}^{(e)} + n_p C_{2s \to 2p}^{(p)} + n_e C_{2s \to 1s}}, \quad (14)$$

where $C_{1s \to 2s}$ and $C_{2s \to 1s}$ are collisional excitation and de-excitation rates given by CHIANTI, $A_{2s} = 8.26 \text{ s}^{-1}$ is the two-photon decay rate, $\Gamma_{2s} = 576 \text{ s}^{-1}$ is the photoionization rate of H(2s), and $\alpha_{2s}$ is the recombination rate to H(2s) as the function of $n_e$ and T. Based on the tabulated rates in Storey & Sochi (2015) with the largest fitting parameter $\kappa$, which applies to a Maxwell-Boltzmann electron distribution, we use the following formula in the calculation:

$$\log_{10} \left( \frac{\alpha_{2s}}{\text{cm}^3 \text{s}^{-1}} \right) = \frac{0.697(\log_{10} n_e)^{2.675}}{T} - 0.269(\log_{10} T)^{1.428} - 11.136. \quad (15)$$

Having calculated $n_{2s}$, and $n_{2p}$, we update their contribution to hydrogen ionization as well as heating. Photoelectric heating from the H(n=2) state by the Balmer continuum as well as collisional de-excitation are included in the heating rates. As stated above, we iterate between the Monte Carlo model and the hydrodynamic model to reach self-consistency.

### 3. THE REFERENCE ATMOSPHERE MODEL

We first run the model under spherical symmetry with the observed system parameters of WASP-121b (see Table 1). Hereafter, we will refer to this model as **Case A** and it represents the globally averaged atmosphere of WASP-121b under the strong irradiation experienced by the planet at its current orbital distance, without Roche lobe effects. Figure 9 shows the temperature and radial velocity profile predicted by the model. The consistency between the temperature profile after two and three iterations demonstrates the efficiency of the radiative transfer and hydrodynamic simulations in achieving rapid convergence. The locations where H atoms are 50% ionized, where $\tau = 1$ for stellar 13.6 eV radiation, where the outflow speed equals the local adiabatic sound speed, and the size of the star are marked with vertical lines. In this case, the predicted mass-loss rate is $3.72 \times 10^{12} \text{ g s}^{-1}$ or $0.052 M_p/\text{Gyr}$.

Figure 10 shows the key heating and cooling rates. The thick black line shows the total photoionization heating rate while the contributions to this heating rate by ground state H, excited H, He, Si, and C are shown by the thin solid lines. Photoionization of ground state H and He is the primary source of heating. Adiabatic cooling is the dominant cooling mechanism at high altitudes

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Temperature and radial velocity simulated by the hydrodynamic model for Case A after three iterations. The blue dotted line that overlaps exactly with the black solid line shows the temperature after two iterations. The black dotted line shows the temperature of the case A' model that simulates an atmosphere with H and He only, excluding metals.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Rate of heating (solid line) and cooling (dashed or dotted line) process for case A. The thin solid lines in the figure represent the contributions of several major species to the total photoionization heating (thick black line).</figcaption></figure>

because of the relatively high mass-loss rate. Radiative cooling dominates at radii below about $2.2 R_p$, around the heating peak and below.

Figure 11 shows the major proton production rates and destruction rates. For a pair of charge exchange processes of H with other species (see Section 2.4.2), if the net effect is ionizing H, the net rates are shown with solid lines. Otherwise, the net rates are shown with dashed lines. In case A, the photoionization is well balanced by local recombination. This is in contrast to the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Rates of major proton production (solid) and destruction (dashed) for case A. As a component of H photoionization, photoionization of H(n = 2) is shown with the black dotted line.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Number density distributions of major species in both the upper atmosphere and lower atmosphere of case A. The region of the lower atmosphere calculated with the photochemical model is marked with green background.</figcaption></figure>

ionosphere of HD 209458b where advection from below is more efficient in replenishing H than recombination (Koskinen et al. 2013a). The difference is mainly due to the much higher ionization fraction of WASP-121b's atmosphere caused by significantly higher stellar XUV flux incident on WASP-121b.

Figure 12 shows the number density distributions of major species that contribute to absorption features calculated by the hydrodynamic model for the upper atmosphere and the hydrostatic atmosphere model for the lower atmosphere. Although the ionization balance is

calculated independently in the two models, the number densities of most atoms and ions from the two models connect smoothly. Metal species are mostly ionized throughout the upper atmosphere.

To analyze the impact of metal species on the model, we construct another H/He-only hydrodynamic upper atmospheric model based on the same lower/middle atmosphere setup and iteration procedure as in case A, which will be referred to as **case A'**. The temperature of case A' is shown as the black dotted line in Figure 9.

Metal species can impact the model through the following four mechanisms: (1) The absence of radiative cooling by metals, mostly by Mg II and Fe II, that dominate cooling between 1.15 and 1.4 $R_p$ (see Figure 10), should make the temperatures higher near the base of the model, (2) Si I, Mg I, Fe I, and Ca I can absorb Lyα photons, leading to a reduction in the Lyα mean intensity and ionization fraction, (3) heating by photoionization of Si I and C I can heat the atmosphere near the transition from the molecule-dominated to the atom-dominated layer, and (4) charge exchange between Fe I and $H^+$ may reduce the H ionization fraction. On WASP-121b, the case A' temperatures are warmer than the reference model temperatures near the base of the model and at higher altitudes. However, the presence of metals does not lead to an appreciable reduction in the Lyα intensity because the metals are mostly ionized throughout the model. For the same reason, heating by metal photoionization is not important either. Also, the charge exchange rate between Fe I and $H^+$ is lower than the H photoionization and recombination rates. We expect that mechanisms (2)-(4) can be more important on planets orbiting stars with lower XUV flux and mass-loss rates but they do not play a significant role on WASP-121b.

## 4. TRANSMISSION SPECTRUM

To calculate the transmission spectrum, we connect the upper atmosphere structure simulated by the hydrodynamic model and the lower atmosphere structure simulated with the photochemical model as shown in Figure 12. The outflow velocity of the lower and middle atmosphere photochemical model is set to conserve the overall mass-loss rate, given by the hydrodynamic model. Thus, we have radial distributions of temperature, outflow velocity, and number densities of the atmosphere at a total of 781 grid points spanning from 100 bar, which radius is referred to as $R_{core}$ in the following, up to ~ $10^{-6} \mu\text{bar}$ at ~19 $R_p$. We place this simulated planetary atmosphere at the center of the stellar disk to calculate absorption by the atmosphere and the transit depths.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Continuum of the simulated transmission spectrum. The horizontal dashed line represents the value of $R_p/R_*$ reported by Delrez et al. (2016). The broadband radius ratios derived from the same observational data by Evans et al. (2016) are also depicted.</figcaption></figure>

### 4.1. Transit Continuum

For the transit continuum that underlies the individual absorption lines, we include extinction by $H^-$ (John 1988) and Rayleigh scattering by H (Lee & Kim 2004), He (Fišák et al. 2017), and $H_2$ (Dalgarno & Williams 1965). Figure 13 shows a section of the simulated continuum profile, which is in agreement with the measured planet radius marked by the horizontal dashed line, or broad-band optical transit depths indicated by the markers within the error bar. This indicates that our choice of $P = 4$ mbar at $R_p$ is appropriate. Rayleigh scattering and extinction by $H^-$ are responsible for the continuum slope of the blue end and the red end, respectively. We note that our model does not include possible molecular absorbers, as our primary focus is on the interpretation of the upper atmosphere observations (see below).

### 4.2. Absorption Lines

For the line profile of each transition in the reference frame comoving with the bulk outflow of the atmosphere, we apply a Voigt profile with thermal broadening based on the local temperature and natural broadening based on each transition. We ignore collisional broadening since our focus is in the planetary upper atmosphere where it is negligible. We note that pressure broadening is obviously important for strong metal lines, such as those of Na and K, probing the middle and lower atmosphere. Since transit depths are not cumulative with altitude, however, we can still accurately represent the

profile of the line core even if our middle atmosphere transit depths are not fully accurate as long as the transit is driven by the upper atmosphere.

In calculating the transmission spectrum, we include the H Balmer, Mg I, Mg II, Ca I, Ca II, Fe I, Fe II, Na I, and K I lines. Here, it is simple to estimate the features of hydrogen-like metal atoms/ions because only a few electric dipole transitions from the ground state need to be considered. The included transitions for these species are Mg II λ2796, λ2804, Ca II λ3935, λ3970, Na I λ3303, λ3304, λ5892, λ5898, and K I λ4045, λ4048, λ7667, λ7701. We apply vacuum wavelengths throughout the simulation. The methods used to calculate the absorption features of other species in the model are discussed below.

#### 4.2.1. Balmer Lines

The Balmer line (Hα, Hβ, and Hγ) absorption profiles are calculated based on the H(2s) and H(2p) distributions (see black solid line in Figure 12). The procedure to estimate these distributions in the hydrodynamic atmosphere model is described in Section 2.7. Due to the high electron number density and low Lyα intensity in the lower atmosphere, the electron collisional processes may be more effective in exciting H atoms compared to radiative processes and thus H atom level population may follow Boltzmann distribution. Therefore, in the lower atmosphere where $P > 1\mu$bár, LTE H($n = 2$) number densities are applied if they are higher than the number density based on the Lyα intensity, obtained using equation 11. Because the lower atmosphere is relatively thin with a low number density of excited H, the contribution of the lower atmosphere to the absorption by the Balmer lines on WASP-121b is negligible.

#### 4.2.2. Mg I

In UHJ upper atmospheres, the electron collisional rates of almost all transitions are much smaller than radiative rates. The first excited state of Mg I at 2.71 eV (3s3p $^3$P$^o$) cannot spontaneously decay to the ground state through an electric dipole transition. We assume that the $^3$P$^o$ state and the ground state $^1$S can be described as a two-level system, using equation (7) with the radiative and collisional transitional rates between these two levels. An electron collisional de-excitation rate of $3 \times 10^{-8}$ cm$^3$ s$^{-1}$ is applied for the temperature range of the model (Osorio et al. 2015). The hydrogen atom collisional rate is always small for this transition (Barklem et al. 2012), and is thus ignored.

The second excited state of Mg I at 4.34 eV (3s3p $^1$P$^o$) can decay to the ground state through a permitted transition with a rate of $4.9 \times 10^8$ s$^{-1}$, which is significantly larger than any collisional rates. Therefore, this excited
state can only be populated through radiative excita-
tion. However, because it is in a deep absorption line in
the stellar spectrum, the number density of Mg I in this
excited state can be neglected. For a similar reason, we
assume all higher energy levels of Mg I are depleted in
our model.

From the NIST atomic spectra database, we retrieve
a total of 31 Mg I transitions with the lower energy level
being the ground state or ³Pº and the radiative decay
rate greater than 10⁶ s⁻¹.

4.2.3. Ca I

Another potentially abundant species that has an in-
tricate radiative network is Ca I. Based on the NIST
atomic spectra database, we supplement and update the
list of radiative transitions of Ca I according to the re-
cent calculation results of Yu & Derevianko (2018). Sim-
ilar to Mg I, we assume that the Ca I 3s4p ³Pº state at
1.88 eV and the ground state ¹S can be described as a
two-level system, with electron collisional de-excitation
rate of 8 × 10⁻⁸ cm³ s⁻¹ (Osorio et al. 2019), hydro-
gen collisional de-excitation rate of 3 × 10⁻¹⁰ cm³ s⁻¹
(Barklem 2016) and A = 2.74 × 10³ s⁻¹. In the lower
atmosphere, the hydrogen atom collisional de-excitation
rate of this transition can be larger than the spontaneous
radiative transition which is rarely seen in models. We
assume that other Ca I triplet states are depleted be-
cause they can spontaneously decay to a lower level.
Among the singlets, the first excited state 4s3d ¹D₂ is
a meta-stable state since it cannot decay to the ground
state by emitting a single photon. However, it can de-
cay to the ³Pº state through a forbidden transition with
A ~ 400 s⁻¹, faster than the production rate by colli-
sional excitation.

In this work, we include 55 Ca I transitions with lower levels no higher than $^3P^o$ and radiative decay rates greater than $5 \times 10^5$ s$^{-1}$. This is in contrast to large number of absorption lines in the Ca I cross-correlation template created based on an isothermal LTE atmosphere (e.g. Hoeijmakers et al. 2020). Although not observed, the model suggests that the Ca I $\lambda$4226 Å absorption feature should be stronger than K I $\lambda$7699 Å feature, and may be accompanied by a broad line wing.

4.2.4. Fe

As described in Section 2.5, we assume that the
$3d^64s4p \ ^5D^o$ state and higher energy levels of Fe I and
the $3d^64p \ ^6D^o$ state and higher of Fe II are depleted.
Combining this with the criteria of spontaneous rates
greater than $10^6$ s$^{-1}$, we include a total of 492 Fe I tran-
sitions and 829 Fe II transitions from the NIST atomic
spectra database in the calculation.

4.3. Line Broadening Due to Bulk Outflow Velocity
and Planet Rotation

Several broadening mechanisms, including natural broadening, collisional broadening, thermal broadening, and velocity broadening, determine the line width of spectral features. In the low-pressure environment of the planetary upper atmosphere, the effects of natural broadening and collisional broadening of metal lines are negligible compared to the observed FWHM of tens of km s⁻¹ of the spectral line probing the upper atmosphere (Borsa et al. 2021).

The width of many observed features cannot be ex-
plained by thermal broadening either. For example,
the required atmospheric temperature to explain the ob-
served ~ 80km s⁻¹ FWHM of Ca II features with ther-
mal broadening is unrealistically high at 2.8 × 10⁶ K.
In addition, because of the higher mass of Na compared
to H, the thermal broadening of the Na spectral lines
should produce significantly narrower lines than those
of hydrogen. This contradicts the observations showing
that the width of Na absorption line profiles is similar
to Hα (Borsa et al. 2021; Cabot et al. 2020). To ac-
count for this, we introduce the line-broadening effect
of radial outflow and planetary rotation. Both of them
can exceed 10 km s⁻¹ and are much higher than the
thermal velocity of the metal species. We assume that
the direction of the outflow (escape) velocity is radially
outward. The line-of-sight component of this outflow
velocity leads to line broadening, as illustrated by Fig-
ure 14.

The rotational velocity of a tidally locked planet can dominate line broadening in the atmosphere if the outflow velocity is small. In order to include rotational broadening in our calculation, we divide each quarter of the atmospheric annulus in the y-z plane into 20 circular sectors to calculate the atmospheric absorption, as illustrated by Figure 15 that shows a simplified model with only 3 sectors. Each sector has the same width in the y-direction — $\Delta y_i = r\Delta\sin(\theta_i)$ is a constant, where $\theta_i$ is the polar angle of the sector in the y-z plane. The line-of-sight rotational velocity is $v_{\text{rot},i} = \Omega y_i$, where $\Omega$ is the planet's rotational angular frequency. Thus, when rotational velocity dominates line broadening, the transit depth at a distance of $\Delta\nu/\nu_0$ away from the line center of an absorption feature corresponds to the absorption of stellar light by a planetary atmospheric slice with constant $y = \frac{\Delta\nu c}{\nu_0 \Omega}$.

This is the case for the spectral profiles with small
transit depths that probe the deeper atmosphere. They
typically show a dip at the line center (e.g. Na, K, and
Hβ in Figure 24) due to the geometry effect: an at-
mospheric slice with a small y has less area in the y-z
<figure><img src="image_15.png" /></figure>

**Figure 14.** Schematic diagram illustrates the velocity broadening mechanism of spectral lines and the line-of-sight optical depth integral. A light ray originating from the star on the right passes through the planet's atmosphere with an impact parameter *b* to reach the observer on the left. The width of the velocity broadening corresponds to the line-of-sight component of the atmospheric radial outflow velocity *v*<sub>r||</sub>. The optical depth τ<sub>λ</sub> along a light trajectory is the sum of the product of the length of the light propagation path within each radial shell of the model, illustrated with red dashed circles, and the opacity of that shell.

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

**Figure 15.** A cross-section view of a quarter of the planetary atmosphere in the y-z plane. The observer is outside the paper. This sketch illustrates the grids in radial and polar angular directions that are used to calculate the optical depth and transmission spectrum. The line-of-sight rotational velocity is Ωy. The three vertical yellow stripes represent three atmospheric slices, each with the same line-of-sight rotational velocity.

plane because the opaque core of the planet that does not contribute to the spectral features fills a large fraction. The slice that is approximately tangent to the planet's surface, shown as the middle vertical stripe in Figure 15, has the largest area. Thus, the split of the two peaks roughly corresponds to the rotational velocity at the planet's equator, which is approximately 7.3 km s⁻¹. A case study on the impact of these two broadening mechanisms on Ca II Kλ3934 Å line is presented in Section 7.3.

4.4. *Algorithm for calculating the transit spectrum*

In the wavelength domain, we divide the simulated spectral range of 2000 – 8000 Å into more than $2 \times 10^5$ wavelength points. Wavelength bands near the center of absorption features are divided into finer spacing, with a spectral resolution of up to $R = 2 \times 10^6$. First, a 2-D opacity array in the Lagrangian frame that does not include the bulk velocity is constructed for each shell and wavelength point according to the local number densities and temperature.

Because of the velocity shifts described in Section 4.3, the atmospheric optical depth $\tau_{\lambda}$ is a function of the impact parameter *b* of the light trajectory and its polar angle. To calculate the absorption in transit, we divide the cross-section annulus of the atmosphere into a 2-D grid in the radial and angular directions, as illustrated in Figure 15, and calculate $\tau_{\lambda}$ in each grid cell. In the angular direction, a quarter of the atmosphere is divided into 20 sectors as described in Section 4.3. In the radial direction, we uniformly sample 2000 grid points between $R_{\text{core}}$ and the stellar radius $R_{*}$. A layout with three radial grid points is shown in Figure 15. The opacity $\kappa_{\lambda}$ in the proper frame of each radial atmospheric shell, which is attached to the bulk outflow gas, is obtained by interpolating the 2-D opacity array to the wavelength points accounting for the blueshift/redshift caused by the bulk outflow and planet rotation. Multiplying the distance that the trajectory passes through each shell and integrating along the line of sight, we obtain the optical depth at a given wavelength, impact parameter, and polar angle $\tau_{\lambda}(b, \theta)$.

Trajectories with large impact parameters may have long segments in the innermost few shells they pass through, which leads to a large dispersion of the line-of-sight component of the atmospheric outflow velocity within the same shell. Because the single line-of-sight velocity values cannot adequately represent the whole segment, these segments are further divided along the trajectory to improve the characterization of spectral line shapes, illustrated with short yellow vertical bars in Figure 14.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Simulated transmission spectrum of case A. The red dashed line indicates the radius of the Roche Lobe at the terminator.</figcaption></figure>

The absorption of light is symmetric between the northern and southern hemispheres. Combining the optical depth of each trajectory with the cross-sectional area it represents, and summing over the impact parameters $b_j$ as well as the polar angles $\theta_i$ in one hemisphere, we obtain the square of the planetary apparent radius or stellar light attenuation fraction:

$$ \left[ \frac{R_p}{R_*}(\lambda) \right]^2 = \frac{1}{R_*^2} \left\{ R_{\text{core}}^2 + 2\Delta b \sum_{i,j} b_j \Delta \theta_i \{1 - \exp[-\tau_{\lambda}(b_j, \theta_i)]\} \right\}. \quad (16) $$

For example, the transmission spectrum calculated based on case A is shown in Figure 16. We note that the cores of the Mg II lines here correspond to an obstacle with a radius of about 3.08 $R_p$ ($R_p/R_* = 0.38$) viewed in transit while the cores of the other strong lines correspond to an obstacle of about 1.65 $R_p$ ($R_p/R_* = 0.21$).

## 4.5. Comparison with the NUV Observations

Since the spectral resolution of HST NUV observations is not high enough to resolve individual spectral features, we processed our simulated spectrum as follows to approximate the treatment of the data by Sing et al. (2019). First, we multiply the simulated stellar light attenuation fraction [$R_p/R_*(\lambda)]^2$ by the stellar spectrum to get the expected spectral energy distribution obscured by the planet during transit. Then, the obscured and stellar spectra are both convolved with a Gaussian smoothing function having an FWHM of 0.09Å, which mimics the HST spectral resolution of R=30,000. The

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Binned simulated NUV transmission spectrum of case A compared to HST observation (Sing et al. 2019). $\chi^2/N = 2.07$.</figcaption></figure>

ratio of these two spectra reproduces the expected transit radius to be observed. Finally, following the data reduction process, the convolved NUV spectra are divided into 187 4Å-wide bins with the same wavelength grid as the processed observed spectrum. The boxcar averaged transit depths for each bin are shown with black dots connected by a black solid line in Figure 17. The red horizontal dashed line marks the equipotential surface that coincides with the first Lagrangian point (L1) i.e., the boundary of the Roche lobe along the terminator that is lower than the distance to the L1 point indicated in Figure 8 of Sing et al. (2019).

Compared to the observations, the model has a reduced $\chi^2$ of 2.08 with 187 degrees of freedom. The spectrum is dominated by Fe II and strong Mg II lines. Although, for example, the size of the planet in the Mg II lines at full model resolution exceeds 35% of the stellar radius, the transit depths of all lines in the binned model spectrum are significantly lower than the observed transit depths. This is because the line widths based on the Case A model are too narrow to match the data.

## 4.6. Comparison with the Optical Observations

The absorption lines of the different species in the optical transmission spectrum provide an additional check on the model and much more robust constraints on the upper atmosphere than the NUV data alone. From Figures 9–12 in Borsa et al. (2021), we use `webplotdigitizer` (Rohatgi 2021) to extract the observed spectral features. We extract the profiles of absorption features as well as 5 additional data points in the continuum region on each side of the spectral lines. When comparing with the model, we normalize the observed features based on the
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Comparison of optical features observed by Borsa et al. (2021) with simulated profiles of case A. The horizontal dashed lines indicate $R_p$ and the radius of the Roche Lobe at the terminator. The black dashed lines in the panels of Hα and Hβ show the profile produced with the H/He only A' model.</figcaption></figure>

outermost 4 points on both sides to match them with the simulated profile. Borsa et al. (2021) noted that all detected absorption features are blueshifted by different degrees from the expected wavelengths of the lines. Due to the 1D nature of the model, our simulated spectra cannot reproduce these shifts. Therefore, in addition to shifting the observed wavelengths to vacuum values by applying air refractive index $n = 1.0002762$, we also shift the observed features according to the blueshift velocities reported in Table 4 in Borsa et al. (2021). This is equivalent to assuming that the blueshift arises in aggregate from horizontal winds in the atmosphere, as suggested by Borsa et al. (2021), that cannot be directly simulated by our approach. Because the spectral resolution is high enough to resolve spectral features, we do not convolve or bin the model spectral lines.

As a result, Figure 18 shows the comparison of several observed spectral features with the simulated transmission spectra obtained based on the Case A model. It is worth noting that although we show the Ca and Na doublet features in single panels, the two observed absorption peaks are normalized separately. The model provides a good fit to the Na and K lines, as they do not probe the escaping atmosphere where Na and K are mostly ionized. In line with the NUV spectrum, however, the model underestimates the strength and broad-

ening of the Ca II lines that probe the escaping part of the atmosphere.

In contrast to the lines of metal ions, the model overestimates absorption in the H Balmer lines, while not matching the broad wing of the Hα line. It is the result of the combination of many factors, including the high stellar Lyα intensity and the close proximity between the peak of the temperature and the peak of the neutral atomic hydrogen number density. This problem would be worse for models that do not include metals. The Hα and Hβ profiles produced with the metal-free A' model are depicted with black dashed lines in Figure 18. Due to the higher temperatures and the resulting more extended atmosphere, the Balmer lines exhibit deeper transit depths.

The comparison of our model spectrum with the NUV and optical observations indicates that a higher density and velocity of the escaping material are required to explain the observations. This implies a significantly higher mass-loss rate at the terminator of the planet than that predicted by our one-dimensional globally averaged model. Below, we argue that Roche lobe overflow (RLOF) can sufficiently enhance the model transit depths to match the observations and thus, the observations of WASP-121b provide a confirmation that the atmospheric escape from this planet is enhanced by
RLOF. We caution the reader that RLOF is an inherently multi-dimensional problem. In the absence of a multi-dimensional model that includes the required photochemistry and radiative transfer, we use results from our one-dimensional model to demonstrate the density and velocity enhancements that are required to match the data and discuss their relationship to what would be expected under RLOF.

## 5. TIDAL POTENTIAL

Because of the small surface gravity and the proximity to its host star, the stellar tide can significantly affect WASP-121b, as it does on WASP-12b (Dwivedi et al. 2019). As mentioned above, a larger atmospheric outflow velocity is necessary to match the observed spectral line profiles. A strong tidal force can reduce the average effective gravity, boost the mass-loss rate, and increase the outflow velocity. Therefore, in the following, we add the tidal (Roche) potential to the atmospheric escape model and simulation of the transmission spectrum to investigate its possible impact on the atmospheric properties and spectral line profiles.

### 5.1. Atmosphere Model with Tidal Potential

In the planetary atmosphere model, we define a Cartesian coordinate $x$ that points towards the host star, $y$ is the direction of the orbital motion, and $z$ points to the north pole. The tidal (Roche) potential of the two-body system in the co-rotating frame, centered at the tidally locked planet, is

$$ \phi(x, y, z) = -\frac{GM_p}{(x^2 + y^2 + z^2)^{1/2}} - \frac{GM_*}{[(a-x)^2 + y^2 + z^2]^{1/2}} - \frac{1}{2}\Omega^2\left[\frac{M_*a}{M_* + M_p} - x\right]^2 + y^2. \quad (17) $$

The planet radius $R_p$ is determined based on the transit depth assuming a circular planetary cross-section. However, because of the tidal potential, the planet $y$-axis size $R_{py}$ is larger than its equal-potential point in the $z$-axis $R_{pz}$. Considering that the surface of the planet is a nearly elliptical equipotential surface, we solve $R_{py}R_{pz} = R_p^2$ to obtain $R_{py}$ and $R_{pz}$, and then calculate the radius $R_{px}$ of the equipotential surface along the $x$-axis. The triaxial radii $R_{px}$, $R_{py}$, and $R_{pz}$ of the planet obtained after considering the tidal potential factor are listed in Table 2. These values agree with the structural parameters corrected for asphericity listed in Delrez et al. (2016), except for the different values of Jupiter radius units adopted, as mentioned in Section 2.1.

The two schematics in Figure 19 visualize, from two angles, the relative size of the planet with respect to the Roche lobe (RL) and the star when using the canonical parameters. The distance to the L1 point $R_{L1}$, and

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Geometry of WASP-121b system and its Roche lobe (RL). (a) The planet viewed from the pole. The host star is towards the left. The disk shows the planet and the blue curve shows the location of the RL. L1 indicates the location of the first Lagrangian point. (b) The planet as viewed by an observer in the middle of transit. The filled black disk shows the planet. The first blue elliptical shape indicates the RL. The outermost yellow circle shows the size of the host star. The two curves in between show the contours of two equipotential surfaces.</figcaption></figure>

<table><caption>**Table 2.** Derived Values and Best-fit Values for the Orbital and Physical Parameters of WASP-121 and WASP-121b</caption><thead><tr><th colspan="2">Derived quantities using the canonical value</th></tr></thead><tbody><tr><td>R<sub>px</sub></td><td>1.945 R<sub>J</sub></td></tr><tr><td>R<sub>py</sub></td><td>1.785 R<sub>J</sub> = 0.1259 R<sub>*</sub></td></tr><tr><td>R<sub>pz</sub></td><td>1.747 R<sub>J</sub> = 0.1232 R<sub>*</sub></td></tr><tr><td>R<sub>L1</sub></td><td>1.748 R<sub>px</sub></td></tr><tr><td>R<sub>RL</sub></td><td>2.233 R<sub>J</sub> = 0.1575 R<sub>*</sub></td></tr></tbody><thead><tr><th colspan="2">Parameters that provide a better match</th></tr></thead><tbody><tr><td>R<sub>*</sub></td><td>1.457 R<sub>☉</sub></td></tr><tr><td>M<sub>*</sub></td><td>1.352 M<sub>☉</sub></td></tr><tr><td>P</td><td>1.2749255 day</td></tr><tr><td>M<sub>p</sub></td><td>1.120 M<sub>J</sub></td></tr><tr><td>R<sub>p</sub></td><td>1.766 R<sub>J</sub></td></tr><tr><td>a</td><td>0.02545 AU</td></tr><tr><td>R<sub>px</sub></td><td>1.959 R<sub>J</sub></td></tr><tr><td>R<sub>py</sub></td><td>1.786 R<sub>J</sub></td></tr><tr><td>R<sub>pz</sub></td><td>1.746 R<sub>J</sub></td></tr><tr><td>R<sub>L1</sub></td><td>1.705 R<sub>px</sub></td></tr><tr><td>R<sub>RL</sub></td><td>2.194 R<sub>J</sub></td></tr></tbody></table>

the radius of the RL in the terminator plane $R_{RL}$ (see Section 5.4) are also listed in Table 2. The modified system parameters that can provide a better fit to the observations, as described in Section 6.1, are listed in the bottom half of Table 2.

In the following, we calculate escape through the L1 point along the substellar streamline by using gravita-
tional potential based on the Roche potential:

$$ \phi(x) = -\frac{GM_p}{x} - \frac{GM_*}{a-x} - \frac{1}{2}\Omega^2\left(\frac{M_*a}{M_*+M_p} - x\right)^2. \quad (18) $$

Because the photochemical model of the lower atmosphere does not include the tidal force, the pressure-radius relationship at pressures higher than 1 µbar is recalibrated to account for it. Using the temperature and mean molecular weight of each pressure level, we integrate the equation of hydrostatic equilibrium along the substellar direction both outward and inward starting from the radius $R_{px}$ where pressure is set to 4 mbar.

In contrast to the case A model, the stellar radiation is incident vertically on the atmosphere instead of at a 60° zenith angle. Accordingly, to represent globally averaged stellar radiation, the incident stellar flux is divided by a factor of 4 instead of a factor of 2. The atmosphere model that includes the tidal potential with canonical system parameters will be referred to as case B.

The temperature and radial velocity profile of case B is shown in Figure 20. At the radius equal to the stellar radius, the radial outflow velocity is over 7 times higher than in case A. The higher outflow velocity results in more effective adiabatic cooling and thus lower atmospheric temperature. Outside of the temperature peak, the hydrogen ionization fraction becomes lower than in case A due to the larger recombination coefficient caused by the lower temperatures. The radius at which $\tau = 1$ for stellar radiation with the energy of 13.6 eV, the ionization threshold of H, is at higher altitudes due to the higher neutral hydrogen column density. To connect smoothly to the lower temperature near the base of the hydrodynamic model when the tidal potential is included, the dotted line in Figure 21 is applied as the temperature for the photochemical lower atmosphere model.

## 5.2. Model inputs and results

We generate a suite of atmospheric models, the key information of which is listed in Table 3. The input parameters listed include whether the tidal potential is considered, the masses and radii of the planet and the star, orbital separations, temperature adjustment to the lower and middle atmosphere, planetary gravity at 1 µbar, multiplier of the stellar XUV flux, Lyα flux, and multiplier of the Fe abundance relative to the default solar value. More details on how these parameters are chosen are given in Section 6.1. The purpose of these models is to explore how different parameter choices and model inputs in general can change the predicted properties of the planetary atmosphere and the transit depths.

Correspondingly, we extract representative information from the result of each atmospheric model and the

<figure><img src="image_21.png" /><figcaption>Figure 20. Temperature and radial velocity along the substellar direction simulated by the hydrodynamic model for case B.</figcaption></figure>

simulated transmission spectrum and also list them in Table 3. The atmospheric properties included in the table are the simulated global mass-loss rate (see below), mass-loss rates according to energy-limited escape, without and with the RL correction factor (see Section 5.3), and the atmospheric pressure at the L1 point.

To represent the features in the simulated transmission spectrum, we select the line center transit depths of Ca II K (λ3935 Å), Hα, Hβ, and Na D$_2$ (λ5890 Å). As a comparison, the depths observed by Borsa et al. (2021) using 4-UT mode are 0.281±0.009, 0.186±0.003, 0.143±0.005, and 0.147±0.002, respectively. In addition, we list the Mg II λ2796 line center transit depth binned to a 4 Å passband, as shown in Figure 17, which can be compared to the observed value of 0.309±0.036 obtained by Sing et al. (2019) using the same passband at the same wavelength. The reduced $\chi^2$ that compares binned NUV simulated spectrum with observation is also listed.

We chose these spectral features because they can be directly compared with high-quality observations that trace species inside and escaping from the RL. Together, these features yield more information about the properties of the escaping atmosphere. The Na D$_2$ line is a good tracer of the lower and middle atmosphere below the RL boundary. In the escaping atmosphere, the most abundant ionization states of the metals, such as Mg$^{2+}$, Na$^{+}$, Ca$^{2+}$, or Fe$^{2+}$, do not have strong detectable spectral features at the observable wavelengths. The line center transit depth of the Ca II K line reflects the extent of the Ca$^{+}$ population, which is the second most abundant ionization level of the element. Combining it with the binned Mg II transit depths can further con-
strain the width of the spectral lines. On the other hand, Hα and Hβ features together reveal the extent and slope of the H(n = 2) population, which are sensitive to the atmospheric temperature and the Lyα intensity in the atmosphere.

### 5.3. Mass-Loss Rates

If the mass-loss rate in each direction is assumed to be the same as the rate along the simulated substellar streamline, then the mass-loss rate of the atmosphere is $\dot{M} = 4\pi r^2 v_r \rho$, where $v_r$ is the bulk radial outflow velocity. The mass-loss rate given by each model is listed in the first column of the bottom half of Table 3. This assumption is reasonable for planets with spherical gravitational potential, as long as the globally averaged stellar flux is used in the simulations. However, for systems with a non-negligible tidal potential, the mass-loss rate estimated using the substellar atmospheric streamline is an upper limit. For example, Guo (2013) compared results from 1D and 2D versions of his model to explore RL effects on a planet like HD209458b. The simulations indicate that a 1D model with substellar gravity overestimates the mass-loss rate by a factor of ~2 compared to a 2D model with uniform heating around the planet. According to the same study, significant day/night temperature differences in the upper atmosphere could reduce the mass-loss rate by a factor of 7 compared to a 1D model with substellar gravity.

We compare the mass-loss rates with the energy-limited mass-loss rates of the system, which are often assumed in the literature for atmospheric escape that is driven by stellar XUV flux. Considering only the gravitational potential energy of the planet itself and assuming the mass-loss efficiency $\epsilon = 1$, we have

$$ \dot{M}_{XUV} = \frac{R_p^3 L_{XUV}}{4a^2GM_p}, \qquad (19) $$

where $L_{XUV}$ is the stellar luminosity at energies higher than the ionization threshold of H at 13.6 eV. To account for the Roche potential, Erkaev et al. (2007) proposed a correction factor for the potential $K(\eta) = 1 - \frac{3}{2\eta} + \frac{1}{2\eta^3}$, where $\eta = R_{L1}/R_p$. The energy-limited $\dot{M}_{XUV}$ and $\dot{M}_{XUV}/K(\eta)$ are also listed in Table 3.

### 5.4. Transmission Spectrum with Tidal Potential

The calculation of the transmission spectrum is complicated by the fact that the geometric layout of the atmospheric cross-section is no longer spherically symmetric and we have to try and approximate it by using a 1-D model. Inspired by the equivalent sphere of the Roche lobe discussed in Eggleton (1983), we map

<figure><img src="image_22.png" /><figcaption>Figure 21. Temperature profiles applied in the lower and middle atmosphere photochemical model. The blue dashed line shows the temperature profile of case A, the same profile as the one shown in Figure 8, which does not include the tidal potential. The blue dotted line shows the temperature profile of case B, which is smoothly connected with an upper atmosphere temperature profile that considers the tidal potential. The blue solid line shows the temperature profile of the best-fit case D model, which uniformly increases the temperature of Case B by 350 K.</figcaption></figure>

our simulated 1-D atmosphere profile along the substellar direction to an equivalent spherical obstacle for the transmission spectrum calculation. As demonstrated in Figure 19, the outflow along the polar direction is suppressed by the tidal force compared to the equatorial direction. As a representative direction of the atmosphere in the terminator plane, we choose the atmosphere properties along the $\theta = 45^\circ$ direction in the $y-z$ plane to construct the equivalent spherically symmetric atmosphere. The radius of the RL in this representative direction is listed as $R_{RL}$ in Table 2.

To map the simulated atmosphere profile from the substellar direction to the $y-z$ plane within the RL, we assume the atmospheric temperature and number densities are uniform across the equipotential surface based on the assumption that the energy is efficiently redistributed around the planet. For regions beyond the L1 point, we conserve the number densities, temperature, and radial velocity, while scaling the radii in the atmosphere by the ratio of $R_{RL}$ to $R_{L1}$ at the RL boundary, which is approximately 2/3.

After this conversion, the reduced mass-loss rate per solid angle $\rho v_r r^2$ outside of the RL is naturally conserved along the radial direction, as in the original substellar solution. The mass-loss rate of our equivalent spherical atmosphere is suppressed by a factor of $\sim (2/3)^2 = 4/9$ compared to the 1D substellar atmosphere model. The
<table><caption>**Table 3.** Input parameters of each model and key result features</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Tidal potential</th><th><i>M</i><sub>p</sub></th><th><i>R</i><sub>p</sub></th><th><i>M</i><sub>*</sub></th><th><i>R</i><sub>*</sub></th><th><i>a</i></th><th>dT</th><th>log <i>g</i></th><th><i>a</i></th><th colspan="3">Multiplier</th></tr><tr><th>(<i>M</i><sub>J</sub>)</th><th>(<i>R</i><sub>J</sub>)</th><th>(<i>M</i><sub>☉</sub>)</th><th>(<i>R</i><sub>☉</sub>)</th><th>(AU)</th><th>(K)</th><th>@ 1μbar</th><th>XUV</th><th>Lyα</th><th>Fe</th></tr></thead><tbody><tr><td>Case A</td><td>N</td><td>1.1824</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>0</td><td>2.84</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Case B</td><td>Y</td><td>1.1824</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>0</td><td>2.70</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Case C</td><td>Y</td><td>1.1204</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>350</td><td>2.62</td><td>1</td><td>1</td><td>4</td></tr><tr><td><b>Case D</b><sup>e</sup></td><td>Y</td><td>1.1204</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>350</td><td>2.62</td><td>0.74</td><td>0.35</td><td>4</td></tr><tr><td>Case E</td><td>Y</td><td>1.1204</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>350</td><td>2.62</td><td>1</td><td>0.25</td><td>5</td></tr><tr><td>Case F</td><td>Y</td><td>1.1204</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>350</td><td>2.62</td><td>0.47</td><td>0.6</td><td>3.5</td></tr><tr><td>Case G</td><td>N</td><td>1.1204</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>350</td><td>2.80</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Case H</td><td>N</td><td>0.9964</td><td>1.8911</td><td>1.6603</td><td>1.5604</td><td>0.02725</td><td>1000</td><td>2.63</td><td>1</td><td>1</td><td>4</td></tr><tr><td>Case I</td><td>Y</td><td>1.152</td><td>1.766</td><td>1.3521</td><td>1.4572</td><td>0.02545</td><td>550</td><td>2.63</td><td>0.74</td><td>0.35</td><td>4</td></tr><tr><td>Case J</td><td>Y</td><td>1.1204</td><td>1.785</td><td>1.396</td><td>1.473</td><td>0.02572</td><td>120</td><td>2.61</td><td>0.21</td><td>0.7</td><td>3.5</td></tr><tr><td>Case K</td><td>Y</td><td>1.1204</td><td>1.785</td><td>1.396</td><td>1.473</td><td>0.02572</td><td>120</td><td>2.61</td><td>1</td><td>0.15</td><td>4</td></tr></tbody></table>

<table><thead><tr><th><i>M</i>&#775;</th><th><i>M&#775;</i><sub>XUV</sub></th><th><i>M&#775;</i><sub>XUV</sub><br><i>K</i>(η)</th><th>P @ L1</th><th>Ca II K <sup>b</sup></th><th>Hα <sup>b</sup></th><th>Hβ <sup>b</sup></th><th>Na D<sub>2</sub> <sup>b</sup></th><th>Mg II <sup>c</sup></th><th><i>χ</i><sup>2</sup>/<i>N</i></th></tr><tr><th colspan="2">(<i>M</i><sub>p</sub>/Gyr)<sup>d</sup></th><th colspan="2">(nbar)</th><th colspan="5">(<i>R</i><sub>p</sub>/<i>R</i><sub>*</sub>)</th><th></th></tr></thead><tbody><tr><td>0.052</td><td>0.16</td><td>0.55</td><td>1.07 × 10<sup>-3</sup></td><td>0.199</td><td>0.201</td><td>0.174</td><td>0.152</td><td>0.182</td><td>2.08</td></tr><tr><td>0.32</td><td>0.16</td><td>0.68</td><td>2.67</td><td>0.200</td><td>0.166</td><td>0.133</td><td>0.141</td><td>0.202</td><td>2.03</td></tr><tr><td>1.20</td><td>0.18</td><td>0.80</td><td>2.52</td><td>0.270</td><td>0.221</td><td>0.145</td><td>0.147</td><td>0.303</td><td>1.18</td></tr><tr><td><b>1.03</b></td><td>0.13</td><td>0.59</td><td>2.02</td><td>0.278</td><td>0.185</td><td>0.135</td><td>0.147</td><td>0.302</td><td>1.20</td></tr><tr><td>1.21</td><td>0.18</td><td>0.80</td><td>2.51</td><td>0.269</td><td>0.191</td><td>0.137</td><td>0.147</td><td>0.303</td><td>1.21</td></tr><tr><td>0.80</td><td>0.084</td><td>0.38</td><td>1.47</td><td>0.268</td><td>0.180</td><td>0.133</td><td>0.143</td><td>0.292</td><td>1.24</td></tr><tr><td>0.073</td><td>0.18</td><td>0.63</td><td>1.31 × 10<sup>-3</sup></td><td>0.217</td><td>0.209</td><td>0.181</td><td>0.155</td><td>0.192</td><td>1.90</td></tr><tr><td>0.188</td><td>0.27</td><td>1.25</td><td>2.69 × 10<sup>-3</sup></td><td>0.256</td><td>0.241</td><td>0.206</td><td>0.166</td><td>0.229</td><td>1.43</td></tr><tr><td>1.08</td><td>0.124</td><td>0.54</td><td>2.07</td><td>0.280</td><td>0.185</td><td>0.135</td><td>0.148</td><td>0.307</td><td>1.22</td></tr><tr><td>0.55</td><td>0.038</td><td>0.179</td><td>0.85</td><td>0.296</td><td>0.183</td><td>0.135</td><td>0.143</td><td>0.318</td><td>1.33</td></tr><tr><td>1.3</td><td>0.183</td><td>0.86</td><td>2.51</td><td>0.265</td><td>0.183</td><td>0.135</td><td>0.144</td><td>0.308</td><td>1.2</td></tr></tbody></table>

<sup>a</sup><i>g</i> = GM<sub>p</sub>/<i>r</i> in cgs unit, where <i>r</i> is the radius where P=1 μbar in the substellar direction.

<sup>b</sup>The line center transit depths of Ca II K (λ3935 Å), Hα, Hβ, and Na D<sub>2</sub> (λ5890 Å) obtained by Borsa et al. (2021) using the 4-UT mode are 0.281±0.009, 0.186±0.003, 0.143±0.005, and 0.147±0.002 respectively.

<sup>c</sup>Mg II λ2796 transit depth binned with 4 Å passband as shown in Figure 17. The transit depth obtained by Sing et al. (2019) using the same passband is 0.309±0.036.

<sup>d</sup><i>M</i><sub>p</sub>/Gyr = 7.1 × 10<sup>13</sup> g s<sup>-1</sup>.

<sup>e</sup>Preferred model.

similarity of this ratio to the factor of 2 difference in mass-loss rates based on 1D and 2D models found by Guo (2013) provides some justification for the conversion that we apply.

Within the RL, however, the mass-loss rate is not conserved along the radial direction if the outflow velocity is assumed to be uniform across the equipotential surface. Therefore, we adjust the outflow velocity below the RL boundary to conserve mass based on the converted number densities and radii. Since the outflow velocity

within the RL makes only a minimal contribution to the spectral line broadening and is smaller than the effect of rotational velocity, this treatment of the radial velocity does not have a noticeable impact on the line profiles in the calculated transmission spectrum.

The use of a 1-D model to simulate the atmospheric structure within strong tidal potential has obvious limitations. The deviation of the escaping cloud of plasma from a spherical shape near the planet, however, may be more limited than first imagined, and the metal line
cores in our best-fit model only extend to a radius of about 4.3 $R_p$ ($R_p/R_* = 0.52$) (see Figure 22). For example, for a planet like HD209458b, Figure 8 in Shaikhislamov et al. (2020) shows the distribution of C II ions responsible for C II λ1336 Å line absorption in the y-z plane, as seen in transit, simulated using their 3D atmospheric hydrodynamic escape model. This distribution is approximately spherical symmetric. Similarly, using the 3D atmospheric hydrodynamic model, Wang & Dai (2021b) showed that the absorption of the He 10830Å line by the atmosphere of WASP-69b, including regions outside the RL, is approximately spherically symmetric.

We repeat the steps in Sections 4.5 and 4.6 to compare the simulated transmission spectra with the NUV and optical observations. Although the atmosphere in the substellar direction becomes more extended with the tidal potential, this effect is offset by the suppression due to the substellar-to-terminator profile mapping described above. As a result, compared to the results of case A that do not include the tidal effect, the line center transit depths of most features in case B are very similar. This is reflected in Table 3, where the Ca II line core transit depths of case A and B are very close. Because of the higher outflow velocity of case B, the absorption profiles are broader. Therefore, although the line center transit depths of Mg II λ2796 are similar, the broader line width causes the larger case B transit depths in the Mg II lines listed in Table 3.

The Balmer lines constitute an exception. As a result of the higher neutral hydrogen column density in case B, less stellar Lyα radiation reaches the atmosphere near the RL. Furthermore, the lower atmospheric temperature within the RL leads to less collisional excitation induced Lyα radiation. The combination of these two factors causes lower Lyα intensity and thus weaker Balmer line absorption.

## 6. A PARAMETER SPACE INVESTIGATION

Although the WASP-121b system is one of the best-observed exoplanet systems, there are still uncertainties in the measured system parameters. We adjust several parameters to explore their effects on the atmospheric properties and transmission spectrum, and look for models that can reproduce all the observed spectral features within the range of their uncertainties. In the following, we first describe the quantities that we adjusted to obtain a better match to the observations, and then conduct a parameter space investigation around the best-fit model and evaluate the impact of each factor.

### 6.1. Matching the observations

We first adjust surface gravity, based on the suggestion of García Muñoz & Schneider (2019) that the mass-loss rates and absorption line profiles of UHJs are sensitive to surface gravity. WASP-121b is closer to the Roche limit than KELT-9b, the subject of their study, and due to the prevalence of RLOF, changes to surface gravity should play a greater role on WASP-121b. Therefore, the mass loss from WASP-121b is expected to be even more sensitive to surface gravity.

Compared to directly observable quantities such as the rotational period, transit duration, and the planet-to-star radius ratio, the uncertainties in the radius and mass of the planet and the star are much larger. Accordingly, when we vary planetary properties to discuss their effect on the model, we keep the rotational period, transit duration, and radius ratio unchanged by varying the radius and mass of the planet and the star as well as the orbital separation. Our intent is to explore the effect of these parameters on the transit fit, and to explore if the uncertainties in the system parameters for WASP-121b allow for a combination of parameters that lead to larger transit depths.

Because case B underestimates the observed transit depths, we search for input parameters that produce a lower planetary surface gravity. First, we discuss a case in which we keep the planet's radius unchanged but decrease the planet's mass by $1\sigma$ (5%). The surface gravity of the planet in this case is $1.4\sigma$ smaller than the canonical value. The physical parameters used by the model are listed in the lower part of Table 2. Lowering the gravity of the hydrodynamic model of the upper atmosphere further can also be achieved by using higher temperatures for our lower and middle atmospheric model, which is retrieved from secondary eclipse observations and is also subject to large uncertainty. The $1\sigma$ uncertainty of the model-dependent MCMC temperature retrieval provided by Mikal-Evans et al. (2019) is ~200 K. Here, we also uniformly increase the lower and middle atmospheric temperature by 350 K, as shown by the blue solid line in Figure 21. The resulting larger atmospheric scale height effectively increases the altitude of the $1\ \mu\text{bar}$ level.

Even with these measures in place to reduce the surface gravity and to increase the outflow velocity and mass-loss rate, the model with a solar abundance of iron falls short of the observed NUV transit depths. For a better agreement with NUV absorption, we increase the abundance of iron by a factor of 4 in the model, while retaining solar abundances for the other elements. Because Fe has almost no impact on the thermal balance of the atmosphere for models with tidal potential (see Section 7.1), this change improves the agreement of
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Simulated transmission spectrum of case D.</figcaption></figure>

the model with the NUV transmission spectrum without significantly changing atmospheric structure. A super-solar metal abundance would be in agreement with the conclusions of Evans et al. (2018) based on optical-NIR HST transmission observation. However, it would be inconsistent with the claim by Gibson et al. (2022) based on the UVES/VLT optical band transmission spectrum that the Fe/Mg abundance ratio is subsolar.

The model with the above adjustments is referred to as the case C model, and the corresponding transit depths for this model are also listed in Table 3. Except for the Balmer lines whose core depth and width are over-estimated, the simulated NUV and optical transmission spectra of case C agree with the observations.

Other relevant uncertainties include the intensity of the stellar XUV flux used in the hydrodynamic model and the Lyα flux used in the Lyα radiative transfer calculation. Huang et al. (2017) showed that the Balmer lines are much more sensitive to the stellar XUV and Lyα flux than, say, the Na lines because of their dependence on Lyα intensity inside the atmosphere. To reduce the simulated Balmer line absorption in case C, we reduce the stellar flux at wavelengths shorter than 1700 Å to 74% of the reference spectrum, shown as the magenta line in Figure 1, and weaken the stellar Lyα radiation intensity to 35% of the canonical value in the case C model. The new model with these adjustments will be referred to as the case D model. The transmission spectrum for case D, shown in Figures 22–24, can broadly reproduce all of the observed spectral features that we consider in this study.

## 6.2. *Parameter variations*

<figure><img src="image_24.png" /><figcaption>**Figure 23.** Binned simulated NUV transmission spectrum of case D compared to HST observation (Sing et al. 2019). $\chi^2/N = 1.20$.</figcaption></figure>

On the basis of the good fit to the observations, together with the fact that all the input parameters to the model are broadly consistent with their uncertainties, we select the case D model as the best-fit model. In the following, we present some other simulations that explore the parameter space more fully to illustrate the effect of the uncertainties on the interpretation of the observations. The input parameters and results of these simulations are summarized in Table 3.

For the purposes of this exercise, we create case E and F models that apply higher and lower stellar XUV fluxes, respectively, than the case D model while the other planetary system parameters remain the same. In these models, the stellar Lyα fluxes are reduced/increased accordingly to maintain a good match to the Balmer lines. The purpose of these models is to illustrate the dependency of the transit depths on the modeled stellar XUV spectrum. Model G also has the same planet parameters as the best-fit case D, but a spherically symmetric gravitational potential is applied. Model H also uses spherically symmetric potential but in this model, the planet mass is $3\sigma$ lower than the canonical value and the radius is larger by $3\sigma$. The lower atmospheric temperature in model H is also increased by 1000 K to raise the altitude at which the atmospheric pressure is $1 \mu\text{bar}$. The purpose of models G and H is to investigate whether a hotter model with lower surface gravity could match the observations without RLOF. Case I is comparable to case D, but the planet mass is only reduced by $0.5\sigma$. In this case, the lower atmospheric temperature has to be increased to match the observed features of the metal ions. Cases J and K as-
<figure><img src="image_25.png" /><figcaption>Figure 24. Comparison of optical features observed by Borsa et al. (2021) with simulated profiles of case D.</figcaption></figure>

assume a larger planet radius by 0.5σ with the same planet mass as in case D, along with adjustments to the lower and middle atmosphere temperature to fit observations.

Comparing cases C–F, we learn that either reducing the XUV flux or reducing the Lyα flux can effectively decrease the transit depths in the Balmer lines. Varying the stellar Lyα flux affects the proton and electron number densities (see Section 7), but has a minimal impact on the distribution of metal species. Ca and Mg are mostly doubly ionized, and Na and K are mostly singly ionized in the upper atmosphere focused in this work. Reducing the ionization fraction due to a lower stellar XUV flux can increase the number densities of Ca II, Mg II, Na I, and K I in the upper atmosphere. However, its effect on the modeled metal lines is offset by the reduction in atmospheric expansion caused by the lower XUV flux.

Comparing the mass-loss rates of case B and C, for a planet close to the Roche limit, a small decrease in planetary surface gravity can significantly increase the mass-loss rate. Comparing the mass-loss rates of cases C–F or case J vs. K, we learn that the reduction in the atmosphere mass-loss rate is smaller than the reduction in the stellar XUV flux. The impact of the surface gravity on the mass-loss rate is much less significant if the tidal potential is not considered, as shown by cases G and H. It is clear that without considering the tidal

potential, the metal species in the atmosphere cannot extend to the altitudes that are probed by the observation even if we significantly reduce the planet's surface gravity. This is in line with the results of Koskinen et al. (2022) who pointed out that for gas giants very close to the host star, the main mechanism of atmospheric loss is Roche Lobe overflow. For WASP-121b, heating of the middle atmosphere by the stellar bolometric flux is an important additional source of energy for atmospheric escape. This also highlights the importance of carefully treating lower and middle atmosphere in atmospheric escape models for such systems.

Although the mass-loss rate in case A is over an order of magnitude lower than in case D, the model can still reproduce the observed Hα transit profile if we reduce the stellar XUV flux or Lyα flux. The dependency of the H Balmer line transit depths on the stellar Lyα and XUV flux makes these lines an ambiguous probe of mass loss from exoplanetary atmospheres. In general, the results indicate that a single transmission feature cannot reliably constrain the properties of escaping atmospheres because of the number of degenerate free parameters in the atmosphere models. By combining many high-quality spectral features with various absorption depths, wavelengths, and absorbers, models can provide useful constraints on the atmospheric properties of UHJ, but care still needs to be taken to avoid drawing overly
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Temperature and radial velocity along the substellar stream line simulated by the hydrodynamic model for case D.</figcaption></figure>

strong conclusions due to the neglect of possible parameter degeneracies. At the same time, while models D–F, and I–K can all deliver a reasonable match to all observed features, the mass-loss rates predicted by these models vary only by a factor of about three. This is because sufficient density and velocity at high altitudes around the planet are required to match the transit observations of the metal lines, regardless of the parameter degeneracies, and the mass-loss rate is based on the density and velocity profiles. Therefore, it is not surprising that the best-fit models yield similar values for the mass-loss rate.

## 7. ATMOSPHERIC PROPERTIES AND SPECTRAL FEATURES OF THE BEST-FIT MODEL

In this section, we discuss the atmospheric properties and dominant physical processes that control the upper atmosphere of WASP-121b based on our results, relying on the case D model. We also discuss the caveat of using a 1-D model to simulate an atmosphere that undergoes Roche lobe overflow, the possible impact of stellar winds, and point to directions for future research.

### 7.1. Atmospheric Properties

The surface gravity in our case D model is $1.4\sigma$ lower than the value based on the Delrez et al. (2016) system parameters. Bourrier et al. (2020), however, reported measurements of planetary parameters for WASP-121b, with $M_p = 1.157 \pm 0.07 M_J$ and $R_p = 1.773_{-0.033}^{+0.041} R_J$. The surface gravity of our best-fit model is only 2.5% (or $0.3\sigma$) lower than the surface gravity based on this set of new parameters. Using the estimated upper limit on the mass-loss rate based on model D, we conclude that

<figure><img src="image_27.png" /><figcaption>**Figure 26.** Rate of heating (solid line) and cooling (dashed or dotted line) process for case D.</figcaption></figure>

the planet can retain its atmosphere for longer than the estimated $1.5 \pm 1$ Gyr age of the star (Delrez et al. 2016).

The temperature and bulk outflow velocity along the substellar streamline are shown in Figure 25. Within the RL, the temperature decreases significantly with altitude due to the effective adiabatic cooling caused by the strong outflow. In agreement with the analytical estimate given by Lubow & Shu (1975), the atmospheric escape velocity exceeds the sound speed slightly above the RL. In this case, the whole thermosphere (i.e., the region primarily heated by the stellar XUV radiation) is located outside of the RL. In addition, because the effective gravity points towards the star and away from the planet above the RL, when the temperature rises rapidly, pressure may increase with radius. At the L1 point, the pressure of the atmosphere is ~2.0 nbar.

This atmospheric structure is noticeably different from the WASP-121b predictions presented by Yan et al. (2021). In addition to the inclusion of a more comprehensive treatment of the chemistry and transport of metals and the inclusion of ionization and heating by excited-state H in our model, an important difference between the two models that drive the difference in atmospheric structure is the treatment of the bottom boundary conditions of the upper atmosphere and escape model. In both studies, the pressure at the bottom boundary of the escape model is 1 $\mu$bar. However, in Yan et al. (2021), the radius corresponding to this pressure is set at the planet transit radius $R_p$. In contrast, the radius of the 1 $\mu$bar level in our case D model is 1.46 $R_p$, due to the inclusion of lower and middle atmosphere structure as described in Sections 2.6 and 5.1.

The rates of the major heating and cooling processes are shown in Figure 26. From the top to the bottom of the upper atmosphere, the atomic species that con-
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Rates of H ionization and recombination of case D. Solid lines indicate processes that ionize H and dashed lines indicate processes that recombine H. As a component of the H photoionization, the rate of photoionization from $H(n = 2)$ state is shown in the black dotted line.</figcaption></figure>

tribute most to heating by photoionization are, in order, H, He, and C. The higher mass-loss rate predicted by the models with the tidal force can transport more neutral H to the upper atmosphere and partly offset photoionization. Because of absorption by neutral H atoms, the EUV radiation in such models is absorbed at a higher altitude compared to case A and thus the peak photoionization heating rate is reduced.

In contrast to the KELT-9b model of García Muñoz & Schneider (2019), photoionization of excited state H is not the main heating source at any altitudes of the atmosphere. The NUV flux of the A-type star KELT-9, which is an order of magnitude stronger than that of WASP-121, can partially explain this difference. In addition, the different choices of $Ly\alpha$ redistribution function as described in Section 2.7 and illustrated by Figure 28 may also play a role.

Due to the high mass-loss rate, adiabatic cooling dominates the cooling rate of almost the entire simulated atmosphere. Because of the steep changes in temperature, heating and cooling by advection that originate from the transport of hotter and colder fluid by the flow can also play an important role. Although radiative cooling by metal lines is negligible in case D, its role in the case A model with a lower atmospheric mass-loss rate, as shown in Figure 10, means that it can be more important for planets with less significant mass loss, for example, on KELT-9b.

The rates of the major proton production and destruction processes are shown in Figure 27. Although absorption of starlight in the Balmer continuum is not an important heating mechanism in our model, it dominates

<figure><img src="image_29.png" /><figcaption>**Figure 28.** Line profile averaged $Ly\alpha$ mean intensity calculated with Monte Carlo radiative transfer for case D. The contribution originated from stellar irradiation, H recombination, and electron collisional excitation of $H(1s)$ are shown in black, blue, and red dotted lines, respectively. The sum of them is shown in the magenta dotted line. Being smoothed using B-spline as mentioned in Section 2.7, the profile of this overall $Ly\alpha$ mean intensity is shown as the green solid line. The result based on the complete frequency redistribution function, instead of the Hummer-IIB partial redistribution function, is shown by the yellow line.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 29.** Number density distributions of major species along the substellar direction from both the case D upper atmosphere and the corresponding lower atmosphere model.</figcaption></figure>

proton production below $2R_{px}$. As shown in Figure 28, stellar $Ly\alpha$ radiation is typically the main source of $Ly\alpha$ intensity near the base of the atmosphere. As a result, the proton number density in this region is sensitive to the stellar $Ly\alpha$ flux.
Proton charge exchange with Fe can be more impor-
tant in the removal of H than recombination near the
base of the hydrodynamic model. Fe I and Mg I number
densities increase with the sharp increase in temperature
between 3 and $4R_p$, as shown in Figure 29. The higher
temperature increases the rates of charge exchange with
hydrogen and drives the ionization fractions of Fe and
Mg closer to the ionization fraction of H, which is more
neutral. It is worth noting that the charge exchange
rate between $Fe^+$ and H used in the calculation orig-
inates from a cursory theoretical estimation. Consid-
ering the large impact of this charge exchange rate on
the model, experimental measurement or careful theo-
retical consideration of this rate would be helpful to the
understanding of the planetary upper atmosphere.

7.2. *Stellar wind and diurnal temperature variations*

Strong stellar wind may modulate atmospheric out-
flow and shape the distribution of the escaping plasma
away from the planet (Wang & Dai 2018). Here we
discuss some aspects of the potential impact of the stel-
lar wind on WASP-121b based on model D. We note
that the stellar wind properties of WASP-121 are uncon-
strained by observations. Instead, we use an estimate of
the stellar mass-loss rate and the scaling proposed by
Preusse et al. (2005) to estimate the thermal, ram, and
magnetic pressures of the wind at the orbital distance
of WASP-121b.

Delrez et al. (2016) noted that WASP-121 has similar
stellar type and RV jitter as τ Bootis A. The Mount-
Wilson S-index of τ Boo is 0.185 (Mengel et al. 2016),
corresponding to a log 10(R'_{HK}) of -4.70. This larger
activity index compared to WASP-121, which has a
log 10(R'_{HK}) of -4.81 (Borsa et al. 2021), indicates that
τ Boo is more active than WASP-121. As an upper
limit, we use the stellar wind parameters of τ Boo
(Nicholson et al. 2016) to estimate its impact on the
planetary atmosphere. We have the stellar wind mass-
loss rate of 2.3 × 10⁻¹² M☉yr⁻¹ with a constant veloc-
ity of u = 220km s⁻¹. This yields the plasma den-
sity of nₚ(0.02545AU) = 2.2 × 10⁶ cm⁻³. Assuming
magnetic field strength varying as r⁻² and B = 0.03
g at 0.049 AU from the host star, as the value at
τ Boo b, we have B(0.02545AU) = 0.12 G. Further-
more, we estimate the electron temperature by scaling as
Tₚ(0.02545AU) = 10⁶ K(r/0.049AU)⁻⁰.⁷ = 1.6 × 10⁶ K
and Tₑ = Tₚ * 2 = 3.2 × 10⁶ K. Combining these, the
thermal pressure of the stellar wind at the distance of
the planet Pth = 1.5 × 10⁻³ µbar, the ram pressure
Pram = ρ(v²ₚₗd + v²ₚₗb) = 3.5 × 10⁻³ µbar, and the
magnetic pressure PB = 6 × 10⁻⁴ µbar. Therefore, we

have the total stellar wind pressure at the planet should
be less than 5.6 × 10⁻³ µbar.

WASP-121b has a near-polar orbit (Delrez et al. 2016) and is therefore subjected to the stellar wind in the polar direction along its orbit, an area that has been studied to a lesser extent. In the absence of any information on the spatial structure of the stellar wind in systems like WASP-121b, we estimate the dependence of the stellar wind pressure on latitude based on observations of the solar wind by Ulysses (Verscharen et al. 2021). During solar minimum, the solar wind velocity is higher in the polar direction compared to the equatorial direction, while the number density is lower in the polar direction. Verscharen et al. (2021) suggested that the momentum flux ($n_p m_p r^2 v^2$) and sonic Mach number ($v^2/T$) exhibit little latitude variation during solar minimum. Therefore, it is expected that the ram pressure and thermal pressure do not change dramatically with solar latitude.

The total stellar wind pressure balances the pressure
of planetary wind in model D at 7.5 $R_{px}$. In compari-
son, as we noted before, the line cores of the Mg II lines,
which probe the highest altitudes, reach the effective ra-
dius of 4.3$R_p$. At 4.3$R_{px}$ along the substellar direction,
the total atmospheric pressure is 8.5 × 10⁻³ µbar. Addi-
tionally, the L1 point, which roughly coincides with the
sonic point along the substellar streamline, is at about
1.7$R_{px}$. Therefore, unless the magnetic pressure of the
stellar wind of WASP-121 is significantly stronger, it
is unlikely to have a significant impact on the dayside
mass-loss rate estimated in this model. This is consis-
tent with the recent conclusion by Mitani et al. (2022).
Along the terminator, the stellar wind would push mate-
rial away from the planet and potentially enhance mass-
loss rates, although its relative importance is limited by
the relatively low stellar wind pressure, as the thermal
escape rates are already high. The stellar wind can also
shape the density and velocity distribution of the escap-
ing material. Modeling this is beyond the scope of the
current work.

Our globally averaged approach and approximations
of the atmosphere under Roche lobe overflow based on
a 1-D model ignore day-night variations of density and
temperature in the atmosphere. Using a phase curve ob-
served with the HST/WFC3 instrument in the near-IR,
Mikal-Evans et al. (2022) inferred diurnal temperature
variations in the lower atmosphere of WASP-121b in the
pressure range of 1–100 mbar. The results are supported
by more recent observations with the JWST NIRSpec
instrument (Mikal-Evans et al. 2023). As shown in Fig-
ure 8, the lower atmosphere temperature applied in this
study is roughly comparable to their dayside temper-
ature, although it lacks the curious steep increase with
altitude and is therefore likely to be more representative
of an average temperature profile.

The much lower night side temperature observed in
the near-IR is below the condensation temperature of
Mg and Fe, suggesting that mineral clouds can form on
the night side. As noted by Mikal-Evans et al. (2022),
however, the clear detection of escaping Mg II and Fe II,
with absorption lines that also trace the middle and up-
per atmosphere, indicates that these species are present
at the day-night terminator in the gas phase. This is
possible, for example, if efficient mixing acts to prevent
day-night cold trapping in the lower atmosphere. We
point out that escape on the dayside and along the ter-
minator must also be strong because the NUV spectrum
cannot be explained without sufficient densities and ve-
locities at high altitudes around the planet.

In addition, we note that the near-IR phase curve is
probing the lower atmosphere whereas our focus is on
the upper atmosphere. In the extended upper atmo-
sphere, radiation can penetrate to the night side and
leave a relatively small region around the anti-stellar
point in the shadow. In addition, WASP-121b is very
close to the host star and stellar rays incident on it are
not all parallel. The maximum angle of rays incident on
the planet's terminator is about 13–17° from the parallel
direction. The diffuse radiation field, together with the
extent of the upper atmosphere, can maintain ionization
and heating along the terminator and on the night side.
Should the nightside mass-loss rate be reduced, however,
the extreme assumption of no mass loss from the night
side would lead to a reduction of the global mass-loss
rate by a factor of two. Since our radiation field is glob-
ally averaged, this geometrical factor is already partly
accounted for. Beyond that, the mass-loss rate on the
dayside can be further reduced by day-night winds that
redirect flows and can extract energy available to power
mass loss, although we suspect that the very strong ef-
fect of RLOF on the dayside of WASP-121b makes this
of limited importance.

7.3. *Spectral Features*

As we have shown, velocity broadening due to planetary rotation and escape makes important contributions to the interpretation of the observations. Figure 30 illustrates the impact of line broadening caused due to the atmospheric bulk outflow velocity and the rotation of the planet, which we discussed in Section 4.3. The red solid line shows the line profile based on the best-fit model with thermal and natural broadening only. It is only about 1/10 of the width of the observation, which indicates that thermal and natural broadening alone cannot explain the width of the observed line profile.

<figure><img src="image_31.png" /><figcaption>**Figure 30.** Comparison of Ca II Kα3934 Å line profiles given by different treatments for velocity broadening caused by atmospheric outflow and planet rotation. The black solid line shows the enlarged profile of the same best-fit model as in Figure 24. The red solid line shows the line profile that only includes thermal broadening. In addition to thermal broadening, the magenta and blue solid lines show the line profile that also includes the effect of rotation and outflow velocity, respectively.</figcaption></figure>

When rotation is included without radial outflow,
the line profile exhibits the double-peaked structure de-
scribed in Section 4.3. In this case, because absorption
at a given wavelength originates in narrow strips, all
at effectively different line-of-sight velocities, shown in
Figure 15, the transit depth is shallower due to the rel-
atively small area of the strips. Introducing broadening
due to the outflow velocity, the line profile is closer to
the observed line profile.

When some of the key number density profiles along
the substellar direction of the case D model shown in
Figure 29 are compared to those of the case A model
shown in Figure 12, the excited state H(2p) profile in
the case D model has a lower peak number density and a
shallower decrease with radius. This means that WASP-
121b has an extended excited-state hydrogen distribu-
tion, which is marginally optically thick at Hα. This
coincides with the physical image reflected by the ob-
served large line center absorption depths ratio between
the Balmer lines, which are close to the ratio of their
spontaneous transition rates. In case D, because the
absorption originates from material at higher altitudes
and with higher outflow velocity, the velocity broaden-
ing leads to a feature that is broader than that of case
A and somewhat broader than the observed line profile.
This difference, combined with the fact that the model
slightly underestimates the absorption depth at Hβ may
<figure><img src="image_32.png" /><figcaption>**Figure 31.** Line profile of Fe II λ5170Å, Fe I λ5169Å, and Mg I λ5169Å, λ5174Å for case D.</figcaption></figure>

imply that the slope of the H(2p) number density profile should be slightly steeper than it is in the case D model.

Sing et al. (2019) did not detect absorption by Mg I. However, as shown in Figure 22, the Mg I λ2853Å line is clearly visible with a line center transit depth of $R_p/R_* = 0.275$. Compared to the Mg II lines, the absorption by Mg I originates from a region at lower altitudes where the outflow velocity and rotational velocity are lower, and Mg I λ2853Å absorption has a much narrower line width. As a result, after binning to 4 Å wavelength band, Figure 23 shows that the non-detection of Mg I λ2853Å line is not statistically significant.

Comparing the simulated NUV spectra with observations, the biggest discrepancy is that the model fails to explain the absorption signal observed at 2438Å and 2485Å. Sing et al. (2019) suggested the absorption feature at 2485 Å may be due to the Fe I λ2484Å line. The Fe I λ2484Å line, which is weaker than the Mg I λ2853Å line, however, is not strong or broad enough to match the observations with a bin width of 4 Å.

Other notable features in the simulated spectrum are the Fe II λ5170Å ($3d^54s^2 \, ^6S \rightarrow 3d^64p \, ^6P^o$) and the Fe I λ5169Å ($3d^74s \, ^3F \rightarrow 3d^64s4p \, ^3D^o$) lines shown in Figure 31, positioned between the Mg I triplets. In particular, the line center transit depth of the Fe II λ5170Å line reaches the RL in case D. The lower levels of the two Fe transitions are at the energy of 2.9 eV for Fe II and 1.5 eV for Fe I. The intensities of these two features are sensitive to the temperature profile and radiative intensity within the atmosphere. For example, case J, where the stellar EUV flux is lower compared to case D, predicts a much weaker Fe II λ5170Å line with a stronger Fe I λ5169Å line. We note that Borsa et al. (2021) detected the Mg I triplet while they did

not mention the Fe I lines. Pai Asnodkar et al. (2022), on the other hand, resolve individual Fe II features in the optical transmission spectrum and include the Fe II λ5170 Å line displayed in Figure 31, but their results are for KELT-9b. If our model significantly overestimates the absorption by these two Fe I and II lines, the most likely explanation is our assumption that the energy levels of Fe satisfy the Boltzmann distribution. If these lines remain undetected, this assumption should be adjusted in future work.

## 8. CONCLUSIONS

In this work, we present a hydrodynamic escape model of the upper atmosphere that includes excited-state H and metals such as Mg and Fe to interpret the NUV and optical transit observations of the ultra-hot Jupiter WASP-121b. The density profiles predicted by the model are combined with density profiles from a photochemical model for the lower and middle atmosphere to provide a consistent, one-dimensional, globally averaged model for the whole atmosphere. Based on the atmospheric structure in this model and adjusting some of the free parameters in the model, we can broadly reproduce the observed NUV spectrum and some of the key optical absorption lines of WASP 121b within uncertainties if we enhance the densities and velocities of the escaping material based on Roche lobe overflow. This allows us to place firmer constraints on the mass-loss rate from one of the most extreme hot Jupiters among the known transiting planets.

As expected, the model demonstrates the importance of the high mass-loss rate caused by Roche lobe (RL) overflow for this extremely close-in UHJ. In particular, the deep NUV Fe II transit depths and the broad Ca II line widths are difficult to explain with thermal broadening alone or without invoking RL overflow. On the other hand, the Na D and K lines probe the atmosphere below the Roche lobe and are not good tracers of escape. At the same time, the transit depths in the $H\alpha$ and $H\beta$ lines are sensitive to the incident stellar radiation field and do not generally provide an unambiguous probe of the properties of the escaping material.

Our model demonstrates that careful consideration of the lower and middle atmosphere that provides the lower boundary conditions for the escape model is important for interpreting the observations and correctly predicting mass-loss rates. For planets such as WASP-121b that lie close to the fundamental stability limit for rapid Roche lobe overflow of the middle atmosphere (Koskinen et al. 2022), this is especially important. Here, we also show that the mass-loss rate from planets like WASP-121b is very sensitive to surface gravity.
According to the best-fit case D model, photoionization of excited H(n=2) plays an important role in determining the overall ionization fraction of H, especially at radii less than about 2 R<sub>p</sub>. In contrast to recent work on KELT-9b by García Muñoz & Schneider (2019), we do not find it to be an important heating mechanism in the atmosphere of WASP-121b. On average, the energy input per ionization event is much lower than, for example, from ionization of ground-state H. In addition, the upper atmospheric structure on WASP-121b differs significantly from KELT-9b due to strong cooling induced by rapid escape, and the input stellar spectrum is very different from KELT-9b. Rapid escape also means that adiabatic cooling overshadows radiative cooling by metal species on WASP-121b. We note, however, that radiative cooling by metal lines can be important on planets with a lower mass-loss rate, and may be constrained by similar observations of such planets.

In addition, our results highlight the fact that a single transmission feature cannot exclusively determine upper atmosphere structure and, therefore, the mass-loss rate from a planet. There are a large number of degenerate free parameters in the models, and different atmospheres can produce the same transmission profile for a given feature in spatially unresolved observations. Combining high-quality spectral features of different species and covering different transit depths and wavelengths simultaneously can effectively break the degeneracy and provide a more reliable constraint on the atmospheric structure. Even so, for many planets, it can still be hard to resolve parameter degeneracy in the fits, even when several features are available. For example, we found multiple combinations that can fit the observations of WASP-121b with our 1D model. Comparing the line shape of high S/N line profile measurements with detailed 3D simulations would help to impose better constraints on the atmosphere from the observations. More accurate measurement of the system parameters can also significantly benefit the effort to interpret observations of atmospheric escape.

The models of the atmosphere presented in this work have their limitations, and the chief among them is that they are one-dimensional. One-dimensional models obviously cannot represent multidimensional physics such

as the Coriolis effect, acceleration by stellar radiation pressure, the interaction with the stellar wind, and the effects of the stellar or planetary magnetic field. Although we show the latter three factors are unlikely to significantly change the estimated mass-loss rate, three-dimensional dynamics is likely necessary to explain the observed blueshift of the optical absorption lines and the asymmetry of the Ca II line profiles that may be due to horizontal winds driven by day-night temperature differences in the planet's atmosphere and/or by radiation pressure or stellar wind in the extended exosphere. It would still be computationally difficult to incorporate the detailed physical processes discussed in this work into 3D atmospheric models without additional simplifications that would require more analysis of the 1D results to justify. While the observations of WASP-121b clearly provide strong motivation for multidimensional modeling, such modeling is therefore beyond the scope of the present work. Instead, our results provide an overview of the relevant physics and chemistry that future multi-dimensional models can use as guidance while in development.

We note that the mass-loss rate of (4–9) × 10<sup>13</sup> g s<sup>-1</sup> based on our 1D models is likely to differ from a similar 3D model by a limited factor of a few, or at most an order of magnitude (Koskinen et al. 2022). Furthermore, despite the degeneracy in the parameters when fitting the observations, the mass-loss rates based on the best-fit combinations that we explored are consistent roughly within a factor of three. In all cases, the models without the RL enhancement on the mass loss, therefore, fail to match the observations. This indicates that WASP-121b is likely to suffer, and to have suffered, significant envelope loss during its main sequence lifetime.

## ACKNOWLEDGEMENTS

We thank Orly Gnat for helping us to set up the Cloudy script to reproduce their results. C.H. thanks Xiao Fang for the useful information about the transitional rates of Ca I. We also thank David Sing for a helpful discussion. T.T.K and C.H. acknowledge support by the NASA Exoplanet Research Program through grant 80NSSC18K0569.

## REFERENCES

Anderson, H., Ballance, C. P., Badnell, N. R., & Summers, H. P. 2002, Journal of Physics B Atomic Molecular Physics, 35, 1613, doi: 10.1088/0953-4075/35/6/701

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
Azevedo Silva, T., Demangeon, O. D. S., Santos, N. C., et al. 2022, A&A, 666, L10, doi: 10.1051/0004-6361/20224489

Ballester, G. E., & Ben-Jaffel, L. 2015, ApJ, 804, 116, doi: 10.1088/0004-637X/804/2/116

Barklem, P. S. 2016, PhRvA, 93, 042705, doi: 10.1103/PhysRevA.93.042705

Barklem, P. S., Belyaev, A. K., Spielfiedel, A., Guitou, M., & Feautrier, N. 2012, A&A, 541, A80, doi: 10.1051/0004-6361/201219081

Ben-Yami, M., Madhusudhan, N., Cabot, S. H. C., et al. 2020, ApJL, 897, L5, doi: 10.3847/2041-8213/ab94aa

Black, J. H. 1981, MNRAS, 197, 553, doi: 10.1093/mnras/197.3.553

Borsa, F., Allart, R., Casasayas-Barris, N., et al. 2021, A&A, 645, A24, doi: 10.1051/0004-6361/202039344

Bourrier, V., Ehrenreich, D., Lendl, M., et al. 2020, A&A, 635, A205, doi: 10.1051/0004-6361/201936640

Cabot, S. H. C., Madhusudhan, N., Welbanks, L., Piette, A., & Gandhi, S. 2020, MNRAS, 494, 363, doi: 10.1093/mnras/staa748

Cecchi-Pestellini, C., Ciaravella, A., Micela, G., & Penz, T. 2009, A&A, 496, 863, doi: 10.1051/0004-6361/200809955

Chenel, A., Mangaud, E., Justum, Y., et al. 2010, Journal of Physics B Atomic Molecular Physics, 43, 245701, doi: 10.1088/0953-4075/43/24/245701

Cunto, W., Mendoza, C., Ochsenbein, F., & Zeippen, C. J. 1993, A&A, 275, L5

Dalgarno, A., & Williams, D. A. 1965, Proceedings of the Physical Society, 85, 685, doi: 10.1088/0370-1328/85/4/308

Del Zanna, G., Dere, K. P., Young, P. R., & Landi, E. 2021, ApJ, 909, 38, doi: 10.3847/1538-4357/abd8ce

Delrez, L., Santerne, A., Almenara, J. M., et al. 2016, MNRAS, 458, 4025, doi: 10.1093/mnras/stw522

Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149, doi: 10.1051/aas:1997368

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium

Dutta, C. M., Nordlander, P., Kimura, M., & Dalgarno, A. 2001, PhRvA, 63, 022709, doi: 10.1103/PhysRevA.63.022709

Dwivedi, N. K., Khodachenko, M. L., Shaikhislamov, I. F., et al. 2019, MNRAS, 487, 4208, doi: 10.1093/mnras/stz1345

Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960

Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459, doi: 10.1038/nature14501

Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329, doi: 10.1051/0004-6361:20066929

Evans, T. M., Sing, D. K., Wakeford, H. R., et al. 2016, ApJL, 822, L4, doi: 10.3847/2041-8205/822/1/L4

Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ, 156, 283, doi: 10.3847/1538-3881/aaebff

Ferland, G. J. 1996, Hazy, A Brief Introduction to Cloudy 90

Fišák, J., Kubát, J., & Krtička, J. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 508, The B[e] Phenomenon: Forty Years of Studies, ed. A. Miroshnichenko, S. Zharikov, D. Korčáková, & M. Wolf, 115

Fossati, L., Koskinen, T., Lothringer, J. D., et al. 2018, ApJL, 868, L30, doi: 10.3847/2041-8213/aaf0a5

Gaia Collaboration. 2020, VizieR Online Data Catalog, I/350

García Muñoz, A., & Schneider, P. C. 2019, ApJL, 884, L43, doi: 10.3847/2041-8213/ab498d

Gibson, N. P., Nugroho, S. K., Lothringer, J., Maguire, C., & Sing, D. K. 2022, MNRAS, 512, 4618, doi: 10.1093/mnras/stac091

Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215, doi: 10.1093/mnras/staa228

Glover, S. C. O., & Jappsen, A. K. 2007, ApJ, 666, 1, doi: 10.1086/519445

Gnat, O., & Ferland, G. J. 2012, ApJS, 199, 20, doi: 10.1088/0067-0049/199/1/20

Gordon, S., & McBride, B. J. 1999, Thermodynamic data to 20,000 K for monatomic gases (Glenn Research Center)

Guo, J. H. 2013, ApJ, 766, 102, doi: 10.1088/0004-637X/766/2/102

Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020, A&A, 641, A123, doi: 10.1051/0004-6361/202038365

Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306, doi: 10.1086/167595

Huang, C., Arras, P., Christie, D., & Li, Z.-Y. 2017, ApJ, 851, 150, doi: 10.3847/1538-4357/aa9b32

Hummer, D. G. 1962, MNRAS, 125, 21, doi: 10.1093/mnras/125.1.21

Jefferies, J. T. 1968, Spectral line formation (Blaisdell)

John, T. L. 1988, A&A, 193, 189

Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2019, ApJ, 885, 67, doi: 10.3847/1538-4357/ab46a4

Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205, doi: 10.1086/192335

Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013a, Icarus, 226, 1678, doi: 10.1016/j.icarus.2012.09.027
Koskinen, T. T., Lavvas, P., Huang, C., et al. 2022, The Astrophysical Journal, 929, 52, doi: 10.3847/1538-4357/ac4f45

Koskinen, T. T., Strobel, D. F., & Brown, Z. 2021, Icarus, 362, 114396, doi: 10.1016/j.icarus.2021.114396

Koskinen, T. T., Yelle, R. V., Harris, M. J., & Lavvas, P. 2013b, Icarus, 226, 1695, doi: 10.1016/j.icarus.2012.09.026

Kramida, A., Yu. Ralchenko, Reader, J., & NIST ASD Team. 2020, NIST Atomic Spectra Database (ver. 5.8), [Online]. Available: https://physics.nist.gov/asd [2017, April 9]. National Institute of Standards and Technology, Gaithersburg, MD.

Kwolek, J. M., Goodman, D. S., Slayton, B., et al. 2019, PhRvA, 99, 052703, doi: 10.1103/PhysRevA.99.052703

Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13, doi: 10.1051/0004-6361/201937175

Landini, M., & Fossi, B. C. M. 1991, A&AS, 91, 183

Lavvas, P., Koskinen, T., & Yelle, R. V. 2014, ApJ, 796, 15, doi: 10.1088/0004-637X/796/1/15

Lee, H.-W., & Kim, H. I. 2004, MNRAS, 347, 802, doi: 10.1111/j.1365-2966.2004.07255.x

Lemaire, P., Emerich, C., Vial, J. C., et al. 2005, Advances in Space Research, 35, 384, doi: 10.1016/j.asr.2004.11.004

Lin, C. Y., Stancil, P. C., Gu, J. P., Buenker, R. J., & Kimura, M. 2005, PhRvA, 71, 062708, doi: 10.1103/PhysRevA.71.062708

Linsky, J. L., Fontenla, J., & France, K. 2014, ApJ, 780, 61, doi: 10.1088/0004-637X/780/1/61

Lubow, S. H., & Shu, F. H. 1975, ApJ, 198, 383, doi: 10.1086/153614

Mamajek, E. E., Prsa, A., Torres, G., et al. 2015, arXiv e-prints, arXiv:1510.07674. https://arxiv.org/abs/1510.07674

Mengel, M. W., Fares, R., Marsden, S. C., et al. 2016, MNRAS, 459, 4325, doi: 10.1093/mnras/stw828

Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2021, MNRAS, 506, 3853, doi: 10.1093/mnras/stab1878

Mikal-Evans, T., Sing, D. K., Goyal, J. M., et al. 2019, MNRAS, 488, 2222, doi: 10.1093/mnras/stz1753

Mikal-Evans, T., Sing, D. K., Barstow, J. K., et al. 2022, Nature Astronomy, 6, 471, doi: 10.1038/s41550-021-01592-w

Mikal-Evans, T., Sing, D. K., Dong, J., et al. 2023, ApJL, 943, L17, doi: 10.3847/2041-8213/acb049

Miroshnichenko, I. B., Shaikhislamov, I. F., Berezutskii, A. G., Rumenskikh, M. S., & Vetrova, E. S. 2021, Astronomy Reports, 65, 61, doi: 10.1134/S1063772921010030

Mitani, H., Nakatani, R., & Yoshida, N. 2022, MNRAS, 512, 855, doi: 10.1093/mnras/stac556

Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23, doi: 10.1088/0004-637X/693/1/23

Nahar, S. N. 1997, PhRvA, 55, 1980, doi: 10.1103/PhysRevA.55.1980

Nahar, S. N., Bautista, M. A., & Pradhan, A. K. 1997, ApJ, 479, 497, doi: 10.1086/303874

Neufeld, D. A., & Dalgarno, A. 1987, PhRvA, 35, 3142, doi: 10.1103/PhysRevA.35.3142

Nicholson, B. A., Vidotto, A. A., Mengel, M., et al. 2016, MNRAS, 459, 1907, doi: 10.1093/mnras/stw731

Oklopčić, A. 2019, ApJ, 881, 133, doi: 10.3847/1538-4357/ab2f7f

Osorio, Y., Barklem, P. S., Lind, K., et al. 2015, A&A, 579, A53, doi: 10.1051/0004-6361/201525846

Osorio, Y., Lind, K., Barklem, P. S., Allende Prieto, C., & Zatsarinny, O. 2019, A&A, 623, A103, doi: 10.1051/0004-6361/201834680

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei

Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29, doi: 10.3847/1538-4357/aa890a

Pai Asnodkar, A., Wang, J., Gaudí, B. S., et al. 2022, AJ, 163, 40, doi: 10.3847/1538-3881/ac32c7

Prasad, S. S., & Huntress, W. T., J. 1980, ApJS, 43, 1, doi: 10.1086/190665

Preusse, S., Kopp, A., Büchner, J., & Motschmann, U. 2005, A&A, 434, 1191, doi: 10.1051/0004-6361:20041680

Rohatgi, A. 2021, Webplotdigitizer: Version 4.5. https://automeris.io/WebPlotDigitizer

Rutherford, J. A., & Vroom, D. A. 1972, JChPh, 57, 3091, doi: 10.1063/1.1678724

Salz, M., Schneider, P. C., Fossati, L., et al. 2019, A&A, 623, A57, doi: 10.1051/0004-6361/201732419

Satta, M., Grassi, T., Gianturco, F. A., Yakovleva, S. A., & Belyaev, A. K. 2013, MNRAS, 436, 2722, doi: 10.1093/mnras/stt1771

Seaton, M. J. 1955, Proceedings of the Physical Society A, 68, 457, doi: 10.1088/0370-1298/68/6/301

Seon, K.-i., & Kim, C.-G. 2020, ApJS, 250, 9, doi: 10.3847/1538-4365/aba2d6

Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2020, MNRAS, 491, 3435, doi: 10.1093/mnras/stz3211

Shull, J. M., & van Steenberg, M. 1982, ApJS, 48, 95, doi: 10.1086/190769

Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, C., & Weiss, W. W. 2004, A&A, 428, 993, doi: 10.1051/0004-6361:20034169

Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91, doi: 10.3847/1538-3881/ab2986
Smith, F. L., I., & Smith, C. 1972, J. Geophys. Res., 77, 3592, doi: 10.1029/JA077i019p03592

Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68, doi: 10.1038/s41586-018-0067-5

Stancil, P. C., Schultz, D. R., Kimura, M., et al. 1999, A&AS, 140, 225, doi: 10.1051/aas:1999419

Stancil, P. C., Havener, C. C., Krstić, P. S., et al. 1998, ApJ, 502, 1006, doi: 10.1086/305937

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41, doi: 10.1093/mnras/272.1.41

Storey, P. J., & Sochi, T. 2015, MNRAS, 446, 1864, doi: 10.1093/mnras/stu2243

Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, The Astrophysical Journal, 621, 1049, doi: 10.1086/427204/FULLTEXT/

van Regemorter, H. 1962, ApJ, 136, 906, doi: 10.1086/147445

Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467, doi: 10.1086/192284

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435

Verner, D. A., & Yakovlev, D. G. 1995, A&AS, 109, 125

Verner, E. M., Verner, D. A., Korista, K. T., et al. 1999, ApJS, 120, 101, doi: 10.1086/313171

Verscharen, D., Bale, S. D., & Velli, M. 2021, MNRAS, 506, 4993, doi: 10.1093/mnras/stab2051

Voronov, G. S. 1997, Atomic Data and Nuclear Data Tables, 65, 1, doi: 10.1006/adnd.1997.0732

Wang, L., & Dai, F. 2018, ApJ, 860, 175, doi: 10.3847/1538-4357/aac1c0

—. 2021a, ApJ, 914, 99, doi: 10.3847/1538-4357/abf1ed

—. 2021b, ApJ, 914, 98, doi: 10.3847/1538-4357/abf1ee

Watanabe, A., Dutta, C. M., Nordlander, P., Kimura, M., & Dalgarno, A. 2002, PhRvA, 66, 044701, doi: 10.1103/PhysRevA.66.044701

Wood, B. E., Redfield, S., Linsky, J. L., Müller, H.-R., & Zank, G. P. 2005, ApJS, 159, 118, doi: 10.1086/430523

Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197, doi: 10.1051/0004-6361:20064981

Yan, D., Guo, J., Huang, C., & Xing, L. 2021, ApJL, 907, L47, doi: 10.3847/2041-8213/abda41

Yan, D., Seon, K.-i., Guo, J., Chen, G., & Li, L. 2022, ApJ, 936, 177, doi: 10.3847/1538-4357/ac8793

Yan, F., & Henning, T. 2018, Nature Astronomy, 2, 714, doi: 10.1038/s41550-018-0503-3

Yelle, R. V. 2004, Icarus, 170, 167, doi: 10.1016/j.icarus.2004.02.008

Yu, Y., & Derevianko, A. 2018, Atomic Data and Nuclear Data Tables, 119, 263, doi: 10.1016/j.adt.2017.03.002

Zatsarinny, O., Bartschat, K., Fernandez-Menchero, L., & Tayal, S. S. 2019, PhRvA, 99, 023430, doi: 10.1103/PhysRevA.99.023430

Zhao, L. B., Stancil, P. C., Gu, J. P., et al. 2005, PhRvA, 71, 062713, doi: 10.1103/PhysRevA.71.062713

—. 2004, ApJ, 615, 1063, doi: 10.1086/424729
# APPENDIX

## A. LIST OF CHARGE EXCHANGE RATES

See Table 4.

<table><caption>**Table 4.** Charge exchange rates</caption><thead><tr><th>Reactants</th><th>Rate<sup>a</sup> (cm<sup>3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>He + H<sup>+</sup> <sup>b</sup></td><td>1.75 × 10<sup>-11</sup>(T/300)<sup>-0.75</sup> exp(-12.75/T<sub>4</sub>)</td></tr><tr><td>He<sup>+</sup> + H <sup>b</sup></td><td>1.25 × 10<sup>-15</sup>(T/300)<sup>0.25</sup></td></tr><tr><td>Mg + H<sup>+</sup> <sup>c</sup></td><td>9.76 × 10<sup>-12</sup>T<sub>4</sub><sup>3.14</sup>[1 + 55.54 exp(-1.12T<sub>4</sub>)]</td></tr><tr><td>Mg<sup>+</sup> + H <sup>c</sup></td><td>2.95 × 10<sup>-12</sup>T<sub>4</sub><sup>3.28</sup>[1 + 55.54 exp(-1.12T<sub>4</sub>)] exp(-6.91/T<sub>4</sub>)</td></tr><tr><td>Mg<sup>+</sup> + H<sup>+</sup> <sup>c</sup></td><td>7.6 × 10<sup>-14</sup>[1 - 1.97 exp(-4.32T<sub>4</sub>)] exp(-1.67/T<sub>4</sub>)</td></tr><tr><td>Mg<sup>2+</sup> + H <sup>c</sup></td><td>8.58 × 10<sup>-14</sup>T<sub>4</sub><sup>2.49×10<sup>-3</sup></sup>[1 + 0.0293 exp(-4.33T<sub>4</sub>)]</td></tr><tr><td>Fe + H<sup>+</sup> <sup>d</sup></td><td>4.0 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup> + H <sup>d</sup></td><td>1.16 × 10<sup>-9</sup>(T/300)<sup>0.072</sup> exp(-6.61/T<sub>4</sub>)</td></tr><tr><td>Fe<sup>+</sup> + H<sup>+</sup> <sup>e</sup></td><td>2.3 × 10<sup>-9</sup> exp(-3.0/T<sub>4</sub>)</td></tr><tr><td>Fe<sup>2+</sup> + H <sup>e</sup></td><td>1.26 × 10<sup>-9</sup>T<sub>4</sub><sup>0.0772</sup>[1 - 0.41 exp(-7.31T<sub>4</sub>)]</td></tr><tr><td>Si + H<sup>+</sup> <sup>b</sup></td><td>7.41 × 10<sup>-11</sup>(T/300)<sup>0.85</sup></td></tr><tr><td>Si<sup>+</sup> + H <sup>b</sup></td><td>4.71 × 10<sup>-11</sup>(T/300)<sup>0.95</sup> exp(-6.32/T<sub>4</sub>)</td></tr><tr><td>Si<sup>+</sup> + H<sup>+</sup> <sup>c</sup></td><td>4.1 × 10<sup>-10</sup>T<sub>4</sub><sup>0.24</sup>[1 + 3.17 exp(-4.18 × 10<sup>-3</sup>T<sub>4</sub>)] exp(-3.18/T<sub>4</sub>)</td></tr><tr><td>Si<sup>2+</sup> + H <sup>c</sup></td><td>1.26 × 10<sup>-9</sup>T<sub>4</sub><sup>0.24</sup>[1 + 3.17 exp(-4.18 × 10<sup>-3</sup>T<sub>4</sub>)]</td></tr><tr><td>Si + He<sup>+</sup> <sup>f</sup></td><td>3.32 × 10<sup>-13</sup>√T + 1.2 × 10<sup>-16</sup>T + 4.2 × 10<sup>-9</sup>/√T - 7.9 × 10<sup>-13</sup></td></tr><tr><td>Si<sup>+</sup> + He <sup>f</sup></td><td>1.415 × (3.32 × 10<sup>-13</sup>√T + 1.2 × 10<sup>-16</sup>T + 4.2 × 10<sup>-9</sup>/√T - 7.9 × 10<sup>-13</sup>)T<sup>0.103</sup> exp(-19.1/T<sub>4</sub>)</td></tr><tr><td>O + H<sup>+</sup> <sup>g</sup></td><td>2.08 × 10<sup>-9</sup>T<sub>4</sub><sup>0.405</sup> + 1.11 × 10<sup>-11</sup>T<sub>4</sub><sup>-0.458</sup></td></tr><tr><td>O<sup>+</sup> + H <sup>g</sup></td><td>(1.26 × 10<sup>-9</sup>T<sub>4</sub><sup>0.517</sup> + 4.25 × 10<sup>-10</sup>T<sub>4</sub><sup>6.69×10<sup>-3</sup></sup>) exp(-227/T)</td></tr><tr><td>C + H<sup>+</sup> <sup>h</sup></td><td>1.31 × 10<sup>-15</sup>(T/300)<sup>0.213</sup></td></tr><tr><td>C<sup>+</sup> + H <sup>h</sup></td><td>6.3 × 10<sup>-17</sup>(T/300)<sup>1.96</sup> exp(-17/T<sub>4</sub>)</td></tr><tr><td>C + He<sup>+</sup> <sup>b</sup></td><td>2.5 × 10<sup>-15</sup>(T/300)<sup>1.597</sup></td></tr><tr><td>C<sup>+</sup> + He <sup>b</sup></td><td>6.75 × 10<sup>-15</sup>(T/300)<sup>1.654</sup> exp(-15.5/T<sub>4</sub>)</td></tr><tr><td>C + Si<sup>+</sup> <sup>f</sup></td><td>(0.724T<sup>0.0463</sup> exp(-3.61/T<sub>4</sub>))/(1.87 × 10<sup>8</sup> + 5.09 × 10<sup>10</sup>T<sup>-0.527</sup>)</td></tr><tr><td>C<sup>+</sup> + Si <sup>f</sup></td><td>1/(1.87 × 10<sup>8</sup> + 5.09 × 10<sup>10</sup>T<sup>-0.527</sup>)</td></tr><tr><td>Mg + Si<sup>+</sup> <sup>i</sup></td><td>2.9 × 10<sup>-9</sup></td></tr><tr><td>Mg<sup>+</sup> + Si <sup>i</sup></td><td>9.8 × 10<sup>-10</sup>(T/300)<sup>-0.0264</sup> exp(-0.59/T<sub>4</sub>)</td></tr><tr><td>Fe + Si<sup>+</sup> <sup>i</sup></td><td>1.9 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup> + Si <sup>i</sup></td><td>1.4 × 10<sup>-9</sup>(T/300)<sup>-0.236</sup> exp(-0.299/T<sub>4</sub>)</td></tr><tr><td>Na + Mg<sup>+</sup> <sup>i</sup></td><td>1.0 × 10<sup>-11</sup></td></tr><tr><td>Na<sup>+</sup> + Mg <sup>i</sup></td><td>3.76 × 10<sup>-11</sup>(T/300)<sup>-0.0302</sup> exp(-2.909/T<sub>4</sub>)</td></tr><tr><td>Mg + C<sup>+</sup> <sup>i</sup></td><td>1.1 × 10<sup>-9</sup></td></tr><tr><td>Mg<sup>+</sup> + C <sup>i</sup></td><td>3.01 × 10<sup>-10</sup>(T/300)<sup>-0.0832</sup> exp(-4.19/T<sub>4</sub>)</td></tr><tr><td>Mg + N<sup>+</sup> <sup>i</sup></td><td>1.2 × 10<sup>-9</sup></td></tr><tr><td>Mg<sup>+</sup> + N <sup>i</sup></td><td>9.48 × 10<sup>-10</sup>(T/300)<sup>0.139</sup> exp(-7.99/T<sub>4</sub>)</td></tr><tr><td>Na + Fe<sup>+</sup> <sup>i</sup></td><td>1.0 × 10<sup>-11</sup></td></tr></tbody></table>

<caption>**Table 4 continued**</caption>
<table><caption>**Table 4 (continued)**</caption><thead><tr><th>Reactants</th><th>Rate<sup>a</sup> (cm<sup>3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Na<sup>+</sup> + Fe<sup>i</sup></td><td>2.4 × 10<sup>-11</sup>(T/300)<sup>-0.0987</sup> exp(-3.21/T<sub>4</sub>)</td></tr><tr><td>Fe + C<sup>+</sup><sup>i</sup></td><td>2.6 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup> + C<sup>i</sup></td><td>1.12 × 10<sup>-9</sup>(T/300)<sup>0.0147</sup> exp(-3.90/T<sub>4</sub>)</td></tr><tr><td>Fe + O<sup>+</sup><sup>d</sup></td><td>1.71 × 10<sup>-9</sup></td></tr><tr><td>Fe<sup>+</sup> + O<sup>i</sup></td><td>5.0 × 10<sup>-10</sup>(T/300)<sup>-0.0337</sup> exp(-6.63/T<sub>4</sub>)</td></tr><tr><td>Mg + S<sup>+</sup><sup>i</sup></td><td>2.8 × 10<sup>-10</sup></td></tr><tr><td>Mg<sup>+</sup> + S<sup>i</sup></td><td>5.35 × 10<sup>-11</sup>(T/300)<sup>0.121</sup> exp(-3.15/T<sub>4</sub>)</td></tr><tr><td>Na + S<sup>+</sup><sup>i</sup></td><td>2.6 × 10<sup>-10</sup></td></tr><tr><td>Na<sup>+</sup> + S<sup>i</sup></td><td>1.86 × 10<sup>-10</sup>(T/300)<sup>0.151</sup> exp(-6.06/T<sub>4</sub>)</td></tr><tr><td>Fe + S<sup>+</sup><sup>i</sup></td><td>1.8 × 10<sup>-10</sup></td></tr><tr><td>Fe<sup>+</sup> + S<sup>i</sup></td><td>5.38 × 10<sup>-11</sup>(T/300)<sup>0.052</sup> exp(-2.85/T<sub>4</sub>)</td></tr><tr><td>S + H<sup>+</sup><sup>j</sup></td><td>exp[-50 + 13.3 ln T - 2.77(ln T)<sup>2</sup> + 0.243(ln T)<sup>3</sup> - 7.24 × 10<sup>-3</sup>(ln T)<sup>4</sup>]</td></tr><tr><td>S<sup>+</sup> + H<sup>j</sup></td><td>exp[-50.14 + 13.3 ln T - 2.77(ln T)<sup>2</sup> + 0.243(ln T)<sup>3</sup> - 7.24 × 10<sup>-3</sup>(ln T)<sup>4</sup> - 3.76/T<sub>4</sub>]</td></tr><tr><td>S + C<sup>+</sup><sup>k</sup></td><td>3.26 × 10<sup>-11</sup>(T/300)<sup>0.289</sup> exp(-338/T)</td></tr><tr><td>S<sup>+</sup> + C<sup>k</sup></td><td>6.97 × 10<sup>-11</sup>(T/300)<sup>0.356</sup> exp(-1.06/T<sub>4</sub>)</td></tr><tr><td>Fe + N<sup>+</sup><sup>d</sup></td><td>9.46 × 10<sup>-10</sup></td></tr><tr><td>Fe<sup>+</sup> + N<sup>d</sup></td><td>1.17 × 10<sup>-9</sup>(T/300)<sup>0.07</sup> exp(-7.696/T<sub>4</sub>)</td></tr><tr><td>O + He<sup>+</sup><sup>l</sup></td><td>4.99 × 10<sup>-15</sup>T<sub>4</sub><sup>0.379</sup> + 2.78 × 10<sup>-15</sup>T<sub>4</sub><sup>-0.216</sup> exp(T<sub>4</sub>/81.97)</td></tr><tr><td>O<sup>+</sup> + He<sup>l</sup></td><td>3.2 × (5.0 × 10<sup>-15</sup>T<sub>4</sub><sup>0.38</sup> + 2.78 × 10<sup>-15</sup>T<sub>4</sub><sup>-0.22</sup> exp(T<sub>4</sub>/81.97))T<sup>0.0377</sup> exp(-12.7/T<sub>4</sub>)</td></tr><tr><td>Na + Ca<sup>+</sup><sup>m</sup></td><td>3.0 × 10<sup>-9</sup></td></tr><tr><td>Na<sup>+</sup> + Ca<sup>m</sup></td><td>1.7 × 10<sup>-8</sup>(T/300)<sup>-0.168</sup> exp(-1.13/T<sub>4</sub>)</td></tr><tr><td>N + H<sup>+</sup><sup>n</sup></td><td>exp[-35.4 + 1.94 ln T - 0.154(ln T)<sup>2</sup> - 6.3 × 10<sup>-3</sup>(ln T)<sup>3</sup> - 1.16 × 10<sup>-3</sup>(ln T)<sup>4</sup>]</td></tr><tr><td>N<sup>+</sup> + H<sup>n</sup></td><td>exp[-40.1 + 6.4 ln T - 1.75(ln T)<sup>2</sup> + 0.18(ln T)<sup>3</sup> - 5.96 × 10<sup>-3</sup>(ln T)<sup>4</sup>]</td></tr><tr><td>Na + H<sup>+</sup><sup>o,p</sup></td><td>exp[48.2 - 43.0 ln T + 8.74(ln T)<sup>2</sup> - 0.77(ln T)<sup>3</sup> - 0.0256(ln T)<sup>4</sup>]</td></tr><tr><td>Na<sup>+</sup> + H<sup>o</sup></td><td>exp[46.6 - 41.6 ln T + 8.02(ln T)<sup>2</sup> - 0.65(ln T)<sup>3</sup> + 0.0197 × 10<sup>-3</sup>(ln T)<sup>4</sup> - 9.817/T<sub>4</sub>]</td></tr><tr><td>K + Na<sup>+</sup><sup>q</sup></td><td>1.0 × 10<sup>-11</sup></td></tr><tr><td>K<sup>+</sup> + Na<sup>q</sup></td><td>6.58 × 10<sup>-12</sup>(T/300)<sup>-0.203</sup> exp(-0.927/T<sub>4</sub>)</td></tr><tr><td>K + H<sup>+</sup><sup>p</sup></td><td>exp[-27.8 + 0.125 ln T + 0.0663(ln T)<sup>2</sup> - 0.0237(ln T)<sup>3</sup> - 1.36 × 10<sup>-3</sup>(ln T)<sup>4</sup>]</td></tr></tbody></table>

NOTE—$^aT_4 = T/10^4 K$, $^bGlover & Jappsen (2007)$, $^cKingdon & Ferland (1996)$, $^dRutherford & Vroom (1972)$, $^eNeufeld & Dalgarno (1987)$, $^fSatta et al. (2013)$, $^gStancil et al. (1999)$, $^hStancil et al. (1998)$, $^iPrasad & Huntress (1980)$, $^jZhao et al. (2005)$, $^kChenel et al. (2010)$, $^lZhao et al. (2004)$, $^mKwolek et al. (2019)$, $^nLin et al. (2005)$, $^oDutta et al. (2001)$, $^pWatanabe et al. (2002)$, $^qLavvas et al. (2014)$.