Survival of Terrestrial N₂–O₂ Atmospheres in Violent XUV Environments through Efficient Atomic Line Radiative Cooling

We integrate the time-dependent equation of continuity,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{j}}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}{n}_{j}(u+{u}_{j})\right]={S}_{j},\end{eqnarray} \tag{ 1 }$$

where t is the time, r is the radial distance measured from the planet's center, u is the bulk velocity, and n_j , u_j , and S_j are the number density, diffusive velocity, and source/sink term (including thermochemical and photochemical reactions) of species j, respectively. The energy equation is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}={Q}^{\mathrm{chem}}+{Q}^{\mathrm{uv}}+{Q}^{\mathrm{pe}}+{Q}^{\mathrm{diff}}+{Q}^{\mathrm{rad}}-{Q}^{\mathrm{hy}},\end{eqnarray} \tag{ 2 }$$

where e is the energy density, Q^chem is the heating rate through chemical reactions, Q^uv is that of photodissociation, Q^pe is that of collision with photoelectrons, Q^diff is that of thermal diffusion, Q^rad is the total rate of molecular radiation and absorption, and atomic line cooling, and Q^hy is the rate of hydrodynamic cooling. Note that no kinetic energy of the bulk flow is included because a steady-state motion of the atmospheric gas is assumed in this study.

We consider a steady-state motion of an inviscid, ideal gas. The equation of motion (or the Euler equation) is written in terms of mass density ρ = ∑_i m_i n_i as

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{dr}}=-\displaystyle \frac{1}{T}\displaystyle \frac{{dT}}{{dr}}+\displaystyle \frac{1}{\bar{m}}\displaystyle \frac{d\bar{m}}{{dr}}-\displaystyle \frac{g}{{u}_{0}^{2}}-\displaystyle \frac{u}{{u}_{0}}\displaystyle \frac{{du}}{{dr}},\end{eqnarray} \tag{ 3 }$$

with the ideal equation of state,

$$\begin{eqnarray}&&P=\displaystyle \frac{{k}_{{\rm{B}}}T}{\bar{m}}\rho ={u}_{0}^{2}\rho ,\end{eqnarray} \tag{ 4 }$$

where P is the pressure, T is the temperature, $\bar{m}$ is the mean mass of a gas particle, k_B is the Boltzmann constant, and ${u}_{0}(=\sqrt{{k}_{{\rm{B}}}T/\bar{m}})$ is the thermal velocity; g is the gravitational acceleration, which is given by the planet's mass M_p and the gravitational constant G as

$$\begin{eqnarray}&&g=\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{r}^{2}}.\end{eqnarray} \tag{ 5 }$$

The temperature is calculated from the energy density and specific heat C_v,j as

$$\begin{eqnarray}&&T=\displaystyle \frac{e}{{\sum }_{j}{n}_{j}{C}_{{\rm{v}},j}}.\end{eqnarray} \tag{ 6 }$$

The majority of the values of C_v,j are taken from NASA CEA (Gordon & McBride 1994). For species unavailable in NASA CEA, we adopt the specific heats for ideal gases with the internal degrees of freedom taken as three for atoms, five for diatomic molecules, and seven for all other molecules.

By use of the steady-state equation of continuity (i.e., r² ρ u = const.), Equation (3) can be written in terms of the bulk velocity uu as

$$\begin{eqnarray}&&\left(1-\displaystyle \frac{{u}^{2}}{{u}_{0}^{2}}\right)\displaystyle \frac{1}{u}\displaystyle \frac{{du}}{{dr}}=\displaystyle \frac{1}{T}\displaystyle \frac{{dT}}{{dr}}-\displaystyle \frac{1}{\bar{m}}\displaystyle \frac{d\bar{m}}{{dr}}+\displaystyle \frac{g}{{u}_{0}^{2}}-\displaystyle \frac{2}{r},\end{eqnarray} \tag{ 7 }$$

which is referred to as the wind equation (e.g., Parker 1958). We integrate Equation (7) inward from the upper boundary, which we locate at the exobase. We adopt the Jeans effusion velocity ${ \mathcal U }$ as the outward velocity at the exobase (e.g., Chamberlain 1963):

$$\begin{eqnarray}&&{{ \mathcal U }}_{j}={u}_{0,j}\displaystyle \frac{(1+{X}_{j})\exp (-{X}_{j})}{2\sqrt{\pi }},\end{eqnarray} \tag{ 8 }$$

where u_0,j and X_j are the thermal velocity and Jeans escape parameter of gas species j at the exobase, respectively; X_j is defined by

$$\begin{eqnarray}&&{X}_{j}=\displaystyle \frac{{r}_{\mathrm{exo}}{g}_{\mathrm{exo}}{m}_{j}}{{k}_{{\rm{B}}}{T}_{\mathrm{exo}}},\end{eqnarray} \tag{ 9 }$$

where the subscript "exo" means the exobase. Then, the bulk velocity is estimated by density-weighted Jeans velocity:

$$\begin{eqnarray}&&u=\displaystyle \frac{1}{\rho }\sum _{j}{n}_{j}{m}_{j}{{ \mathcal U }}_{j}.\end{eqnarray} \tag{ 10 }$$

At the upper boundary of the exobase, the bulk velocity is replaced with the Jeans effusion velocities of the individual species in Equation (1).

We consider both thermo- and photochemical reactions. We adopt the chemical network provided in Johnstone et al. (2018), which is validated to reproduce the properties of the Earth's upper atmosphere and has been used in several previously published atmospheric models. The chemical network, however, does not include some reactions that are expected to play a crucial role in high-XUV environments of interest in this study.

The total heating rate due to chemical reactions, Q^chem, is calculated as

$$\begin{eqnarray}&&{Q}^{\mathrm{chem}}=\sum _{{\ell }}{E}_{{\ell }}{R}_{{\ell }},\end{eqnarray} \tag{ 11 }$$

where E_ℓ is the enthalpy of reaction for the ℓth reaction and R_ℓ is the reaction rate, which is calculated as

$$\begin{eqnarray}&&{R}_{{\ell }}={k}_{{\ell }}\prod _{j}{n}_{j};\end{eqnarray} \tag{ 12 }$$

the RHS means the reaction coefficient, k_ℓ , times the product of the number densities of all of the reactants. We calculate the source/sink term of Equation (1) by summing all reactions, including the chemical and photochemical reactions.

2.2.1. Thermochemistry

