Global 3D Simulation of the Upper Atmosphere of HD189733b and Absorption in Metastable He i and Lyα Lines

Gas giants are numerous and are the most studied among the known extrasolar planets, because they produce strong, easily detectable transit signals. These so-called 'hot Jupiters' are exposed to an intense X-ray and ultraviolet (XUV) radiation coming from the host stars, and the planetary upper atmospheres are subject to strong radiative heating and ionization. This results in the expansion and hydrodynamic escape of the upper atmosphere and thus mass loss (Lammer et al. 2003). Furthermore, the presence of a transit gives one the possibility of detecting absorption by the upper atmosphere, probing the size, composition, temperature, and velocity of the surrounding media.

The first exoplanet atmosphere detected with the transit method was HD209458b (Charbonneau et al. 2000). Its transits were also observed in the Lyα line (Vidal-Madjar et al. 2003) and revealed absorption several times higher than that in optical photometry (Ben-Jaffel 2007; Vidal-Madjar et al. 2008). Later, Vidal-Madjar et al. (2004), Linsky et al. (2010), and Cubillos et al. (2020) measured the comparable transit absorption depths of 6%–9% for this exoplanet in the C ii, O i, Si iii, and Fe ii lines. These observations of HD209458b were followed by similar spectral measurements for HD189733b, which also transits a bright star. Due to the greater planet/star radii ratio, it is even more suitable for transit measurements. However, in other aspects, HD189733b is quite different from HD209458b. Most notably, it is smaller, but more massive, and it orbits significantly closer to its active K1.5V host star (Bouchy et al. 2005).

Lecavelier des Etangs et al. (2010) and Lecavelier des Etangs et al. (2012) reported one spectroscopically unresolved absorption in the Lyα line at the level of 5.05% ± 0.75% and one resolved absorption line with almost the same value of the total absorption of 5.0% ± 1.3% over the whole line, but with an essentially enhanced absorption of 14.4% ± 3.6% in the range of blueshifted velocities [−230, −140] km s⁻¹. Lecavelier des Etangs et al. (2012) pointed out that this effect might be related to the upper atmosphere of HD189733b, which overflows the Roche lobe and supersonically streams away. There was also one spectroscopically resolved measurement, which yielded a very small excess absorption in the Lyα line (Lecavelier des Etangs et al. 2012). In addition to that, the far ultraviolet (FUV) absorption lines of ionized C, N, and Si have shown significant short-time variability, most probably caused by flare activity of the host star (Pillitteri et al. 2015). This makes the analysis and interpretation of the transit absorption by this exoplanet in FUV lines particularly difficult and ambiguous.

An important feature of HD189733b consists of a severe day-night heating contrast. The simulations with a 3D global circulation model (GCM) predicted the zonal winds with velocities of [−5, 10] km s⁻¹ in the lower layers of the atmosphere (Showman & Polvani 2011). The upper atmosphere escape on HD189733b was modeled with 1D hydrodynamic (HD) codes. In particular, for the ionizing (λ < 912 Å) XUV radiation flux F_XUV = 2 × 10⁴ erg cm⁻² s⁻¹ at the planet orbit, Guo (2011) obtained the planetary mass-loss rate of 5 × 10¹⁰ g s⁻¹, the upper atmosphere maximum temperature slightly above 10⁴ K, and the location of half-ionization H/H⁺ = 1 level at ∼1.1R_p, indicating a very fast photoionization of the atmospheric hydrogen. Menager et al. (2013), using an aeronomy code, found that the H/H⁺ = 1 level appears at ∼1.5R_p. The most detailed and recent 1D HD simulation by Salz et al. (2016) predicts for F_XUV = 1.8 × 10⁴ erg cm⁻² s⁻¹ a somewhat lower mass-loss rate of about (1−2) × 10¹⁰ g s⁻¹, the maximum temperature of 1.2 × 10⁴ K, and the location of the H/H⁺ = 1 level at 1.2R_p. A reason for the lower mass-loss rate was found to be the cooling effect, caused by the excitation of hydrogen and emission of Lyα photons.

There are only a few works dedicated to the interpretation of the Lyα observations of the HD189733b. Bourrier & Lecavelier des Etangs (2013) used a 3D Monte Carlo code, tracing the hydrogen atoms launched from an inner spherical boundary surrounding the planet at 3R_p. They found that acceleration of the atoms by radiation pressure alone is not sufficient for the generation of a population of energetic neutral atoms (ENAs) with velocities of up to 200 km s⁻¹, and that an additional mechanism should operate, namely charge exchange between the atmospheric atomic hydrogen and stellar protons. In Odert et al. (2020), a 1D HD model of the planetary wind was combined with a 3D MHD model of the magnetized stellar wind (SW). For the XUV flux of 1.8 × 10⁴ erg cm⁻² s⁻¹, they obtained the upper atmosphere maximum temperature of ≈1.1 × 10⁴ K, the mass-loss rate of 5 × 10¹⁰ g s⁻¹, and the location of the H/H⁺ = 1 level at ≈1.7R_p. Odert et al. (2020) also found that the absorption in the Lyα line is mostly due to the natural line broadening mechanism, and is of order ∼11% in the blue wing of the line [−400, −40] km s⁻¹, and 13% in the red wing [40, 400] km s⁻¹. Both numbers are higher than the observed value. Odert et al. (2020) have shown that even an extremely strong SW would not generate enough ENAs to produce large blueshifted absorption in excess of 140 km s⁻¹ as reported in Lecavelier des Etangs et al. (2010, 2012). The parameters of the considered strongest SW corresponded to the conditions in coronal mass ejections (CMEs) with the velocity of 1000 km s⁻¹ and the integral stellar mass loss 10³ times higher than the average solar value (≈2.5 × 10¹² g s⁻¹). Under these conditions, the magnetopause was pushed by the SW ram pressure toward the planet as close as 2.3R_p.

In 2018 a new opportunity for probing the expanding and escaping planetary atmospheres was opened (Seager & Sasselov 2000; Oklopčić & Hirata 2018), based on the measurements of the transit absorption at the position of the metastable helium 2³S triplet line at 10830 Å (hereafter He i(2³S)). The metastable 2³S level of the helium atom has a high excitation energy of 19.8 eV and a radiative lifetime of about 2.2 hr. It is efficiently pumped by the recombination of ionized helium and results in the absorption in the infrared band due to the triplet transition 2³S−2³P at 10830.34 Å, 10830.25 Å, and 10829.09 Å. The He i(2³S) 10830 Å triplet line has been used to estimate the magnetic field in solar filaments (Lin et al. 1998) and to determine the helium abundance in extragalactic H ii regions (Aver et al. 2015), as well as to probe the material flows in active galactic nuclei (Leighly et al. 2011). Seager & Sasselov (2000) proposed using this line to observe the transits of HD209458b. Determination of the population of the He i(2³S) state in the upper atmospheres of hot exoplanets is a rather complex problem tackled in Oklopčić & Hirata (2018) by means of a 1D HD model. They predicted the absorption at the line center at a level of 8% and 2% for GJ436b and HD209458b, respectively.

