Global 3D hydrodynamic modeling of in-transit Ly{\alpha} absorption of GJ436b

Using a global 3D, fully self-consistent, multi-fluid hydrodynamic model, we simulate the escaping upper atmosphere of the warm Neptune GJ436b, driven by the stellar XUV radiation impact and gravitational forces and interacting with the stellar wind. Under the typical parameters of XUV flux and stellar wind plasma expected for GJ436, we calculate in-transit absorption in Ly{\alpha} and find that it is produced mostly by Energetic Neutral Atoms outside of the planetary Roche lobe, due to the resonant thermal line broadening. At the same time, the influence of radiation pressure has been shown to be insignificant. The modelled absorption is in good agreement with the observations and reveals such features as strong asymmetry between blue and red wings of the absorbed Ly{\alpha} line profile, deep transit depth in the high velocity blue part of the line reaching more than 70%, and the timing of early ingress. On the other hand, the model produces significantly deeper and longer egress than in observations, indicating that there might be other processes and factors, still not accounted, that affect the interaction between the planetary escaping material and the stellar wind. At the same time, it is possible that the observational data, collected in different measurement campaigns, are affected by strong variations of the stellar wind parameters between the visits, and therefore, they cannot be reproduced altogether with the single set of model parameters.


Introduction
A series of observations with the Hubble Space Telescope/ Space Telescope Imaging Spectrograph (HST/STIS; Kulow et al. 2014;Ehrenreich et al. 2015;Lavie et al. 2017) revealed that the warm Neptune GJ 436b has a very deep transit in the Lyα line, with up to 60% of the stellar flux in this line being absorbed. Moreover, this strong absorption takes place mostly in the blue wing of the line in the range of Doppler-shifted velocities of [−120, −40] km s −1 . The good signal-to-noise ratio clearly revealed the first time in visual ultraviolet (VUV) observations of close-orbit hot/warm exoplanets such details of the transit light curve as early ingress  and extended egress (Lavie et al. 2017). The quality of obtained data enables quantitative testing of existing theoretical concepts and numerical models proposed hitherto for GJ 436b and other similar closeorbit exoplanets.
The basic physical concept behind the measured absorption in the Lyα line greatly exceeding the optical transit depth and duration is related to the expanding hydrogen-dominated upper atmospheres of close-orbiting exoplanets. According to the energy-limited estimates (Lammer et al. 2003), the ionizing radiation of a host star leads at orbital distances <0.2au to the intensive thermal escape and mass loss of hydrogen-dominated upper atmospheres of hot gas giants. Since the first detection of an excess absorption in the Lyα of the hot Jupiter HD 209458b at the level of 10% (Vidal-Madjar et al. 2003), this concept has been rapidly developed, with continuously increasing complexity of numerical models. The first generation of 1D aeronomy codes (Yelle 2004;García Muñoz 2007;Koskinen et al. 2007) clarified the basic physics of the escaping upper atmosphere in the form of planetary wind (PW), which includes the X-ray and ultraviolet (XUV) heating, hydrogen plasma photochemistry, radiation cooling, and gravitational and thermal pressure forces. They helped to explain some of the in-transit spectral observations by the presence of expanded partially ionized upper atmospheres, which fill the Roche lobes of hot giant exoplanets, such as HD 209458b and HD 189733b (Ben-Jaffel 2007; Koskinen et al. 2007Koskinen et al. , 2010Ben-Jaffel & Sona Hosseini 2010). These expanding atmospheres were shown to be sufficiently dense to produce the absorption in Lyα owing to the natural line broadening mechanism.
However, the detection of absorption in the resonant lines of heavy elements such as O I, C II, and Si III Linsky et al. 2010) has shown that the absorbing material of planetary origin far beyond the Roche lobe has to be considered as well (Ben-Jaffel & Sona Hosseini 2010;Shaikhislamov et al. 2018a). The presence of a huge hydrogen corona also is a prerequisite for the explanation of strong in-transit Lyα absorption of GJ 436b. By this, there is another crucial factor, besides the Roche lobe effect, that has to be properly taken into account in the modeling of large-scale plasma dynamics around the close-orbit exoplanets-the stellar wind (SW) plasma. Self-consistent description of the escaping multicomponent PW flow, accelerated by the pressure gradients and stellar gravity, and its interaction with the SW required an upgrade of the first generation of 1D models to 2D and 3D ones. The corresponding effort has been undertaken in Shaikhislamov et al. (2016Shaikhislamov et al. ( , 2018aShaikhislamov et al. ( , 2018b and Khodachenko et al. (2017).
Outside the Roche lobe, under the conditions of escaping planetary upper atmospheric material (hydrogen and/or other elements), the dominating absorption mechanism is the thermal line broadening, provided by the resonant particles moving at the corresponding matching speeds. Thus, for example, for Lyα absorption at Doppler-shifted velocities up to a hundred kilometers per second, the presence of hydrogen energetic neutral atoms (ENAs) is required, as the atoms of PW are two slow and cold and they cannot contribute to this process. Two mechanisms for the generation of ENAs were proposed: (1) acceleration by the stellar radiation pressure, which in the case of hydrogen is provided by the stellar Lyα flux (Vidal-Madjar et al. 2003), and (2) charge exchange between fast stellar protons and planetary atoms (Holmström et al. 2008).
Taking the appropriate values for the parameters of PW, passing the Roche lobe, the generation and further dynamics of ENAs can be simulated with the kinetic Monte Carlo models, to interpret the transit observations. Such a 3D modeling approach, which covers the whole scale of the planet-star system, has been extensively applied to a variety of hot giant exoplanets, for which the absorption in Lyα and other resonant lines were measured. Having included initially only the radiation pressure force as the ENAs' generation mechanism (Vidal-Madjar et al. 2003, 2004Lecavelier des Etangs et al. 2004Bourrier et al. 2015), later, these models also incorporated the charge exchange process (Holmström et al. 2008;Kislyakova et al. 2014;Bourrier et al. 2016;Lavie et al. 2017). In particular, for the hot Neptune GJ 436b, considered in this paper, such modeling was able to reproduce not only the strongly asymmetric absorption profile and transit depth Bourrier et al. 2016) but also the observed early ingress and delayed egress of the transit light curves (Lavie et al. 2017). However, despite the good agreement of numerical Monte Carlo simulations of GJ 436b with observations, the inferred best-fit parameters of PW, used as the model input, appear in obvious disagreement with the predictions of aeronomy models for this exoplanet. The chemically complex multispecies 1D model by Loyd et al. (2017) and 3D simulations of the escaping hydrogen-helium upper atmosphere of GJ 436b by Shaikhislamov et al. (2018b) both show very similar results. Specifically, the self-consistently derived atmospheric mass-loss rate in these models is ∼3× 10 9 gs −1 , with the outflow velocity not exceeding 10 km s −1 , and the maximum temperature is ∼4500 K. On the other hand, the values required for the Monte Carlo modeling in order to fit the observations, derived in Bourrier et al. (2016) and Lavie et al. (2017), are 3×10 8 gs −1 and 70 km s −1 , respectively. To justify such discrepancy, Bourrier et al. (2016) and Lavie et al. (2017) assumed that such a low mass-loss rate is due to a very small heating efficiency by stellar XUV radiation, whereas the high outflow velocity of PW is due to the acceleration by hypothetical Alfvén waves.
In another recent simulation of the Lyα in-transit signatures of GJ 436b with the Monte Carlo model by Kislyakova et al. (2019), somewhat different conclusions regarding the parameters of escaping PW and ENAs' generation mechanism have been made. A good agreement of the modeling results with the observations in both the blue and red wings of the Lyα line was achieved, but at the essentially smaller velocity of the escaping planetary upper atmospheric material of only ∼8 km s −1 . Note that such a velocity value agrees well with the aeronomy simulations. Moreover, the modeling by Kislyakova et al. (2019) confirms the conclusion of the present paper that atomic hydrogen acceleration by the radiation pressure is insignificant for the observed in-transit Lyα signatures of GJ 436b, and the key role in the production of ENAs belongs to charge exchange.
There are several other important physical aspects that cannot be modeled with the Monte Carlo approach and that can be adequately described only by HD/MHD models. Among those are the effects of thermal pressure and shock waves. In the case of a supersonic interaction between the PW and SW, one can naturally expect the formation of shocks, or thermally compressed regions. Moreover, according to 1D, 2D, and 3D aeronomy modeling of several hot exoplanets (Koskinen et al. 2013;Shaikhislamov et al. 2016Shaikhislamov et al. , 2018aShaikhislamov et al. , 2018bKhodachenko et al. 2017;Loyd et al. 2017), the PW remains collisional up to rather far distances (of several tens of R p ) from the planet. The exobase is located usually far beyond the Roche lobe, as well as beyond the shocked region of the interacting PW and SW. This makes the hydrodynamic approach to the modeling of plasma environments of hot exoplanets a necessary and efficient way to understand the physical processes there and simulation of the related transit observations.
In the present paper, a fully self-consistent 3D hydrodynamic model is applied for the first time to simulate the in-transit Lyα line absorption features of GJ 436b as observed with the HST/ STIS. The multifluid aeronomic code of the model includes the hydrogen plasma photochemistry, the self-consistent stellar radiation energy input to the upper atmosphere of the planet that drives its expansion, and the effects of stellar and planetary gravity, as well as the SW plasma flow. The HD/MHD modeling of hot close-orbit exoplanets has been steadily progressing nowadays from 1D to 3D codes (Bisikalo et al. 2013;Tremblin & Chiang 2013;Owen & Adams 2014;Trammell et al. 2014;Khodachenko et al. 2015Khodachenko et al. , 2017Matsakos et al. 2015;Tripathi et al. 2015;Bisikalo & Cerenkov 2016;Shaikhislamov et al. 2016Shaikhislamov et al. , 2018aShaikhislamov et al. , 2018bErkaev et al. 2017). At the same time, the majority of existing 3D models have not yet reached the same level of physics-andchemistry complexity as that of the first generation of 1D aeronomy models, which would allow self-consistent simulation of the outflowing PW and its interaction with the surrounding SW. In that respect, the model presented in this paper is the first one that includes the aeronomy part, comparable to the existing 1D simulations, while covering, at the same time, the global 3D plasma environment of the whole stellar-planetary system. It is also the only model among the existing ones that incorporates the key processes considered in the Monte Carlo simulations, such as radiation pressure and charge exchange. Therefore, it can be directly compared to this kind of simulation as well. The proposed 3D hydrodynamic model not only reveals the parameters of PW, stellar XUV flux, and SW, at which good agreement of the calculated (synthetic) and experimental line absorption profiles and transit light curves is achieved, but also describes how different physical processes affect the observations.
Our previous paper on 3D simulation of GJ 436b (Shaikhislamov et al. 2018b) was dedicated to upper atmosphere escape and did not take the stellar plasma into account. Its primary goal was to compare the results of the selfconsistent 3D modeling to those of the previous 1D one and to investigate the role of helium abundance. It was found, in particular, that while the total influence of helium abundance on the mass loss is relatively small, the presence of helium affects the structuring of the nearby planetary plasmasphere. Because of its higher mass, the helium component decreases the gravitational scale height of the planetary atmosphere, resulting in a sharper decrease of density with height. This shifts the H 2 dissociation front closer to the planet, and for GJ 436b we saw two qualitatively different kinds of thermosphere -one with an extended molecular hydrogen envelope, and another (for He/H<0.1) with a restricted H 2 area, closely localized around the planet. The presence of helium also affects the photoionization of hydrogen. While the study of helium's influence on the transit observations of GJ 436b is not the major goal of the present paper, we nevertheless take it into account in the modeling.
It has to be noted that the results of 3D modeling of GJ 436b in Shaikhislamov et al. (2018b) are in good agreement with the chemically more complex aeronomy simulation of the planet by Loyd et al. (2017), which was just a 1D model but included the oxygen and carbon components in the atmosphere. In particular, such features as temperature maximum, escaping PW velocity, H 2 half-dissociation height, and the mass-loss rate, predicted by both models, are quantitatively similar, and the differences of up to 25% can be attributed to the specifics of 3D modeling, not reproduced in 1D. The model in Shaikhislamov et al. (2018b) confirmed the conclusion of the previous estimates that the exosphere of GJ 436b is indeed very extended and obscures most of the stellar disk, even in spite of the high impact parameter of the orbit, making the planet transit close to the disk edge. Among the newly simulated essentially 3D features are the Coriolis twisting of the PW streams in the noninertial planetbased reference frame and their compression toward the ecliptic plane by the stellar gravity, which in the tidally locked rotating system is not compensated by the centrifugal force in the polar directions. This effect results in the formation of sharp plasmasphere boundaries at distances of (10-20)R p below and above the ecliptic plane.
The modeling results reported in the present paper demonstrate that the huge weakly ionized plasmasphere built around GJ 436b causes the strong absorption in the Lyα line, which can be as large as 70%-90%. By this, the absorption at high Doppler-shifted velocities is produced by the ENAs generated by charge exchange between the planetary atoms and fast protons of the SW (in case the latter is prominently present). Therefore, the Lyα absorption is strongly asymmetric and shifted to the blue wing of the line. The interaction between the outflowing PW and SW is characterized by the presence of a clearly pronounced shocked, or compressed, region. This is the region where the ENAs are produced because the planetary atoms penetrate there through the ionopause and charge exchange with the compressed proton fluid of the SW. The sharp early ingress starts at the position of the shock ahead of GJ 436b, which depends on the total pressure balance between the PW and SW. The egress, on the other hand, tends to be long-lasting in duration as a part of the escaping PW material trails behind the planet. With the dedicated model runs we determine the parameters of PW and SW, as well as the stellar XUV flux, at which good agreement of the calculated and the measured Lyα absorption profiles and the corresponding transit light curves is achieved, and describe how different physical processes affect the observations. It has been also found that at the Lyα radiation flux expected for GJ 436, based on actual measurements and reconstruction of the line core , the radiation pressure produces too small of an effect in the acceleration of hydrogen atoms to be seen in observations. In this respect it should be noted that the acceleration of the escaping upper atmospheric planetary atoms by the radiation pressure has been recently reanalyzed with simple analytical estimations and hydrodynamic modeling. By this, the coupling between protons and hydrogen atoms, is crucial for the acceleration of the last, which was not considered in the Monte Carlo approach, was properly taken into account (Shaikhislamov et al. 2016;Cerenkov et al. 2017;Khodachenko et al. 2017;Debrecht et al. 2019). In the present paper we recall further the question about the role of the radiation pressure. In the Appendix we discuss why the radiation pressure force has so little influence on the Lyα absorption features of GJ 436b.
The Lyα absorption modeled here appears in good agreement with the measurements of major observational features, such as strong asymmetry between blue and red wings of the absorbed Lyα line and deep (>70%) in-transit depth in the high-velocity blue part of the line, as well as the early ingress and extended egress in parts of the transit light curve. However, the exact simultaneous fitting of all these observed features appears to be difficult. The transit light curve of GJ 436b, built on the basis of all observational visits, has a complex egress part, which consists of a sharp initial drop followed by a gradual decrease, as described in Section 2. This complexity can be either a result of action of specific physical processes leading to such an ambivalent behavior of the light curve or a consequence of overlaying in one light curve of the measurements performed at different conditions of different epochs, for which the data were obtained. The major discrepancy here is that the modeling shows the absorption maximum at 1-3hr after midtransit, followed by a long egress part, which is significantly deeper than the one in the observations. The egress decline can be made sharper by switching to a denser and faster SW, which blows away the trailing planetary tail. However, in this case the early ingress will also be affected, and the shallow and long part of the egress will disappear.
It is worth noting that there are physical aspects, still not included in our model, that certainly can affect (qualitatively and quantitatively) the interaction between the PW and SW. These are, first of all, the magnetic fields of the star and SW, on one hand, and the planetary magnetic field, on the other. The presence of such large-scale (stellar and planetary) magnetic fields modifies the geometry of the regions where the ENAs are generated, thus resulting in an additional physical factor that influences the ingress and egress dynamics of the transit light curve. Nevertheless, it should not be forgotten that the available observational data for none of the single visits properly cover the whole transit light curve with all of its phases, i.e., the ingress, midtransit, and extended egress parts. This fact leaves the possibility of an influence of the varying stellar radiation and plasma conditions from visit to visit on the observed behavior of the transit light curve.
The paper is organized as follows. In Section 2, the available relevant observational data on GJ 436b are briefly reviewed. Section 3 describes the model details. In Section 4, the results of simulations for the different SW conditions and XUV radiation fluxes, as well as other modeling parameters, are reported, and the role of the radiation pressure is quantitatively investigated. Section 5 presents the discussion of the obtained results and main conclusions.

Summary of the Lyα Absorption Measurements for GJ 436b
In this section, we briefly summarize the observational data regarding Lyα absorption during the transits of GJ 436b. There are in total eight data measurement campaigns, so-called visits, when the Lyα flux of GJ 436b was recorded with HST, as described in Ehrenreich et al. (2015), Bourrier et al. (2016), and Lavie et al. (2017). Visits 0, 4, and 5, which cover the time interval of 26-37hr after the midtransit, are used to derive the out-of-transit Lyα flux of GJ 436, as well as to infer its shortterm and long-term variability. In-transit light curves are constructed with data from two other sets of visits. Visits 1-3 cover about±3.3hr around the midtransit and show more or less the same details of the transit light curve, namely, a sharp ingress, starting at about −2hr; deep maximum of absorption, with a depth of up to 70%, continued for about 0.75hr; and a sharp fall at the interval of [1, 3]hr. Here and further on the phase of transit light curve is referred to as a time relative to the ephemerid midtransit time of the planet, with the positive and negative values corresponding to the phases before and after the midtransit, respectively. Visits 6 and 7 were aimed at measuring the duration of egress. They cover the interval of [4,9]hr after the midtransit and, as an out-of-transit reference, the interval well before the midtransit [−8.5, −3.3]. Combined together, visits 6 and 7 show a different transit light curve, as compared to one revealed at visits 1−3, which is much shallower (only 25%) but smoother and more extended in time. This difference in data can be interpreted either as evidence of two independent transit scenarios-(1) sharp and deep, and (2) long and shallow-or as a combined scenario with a sharp increase in the absorption depth during the interval [−2, 2]hr, when the planet crosses the midtransit, followed then by a long shallow egress.
For the comparison with simulations, we use the data of all visits from Lavie et al. (2017). The left panel of Figure 1 shows the transit light curve integrated over the blue wing of the Lyα line in the interval of Doppler-shifted velocities [−120, −40] km s −1 , composed of all visits. The right panel of Figure 1 shows the absorption profiles over the whole Lyα line at the midtransit (0hr) and post-transit (sum of measurements at 4.2 and 5.8hr) averaged over all visits, as well as the out-of-transit profile for reference.
The absorption depth, averaged over the blue wing interval [−120, −40] km s −1 of the Lyα line, reaches in some visits up to 70% of the stellar flux at the midtransit and 32% at the posttransit. In the red wing interval [30, 110] km s −1 of the line it is 16% at the midtransit and 22% at the post-transit. Thus, the Lyα absorption of GJ 436b at the midtransit is very asymmetric, being much stronger in the blue wing of the line, while in the post-transit it is more symmetric and reaches moderate values in both the blue and red wing intervals. These features of the Lyα absorption are related to the dynamical conditions realized in the interacting PW and SW around GJ 436b. The detailed modeling of this interaction and the related distribution and motion of the absorbing hydrogen shed light on the key physical processes in the system, transit scenarios, and dominating absorption mechanisms.
Further on, for the purpose of comparison of the modeled and observed transit light curves and the Lyα absorption profiles, we use the same observational data as those displayed in Figure 1, if not specified otherwise.

The Model
A 3D hydrodynamic multifluid model, used here, further extends the 3D model of the expanding hydrogen-dominated upper atmosphere of GJ 436b, employed in Shaikhislamov et al. (2018b). The present simulations of GJ 436b incorporate also the stellar plasma flow. Similar to Shaikhislamov et al. (2018b), we consider, in addition to the dominating hydrogen of the planetary and stellar origin, the helium (He) component, which is taken with an abundance He/H=0.1 (if not specified otherwise). Therefore, among the major species of the simulated fluid are the hydrogen and helium particles and electrons. The model code solves numerically the hydrodynamic equations of continuity, momentum, and energy for all species (marked with index j) of the simulated multicomponent flow, written in the following form (Shaikhislamov et al. 2016): Equations (1)-(3) are solved in a noninertial spherical frame of reference fixed at the planet center, which orbits together with the planet and rotates with the same rate Ω so that one of the axes is continuously pointed to the star. The Z-axis is directed perpendicular to the ecliptic plane. This is a so-called tidally locked frame of reference. Though, in the general case, the planet itself can rotate around its axis with its own (different) rate, in this work GJ 436b is taken to be also tidally locked to the star. Equations (1)- (3) properly account for the noninertial terms, i.e., the generalized gravity potential U and Coriolis force. The particular expressions for these terms are well known, and they were also given in our previous papers, e.g., in Khodachenko et al. (2015), Shaikhislamov et al. (2016), and Dwivedi et al. (2019). Other basic notations used in Equations (1)-(3) are the standard ones, with m j , z j , and n j standing for the mass, charge, and number density of species j, respectively (for the neutrals z j is taken to be zero). Note that the electrons are considered as a massless fluid, propagated under the action of its pressure gradient and Coulomb force. Magnetic field is not considered, while the electric field is expressed using the momentum equation for electrons. Altogether, this results in the appearance of a specific term in Equation (2) proportional to the electron pressure gradient. The atomic hydrogen component is not originally present in the SW. It is produced by charge exchange between the SW protons and planetary atomic hydrogen. The minor species, which in the considered here case are represented by helium, are also described as separate fluids with the corresponding energy, momentum, and ionization-recombination equations. An account of molecular hydrogen (H 2 ) and the corresponding molecular ions + H 2 and + H 3 , as, e.g., in Shaikhislamov et al. (2018aShaikhislamov et al. ( , 2018b and Khodachenko et al. (2017), enables more accurate treatment of the inner regions of the planetary thermosphere.
The basic physical parameters of the model and the simulation results are scaled in units of the characteristic values of the problem: temperature, in T 0 =10 4 K; speed, in the corresponding proton thermal speed, V 0 =9.07 km s −1 ; distance, in radius R p =2.5×10 9 cm of GJ 436b; and time, in τ=R p /V 0 ∼0.75hr. These units are also used for scaling of axes in the appropriate figures below.
The main processes, responsible for the transformation between neutral and ionized particles, are photoionization, electron impact ionization, and dielectronic recombination. These are included in the term XUV,j j XUV,j XUV j e Te ion,j j e Te rec,j in the continuity Equation (1) and are applied for all species with indexes j and j+, where the latter indicates ions of the particles of sort j. Photoionization also results in a strong heating of the planetary atmospheric material by the produced photoelectrons. The corresponding term W XUV,j in the energy Equation (3) (addressed in Trammell et al. 2011;Shaikhislamov et al. 2014;Khodachenko et al. 2015) is derived by integration of the stellar XUV spectrum. For the M dwarf GJ 436 we use the XUV spectrum compiled by the MUSCLES survey (France et al. 2016), which is based on actually measured fluxes in VUV and X-ray bands. The corresponding integrated XUV flux (λ<91.2nm) is F XUV ∼0.86 erg s −1 cm −2 at 1au. We also applied the regression model by Arkhypov et al. (2018) based on the statistical analysis of available stellar parameters (rotation periods and effective temperatures) and observed activity rates, which predicts for GJ 436 the flux F XUV ∼2 erg cm −2 s −1 at 1au. Note that the same values of 1-2 erg cm −2 s −1 at 1au are estimated also in Bourrier et al. (2016) and Lavie et al. (2017). The XUV photons ionize hydrogen atoms, as well as H 2 and He components, according to the wavelength-dependent cross sections. The model assumes that the energy, released in the form of photoelectrons, is rapidly and equally redistributed between all locally present particles with an efficiency of 50%. This is a commonly used assumption, which we adopted on the basis of qualitative analysis (Shaikhislamov et al. 2014). The attenuation of the XUV flux inside the planetary atmosphere is calculated for each spectral bin according to the wavelength-dependent cross sections. The heating term includes also energy losses due to excitation and ionization of hydrogen atoms, and in a simplified form it can be written as follows: Of certain importance are also the energy losses due to infrared radiation of the + H 3 molecule with a rate C H3+ calculated by Miller et al. (2013) and those due to dissociation of the H 2 molecule, W cool =(γ −1)(n H3+ × C H3+ +n H2 × E diss × C diss )/N tot . Here N tot is the total number density of all particles, including electrons. The specific photochemistry reactions and their rates used in the model, including those of the dissociation, are listed in Khodachenko et al. (2015).
Another kind of exchange between the considered particle populations are the resonant charge exchange collisions. Indeed, charge exchange has a typical cross section of σ exh = 6×10 −15 cm 2 at low energies, which is an order of magnitude larger than the elastic collision cross section. The charge exchange process is included in the term N exh,j =n j+ n j σ exh v j+,j in Equation (1). Also, when the planetary atoms and protons slip relative to each other, because they have different thermal pressure profiles, and the protons additionally feel the electron pressure while the atoms do not, the charge exchange between them leads to their velocity and temperature interchange. We describe this process with the collision rate coefficients u C ji and C ji T , where the upper index indicates the physical quantity being interchanged. For example, in the momentum Equation ( H H 2 depends in general on thermal speeds and relative velocities of the interacting fluids (in the considered example-of protons and hydrogen atoms). Besides charge exchange, ordinary elastic and Coulomb collisions (∝ n i σ ji v j,i ) are included in the u C ji and C ji T terms.
For the typical parameters of the planetary plasmasphere, Coulomb collisions effectively couple the ionized species between each other and with electrons. For example, at T<10 4 K and n H+ >10 6 cm −3 the collisional equalization time for temperature and momentum for protons is less than 2.4 s (Braginskii 1965). The analogous value for He + ions is the factor of (m He /m H ) 1/2 ∼ 1.8 higher (i.e., 4.4s). These times are several orders of magnitude less than the typical gasdynamic timescale of the problem τ=R p /V 0 ∼3×10 3 s. There is also another physical reason for the strong coupling of charged particles within the typical spatial scale of the problem, which is about R p (∼10 9 cm). The chaotic/sporadic magnetic field, which is naturally present in the considered interacting plasma flows, affects the relative motion of ions, so that they become coupled via the magnetic field owing to the Lorentz force and exchange momentum on a timescale of the Larmor period. Even for a very weak magnetic field of ∼1 nT the proton gyroradius remains smaller than the typical spatial scale. Thus, for the expected sporadic magnetic fields of the order of 10 3 nT, the charged particles can be safely assumed to be effectively coupled with a sufficiently small mean free path (of about a gyroradius). For the same reason, the protons of hot and rarefied SW can be treated as strongly coupled ones, even in spite of the fact that Coulomb collisions are negligible there. By this, the interpenetration of planetary and stellar protons in the region of counterstreaming PW and SW plasma flows would be microscopically small. Therefore, a sharp ionopause and bow shock (in the case of a supersonic relative velocity) have to be formed in the region where the expanding PW meets the SW. The location of ionopause would be defined by the total pressure balance of the counterstreaming PW and SW flows. In view of the strong coupling of the charged components of the simulated material flows, there is no need to calculate separately the dynamics of every charged species, and we assume in the simulations that all of them (i.e., ) have the same temperature and velocity. On the other hand, the temperature and velocity of each neutral component are calculated individually by solving the corresponding energy and momentum equations. Note that the neutral hydrogen atoms are to a certain extent coupled to the main material flow also by elastic atomic-atomic and proton-atomic collisions. With a typical cross section σ HH ≈σ H+H ≈ 10 −16 cm 2 , the mean free path at a density of 10 6 cm −3 is comparable to R p . At the same time, the above-discussed charge exchange ensures more efficient coupling between hydrogen atoms and protons (Shaikhislamov et al. 2016;Khodachenko et al. 2017).
The total fluid velocity at the planetary surface is taken to be zero, whereas at the stellar surface the normal component of the plasma flow velocity is defined by the polytropic solution of the SW model, described below. The tangential component of SW velocity is connected to the stellar surface rotation in the chosen noninertial frame of reference. It is worth noting that those velocities at the stellar surface are small, as compared to the typical SW velocity. To keep the number of points that operated in the numerical code, tractable for processing, while enabling sufficient resolution of the model, a nonuniform radial mesh was used with the grid step increasing proportionally with distance from the planet surface. This allows resolving of the highly stratified upper atmosphere of the planet, where the required grid step is as small as Δr=R p /400. Taking into account that the planet-star distance is usually tens, or even hundreds, of R p , the corresponding grid step of ∼R p , at the distances of several tens of planetary radii, R p , ensures sufficient spatial resolution of the model also in these moderately remote regions. The applied radial spacing in the spherical coordinate system also keeps the same resolution in all three dimensions, if the azimuthal and latitudinal mesh is chosen so that Δj≈Δθ≈Δr/r. Further details regarding the numerical code and the applied grid are provided in Appendix B. For the analysis of the simulation results, a related Cartesian coordinate system with the X-axis, pointed from the planet to the star, is also used further.
For the self-consistent calculations on the scale of the whole star-planet system, we incorporate in the same code the SW plasma dynamics as well. Besides the upper planetary atmosphere, the stellar corona appears as another boundary of the simulation domain, at which the corresponding coronal values of the plasma parameters are fixed. Outside this boundary, taken for simplicity at the stellar surface with the radius R star , the SW, i.e., the proton fluid, is calculated by the same code as that used for the protons of PW. The simulation of SW is a complex problem, which has been extensively tackled in recent decades (see, e.g., Usmanov et al. 2011 and references therein). Its detailed treatment, especially the processes of the SW heating and acceleration up to supersonic velocities (still not yet fully understood), appears beyond the scope of the present paper. A common approach of the majority of the global astrophysical codes, adopted for the simulation of exoplanetary environments, is to model the SW with a simplified approach, using a polytropic or even an isothermal specific heat ratio 1γ p 1.2 (e.g., in Bisikalo et al. 2013;Matsakos et al. 2015;Christie et al. 2016;Daley-Yates & Stevens 2018). Any polytropic index lower than the adiabatic one, γ a =5/3, introduces in the energy equation an effective source term that replaces the variety of complex SW energizing processes. However, simulation with a nonadiabatic index makes unrealistic the treatment of shocks, developed in the region of the colliding SW and PW flows. To avoid this difficulty, we introduce in the energy equation an empirical heating source to drive the SW. The corresponding source term is found from a semianalytical solution of the 1D polytropic Parker-like model, which yields for a given stellar gravity the base temperature T cor , the asymptotic supersonic velocity ¥ V sw, , a unique value of γ p , and the radial profiles of all SW parameters (see, e.g., in Keppens & Goedbloed 1999). To fix the SW plasma density at the stellar surface boundary, we use an integral characteristic parameter-stellar mass-loss rate, ¢ M sw . The heating term, therefore, is found as Here T p and V p are obtained from the polytropic solution, and the distance R is measured from the center of star. It was specially checked in the simulations that without PW the solution obtained for SW is in a good agreement with the analytical model. Note that Equation (5) specifies a spatially distributed, continuous, and stationary-in-time heat source, which is defined by the above-specified model parameters only.
To avoid unphysical artifacts, the heating source is switched off in the regions where a significant amount of particles of planetary origin penetrate into the SW. The modeling approach, described above, allows a physically self-consistent global simulation of the generation of supersonic PW and SW, as well as their interaction. For the sake of generality and for comparison with other models, we also include in the code the radiation pressure force, acting on the hydrogen atoms owing to the Lyα flux. The reconstructed Lyα line profile, based on the actual measurements for GJ 436, can be found in Bourrier et al. (2015). The total flux of ≈1 erg cm 2 s −1 at the reference distance of 1au generates at the center of the emission line a force equal to ≈0.7 of the stellar gravity pull. The radial dependence of the Lyα flux and the particle instant velocity relative to the line center are properly taken into account during the calculation of the acceleration of the hydrogen atoms. An important factor is the Lyα flux self-shielding. This we compute by Voigt convolution of the Lorentz line shape with the natural width and the Maxwellian distribution of atoms with the temperature and density calculated by the hydrodynamic model. The convolution integral is expressed via the empirical formula from Tasitsiomi (2006), as described in Shaikhislamov et al. (2016) and Khodachenko et al. (2017). The self-shielding is calculated in the code in 30 km s −1 bins in the range of±120 km s −1 .
The synthetic absorption in the Lyα line by the transiting GJ 436b is calculated with a dedicated data processing program, on the basis of the numerical simulation results, as follows. We use the available observational data on the spectrally resolved out-of-transit and in-transit fluxes, F out (λ) and F in (λ), respectively. By this, the integrated absorption in a particular part [λ 1 , λ 2 ] of the line profile is calculated as The physical parameters of the dynamic planetary environment obtained in the course of the simulations are used to calculate the absorption A abs (λ) of a given atomic line (in the present study, of the Lyα line) in dependence on λ, while the emissionline profile itself is not specified. Since, by definition, F in (λ)=[1−A abs (λ)]×F out (λ), the integrated synthetic absorption is obtained, using Equation (6), by averaging of A abs (λ) over the emission-line profile F out (λ) measured in observations: ( ) A summary of the parameters of the star-planet system used in the modeling is given in Table 1. Further on, the following parameters, mostly based on the analogy with the Sun and general physical reasoning, will be assumed as the so-called "standard" ones: upper atmospheric temperature and pressure of GJ 436b, i.e., at the inner boundary of the simulation domain (at the conventional planet surface), T base =750 K and P base = 0.05 bars, respectively; helium abundance He/H=0.1; ionizing radiation of GJ 436 at 1au F XUV =0.86 erg cm −2 s −1 ; stellar coronal temperature T cor =2×10 6 K; the terminal SW velocity = ¥ V 400 sw, km s −1 ; and the stellar mass-loss rate ¢ = M sẃ 2.5 10 11 gs −1 . Note that the assumed stellar radiation and mass-loss rate are about 5-10 times less than those for the Sun, which is mostly due to the smaller size of GJ 436.

Simulation Results
About 100 simulations were performed with the varying of most important parameters of the system in relatively broad ranges, specifically with the stellar XUV flux scaled to 1au 13 gs −1 . The total mass loss of the star is an unconstrained parameter, so it was varied from the small values, at which SW does not affect the planetary atmospheric material outflow, up to extremely large figures, exceeding those of the quiet Sun by more than an order of magnitude.
Besides these main parameters of the system, there are several others that to a lesser extent also affect the absorption in Lyα. Among those is the stellar coronal temperature, T cor , which may influence the formation of the shock, and which was varied for different modeling runs in the range from 1.5 to 4.5 MK. Also, the changes of cooling rate due to the + H 3 molecule were investigated, which, being varied by an order of magnitude, correspondingly affects the planetary atmospheric mass-loss rate up to ∼50%. Note that this cooling rate so far has not been directly measured. It is based on complicated theoretical calculations and estimates (Miller et al. 2013). With a series of dedicated simulation runs, it has also been found that variation by a factor of two of such parameters of the upper planetary atmosphere (i.e., of the inner boundary conditions) as T base , P base , and helium abundance does not produce any significant difference in the modeling results.

Weak Stellar Wind and the Role of Radiation Pressure Force
First, we consider the so-called "captured by the star" regime of the PW and SW interaction (Shaikhislamov et al. 2016). It corresponds to the case of a sufficiently weak SW, when the escaping PW material is captured by the stellar gravity and accumulates around the star. For the standard parameters of the system specified above, the captured by the star regime takes place for a small stellar mass-loss rate ¢~Ḿ 5 10 sw 9 gs −1 , which is ∼500 times less than that of the present Sun. Note that at four times higher ¢ =Ḿ 2 sw 10 10 gs −1 , the SW is already strong enough (i.e., dense and fast) to blow away the escaping PW material, which corresponds to the "blown by the wind" regime of the PW and SW interaction according to the terminology used in Shaikhislamov et al. (2016). Figure 2 shows the color plots of H and H+ density distributions and the typical streamlines, achieved after about 25 orbital revolutions of the planet after the start of the simulation run.
One can see that the PW stream spirals and quickly (in several tens of revolutions) falls on the star. The calculation of long-term accretion and/or accumulation of the material , ¢ =Ḿ 5 10 sw 9 gs −1 . Here and in other similar figures below the distance is scaled in units of GJ 436b radius, R p . The planet is located at the center of the coordinate system. The star is located to the right from the planet at X=158 and is marked with the scaled black circle. Black lines in both panels show the streamlines of the material flow.  Figure 2. The planet is located at the center of the coordinate system. The star is located to the right from the planet at X=158 and is marked with the scaled black circle. Black lines show the streamlines of the atomic hydrogen flow.
around the star requires the inclusion of more sophisticated physics than that in the present version of the code. The performed simulations show that the captured by the star regime of the PW and SW interaction is unlikely for GJ 436b, because under the typical standard parameters of the system it takes place only at extremely low stellar mass-loss rates. Therefore, in further modeling and the study of the Lyα absorption features of GJ 436b only the blown by the wind regime is considered. At the same time, the captured by the star regime is the most suitable one for the demonstration of the role of radiation pressure.
The left panel of Figure 3 shows the same plots as those in Figure 2 (left panel), but obtained with the inclusion of the radiation pressure force, produced by the stellar Lyα flux of F Lyα =2 erg cm −2 s −1 (scaled to 1au). At this value the acceleration of hydrogen atoms due to the radiation pressure force is about the same as that caused by the stellar gravity. Although the PW material also accumulates around the star as in the previous case, one can clearly see that the hydrogen atoms are pushed a bit farther away from the star and the spiraling of the PW stream lasts significantly longer than in the simulations without radiation pressure. The extended helmetshape structure, seen in the meridional cut across the ecliptic plane in the right panel of Figure 3, is produced mainly as a result of the radiation pressure force, and it is practically absent without it. The latter is clearly indicated also by the difference between green halos (n a level ∼10 2 cm −3 ) in the ecliptic plane in the left panels of Figures 2 and 3.
For the quantitative characterization of the role of the radiation pressure force in terms of its measurable manifestation in the Lyα absorption profiles, these, as well as the corresponding transit light curves, were calculated for both cases (i.e., with and without inclusion of the radiation pressure). The results are shown in Figure 4. In the case of weak SW and without the radiation pressure force, the absorption in the high-velocity blue wing of the Lyα line is just a few percent (as expected). Note that, besides the spike at the midtransit, the absorption in this case takes place over an extended time interval before and after the transit. This happens as a result of formation of an extended belt of atomic hydrogen around the star supplied mainly by the escaping PW material. The inclusion of radiation pressure force does not make a significant difference, although additional acceleration of atoms results in a pushing of the accreting PW stream away from the star. The consequences of this effect can be seen in the left panel of Figure 4 before the transit (i.e., ahead of the planet). The Lyα absorption profiles provided in the right panel of Figure 4 show that the velocities of hydrogen atoms are shifted, due to the action of the radiation pressure force, by about −20 km s −1 .
The effect of atomic hydrogen acceleration by the radiation pressure force is further demonstrated in Figure 5 in the velocity profiles along the planet-star line (Z=0, Y=0) in the star-centered coordinate system. As can be seen, if the radiation pressure is not included, the hydrogen atoms move toward the star, while being accelerated by the stellar gravity up to velocities of ∼100 km s −1 . The radiation pressure force, as expected from the estimations, nearly balances the stellar gravity, so that the hydrogen atoms do not fall to the star, except for the small region at R<50, where atoms are coupled to the protons and a complex mixing of the PW and SW material takes place. Note also that far from the planet, where  Figure 2) and with (also in Figure 3) the inclusion of the radiation pressure force, respectively. Figure 5. Velocity of hydrogen atoms along the planet-star line (Z=0, Y=0) in the star-centered coordinate system. The distance is scaled in units of GJ 436b radius, R p . The planet is located at X=−158. Positive/negative values of the velocity mean the motion away/toward the star. Dashed and solid lines correspond to the cases without (as in Figure 2) and with (as in Figure 3) the inclusion of the radiation pressure force, respectively. the material is rarified and hydrogen atoms become uncoupled from the protons, the streamlines of hydrogen atoms (see in Figure 3, left panel) are almost circular, as, being additionally supported by the radiation pressure force, atoms do not fall toward the star.
Finally, one of the reasons why the effect of the radiation pressure force is not prominent in the case of GJ 436b under study is the self-shielding. Figure 6 shows the optical depth for the Lyα photons ( ) as a function of distance along the line of sight, passing parallel to the X-axis and sufficiently close (Y=−20, Z=0) to the planet, for a set of different Doppler-shifted velocities V=c·(1-λ/ λ o ). As can be concluded from the plots in Figure 6, the central part of the Lyα line in the range [−30, 30] km s −1 is strongly absorbed in the close vicinity of the planet, e.g., for | | < Y 20. This means that close to the planet the acceleration of slow atoms (<30 km s −1 ) cannot take place because of the Lyα shielding.
Altogether, the numerical simulations and tests, performed with the inclusion of all the key physics, demonstrate that the radiation pressure force produces a relatively small effect on the Lyα absorption for GJ 436b. Therefore, for simplicity's sake, it will not be included in further modeling calculations. Nevertheless, to stay on the safe side with this kind of simplification, in several cases considered below we continued the comparison of the simulation results with and without the radiation pressure force, finding always very little difference between them. The question on the role of the radiation pressure force in the Lyα absorption by GJ 436b is further discussed in Appendix A, where the results of numerical simulations are supported and justified with a theoretical analysis and analytic estimates.

Material Escape on GJ 436b; Interaction with SW-Modeling under Different Conditions
The parameters of the discussed further simulation runs are listed in Table 2. For referencing purposes, a certain number (1...14) is assigned to each set of parameters. Figures 7 and 8 show the color plots of atomic hydrogen and proton density distributions, obtained in the modeling run with the parameter set No. 1. They reveal that the SW flow sweeps away the escaping PW material and redirects it to the trailing tail. A compressed layer is formed in front of the planet, well seen in the ENAs and temperature distributions (Figures 7 and 8, right panels), as well as in the temperature plot along the planet-star line in Figure 9.
However, due to the high temperature of SW plasma, the relative SW and PW velocity is not enough to generate a strong shock wave (the thermal speed in the SW plasma at the orbit of the planet is ∼130 km s −1 ). The presence of a weak shock is seen in Figures 7 and 8 as a slight discontinuity of the proton streamlines at the bow shock. The ENAs are mainly generated in the region between the ionopause (∼20R p ) and the bow shock (∼40R p ), both clearly seen in Figure 8 (right panel) as sharp boundaries in the ENAs' density color plot. The neutral atoms of planetary origin, which approach closely to the ionopause, penetrate into the shocked region, where they become uncoupled from the protons. The structure of different regions along the planet-star line (i.e., at Z=0, Y=0), from the planetary upper atmosphere up to the star, is shown in more detail in Figure 9. All these features have been already addressed and discussed with respect to our previous 2D simulations (Shaikhislamov et al. 2016;Khodachenko et al. 2017). They are qualitatively similar to those revealed in the 2D case during the "captured by the star" regime of SW and PW interaction, though quantitatively there are differences, related to the specifics of 3D, that cannot be reproduced in 2D modeling.
Next, we describe the absorption in Lyα derived in the simulation run with the parameter set No. 1. Figure 10 shows the transit light curve in the blue wing of the Lyα line and the modeled full line profiles at the midtransit (t=0 hr) and posttransit (t=5 hr), respectively. For comparison, the measured data from Figure 1 (according to Lavie et al. 2017) are reproduced as well.
As can be seen in Figure 10, our self-consistent 3D modeling of the dynamical environment of GJ 436b reproduces the main observational feature, namely, the strong absorption in the blue wing of the Lyα line. Under the conditions of the modeling parameter set No. 1, it reaches almost 70%. We note that these parameters have not been specifically fitted, but rather assumed as the most relevant and typical ones. The simulated Lyα absorption is strongly asymmetric, and in the red wing of the line it does not exceed 15%. This is because the absorption at velocities >30 km s −1 is to a significant extent due to the ENAs, created by charge exchange between the planetary atoms and SW protons. This fact is illustrated by the relative closeness of the blue and violet lines in the left panel of Figure 10. Note that the modeled transit light curve possesses quite an extended egress part that gradually decreases for as long as 10hr after the midtransit. This duration corresponds to about 1/6 of the whole orbit of GJ 436b. Such an effect is caused by the long trailing planetary tail. The line profiles in the right panel of Figure 10 show that the absorption takes place mostly in the blue wing and extends up to −150 km s −1 , which is in good agreement with the observations. A view of the geometry of the Lyα absorbing region is presented in Figure 11, where the color map of the line-of-sight absorption, expressed via the optical depth as 1−exp(−τ) and Figure 6. Logarithm of the integral optical depth for the Lyα photons versus distance from the planet, scaled in units of R p , along the line-of-sight parallel to the X-axis, for a set of different Doppler-shifted velocities. Negative velocity values correspond to the motion away from the star. The integration goes along the X-axis starting from the star position at X=158, at fixed Y=−20 and Z=0 in the planet reference frame.

Table 2
The Parameter Sets for the Discussed Simulations Notes.The XUV column contains the stellar flux in the range of wavelengths 1 nm<λ<91.2 nm, expressed in erg cm −2 s −1 and scaled to 1au. The mass-loss rates of the planet ¢ M pw and the star ¢ M sw are expressed in 10 10 gs −1 . If not specified otherwise, other parameters of the simulations are P base =0.05 bars, T base =750 K, and He/H=0.1. Index "swp" denotes SW parameters at the orbit of the planet. The  averaged over the blue wing of the Lyα line, is shown. Integration of this distribution over the whole stellar disk gives the absorption ≈60%, which corresponds the midtransit point (t=0 hr) on the transit light curve in Figure 10 (left panel).
One can clearly see in Figure 11 that the strongest absorption comes from the regions associated with the shock, located relatively far from the planet (up to ∼11R p ). The compressed layer between the planetary ionopause and bow shock, where the ENAs are generated, covers most of the central part of the stellar disk. Because of that, the grazing trajectory of the planet, as seen by an observer (i.e., from Earth), is actually favorable for the absorption. If the planet would cross the stellar disk close to the star's center, the maximum absorption would be only 43%, and the egress part of the transit light curve would be twice shorter in duration. The details of the plasma parameters' distribution are shown in Figure 12, where the profiles of proton, atomic hydrogen, and ENA densities along the line of sight, parallel to the X-axis, are presented for the fixed Z=11, and Y=0 (the planet-centered coordinate system is used here). One can see that along this particular line of sight, which in the observer reference frame in Figure 11 has the coordinates Z O =2, Y O =0, the major part of absorption should happen in the interval of X= [−30, 20], mostly populated by the atomic hydrogen fraction. Since the considered line of sight crosses the shocked region, the corresponding temperature and velocity there are relatively high. Figure 12 also shows the component of velocity, V s , along the line of sight, which is positive everywhere (directed to the observer), except a small region near the planet. It confirms that the neutrals move away from the star with speeds corresponding to the Doppler-shifted velocities of the blue wing of the Lyα line, from −10 to −200 km s −1 (in the planet reference plane).
The next series of simulation runs with the parameter sets No. 1-6 demonstrates how the modeling results are changed owing to a factor of two variation of the major two parameters of the problem, specifically the stellar XUV flux F XUV and the mass-loss rate of the star ¢ M sw . Note that the intensity of SW (i.e., the mass flux) in all these cases is varied by changing the SW density, while keeping its velocity and temperature the same. The corresponding transit light curves and absorption profiles are shown in Figure 13. Since the maxima of the simulated transit light curves are shifted relative to the ephemerid center of transit, the absorption profiles in the right panel of Figure 13 are shown for the time t=1 hr, i.e., slightly after the midtransit. One can see that the decrease/increase of the XUV flux and SW intensity directly influences the starting time of the early ingress and the duration of egress. The  maximum depth depends mostly on the density of SW, which influences the efficiency of ENAs' generation. In general, the variation of F XUV and ¢ M sw by a factor of two results in variation of the absorption in the blue wing of the Lyα line by a factor of 1.5, while the red wing remains unaffected.
Another parameter that affects the absorption is the helium abundance He/H. When it is very small, the PW outflow velocity is higher, resulting in a denser hydrogen cloud around the planet (Shaikhislamov et al. 2018b) and an extremely deep and broad transit light curve. Figure 14 shows the simulated transit light curves of GJ 436b in the blue wing of the Lyα line and the corresponding modeled Lyα line profiles at t=1 hr (i.e., close to the midtransit), obtained for different helium abundances with the model parameter sets No. 1 and No. 7-10. When He/H increases to about unity, the PW velocity significantly decreases because of the mass loading by heavier helium atoms, resulting in less hydrogen outflow and, correspondingly, less absorption in Lyα.
Despite the general good agreement between the modeling results and observations, some of the features, obtained in the simulations and shown above, do not fit the observations in all details. In particular, the ingress starts ∼1 hr later; the maximum of absorption is also shifted by ∼1hr; and the post-transit phase of the light curve looks like a simple decline, instead of the two-stage decreasing with a sharp fall, followed by a slow decay phase. Also, the depth of egress in general is much higher. While in the simulations with the parameter set No. 1 the absorption profile at the midtransit in the blue wing of the Lyα line is in a good general agreement with the observations, i.e., quantitative differences are not crucial (see in Figure 10), in the red wing there is a noticeable qualitative difference. The modeled absorption in the Doppler shift  Figure 11. Distribution of the line-of-sight absorption, averaged over the blue wing of the Lyα line, as seen by a remote observer at the midtransit. The stellar disk is shown with a black-line circumference, and GJ 436b with a black circle. Note that this picture uses the observer reference frame, in which the planet is shifted to Z O =−9 relative to the observer-star line (Z O =0, Y=0). velocity interval [30, 110] km s −1 at the midtransit is larger than at the post-transit, while in the observations the situation is the opposite.
The specially performed modeling study reveals that it is difficult to fit all the features of the observed Lyα absorption just by varying of the available modeling parameters. The start time of the ingress is easily shifted by a factor of two variation of the SW pressure (e.g., by manipulation with the SW density and, to a certain extent, with the velocity). It can be also adjusted by modifying the planetary mass-loss rate via the variation of the stellar XUV flux (see, e.g., Figure 13). However, fitting of the early ingress observations leads to an extended egress, resulting in an even stronger disagreement with the observations in the post-transit part of the light curve. On the other hand, the egress fall time can be adjusted by increasing the SW velocity and pressure. In this case, the ingress starts too late, and the absorption in the red wing of the Lyα line becomes too small. The absorption in the red wing, in its turn, is sensitive to the heating of planetary atoms due to the interaction with SW, and it can be increased by decreasing the SW velocity while increasing its temperature. The post-transit absorption in the red wing of the Lyα line at the level of about 20% is achieved in the simulations under the conditions of a rather small SW velocity, V swp <100 km s −1 , at the planet orbit.
In order to fit the transit light curves and the absorption profiles more precisely, checking of other physical parameters and possible processes was done. Specifically under consideration were (1) the heating efficiency by the XUV radiation, which influences the intensity of the PW outflow, while keeping the photoionization the same; (2) the cooling rate due to H3+, which affects the thermosphere structure; (3) the base temperature and density, which in a small degree control the mass-loss rate; and (4) the planet rotation rate around its axis, which to a certain degree affects the zonal flows in the thermosphere. However, the resulting manifestations of the considered parameter variations were found to be either small or not suitable for the intended fitting.
As was already pointed out in Section 2, the available observations can be interpreted as two different scenarios, realized under different conditions. Below we consider them both in more detail.  The first scenario with a sharp ingress, starting at −2 hr, and deep midtransit ∼70% phase, lasting for about 2hr, is followed by a sharp egress. In the course of the performed simulation runs, it has been found that this kind of transit light curve can be realized under the conditions of a very fast SW with an asymptotic terminal velocity reaching the values of = ¥ V sw, ¼ 700 2000 km s −1 . At such velocities, the spiraling of material flow, due to Coriolis force, in the planet-fixed frame of reference is less significant, resulting in a more pronounced sweeping of the PW tail. This can be seen in Figure 15, obtained with the modeling parameter set No. 12. In this case, the inclination of the planetary tail to the orbit trajectory is about 32°, whereas under the conditions of parameter set No. 1 (for comparison), it is only 16°.
As demonstrated in Figure 16, obtained in the simulations with the modeling parameter sets No. 11 and 12 (i.e., fast SW and high F XUV ), a lesser inclination of the planetary tail relative to the line of sight results in a significantly sharper post-transit decline of the light curve. To compensate for the higher pressure of the fast SW (proportional to the SW temperature and squared velocity) and to reproduce the effect of early ingress, one has to increase the pressure of PW as well. This is why the higher values of XUV flux are taken in the applied simulation parameter sets No. 11 and 12. Note that the Lyα line profile simulated with the parameter set No. 11 at the midtransit nicely fits the observations in the blue wing and remains also in good agreement with the integral depth in the red wing.
The right panel of Figure 15 shows the features that do not affect the observations but demonstrate the capacity of the multifluid model to reproduce the fine details of the complex PW and SW interaction. Specifically, the neutral particles, such as helium and hydrogen atoms, are coupled to ions in the PW flow but decouple beyond the ionopause, in the shocked region, where the density of the stellar protons is not high enough for the collisional mixing. As demonstrated in the right panel of Figure 15 for helium, as an example, in these regions the temperature of planetary neutrals is distinctly different from the temperature of SW. Closer to the star, the proton density becomes again sufficiently high for elastic collisions to stop helium atoms falling on the star. For a supersonic helium flow this results in the generation of a shock, positioned at the distance of ∼60R p from the star. At this distance, the density of atomic helium is orders of magnitude lower than that of the SW protons and the shock does not influence the parameters of SW. Correspondingly, the shock generated in the atomic hydrogen fluid is positioned farther from the star because hydrogen atoms and protons intermix with each other via charge exchange, which has a larger cross section than the elastic collisions.
The second scenario, with extended and shallow ingress and egress phases and a relatively small maximal depth of ∼25%, is suggested for the transits, revealed in visits 6 and 7 (Lavie et al. 2017). In our simulations such a scenario was realized under the conditions of slow SW with a low pressure. The slow SW velocity is needed to increase the absorption in the red wing of the Lyα line, while the low pressure enables expanding of the PW stream ahead of the planet, resulting in an earlier ingress. To compensate the overall decrease of the absorption due to the low density of SW and the related less amount of the generated ENAs, a slightly increased XUV flux is necessary. The transit light curves and Lyα line absorption profile at the post-transit, obtained with the modeling parameter set No. 13, are shown in Figure 17. One can see that the transit starts sufficiently early, which is in agreement with the data from visits 6 and 7, and lasts for a long time of ∼20hr at approximately the same level of ∼25%. The simulated Lyα absorption profile in the red wing agrees with the measured profile in the range [40, 65] km s −1 , though it differs at larger velocities.

Discussion and Conclusions
Let's estimate how high is the absorption by the expanded upper atmosphere of GJ 436b inside the Roche lobe, due to the natural line broadening mechanism. At the velocity of 60 km s −1 , the corresponding cross section of the Lyα line is ∼4× 10 −19 cm 2 . Taking the typical density n H ∼10 7 cm −3 (see, e.g., in Figure 9), the optical depth can be estimated as follows: τ∼n HI σ nat R Lag1 ∼0.06, where R Lag1 =5.8R p is the distance to the first Lagrange point. The integral absorption should be further reduced by a factor of (R Lag1 /R star ) 2 . Thus, as expected, the natural line broadening is too small for GJ 436b to play any significant role. The thermal broadening for the resonant atoms has a much larger cross section, · · -K T 6 10 10 14 4 cm 2 . As follows from the simulations in Figure 12, the bulk velocity V bulk in the shocked region varies in the range [−10, −150] km s −1 , whereas the thermal velocity of hot hydrogen atoms with the Maxwellian distribution at the temperature T∼4×10 5 K has a comparable value V thermal ∼V bulk . The integrated column density of ENAs and planetary atoms is about NL∼6×10 2 × 50×R p ∼10 14 cm −2 . Thus, the optical depth for the thermal line broadening absorption is about unity, NL res ∼1. Since the distributions of physical parameters, shown in Figure 12 for the particularly chosen line of sight, almost do not change over the significant part of the integration region, as can be seen in Figure 11, the absorption due to the thermal broadening mechanism should be significant in the blue wing of the Lyα line, in agreement with the observations. A typical product of the PW and SW interaction is the formation of a bow shock, or a thermally compressed region ahead of the planet, whereas the escaping PW material trails the planet as a tail, due to momentum conservation. The pressure exerted by the SW pushes the tail away from the star. The propagation of PW material above and below the ecliptic plane is restricted by the stellar gravity. One of the distinct features of the observed transit light curves is an early ingress, which starts at about t=−2hr. With respect to the planet, in view of its orbital motion, this corresponds to the position of about 30R p ahead. The supersonic interaction of PW with SW naturally creates such a sharp boundary in the form of ionopause and compression layer/bow shock. The position (so-called standoff distance) of the ionopause is determined by the pressure balance condition between the PW and SW. To estimate it, one can assume that the largest part of the SW total pressure is contributed by the dynamic pressure component related to the orbital motion of GJ 436b, because the background SW at the orbital distance of ∼0.028au is still less than the Keplerian speed and is approximately perpendicular to the direction of the orbital motion. Therefore, expressing the total pressure of PW via the planetary massloss rate as case, the duration of the absorption at the egress phase would be determined by the width of the PW stream. However, the performed simulations show that a very high SW velocity of ∼1000 km s −1 is needed for such a scenario to take place, which in the case of the Sun is achieved only in the coronal mass ejections. It has to be noted that reducing the egress duration by means of increasing the SW density, instead of velocity, results in a simultaneous cut of the early ingress part, which would contradict the observations. Another observational feature, detected in visits 1 and 6, that is not reproduced in simulations is the absorption in the red wing of the Lyα line at relatively high velocities [60, 110] km s −1 . However, as was pointed out in Lavie et al. (2017), this challenging-forinterpretation feature may be due to stellar variability.
Finally, the following physical picture, consistent also with the basic general estimates, follows from the undertaken numerical simulations. The supersonic interaction of SW with PW forms a bow shock ahead of GJ 436b. In the transition region between the ionopause and the bow shock the hydrogen atoms, carried in the PW stream, interact with the SW protons and generate ENAs. The absorption in Lyα observed for GJ 436b is produced mostly by these ENAs, due to the resonant, or thermal, line broadening mechanism, augmented by the planetary hydrogen atoms, energized in the hot transition layer. The velocity of ENAs, like that of protons, is directed approximately away from the star. This results in a strong asymmetry of the absorption between the blue and the red wings of the Lyα line. It has to be noted that no significant absorption comes from the regions beyond the compressed layer, because the hydrogen atoms are either too slow or too rare there. Since the rate of ENA production is proportional to the density of the planetary hydrogen atoms and stellar protons, the depth of the transit light curve directly depends on the density of SW and intensity of the PW outflow. The PW intensity, in its turn, is proportional to the XUV flux. Because of the bow shock formation, the grazing transit due to orbital inclination is actually favorable to the absorption in Lyα, since the compressed region in this case transits across the center of the stellar disk. The timing of early ingress is determined by the position of the bow shock ahead of the planet and depends mostly on the SW density and PW intensity. At the same time, the shape of the egress part of the transit light curve is formed by the SW velocity.
The modeling, reported in this paper, quantitatively reproduces such observed features as the Lyα absorption line profiles at the midtransit and post-transit phases (right panels of Figures 10 and 13, blue line), as well as the timing of early ingress and the maximum depth of the transit light curve in the blue wing (Figures 14, 16, and 17). There is also an agreement between the simulation results and observations regarding the strong asymmetry of the absorption in the blue and red wings of the Lyα line and in the range of Doppler-shifted velocities (<150 km s −1 ) at which the absorption in the blue wing takes place. However, the model could not reproduce successfully the egress part of the transit light curve. The suggested conclusion in that respect is twofold. It is possible that the available data are affected by strong variations of the SW parameters between the visits. At the same time, the effect of other physical factors, not yet included in the model, that may influence the interaction between the SW and PW around GJ 436b cannot be excluded as well. In particular, the magnetic fields of stellar and planetary origin can be important. As demonstrated by the existing generic MHD simulations (e.g., Bisikalo et al. 2013;Matsakos et al. 2015;Erkaev et al. 2017), the presence of magnetic field extends the complexity of physics of the planetary environment and directly influences the interaction between PW and SW. It therefore brings an additional parameter for fitting of the modeling results to observations.
Altogether, the first application of the 3D hydrodynamic self-consistent multifluid model for the interpretation of transit observations of GJ 436b, presented here, has shown that the concept of expanding and escaping upper atmospheres of hot exoplanets, developed earlier on the basis of 1D aeronomy models, gives adequate values for the density, velocity, and ionization degree of the outflowing PW material and explains reasonably well the Lyα absorption measurements for GJ 436b. The performed more realistic 3D simulations confirm the conclusion of previous 1D and 2D models that the generation of ENAs during the interaction between the PW and SW plasmas is the main process responsible for the observed absorption in the Lyα line, whereas the radiation pressure effects play a secondary and less important role. It is also demonstrated that the Monte Carlo modeling, used for the interpretation of observations in a number of cases before, has to be refined to include the parameters and geometry of the bow shock and/or ionopause, because the ENAs are generated in sufficient amount only in the compressed transition layer between them.
The remaining discrepancy between the simulated and observed Lyα transit light curves and the line profiles indicates, on one hand, that the modeling, being basically correct and promising, admits further refining by taking into account of additional factors. On the other hand, the model-based fitting of observations constrains the parameters of SW, which, at the same time, may vary from one observational campaign to another. This work was supported by grant No. 18-12-00080 of the Russian Science Foundation. M.L.K. and H.L. acknowledge the support from projects I2939-N27, S11606-N16, S11607-N16, and P25256-N27 of the Austrian Science Fund (FWF). Parallel computing simulations, key for this study, have been performed at the Computation Center of Novosibirsk State University, the SB RAS Siberian Supercomputer Center, the Joint Supercomputer Center of RAS, and the Supercomputing Center of the Lomonosov Moscow State University.

Appendix A Analytic Estimations of the Role of Radiation Pressure
Close to the planet the escaping PW material is sufficiently dense to absorb efficiently the Lyα flux. If the self-shielding is strong, then it acts as a surface pressure in the same manner as the SW ram pressure. The comparative role of both effects can be characterized by their ratio: where L Lyα is the total luminosity in the Lyα line and ¢ M sw is the total mass-loss rate of the star due to SW. Assuming typical values of F Lyα =1 erg cm −2 s −1 , at 1au, · ¢M 2.5 10 sw 11 gs −1 , and V sw =10 7 cm s −1 , one obtains that the radiation pressure is about one order of magnitude less than the SW ram pressure, p Lyα /p SW ∼0.1:ã P P 0.1.

Ly sw
Let us consider now the acceleration of neutrals by the radiation pressure only, i.e., disregarding the SW and gasdynamic interactions. The motion of a single atom in the stellar gravity field under the radiation pressure action can be easily solved, assuming the orbital momentum conservation. The equation for the radial velocity in the equatorial plane is The azimuthal velocity can be expressed as V j R=const., = j V GR R 2 orb 2 . Here we took into account that in the planet's frame of reference, from which the atom is launched, the centrifugal force is equal to the gravity force. The factor β in Equation (10) is the relation between the radiation pressure and gravity forces, and g(V ) is the stellar Lyα profile for which we take, for simplicity, but without loss of generality, the Gaussian profile, . For GJ 436, which radiates in the Lyα line at the reference distance of 1au F Lyα about 1 erg cm −2 s −1 , the radiation force is a fraction of β=0.7 of the stellar gravity, while the line width V Lyα is about 80 km s −1 . In the dimensionless form for GJ 436b, with v=V R /V Lyα , r=R/R orb , where R orb =4.35×10 11 cm and V orb =120 km s −1 , Equation (10) can be rewritten as follows: Maximum terminal velocity, due to the finite line width, can be found by direct integration of equation a r 2 with the result: Ly orb 2 Ly 2 -95 km s 1 . However, the actual acceleration of an atom by the radiation pressure force is restricted by the decrease of centrifugal force with distance, which compensates for the gravity pull. An approximate analytic solution, which is sufficiently precise up to the maximum velocity, can be found by the series expansion of the O(1−1/r) term in Equation (11): ln 1 1 1 0.5 1 1 0.25 1 1 . 12 2 2 3 Figure 18 shows the results of numerical integration of Equation (11) for the cases of β=0.7, 05, and 1. One can see that the maximum velocity does not exceed 60 km s −1 . It is reached rather far from the planet, at ∼0.5R orb ∼ 7R star . Another factor that strongly decreases the impact of radiation pressure is the gasdynamic coupling between the neutrals and protons (Shaikhislamov et al. 2016). Due to charge exchange with a cross section σ∼3×10 −15 cm −2 , the hydrogen atoms effectively exchange momentum with protons over the distance of ∼R p at densities as low as n=10 5 cm −3 . Thus, when the density of protons exceeds that of atoms, the propulsion effect of the radiation pressure force on the plasma as a whole correspondingly decreases. This is especially true at small velocities. This effect can be taken into account by decreasing Figure 18. Radial velocity of a hydrogen atom as a function of distance, calculated by integration of Equation (11), for different values of the radiation pressure force. The dotted line shows the analytic solution (12) for β=0.7. The dashed line shows the solution of an extended version of Equation (11) with the inclusion of ionization of hydrogen for the case of β=0.7. the radiation pressure force (i.e., the parameter β) by a factor of n a /(n a +n p )≈exp (−t/τ ion ), which depends on the ionization time of hydrogen atoms. The time of particle's flight is found as t=(R orb /V Lya ) ò dr/v. For the particular spectrum of GJ 436b (France et al. 2016) with τ ion =29 hr, the dimensionless parameter R orb /(V Lya τ ion )≈0.5 is of the order of unity. The numerical solution of the extended version of Equation (11) with the inclusion of ionization is shown in Figure 18. One can see that photoionization reduces the velocity range of the accelerated by the radiation pressure hydrogen atoms down to values that appeared within the interstellar medium absorption and geocoronal contamination window. Therefore, these atoms will not contribute to the observed Lyα line absorption profile.
Altogether, the quantitative arguments presented above show that the radiation pressure force of GJ 436 is relatively insignificant, as compared to other forces, and alone it is not able to accelerate hydrogen atoms up to velocities necessary to produce the measured absorption in the Lyα line. This explains why in our simulations we see that the effect of the radiation pressure is very small and the main observational feature-the absorption of the Lyα line-is practically independent of it. To verify the performed estimates, we did also a simulation run with a very weak SW ( ¢M 10 sw 10 gs −1 ) and 20 times increased Lyα flux (∼ 20 erg cm −2 s −1 at 1au). Only in this extreme case did we see the increased manifestation of the radiation pressure effect, resulting in a strong absorption with a depth comparable to that measured in observations.

Appendix B Technical Details Regarding the Numerical Code of Model
Our code employs an explicit numerical scheme with an upwind donor cell method for flux calculations. To achieve a second-order spatial accuracy for the differentials, two grids are used, shifted relative to each other on a half step along each dimension. One grid is reserved for the densities, temperatures, and gravity potential, whereas another grid is reserved for velocities. To achieve the second-order accuracy in time, the code at each time step calculates at first n and T values using velocity field V from the previous time step and then recalculates V, using the new n and T values. The applied numerical scheme fully conserves fluxes and the total mass. Also, the Bernoulli constant along the characteristics is conserved. For the energy, a simple Equation (3) in a nonconservative form is used. The conservation of energy is checked by global integration, and the result is used to evaluate the accuracy of the simulations. Usually, the energy is balanced within the range of±25%. The code does not use any particular method to capture the shocks. For the problems under consideration, an accurate revealing of the position of shock and high resolution of its front are not crucial.
The spatial spherical grid uses a uniform step for the azimuthal angle. In the present modeling of GJ 436b the azimuthal step is Δj=0.065, which corresponds to 96 equidistant points over the whole circumference. In the radial direction the grid step varies linearly with the radius as follows: Δr=Δr min +(Δr max -Δr min )×(r-R p )/(R max -R p ), where Δr min is taken at the planet surface as small as R p /400, and for Δr max , a value equal to Δj×R max is applied. Therefore, in the shock region at ∼20R p the resolution of grid is about R p . For the polar angle, the grid step is defined as a quadratic function, Δθ=Δθ min ×i+α×i 2 , with the minimal step Δθ min at the equatorial plane. The parameters i and α usually are selected so that Δθ=Δθ min =Δj at the equator, whereas at the polar axis Δθ=2Δj. Note that the applied radial spacing in the spherical coordinate system allows keeping the same resolution in all three dimensions, if the azimuthal and latitudinal steps are chosen so that Δj≈Δθ≈Δr/r. The influence of spatial resolution was checked by doubling the number of grid steps in each dimension. By this, the difference in results of the test runs appeared sufficiently small.
For the calculation of column density along stellar rays in the star-centered coordinate system, we use the numerically consuming but straightforward procedure of integration from the star to each elementary cell in the planet-centered spherical frame, where the modeling was performed. The density values (but not the column densities) interpolated from the data in the nearby pixels are used along the integration path. The calculated column density is used to determine attenuation of XUV flux in each spectral bin (0.1 nm). For the Lya, the opacity is calculated in 30 km s −1 bins (in Doppler-shifted velocity units). The scattering of primary photons is not taken into account, whereas the re-emitted photons are not considered because they do not produce a net radiation force anymore. To save numerical time, the radiation transfer is calculated usually each fourth step of fluid dynamics and chemistry. We also apply an optically thin approximation for the photons generated by proton recombination to the ground state, so the standard recombination coefficients are used (Verner & Ferland 1996).
The photochemistry reactions (the same list as those used in Khodachenko et al. 2015) are calculated by direct conversion of matrix dn i =R ji (t, r)×n j ×n i at each time step and at each pixel. This is not efficient numerically, but it eliminates the convergence problem caused by rather different reaction rates R ji .
The boundary of the simulation domain connected to the star is treated in the same manner as the planet-related boundary, i.e., by imposing and keeping particular values. To launch the stellar wind in a more accurate way, the fixed values are prescribed at 1.5R star according to a predefined analytic solution. In fact, the spatial resolution of the grid around the star is rather coarse, just several points per stellar diameter scale, but it is sufficient for the SW simulation within the goals of the present study.
In the course of the present study the code was executed over 500 characteristic time units τ=R p /V 0 ∼0.75 hr, which correspond to approximately six orbital periods of the planet. As the main criterion for achieving a steady state, the convergence of the planetary mass-loss rate to a quasi-constant value has been applied. In fact, already after about three orbits the mass-loss rate within an accuracy of 5% reached an asymptotic quasi-stationary value. At the same time, no disruption-type variabilities of the planetary material flow or any instabilities of the shock interface during the interaction between PW and SW have been observed in the simulations.