We consider thermochemical reactions between 47 chemical species composed of the four elements, H, C, O, and N, including ions and electrons. We adopt all of the reactions and species from Johnstone et al. (2018). A major difference from Johnstone et al. (2018) and other previous models is that we also include chemical reactions involving internal excitation and ionization. The former is detailed in Section 2.2.1.1, while for the latter, we take into account ionization via collision with electrons, following Voronov (1997), who presented fitting formulae for the reaction coefficients. In addition, unlike previous models, we consider the effects of endothermic reactions (Section 2.2.1.2) and recombination and emissive de-excitation (Section 2.2.1.3).

2.2.1.1. Excitation and De-excitation

We consider the excitation and de-excitation of N, O, and O⁺ via collisions with electrons, following Tayal (2006) and Ito & Ikoma (2021). Those rate coefficients are estimated with the effective collision strength, γ. As for excitation, the rate coefficient, k_lu , is given by (e.g., Tielens 2005)

$$\begin{eqnarray}&&{k}_{{lu}}=\gamma \displaystyle \frac{8.629\times {10}^{-6}}{{s}_{l}\sqrt{T}}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{E}_{{lu}}}{{k}_{{\rm{B}}}T}\right),\end{eqnarray} \tag{ 13 }$$

where the subscripts l and u denote the lower and upper states, respectively, s_l is the statistical weight of state l, and ΔE_lu is the energy difference between the upper and lower states. As for de-excitation, the rate coefficient is given by

$$\begin{eqnarray}&&{k}_{{ul}}=\displaystyle \frac{{s}_{l}}{{s}_{u}}{k}_{{lu}}\exp \left(\displaystyle \frac{{\rm{\Delta }}{E}_{{lu}}}{{k}_{{\rm{B}}}T}\right).\end{eqnarray} \tag{ 14 }$$

We ignore the small dependence of γ on T and adopt the value of γ at T = 1 × 10⁴ K, which is typical for highly irradiated upper atmospheres, following Ito & Ikoma (2021). Also, we include transitions between states with different electronic configurations, namely, ⁴S-²D-²P for N, ³P-¹D-¹S for O, and ³P-¹D-¹S for O⁺. The detailed level parameters of the statistical weight and energy and transition parameters are summarized in Appendix C.

2.2.1.2. Endothermic Reactions

We consider the reverse reactions for the reactions in the chemical network of Johnstone et al. (2018). Unlike in the atmospheres of the present-day solar system planets, the upper atmosphere of interest here is sufficiently hot that endothermic reactions are expected to occur efficiently. The rate coefficients of reverse reactions, k_r , are derived according to the principle of microscopic reversibility, which is discussed in several studies focusing on the atmospheric composition of hot Jupiters (e.g., Visscher & Moses 2011; Heng et al. 2016), and are given by

$$\begin{eqnarray}&&{k}_{r}={k}_{f}\exp \left(-\displaystyle \frac{{\rm{\Delta }}G}{{k}_{{\rm{B}}}T}\right){\left(\displaystyle \frac{{k}_{{\rm{B}}}T}{{P}_{0}}\right)}^{{q}_{p}-{q}_{r}},\end{eqnarray} \tag{ 15 }$$

where k_f is the rate coefficient of the forward reaction, ΔG is the Gibbs free energy change for the forward reaction at the standard pressure (P₀ = 1 × 10⁶ dyn/cm²), and q_p and q_r are the stoichiometric coefficients of products and reactants in the forward reaction, respectively. We calculate ΔG for a given temperature from the enthalpy and entropy taken primarily from NASA CEA. For the values unavailable in NASA CEA, we adopt the values for ideal gases with the internal degrees of freedom taken as three for atoms, five for diatomic molecules, and seven for other molecules. For their specific entropy at the standard state, we use the values calculated by the Gaussian-2 composite method, which are available in the NIST Computational Chemistry Comparison and Benchmark Database (Johnson 1999). We calculate the enthalpy for excited states simply by calculating the sum of the excitation energy and ground-state enthalpy, for which no measurements are available.

2.2.1.3. Recombination and Emissive De-excitation

Our chemical network also includes radiative recombination and de-excitation via spontaneous emission. We assume that the radiation emitted upon recombination directly escapes to space without extinction by the atmospheric gas.

As for the spontaneous emission, we approximate the effects of radiative transfer by estimating the frequency-integrated probability for photons to escape from the atmosphere (known as the escape probability method; Irons 1978), instead of using a detailed treatment of radiative transfer. The escape probability P_e is given as a function of optical thickness τ (Hollenbach & McKee 1979; Kwan & Krolik 1981):