Recent infrared observations of exoplanetary transits with ground telescopes, mostly with the CARMENES instrument at the Calar Alto Observatory and with NIRSPEC on Keck II, provided a series of valuable data for the analysis and interpretation. The He i(2³S) absorption depths of ∼8% for Wasp107b (Allart et al. 2019; Kirk et al. 2020) and ∼1% for HATP11b (Allart et al. 2018; Mansfield et al. 2018), GJ3470b (Ninan et al. 2020, Pallé et al. 2020), and HD209458b (Alonso-Floriano et al. 2019), as well as ∼3% for Wasp69b (Nortmann et al. 2018), were measured. Most of the absorption profiles have a typical spectral width of ±10 km s⁻¹ in Doppler velocity units. The reproducibility of the experimental data, the typical signal-to-noise ratio, and the statistical significance ≥10σ noticeably exceed those achieved for the measurements in the FUV band. As to the hot Jupiter HD189733b, considered in the present paper, there are two independent observations described by Salz et al. (2018) and Guilluy et al. (2020). They are similar in the half-width of the He i(2³S) absorption profile, but differ in depths by a factor of 1.5 (i.e., 1% and 0.7%, respectively). There is also a net blueshift of the whole absorption profile by ≈2.5 km s⁻¹ reported in Salz et al. (2018), but absent in the Guilluy et al. (2020) data. We note that signatures of a small blueshift are also reported for other planets. However, Salz et al. (2018) remarks that these blueshift velocities could be potentially caused by the stellar pseudo-signals, i.e., due to stellar surface regions with stronger He i absorption. To simulate the absorption profiles at 10830 Å for different exoplanets a number of 1D HD models were used. These 1D models calculate temperature, velocity, and electron, hydrogen, and helium density distribution profiles in the upper atmosphere, as well as the kinetics of the He i(2³S) state population. Ninan et al. (2019) and Palle et al. (2020) found for GJ3470b that the area in which He i(2³S) is present extends up to at least 10R_p, and that the observed transit depths can be fitted in the modeling with a relatively low helium abundance of about He/H ∼ 0.01. The most sophisticated simulation of the He i(2³S) absorption, based on 1D HD modeling, is described by Lampón et al. (2021). Because we explored GJ3470b in our previous work (Shaikhislamov et al. 2021) and consider HD189733b in this paper, we discuss the approach of Lampón et al. (2021) in more detail. It is based on a Parker-like solution, where instead of temperature, the sound speed is assumed to be constant. In addition to temperature, other fitting parameters are the mass-loss rate and the helium abundance. The fixed XUV flux of 5.7 × 10⁴ erg cm⁻² s⁻¹ is used in the model to calculate the ionization balance and the kinetics of the He i(2³S) state. To fit the measured absorption profiles, Lampón et al. (2021) used up to three velocities (that is, line-of-sight components) with different coverage of stellar disk for HD189733b (altogether, six empirical fitting parameters), whereas for GJ3470b only the net blueshift appeared to be sufficient. However, the fitting of the He i(2³S) absorption profiles was strongly dependent on the temperature and mass loss considered as the modeling free parameters, but was equally good for the broad range of helium abundances, from 1% to 10%. To constrain this parametric degeneracy, more self-consistent modeling approaches have to be used by Lampón et al. (2021). They used the temperature and the mass loss derived from the best fit of the Parker-like solution profiles to the profiles from Salz et al. (2016, 2018) for GJ3470b, and from Salz et al. (2016) and Odert at al. (2020) for HD189733b. It is argued in Lampón et al. (2021) that to derive the He/H ratio, the parameters of planetary outflow should be constrained by the Lyα measurements. At the same time, the detected Lyα absorption of GJ3470b is significantly blueshifted (35% ± 7% in the range [−94, −41] km s⁻¹; Bourrier et al. 2018) and it cannot be generated just by the planetary atmosphere outflow. Therefore, as in the case of GJ3470b, the contribution of SW and the generation of ENAs are essential. These factors, however, cannot be properly included in the 1D model. Nevertheless, to give at least a hint of them, Lampón et al. (2021) performed one of the fits of the He(2³S) absorption profile with the assumption that the escaping planetary material beyond the Roche lobe moves away from the star at a constant prescribed velocity.

Thus, combining different 1D models still requires, in order to address the essentially 3D processes and observational features, the artificial addition of nonspherically symmetric effects. Only a global 3D multicomponent model can simulate such effects and predict realistic density distributions of the major species and the structure of the interacting planetary and SW material flows. This is especially true for such objects as HD209458b and HD189733b, for which the spectral absorption measurements of lines of different species have been made, e.g., H i (Lyα), O i, C ii, and He i(2³S). Ideally, the interpretation of these measurements should be based on the simulation of dynamical behavior of all species and respective lines within the frame of a global self-consistent model. In the present paper we make a step forward in this direction and perform a multifluid 3D HD self-consistent modeling of the escaping upper atmosphere of HD189733b surrounded by the stellar wind, simulated with the same model, and calculate the synthetic transit absorption features. In addition to presenting the simulation and the interpretation of observations of HD189733b, we compare our results with those of other modeling approaches.

The model we use in our study here has been already applied for the simulation and interpretation of observations of GJ436b (Khodachenko et al. 2019), π Men c (Shaikhislamov et al. 2020a), HD209458b (Shaikhislamov et al. 2020c; Khodachenko et al. 2021b), GJ3470b (Shaikhislamov et al. 2020b), and Wasp107b (Khodachenko et al. 2021a). The last three most recent papers also modeled the absorption by metastable helium at 10830 Å. The interpretation of observations of GJ3470b in both the Lyα and He i(2³S) lines revealed that the bowshock layer formed by SW around the planet plays a significant role. They found the best-fit parameters for both lines to be F_XUV = 8 erg cm⁻² s⁻¹, He/H = 0.013, and the SW with a density and velocity similar to the fast solar wind, strong enough to form a bowshock at about 20R_p (Shaikhislamov et al. 2021).

The simulation of absorption at 10830 Å during the transits of Wasp107b revealed another feature that is beyond the capacity of 1D models; it concerns the radiation pressure, acting on metastable He i(2³S) atoms. As pointed out in Allart et al. (2019), the acceleration of helium atoms by infrared (10830 Å) radiation flux of the star Wasp107 is ∼75 times higher than that due to the stellar gravity. At the same time, as shown in Khodachenko et al. (2021a) the observed He i(2³S) absorption at Wasp107b can be well reproduced with the typical expected parameters of the system and by taking account of the electron and atom impact processes, without invoking of any specific additional assumptions of the stellar radiation and helium abundance.

