Global 3D simulation of the upper atmosphere of HD189733b and absorption in metastable HeI and Ly{\alpha} lines

A 3D fully self-consistent multi-fluid hydrodynamic aeronomy model is applied to simulate the hydrogen-helium expanding upper atmosphere of the hot Jupiter HD189733b, and related absorption in the Lya line and the 10830 A line of metastable helium. We studied the influence of a high-energy stellar flux, stellar wind, and Lya cooling to reproduce the available observations. We found that to fit the width of the absorption profile in 10830 A line the escaping upper atmosphere of planet should be close to the energy limited escape achieved with a significantly reduced Lya cooling at the altitudes with HI density higher than 3*10^6 cm^-3. Based on the preformed simulations, we constrain the helium abundance in the upper atmosphere of HD189733b by a rather low value of He/H~0.005. We show that under conditions of a moderate stellar wind similar to that of the Sun the absorption of Lya line takes place mostly within the Roche lobe due to thermal broadening at a level of about 7%. At an order of magnitude stronger wind, a significant absorption of about 15% at high blue shifted velocities of up to 100 km/s is generated in the bowshock region, due to Doppler broadening. These blue shifted velocities are still lower than those (~200 km/s) detected in one of the observations. We explain the differences between performed observations, though not in all the details, by the stellar activity and the related fluctuations of the ionizing radiation (in case of 10830 A line), and stellar wind (in case of Lya line).

