© The Authors 2025

1 Introduction

Hot Jupiters (HJs) constitute a significant portion of all confirmed exoplanets. They show remarkable diversity in the conditions of their upper atmospheres. Such diversity is caused by an interplay of the planetary mass, proximity to the star, and the type of host star. Due to their large size and short orbital period, HJs are good targets for spectrally resolved observations with modern space- and ground-based instruments. A unique feature of HJs is an extended upper atmosphere heated by the short-wave radiation from a host star. Transit absorption in the Lyα line measured with the Hubble Space Telescope revealed massive hydrodynamic outflows, as well as dense partially ionized plasma that overflows the Roche lobe for a number of typical HJs, such as HD 209458 b (Vidal-Madjar et al. 2003) and HD 189733 b (Lecavelier des Etangs et al. 2010). The depth and spectral width of absorption in the vacuum ultra-violet (VUV) resonant lines of O I, C II, and Si III measured for HD 209458 b (Vidal-Madjar et al. 2004; Linsky et al. 2010) provided the evidence of a supersonic outflow of these trace elements carried by the escaping planetary wind. Recently a wealth of information was obtained in observations of the metastable helium line at 1083 nm. Out of approximately 36 exoplanets, half yielded a positive detection of absorption (Fossati et al. 2023). Numerical simulation shows that typically this absorption is generated by an upper atmosphere inflated to about three planetary radii and heated to about 10⁴ K (Khodachenko et al. 2021; Rumenskikh et al. 2022).

The planet KELT-9 b, a so-called ultra-hot Jupiter (UHJ), is unique in several aspects. It is a massive planet with an escape speed of about 50 km s⁻¹, which prohibits the significant outflow velocity inside the Roche lobe expected for most exoplanets at orbits as close as 0.35 au, and a related equilibrium temperature of more than 4000 K. It orbits a massive A0-type star with a spectral energy density (SED) that is very different from that of late-type stars. In particular, the radiation flux in the X-ray band (λ < 50.4 nm) is practically absent, whereas the moderate value of the Extreme Ultraviolet (EUV) flux (λ < 91.2 nm, with a value estimated via stellar modeling to be of the order of ∼5 erg cm⁻² s⁻¹ at a reference distance of 1 AU), is followed by a steeply increasing VUV flux of up to as much as 10⁷ erg cm⁻² s⁻¹ at λ < 300 nm. In Fossati et al. (2018), the SEDs of intermediate-mass stars were reviewed, and the authors showed that, in general, the X-ray and EUV (XUV) fluxes relevant for upper atmospheric heating are higher for cooler stars and lower for hotter ones, while the VUV fluxes increase with the increasing stellar temperature. This is explained by the supposed suppression of convection and the degeneration of corona at the stars with a surface temperature higher than 8500 K. As the number of extremely irradiated UHJs discovered around intermediate-mass F5- to B5-type stars rapidly increases (e.g., WASP-33 b, KELT-20 b, WASP-189 b, MASCARA-1 b, KELT-9 b, TOI-2046 b and TOI-2046 b. See Cameron et al. 2010; Gaudi et al. 2017; Lund et al. 2017; Talens et al. 2017; Anderson et al. 2018; Wong et al. 2021; Kabáth et al. 2022), it is important to study the specific differences of such stellar-planetary systems in comparison to the typical well-known HJs around solar type G-stars.

Since its recent discovery (Gaudi et al. 2017), an outstanding volume of observational data have been obtained for KELT-9 b. Besides the detection of a number of heavy elements (e.g., Fe I, Fe II, Ti I, Ti II, Sc II, Cr II; Hoeijmakers et al. 2019; D’Arpa et al. 2024), the population of excited hydrogen has been observed in Balmer lines with several different instruments (Yan & Henning 2018; Wyttenbach et al. 2020; Cauley et al. 2019; Turner et al. 2020; D’Arpa et al. 2024). The analysis in Fossati et al. (2020) gives the following absorption parameters averaged over all observations for the most important Hα line: depth A = 1.1% and full width at half maximum (FWHW) = 39 km s⁻¹.

Recently absorption by the higher level Paschen line was also detected (Sánchez-López et al. 2022). Moreover, the frist detection of absorption by excited O I at 777.4 nm was reported (Borsa et al. 2022) with a depth of A = 0.26% and a width of FWHW = 21 km s⁻¹, which may become a new channel for probing exoplanetary atmospheres. The measurements of transit absorption at the lines of two different elements, for example, H I and O, provide a challenge for the application of numerical models that simulate the stellar-planetary environment and the corresponding spectral observations, while the unique features of particular stellar-planetary systems may lead to the discovery of novel processes and features that have not yet been explored.

The detection of an optically thick Hα absorption in the atmosphere of KELT-9 b indicates that it takes place at heights of up to ∼1.6 R_p at a temperature of T ∼ 10⁴ K; it also indicates that hydrogen fills up the planetary Roche lobe (≈2.0 R_p) and escapes at a rate of ∼10¹² g s⁻¹ (Wyttenbach et al. 2020; Fossati et al. 2020). The crucial question is how and to what extent is this escape related to the population of excited hydrogen that we see in the observations.

García Muñoz & Schneider (2019) have shown that the specific SED of an A-type star dramatically changes the heating mechanism of the upper atmosphere of a close-orbit exoplanet. Photoionization from excited levels, rather than from the ground state, becomes dominant. Thus, for the hydrogen dominated atmospheres, the excitation of H I(2s, 2p) states by stellar radiation, electron collisions, and recombination drives the heating of the upper atmosphere. Simulations including metals revealed that elements with low-ionization potential, such as Mg and Fe, may significantly contribute to the overall heating, as well as cooling, processes (Huang et al. 2017; Fossati et al. 2021; Fossati et al. 2022a). Recently, Nakayama et al. (2022) demonstrated that the metastable level of O I(2s² 2p⁴) with a low excitation potential can be very important for the energy budget of the upper atmospheres of hot oxygen-rich exoplanets.

It appears that the excited levels related to atomic and ionic lines whose absorption in upper atmospheres was detected (e.g., Hα, He I(2³S), O I(5S), Fe I) cannot be in general described within the local thermal equilibrium (LTE). In this regard, we also note that the partially ionized plasma of hot exoplanets does not follow the Saha equilibrium as well. The reasons for this are the relatively low densities of the electrons and neutrals, which are insufficient for collisional as well as radiation equilibrium; as a consequence, the dominating ionization mechanism is by photoionization, whereas radiative recombination proceeds with the reabsorption of radiation and the excitation of levels. However, the opacity effects do play an important role for some of the transitions, such as the trapping of Lyα photons in dense layers (García Muñoz & Schneider 2019; Huang et al. 2017; Miroshnichenko et al. 2021).

To simulate the kinetics of the population of levels of different elements along with an account of radiation transfer is a formidable task. Lothringer et al. (2018) used the stellar and planetary atmosphere code PHOENIX to compute the temperature profile of KELT-9 b and showed that, for UHJs, the temperature is formed due to the absorption of intense UV and optical stellar radiation. The numerical package CLOUDY is a widely used platform to calculate chemistry, levels populations, and radiation transfer in mutlicomponent astrophysical plasmas. Fossati et al. (2020) employed CLOUDY with a set of fixed temperature-density (TP) profiles to calculate absorption in the Hα and Hβ lines and to further constrain the TP profiles of KELT-9 b. Significantly, a deviation from the LTE population of the H I(n2) level of several orders of magnitude was observed in the simulations, which fitted observations much better than that obtained with the LTE assumption. In a follow-up paper, Fossati et al. (2021) used a HELIOS radiative-convective equilibrium code to obtain the TP and density profiles at high pressures (>10⁻⁴ bar), combining it with CLOUDY non-LTE at lower pressures and calculating the TP profiles self-consistently. They found that the non-LTE treatment by CLOUDY generates a significantly hotter upper atmosphere, that is, up to 8500 K at pressures of 10⁻⁶−10⁻¹⁰ bar, versus only 6500 K obtained with the LTE treatment. Interestingly, this is mostly due to the contribution of Fe II ions and their low-lying levels being overpopulated, while other species, including excited hydrogen and Mg, provide a rather moderate input. At the same time, at pressures below 10⁻¹⁰ bar the excited hydrogen becomes the dominant heating driver. A higher temperature of the upper thermosphere results in a much better fit with observations of the simulated absorption in Hα and Hβ, as well as other Balmer lines (D’Arpa et al. 2024).

It appears that KELT-9 b is one of the few exoplanets in whose atmosphere the negative hydrogen ion H I⁻ can play a significant role (Arcangeli et al. 2018). The only element that is capable of efficiently absorbing the SED in the VUV range at a sufficiently dense atmosphere is a well-known bound-free continuum opacity of H I⁻(Wyttenbach et al. 2020). At the relative population of H I⁻/H I ∼ 10⁻⁸ and the scattering cross section of ∼10⁻¹⁷ cm⁻², it makes the atmosphere at pressures >0.01 bar optically thick for photons with λ < 300 nm. However, at such pressure, the population of H I⁻ is not in the LTE. Therefore, it has to be calculated using an account of the processes of H I⁻ creation and annihilation (Lenzuni et al. 1991).

The extreme VUV flux of the KELT-9 b host star makes this planet unique in another aspect. It is well known that, depending on the optical thickness of a particular atmospheric layer, a type A or type B recombination coefficient should be used. The populating of excited states due to type B recombination is usually several times higher than that in the case of direct recombination into the ground state. In a dense gas, where the resonant photons are trapped, the excited states relax into the ground state via electron collisions, thus finalizing the recombination act. However, when the process of the photoionization of excited states is faster than the electron collisions, the type B recombination rate decreases. This leads to an ionization-recombination balance at a higher degree of ionization. We found that this factor facilities much higher electron densities and, consequently, higher densities of excited hydrogen, as well as higher temperatures. This feedback process, pointed out in García Muñoz & Schneider (2019), remains overlooked in most studies and simulations of HJs.

García Muñoz & Schneider (2019) is the only work so far devoted to the upper atmosphere of KELT-9 b, where the level kinetics of the hydrogen atom is combined with the detailed radiation transfer of Lyα photons. The trapping and diffusion of Ly α is crucial for the calculation of heating and ionization via excited hydrogen, because at pressures above 10⁻¹⁰ bar, it increases the population of the H I(2p) state by orders of magnitude. This feature was not included in the simulations by Fossati et al. (2020); Fossati et al. (2021), so the overall heating they obtained could be significantly underestimated. So far, an account of the hydrogen level kinetics and Lyα diffusion has only been done in 1D models, such as Christie et al. (2013) (for HD 189733 b and HD 209458 b), Huang et al. (2017) (for HD 189733 b), and García Muñoz & Schneider (2019) (for KELT-9 b), and in an axis-symmetric 2D model by Miroshnichenko et al. (2021) (for HD 189733 b and HD 209458 b) and Sharipov et al. (2023) (for WASP-52 b).