Regarding the results reported in the present paper, we note that the 3D aeronomy simulations and global HD modeling are performed for the first time in the study of HD189733b. It is also the first attempt to interpret simultaneously both Lyα and He i(2³S) measurements at transits of HD189733b within a single model. In this respect it is worth mentioning that Odert et al. (2020), aiming at the interpretation of Lyα observations during the transits of HD189733b, considered a more sophisticated magnetized SW case in 3D, but they took the escaping planetary atmospheric flow from the 1D simulations and did not take helium into account at all. In the present paper, we investigate, as in Odert et al. (2020), whether a strong SW can generate a sufficient amount of ENAs to affect the Lyα blue wing absorption and to explain the significant difference between the measured Lyα absorption profiles in different observations. Based on the performed simulations and fitting to the measured He i(2³S) absorption profile, we also constrain the helium abundance value in the upper atmosphere of HD189733b and show that the differences in observations can be explained by the temporal variations of the stellar XUV flux.

The paper is organized as follows. Section 2 contains a short description of the model and provides information about initial parameters. Section 3 presents the obtained results, which are discussed in Section 4.

The 3D multifluid global model used in this paper has been developed since Shaikhislamov et al. (2014). Its most important parts are described in the following papers: the details of gas dynamic equations and terms related to the generation of a planetary wind, the collisional interaction of species, and the applicability of a hydrodynamic approach are given in Shaikhislamov et al. (2016) and in the Appendix of Dwivedi et al. (2019). The details of the 3D code and of the stellar plasma wind simulation are given in Khodachenko et al. (2019); procedures of calculating transit absorption in lines are given in Shaikhislamov et al. (2018b); and the simulation of the metastable helium population and absorption are given in Shaikhislamov et al. (2021). Here, we only briefly outline the model. It solves numerically for the system of mass, momentum, and energy equations in spherical geometry, taking into account the gravity and Coriolis forces, as well as the radiation pressure of the stellar Lyα and IR flux acting on H i and He i(2³S) atoms, respectively. The model includes also the processes of hydrogen and helium photochemistry, in particular the reactions of photoionizaton, recombination, and electron impact. For the upper atmosphere of HD189733b, we neglect the molecular species that are present in the lower atmosphere, which are quickly dissociated with height. Thus, in this paper we deal with the following species: H, H⁺, He, He⁺, He²⁺. An additional argument in support of the atomic upper atmosphere assumption comes from the modeling of He i(2³S) absorption discussed later in the paper. However, in several dedicated model runs we studied the effects of molecules and included H₂, H₂ ⁺, and H₃ ⁺.

The metastable He i(2³S) atoms are treated in the model as a separate fluid with its own velocity and temperature, determined by those of the species from which it originates, i.e., He⁺ in the case of recombination, or He i(1¹S), in the case of excitation from the ground state. The reactions involved in the processes of population and depopulation of the He i(2³S) triplet level are similar to those considered in Oklopčić & Hirata (2018) and described in Shaikhislamov et al. (2021).

The code uses an explicit numerical scheme of second-order spatial and temporal accuracy, realized on a spherical grid in the planet-centered reference frame with the polar axis, Z, perpendicular to the orbital plane. More details about the model implementation can be found in Shaikhislamov et al. (2018a, 2020a, 2021) and in Khodachenko et al. (2019) and Khodachenko et al. (2021). At the inner boundary of the simulation domain, i.e., the planet surface at r = R_p, we fix the radial velocity of the atmospheric material to zero, the base atmosphere temperature to 1000 K and the pressure to 0.05 bar. The photoionization and heating of gas by photoelectrons, as well as the radiation attenuation in the planetary atmosphere, are calculated over the XUV and near-ultraviolet (NUV) spectra binned by 1 Å according to wavelength-dependent cross sections of various elements involved in the simulation. For the particular spectra we employ the rescaled spectral energy distribution (SED) of epsilon Eri to emulate the emission of HD189733. This choice is motivated by the fact that the two stars have the same spectral type and comparable activity levels (see, for example, Chadney et al. 2017). In the XUV band (10 < λ < 912 Å) the radiation flux is F_XUV ≈ 25 erg cm⁻² s⁻¹ at the reference distance of 1 au. The latest reanalysis of available X-ray and FUV observations and modeling of synthetic extreme ultraviolet (EUV) spectra of HD189733 (Bourrier et al. 2020) produced stellar XUV fluxes varying from visit to visit of 23 to 31 erg cm⁻² s⁻¹ at 1 au. At the same time, we distinguish in our model a hard energy part of the flux (XUVH) at λ < 510 Å with F_XUVH = 16 erg cm⁻² s⁻¹ at 1 au, which influences the metastable helium population. The ratio of F_XUVH/F_XUV = 2/3 was kept unchanged in course of the modeling (except for specially addressed cases), whereas the value of the total integral XUV flux F_XUV was varied for different model runs. The employed SED is characterized in the NUV (912 < λ < 3000 Å, which affects the He i(2³S) population by its photoionization) by F_NUV = 1830 erg cm⁻² s⁻¹ at 1 au, while in the near-infrared part (NIR, around 10830 Å, which is used to calculate the radiation pressure force acting on He i(2³S)) it is F₁₀₈₃₀ =20 erg cm⁻² s⁻¹ Å ⁻¹ at 1 au.

The absorption is calculated post-simulation by a special procedure described in Khodachenko et al. (2017) and Shaikhislamov et al. (2018b) for Lyα line and in Shaikhislamov et al. (2021) for the He i 2³P−2³S transitions. For the latter, the synthetic spectral profiles are obtained such as those for the observational data of Salz et al. (2018) and Guilluy et al. (2020); that is, by averaging the absorption between the second and the third contact transit points in a planet reference frame.

3.1. The Role of Lyα Cooling of the Upper Atmosphere

The cooling of the expanding upper atmospheric material due to the excitation of hydrogen atoms and the emission of Lyα photons is an important effect in the partially ionized upper atmospheres of hot Jupiters. However, its quantitative description is nontrivial because of the trapping and diffusion of Lyα photons in dense regions of the planetary atmosphere. The first aeronomy codes (Yelle 2004; Koskinen et al. 2007; García Muñoz 2007) ignored this Lyα cooling effect, assuming the eventual loss of Lyα photons to photoionization of the excited states of trace elements.