$$\begin{eqnarray}{P}_{e}=\left\{\begin{array}{ll}\sqrt{\pi }\tau {\left({a}_{e}+\displaystyle \frac{\sqrt{\mathrm{ln}\tau }}{1+{b}_{e}\tau }\right)}^{-1} & (\tau \geqslant 1),\\ \displaystyle \frac{1-\exp (-2\tau )}{2\tau } & (\tau \lt 1),\end{array}\right.\end{eqnarray} \tag{ 16 }$$

where a_e = 1.2 and b_e = 1.0 × 10⁻⁵. The upper atmosphere of interest here is sufficiently tenuous that the line profile is determined by the Doppler (thermal) broadening. The frequency-integrated optical thickness is thus given by $\tau =\sqrt{\pi }{\tau }_{0}$, with the optical thickness at the line center τ₀, which we estimate using the Einstein A coefficients. In this method, it is assumed that the line profile of radiation, which is determined by the temperature of the source gas, never varies during propagation.

We assume that the radiation propagates equally upward and downward (i.e., a two-stream approximation). The total escape probability P_tot is, thus, given by

$$\begin{eqnarray}&&{P}_{\mathrm{tot}}=\displaystyle \frac{{P}_{e}(\tau )+{P}_{e}({\tau }_{\mathrm{tot}}-\tau )}{2},\end{eqnarray} \tag{ 17 }$$

where τ and τ_tot are the optical thicknesses at the emission altitude and the lower boundary, respectively. We regard the photons absorbed below the lower boundary to be escaped. This is because those absorbed photons are quickly partitioned into the internal energy of the atmospheric gas and are, then, re-radiated by a blackbody emission, since the lower atmosphere is in the local thermal equilibrium. Thus, we calculate the transition rate via spontaneous emission with the Einstein A coefficient multiplied by the total escape probability.

2.2.2. Photochemistry

We use the photochemical network presented in Johnstone et al. (2018). We exclude the reactions numbered 446, 447, 454, 455, 459, and 468 in the appendix of Johnstone et al. (2018), for which the absorption cross sections or quantum yields are unavailable in their references. For all of the reactions, we take the wavelength-dependent absorption cross sections from the PHIDRATES database (Huebner & Mukherjee 2015). Our model covers the wavelength range between 0.5 nm and 400 nm to allow us to resolve the oxygen chemistry.

The reaction rate coefficient of the ℓth photoreaction is given by

$$\begin{eqnarray}&&{k}_{{\ell }}={\int }_{{\lambda }_{{\ell }}^{\mathrm{cr}}}^{\infty }{\sigma }_{{\ell }}(\lambda ){I}_{\lambda }(r)\,{\rm{d}}\lambda ,\end{eqnarray} \tag{ 18 }$$

where λ is the wavelength, ${\lambda }_{{\ell }}^{\mathrm{cr}}$ is the threshold wavelength for the ℓth reaction, σ_ℓ (λ) is the absorption coefficient for a given wavelength, and I_λ (r) is the irradience per photon, quadratic centimeter, wavelength, and second for a given wavelength at altitude r. I_λ (r) is calculated by the radiative transfer equation:

$$\begin{eqnarray}&&{I}_{\lambda }(r)={I}_{\lambda }(\infty )\exp \left\{-\displaystyle \frac{1}{\mu }\sum _{{\ell }}\left({\int }_{r}^{\infty }{\sigma }_{{\ell }}(\lambda ){n}_{{\ell }}\,{\rm{d}}r\right)\right\},\end{eqnarray} \tag{ 19 }$$

where I_λ ( ∞ ) is the irradiance at the top of the atmosphere, μ is the cosine of the stellar zenith angle, and n_ℓ is the number density of the absorbed species for the ℓth reaction. As in Johnstone et al. (2018), we assume $\mu =\cos (66^\circ )$ to evaluate the global averaged structure of the upper atmosphere.

The absorbed photon energy is largely consumed by the chemical reactions. For photodissociative reactions, the remaining energy is converted into the kinetic energy of molecules and atoms, which is subsequently dissipated as heat through collisional relaxation. We assume that these conversions occur instantaneously. Thus, the heating rate due to photodissociation is given by

$$\begin{eqnarray}&&{Q}^{\mathrm{uv}}=\sum _{{\ell }}{\int }_{0}^{{\lambda }_{{\ell }}^{\mathrm{cr}}}\left(\displaystyle \frac{{hc}}{\lambda }-\displaystyle \frac{{hc}}{{\lambda }_{{\ell }}^{\mathrm{cr}}}\right){n}_{{\ell }}{\sigma }_{{\ell }}(\lambda ){I}_{\lambda }(r)\,{\rm{d}}\lambda ,\end{eqnarray} \tag{ 20 }$$

where h is the Planck constant, and c is the speed of light. The terms with (hc/λ) and $({hc}/{\lambda }_{{\ell }}^{\mathrm{cr}})$ represent the absorbed photon energy for a given wavelength and the energy required for the ℓth reaction to occur, respectively.

For photoionization reactions, the remaining energy is consumed in a complex manner. The energy is first transformed into the kinetic energy of electrons. Those electrons, which are termed photoelectrons, have a higher energy than the ambient thermal electrons. Due to their high energy, collisions with photoelectrons result in chemical reactions of neutral species, including secondary photoionization, and heating of thermal electrons. For such photoelectron-induced processes, we assume a local approximation such that the produced photoelectrons are consumed locally, as assumed in Johnstone et al. (2018). The outline of the calculation method is as follows (see Johnstone et al. 2018, for the details): First we obtain the initial energy distribution of photoelectrons, assuming the energy of each photoelectron as the remaining energy from each photoionization process. Then, the degradation of the energies of the photoelectrons is calculated by inelastic collisions between photoelectrons and neutral species, including chemical reactions and excitation, from higher energy to lower energy. Finally, we estimate the heating rate for the photoelectrons, Q^pe, using the energy distribution calculated above in the expression of energy transfer from photoelectrons to ambient thermal electrons given by Schunk & Nagy (1978). Here we consider the inelastic collision of the photoelectrons with H, C, O, N, O₂, N₂, CO, and CO₂. The cross sections, as a function of the photoelectron energy, are taken from Voronov (1997) for H, Suno & Kato (2006) for C, Jackman et al. (1977) for O, Tian et al. (2008b) for N, Green & Sawada (1972) and Jackman et al. (1977) for O₂, Green & Barth (1965) and Green & Sawada (1972) for N₂, Sawada et al. (1972) and Jackman et al. (1977) for CO, and Jackman et al. (1977) and Bhardwaj & Jain (2009) for CO₂. Excitation states that are not included in our chemical model are assumed to be quenched immediately via radiative de-excitation.

2.3.1. Diffusion

We consider multicomponent molecular and eddy diffusion. The total diffusion velocity of species j is expressed as the sum:

$$\begin{eqnarray}&&{u}_{j}={u}_{j}^{\mathrm{mol}}+{u}_{j}^{\mathrm{eddy}},\end{eqnarray} \tag{ 21 }$$

where ${u}_{j}^{\mathrm{mol}}$ and ${u}_{j}^{\mathrm{eddy}}$ are the velocities caused by molecular and eddy diffusion, respectively. Although a minor-component approximation is often used in upper-atmospheric modeling (e.g., Banks & Kockarts 1973), we should use a diffusion model applicable to a wide range of conditions because the dominant species and temperature change with XUV intensity and altitude in the upper atmosphere.

For molecular diffusion, we adopt the formula derived by García Muñoz (2007a), which is also used in hydrodynamic simulations for hot Jupiters (García Muñoz 2007b) and for hot rocky super-Earths (Ito & Ikoma 2021). The formula is based on the momentum equations for a multicomponent gas derived by Burgers (1969), assuming no momentum transfer via diffusion, local gas neutrality (i.e., the ambipolar constraint), and no external electromagnetic forces. The diffusion model explicitly includes the collisions among all species and the resultant dragging effect. Although not repeating the explicit expression of the diffusion model here, we adopt the second-order approximation of the diffusion matrix and their binary diffusion coefficients for evaluating diffusion velocity (see García Muñoz 2007a, for the details). Empirical coefficients of binary diffusion are preferentially adopted. When no reliable data are available, however, the binary diffusion coefficients are estimated based on the hard-sphere model for neutral–neutral and neutral–electron pairs, the Coulomb interaction model for ion–ion and ion–electron pairs, and the interaction via induced polarization potential for neutral–ion pairs. For neutral–ion pairs, we adopt values of the polarizability of neutral species taken primarily from the CRC Handbook (William 2016) or the typical value of 1.0 × 10⁻²⁴ cm³ is used when no data is available. In addition, we include the binary diffusion coefficients for charge exchange given in Table 4.4 of Schunk & Nagy (2000).

For eddy diffusion, the diffusion speed is given by (e.g., Banks & Kockarts 1973)

$$\begin{eqnarray}&&{u}_{j}^{\mathrm{eddy}}=-{K}_{{\rm{E}}}\left(\displaystyle \frac{1}{{n}_{j}}\displaystyle \frac{\partial {n}_{j}}{\partial r}+\displaystyle \frac{1}{\bar{H}}+\displaystyle \frac{1}{T}\displaystyle \frac{\partial T}{\partial r}\right),\end{eqnarray} \tag{ 22 }$$

where K_E is the eddy diffusion coefficient, and $\bar{H}$ is the pressure scale-height. We adopt the conventional form of the eddy diffusion coefficient,

$$\begin{eqnarray}&&{K}_{{\rm{E}}}={A}_{{\rm{E}}}{N}^{{B}_{{\rm{E}}}},\end{eqnarray} \tag{ 23 }$$

where N is the particle number density of the entire gas, and A_E and B_E are empirical constants. We use A_E = 10⁸ and B_E = − 0.1, which are valid for the Earth's current atmosphere (Johnstone et al. 2018).

For the diffusion of electrons, we assume that the atmosphere maintains local electrical neutrality. Assuming no magnetic field, the electron diffusive flux is equal to the sum of the ion's diffusive fluxes (e.g., Shinagawa & Cravens 1989):

$$\begin{eqnarray}&&{u}_{{\rm{e}}}{n}_{{\rm{e}}}=\sum _{j}{u}_{j}{n}_{j}.\end{eqnarray} \tag{ 24 }$$

This condition is also applied to the Jeans effusion velocity of electrons.

2.3.2. Conduction

Heat transport occurs via thermal conduction. The rate of heating via conduction due to both the molecular and eddy diffusion (or molecular and eddy conduction) is given by (Banks & Kockarts 1973; Hunten 1974)

$$\begin{eqnarray}&&{Q}^{\mathrm{diff}}=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}{\kappa }_{\mathrm{mol}}\displaystyle \frac{\partial T}{\partial r}+{r}^{2}{\kappa }_{\mathrm{eddy}}\left(\displaystyle \frac{\partial T}{\partial r}+\displaystyle \frac{g}{{C}_{P}}\right)\right],\end{eqnarray} \tag{ 25 }$$