The aim of the present study is to use a single unified self-consistent code to model the absorption in lines with the available observations, which reflect the parameters of the thermosphere and the upper atmosphere, specifically in Hα (656.3 nm), Paβ (1282 nm), and O I (777.4 nm). We note that the Hα and Paβ lines together probe two excited levels of hydrogen, H I(n2) and H I(n3).

We use our global 3D multifluid hydrodynamic (HD) model “exo3D”, which enables a calculation of the synthetic absorption without any geometrical approximations. This is important because the structure of the atmosphere where absorption takes place is influenced by tidal forces, anisotropic radiation flux, and stellar wind (SW) plasma flow.

We explicitly model the population of excited hydrogen atoms based on all the reactions of excitation and de-excitation, we restrict ourselves to the n = 2, 3, 4 levels, and we repeat the approach and employ the rate values of García Muñoz & Schneider (2019). The 2s and 2p sublevels are calculated separately taking into account the fast mixing between them induced by atom-proton collisions. The effect of resonant photon trapping (reabsorption) is included via semi-empirical analytic expression derived from the comparative analysis of direct Monte-Carlo simulations of Lyα scattering and diffusion. We note that suchan approach has not been attempted before, but appears to be rather efficient when applied to the exponentially scaled atmospheres. We note that previously we simulated the absorption in the metastable helium line at 1083 nm with the same code and exactly in the same generic way, by calculating all excitation and de-excitation reactions (e.g., in Khodachenko et al. 2021; Shaikhislamov et al. 2021; Rumenskikh et al. 2022).

The population of levels of minor trace elements is simulated by the Single Level Population Model (SLPM), which excludes interaction between levels. This approach is useful for the resonant transitions that are mostly responsible for line cooling and heating by photoionization. For KELT-9 b, it is also valid in view of the fast photoionization of the excited states. To assess the role of heating and/or cooling due to trace metals and to include this effect directly into the 3D simulation, we use the SLPM to model the population of the lowest levels of O I, C I, Mg I, as well as C II and Mg II ions. Those elements are considered to be most important, especially for cooling (Huang et al. 2017; Fossati et al. 2021; Nakayama et al. 2022). At the same time, Fe I and Fe II particles, which Fossati et al. (2021) found to be the main heating agents, have energy level systems that are too complex for such an approach.

To evaluate the validity of the SLPM approximation, we employ a dedicated platform, ASTERA (analogous to CLOUDY), which we are developing to simulate multilevel populations of elements. This platform has been constructed on the databases developed previously for Sun and stellar studies (Sitnova et al. 2013). It was also benchmarked versus CLOUDY in application to a simpler problem of an optically thin gas cooling function (Gnat & Ferland 2012). With regard to KELT-9 b, using the ASTERA platform we elucidate the fact that the population of O I(3s 5s) level, which is responsible for the absorption at 777.4 nm, can be calculated independently of other excited levels of O I, which thereby validates the SLPM approach.

The paper is organized as follows. In Sect. 2, we briefly describe our numerical model and its novel features. In Sect. 3, the simulation results are presented, subdivided into subsections devoted to particular aspects. Section 4 contains the discussion and conclusions.

2 Numerical model and basic issues

The global 3D multifluid HD model used in the present work has been already described in our earlier papers (Shaikhislamov et al. 2018, 2020b,a, 2021; Khodachenko et al. 2019; Khodachenko et al. 2021). It was developed as an upgrade of the previous 1D (Shaikhislamov et al. 2014) and 2D (Khodachenko et al. 2015, 2017; Shaikhislamov et al. 2016) models. The model code numerically solves the hydrodynamic equations of continuity, momentum, and energy for all species of the simulated multicomponent plasmas. Among the considered species, the model includes hydrogen and helium particles (H, H⁺, He, He⁺, He²⁺) as well as the heavier particles, O, C, Mg, and Fe in the atomic and up to Z = +2 ion stages. The model also includes molecular hydrogen species, such as ${\mathrm{H}}_2,{\mathrm{H}}_2^{+}$ , and ${\mathrm{H}}_3^{+}$ (Shaikhislamov et al. 2018; Khodachenko et al. 2019), however, in the present study we do not take them into account, because of the relatively high base temperature of the KELT-9 b atmosphere. Without accounting for molecular species, the photochemistry of H, He, and trace elements is driven by photoionization, dielectronic recombination, electron impact excitation, and ionization.

We describe the calculation of the population of the selected levels of elements with the SLPM in Appendix B. We take into account all the population and depopulation processes, except for the transitions between the excited levels. For hydrogen, this model is generalized to account for all possible transitions between its 1s, 2s, 2p, 3, and 4 levels. Within our 3D multifluid HD code, the populations of levels are treated as separate fluids with the same velocity and temperature as ground-state atomic hydrogen. Below, we describe some of the most important physical aspects and the related modeling issues.

2.1 Photoionization heating

Besides its effect on the degree of ionization, photoionization also results in strong heating of the gas by the photoelectrons produced, which in turn drives the hydrodynamic outflow of the planetary upper atmosphere. The novel aspect implemented in the present study, as compared to previous extensive multispecies modeling of, for example, HD 209458 b (Shaikhislamov et al. 2020b; Khodachenko et al. 2021), is that all photochemistry processes are taken into account, not only for the ground state of species, but also for their lowest excited states as well. In particular, in the case of hydrogen, its next most important state is the H I(n2) level with an ionization threshold of E_thr = E_ion − E_n2 = 3.4 eV. Its photoionization is driven by Balmer continuum absorption (n = 2 → ∞) at λ < 364.6 nm. To ensure a weak dependence of the atmospheric escape solutions on the inner boundary conditions, the corresponding stellar flux should be mostly absorbed at the base pressure. As is known from the stellar models, effective continuum absorption at these wavelengths (and up to an IR range <1600 nm) is provided by the negative hydrogen ion H I⁻ (Wildt 1939). Fossati et al. (2021) accounted for this within the CLOUDY model and H I⁻ provided the largest absorption at pressures >10⁻⁴ bar. In our model, we calculate the total heating term by integrating the absorption with photoionization by all considered species (in the ground and excited states), including H I⁻, over the wide range of the XUV+VUV+Visible spectrum according to wavelength-dependent cross sections. The attenuation of flux in the atmosphere is calculated based on the bound-free transitions, but only from the ground state of elements. For excited levels of hydrogen, the atmosphere remains optically thin up to very large pressures where the absorption by H I⁻becomes dominant. As a proxy for the SED function F_St(λ) of the host star, KELT-9, we use the PHOENIX spectrum (Husser et al. 2013). The corresponding stellar radiation flux scaled at 1 AU contains ≈4 erg cm⁻² s⁻¹ in XUV (λ < 91.2 nm) and ≈2 × 10⁷ erg cm⁻² s⁻¹ in VUV (100 < λ < 300 nm). The equations for SED absorption, ionization, and heating by photoelectrons are $$F_{\lambda }=F_{\mathrm{S}\mathrm{t}}(\lambda )\exp [-\int {\displaystyle \sum _k}{n}_k{\sigma }_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}(\lambda )dL],$$(1) $$R_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}=\int {\sigma }_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}(\lambda )F_{\lambda }d\lambda ,\lambda <\hslash c/E_{\mathrm{t}\mathrm{h}\mathrm{r},k},$$(2) $$W_{\text{heating\hspace{0.17em}}}=\int \left[{\displaystyle \sum _k}{n}_k{\sigma }_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}(\hslash c/\lambda -E_{\mathrm{t}\mathrm{h}\mathrm{r},k}-1.5k{T}_{\mathrm{e}})\right]F_{\lambda }d\lambda .$$(3)

We note that in Eq. (3) we included the thermal energy of background electrons, because the photoelectrons only lose available excess energy above this level. This simple fact, which has been mostly overlooked in aeronomy studies, is crucial for KELT-9 b. The very steep decrease towards larger photon energies in the SED of KELT-9 has an important consequence for the heating of the atmosphere. The mean energy of released photoelectrons, produced by photoionization of a particular specie, is given by $${\overline{E}}_{\mathrm{e}\mathrm{l},\mathrm{p}\mathrm{h}}=R_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}^{-1}\int {\sigma }_{\mathrm{p}\mathrm{h}-\mathrm{i}\mathrm{o}\mathrm{n},k}(\hslash c/\lambda -E_{\mathrm{t}\mathrm{h}\mathrm{r}})F_{\lambda }d\lambda .$$(4)

This energy equalizes with the background electrons, and the energy gain per photoelectron is Ē_el,ph − 1.5 kT_e. It appears that Ē_el,ph for the KELT-9 SED is rather small, which thus restricts the maximal possible heating, as shown in Fig. 1. For example, photoionization from the ground state of hydrogen, as the usual main channel of heating, produces photoelectrons with a mean energy of ∼0.6 eV capable of heating gas only up to ∼5000 K. This temperature is actually close to the supposed equilibrium temperature of KELT-9 b. In comparison, for the solar SED, the mean energy of photoelectrons is much higher: 2.6 eV. Altogether, this means that if the actual temperature of the planetary atmosphere is higher than 5000 K, the photoionization of hydrogen H I(n1) by the radiation of the host star, KELT-9, will ionize and cool the atmospheric gas at the same time. Some values that characterize the processes in the KELT-9 b atmosphere are listed in Appendix A.

2.2 Recombination-ionization balance

While photoionization of the upper atmosphere of KELT-9 b results in its partial cooling, the recombination process, on the other hand, leads to effective heating of plasma. The recombination decreases the number of particles, which results in a decrease of the total thermal energy, but the average thermal energy of the electrons actually increases. This is because the recombination cross section is larger for cooler electrons. The average energy of the recombining electrons appears to be about 0.7 kT_e (instead of 1.5 kT_e), so that each act of recombination increases the thermal energy of the remaining electrons by about 0.8 kT_e (Draine 2010). Accounting for the corresponding term in the energy equation becomes important, especially for nearly static atmospheres with a dominating recombination-ionization balance, such as that of KELT-9 b, rather than advection-ionization balance.