However, Murray-Clay et al. (2009) simulated the escaping planetary wind accounting for Lyα cooling and have shown that it changes the energy budget of a typical hot Jupiter by of order 10%. Later, Shaikhislamov et al. (2014), using a simple model of the diffusion of Lyα photons, have shown that the collisional electron de-excitation of the excited hydrogen 2p states reduces the effect of Lyα cooling and influences the temperature profile of the upper atmospheric material near the inner boundary of the simulation domain at low altitudes. For sufficiently massive planets the increasing effect of Lyα cooling suppresses the expansion of the planetary atmosphere, until the regime of the material dynamical escape is replaced by the quasi-static radiatively cooled atmosphere (Weber et al. 2018). Salz et al. (2016) obtained for HD189733b exactly this regime of a quasi-stationary upper atmosphere with a low mass-loss rate of (1−2) × 10¹⁰ g s⁻¹.

At the same time, there exist much more sophisticated simulations of the stationary upper atmosphere of HD189733b based on direct Monte Carlo modeling (Christie et al. 2013 and Huang et al. 2017). This approach takes into account a number of processes that terminate the Lyα scattering cycle. Among them are photoionization from the n = 2 state, collisional de-excitation, two photon decay, and photoionization of metals, e.g., Na and Mg. Altogether, such Monte Carlo simulations led to the conclusion that the major cooling processes involve metals and other trace elements, while the contribution of Lyα cooling does not exceed 10%, as compared with the total XUV heating.

Therefore, so far there is no model that combines all necessary physical processes needed for the detailed and self-consistent simulation of Lyα cooling of the expanding upper atmosphere of a hot Jupiter. For a reason that will be specified later, we consider in this work the effectiveness of Lyα cooling depending on an empirical parameter defined via a so-called density cutoff in the following way. We reduce the Lyα cooling rate by a scaling factor exp(−σ H n_H/α), where n_H is the atomic hydrogen density, σ ≈ 6 × 10⁻¹⁴ cm² is the cross section of Lyα photon absorption at a typical temperature of 10⁴ K, and H is the barometric scale height at a given gas temperature. The empirical parameter α defines the number of scatterings after which the Lyα cycle is terminated. Thus, the typical value n_cut = α/σ H corresponds to the density layer below which Lyα cooling rapidly decreases in the model.

We found out that in the case of Lyα cooling acting over all altitudes the hydrodynamic escape in our model gets very low and the upper atmosphere approaches radiative hydrostatic equilibrium. In this regime, the total rates of Lyα-emitted and XUV-absorbed energy are related as W_Lyα /W_XUV = 0.76, with W_XUV ≈ 4 × 10²⁴ erg s⁻¹. The corresponding mass-loss rate ≈ 1.5 × 10¹⁰ g s⁻¹, the maximum temperature ≈10⁴ K, and the velocity at 2R_p is about 1.5 km s⁻¹. All of these values are rather close to those obtained in Salz et al. (2016), despite the differences between the 1D and 3D models. However, without accounting for Lyα cooling the mass-loss rate becomes an order of magnitude higher, 2.2 × 10¹¹ g s⁻¹, with the temperature maximum of 1.5 × 10⁴ K, and the velocity at 2R_p of about 6 km s⁻¹. The differences between these two cases are shown in Figure 1. Altogether, the performed comparison reveals that without Lyα cooling the mass-loss rate and the velocity of material escape are higher, mostly because of the higher temperature.

Zoom In Zoom Out Reset image size

The most important point here is that the regimes with and without Lyα cooling are very different with respect to the He i(2³S) absorption, which prompted us to reconsider the problem of Lyα cooling. Under the conditions of full Lyα cooling, the velocity of escaping material (PW) is low, and the width of the corresponding He i(2³S) absorption profile appears to be too small to fit the observations, whereas without Lyα cooling, the absorption width agrees with the measurements. Figure 2 (left panel) shows the simulated absorption profiles at the position around 10830 Å for different levels of Lyα cooling (including full- and no-cooling cases), obtained with different values of the above-defined cutoff density, which controls the effective trapping of Lyα photons. It shows that to make the shape of the synthetic He i(2³S) absorption profiles closer to the measurements, we should assume a significant reduction of Lyα cooling, with the cutoff density of at least n_cut =3 × 10⁶ cm⁻³, which corresponds to the trapping parameter σ H n_H > 10². In this case, the energy rates are related as W_Lyα /W_XUV = 0.37. According to our approach, the case of practically no Lyα cooling corresponds to n_cut = 3 × 10⁵ cm⁻³, whereas for full Lyα cooling, n_cut = 10¹¹ cm⁻³.

Zoom In Zoom Out Reset image size

It should be noted that our model does not show the increased absorption at the weakest component of the triplet between [−50, −20] km s⁻¹, nor the small net blueshift observed by Salz et al. (2018; though it is absent in the observation of Guilluy et al. 2020). These features we discuss later.

With a series of dedicated simulation runs we also checked another possibility that can produce the broadening of the He i(2³S) absorption profile, namely the zonal wind, whose speed for HD189733b might reach up to 10 km s⁻¹ at altitudes with pressures of ≤30 mbar (Showman & Polvani 2011). For that we emulated the equatorial jet as a prescribed additional rotation of the planet by imposing the corresponding azimuthal velocity V_ϕ = R_p × Ω × sinθ at the inner boundary of the simulation domain at r = R_p. Here, Ω is the cyclic frequency of the prescribed rotation. The right panel of Figure 2 shows the simulated absorption profiles obtained in the case of full Lyα cooling (n_cut = 3 × 10¹⁶ cm⁻³) for three different values of the azimuthal velocity V_ϕ = {19, 11, 5.5} km s⁻¹ at the equator, which correspond to the additional rotation periods P = 2π/Ω = {8, 12, 24} hr. In particular, it can be seen that the close-to-observations width of the absorption profile is achieved at a rather high equatorial velocity in the range V_ϕ = (10–20) km s⁻¹.

Based on this preliminary study, further on we consider only the scenario of the reduced Lyα cooling with the cutoff density n_cut = 3 × 10⁶ cm⁻³. This value allows us to obtain the most plausible widths for the spectral profiles that fit the observations of the He i(2³S) line.

3.2. Constraining the XUV Flux and Helium Abundance

At first, we consider the case of a weak SW with the following parameters at the planetary orbit: V_sw = 210 km s⁻¹, T_sw = 0.7 MK, and n_sw = 10³ cm⁻³. It is worth noting that, at the orbital distance of HD189733b, the SW is still in the process of accelerating, and an orbital speed of ≈150 km s⁻¹ provides a significant component to the total velocity of the planet relative to the background SW plasma. Figure 3 shows the global view of the escaping upper atmosphere of HD189733b interacting with the flow of SW plasma simulated in the model run N11 under the typical conditions specified above (see also Table 1).