There are only a few works dedicated to the interpretation of the Lyα observations of the HD189733b. Bourrier et al. (2013) used a 3D Monte Carlo code, tracing the hydrogen atoms launched from an inner spherical boundary, surrounding the planet at 3Rp. They found that acceleration of the atoms by radiation pressure, alone is not sufficient for the generation of energetic neutral atoms (ENA) population with the velocities of up to 200 km/s, and 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 planetary wind was combined with a 3D MHD model of the magnetized stellar wind (SW). For the XUV flux of 1.8·10 4 erg cm -2 /s, they obtained the upper atmosphere maximum temperature of ≈1.1·10 4 K, the mass loss rate of 5·10 10 g/s, and the location of H/H + =1 level at ≈1.7Rp. Odert et al. (2020) also found that the absorption in Lyα line is mostly due to the natural line broadening mechanism, and is of ~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 blue shifted absorption in excess of 140 km/s reported in Lecavelier des Etangs et al. (2010Etangs et al. ( , 2012. The parameters of the considered strongest SW corresponded to the conditions in CMEs with the velocity of 1000 km/s and the integral stellar mass loss 10 3 times higher than the average solar value (≈2.5·10 12 g/s). Under these conditions the magnetopause was pushed by the SW ram pressure towards the planet as close as 2.3Rp.
Since 2018 a new opportunity for probing of 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 metastable helium 2 3 S triplet line at 10830 Å (hereafter HeI(2 3 S)). The metastable 2 3 S level of helium atom has a high excitation energy of 19.8 eV and a radiative life-time of about 2.2 hours. It is efficiently pumped by the recombination of ionized helium and results in the absorption in infrared band due to the triplet transition 2 3 S-2 3 P at 10830.34 Å, 10830.25 Å, and 10829.09 Å. The HeI(2 3 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 HII regions (Aver et al. 2015), as well as to probe the material flows in active galactic nuclei (Leighly et al. 2011). Seager & Susselov (2000) proposed using this line to observe the transits of HD209458b. Determination of the population of HeI(2 3 S) state in the upper atmospheres of hot exoplanets is a rather complex problem tackled in Oklopčić & Hirata (2018) by means of 1D HD model. They predicted the absorption at the line center at a level of 8% and 2% for GJ436b and HD209458b, respectively.
Recent IR observations of exoplanetary transits with ground telescopes, mostly with CARMENES instrument on the Calar Alto and NIRSPEC on the Keck II, provided a series of valuable data for the analysis and interpretation. The HeI(2 3 S) absorption depths of ~8% for the Wasp107b (Allart et al. 2019;Kirk et al. 2020), ~1% for the HATP11b , 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 the Wasp69b ), 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 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 half-width of the HeI(2 3 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 blue-shift of the whole absorption profile by ≈2.5 km/s reported in Salz et al. (2018), but absent in Guilluy et al. (2020) data. We note that signatures of a small blue-shift is also reported for other planets. However, Salz et al. (2018) remarks that these blue-shift velocities could be potentially caused by the stellar pseudo-signals, i.e. due to the stellar surface regions with a stronger HeI 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, electron, hydrogen and helium density distribution profiles in the upper atmosphere and the kinetics of HeI(2 3 S) state population. Ninan et al. (2019) and Palle et al. (2020) found for the GJ3470b that the area of HeI(2 3 S) presence extends up to at least 10Rp, and 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 HeI(2 3 S) absorption, based on 1D HD modeling is described by Lampón et al. (2021). Because we explored the GJ3470b in our previous work (Shaikhislamov et al. 2020), 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. Besides of temperature, other fitting parameters are the mass loss rate and the helium abundance. The fixed XUV flux of 5.7·10 4 erg cm -2 /s is used in the model to calculate the ionization balance and the kinetics of HeI(2 3 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 the HD189733b (altogether, six empirical fitting parameters), whereas for the GJ3470b only the net blue shift appeared to be sufficient. However, the fitting of the HeI(2 3 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 temperature and the mass loss derived from the best fit of the Parker-like solution profiles to the profiles from Salz et al. (2016 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 blue-shifted (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, like 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) perform one of the fits of the He(2 3 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, artificially adding of the non-spherically symmetric effects, which only a global 3D multi-component model can simulate 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 in lines of different species have been made, like e.g., in HI (Lyα), OI, CII, HeI(2 3 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 on this way and perform a multi-fluid 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.
Besides of the simulation and interpretation of observations of HD189733b, we compare our results with other modeling approaches.
The model we use in our study here has been already applied for the simulation and interpretation of observations of GJ436b , π Men c (Shaikhislamov et al. 2020a), HD209458b , 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 Lyα and HeI(2 3 S) lines revealed that bowshock layer formed by SW around the planet plays significant role. They found the best fit parameters for both lines as FXUV=8 erg cm -2 /s, He/H=0.013, and the SW with density and velocity similar to the fast solar wind, strong enough to form a bowshock at about 20Rp (Shaikhislamov et al. 2021).
The simulation of absorption at 10830 Å during the transits of Wasp107b revealed another feature which is beyond the capacity of 1D models; it concerns the radiation pressure, acting on metastable HeI(2 3 S) atoms. As pointed out in Allart et al. (2019), the acceleration of helium atoms by IR (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. (2021) the observed HeI(2 3 S) absorption at Wasp107b can be well reproduced with the typical expected parameters of the system and appropriate account of the electron and atom impact processes, without invoking of any specific additional assumptions on 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 HeI(2 3 S) measurements at transits of HD189733b within a single model. In this respect it is worth to mention 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 well, like in Odert et al. (2020), weather a strong SW can generate sufficient amount of ENAs to affect the Lyα blue wing absorption and to explain the significant difference of the measured Lyα absorption profiles in different observations. Based on the performed simulations and fitting to the measured HeI(2 3 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 short description of the model and provides information about initial parameters; Section 3 presents the obtained results, which are discussed in Section 4.

Modeling approach
The 3D multi-fluid global model used in this paper is being developed since Shaikhislamov et al. (2014). Its most important parts are described in the following papers: the details of gasdynamic equations and terms related to generation of planetary wind, collisional interaction of species and applicability of hydrodynamic approach are given in the Shaikhislamov et al. (2016) and in the appendix of Dwivedi et al. (2019); the details of 3D code and of stellar plasma wind simulationin Khodachenko et al. (2019); procedures of calculating transit absorption in linesin Shaikhislamov et al. (2018b); simulation of metastable helium population and absorptionin Shaikhislamov et al. (2021a). Here we only briefly outline the model. It solves numerically 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 HI and HeI(2 3 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+, He2+. An additional argument in support of the atomic upper atmosphere assumption comes from the modeling of HeI(2 3 S) absorption discussed later in the paper. However, in several dedicated model runs we studied the effects of molecules and included H2, H2 + , and H3 + into consideration.
The metastable HeI(2 3 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 they originate, i.e., He + in case of recombination, or HeI(1 1 S), in case of excitation from the ground state. The reactions involved in the processes of population and depopulation of the HeI(2 3 S) triplet level are similar to those considered in Oklopčić & Hirata (2018) and described in Shaikhislamov et al. (2021a).
The code uses an explicit numerical scheme of the second order spatial and temporal accuracy, realized on a spherical grid in the planet-centered reference frame with polar axis Z, perpendicular to the orbital plane. More details about the model implementation can be found in Shaikhislamov et al. (2018aShaikhislamov et al. ( , 2020aShaikhislamov et al. ( , 2021 and Khodachenko et al. (2019Khodachenko et al. ( , 2021. At the inner boundary of the simulation domain, i.e., the planet surface at r=Rp, we fix zero radial velocity of the atmospheric material, the base atmosphere temperature of 1000 K and pressure of 0.05 bar. The photoionization and heating of gas by photo-electrons, as well as the radiation attenuation in the planetary atmosphere, is calculated over the XUV and NUV spectra binned by 1 Å according to wavelength-dependent cross-sections of various elements involved in 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 of this spectra is FXUV ≈25 erg cm -2 /s at the reference distance of 1 a.u. The latest reanalysis of available X-ray and FUV observations and modeling of synthetic EUV spectra of HD189733 (Bourrier et al. 2020) produced stellar XUV fluxes varying from visit to visit from 23 to 31 erg cm -2 /s at 1 a.u. At the same time, we distinguish in our model a hard energy part of the flux (XUVH) at λ<510 Å with FXUVH =16 erg cm -2 /s at 1 a.u., which influences the metastable helium population. The ratio of FXUVH/FXUV =2/3 was kept unchanged in course of the modeling (except of specially addressed cases), whereas the value of total integral XUV flux FXUV was varied for different model runs. The employed SED is characterized in the near-ultraviolet part (NUV, 912<λ<3000 Å, effects HeI(2 3 S) population by its photo-ionization) by FNUV=1830 erg cm -2 /s at 1 a.u., while in the and near-infrared part (NIR, around 10830 Å, is used to calculate the radiation pressure force acting on HeI(2 3 S)) by F10830 =20 erg cm -2 /s Å -1 at 1 a.u.
The absorption is calculated post-simulation by a special procedure described in Khodachenko et al. (2017) and Shaikhislamov et al. (2018bShaikhislamov et al. ( , 2020a for Lyα line and in Shaikhislamov et al. (2021a) for HeI 2 3 P-2 3 S transitions. For the latter, the synthetic spectral profiles are obtained, like those in 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.

The role of the Lyα cooling of upper atmosphere
The cooling of expanding upper atmospheric material due to excitation of hydrogen atoms and emission of the Lyα photons is an important effect in partially ionized upper atmospheres of hot Jupiters. However, its quantitative description is non-trivial because of the trapping and diffusion of the Lyα photons in dense regions of the planetary atmosphere. The first aeronomy codes (Yelle et al. 2004, García Muñoz et al. 2007, Koskinen et al. 2007) ignored this Lyα cooling effect, assuming eventual loss of Lyα photons to photo-ionization of the excited states of trace elements.
However, Murray-Clay et al. (2009) simulated the escaping planetary wind with the account of Lyα cooling and have shown that it changes the energy budget of a typical hot Jupiter by an order of 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 the 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 the HD189733b exactly this regime of a quasi-stationary upper atmosphere with low mass-loss rate of (1-2)·10 10 g/s. At the same time, there exist much more sophisticated simulations of the stationary upper atmosphere of HD189733b based on the direct Monte-Carlo modeling (Christie et al. 2013 andHuang et al. 2017). This approach takes into account a number of processes, which terminate the Lyα scattering cycle. Among them are photoionization from the n=2 state, collisional deexcitation, two photon decay, and photoionization of metals, e.g., Na and Mg. Altogether, such Monte-Carlo simulation 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, which 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, which will be specified later, we consider in this work the effectiveness of the Lyα cooling depending on an empirical parameter defined via so-called density cut-off in the following way. We reduce the Lyα cooling rate by a scaling factor exp(-σHnH/α), where nH is atomic hydrogen density, σ≈6·10 -14 cm 2 is cross-section of Lyα photon absorption at a typical temperature of 10 4 K, and H is 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 ncut=α/σH corresponds to the density layer below which the 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 becomes close to the radiative hydrostatic equilibrium. In this regime, the total rates of Lyα-emitted and XUV-absorbed energy are related as WLyα/WXUV=0.76 with WXUV≈4·10 24 erg/s. The corresponding mass loss rate is ≈1.5·10 10 g/s, maximum temperature ≈10 4 K, and velocity at 2Rp about 1.5 km/s. All these values are rather close to those obtained in Salz et al. (2016), despite the differences between the 1D and 3D models. However, without account of the Lyα cooling the mass loss rate becomes an order of magnitude higher, 2.2·10 11 g/s, with the temperature maximum 1.5·10 4 K, and velocity at 2 Rp 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 velocity of material escape are higher, mostly because of the higher temperature.
The most important point here is, that the regimes with and without Lyα cooling are very much different with respect to the HeI(2 3 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 HeI(2 3 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 Lyα cooling (including fulland no-cooling cases), obtained with different values of the defined above cut-off density, which controls the effective trapping of Lyα photons. It shows that to make the shape of the synthetic HeI(2 3 S) absorption profiles closer to the measurements, we should assume significant reduction of the Lyα cooling, with the cut-off density at least ncut=3·10 6 cm -3 , which corresponds to the trapping parameter σHnH >10 2 . In this case, the energy rates are related as WLyα/WXUV =0.37. According to our approach, the case of practically no Lyα cooling corresponds to ncut=3·10 5 cm -3 , whereas full Lyα cooling, to ncut=10 11 cm -3 . It should be noted that our model does not show the increased absorption at the weakest component of triplet between [-50; -20] km/s, as well as the small net blue shift observed by Salz et al. (2018) (though 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 broadening of the HeI(2 3 S) absorption profile, namely the zonal wind, whose speed for HD189733b might reach up to 10 km/s at altitudes with pressure 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φ =Rp·Ω·sinθ at the inner boundary of the simulation domain at r=Rp. Here Ω is the cyclic frequency of the prescribed rotation. Right panel in Figure 2 shows the simulated absorption profiles obtained in the case of full Lyα cooling (ncut =3·10 16 cm -3 ) 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} hours. In particular, it can be seen that the close to observations width of the absorption profile is achieved at 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 cut-off density ncut =3·10 6 cm -3 . This value allows obtaining the most plausible widths of spectral profiles, which fit the observations at HeI(2 3 S) line.  Table 1 given below), the corresponding cut-off densities and widths of the absorption profile.
Right panel: HeI(2 3 S) triplet absorption profiles simulated in the case of full Lyα cooling for different planet rotation velocities (provided in legend) emulating the effect of global circulation and corresponding to the diurnal periods of 8, 12 and 24 hours.

Constraining XUV flux and helium abundance
At first, we consider the case of a weak SW with the following parameters at the planetary orbit: Vsw =210 km/s, Tsw =0.7 MK, nsw =10 3 cm -3 . It is worth to note that at the orbital distance of HD189733b, the SW is still in the process of acceleration, and orbital speed of ≈150 km/s provides a significant component to the total velocity of the planet relative the background SW plasma. Figure 3 shows the global view of the escaping upper atmosphere of the HD189733b interacting with the flow of SW plasma simulated in the model run N11 under the typical conditions specified above (see also Table 1).
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 non-inertial 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 the stellar gravity, the stream is confined within the range ±8Rp across the orbital plane.   Table 1), which corresponds to the case of weak SW, no zonal wind, and reduced Lyα cooling. The planet is located at the center of coordinates and rotates anti-clockwise around the star marked as black circle at X=55. Proton fluid streamlines originated from star are shown with black lines, whereas gray lines indicate the streamlines of planetary material flow.
Right panel: The distribution of the line-of-sight absorption at the position of HeI(2 3 S) 10830 Å line integrated over ±15 km/s, as seen by a remote observer at the mid-transit. The stellar limb is the picture circumference; the planet is marked by black circle.
Left panel in Figure 4 shows detailed information on the distribution of densities of species, velocity and temperature of planetary material along the planet-star line. Due to the high gravity, the outflow velocity is relatively low in comparison to other well studied exoplanets HD209458b and GJ436b. The temperature reaches 15000 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 2Rp. The H + /H=1 level appears at r=1.6Rp, which is close to that predicted in aeronomy simulations in Menager et al. (2013), but perceptibly farther than in 1D models of Guo (2011) and Salz et al. (2016). With the same value of XUV flux, we obtained for the HD189733b the mass loss rate of about three times of that reported in Guo (2011) and Odert et al. (2020).
From Figure 4 we can also surmise that PW is sufficiently collisional. For example, at a distance of 5Rp the helium atoms have the Knudsen number of about 0.1, calculated by collisions with protons with a cross section of about 4·10 -16 cm 2 . Right panel shows the rates of chemical processes increasing and decreasing the metastable HeI(2 3 S) population. One can see that the recombination of ionized helium HeII into HeI(2 3 S) state is counteracted by the auto-ionization collisions with H atoms at low altitudes r<1.3Rp, and by the electron collisional depopulation above this height. Note that photo-ionization of HeI(2 3 S) state becomes important relatively far from the planet where the density of HeI(2 3 S) atoms becomes too low to affect the absorption. Thus, the ionizing stellar flux in the NUV band is not crucial for the HeI(2 3 S) absorption at 10830 Å for the HD189733b. In the course of modeling, by varying flux at around 10830 Å and comparing velocities, we found out that the radiation pressure causes no significant acceleration of HeI(2 3 S) atoms for this exoplanet. To find the model run which provides a best fit of observations in HeI(2 3 S) and Lyα spectral lines, we vary the value of integral stellar XUV flux, which is a 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 doesn't affect the expansion of the escaping planetary atmospheric material. Then we will compare the obtained results with the cases of more intense SW characterized by 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 the available for HD189733b two independent measurements at HeI(2 3 S) line have significantly higher S/N ratio and similarity between each other, than the Lyα observations. Therefore, by comparison of the simulated HeI(2 3 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. 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 -2 s -1 at 1 a.u. (the hard energy part λ<510 Å of the flux is always taken as FXUVH =2/3·FXUV except the run N13); helium abundance in the base atmosphere; mass loss rate of the star which specifies the intensity of SW; cut-off density value which parametrizes the reduction of the Lyα cooling; calculated mass loss rate of the planet; calculated absorption maximum and half-maximum width of the simulated HeI ( Figure 5 shows the excess absorption profiles at the position around 10830 Å caused by the HeI(2 3 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 than that. In the model run N13 we took the same features as in Lampón et al. (2021): XUV flux of FXUV =55 erg cm -2 /s at 1 a.u., with the hard flux at λ<510 Å three times less than the total (i.e., FXUVH/FXUV =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, in lesser degree, the width of the absorption profile, increase with the increasing XUV flux. This is because of more efficient ionization of helium and 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 HeI(2 3 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 bit higher helium abundance of He/H=0.005.
Altogether, the model run N16 provides relatively good correspondence with the measurements of Guilluy et al. (2020). Its reduced χ 2 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 χ 2 we used for the degrees of freedom the number of observational data points minus the number of modeling parameters (FXUV 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 18 cm -3 ) 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.2Rp, the HeI(2 3 S) atoms are much more efficiently depopulated rather than produced. It is worth to note that to include the optically thick region, we adopt in our model sufficiently high density of 10 17 cm -3 at the level of r=Rp. Another feature, which our model does not reproduce, is a net blue shift of about 2.5 km/s in the absorption data set in Salz et al. (2018), though such 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 HeI(2 3 S) atoms can accelerate them up to velocities necessary for the observed blue shift. 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 acceleration during the average life-time of metastable atom.
An alternative way to explain the net blue shift in the HeI(2 3 S) absorption profile is related to the zonal flows from the dayside 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 HeI(2 3 S) absorption profile (in the model runs N15, N17), which results in the remarkably good correspondence with observations, characterized by sufficiently low χ 2 =0.01.
To demonstrate the HeI(2 3 S) absorption in more details, 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; 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 ~3Rp. The absorption in 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 Coriolis force. Therefore, the calculation of spectral absorption over the whole line, based on the predictions of simplified 1D models leads to the physically erroneous estimates and conclusions. Altogether, based on the performed modeling results, presented above, and 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 -2 /s at 1 a.u., whereas helium abundance should be He/H=(3-5)·10 -3 . We explain the difference between peak absorption values measured in two different observation campaigns by possible variations of XUV flux, related with the activity of stars like HD189733.
Despite of the lower atmosphere of HD189733b dominantly consisting of the molecular hydrogen, it is rapidly dissociated at heights exposed to the stellar VUV and XUV fluxes. In the presence of even small sub-solar amounts of oxygen this dissociation proceeds much faster, facilitated by the molecular reactions with O and H2O (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 assumed molecular hydrogen base atmosphere, resulted in an order of magnitude smaller HeI(2 3 S) absorption and its significantly narrower profile. This is a consequence of the twice reduced in this case scale-height and more compact upper atmosphere. Therefore, our modeling results and their comparison with the measurements of HeI(2 3 S) absorption confirm the absence of the molecular hydrogen in the upper atmosphere of HD189733b, which has to be dissociated due to the complex molecular chemistry.

Strong stellar wind and the role of ENAs in Lyα absorptions
Now let us discuss the Lyα absorption. In that respect, our model produces the 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 (<3Rp) in the region populated with the sufficiently high amount of the 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. 2020, Khodachenko et al. 2017) and GJ436b .
Regarding the 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 blue-shifted 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 and GJ3470b , Shaikhislamov et al. 2021). Figure 7 shows several Lyα absorption profiles simulated for different SW densities and stellar XUV fluxes, taking the helium abundances, which correspond to best fits of the HeI(2 3 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 much symmetric. However, 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 with the different character of interaction of the escaping planetary atmosphere and 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 ~7Rp is formed, resulting in the production of ENAs and 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 the CMEs, the bowshock region approaches extremely close to the planet, down to ~3Rp. This results in a qualitative change of the Lyα absorption character. First, the absorption depth at the line center drops up to several times, since the region of strong atmospheric absorption by the planetary hydrogen atoms shrinks down to an area with size of ~3Rp around the planet, whereas under a weak SW conditions the full size of this region was ~5Rp. 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 part 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 6 cm -3 , and the planetary atoms are picked-up by the SW flow. Similar effect was also found in Khodachenko et al. (2017) while studying with a 2D HD model the interaction of hot Jupiter planetary wind with different stellar winds. Figure 8 shows the proton temperature and ENA density distributions around the planet simulated in the model run N21 under the conditions of 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 the most ENAs are produced in the tail region, which is absent in the simulations of Odert et al. (2020).
To shed more light on the origin of the simulated and measured Lyα absorption in HD189733b, we present in Figure 9 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 of [-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 ≈3Rp. However, it is worth to note 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). Finally, in Figure 10 we present how the strong SW might affect the HeI(2 3 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 discussed in this paper measurements of the HeI(2 3 S) absorption, the SW could not be strong.

Conclusions
Despite a number of FUV spectral observations performed with 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 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 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 et al. (2013) reports that the Monte Carlo model is able to fit this particular observation by the radiation pressure acceleration of the atomic hydrogen and charge exchange with SW protons. It is obtained at planetary mass losses 4·10 8 -4·10 9 g/s, which are significantly lower than those predicted by gasdynamic models at similar XUV fluxes of 10-30 erg cm -2 /s at 1 a.u. At the same time, according to Bourrier et al. (2013), the stellar wind at the planet orbit should have velocity of 200 km/s and rather low temperature of 3·10 4 K, while density can vary in a wide range around 10 5 cm -3 . 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 the velocities in excess of 140 km/s by either of the 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 up to 16% under the conditions of a strong SW, characterized by velocity of 240 km/s and density of >10 5 cm -3 , while the temperature, as dictated by our empirical heating model, should be around 10 6 K. Note that the measured by HST the absorption of 13.2% at velocities of [-117; -156] km/s at the 3.11.2013 visit (Bourrier et al. 2020), agrees much better with our simulation. One of the physical features, which can significantly influence the interaction of planetary and stellar particles and 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 HeI(2 3 S) triplet line, performed with the 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 HeI(2 3 S) absorption appeared significantly narrower than that in the observations. This is because the upper atmosphere of the 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 planetary wind. Within the frame of our model we could reproduce the width of the measured HeI(2 3 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 Lyα line is only slightly affected by the cooling.
With the reduced Lyα cooling, our simulations indicates 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 HeI(2 3 S) absorption is matched in the range of stellar XUV fluxes (7.5-55) erg cm -2 /s at 1 a.u., whereas the helium abundance is constrained within the range of (3-5)·10 -3 . Our results are in better agreement with the reconstructed SED of Bourrier et al. 2020 at lower He/H ratio (run N16). We explain the factor of 1.5 difference between two available measurements of the HeI(2 3 S) absorption by 2-2.5 times natural fluctuations of the stellar XUV flux. Another constrain that our modeling enables to put, based on the comparison of the measured and simulated HeI(2 3 S) absorption profiles, is that the SW could not be very strong during the observations.
Remarkably, the best-fit result regarding the helium abundance, (5-8)·10 -3 in Lampón et al. (2021), obtained at the mass loss rate of 1.1·10 11 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 the 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 a subject to temporal variability, and helium abundance, which is not known and has to be constrained. It turned out that the Lyα absorption, which is to a significant extent is produced by the planetary hydrogen due to the natural line broadening, is quite insensitive to XUV flux and the related with it 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 it was done in Lampón et al. (2021).
In general, our simulations show that the global 3D fully self-consistent multi-component modeling produces the results consistent with observations at transits of HD189733b in Lyα and HeI(2 3 S) lines at physical parameters which are well within the expected ranges. The synthetic profiles simulated with our model for the HeI(2 3 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 FXUXH/FXUV, elements abundance, as well as SW parameters and 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 blue shift of the HeI(2 3 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 there exists in the atmosphere a strong eastward zonal jet (Showman et al. 2015). We checked this possibility versus the observations and obtained 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, detailed conclusion on the role of the zonal winds is a subject for a special study, because the parameters of an atmospheric winds cannot be formally inserted in a self-consistent model to produce the desired effect. All the assumptions and prescribed features should have firm physical grounds, or be provided as an outcome of other dedicated and integrated models.