The planet mass of KELT-9 b is very high, and its atmosphere cannot accelerate to a significant outflow velocity inside the Roche lobe, despite the integral mass-loss rate appearing to be very high under the impact of radiation heating. This means that we can assume a barometric equilibrium between thermal pressure and gravity for the KELT-9 b upper atmosphere. This implies that there is a local balance between total ionization and recombination. We note that this is not a trivial assumption since for many exoplanets, for example HD 209458 b, the upper atmospheric outflow disrupts such a balance making recombination insignificant in comparison with ionization. The fact that the atmospheres of massive planets are in recombination-ionization balance, in contrast to the atmospheres of moderate mass planets being in advection-ionization balance, has been analysed in Erkaev et al. (2022); Fossati et al. (2018).

For the barometric exponential atmosphere with scale height H, the LOS (line-of-sight) density integral along the planet-star direction can be approximated as $$\int ndl=n(r)H,H\approx kT{R}_{\mathrm{p}}^2/\left(m_{\mathrm{i}}{M}_{\mathrm{p}}G\right).$$(5)

For simplicity, let us consider hydrogen alone, with photoionization of the ground and the excited states by the stellar radiation. The photoionization of the ground level decreases at the density n_a>(σ_thr H)⁻¹ ≈ 3 × 10⁸ cm⁻³ as the gas becomes optically thick for photons with energy E_thr > 13.6 eV. Here, σ_thr ≈ 7 × 10⁻¹⁸ cm² is the cross section at the threshold, and H ≈ 5 × 10⁸ cm is the scale height of KELT-9 b. At the same time, the photo-ionization of the excited states proceeds in an optically thin regime, down to much deeper layers of the atmosphere. For the photons generated by radiative recombination to the H I(n1) state, the gas becomes optically thick at the same densities, n_a > (σ_thr H)⁻¹. The loss of such trapped photons in the stratified atmosphere proceeds due to their diffusion in the direction of a negative gradient, or absorption by metals, if present. The effective reduction of the recombination by diffusion can be estimated as a factor of ∼(2n_a σ_thr H)⁻². Therefore, the local recombination rate to the H I(n1) state is $$R_{n1}^{\mathrm{r}\mathrm{e}\mathrm{c}}\approx R_{n1}^{\mathrm{r}\mathrm{e}\mathrm{c},0}/[1+{\left(2n_{\mathrm{a}}{\sigma }_{\mathrm{t}\mathrm{h}\mathrm{r}}H\right)}^2].$$(6)

Here $R_{n1}^{\text{rec},0}$ is the recombination rate to the ground state in an optically thin plasma The index notations used further on are specified in Appendix B. For a hydrogen atom, $R_{n1}^{\text{rec},0}$ is about 35% of the total rate. We note that with increasing density the local recombination rate to the ground state rapidly becomes insignificant as compared to photoionization by SED. The latter decreases much more slowly due to LOS absorption. We note that expression (6) does not provide a formula for accurate calculations, but rather gives a switch-off criterion due to the strong dependence of the recombination process on the density in the exponential atmosphere.

Let us now consider type B recombination into the bound states with n > 1. The ns states do not have direct radiative transitions to the ground state, 1s. However, they rapidly transit to the 2p state, while emitting a photon that escapes the system. The np states have resonant transitions to the ground state, but in a relatively dense plasma all resonant photons are trapped and undergo a cycle of scattering (emission and reabsorption).

As the average most probable cross section of absorption of the resonant photons is orders of magnitude larger than the photoionization cross section, the overall time it takes to radiatively decay to the ground state, in the regime in which n_a σ_thr H > 1, is also very large. At the same time, all states with n > 2 have a rapid mixing between sublevels and radiative transitions to 2s and 2p states. This leads to the accumulation of all recombined atoms in these two states (for a more detailed description, see Draine (2010).

When the fast radiative decay from the 2 p state is restricted due to the trapping of resonant photons, depopulation of the excited states takes place via electron collisional transition to the ground state, $R_{21}^{\text{coll\hspace{0.17em}}}$ , or by photoionization, $R_{2f}^{\mathrm{p}\mathrm{h}}$ . The first process leads to a finalization of the act of recombination, as the electron transfers from the free state to the ground state. On the other hand, the photoionization transfers the electron back to the free state, thus canceling the act of recombination.

Altogether, it turns out that the recombination into the excited states is also reduced in an optically thick plasma. This reduction can be estimated as $R_{2+}^{\mathrm{r}\mathrm{e}\mathrm{c}}\approx R_{2+}^{\mathrm{r}\mathrm{e}\mathrm{c},0}{R}_{21}^{\text{coll\hspace{0.17em}}}/(R_{21}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}+R_{2f}^{\mathrm{p}\mathrm{h}})$. The specific feature of the KELT-9 system is that photoionization from the excited states of hydrogen goes significantly faster than the collisional processes, resulting in $R_{2+}^{\text{rec\hspace{0.17em}}}\ll R_{2+}^{\text{rec,\hspace{0.17em}}0}$ . At the same time, under these conditions, the population of the H I(n2) level is close to LTE (see Eq. (D.8) in Appendix D). We note that for photons with an energy sufficient for photoionization from the excited states of hydrogen, the entire gas remains optically thin up to very large densities (for KELT-9 b up to at least ∼10¹⁷ cm⁻³). This is because the absolute density of excited hydrogen remains rather low.

From the analyses above, it follows that under conditions of an optically thick gas with n_a σ_thrH ≫ 1 and static atmosphere, the ionization-recombination balance requires a higher degree of ionization to compensate for the decrease in the effective recombination rate. At a sufficiently large $R_{2f}^{\mathrm{p}\mathrm{h}}$ , the gas can achieve the largest degree of ionization, which corresponds to the Saha balance. The performed simulations show that under the conditions of KELT-9 b, the dense layers of its atmosphere approach the Saha equilibrium.

A more detailed description of the model and further analytical analyses are provided in Appendix D. To verify the analytical estimates using the numerical model, in the Fig. 2 we demonstrate the results of the global 3D simulations of an H – He atmosphere of KELT-9 b. The profiles of various quantities are shown along the y axis, which is directed opposite to the orbital motion of the planet. One can see that the population of the excited state of hydrogen remains in LTE until the gas becomes transparent for type A recombination photons at pressures of ∼10⁻⁶ bar. The ionization degree in the same pressure range is very large and close to the Saha equilibrium, as well as the analytical formula based on type B recombination.

Let us undertake the estimate to demonstrate the effect of Lyα trapping and reduction of the effective radiative lifetime as described by Eq. (B.5), $R_{21}^{\text{rad,trap\hspace{0.17em}}}=R_{21}^{\text{rad\hspace{0.17em}}}{R}_{\text{loss\hspace{0.17em}}}/(R_{\text{loss\hspace{0.17em}}}+{\sigma }_{12}{n}_{0}c)$ . At a pressure of 10⁻⁶ bar, we see from Fig. 2 that T ≈ 7000 K, n_e ≈ 0.3n_a, and n_a ≈ 10¹² cm⁻³. The cross section of absorption of Lyα photons is ≈3 × 10⁻¹⁴ cm⁻³. Even for a direct time-flight loss, R_loss = c/H ≈ 10² s⁻¹, the reduction factor is 10⁻⁷, and $R_{21}^{\text{rad,trap\hspace{0.17em}}}\approx 63{\mathrm{s}}^{-1}$ . We note that an analysis of performed MonteCarlo simulations (see Appendix B) gave us a much smaller value of escape rate, R_loss ≈ 2 s⁻¹. At the same time, the rate of electron de-excitation is $R_{21}^{\text{coll\hspace{0.17em}}}\approx 1.5\times {10}^{-8}{n}_{\mathrm{e}}\approx 0.5\times {10}^4{\mathrm{s}}^{-1}$ .

The bottom panel in Fig. 2 demonstrates that the analytical estimate of the heating rate due to the H I(n2) state is also in good agreement with the simulation. At pressures above ∼10⁻² bar, the heating decreases due to absorption of the stellar radiation by H I⁻ ions.

Among hydrogen species, we also include negative ion H I⁻. Its main role is continuum absorption close to the photometric radius at pressures of p > 0.1 bar. A detailed description is given in Appendix C.

Fig. 1 Mean energy of photoelectrons (asterisks, in eV) calculated with Eq. (4) and a corresponding thermal equilibrium temperature of 1.5 kTe = Ēel, ph (dashed line) versus the threshold of photoionization calculated for the SED of KELT-9. For the cross section, a simple formula is used: σph-ion = σthr(Ethr/E)3, E > Ethr. The labeled vertical lines indicate the photoionization thresholds of several particular elements, or transitions, of interest.

Fig. 2 Profiles of various quantities of the H–He atmosphere of KELT-9 b along the y axis directed opposite to the orbital motion of the planet, obtained in simulation with the parameter set N 1 (Table 2 below) and expressed via pressure. Top panel: Population of H I(n2) level in terms of LTE value (red, left axis); ionization degree Ne/Na (orange) in comparison to Saha equilibrium (turquoise) and type B recombination (blue, formula (D.8)); temperature (gray, upper horizontal axis). Bottom panel: Heating (>0, solid lines) or cooling (<0, dashed lines) terms in the energy equation associated with H I(n1) and H I(n2) levels. Photoionization heating from the H I(n2) level given by the analytical formula $W=n_{{\mathrm{H}\mathrm{I}}_{\mathrm{I}}}(n2)R_{2f}^{\mathrm{p}\mathrm{h}}{\overline{E}}_{2f}^{\mathrm{p}\mathrm{h}}$ is shown with a green line.

Fig. 3 Cooling rates for Mg and O calculated with a unit particle density (i.e., 1 cm−3) for all species, including electrons, and with zero radiation field. Black lines show cooling due to the excited levels listed in the table, red lines show full cooling including ionization from the ground state, and, for comparison, yellow lines show the analogous simulation result with the code CLOUDY from Gnat & Ferland (2012).

Fig. 4 Cooling or heating rates of elements, calculated for KELT-9 b with an account of the lines listed in Table 1 at ne = 1011 cm−3 under the impact of an unshielded radiation flux of the host star. Depending on temperature, the rates are either negative (i.e., cooling, dashed lines), or positive (i.e., heating, solid lines). The vertical axis is in log scale.

2.3 Effect of minor trace elements

Let us now consider species other than hydrogen. Helium, with its large ionization potential, remains practically neutral in the upper atmosphere of KELT-9 b and it does not affect the energy balance there. Its only role is in influencing the total local scale height. Further on, we consider a fiducial abundance to be He/H = 0.03, which is somewhat less than the standard Solar figure and appears typical for hot Jupiters, judging by the absorption measurements of the He I 1083 nm line (Fossati et al. 2023).