where κ_mol and κ_eddy are the thermal conductivity due to molecular and eddy diffusion, respectively, and C_P is the specific heat at constant pressure, which is derived in the same manner as C_v.

The molecular conduction occurs through neutral particles, ions, and electrons. We adopt a conduction model for each species almost identical to that used in Johnstone et al. (2018), which is outlined as follows: For neutral particles, we consider N₂, O₂, CO₂, CO, O, and H and take their conductivities from Schunk & Nagy (2000). We use the expression of the total conductivity given by Banks & Kockarts (1973) as the mixture of neutral species. For ions and electrons, we adopt the models given in Banks & Kockarts (1973). We use the number density-weighted average value of the ion conductivity. The electron conductivity includes a reduction in conductivity caused by collisions with neutral species for which we only consider the effects of N₂, O₂, O, and H.

The eddy conductivity is related to the eddy diffusion coefficient as κ_eddy = ρ C_P K_E (Hunten 1974). As found in Equation (25), the eddy conduction includes the convective turbulence term expressed as the dry adiabatic lapse rate, g/C_P .

We consider the radiative cooling via atomic electronic transitions and molecular vibrational-rotational transitions. In addition, we consider heating by the absorption of stellar near-IR radiation by molecular species such as CO₂ and H₂O. For both of the atomic and molecular transitions, we assume the statistical equilibria (e.g., Manuel et al. 2001). The populations (or number densities) of respective excited states in the atomic and molecular species, hereafter simply termed the level populations, are determined via collisional and radiative excitation/de-excitation. In the statistical equilibrium, the level populations achieve a steady state such that

$$\begin{eqnarray}&&{\chi }_{{ij}}{n}_{j}^{(x)}=0\end{eqnarray} \tag{ 26 }$$

where χ_ij is the matrix of excitation/de-excitation rate coefficients and ${n}_{j}^{(x)}$ is the population of excited state j in species x. Combining the processes mentioned above together, we can express χ_ij as

$$\begin{eqnarray}&&{\chi }_{{ij}}={C}_{{\rm{e}},{ij}}-{C}_{{\rm{d}},{ij}}+{R}_{{\rm{e}},{ij}}-{\left({{AP}}_{\mathrm{tot}}\right)}_{{ij}},\end{eqnarray} \tag{ 27 }$$

where C_e,j and C_d,j are the matrices of collisional excitation and de-excitation coefficients, respectively, R_e,j is that of radiative excitation rate coefficients, A represents the Einstein coefficients, and P_tot is the escape probability. We solve Equation (26) with mass conservation and estimate the radiative cooling/heating rate Q^rad as

$$\begin{eqnarray}&&{Q}^{\mathrm{rad}}=\sum _{j}{n}_{j}\sum _{i}{E}_{{ij}}\left\{{R}_{{\rm{e}},{ij}}-{({{AP}}_{\mathrm{tot}})}_{{ij}}\right\},\end{eqnarray} \tag{ 28 }$$

where E_ij is the energy difference between the i and j states. Below we describe the detailed treatments and assumptions for evaluating the transitions of atomic and molecular species.

2.4.1. Atomic Electronic Transitions

2.4.2. Molecular Vibrational-rotational Transition

Hydrodynamic cooling, which is the combined effects of work done by pressure, internal advective transport plus kinetic energy, and conversion to the gravitational potential, is known to dominate the cooling process in the upper parts of escaping atmospheres of highly irradiated planets (e.g., Yelle 2004; Tian et al. 2008a). Instead of the approximate formula used in Tian et al. (2008a), who only considered adiabatic cooling, we adopt the standard energy equation of an inviscid fluid under the influence of planetary gravity (e.g., García Muñoz 2007a):

$$\begin{eqnarray}&&{Q}^{\mathrm{hy}}=\displaystyle \frac{\partial \left[\left(e+\tfrac{1}{2}\rho {u}^{2}+P\right){{ur}}^{2}\right]}{{r}^{2}\partial r}+\rho {ug}.\end{eqnarray} \tag{ 29 }$$