Zoom In Zoom Out Reset image size

Table 1. The List of Simulation Scenarios with Modeling Parameters and Calculated Absorption Features

N	FXUV	He/H	Msw (×1010 g s−1)	ncut (cm−3)	Mpw (×1010 g s−1)	Abs. He i max,%	Half-width He i (km s−1)	Abs. Lyα, blue wing [−230, −140] km s−1	Abs. Lyα, red wing [60, 110] km s−1	Abs. Lyα, total	Others
1	15	0.006	10	3.3 × 105	17.3	0.95	22	5.8	7.6	6.3
2	15	0.006	10	3.3 × 106	17.1	1.1	21	5.9	8.0	6.8
3	15	0.006	10	3.3 × 107	9.8	1.2	16.5	6.3	8.5	7.1
4	15	0.006	10	3.3 × 108	5.8	1.2	14	6.4	7.8	6.9
5	15	0.006	10	3.3 × 1010	1.9	0.8	13.2	4.9	5.7	5.1
6	15	0.006	10	3.3 × 1016	1.6	0.73	12.8	4.7	5.4	4.9
7	15	0.006	10	3.3 × 1016	9.0	1.0	26.6	6.8	8.4	7.5	Vϕ = 19 km s−1
8	15	0.006	10	3.3 × 1016	8.0	1.25	19.5	3.9	4.6	7.7	Vϕ = 11 km s−1
9	15	0.006	10	3.3 × 1016	4.7	1.15	15.6	6.2	7.4	6.6	Vϕ = 5.5 km s−1
10	7.5	0.006	10	3.3 × 106	10.6	0.9	19.1	6.0	8.3	6.9
11	25	0.006	10	3.3 × 106	24.7	1.35	21.2	5.9	7.8	6.6
12	40	0.006	10	3.3 × 106	35.8	1.6	21.9	6.0	7.8	6.7
13	55	0.006	10	3.3 × 106	46.3	1.75	22.2	6.0	7.9	6.2	FXUVH/FXUV = 1/3
14	7.5	0.005	10	3.3 × 106	7.5	0.77	19.2	6.1	7.4	6.9
15	20	0.005	10	3.3 × 106	16.2	1.1	22.4	6.3	7.7	6.6
16	25	0.003	10	3.3 × 106	21.6	0.8	22	6.6	8.2	6.9
17	55	0.003	10	3.3 × 106	40.4	1.05	22.4	7.1	8.5	7.1
18	25	0.003	250	3.3 × 106	26.6	0.81	20.7	6.2	9.0	7.4
19	7.5	0.005	1000	3.3 × 106	7.5	0.55	16.4	6.1	7.4	6.5
20	20	0.005	2000	3.3 × 106	16.2	0.82	17.0	6.3	7.7	7.0
21	25	0.003	2000	3.3 × 106	21.6	0.69	16.9	6.6	8.2	7.5
22	55	0.003	4000	3.3 × 106	40.4	0.90	18.0	7.1	8.5	3.3
Observations						1.1 ± 0.1 a	22.5 a	14.4 ± 3.6 c	7.7 ± 2.7 c	5.0 ± 1.3 c
						0.7 ± 0.1 b	27.2 b			2.9 ± 1.4 d

Notes. Columns from left to right: number of the model run; assumed value of the integrated stellar XUV flux in the range of wavelengths 10 Å < λ < 910 Å in erg cm⁻² s⁻¹ at 1 au (the hard energy part λ < 510 Å of the flux is always taken to be F_XUVH = 2/3 × F_XUV except in the run N13); helium abundance in the base atmosphere; mass-loss rate of the star, which specifies the intensity of the SW; cutoff density value, which parameterizes the reduction of the Lyα cooling; calculated mass-loss rate of the planet; calculated absorption maximum and half-maximum width of the simulated He i(2³S) absorption profile at midtransit; simulated midtransit Lyα absorption averaged over the blue wing [−230, −140] km s⁻¹ and over the red wing [50, 230] km s⁻¹ of the line; simulated total midtransit Lyα absorption averaged over the whole line [−400, 400] km s⁻¹, with the exception of the [−50, 50] km s⁻¹ interval. The last column shows other modeling parameters varied in addition and described in the text. We chose specific absorption line averaging intervals to enable the comparison between simulations and observations. The bottom row indicates the measured absorption values provided in:

^a Salz et al. (2018). ^b Guilluy et al. (2020). ^c Lecavelier des Etangs et al. (2012; the measurements in 2011 September, supported also with 2007/2008 measurements). ^d Lecavelier des Etangs et al. (2010; the measurements in 2010 August).

Download table as:  ASCIITypeset image

The flow of atmospheric material occupies a significant area along and around the orbit, far ahead and well behind the planet. The realized, in this case, regime of interaction between the escaping atmospheric flow of HD189733b and the SW corresponds to the "captured by the star" regime, according to the classification introduced in Shaikhislamov et al. (2016). Because of the strong tidal and noninertial forces, which turn the planetary material stream clockwise and align it with the orbital plane, the flow is relatively thin across the planet−star direction. Due to stellar gravity, the stream is confined within the range ±8R_p across the orbital plane.

Figure 3 (right panel) demonstrates how the He i(2³S) absorption is distributed over the stellar disk as seen by the observer, i.e., the line-of-sight absorption. It takes place relatively close to the planet within a sphere of radius ∼3R_p.

The left panel in Figure 4 shows detailed information of the distributions of the densities of species, the velocity, and the temperature of planetary material along the planet−star line. Due to the high gravity, the outflow velocity is relatively low in comparison to that of the well-studied exoplanets HD209458b and GJ436b. The temperature reaches 15,000 K rather close to the planet and is higher than that obtained in the 1D simulations (for example, Salz et al. 2016), which is due to the reduced Lyα cooling adopted in our modeling. Half-ionization levels for H and He are both below 2R_p. The H⁺/H = 1 level appears at r = 1.6R_p, which is close to that predicted in aeronomy simulations in Menager et al. (2013), but is perceptibly farther than in 1D models of Guo (2011) and Salz et al. (2016). With the same value of XUV flux, we obtained for HD189733b the mass-loss rate of about three times that reported in Guo (2011) and Odert et al. (2020).

Zoom In Zoom Out Reset image size

From Figure 4 we can also surmise that PW is sufficiently collisional. For example, at a distance of 5R_p the helium atoms have the Knudsen number of about 0.1, calculated by collisions with protons with a cross section of about 4 × 10⁻¹⁶ cm².