In the simulation, we include C, O, Mg, as well as Fe in the ionization stages of up to 2+ for all elements. We also include several lowest excited states of the neutral and single ionized particles of C, O, and Mg, as listed in the table below. The population of these levels is simulated within a SLPM approximation by the same generic routine. This enables a more accurate calculation of the ionization and contribution of these elements to the energy balance. Since the maximum temperature of the upper atmosphere of KELT-9 b is below 1.5 eV, an account of only a few levels with a relatively simple structure and lowest excitation energy is sufficient for C, O, and Mg. This simplified approach can be justified by the smallness of the rates of collisional transitions between the excited levels (whether considered or not) in comparison to the rates of photoionization from the excited states. We note that here we do not consider any kinetics for O II, because the low energy levels of this ion (3.32 and 5.02 eV), which might affect the energy balance, are metastable, (one) and they have too large (>30 eV) an ionization energy gap. It means that these levels, which are not affected by photoionization and radiative decay, should be in LTE. The species of Fe I and Fe II play some a role by absorbing stellar radiation and providing additional electrons. Therefore, we take them into account in the ground states, but do not calculate the level kinetics because it is too complex for the SLPM approach that is implemented in our global 3D HD code.

As a first step and test illustration, in Fig. 3 we show the cooling rates of Mg I, Mg II, and O I per particle per electron obtained with our model in a zero radiation field versus the results of analogous calculations with the code CLOUDY (Gnat & Ferland 2012). As a next step, to approach the real situation, we calculate the cooling or heating by elements, based on the levels’ data listed in Table 1 (including the ground states) at a realistic density of electrons n_e = 10¹¹ cm⁻³ and with the radiation flux of the host star of KELT-9 b included without any absorption (self-shielding) and trapping of resonant photons. In this case, all processes, except recombination, are “switched-on”. The recombination is put to zero for this particular demonstration. One can see in Fig. 4 that the heating due to photoionization is rather negligible and it only takes place at relatively low temperatures. For the majority of the analyzed temperature range, the considered elements mainly provide cooling of the atmosphere. This is because photoionization from the excited levels, while producing photoelectrons, consumes the energy that is spent for the excitation of these levels. We note that hydrogen is not particularly different in that respect from other atmospheric species.

In Huang et al. (2017); Fossati et al. (2021); Fossati et al. (2023), magnesium is shown to be important for the cooling of hot exoplanet atmospheres. As the Mg I atom is very rapidly photoionized, which supplies most of the electrons in a dense atmosphere, we focus our attention on magnesium ions. The first level of Mg II has a rapid resonant transition to the ground state. The difference, as compared to the H I(n2) level, is that the recombination pumping of Mg II(n2) level can be ignored (due to low density of Mg III), as well as photoionization (due to the large energy of photons required, λ < 100 nm). Thus, we can expect that the population will change from LTE in dense atmospheres where resonant photons are trapped to significantly below LTE as plasma becomes optically thin. The next level, Mg II(n3), is metastable. It is excited collisionally, with subsequent rapid photoionization by photons with energy λ < 168 nm, and therefore this level is significantly underpopulated. Thus, both levels will contribute to cooling.

3 Simulation results

Using the exo3D code we adopted to study KELT-9 b, we calculate the transit absorption profiles of the Hα, Paβ, and O I 777.4 nm lines and compare them with observations. Table 2 lists the modeling parameters that have been varied, and the values that have been derived in a series of simulations. We note that the parameter set N1 was already discussed with respect to Fig. 2. For the fiducial parameters (set N5), we take the XUV (λ < 91.2 nm) and VUV (100 < λ < 360 nm) flux values at 1 a.u., respectively, F_XUV = 5 erg cm⁻² s⁻¹ and F_VUV = 2 × 10⁷ erg cm⁻² s⁻¹; the atmospheric base temperature and pressure as T_base = 4500 K, P_base = 100 bar; and helium content as He/H = 0.03. Since the SW is expected to be negligible for the KELT-9 star, we take it to be relatively weak with a small integral mass loss of 10¹⁰ g s⁻¹, which corresponds to the following model parameters at the planet’s orbit: V_sw = 220 km s⁻¹, T_sw = 8 × 10⁵ K, and n_sw = 10² cm⁻³. We note that SW does not affect the outflow of planetary atmosphere as well as the absorption in the lines of interest. Therefore, we do not discuss it further on.

The top panel in Fig. 5 shows an overall global picture of the escaping atmosphere of KELT-9 b interacting with the SW plasma flow at the scale of the whole stellar system. The streams of upper atmospheric material accelerated away from the planet, beyond the Roche lobe, are forced to flow along the planet orbit due to momentum conservation. At the same time, because of the interaction with the SW plasma, part of the stream ahead the planet loses momentum and accretes onto the star. However, due to the low density, plasma outside the Roche lobe does not absorb in the lines of interest. For a better view of the escaping atmospheric material structure in closer vicinity to the planet, comparable to the size of the Roche lobe, which is important for transit absorption, we “zoom” the simulation results to this scale and show the distribution of species there. In particular, the bottom panel of Fig. 5 shows the meridional plane cut of the density of the excited hydrogen atoms H I(n2).

In Fig. 6, we show the temperature and density profiles of various elements, including the population of excited hydrogen, along the planet-star line. One can see that temperature rises gradually, starting from a pressure of ∼0.1 bar (at h = 0.1 R_p) until the atmosphere spills over the Roche lobe (at h ≈ R_p). The ionization of hydrogen is very high, with a proton density of up to 10¹² cm⁻³. Close to the planet, the density of electrons, supplied by ionization of Mg and Fe, is even higher.

Figure 7 shows the distributions along the planet-star line of the population of different levels in comparison with those calculated by the LTE formulae. The hydrogen levels (Fig. 7, left panel) within the Roche lobe remain rather close to the LTE values and become different in the tail region. A comparison of the population profiles in the day and night sides leads to the following important conclusions. For the generation of excited H I, the heating and photoionization of that same excited H I are crucial, and constitute a feedback process. Without a high level of ionization supported by photoionization, the levels with n > 2 are significantly below LTE due to radiative de-excitation. The SED of KELT-9 in far ultraviolet (FUV) is very strong and can populate the excited levels by resonant photons. For example, while the SED of KELT-9 has no Lyα emission line, the continuum in the region of Lyα is strong. However, due to the absorption (self-shielding), direct pumping by the stellar resonant photons does not affect the population of most resonant levels of atomic species, for example H I(n2), Mg I(n2) and Mg II(n2).

However, one can see that the O I(n4) level population is different. In particular, at the day side it exceeds the LTE level, which is obviously due to excitation by the stellar continuum. At the same time, the collisional excitation that supports the LTE level is small due to relatively large excitation potential. Excitation by stellar radiation is rather strong in contrast to the case of strong resonant transitions due to a relatively low radiative decay rate and low self-shielding. The O I(n4) level is important because it causes absorption at the 777.4 nm line, and below we estimate its population at x ≈ 1.25 using the profiles in Fig. 6. The corresponding relevant rates and beta-factor (level population normalized by its LTE value) are $$\begin{array}{rl}& R_{0i}^{\mathrm{p}\mathrm{h}}\approx 4.8\times {10}^{-3}{\mathrm{s}}^{-1},R_{i0}^{\mathrm{p}\mathrm{h}}=2.8\times {10}^3{\mathrm{s}}^{-1},R_{if}^{\mathrm{p}\mathrm{h}}\approx 0.29\times {10}^3{\mathrm{s}}^{-1},\\ & R_{i0}^{\text{coll\hspace{0.17em}}}\sim {10}^3{\mathrm{s}}^{-1},R_{0i}^{\text{coll\hspace{0.17em}}}\sim {10}^{-4}{\mathrm{s}}^{-1},\text{and}\\ & \beta \approx \frac{R_{0i}^{\mathrm{p}\mathrm{h}}}{R_{i0}^{\mathrm{p}\mathrm{h}}}\frac{\exp \left(10.61/T_{\mathrm{e}}\right)}{g_{i0}}=3.1\times {10}^{-6}\exp \left(10.61/T_{\mathrm{e}}\right)\text{\hspace{0.17em}where}\\ & \beta >1\text{\hspace{0.17em}at\hspace{0.17em}}{T}_{\mathrm{e}}<8300\mathrm{K}.\end{array}$$

Thus, the population of the O I(n4) level at the day side should be above the LTE level due to excitation by the stellar continuum.

Figure 8 shows the actual heating or cooling rates due to elements and levels. It demonstrates that the heating due to excited H I remains by far the main energy input. It is balanced by cooling due to excitation and ionization from the ground state of H I. The channel that involves H I(n1) also includes the energy of 1.5kT_e, assigned to the electrons newly created in the course of the total ionization rate of H I, which is an effective energy loss.

In a dense atmosphere, heating is provided via other elements. In particular, above pressures of 10⁻² bar, the main contributors are C I and O I, whereas at very large pressures of 10⁻¹ bar, Mg I plays the main role. At pressures below 10⁻⁶ bar, cooling by Mg II lines makes some difference. The plots for individual levels show which of them makes the major contribution.

As already shown in Fig. C.1, at fiducial parameters we obtain a good agreement with the observed absorption profile of the Hα line, both in the amplitude and in the half-width. Next, we demonstrate how this absorption depends, if indeed it does, on some important processes. In Fig. 9 we show a number of profiles obtained in a series of test simulations. In particular, if VUV heating is excluded, the calculated absorption is significantly smaller, because the atmosphere becomes much cooler in this case. Therefore, the VUV part of the stellar radiation flux is indeed an essential factor in our particular study. When ignoring the average energy, 1.5kT_e, that should be assigned to the newly created electrons, the calculated absorption due to hotter atmosphere becomes significantly higher than the measured value. At the same time, without effective heating caused by recombination (0.7kT_e), the absorption is smaller due to cooler atmosphere. Switching off the XUV heating (not shown in Fig. 9) does not produce any effect on Hα absorption.

Figure 10 demonstrates the absorption profiles for three lines of interest obtained at several parameter sets that deviate from the fiducial set N5. One can see that absorption in the H I lines is strongly dependent on the VUV flux, as well as on helium abundance. At the same time, absorption in the O I line does not depend on heating, but it is sensitive to the scale height of the atmosphere. Independent variation of the VUV and He/H ratio leads to degeneracy, as the same observations can be fitted with higher or lower VUV and He/H as compared to the fiducial set. This is demonstrated in Fig. 11.