We neglect the stellar tide and the centrifugal force due to the orbital motion because the planets of interest are so far from their central star, and their exospheric radius is so small relative to their Hill radius, that those effects are negligible.

We integrate the thermal and compositional structure of the upper atmosphere for a given XUV irradiation spectrum. The atmosphere is divided into spherical layers. The lowermost layer is 2 km thick, and the layer thickness increases with altitude in such a way that the thickness ratio of two neighboring layers is 1.02. Physical quantities for each layer are defined at the midpoint of the layer. To determine the outer edge of the atmosphere (i.e., the exobase), we add layers upwards until the Knudsen number Kn, which is the ratio of the mean free path length of gas particles to the local scale-height, approaches unity (exactly, 0.7 ≤ Kn ≤ 1.2).

We find steady-state solutions of the upper-atmospheric structure, integrating from the lower boundary to the exobase. We adopt the present-day Earth's lower atmosphere as the lower boundary condition; we use the values of the number density of O, O₂, and N₂ and the temperature at the altitude of 75 km taken from the empirical NRLMSISE-00 model (Picone et al. 2002) on 1990 January 1. In addition, we assume Earth-like CO₂ and H₂O mixing ratios of 4.0 × 10⁻⁴ and 6.0 × 10⁻⁶, respectively. The detailed lower boundary conditions used in our model are given in Table 1.

We adopt an implicit solver of DLSODE with the backward differential formula (Hindmarsh 1983) for the time integration of Equations (1) and (2). The solver is suitable for stiff ordinary-differential-equation systems such as the chemical network (e.g., Grassi et al. 2014). We adopt 10⁻⁴ and 10⁻⁷, respectively, as the values of the relative and absolute tolerances for the solver. From the lower boundary, where mean molecular weight, temperature, and wind profiles are known, we radially integrate Equation (3). Allowing the steady-state density structure, we recalculate the number density of each species at each location for every time step, keeping the mixing ratio derived from the continuous equation. For the initial condition, we assume that the atmosphere is perfectly mixed and the temperature and mixing ratios of all layers are the same as those at the lower boundary. For the criterion of convergence, we introduce the variable defined as

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Theta }}=\displaystyle \frac{1}{{\rm{\Theta }}}\displaystyle \frac{d{\rm{\Theta }}}{{dt}},\end{eqnarray} \tag{ 30 }$$

where Θ is the physical quantity. When ΔΘ of all of the physical quantities become less than 10⁻⁶ in all of the layers, we judge the calculation to be converged.

For the present-day solar spectrum, we adopt the spectrum model that Claire et al. (2012) developed based on the ATLAS 1 measurements obtained near the solar maximum (Thuillier et al. 2004). To investigate the response of the atmospheric structure to an increase in the amount of XUV radiation, we regard the present-day solar spectrum as the reference and simply multiply the XUV intensity at each wavelength by N for which we consider values of 1 to 1000 in this study. This approach is similar to that in Tian et al. (2008a), who extrapolated the XUV spectrum to that in more active conditions from the observed solar minimum and maximum spectra. We note that the range of XUV wavelength in this method (≤105 nm) is different from the hydrogen ionization edge (≤91 nm), following the definition of Tian et al. (2008a). The actual shape of the XUV spectrum is not scaled by the amount of XUV radiation, and far-UV (FUV) and near-UV (NUV) spectra also depend on the specific stellar activity (e.g., Claire et al. 2012). To discuss the sensitivity to the UV spectrum, we also calculate the thermospheric structure using spectra emitted by young Sun-like stars (Claire et al. 2012). We investigate the impact of the UV spectrum in Section 3.4.

Figure 1 shows the vertical temperature profiles in the upper atmosphere for three different levels of XUV irradiation (F_XUV), one, three, and five times the present-day Earth's irradiation (F_XUV,⊕), which are indicated with blue, green, and red solid lines, respectively. For comparison, the profiles obtained without atomic line cooling are also shown by dashed lines. In every case, temperature is found to increase with increasing altitude because of heating at high altitudes (see below for details). In the cases of 1F_XUV,⊕ and 3F_XUV,⊕, the solid and dashed lines overlap each other, indicating that the atomic line cooling is ineffective. By contrast, in the case of 5F_XUV,⊕, an obvious difference is found between the solid and dashed lines; namely, the atomic line cooling results in significant temperature. Although not shown, the difference increases with increased XUV irradiation.

Zoom In Zoom Out Reset image size

To understand the effects of atomic line cooling, we show the profiles of the energy budget for F_XUV = 5 F_XUV,⊕ in Figure 2. The upper and lower panels show the results obtained with and without atomic line cooling, respectively. In both cases, incident stellar XUV (≲91 nm) is absorbed and, thereby, photoelectrons are produced at high altitudes (see the dark blue lines). Such photoelectrons collide with and excite the ambient atoms. When atomic line cooling is omitted (see Figure 2(b)), almost all of the energy thus absorbed is transferred downward by conduction (see the purple line), followed by radiative cooling via vibration of molecules such as NO and CO₂ at low altitudes (∼100 km; see the red line). The upper atmosphere is in radiative equilibrium, as previously understood (e.g., Johnstone et al. 2018).

Zoom In Zoom Out Reset image size

By contrast, in Figure 2(a), the atomic line cooling (brown line) produces a major contribution to the energy loss, in addition to conduction (purple line). Under this condition, the radiative emission from O(¹D) associated with excitation via collisions with N₂, O, and e⁻ mainly contributes to atomic line cooling. The dominant cooling atom and species that play a major role in collisional excitation depend on the composition and temperature of the atmospheric gas at each altitude. For instance, emission from excited N, N⁺, O, and O⁺ associated with collisional excitation by e⁻ is dominant under highly irradiated conditions. Although conduction also transfers the absorbed energy downward, similarly to the case without atomic line cooling, its contribution becomes smaller at higher altitudes, where the temperature is sufficient to excite electric transitions. As such, the local energy loss associated with the atomic line cooling results in a reduction in the atmospheric temperature, as illustrated in Figure 1. Until F_XUV = 5 F_XUV,⊕, the atoms are insufficiently excited, so that the energy budget is similar to that of the case without atomic line cooling.