The right panel shows the rates of chemical processes increasing and decreasing the metastable He i(2³S) population. One can see that the recombination of ionized helium He ii into the He i(2³S) state is counteracted by the autoionization collisions with H atoms at low altitudes r < 1.3R_p, and by the electron collisional depopulation above this height. Note that photoionization of the He i(2³S) state becomes important relatively far from the planet where the density of He i(2³S) atoms becomes too low to affect the absorption. Thus, the ionizing stellar flux in the NUV band is not crucial for the He i(2³S) absorption at 10830 Å for HD189733b. In the course of modeling, by varying the flux at around 10830 Å and comparing velocities, we found out that the radiation pressure causes no significant acceleration of He i(2³S) atoms for this exoplanet.

To find the model run that provides the best fit to observations of the He i(2³S) and Lyα spectral lines, we vary the value of the integral stellar XUV flux, which is subject to temporal fluctuations (Bourrier et al. 2020), and the helium abundance, which is a constant value to be constrained. Initially, we assume a very weak SW, which does not affect the expansion of the escaping planetary atmospheric material. Then, we compare the obtained results with the cases with a more intense SW characterized by a stellar mass-loss rate (M_sw).

A list of the model parameters and the corresponding absorption values calculated in different simulation runs is given in Table 1. We note that there are available for HD189733b two independent measurements of the He i(2³S) line that have significantly higher signal-to-noise ratios (S/N) and a greater similarity to each other than the Lyα observations have. Therefore, by comparison of the simulated He i(2³S) absorption profiles with the measurements, we constrain at first the helium abundance. The presence and amount of helium are important for the atmospheric escape of the relatively massive HD189733b, because they influence the hydrostatic scale height of the planetary atmosphere.

Figure 5 shows the excess absorption profiles at the position around 10830 Å caused by the He i(2³S) triplet, simulated for different values of the stellar XUV flux. The considered XUV flux varies from about an average solar value up to an order of magnitude higher. In the model run N13 we took the same features as in Lampón et al. (2021): XUV flux F_XUV = 55 erg cm⁻² s⁻¹ at 1 au, with the hard flux at λ < 510 Å three times less than the total (i.e., F_XUVH/F_XUV = 1/3). The helium abundance in the model runs presented in Figure 5 is fixed at He/H = 0.006. One can see that the simulated absorption depth and, to a lesser degree, the width of the absorption profile, increases with the increasing XUV flux. This is because of the more efficient ionization of helium and the faster expansion of the planetary atmosphere. Therefore, the factor of 1.5 difference between the measurements by Salz et al. (2018) and Guilluy et al. (2020) can be explained by an expected variability of the stellar XUV emission. In the right panel of Figure 5 we plot the results of two modeling scenarios aimed to fit these two different measurements. The first scenario assumes a very high XUV flux needed to reproduce the absorption of about 1% at the He i(2³S) line center, as measured in Salz et al. (2018), and a relatively low helium abundance He/H = 0.003. The second scenario takes a lower XUV flux, to fit absorption of about 0.7% at the line center, measured in Guilluy et al. (2020), and a slightly larger helium abundance of He/H = 0.005.

Zoom In Zoom Out Reset image size

Altogether, the model run N16 provides a relatively good correspondence with the measurements of Guilluy et al. (2020). Its reduced χ² statistics value, calculated in the interval [−40, 20] km s⁻¹ is 0.8, and it is the minimal value among all considered modeling runs. In calculating the χ² we used for the degrees of freedom the number of observational data points minus the number of modeling parameters (F_XUV and He/H).

Regarding the observations reported in Salz et al. (2018), the model does not reproduce the depth of the third line in the triplet, which appears to be larger than the stat-weight value of 1/8 relative to the sum of the first and second lines. Lampón et al. (2021) argued in that respect that the increasing atmosphere density (up to 10¹⁸ cm⁻³) at the level close to the optical radius of the planet results in a stronger absorption at the optically thick region and, consequently, the smaller ratio of the 10832/10833Å lines. However, from Figures 3 and 4 one can see that according to our modeling, in the dense region of the atmosphere of HD189733b at the altitudes r < 1.2R_p, the He i(2³S) atoms are much more efficiently depopulated rather than produced. It is worth noting that to include the optically thick region, we adopt in our model a sufficiently high density of 10¹⁷ cm⁻³ at the level of r = R_p.

Another feature that our model does not reproduce is a net blueshift of about 2.5 km s⁻¹ in the absorption data set in Salz et al. (2018), though such a shift is absent in the data from Guilluy et al. (2020). While the escaping planetary material flow has a complex spiraling structure and asymmetry, it does not result in any significant shift of the line core, because the absorption comes from the region where velocity of the flow is still rather small. With several dedicated model runs, we checked if there might be any reasonable conditions at which the radiation pressure acting on He i(2³S) atoms can accelerate them up to velocities necessary for the observed blueshift. However, the stellar flux at 10830 Å appears to be well constrained by the SED computed for the given stellar parameters, and it appears too low to produce such an acceleration during the average lifetime of a metastable atom.

An alternative way to explain the net blueshift in the He i(2³S) absorption profile is related to the zonal flows from the day side to the night side generated by the so-called hot spot in the atmosphere (Showman et al. 2015). Precise modeling of such flows requires a combination of our model of the escaping upper atmosphere of the planet with a GCM of its lower atmosphere. This challenging task is a subject of future work. At the moment, we can only try an empirical shift of the whole simulated He i(2³S) absorption profile (in the model runs N15 and N17), which results in the remarkably good correspondence with observations, characterized by a sufficiently low χ² = 0.01.

To demonstrate the He i(2³S) absorption in more detail, Figure 6 shows the distribution over the stellar disk, as seen by the observer, calculated for three Doppler velocity intervals: the blue wing [−20, −7] km s⁻¹, the red wing [7, 20] km s⁻¹, and the line core [−7, 7] km s⁻¹. One can see that the line core excess absorption is rather symmetric with a maximum value of about 15% outside the planet. Its appearance is restricted to the area within ∼3R_p. The absorption in the blue and red wings predominantly takes place behind and ahead the planet, respectively. This is related with the clockwise twist of the escaping stream of upper atmospheric material under the action of the Coriolis force. Therefore, the calculation of spectral absorption over the whole line, based on the predictions of simplified 1D models, leads to physically erroneous estimates and conclusions.

Zoom In Zoom Out Reset image size

Altogether, based on the performed modeling results presented above and the comparison with observations, we conclude that in the considered case of HD189733b, the stellar XUV flux is constrained to the range of 7.5–55 erg cm⁻² s⁻¹ at 1 au, whereas the helium abundance should be He/H =(3–5) × 10⁻³. We explain the difference between peak absorption values measured in two different observation campaigns by possible variations of the XUV flux, related to the activity of stars like HD189733.