Besides the maximum of absorption, the full-width at half-maximum (FWHM) of the absorption profile is another important characteristic. For the Hα line measured at KELT-9 b, it is 44 km s⁻¹ and 52 km s⁻¹, according to Wyttenbach et al. (2020) and Yan & Henning (2018), respectively. The average weighted value of FWHM over all the Hα observations provided in Fossati et al. (2020) is 39 km s⁻¹. In our simulations, the FWHM of the Hα line is 48 km s⁻¹. For the O I triplet line at 777.4 nm, Borsa et al. (2022) report a FWHM of 21 km s⁻¹. We note that the 1D modeling performed in Borsa et al. (2022) required the inclusion of empirical micro- and macro-turbulence broadening of 3.0 ± 0.7 and 13 ± 5 km s⁻¹, respectively, to fit the observations. In contrast, our synthetic profile of the O I triplet has a FWHM of 20.5 km s⁻¹. Finally, for the Paschen line, the measured width is 32 [+12; −8] km s⁻¹ (Sánchez-López et al. 2022), while the calculated figure is 42 km s⁻¹. Therefore, our model fits the observed FWHM of all three lines rather well.

As demonstrated in Fig. 12, there is a direct linear relation between the absorption in hydrogen lines and mass loss. In fact, the mass loss obtained for KELT-9 b is one of the highest among the close orbit exoplanets.

Fig. 5 Top panel: distribution of the proton density in the equatorial plane of the KELT-9 stellar system (i.e., in the whole simulation domain) simulated with the fiducial parameter set N5. Black circle indicates the star in scale. The planet is located at the coordinate center. Streamlines are shown originating at the star (black lines), and at the planet (blue lines). Bottom panel: distribution of the excited hydrogen density in the meridional plane in close vicinity to KELT-9 b simulated with the parameter set N5. The scale of plots is in units of Rp. The values outside the specified variation ranges are colored either in red, if smaller than minimum, or in blue, if higher than maximum.

Fig. 6 Density profiles of different elements along the planet-star line (x > 0 spatial axis, expressed in units of pressure). The gray line with the upper horizontal axis shows the temperature (top panel) and the height above the planet h = R/Rp − 1 (bottom panel). Results are obtained with parameter set N5.

Fig. 7 Population profiles of different levels along the planet-star line at the day side (x > 0) and at the night side (x < 0). The left panel shows the hydrogen densities calculated by model (solid) lines and the corresponding LTE values (dashed). The center and right panels show the levels’ profiles as beta-factors (i.e., density values normalized to LTE values) for Mg I, Mg II, and O I, C I, C II.

Fig. 8 Heating (solid lines) or cooling (dashed lines) due to various elements and levels along the planet-star line (x > 0, expressed as a pressure). The upper panel shows the main considered elements of interest. The ground and excited levels of hydrogen are distinguished, whereas for other elements, all simulated levels are combined. The “sum” line indicates the accuracy of the overall balance between heating and cooling achieved in simulation. The gray line shows the temperature profile. The center and bottom panels show the heating or cooling contribution of excited magnesium (top), oxygen, and carbon (bottom), respectively.

Fig. 9 Absorption profiles in Hα line at mid-transit in units of Doppler velocity obtained in test simulations: when the stellar flux in VUV range is put to zero (red); when the residual energy of newly created electrons is put to zero (instead of 1.5kTe, green); when the recombination heating is put to zero (blue). For comparison, the black line shows the Hα absorption profile calculated with the fiducial set of the model parameters N5.

4 Discussion and conclusions

To conclude the reported study, we compare the results obtained with similar previous studies of KELT-9 b in more detail. Table 3 gives some characteristic values, while Fig. 13 shows the temperature and total heating profiles versus pressure. First of all, we list the conclusions that are common to all works. These are: 1) the main heating of the atmosphere of KELT-9 b occurs due to the VUV part of spectrum of the host star; and 2) the heating and mass loss of the planet driven by the XUV part of the spectrum only, related to the photoionization of atoms from the ground state, are very small.

Now let us consider the details that differ, sometimes significantly, but do not change the main results. The characteristic maximum temperature of the thermosphere varies from 0.85 × 10⁴ K to 1.5 × 10⁴ K and is achieved in a relatively narrow pressure range of 3 × 10⁻⁷−10⁻⁹ bar. We note that there is a good qualitative agreement between this work and Fossati et al. (2021). The point of half-ionization of hydrogen also lies in a relatively narrow pressure range of 3 × 10⁻⁶−10⁻⁸ bar. The total mass loss rate is estimated in Borsa et al. (2022) as 5 × 10¹¹ g s⁻¹. In the present work, it reaches much higher values of 3 × 10¹² g s⁻¹.

A significant difference in the work of García Muñoz & Schneider (2019) is that the calculation begins with a low inner boundary pressure of 10⁻⁷ bar and a rather high temperature of 8000 K. The calculated population of excited hydrogen, β ≈ 3 × 10⁻⁴, is much less than the LTE level. In other works, although much less than unity, the β-factor at similar pressures is nevertheless two orders of magnitude higher than that in García Muñoz & Schneider (2019). Probably for this reason, the absorption in the Hα line calculated in García Muñoz & Schneider (2019) is also small. Indeed, in Fossati et al. (2020); Fossati et al. (2021) it is indicated that the absorption region of this line lies in the pressure range of 10⁻⁵−10⁻⁸ bar. In the present work, we can specify this range more precisely. As shown in Fig. 14, half of the Hα absorption forms at the heights that correspond to the pressure range of 10⁻⁷−10⁻⁸ bar, whereas the heights with pressures of 5 × 10⁻⁵−10⁻⁸ bar already provide 80% of the absorption. The total heating of the atmosphere of KELT-9 b appears to be in surprisingly good agreement between the works, despite the differences of the models employed. However, in Fossati et al. (2020); Fossati et al. (2021) the main contributor to heating are excited levels of Fe II ions, while in other similar studies, excited levels of H I are identified as the main contributor.

It follows from Fig. 14 that absorption in the Hα line takes place in the optically thick regime. Indeed, the typical pressure where the major part of absorption takes place is above ∼10⁻⁷ bar, which corresponds to a density of H I(n2) > 10⁴ cm⁻³ (see Fig. 6). At a typical cross section of ∼3 × 10⁻¹³ cm² at the line center, the optical thickness is <3 × 10⁸ cm, which is much smaller than the extent of the atmosphere. At the same time, the optical thickness of the atmosphere for the Paschen Paβ line is much less than for the Hα line, though still larger than unity. To check whether the inferred absorption in the Hα line, provided by the simulated population of the H I(n2) level, is affected by the large optical thickness, we calculated the absorption in the Hβ line, which is also based on the population of the H I(n2) level, but has much less optical thickness due to the smaller cross section. The result is shown in Fig. 15. A sufficiently good agreement with the measurements provides an additional verification for the consistency of the modeling performed.

Regarding the processes that populate excited states related to observed absorptions, the H I levels are populated by collisional excitation and recombination, which strongly depend on temperature and electron density. The electron density in the upper atmosphere of KELT-9 b is significantly increased due to strong photoionization by VUV spectra. The other crucial factor is the trapping of resonant Lyα photons. On the other hand, the O I excited level related to the 777.4 nm line is mostly populated by stellar radiation pumping and only partially be collisions.

Altogether, the results of the present work confirm the conclusion by García Muñoz & Schneider (2019) that the VUV spectrum of an A-class star creates a fundamentally new source of heating of hydrogen atmospheres, which is related to photoionization from excited levels. Our study goes a step further than the previously reported studies and shows that this heating mechanism is sufficient to produce the absorption in the H I lines observed at transits of KELT-9 b. For the first time, we demonstrated that the observed absorption in the lines generated from the excited hydrogen is a direct consequence of the heating process associated with the population of the lower levels of these lines. Also, in simulations and analytically, we demonstrated that the population of excited hydrogen is in LTE before the temperature maximum, then after the maximum it becomes significantly lower. This is because Lyα photons are trapped in a dense atmosphere, which drastically reduces the radiative de-excitation of the excited H I. We have also shown that the heating or cooling due to such elements as O, C, and Mg is overall significantly less than by hydrogen, though at very large pressures it can dominate. Thus, the conclusion based on the presented modeling is that the excited states of hydrogen are able to provide sufficient heating rates to fit the observations, while the additionally considered heavy elements C, O, and Mg make the overall energy budget picture more complex, but essentially not too different. A detailed comparison of the simulations with the model parameter sets N1–N5 reveals that without metals, or without cooling due to metals, the absorption in H I lines significantly increases, and it does not fit observations at the fiducial VUV flux and helium abundance. Overall, the most cooling is produced by Mg II. Our modeling-based conclusion on the cooling efficiency of Mg is in agreement with other works, such as Huang et al. (2017), Fossati et al. (2021); Fossati et al. (2022a).

However, the main conclusion of the present paper, that heating by H I is sufficient to fit the observations, is still in disagreement with the results of Fossati et al. (2021); Fossati et al. (2022a), where Fe II was found to be the main agent of heating, while H I played a secondary role. The final judgment on this controversy remains open, as we did not simulate the kinetics of Fe. The corresponding extension of our model is a subject for future work.

By varying the helium abundance and the base temperature, which change the atmosphere scale height without affecting the heating, as well as the VUV flux, which is directly related to the amount of heat absorbed in the atmosphere, we found that there is a monotonic dependence of the absorption amplitude in the Hα line with the overall mass loss by the atmosphere. It appears that the observed absorption amplitude of ∼1% at KELT-9 b implies that the mass loss of the planet exceeds 10¹² g s⁻¹, which is a rather high figure.

Fig. 10 Absorption profiles in Hα, Paβ, and O I 777.4 nm lines at mid-transit in units of Doppler velocity for the simulations with different model parameter sets specified in Table 2.

Fig. 11 Absorption profiles in Hα, Paβ, and O I 777.4 nm lines calculated with different combinations of VUV flux and He/H abundance ratio, which enable a best fit of the observations.

Fig. 12 Absorption in Hα line averaged in the velocity interval ± 20 km s−1 around the line center versus the integral mass loss of atmosphere. Different points were obtained by independently varying the He/H abundance ratio within the range [0.03; 0.1], VUV integral flux within [X0.5; X2], and base temperature Tb within [4000; 5000] K. The gray circle highlights the points that best correspond to measurements.

Fig. 13 Total heating (red, bottom axis) and temperature (black, top axis) profiles obtained in the present work (solid lines), in Fossati et al. (2021, dotted lines), and in García Muñoz & Schneider (2019, dashed lines).