Figure 3 shows the calculated temperature profiles for the XUV irradiation level, F_XUV, of up to 1000 F_XUV,⊕. For F_XUV < 300F_XUV,⊕ (from black to yellow lines), the temperature profiles are qualitatively similar to one another; the temperature rises monotonically with increasing altitude above ∼100 km. An increase in F_XUV results in an overall temperature rise, including the exospheric temperature. By contrast, for F_XUV ≥ 300F_XUV,⊕ (shown by reddish lines), the temperature profiles are not monotonic. Between the two stratified regions (i.e., between ∼170 and ∼350 km), the temperature falls with increasing altitude. Also, in the uppermost region (≳2000 km), temperature falls with increasing altitude. In this high-F_XUV regime, the exospheric temperature decreases, as the XUV irradiation increases. This qualitative difference in the temperature profile between the low-F_XUV and high-F_XUV regimes can be interpreted as follows.

Zoom In Zoom Out Reset image size

When F_XUV < 300F_XUV,⊕, as shown in Figure 2(a), the absorption of the incident stellar XUV (≲91 nm), which causes photoionization, takes place at all altitudes, whereas the stellar FUV and NUV (∼91–400 nm) are absorbed only below 100 km. The temperature above 100 km is thus controlled by heating via photoionization-driven chemical reactions caused by XUV absorption. The specific heating rate increases with increasing altitude (see Figure 2(a)), since the XUV energy density is lower at lower altitudes because of larger optical depths. This is why temperature rises monotonically with increasing altitude in the low-F_XUV regime.

When F_XUV ≥ 300F_XUV,⊕, the absorption of FUV and NUV also takes place above 100 km, in contrast to the low-F_XUV regime. This is confirmed in Figure 4(a), which shows the altitude profiles of the absorption rate for XUV (red) and FUV+NUV (blue). The contribution of FUV+NUV is much larger than that of XUV below 200 km. This is because neutral molecules, which absorb FUV+NUV, exist below ∼200 km. This is illustrated in Figure 4(b), where the compositional fractions of neutral molecules, neutral atoms, and ion species are shown. Unlike in the low-F_XUV regime, because of relatively high temperatures, and thereby large scale-heights, molecules exist above 100 km in the high-F_XUV regime. Two peaks of the absorption rate for FUV+NUV, found at ∼200 km and ∼100 km in Figure 4(a), are for N₂ and O₂, respectively. Since the FUV+NUV energy density decreases with decreasing altitude below this, the temperature peaks at ∼200 km.

Zoom In Zoom Out Reset image size

Temperature also peaks at ∼2000 km and then falls with increasing altitude for F_XUV ≥ 300F_XUV,⊕ (see Figure 3). At 3000 km, as illustrated in Figure 4(b), the ionized fraction is about 50. Since ions absorb UV poorly, the absorbed energy per unit mass increases with decreasing altitude. The optical depth for XUV is approximately unity at ∼2000 km, below which it increases rapidly with decreasing altitude.

Although thermal conduction and advective energy transport tend to relax such nonmonotonic temperature profiles, the strong dependence of atomic line cooling on temperature inhibits temperature gradients large enough for efficient energy transport. Indeed, at ∼100 km, where atomic line cooling is inefficient, while efficient absorption of long-wavelength radiation by O₂ occurs, eddy thermal diffusion dominates, resulting in no temperature peak. Thus, in systems where advection and thermal conduction work inefficiently, the atmospheric composition and the incident UV spectrum determine the profiles of the heating rate and temperature in the upper atmosphere.

Figure 5 shows the rates of energy deposition/loss integrated over the entire atmosphere (i.e., ∫4π r² Qdr) as functions of the XUV irradiation level. The energy deposition via XUV absorption (red line) increases with the increased XUV irradiation level; the dependence is less than linear because X-ray (< 10 nm) and FUV (≳100 nm) are incompletely absorbed in the atmosphere. The energy deposition rate via molecular absorption (green line) is independent of the XUV irradiation, because the abundances of H₂O and CO₂ are also insensitive to the latter.

Zoom In Zoom Out Reset image size

For energy loss, while molecular rotation-vibration cooling (blue) dominates for F_XUV < 10F_XUV,⊕, atomic line cooling (purple) and radiative recombination (orange) become more effective with increasing F_XUV and dominate for F_XUV > 10F_XUV,⊕. The radiative recombination rate depends more strongly on the XUV irradiation level than on the atomic line cooling, and, consequently, the former dominates over the latter for F_XUV ≥ 30F_XUV,⊕. This can be explained as follows: First, the atomic line cooling rate is proportional to the number density of electrons (produced via photoionization of atoms). This is because under non-LTE conditions, where collisional excitation (mainly with electrons) occurs less frequently than radiative de-excitation, the former controls the atomic line cooling. Furthermore, the radiative de-excitation rate is largely insensitive to the gas species, including ionized states. Meanwhile, the radiative recombination rate is proportional to the square of the number density of electrons, because, by definition, recombination requires a collision between an electron and an ion. Since more electrons are produced at higher XUV levels, the radiative recombination is more sensitive to F_XUV. As such, the temperature rises, and thus expansion of the thermosphere is suppressed under high-XUV conditions.

Finally we make a few comments on the above results. First, although radiative recombination dominates over atomic cooling at high-XUV irradiation levels, this does not mean that the latter is less important. It is because of the suppression of temperature rise via atomic cooling that the radiative recombination becomes effective. Second, for F_XUV = 5F_XUV,⊕, although the atomic line cooling yields an energy-loss rate smaller by a factor of five than the molecular cooling, the former has a significant effect on the thermal structure of the upper atmosphere (see Figure 1). As described in Section 3.1, this is because the local energy balance between XUV absorption and atomic line cooling hinders a rise in temperature in the upper atmosphere. Finally, the hydrodynamic cooling always makes a limited contribution to the energy loss, as discussed in the following subsection.

Figure 6 shows the exobase temperature T_exo as a function of the XUV irradiation level F_XUV with the effects of atomic line cooling (red solid line). For comparison, we also show the solutions without the effects of atomic line cooling (red dashed line). First, without atomic line cooling, the exobase temperature rises rapidly with the increased XUV irradiation level for F_XUV ≤ 5F_XUV,⊕. Beyond this level, although not shown here, an increase in XUV irradiation results in a relatively gradual decrease in exobase temperature, as found in Tian et al. (2008a; see their Figure 2). In this regime (called the hydrodynamic regime in Tian et al. 2008a), high temperature leads to large Jeans effusion velocity. Thus, such a decrease in exobase temperature is a hydrodynamic effect (i.e., the effect of advective cooling). When advection is ignored (i.e., completely hydrostatic equilibrium being assumed), no solution is found for F_XUV > 5F_XUV,⊕ (see Tian et al. 2008a), which is often called the atmospheric blow-off.