Despite the lower atmosphere of HD189733b predominantly consisting of the molecular hydrogen, it is rapidly dissociated at heights exposed to the stellar vacuum ultraviolet (VUV) and XUV fluxes. In the presence of even a small subsolar amount of oxygen this dissociation proceeds much faster, facilitated by the molecular reactions with O and H₂O (García Muñoz 2007), which are not included in our model. This is why we start our simulations already with the atomic hydrogen atmosphere imposed at the inner boundary base. However, the trial run with an assumed molecular hydrogen base atmosphere resulted in an order of magnitude smaller He i(2³S) absorption and a significantly narrower line profile. This is a consequence of, in this case, the twice-reduced scale height and more compact upper atmosphere. Therefore, our modeling results and their comparison with the measurements of He i(2³S) absorption confirm the absence of molecular hydrogen in the upper atmosphere of HD189733b, which has to be dissociated due to the complex molecular chemistry.

3.3. Strong Stellar Wind and the Role of ENAs in Lyα Absorption

Now let us discuss the Lyα absorption. In this respect, our model produces results similar to those of Odert et al. (2020). The absorption is generated by natural line broadening at the level of 5%–7% relatively close to the planet (<3R_p) in the region populated with a sufficiently high amount of neutral atomic hydrogen. In the course of modeling we checked that the Lyα radiation pressure, acting on H atoms, does not accelerate them significantly up to the energies affecting the absorption profile. This is consistent with our previous studies dedicated to HD209458b (Shaikhislamov et al. 2020c; Khodachenko et al. 2017) and GJ436b (Khodachenko et al. 2019). Regarding HD189733b, an important question is whether the ENAs generated by charge exchange of planetary atoms with the SW protons can significantly add to the absorption at high blueshifted velocities of the Lyα line, similar to that seen in the simulations of HD209458b, as well as in the observations and modeling of GJ436b (Lavie et al. 2017; Khodachenko et al. 2019) and GJ3470b (Bourrier et al. 2018; Shaikhislamov et al. 2021).

Figure 7 shows several Lyα absorption profiles simulated for different SW densities and stellar XUV fluxes, taking the helium abundances that correspond to the best fits of the He i(2³S) absorption achieved in the model runs N14–N17. Under the conditions of a weak SW assumed in the model run N16, the simulated Lyα absorption profile appears to be very symmetric. However, the increasing intensity of the SW in the model runs N18–N22 results in a gradual increase of the absorption in the blue wing of the profile with respect to the red wing. This difference is related to the different character of interaction of the escaping planetary atmosphere and the SW and the corresponding structure of material flows and the distribution of the absorbing atomic hydrogen. In particular, already in the model run N18 with the parameters of a moderate SW similar to those in an average solar wind, a well-pronounced bowshock located upstream at ∼7R_p is formed, resulting in the production of ENAs and a significant (up to 1.5 times) increase of the absorption at the range of Doppler-shifted velocities around −75 km s⁻¹. However, in the range of higher velocities, [−230, −140] km s⁻¹, reported in Lecavelier des Etangs et al. (2012), the input of ENAs is still negligible. Under the conditions of more intense SWs with higher densities and velocities, which might take place, e.g., in CMEs, the bowshock region approaches extremely close to the planet, down to ∼3R_p. This results in a qualitative change of the Lyα absorption character. First, the absorption depth at the line center drops as much as several times, since the region of strong atmospheric absorption by planetary hydrogen atoms shrinks down to an area with a size of ∼3R_p around the planet, whereas under weak SW conditions the full size of this region was ∼5R_p. At the same time, significant (≈15%) Lyα absorption at the blue range of Doppler-shifted velocities around −100 km s⁻¹ appears. Only about half of this absorption is due to ENAs. The rest comes from the planetary hydrogen atoms, which are not involved in the charge exchange reaction with the SW protons, but are collisionally energized in the shock region where the density of SW protons is ∼10⁶ cm⁻³, and the planetary atoms are picked up by the SW flow. A similar effect was also found in Khodachenko et al. (2017) while studying with a 2D HD model the interaction of a hot Jupiter planetary wind with different stellar winds.

Zoom In Zoom Out Reset image size

Figure 8 shows the proton temperature and ENA density distributions around the planet simulated in the model run N21 under the conditions of a strong SW, which are similar to the case of CMEs considered in Odert et al. (2020). These plots clearly show the bowshock formed around the planet and details of the planetary and stellar material motion. According to Figure 7, the effect of ENAs is well pronounced in the simulated Lyα absorption, which was not the case in Odert et al. (2020). One of the reasons for this difference is that we model the interaction of the planetary and stellar material flows in full 3D, and most ENAs are produced in the tail region, which is absent in the simulations of Odert et al. (2020).

Zoom In Zoom Out Reset image size

To shed more light on the origin of the simulated and measured Lyα absorption in HD189733b, we present in Figure 9 the distribution of this absorption across the stellar disk in different Doppler velocity intervals and produced by different populations of particles. The absorption over the whole line in the range [−150, 150] km s⁻¹, except of the geo-contaminated interval (±50 km s⁻¹), is more or less symmetric. It originates from a relatively close projected region around the planet within ≈3R_p. However, it is worth noting that this projected region contains also an additional input from the hydrogen atoms in the tail. The input of tail particles becomes more evident at the blue wing of the Lyα line, especially in the part of absorption contributed by the ENAs generated by the charge exchange reaction. These tail particles are clearly seen in Figure 8 (right panel).

Zoom In Zoom Out Reset image size

Finally, in Figure 10 we present how the strong SW might affect the He i(2³S) line absorption. One can see that while the moderate SW does not influence the absorption profile, the strong SW significantly reduces the absorption in the red wing, as well as the overall width and depth of the profile. This result leads us to the conclusion that at the time of the measurements of the He i(2³S) absorption discussed in this paper, the SW could not be strong.

Zoom In Zoom Out Reset image size