Fig. 14 Absorption at the center of the Hα line (at V = 0, black circles) and at the distant wing (V = 80 km s−1, white circles), manually restricted by calculation within region Rp < r < Rmax around the planet, in dependence on Rmax. The vertical dashed lines show the region of formation of the corresponding absorption at V = 0 at the level of 0.5 to 0.9 of the maximum value. The gray curve shows the pressure in the atmosphere at the corresponding altitudes.

Fig. 15 Absorption profiles in Hβ 486.1 nm line at the same model runs as those in Fig. 11 (N5 – red, N12 – olive, N13 – blue). The gray circles with error bars show the average measurements over the available observations (from Fossati et al. 2020).

Acknowledgements

This work is supported by project No. 23-12-00134 of the Russian Science Foundation. Parallel computing simulations, key for this study, have been partially performed at Computation Centre of Novosibirsk State University and SB RAS Siberian Supercomputer Centre. SSS and BAG acknowledge the support of RSF 23-72-10060 grant under which the Monte-Carlo simulation was performed. Authors are grateful to Luca Fossati for inspiring this work, valuable discussions and comments.

Appendix A Typical values and rates for hydrogen at conditions of KELT-9 b

Appendix B Equations of level population

Let’s consider the rate equations for the populations of two quantum levels of a given element (j): the ground state (0) and an excited state (i), with the energy E_i0 separating them, also involving the ionization state (j+1): $$\frac{d}{dt}{n}_i=-R_{i0}^{\mathrm{r}\mathrm{a}\mathrm{d}}{n}_i+R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_0-R_{i0}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_i+R_{0i}^{\mathrm{p}\mathrm{h}}{n}_0+c{\sigma }_{0i}{n}_{\mathrm{p}\mathrm{h}}{n}_0+R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}{n}_{\mathrm{e}}{n}_0^{j+1}-R_{if}^{\mathrm{p}\mathrm{h}}{n}_i-R_{if}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_i,$$(B.1) $$\frac{d}{dt}{n}_0=R_{i0}^{\mathrm{r}\mathrm{a}\mathrm{d}}{n}_i-R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_0+R_{i0}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_i-R_{0i}^{\mathrm{p}\mathrm{h}}{n}_0-c{\sigma }_{0i}{n}_{\mathrm{p}\mathrm{h}}{n}_0+R_{f0}^{\mathrm{r}\mathrm{e}\mathrm{c}}{n}_{\mathrm{e}}{n}_0^{j+1}-R_{0f}^{\mathrm{p}\mathrm{h}}{n}_0-R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_0,$$(B.2) $$\frac{d}{dt}{n}_{\mathrm{p}\mathrm{h}}=R_{i0}^{\mathrm{r}\mathrm{a}\mathrm{d}}{n}_i-c{\sigma }_{0i}{n}_{\mathrm{p}\mathrm{h}}{n}_0-{\hat{R}}_{\text{loss\hspace{0.17em}}}{n}_{\mathrm{p}\mathrm{h}},\text{\hspace{0.17em}and}$$(B.3) $$\frac{d}{dt}{n}_0^{j+1}=-(R_{f0}^{\mathrm{r}\mathrm{e}\mathrm{c}}+R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}})n_{\mathrm{e}}{n}_0^{j+1}+R_{0f}^{\mathrm{p}\mathrm{h}}{n}_0+R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_0+R_{if}^{\mathrm{p}\mathrm{h}}{n}_i+R_{if}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}{n}_i.$$(B.4)

Here,

n₀, n_i represent the densities of particles j in the ground and excited states, respectively,

$n_0^{j+1}$ is the density of particles j + 1 in the ground state,

n_ph is the density of photons resonant with 0 → i transition, and

σ_0i is the average cross section for absorption of resonant photons by the 0 → i transition.

The superscript on the rates R denotes the type of transition: radiative (rad), photoinduced (ph), electron collisional (coll), or recombination (rec). We note that the rate of collisional reactions is proportional to the electron density. The subscripts indicate the initial and final states of the electron during the transition: ground state (0), excited state (i), and free state (f).

${\hat{R}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}$ is the loss rate of resonant photons due to non-resonant absorption. In general, this is an operator that includes spatial diffusion. In this model, we include losses due to photoionization of H I(n2), H I⁻, Mg I(n1), and Fe I(n1). It is well established (Dijkstra 2019) that, due to diffusion in the frequency space, the spatial diffusion of resonant photons in an optically thick regime does not follow simple random walk statistics and proceeds much faster. Since we are dealing with the gas with an exponential density profile, it has a characteristic scale – the barometric scale height (H < 0.1 R_p) – which is much smaller than the size of the entire region (∼R_p) where the equations are solved. Moreover, due to the exponential density profile, the resonant photons, on average, move away from the planet because flight time before reabsorption is longer towards the less dense layers. We note, that there are two simple asymptotes. The fastest loss rate is approximately the photons’ inverse time-of-flight through the scale height, c/H. The slowest is the purely spatial diffusion occurring in the course of scattering with a drift due to the gradient in the mean free path in the stratified gas, R_loss ≈(c/H)(4 σ_0i n₀ H)⁻¹. In any case, the Eq. (B.3) can be eliminated by introducing an effective radiation decay rate $$R_{i0}^{\text{rad,trap\hspace{0.17em}}}=R_{i0}^{\text{rad\hspace{0.17em}}}\frac{R_{\text{loss\hspace{0.17em}}}}{R_{\text{loss\hspace{0.17em}}}+{\sigma }_{0i}{n}_{0}c}.$$(B.5)

To validate the approach used, we compare the result obtained by solving Eqs. (B.1–B.5) using a local approximation with the calculation of the H I(n2) population obtained using the Monte-Carlo code (Miroshnichenko et al. 2021; Sharipov et al. 2023), in which the frequency-spatial trajectory of Lyα photons is traced. In this case, the Monte-Carlo calculation is performed for the distribution of hydrogen atoms, protons, and electrons obtained from the 3D gas-dynamic model with the parameter set N1 (Table 2).

Photoionization of elements is the dominant process that alters the stellar SED as it penetrates into deeper layers of the atmosphere. We calculate this radiation transfer taking into account all elements, including H I⁻ ions. The photoionization cross sections of elements from the ground state are taken from Verner et al. (1996). The cross sections of electron collisional ionization of elements from the ground state are taken from Voronov (1997). Bound-bound electron collision rates are related through the LTE condition and calculated using the effective collision oscillator strength: $$R_{0i}^{\text{coll\hspace{0.17em}}}=R_{i0}^{\text{coll\hspace{0.17em}}}\left(g_i/g_0\right)\exp (-E_{i0}/kT)\text{\hspace{0.17em}and}$$(B.6) $$R_{i0}^{\text{coll\hspace{0.17em}}}=8.63\times {10}^{-8}\left(f_{i0}^{\text{coll\hspace{0.17em}}}{g}_i^{-1}\right)n_{\mathrm{e}}\sqrt{{10}^4/T_{\mathrm{e}}}.$$(B.7)

The photoinduced rates are calculated using a given stellar SED F_st (λ) – the photon flux per unit area per unit time per unit wavelength interval at a given λ: $$R_{0i}^{\mathrm{p}\mathrm{h}}=\int {\sigma }_{0i}(\lambda )F_{\mathrm{s}\mathrm{t}}(\lambda )d\lambda .$$(B.8)

The cross section of absorption by a single atom is given by Lorentz profile: $${\sigma }_{ik}=f_{ik}\left(\frac{e^2}{m_{\mathrm{e}}{c}^2}\right)c\frac{\mathrm{\Delta }{\nu }_{ik}/2}{{(\nu -{\nu }_0)}^2+{\left(\mathrm{\Delta }{\nu }_{ik}/2\right)}^2}={\sigma }_{0,ik}\frac{{\nu }_0}{\pi }\frac{{\left(\mathrm{\Delta }{\nu }_{ik}/2\right)}^2}{{(\nu -{\nu }_0)}^2+{\left(\mathrm{\Delta }{\nu }_{ik}/2\right)}^2}\text{\hspace{0.17em}with}$$(B.9) $${\sigma }_{0,ik}=f_{ik}{\lambda }_{ik}\left(\frac{\pi e^2}{m_{\mathrm{e}}{c}^2}\right)=8.86\times {10}^{-21}{f}_{ik}{\lambda }_{ik}({\mathrm{A}}^{{}^{\circ }}){\mathrm{c}\mathrm{m}}^2=1.34\times {10}^{-36}{A}_{ik}{\lambda }_{ik}^3({\mathrm{A}}^{{}^{\circ }}){\mathrm{c}\mathrm{m}}^2,$$(B.10)

where $f_{ik}=1.5\times {10}^{-16}{\lambda }_{ik}^2\left(g_k/g_i\right)A_{ki}=1.5\times {10}^{-16}{\lambda }_{ik}^2{A}_{ik}$ is the oscillator strength (wavelength expressed in Angstroms) and A_ki – radiative decay rate from the upper level k.

Fig. B.1 The density profile of hydrogen atoms excited to the n = 2 level, obtained with the present model based on Eqs. (B.1–B.5) (red line) and with the Monte-Carlo code (blue dots) and the corresponding interpolation line. The green dashed line shows the level population given by LTE formula. The simulation parameters are those of the set N1 from Table 2.

Let us assume that the spectral flux of photons is broad compared to the line width. Then the spectral integral is given simply by $$R_{0i}^{\mathrm{p}\mathrm{h}}=F_{\lambda 0i}\int {\sigma }_{0i}d\lambda =F_{\lambda 0i}{f}_{0i}\left(\frac{e^2}{m_{\mathrm{e}}{c}^2}\right){\lambda }_{ik}^2\int \frac{\mathrm{\Delta }{\nu }_{0i}/2}{{(\nu -{\nu }_0)}^2+{\left(\mathrm{\Delta }{\nu }_{0i}/2\right)}^2}d\nu =F_{\lambda 0i}{\lambda }_{0i}{\sigma }_0.$$(B.11)

Here F_λ0i = F_st(λ_0i). Using CGS units, with wavelength in Angstroms, this gives: $$R_{0i}^{\mathrm{p}\mathrm{h}}\left(s^{-1}\right)=0.88\times {10}^{-20}{f}_{0i}{F}_{\lambda 0i}{\lambda }_{0i}^2({\mathrm{A}}^{{}^{\circ }})=1.32\times {10}^{-36}\left(g_i/g_0\right)R_{i0}^{\text{rad\hspace{0.17em}}}{\lambda }_{0i}^4({\mathrm{A}}^{{}^{\circ }})F_{\lambda 0i}\left({\mathrm{c}\mathrm{m}}^{-2}{\mathrm{s}}^{-1}{{\mathrm{A}}^{{}^{\circ }}}^{-1}\right).$$(B.12)