Zoom In Zoom Out Reset image size

When the atomic line cooling is incorporated (solid red line), the behavior of the exobase temperature is in clear contrast to that obtained without atomic line cooling. The exobase temperature increases gradually with F_XUV, then levels off at F_XUV ≃ 100F_XUV,⊕, and, thereafter, decreases slightly with increasing F_XUV for F_XUV > 300F_XUV,⊕. Of particular interest is that the atmosphere does not enter the hydrodynamic regime until F_XUV = 1000 F_XUV,⊕. As F_XUV increases from 1 F_XUV,⊕ to 1000 F_XUV,⊕, the number density-weighted escape parameter (see Equation (9)) decreases from 106 to 8.6, indicating that advection becomes increasingly significant. Nevertheless, the value is still larger than the blow-off criterion (∼2.5; Watson et al. 1981). The decrease in exobase temperature found for F_XUV > 300F_XUV,⊕ is not due to advective cooling, but due to an increase in the degree of ionization.

Young Sun-like stars emit much larger amounts of XUV radiation than the present Sun and become less active with age. In addition to the total amount of radiation, the emission intensity at each wavelength (or emission spectrum) varies with age. Thus, taking the age change in stellar emission spectrum (Claire et al. 2012) into account, we calculate and show the exobase temperature as a function of the XUV irradiation level, shown by the blue line in Figure 6. As found by comparison between the blue and red solid lines, the age change in stellar spectrum has a limited effect on the T_exo–F_XUV relationship.

Figure 7 shows temperature profiles for solar irradiation at different ages. Under the same amount of XUV irradiation, the temperature profiles are similar to those obtained by linear-scaling XUV spectra shown in Figure 3, as well as the exobase temperature shown in Figure 6. Small differences are, however, found in the lower thermosphere (at ∼100 km), especially for younger stellar spectra, because young Sun-like stars emit larger amounts of FUV radiation, and FUV absorption mainly contributes to heating of the lower thermosphere for high-XUV conditions (see Figure 4). As such, the detailed shape of stellar UV spectrum has a small impact on the upper thermosphere where atmospheric escape rates are determined.

Zoom In Zoom Out Reset image size

The question of the vulnerability of the present-day Earth's atmosphere to strong XUV irradiation has long been explored, because the Sun is known to have emitted much stronger XUV radiation in the past than at present (e.g., Ribas et al. 2005). Öpik (1963) predicted that the Jeans escape parameter exceeds the critical value of 1.5 at a high level of XUV irradiation, leading to an uncontrolled process of atmospheric escape, and defined such a state as the atmospheric blow-off. With detailed treatments of photochemistry and radiative transfer and, thereby, self-consistent computations of heating efficiencies, Kulikov et al. (2007) modeled the hydrostatic equilibrium structure of the thermosphere of the terrestrial planets with enhanced XUV irradiation. Their model predicted a thermospheric temperature of ∼10,000 K for F_XUV ∼ 10F_XUV,⊕, for which the escape parameter is over the blow-off threshold.

Also self-consistently computing heating efficiencies, Tian et al. (2008a) investigated the response of an Earth-like N₂–O₂ atmosphere to an increase in the XUV irradiation level. First they demonstrated that no hydrostatic equilibrium solution is found beyond F_XUV ≃ 5F_XUV,⊕ and confirmed that the atmospheric blow-off state predicted by Öpik (1963) is achieved. Meanwhile, including the effect of advection, they did not find any uncontrolled process such as that imagined by Öpik (1963). Beyond the hydrostatic blow-off threshold, the adiabatic cooling due to advection (or outflow) becomes effective in keeping the thermosphere at relatively low temperatures. Then, as the XUV irradiation increases, advection becomes strong, and the thermospheric temperature decreases. Consequently the atmospheric escape occurs in an energy-limited fashion (Tian 2013). Recently, Johnstone et al. (2019) developed fully hydrodynamic models of planetary upper atmospheres and confirmed that the atmospheric escape occurs in the form of the transonic hydrodynamic escape (or the Parker wind) for the solar spectrum supposed at 4.4 Gyr ago.

By contrast, our results show that an N₂–O₂ Earth-like atmosphere is much more resistant to high-XUV irradiation than predicted by previous studies (Kulikov et al. 2007; Tian et al. 2008a; Johnstone et al. 2019). This is due to the effect of atomic line cooling, as shown in Section 3. In our model, the thermosphere never enters such a hydrodynamic regime and remains almost in hydrostatic equilibrium until F_XUV = 1000 F_XUV,⊕. This has a great impact on the loss of the atmosphere. For example, while Johnstone et al. (2019) estimated that the mass-loss rate is 1.8 × 10⁹ g s⁻¹ for an N₂–O₂ Earth-like atmosphere irradiated by the 4.4 Ga Sun, our model for the same stellar spectrum estimates the mass-loss rate at 4.8 × 10⁴ g s⁻¹, which is smaller by approximately four orders of magnitude than the previous value.

A fundamental question associated with the Earth's evolution is whether the modern atmosphere would be stable in high-XUV environments of the early Earth. Dividing the mass of the Earth's atmosphere (5.3 × 10²¹ g) by the mass-loss rate (1.8 × 10⁹ g s⁻¹) that Johnstone et al. (2019) obtained using Claire et al.'s (2012) solar XUV spectrum at 4.4 Ga, one can estimate the lifetime of a 1 bar N₂–O₂ Earth-like atmosphere to be 0.09 Myr. Such an extremely short lifetime relative to planetary evolution timescales suggests that modern Earth's atmosphere could never have survived the active phases of the Sun and must have been formed recently. Therefore Johnstone et al. (2019) concluded that Earth's early atmosphere was significantly different in amount and composition from the current atmosphere.

Our model, however, predicts much longer lifetimes. In Figure 8, we show the estimated lifetime of a 1 bar N₂–O₂ atmosphere as a function of XUV irradiation level. The lifetime τ_life is found to depend on the XUV irradiation level in a somewhat complicated way: for F_XUV ≲ 10F_XUV,⊕, τ_life increases slightly with increasing F_XUV. Then, τ_life becomes sensitive to and decreases rapidly with increasing F_XUV until F_XUV ≃ 100F_XUV,⊕. Thereafter, τ_life is almost constant for F_XUV ≳ 100F_XUV,⊕. This dependence can be interpreted as follows, according to our simulation results: For low-level XUV irradiation (F_XUV ≲ 10F_XUV,⊕), only the H species such as H and H⁺ are escaping, while the major components, O and N, are not because their escape parameters are sufficiently large. For sufficiently small values of the escape parameter, the H effusion velocity itself is insensitive to an increase in exobase temperature (see Equation (8)), resulting in almost constant mass-loss rates. A slight decrease in mass-loss rate appears because the mass-loss rate is determined solely by the number density of H at the exobase, which is, by definition of the exobase, inversely proportional to the scale-height or the exospheric temperature.