Despite a number of FUV spectral observations performed with the Hubble Space Telescope (HST), the existing stellar activity and variability prevent an accurate constraining of the parameters of the upper atmosphere of the hot Jupiter HD189733b. While the measurements of the 5% depth in the whole Lyα line are interpreted within the frame of 1D HD models as the natural line broadening absorption by the planetary hydrogen in the expanding atmosphere confined in the region of a few planetary radii, the measurement of the 14% depth at the high velocity >140 km s⁻¹ blue wing of the Lyα line obtained in a single observational session remains enigmatic for hydrodynamic models. Bourrier & Lecavelier des Etangs (2013) report that the Monte Carlo model is able to fit this particular observation with the radiation pressure acceleration of atomic hydrogen and charge exchange with SW protons. It is obtained at planetary mass-loss values of 4 × 10⁸–4 × 10⁹ g s⁻¹, which are significantly lower than those predicted by gas dynamic models at similar XUV fluxes of 10–30 erg cm⁻² s⁻¹ at 1 au. At the same time, according to Bourrier & Lecavelier des Etangs (2013), the stellar wind at the planet orbit should have a velocity of 200 km s⁻¹ and a rather low temperature of 3 × 10⁴ K, while the density can vary in a wide range around 10⁵ cm⁻³. However, the 3D HD modeling in Odert et al. (2020), as well as that presented in this paper, do not confirm the generation of ENAs with velocities in excess of 140 km s⁻¹ by either of these mechanisms. At the same time, we found a significant increase of the Lyα absorption in the range [−50, −150] km s⁻¹, with an average absorption value of up to 16% under the conditions of a strong SW, characterized by a velocity of 240 km s⁻¹ and a density of >10⁵ cm⁻³, while the temperature, as dictated by our empirical heating model, should be around 10⁶ K. Note that the absorption of 13.2% at velocities [−117, −156] km s⁻¹ measured by HST at the 2013 November 3 visit (Bourrier et al. 2020) agrees much better with our simulation. One of the physical features that can significantly influence the interaction of planetary and stellar particles and the interpretation of observations is the planetary and stellar magnetic fields (Pillitteri et al. 2015; Khodachenko et al. 2021b). Investigation of this possibility remains a subject of future studies.

The measurements of absorption at the position of the metastable He i(2³S) triplet line, performed with ground-based telescopes, offered a new possibility to verify the existing models and physical scenarios and to constrain the parameters of the HD189733b system. With the assumed cooling of the planetary upper atmosphere by excitation of the hydrogen atoms, the simulated width of the He i(2³S) absorption appeared significantly narrower than that in the observations. This is because the upper atmosphere of HD189733b under such an assumed cooling is in a radiative balance rather than in an outflow regime. Salz et al. (2016) reported the same conclusion and approximately similar parameters of the planetary wind. Within the frame of our model we could reproduce the width of the measured He i(2³S) absorption profile by reducing the overall effect of Lyα cooling in the energy balance by about two times. This increases by about 2.5 times the part of energy going to the acceleration of the escaping atmosphere flow and increases the mass-loss rate by an order of magnitude. Note that, at the same time, the absorption in the Lyα line is only slightly affected by the cooling.

With the reduced Lyα cooling, our simulations indicate that the helium abundance could not be larger than He/H = 0.01, otherwise the synthetic absorption profiles become too narrow to match the observations. The peak He i(2³S) absorption is matched in the range of stellar XUV fluxes (7.5−55) erg cm⁻² s⁻¹ at 1 au, whereas the helium abundance is constrained to be within the range of (3−5) × 10⁻³. Our results are in better agreement with the reconstructed SED of Bourrier et al. (2020) at a lower He/H ratio (run N16). We explain the factor of 1.5 difference between two available measurements of the He i(2³S) absorption by a 2−2.5 times natural fluctuation of the stellar XUV flux. Another constraint that our modeling enables us to place, based on the comparison of the measured and simulated He i(2³S) absorption profiles, is that the SW could not be very strong during the observations.

Remarkably, the best-fitting result regarding the helium abundance, (5–8) × 10⁻³ in Lampón et al. (2021), obtained at the mass-loss rate of 1.1 × 10¹¹ g s⁻¹, is quite close to the predictions of our global 3D HD modeling. However, the approaches to derive these values are principally different. In fact, the mass-loss rate and temperature of the escaping planetary flow are not free parameters of our model. They are calculated in the simulations in dependence on more fundamental and independent physical parameters of the system, namely, the integral stellar XUV flux, which is subject to temporal variability, and the helium abundance, which is not known and has to be constrained. It turned out that the Lyα absorption, which is to a significant extent produced by the planetary hydrogen due to the natural line broadening, is quite insensitive to the XUV flux and the related atmospheric mass loss. Altogether, our simulations show that the Lyα absorption cannot be used to derive the planetary mass loss and the helium abundance, as was done in Lampón et al. (2021).

In general, our simulations show that global 3D fully self-consistent multicomponent modeling produces results consistent with observations at transits of HD189733b in the Lyα and He i(2³S) lines at physical parameters that are well within the expected ranges. The synthetic profiles simulated with our model for the He i(2³S) absorption do not fit the observation data points as close as they do in Lampón et al. (2021). This is because our self-consistent modeling approach has a limited number of free parameters, specifically the XUV radiation flux with its energy partition F_XUXH/F_XUV and element abundances, as well as SW parameters and the efficiency of Lyα cooling. Therefore, the insignificant difference between the simulations and measurements in our case is more likely because of additional physical effects and processes not covered by the model, which nevertheless well reproduces the major observational phenomena. Also, we do not reproduce the net blueshift of the He i(2³S) absorption profile by a few km s⁻¹, though it greatly improves the fit between the modeling and observations.

At the same time, we can add in our model one more free parameter, assuming that the planet is not completely tidally locked, but rotates, or that there exists in the atmosphere a strong eastward zonal jet (Showman et al. 2015). We checked this possibility versus the observations and determined that the rotational azimuthal velocities at the inner boundary of the simulation domain, i.e., at the base of the modeled upper atmosphere, might influence the width of the absorption profile, given they are of sufficiently large values of more than 10 km s⁻¹. However, a detailed conclusion regarding the role of the zonal winds is a subject for a special study, because the parameters of an atmospheric wind cannot be formally inserted in a self-consistent model to produce the desired effect. All of the assumptions and prescribed features should have firm physical grounds, or be provided as an outcome of other dedicated and integrated models.

This work was performed under the grant No. 18-12-00080 of the Russian Science Foundation. Analysis of the Lyα cooling effect was accomplished in the frame of the of grant No. 21-72-00129 of the Russian Science Foundation. M.S.R. and A.G.B. are thankful to the project "Study of stars with exoplanets" under a grant from the Government of the Russian Federation for scientific research conducted under the guidance of leading scientists (agreement No. 075-15-2019-1875) for supporting the investigation of Lyα absorption properties and additional rotation affects on He i(2³S) absorption properties. M.L.K. acknowledges the projects I2939-N27 and S11606-N16 of the Austrian Science Fund (FWF). H.L. acknowledges the projects P25256-N27 and S11607-N16 of the Austrian Science Foundation (FWF). I.F.S. is thankful to the program "Astrophysical Origins: Pathways from Star Formation to Habitable Planets" (July 2019, ESI, Wien, Austria) for productive discussions. I.B.M. and A.G.B. acknowledge the support of Russian Fund of Basic Research (RFBR) grant No. 20-02-00520 for calculating the stellar wind parameters. Parallel computing simulations, key for this study, have been performed at the Computation Center of Novosibirsk State University, the SB RAS Siberian Supercomputer Center, and the Joint Supercomputer Center of RAS.