Thus, in an optically thin plasma, the radiative pumping and decay rates of the excited level can be directly compared for a given SED. This is shown in Fig. B.2. We note that Eq. (B.12) can be applied to transitions between any two levels.

Pumping of the excited state by resonant photons from stellar radiation must be calculated taking into account self-shielding. For example, calculating the radiation pressure of Lyα stellar photons under typical hot exoplanet conditions without accounting for the attenuation of the flux as it penetrates optically thick layers yields erroneous results, leading to an extremely excessive acceleration of atoms due to radiation pressure. In the case where the spectral features of the SED are smooth in the vicinity of the line under consideration, the flux attenuation can be calculated analytically. $$F_{\lambda 0i}=F_{\lambda 0i}\exp (-\int n_0{\sigma }_{0i}dl)\approx F_{\lambda 0i}\exp (-n_0{\sigma }_{0i}H)$$(B.13)

Fig. B.2 The relationship between pumping and emission rates $(R_{ki}^{\mathrm{p}\mathrm{h}}/R_{ik}^{\mathrm{r}\mathrm{a}\mathrm{d}})$, according Eq. (B.12) for a transition characterized by wavelength λ for the SED of KELT-9. A fixed value of gk/gi = 10 is assumed for the statistical weights.

Here, we adopted an exponential atmosphere with a scale height H and approximated the exponential power integral. The absorption cross section σ_0i(V, T) of an ensemble of atoms is a function of plasma parameters and wavelength, and it is given by the Voigt integral. We note that this differs from the situation when calculating the absorption rate of a single atom (B.9). In terms of Doppler velocity, instead of λ, the corresponding expression is $$R_{0i}^{\mathrm{p}\mathrm{h}}={\lambda }_{0i}{F}_{\lambda 0i}\frac{1}{c}{\displaystyle \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}}\exp [-n_{0}H{\sigma }_{0i}(V)]{\sigma }_{0i}dV.$$(B.14)

For the Doppler broadening, the integral can be expressed numerically as $$R_{0i}^{\mathrm{p}\mathrm{h}}={\lambda }_{0i}{F}_{\lambda 0i}\frac{1}{n_{0}H}\sqrt{\frac{2kT}{m_i{c}^2}}g\left(x_{max}\right)$$(B.15)

where $$g={\int }_0^{x_{max}}\frac{e^{-x}}{\sqrt{\ln (x_{max}/x)}}dx$$

with $x_{max}=n_{0}H{\sigma }_{0i}(V=0)=n_{0}H{\sigma }_{max},{\sigma }_{max}={\sigma }_{\text{abs\hspace{0.17em}}}(V=0)={\sigma }_0\sqrt{m_i{c}^2/2\pi kT}$ and σ₀ = f_0i λ_0i(πe²)/(m_ec²). Finally, $$R_{0i}^{\mathrm{p}\mathrm{h}}\approx F_{\lambda 0i}{\lambda }_{0i}{\sigma }_{0}g\left(x_{max}\right)/x_{max}$$(B.16)

where x_max = n₀Hσ_max, g(x) ≈ x exp (−0.58 x) at x < 2.62, and $g(x)\approx 1/\sqrt{\pi \ln x}$ at x ≥ 2.62.

Accounting for the trapping of resonant photons (B.5) and self-shielding (B.16) changes the excited-state population in a dense gas by orders of magnitude, and therefore cannot be ignored.

Appendix C Negative hydrogen ion

The H I⁻ion population is not in LTE, so we calculate it using the corresponding production and annihilation reactions, listed in the table below. In general, the primary role of H I⁻ is to reduce the opacity of the dense gas in the VUV and optical bands where the radiation absorption by other elements is insignificant. According to Fig. 2, photoionization of H I⁻ does not contribute significantly to the overall energy balance. Another important effect of the H I⁻ ion is its direct influence on the exited states of hydrogen: $${\mathrm{H}\mathrm{I}}^{-}+p\leftrightarrow \mathrm{H}\mathrm{I}(n3)+\mathrm{H}\mathrm{I}(n1).$$(C.1)

This is due to the reaction of neutralization of H I⁻ forming excited hydrogen atoms. This reaction has a relatively large cross section of ∼10⁻¹⁵ cm² (Janev et al. 1987) because, for n = 3, it proceeds with a very small energy change. The rate of the corresponding reverse reaction, which effectively destroys the excited H I by the most abundant particle (H I in the ground state), can be found using LTE relations. We note that even if the H I⁻population is not significant on its own, the reaction (C.1) can be important as a channel for the destruction of excited hydrogen.

Continuum absorption by H I⁻actually affects the photometric radius of the planet as inferred from observations. We are modeling the atmosphere of KELT-9 b, at a relatively high base pressure, to calculate the absorption more accurately, we must define the pressure corresponding to the photometric radius used in our analysis. This issue was analyzed by Fossati et al. (2021); Fossati et al. (2023b) for two similar planets orbiting similar stars: KELT-9 b and KELT-20 b. In both cases, according to their estimates, the photometric radius appears to match the pressure range of 0.1–0.01 bar. For the particular pressure profile of KELT-9 b, obtained in our model simulations based on P_base = 100 bar at R = R_p, this yields R_phm = 1.1 R_p.

Fig. C.1 Absorption profiles in the Hα line at mid-transit of KELT-9 b in units of Doppler velocity, are shown for different assumed photometric radii of the planet, disregarding the absorption by H I−(blue, violet, and turquoise lines). The red line shows similar calculations, but with the inclusion of absorption by H I−and assuming Rphm = Rp. Here and further on, semi-open circles with error bars represent the measurements (Yan & Henning 2018; Wyttenbach et al. 2020).

Figure C.1 shows how the absorption in the Hα line changes with R_phm and how it compares with continuum absorption by H I⁻. It can be seen that at R_phm = R_p and without H I⁻ the absorption line profile has wings caused by natural line broadening in the dense atmosphere, which are absent in the observations. These wings are suppressed by H I⁻ absorption. Thus, taking into account H I⁻ absorption with R_phm = R_p is approximately equivalent to assuming R_phm=1.1 R_p but without H I absorption.

Appendix D Analytical approximations

In the simplest case, considering only collisional excitation and de-excitation of the energy level, the system is described by Eqs. (B.1–B.5). This results in a LTE population of the excited level (beta factor β = 1) with zero net heating or cooling. The next step is to include the radiative pumping and decay of the level. In a stationary case, the corresponding relation between level populations is $$n_i=n_k\frac{R_{ki}^{\text{coll\hspace{0.17em}}}+R_{ki}^{\mathrm{p}\mathrm{h}}}{R_{ik}^{\text{coll\hspace{0.17em}}}+R_{ik}^{\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{a}}}.$$(D.1)

Therefore, as long as $R_{ki}^{\mathrm{p}\mathrm{h}}{R}_{ik}^{\text{coll\hspace{0.17em}}}<R_{ik}^{\mathrm{r}\mathrm{a}\mathrm{d}}{R}_{ki}^{\text{coll\hspace{0.17em}}}$ the population of the upper level is below the LTE value, generating the line cooling.

A relevant asymptotic case is a gas that is optically thick for transitions from the ground state, n₀ σ_max H ≫ 1, that $R_{i0}^{\text{rad,trap\hspace{0.17em}}}$ (B.5) becomes negligible compared to other rates. This also implies that the stellar pumping of the level is negligible (the so-called self-shielding). Under these conditions the level’s population is determined by collisions, recombination and, especially important for KELT-9 b, photoionization. The relevant equations are $$n_i=n_0\frac{R_{0i}^{\text{coll\hspace{0.17em}}}+R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}{n}_{\mathrm{e}}{\mathrm{n}}_0^{j+1}/n_0}{R_{i0}^{\text{coll\hspace{0.17em}}}+R_{if}^{\mathrm{p}\mathrm{h}}+R_{if}^{\text{coll\hspace{0.17em}}}}\text{\hspace{0.17em}and}$$(D.2) $$(R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}+R_{f0}^{\mathrm{r}\mathrm{e}\mathrm{c}})n_{\mathrm{e}}{n}_0^{j+1}=(R_{if}^{\mathrm{p}\mathrm{h}}+R_{if}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}})n_i+n_0(R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}+R_{0f}^{\mathrm{p}\mathrm{h}}).$$(D.3)

For the hydrogen atom, the photoionization from the ground state in (D.3) can be neglected, i.e., $$n_i/n_0=\frac{R_{0i}^{\text{coll\hspace{0.17em}}}+\alpha R_{0f}^{\text{coll\hspace{0.17em}}}}{R_{i0}^{\text{coll\hspace{0.17em}}}+(1-\alpha )(R_{if}^{\mathrm{p}\mathrm{h}}+R_{if}^{\text{coll\hspace{0.17em}}})},$$(D.4) $$n_{\mathrm{e}}{n}_0^{j+1}/n_0=\frac{1}{R_{fi}^{\text{rec\hspace{0.17em}}}+R_{f0}^{\text{rec\hspace{0.17em}}}}[R_{0f}^{\text{coll\hspace{0.17em}}}+\frac{(R_{0i}^{\text{coll\hspace{0.17em}}}+\alpha R_{0f}^{\text{coll\hspace{0.17em}}})(R_{if}^{\mathrm{p}\mathrm{h}}+R_{if}^{\text{coll\hspace{0.17em}}})}{R_{i0}^{\text{coll\hspace{0.17em}}}+(1-\alpha )(R_{if}^{\mathrm{p}\mathrm{h}}+R_{if}^{\text{coll\hspace{0.17em}}})}]\text{\hspace{0.17em}with}$$(D.5) $$\alpha =\frac{R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}}{R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}+R_{f0}^{\mathrm{r}\mathrm{e}\mathrm{c}}}.$$(D.6)

For KELT-9 b, the largest term in Eq. (D.4) is photoionization from the excited level. For type A recombination, when gas is not optically thick for ionizing photons (i.e., n₀ σ_0f H< 1), α ∼ 1, and 1 − α ∼ 1, the population is much smaller than in LTE, while the ionization degree is relatively small. $$\beta \approx R_{i0}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}/R_{if}^{\mathrm{p}\mathrm{h}}(1-\alpha {)}^{-1}\ll 1$$ $$n_{\mathrm{e}}{n}_0^{j+1}/n_0=R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}/(R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}+R_{f0}^{\mathrm{r}\mathrm{e}\mathrm{c}})(1-\alpha {)}^{-1}$$(D.7)