Zoom In Zoom Out Reset image size

For 10F_XUV,⊕ ≲ F_XUV ≲ 100F_XUV,⊕, by contrast, N and O are dominantly escaping, which is shown in Figure 10 of Appendix B. For such relatively high-XUV irradiation, the exobase becomes comparable to the planetary radius; then, a rise in F_XUV results in an increasing escape surface area and decreases the gravitational potential at the exobase. Such feedback brings about the strong dependence of the mass-loss rate on F_XUV. (Note: The response of the escape velocity to the exospheric temperature and the expansion of the thermosphere is also discussed in Tian et al. 2008a). For F_XUV ≳ 100F_XUV,⊕, the escape rate is kept almost constant, because the atomic line cooling and radiative recombination suppress the expansion of the thermosphere, and, consequently, the effect of expansion is canceled out by the exospheric temperature due to the high ionization fraction, as discussed in Section 3.2.

For the same XUV spectrum as adopted by Johnstone et al. (2019; F_XUV ≈ 60 F_XUV,⊕), the lifetime is estimated to be as long as 3.5 Gyr. Also, for F_XUV = 50 F_XUV,⊕, the calculated mass-loss rate is 2.3 × 10⁴ g s⁻¹, corresponding to a lifetime of 7.4 Gyr, longer than the age of the Earth. Thus, based on the vulnerability to XUV radiation, we cannot rule out the possibility that the early Earth had the same atmosphere as today, although the early atmosphere likely differed from that of today, as suggested by other evidence (e.g., the faint young Sun problem; Goldblatt et al. 2009).

To discuss the resistance of exoplanetary atmospheres to stellar XUV radiation, we would have to factor in the evolution of the latter. Observations of stellar X-ray emission suggest that main-sequence stars become less active with age, and the X-ray evolution differs from star to star. The difference is likely due to initial stellar rotation. According to Tu et al.'s (2015) calculations for solar-mass stars, the "fast rotators" younger than a few hundred megayears emit X-ray radiation hundreds of times as strong as the present Sun, while X-ray emission from the "slow rotators" decays rapidly, and its flux decreases to tens of times the present Sun's emission in a few tens of megayears (see Figure 2 of Tu et al. 2015). In our model for F_XUV = 1000 F_XUV,⊕, the mass-loss rate is estimated at 8.5 × 10⁴ g s⁻¹, corresponding to the lifetime of 2.0 Gyr for a 1 bar N₂–O₂ Earth-like atmosphere (see Figure 8). For such significantly long lifetimes relative to the stellar active periods, even for extremely high-XUV irradiation, our results suggest that exo-Earths orbiting in the habitable zone around Sun-like stars can retain atmospheres like the Earth's present atmosphere, regardless of stellar XUV activities.

As mentioned in the introduction, the initial active phases (or the saturation phases) of main-sequence stars increase with decreasing stellar mass. Especially for stars lighter than 0.4 M_⊙, the saturation phase lasts on geological timescales of gigayears and continues emitting XUV hundreds of times the solar XUV (e.g., see Johnstone et al. 2021). Thus, planets currently located in the habitable zone around M dwarfs, targeted by recent and near-future exoplanet surveys, have been exposed to strong XUV radiation. A 1 bar N₂–O₂ Earth-like atmosphere would be lost in the saturation phases. Nevertheless, N₂–O₂ atmospheres of a few bars may survive, because the saturation phase ends no later than a few gigayears, even for very-low-mass stars of ≲0.1 M_⊙ (see Figure 10 of Johnstone et al. 2021).

In reality, diverse atmospheric composition and pressure must be considered. Many pieces of geological evidence suggest that the Earth's atmosphere contained a significant amount of CO₂ or other greenhouse molecules in the past (e.g., Catling & Zahnle 2020). Beyond the solar system, diverse water contents in terrestrial planets predicted by planet formation theories (e.g., see O'Brien et al. 2018; Ikoma et al. 2018, for recent reviews) would result in diverse amounts of atmospheric CO₂ for terrestrial planets in the habitable zone (e.g., Nakayama et al. 2019; Krissansen-Totton et al. 2021). Since molecules are capable of cooling the atmosphere through the rotation-vibration transitions, oxidizing atmospheres with larger amounts of CO₂ are more resistant to XUV radiation (e.g., Johnstone et al. 2018). As seen in Section 3, however, the molecular IR cooling makes a smaller contribution to the thermospheric structure in our model than in previous models for highly irradiated atmospheres. Instead, enhanced atomic line cooling from carbon atoms could lead to more efficient cooling; quantitative investigation will be explored in forthcoming studies.

4.3.1. Stability of Water Vapor Atmosphere

4.3.2. Nonthermal Escape

Observational constraints on physical processes occurring in the thermosphere are crucial in further understanding of the evolution of planetary atmospheres. The compact thermosphere resulting from atomic line cooling will be confirmed through near-future exoplanet observations. For instance, Tavrov et al. (2018) proposed a strategy for constraining the thermospheric structure for strongly XUV-irradiated terrestrial planets: The resonant scattering of the OI triplet at a wavelength of ∼130 nm in O-rich planetary exospheres produces large transit depths during primary transits. Since the transit depth reflects the size of the thermosphere plus exosphere, future UV transit observations such as WSO-UV help us understand whether atomic line cooling is effective in actual systems.

In addition, our model predicts strong line emission in the NUV to optical wavelength range from the thermosphere of O-rich atmospheres in violent XUV environments. For instance, assuming a globally uniform structure, we estimate that such planetary atmospheres emit radiation with an intensity of up to 2 × 10²⁰ erg s⁻¹ at wavelengths of 630.0 nm and 636.4 nm via O(¹D) → O(³P) and at 557.7 nm via O(¹S) → O(¹D). Such strong emission would be detectable during secondary transits for M dwarf systems because M dwarfs are less bright in the optical wavelength, leading to high planet–star contrast ratios. Quantitative investigation of the observational feasibility is beyond the scope of this paper and will be the subject of future work.