For type B recombination, when recombination directly to the ground state is blocked, and 1 − α ≪ 1, the population of the level approaches LTE. Simultaneously, recombination into the excited level is balanced by photoionization from this level, and the ionization degree reaches a relatively high value. $$\beta =1$$ $$n_{\mathrm{e}}{n}_0^{j+1}/n_0=\left(R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}/R_{fi}^{\mathrm{r}\mathrm{e}\mathrm{c}}\right)\left(R_{if}^{\mathrm{p}\mathrm{h}}/R_{i0}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}\right)$$(D.8)

The heating or cooling, W, is produced by electron ionization from the ground state, excitation-de-excitation, and photoionization from the excited level. It also includes the thermal energy of the newly produced electrons. We note that at β ≤ 1 the cooling due to electron ionization from the excited state can be ignored compared to ionization from the ground state. $$W/n_0=E_{i0}(\beta -1)R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}-(E_{\text{ion\hspace{0.17em}}}+1.5k{T}_{\mathrm{e}})R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}++({\overline{E}}_{0f}^{\mathrm{p}\mathrm{h}}-1.5k{T}_{\mathrm{e}})R_{0f}^{\mathrm{p}\mathrm{h}}+({\overline{E}}_{if}^{\mathrm{p}\mathrm{h}}-1.5k{T}_{\mathrm{e}})R_{if}^{\mathrm{p}\mathrm{h}}\beta R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}/R_{i0}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}$$(D.9)

To summarize, at type A recombination, $$W/n_0\approx ({\overline{E}}_{0f}^{\mathrm{p}\mathrm{h}}-1.5k{T}_{\mathrm{e}})R_{0f}^{\mathrm{p}\mathrm{h}}+[({\overline{E}}_{if}^{\mathrm{p}\mathrm{h}}-1.5kT)/(1-\alpha )-E_{i0}]R_{0i}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}$$(D.10)

and at type B recombination: $$W/n_0=-(E_{\text{ion\hspace{0.17em}}}+1.5k{T}_{\mathrm{e}})R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}+({\overline{E}}_{if}^{\mathrm{p}\mathrm{h}}-1.5k{T}_{\mathrm{e}})R_{if}^{\mathrm{p}\mathrm{h}}\left(g_i/g_0\right)\exp (-E_i/k{T}_{\mathrm{e}}).$$(D.11)

We note that the last term in Eq. (D.10) is usually negative, so only photoionization from the ground state and related heating can balance the cooling. This can be achieved at relatively low temperatures when $R_{0i}^{\text{coll\hspace{0.17em}}}\sim \exp (-E_i/k{T}_e)$ is sufficiently small.

At type B recombination: $$W/n_0=-(E_{\text{ion\hspace{0.17em}}}+1.5k{T}_{\mathrm{e}})R_{0f}^{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}+({\overline{E}}_{if}^{\mathrm{p}\mathrm{h}}-1.5k{T}_{\mathrm{e}})R_{if}^{\mathrm{p}\mathrm{h}}\left(g_i/g_0\right)\exp (-E_i/k{T}_{\mathrm{e}}).$$(D.12)

The first term in Eq. (D.12) is small, and the balance is achieved at a relatively high temperature approximately 1.5 kT_e ≈ Ē_ph.

References

All Tables

All Figures

Fig. 1 Mean energy of photoelectrons (asterisks, in eV) calculated with Eq. (4) and a corresponding thermal equilibrium temperature of 1.5 kTe = Ēel, ph (dashed line) versus the threshold of photoionization calculated for the SED of KELT-9. For the cross section, a simple formula is used: σph-ion = σthr(Ethr/E)3, E > Ethr. The labeled vertical lines indicate the photoionization thresholds of several particular elements, or transitions, of interest.

In the text

Fig. 2 Profiles of various quantities of the H–He atmosphere of KELT-9 b along the y axis directed opposite to the orbital motion of the planet, obtained in simulation with the parameter set N 1 (Table 2 below) and expressed via pressure. Top panel: Population of H I(n2) level in terms of LTE value (red, left axis); ionization degree Ne/Na (orange) in comparison to Saha equilibrium (turquoise) and type B recombination (blue, formula (D.8)); temperature (gray, upper horizontal axis). Bottom panel: Heating (>0, solid lines) or cooling (<0, dashed lines) terms in the energy equation associated with H I(n1) and H I(n2) levels. Photoionization heating from the H I(n2) level given by the analytical formula $W=n_{{\mathrm{H}\mathrm{I}}_{\mathrm{I}}}(n2)R_{2f}^{\mathrm{p}\mathrm{h}}{\overline{E}}_{2f}^{\mathrm{p}\mathrm{h}}$ is shown with a green line.

In the text

Fig. 3 Cooling rates for Mg and O calculated with a unit particle density (i.e., 1 cm−3) for all species, including electrons, and with zero radiation field. Black lines show cooling due to the excited levels listed in the table, red lines show full cooling including ionization from the ground state, and, for comparison, yellow lines show the analogous simulation result with the code CLOUDY from Gnat & Ferland (2012).

In the text

	Fig. 4 Cooling or heating rates of elements, calculated for KELT-9 b with an account of the lines listed in Table 1 at ne = 1011 cm−3 under the impact of an unshielded radiation flux of the host star. Depending on temperature, the rates are either negative (i.e., cooling, dashed lines), or positive (i.e., heating, solid lines). The vertical axis is in log scale.
In the text

Fig. 5 Top panel: distribution of the proton density in the equatorial plane of the KELT-9 stellar system (i.e., in the whole simulation domain) simulated with the fiducial parameter set N5. Black circle indicates the star in scale. The planet is located at the coordinate center. Streamlines are shown originating at the star (black lines), and at the planet (blue lines). Bottom panel: distribution of the excited hydrogen density in the meridional plane in close vicinity to KELT-9 b simulated with the parameter set N5. The scale of plots is in units of Rp. The values outside the specified variation ranges are colored either in red, if smaller than minimum, or in blue, if higher than maximum.

In the text

	Fig. 6 Density profiles of different elements along the planet-star line (x > 0 spatial axis, expressed in units of pressure). The gray line with the upper horizontal axis shows the temperature (top panel) and the height above the planet h = R/Rp − 1 (bottom panel). Results are obtained with parameter set N5.
In the text

Fig. 7 Population profiles of different levels along the planet-star line at the day side (x > 0) and at the night side (x < 0). The left panel shows the hydrogen densities calculated by model (solid) lines and the corresponding LTE values (dashed). The center and right panels show the levels’ profiles as beta-factors (i.e., density values normalized to LTE values) for Mg I, Mg II, and O I, C I, C II.

In the text

Fig. 8 Heating (solid lines) or cooling (dashed lines) due to various elements and levels along the planet-star line (x > 0, expressed as a pressure). The upper panel shows the main considered elements of interest. The ground and excited levels of hydrogen are distinguished, whereas for other elements, all simulated levels are combined. The “sum” line indicates the accuracy of the overall balance between heating and cooling achieved in simulation. The gray line shows the temperature profile. The center and bottom panels show the heating or cooling contribution of excited magnesium (top), oxygen, and carbon (bottom), respectively.

In the text

Fig. 9 Absorption profiles in Hα line at mid-transit in units of Doppler velocity obtained in test simulations: when the stellar flux in VUV range is put to zero (red); when the residual energy of newly created electrons is put to zero (instead of 1.5kTe, green); when the recombination heating is put to zero (blue). For comparison, the black line shows the Hα absorption profile calculated with the fiducial set of the model parameters N5.

In the text

	Fig. 10 Absorption profiles in Hα, Paβ, and O I 777.4 nm lines at mid-transit in units of Doppler velocity for the simulations with different model parameter sets specified in Table 2.
In the text

	Fig. 11 Absorption profiles in Hα, Paβ, and O I 777.4 nm lines calculated with different combinations of VUV flux and He/H abundance ratio, which enable a best fit of the observations.
In the text

Fig. 12 Absorption in Hα line averaged in the velocity interval ± 20 km s−1 around the line center versus the integral mass loss of atmosphere. Different points were obtained by independently varying the He/H abundance ratio within the range [0.03; 0.1], VUV integral flux within [X0.5; X2], and base temperature Tb within [4000; 5000] K. The gray circle highlights the points that best correspond to measurements.

In the text

	Fig. 13 Total heating (red, bottom axis) and temperature (black, top axis) profiles obtained in the present work (solid lines), in Fossati et al. (2021, dotted lines), and in García Muñoz & Schneider (2019, dashed lines).
In the text

Fig. 14 Absorption at the center of the Hα line (at V = 0, black circles) and at the distant wing (V = 80 km s−1, white circles), manually restricted by calculation within region Rp < r < Rmax around the planet, in dependence on Rmax. The vertical dashed lines show the region of formation of the corresponding absorption at V = 0 at the level of 0.5 to 0.9 of the maximum value. The gray curve shows the pressure in the atmosphere at the corresponding altitudes.

In the text

	Fig. 15 Absorption profiles in Hβ 486.1 nm line at the same model runs as those in Fig. 11 (N5 – red, N12 – olive, N13 – blue). The gray circles with error bars show the average measurements over the available observations (from Fossati et al. 2020).
In the text

	Fig. B.1 The density profile of hydrogen atoms excited to the n = 2 level, obtained with the present model based on Eqs. (B.1–B.5) (red line) and with the Monte-Carlo code (blue dots) and the corresponding interpolation line. The green dashed line shows the level population given by LTE formula. The simulation parameters are those of the set N1 from Table 2.
In the text

	Fig. B.2 The relationship between pumping and emission rates $(R_{ki}^{\mathrm{p}\mathrm{h}}/R_{ik}^{\mathrm{r}\mathrm{a}\mathrm{d}})$, according Eq. (B.12) for a transition characterized by wavelength λ for the SED of KELT-9. A fixed value of gk/gi = 10 is assumed for the statistical weights.
In the text

Fig. C.1 Absorption profiles in the Hα line at mid-transit of KELT-9 b in units of Doppler velocity, are shown for different assumed photometric radii of the planet, disregarding the absorption by H I−(blue, violet, and turquoise lines). The red line shows similar calculations, but with the inclusion of absorption by H I−and assuming Rphm = Rp. Here and further on, semi-open circles with error bars represent the measurements (Yan & Henning 2018; Wyttenbach et al. 2020).

In the text