Transit Ly-$\alpha$ signatures of terrestrial planets in the habitable zones of M dwarfs

TL;DR

This study models the Lyman-alpha transit signatures of Earth-sized planets in the habitable zones of M dwarfs, assessing their detectability with current telescopes and validating the models against known exoplanet observations.

Contribution

It introduces a detailed modeling approach for Lyman-alpha signatures of terrestrial exoplanets and evaluates their observability with existing instruments, expanding the potential for detecting Earth-like planets.

Findings

Hydrogen-dominated atmospheres are marginally detectable around GJ 436.

Nitrogen atmospheres are not detectable in Lyman-alpha transit observations.

Model successfully reproduces observed signatures of GJ 436b.

Abstract

We modeled the transit signatures in the Lya line of a putative Earth-sized planet orbiting in the HZ of the M dwarf GJ436. We estimated the transit depth in the Lya line for an exo-Earth with three types of atmospheres: a hydrogen-dominated atmosphere, a nitrogen-dominated atmosphere, and a nitrogen-dominated atmosphere with an amount of hydrogen equal to that of the Earth. We calculated the in-transit absorption they would produce in the Lya line. We applied it to the out-of-transit Lya observations of GJ 436 obtained by the HST and compared the calculated in-transit absorption with observational uncertainties to determine if it would be detectable. To validate the model, we also used our method to simulate the deep absorption signature observed during the transit of GJ 436b and showed that our model is capable of reproducing the observations. We used a DSMC code to model the planetary exospheres. The code includes several species and traces neutral particles and ions. At the lower boundary of the DSMC model we assumed an atmosphere density, temperature, and velocity obtained with a hydrodynamic model for the lower atmosphere. We showed that for a small rocky Earth-like planet orbiting in the HZ of GJ436 only the hydrogen-dominated atmosphere is marginally detectable with the STIS/HST. Neither a pure nitrogen atmosphere nor a nitrogen-dominated atmosphere with an Earth-like hydrogen concentration in the upper atmosphere are detectable. We also showed that the Lya observations of GJ436b can be reproduced reasonably well assuming a hydrogen-dominated atmosphere, both in the blue and red wings of the Lya line, which indicates that warm Neptune-like planets are a suitable target for Lya observations. Terrestrial planets can be observed in the Lya line if they orbit very nearby stars, or if several observational visits are available.

Tables

Table 1. Table 1: Stellar, planetary, and stellar wind parameters used in our test simulation for GJ 436b. The stellar and planetary data are from www.exoplanet.eu and Bourrier et al. ( 2016 ) , stellar wind parameters are adopted from Vidotto & Bourrier ( 2017 ) , except of a slightly faster stellar wind speed of 110 km/s instead of 85 km/s. Simulation inner boundary for the DSMC code has been adopted close to the Roche lobe @ ​ 5.83 @ 5.83 @5.83  R pl . The parameters at the inner boundary have been calculated by the hydrodynamic code. The shift of the planet relative to the stellar center in the direction of the z 𝑧 z axis has been calculated according to the inclination. Assumed collision cross sections can be found in Table  2 .

Parameter	Value
Star mass $M_{\mathrm{st}}$ [kg]	$8.98\times {10}^{29}$
Star radius $R_{\mathrm{st}}$ [m]	$3.22\times {10}^8$
Planet’s mass $M_{\mathrm{pl}}$ [kg]	$1.32\times {10}^{25}$
Planet’s radius $R_{\mathrm{pl}}$ [m]	$2.65\times {10}^7$
Semi-major axis $d_{\mathrm{HZ}}$ [AU]	0.0288
Inclination [degree]	85.8
SW number density $n_{\mathrm{sw}}$ [m-3]	$2.0\times {10}^9$
SW velocity $v_{\mathrm{sw}}$ [m$\cdot$s-1]	$110\times {10}^3$
SW temperature $T_{\mathrm{sw}}$ [K]	$0.41\times {10}^6$
Inner boundary radius [m]	$1.54\times {10}^8$
Inner boundary density [m-3]	$2.5\times {10}^{12}$
Inner boundary temperature [K]	2693
Electron impact ionizaion rate (hydrogen) ${\tau }_{\mathrm{ei}}$ [s-1]	$5.6\times {10}^{-5}$
Photoionizaion rate (hydrogen) ${\tau }_{\mathrm{ph}}$ [s-1]	$8.8\times {10}^{-6}$
Bulk atmospheric outflow velocity [m s-1]	$7.9\times {10}^3$
Magnetospheric obstacle substellar point $R_s$ [m]	$1.6\times {10}^8$
Magnetospheric obstacle width $R_t$ [m]	$1.6\times {10}^8$

Table 2. Table 2: Stellar, planetary, and stellar wind parameters used in all simulations. Stellar wind and habitable zone parameters of GJ436 have been adopted from Kislyakova et al. ( 2013 ) , electron impact ionization rates were calculated using the NIST databases. Charge exchange cross sections are from Lindsay & Stebbings ( 2005 ) for H + -H and Rinaldi et al. ( 2011 ) for N + -H. Elastic collision crossections were calculated as described in section  2.2 .

Parameter	Value
Planet’s mass $M_{pl}$ [kg]	$5.97\times {10}^{24}$
Planet’s radius $R_{pl}$ [m]	$6.371\times {10}^6$
HZ distance $d_{\mathrm{HZ}}$ [AU]	0.12 – 0.36
Orbital distance $d_{\mathrm{orb}}$, [AU]	0.24
Electron impact ionizaion rate (N) ${\tau }_{\mathrm{ei}}$ [s-1]	$6.3\times {10}^{-6}$
Electron impact ionizaion rate (H) ${\tau }_{\mathrm{ei}}$ [s-1]	$1.0\times {10}^{-5}$
SW number density $n_{\mathrm{sw}}$ [m-3]	$2.5\times {10}^8$
SW velocity $v_{\mathrm{sw}}$ [m$\cdot$s-1]	$330\times {10}^3$
SW temperature $T_{\mathrm{sw}}$, [K]	$2\times {10}^6$
Charge exchange crossection H+-H [m2]	$2.0\times {10}^{-19}$
Charge exchange crossection H+-N [m2]	$1.0\times {10}^{-19}$
H-H elastic collision crossection [m2]	$3.5\times {10}^{-20}$
H-N elastic collision crossection [m2]	$3.7\times {10}^{-20}$
H-N elastic collision crossection [m2]	$3.9\times {10}^{-20}$

Table 3. Table 3: Atmospheric parameters used in the simulations for nitrogen- and hydrogen-dominated atmospheres for 3, 7, 10, and 20 XUV. Adopted from Tian et al. ( 2008 ) (3 and 20 XUV) and Lammer et al. ( 2008 ) (7 and 10 XUV). The value of the hydrogen density at the modern terrestrial exobase has been adopted from Kislyakova et al. ( 2013 ) . The value of R s subscript 𝑅 s R_{\rm s} is chosen arbitrarily slightly above the exobase.

Parameter	3 XUV	7 XUV	10 XUV	20 XUV
Nitrogen-dominated atmosphere with Earth-like hydrogen content
$R_{\mathrm{ib}}$ [m]	$7.4\times {10}^6$	–	–	–
$n_{\mathrm{ib}}$ [m-3] (N)	$1.2\times {10}^{13}$
$n_{\mathrm{ib}}$ [m-3] (H)	$7.0\times {10}^{10}$
$R_s$ [m]	$7.4\times {10}^6$
$T_{\mathrm{ib}}$ [K]	$2.2\times {10}^3$
Nitrogen-dominated atmosphere
$R_{\mathrm{ib}}$ [m]	$7.4\times {10}^6$	$2.15\times {10}^7$	$2.7\times {10}^7$	$8.6\times {10}^7$
$n_{\mathrm{ib}}$ [m-3] (N)	$1.2\times {10}^{13}$	$2.2\times {10}^{12}$	$4.6\times {10}^{11}$	$4.4\times {10}^{11}$
$R_s$ [m]	$7.4\times {10}^6$	$2.23\times {10}^7$	$2.8\times {10}^7$	$8.64\times {10}^7$
$T_{\mathrm{ib}}$ [K]	$2.2\times {10}^3$	$7.5\times {10}^3$	$5.0\times {10}^3$	$2.53\times {10}^3$
Hydrogen-dominated atmosphere
$R_{\mathrm{ib}}$ [m]	$4.25\times {10}^7$	$5.65\times {10}^7$	$6.3\times {10}^7$	$8.2\times {10}^7$
$n_{\mathrm{ib}}$ [m-3] (H)	$2.5\times {10}^{11}$	$1.74\times {10}^{11}$	$1.54\times {10}^{11}$	$1.07\times {10}^{11}$
$R_s$ [m]	$4.31\times {10}^7$	$5.73\times {10}^7$	$6.37\times {10}^7$	$8.28\times {10}^7$
$T_{\mathrm{ib}}$ [K]	$1.14\times {10}^3$	$2.4\times {10}^3$	$3.04\times {10}^3$	$4.69\times {10}^3$
Species photoionization rates
H photoion. rate [s-1]	$1.43\times {10}^{-7}$	$3.17\times {10}^{-7}$	$4.82\times {10}^{-7}$	$1.03\times {10}^{-6}$
N photoion. rate [s-1]	$1.15\times {10}^{-6}$	$2.14\times {10}^{-6}$	$2.97\times {10}^{-6}$	$5.43\times {10}^{-6}$

Equations

$$n_{\mathrm{exo}}=n_{\mathrm{ib}}=\frac{1}{\sigma H},$$

$${\sigma }_{a-b}=\pi (R_a+R_b{)}^2,$$

$$\begin{array}{r}{\sigma }_{N-N}=3.9\times 10^{-20},\\ {\sigma }_{H-H}=3.5\times 10^{-20},\\ {\sigma }_{H-N}=3.7\times 10^{-20}.\end{array}$$

$$x=R_{\mathrm{s}}\left(1-\frac{y^2+z^2}{R_{\mathrm{t}}^2}\right).$$

$$T=I/I_0=e^{-\tau (f)}=e^{-\sigma (f)Q}.$$

$$\bar{T}(f)=\frac{1}{N_c}\sum _{i=1}^{N_c}{T}_i(f),$$

$$I_m(f)=I_0(f)\bar{T}(f).$$

$$\sigma (f)={\int }_{-\infty }^{+\infty }\check{u}(v_x)\text{leavevmode}{\sigma }_N[(1+v_x/c)f]dv\text{leavevmode}[m^2],$$

$${\sigma }_N(f)=\frac{3{\lambda }_{\alpha }^2{A}_{21}^2}{8\pi }\frac{(f/f_{\alpha }{)}^4}{4{\pi }^2(f-f_{\alpha }{)}^2+(A_{21}^2/4)(f/f_{\alpha }{)}^6}\text{leavevmode}[m^2].$$

$$\tau (f)={\int }_{-\infty }^{+\infty }u(v_x)\text{leavevmode}{\sigma }_N[(1+v_x/c)f]dv.$$

$${\sigma }_N(f)=f_{12}\frac{\pi e^2}{m_{e}c}\frac{\Delta f_L/2\pi }{(f-f_{\alpha }{)}^2+(\Delta f_L/2{)}^2},$$

Full text

*Methods. *We used a direct simulation Monte Carlo (DSMC) code to model the planetary exospheres. The code includes several species and traces neutral particles and ions. It includes several ionization mechanisms, such as charge exchange with the stellar wind, photo- and electron impact ionization, and allows to trace particles collisions. At the lower boundary of the DSMC model we assumed an atmosphere density, temperature, and velocity obtained with a hydrodynamic model for the lower atmosphere.

*Results. *We showed that for a small rocky Earth-like planet orbiting in the HZ of GJ 436 only the hydrogen-dominated atmosphere is marginally detectable with the Space Telescope Imaging Spectrograph (STIS) on board the Hubble Space Telescope (HST). Neither a pure nitrogen atmosphere nor a nitrogen-dominated atmosphere with an Earth-like hydrogen concentration in the upper atmosphere are detectable. We also showed that the Ly$\alpha$ observations of GJ 436b can be reproduced reasonably well assuming a hydrogen-dominated atmosphere, both in the blue and red wings of the Ly$\alpha$ line, which indicates that warm Neptune-like planets are a suitable target for Ly$\alpha$ observations. Terrestrial planets, on the other hand, can be observed in the Ly$\alpha$ line if they orbit very nearby stars, or if several observational visits are available.

1 Introduction

Low-mass M dwarfs are the most numerous stars in the Universe and are often hosts to small rocky planets of a size and mass similar to the Earth. These planets are interesting targets for modern planetology, because conditions in the systems of M dwarfs differ significantly from those in the planetary systems around G stars, thus raising questions about the evolution of such planets (e.g., Scalo et al. 2007).

Lyman-alpha (Ly$\alpha$) transit observations are a powerful tool of characterizing the planetary atmospheres. For hot Jupiters, in a number of cases these observations can be used as an indicator of escape rates from hydrogen-dominated atmospheres (e.g., Vidal-Madjar et al. 2003; Lecavelier des Etangs et al. 2012; Bourrier & Lecavelier des Etangs 2013; Khodachenko et al. 2017), atmosphere properties (Ben-Jaffel 2007) and planetary magnetic moments and stellar winds (Kislyakova et al. 2014a; Vidotto & Bourrier 2017). Besides hot Jupiters, a variety of smaller exoplanets has also been observed in Ly$\alpha$, including planets for which a deep absorption signature indicating a presence of a dense and extensive hydrogen envelope has been detected, such as GJ 436 (Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017) and GJ 4370b (Bourrier et al. 2018), and planets for which no excess absorption in comparison to the visible light has been found (Ehrenreich et al. 2012; Bourrier et al. 2017a). Different character of observed in-transit signals of exoplanets indicates that Ly$\alpha$ observations can be used as a powerful tool not only to detect a presence of a hydrogen envelope around an exoplanet, but also its shape and density. Ly$\alpha$ observations can also give us hints on other parameters such as planetary magnetic moment and stellar wind density and velocity, and the type of interaction between planetary atmosphere and stellar wind (Khodachenko et al. 2015, 2017; Shaikhislamov et al. 2016). However, one should keep in mind that interpreting the Hubble Space Telescope (HST) Ly$\alpha$ observations can be very challenging (for instance, see discussion in Ben-Jaffel 2007 and Vidal-Madjar et al. 2008, and in other papers cited above).

At present, Ly$\alpha$ observations can be obtained only by the Hubble Space Telescope, as this is the only space observatory equipped with an UV instrument. In the future, other space observatories such as LUVOIR (France et al. 2017) and WSO UV (Boyarchuk et al. 2016) may provide us with more precise observations.

In this paper, we estimated the transit depth in the Ly$\alpha$ line of a terrestrial planet orbiting in a habitable zone (HZ) of an M star. Our ultimate goal was to examine if planets with nitrogen-dominated atmospheres showed an excess absorption in the Ly$\alpha$ line due to hydrogen in their exospheres originating from the charge exchange between incoming stellar wind protons and exospheric neutrals. If these signatures would be detectable, the observations in the Ly$\alpha$ line could be used as a characterization tool not only for hydrogen-dominated planets, but also for exoplanets with secondary atmospheres. We have used Ly$\alpha$ observations of an M dwarf GJ 436 of a spectral type M2.5, which hosts a transiting hot Neptune GJ 436b (Ehrenreich et al. 2015; Lavie et al. 2017). This planet has been observed to produce an immense in-transit Ly$\alpha$ absorption, thus allowing for a sophisticated modeling of its hydrogen envelope and other parameters of the system (Bourrier et al. 2015, 2016; Vidotto & Bourrier 2017; Loyd et al. 2017). For this study, we have selected out-of-transit observations of GJ 436 published by Lavie et al. (2017), which we used as stellar background spectrum for our modeling. We modeled atomic coronae surrounding an Earth-mass planet orbiting in the habitable zone (HZ) of GJ 436 with a nitrogen and a hydrogen-dominated atmosphere and then calculated absorption in the Ly$\alpha$ line which would be produced by a transit of such an exospheric cloud. Then, we imposed this calculated absorption on the observed Ly$\alpha$ profile of GJ 436 and compared it with the observational uncertainties. This method allowed us to conclude if planets with different types of atmospheres could be observable in the Ly$\alpha$ line by the HST. We performed the modeling for several XUV (X-ray and extreme ultraviolet) enhancement factors: three times the modern solar flux at Earth (which is the nominal XUV flux in the middle of the HZ of GJ 436 according to observations by Ehrenreich et al. 2011), seven, ten, and twenty times. For the case of 3 XUV, we include an additional case of a nitrogen-dominated atmosphere with a hydrogen content equal to the one at the Earth’s exobase.

This paper is organized as follows: in Section 2, we describe our models and available codes we use for the simulations. To validate our method, in Section 3 we apply our model to GJ 436b and show that our model can reproduce a deep observed transit signature of this planet. In Section 4, we describe simulation parameters and inputs. Section 5 summarizes our results. Section 6 describes possible applications of the results to observation planning and compares our results with the literature, and, finally, Section 7 summarizes our conclusions.

2 Model description

2.1 The star. GJ 436

GJ 436 is an M dwarf of spectral class M2.5 with a mass of 0.452 M⊙ and a radius of 0.464 R⊙. We have selected GJ 436 as a target because for this star Ly$\alpha$ in-transit and out-of-transit observations with the Space Telescope Imaging Spectrograph (STIS) on board the Hubble Space Telescope (HST) are available. Due to its relatively close distance to the Earth of only 10.2 pc, this observation is of a very high quality and has a high signal-to-noise ratio. Transit observations of GJ 436b by the HST revealed an immense absorption in the stellar Ly$\alpha$ line, especially in its blue wing (Ehrenreich et al. 2015). Later visits with longer exposure times confirmed the presence of a deep in-transit, egress, and ingress absorption, with the mid-transit signature being the deepest (Lavie et al. 2017). This indicates a presence of a massive hydrogen cloud around this planet, which obscures almost half of the emission in the stellar Ly$\alpha$ line (the absorption during a transit is $\approx$56%). Numerical modeling confirms that GJ 436b is surrounded by a big cloud of atomic hydrogen with a complex structure (Bourrier et al. 2015, 2016; Loyd et al. 2017; Shaikhislamov et al. 2018).

In this article, we aimed at investigating if a smaller planet located further from GJ 436 than GJ 436b, namely, an Earth-like planet in the habitable zone of this star with the orbital distance of 0.24 au, could also produce a detectable Ly$\alpha$ signature. For this purpose, we used out-of-transit observations of GJ 436b and calculated the Ly$\alpha$ absorption caused by a transit of such a planet in the HZ. M dwarfs such as GJ 436 are very numerous in the solar neighborhood, and therefore this study can be used for predicting the amount of Ly$\alpha$ absorption small rocky planets orbiting them can produce.

2.2 Atmosphere parameters

Even if a planet is similar to Earth in terms of mass and size, it does not necessarily have an atmosphere similar to our planet, in other words, dominated by nitrogen. The amount of hydrogen the planet accumulates depends on its formation history. If a planet grows to a mass close to Earth’s mass before the protoplanetary nebula dissipates, it can accumulate a dense hydrogen envelope (Stökl et al. 2016). After the nebula dissipates, this envelope will partly escape due to enhanced thermal escape (Owen & Wu 2016; Fossati et al. 2017), however, Earth-like planets in the HZs can still keep a significant fraction of it. If the hydrogen envelope does not disappear during the first tens of Myr after the planet has been released from the nebula, it can still escape afterwards, both by thermal (e.g., Erkaev et al. 2013; Lammer et al. 2014; Luger et al. 2015; Johnstone et al. 2015; Kubyshkina et al. 2018) and non-thermal (e.g., Kislyakova et al. 2013; Airapetian et al. 2017) escape mechanisms. Of course, the exact evolutionary path of a planet strongly depends on the stellar activity evolution and on the initial amount of hydrogen the planet accumulates.

On the other hand, if a planet never accumulates enough hydrogen, or if the star is active enough early on, it is possible that a secondary atmosphere forms. In this article, we focus on nitrogen-dominated atmospheres. One should keep in mind that these atmospheres may not be stable around active M dwarfs, because a nitrogen-dominated atmosphere readily expands if exposed to strong XUV radiation (Tian et al. 2008), which leads to strong non-thermal losses and, consequently, to a rapid loss of the whole atmosphere (Lichtenegger et al. 2010; Airapetian et al. 2017). Due to their slow evolution, M stars can keep high levels of XUV radiation in their HZs for a longer time in comparison to solar-like stars (West et al. 2008), which can further enhance non-thermal losses. The atmospheres of planets orbiting active low mass M dwarfs may be replenished by extreme volcanic activity due to tidal or induction heating (Driscoll & Barnes 2015; Kislyakova et al. 2017, 2018). A volcanically active planet would eject gases composed mainly of SO2, CO2, H2O, and S2. The exact composition depends on the redox state and volatile abundances in planetary interiors, and also on the surface pressure (Gaillard & Scaillet 2014), but a volcanically produced nitrogen-dominated atmosphere seems unlikely. At Earth, both life and plate tectonics support the nitrogen-dominated atmosphere (Lammer et al. 2018).

It seems reasonable to assume that Earth-like planets in the HZs of M stars can have both types, hydrogen-dominated and secondary, atmospheres. Observations confirm the presence of planets with and without a hydrogen envelope, which can be concluded about based on their average density (e.g., Lopez et al. 2012; Ginzburg et al. 2016; Fossati et al. 2017).

Hydrogen atmospheres. To obtain the parameters at the exobase of a hydrogen-dominated atmosphere, we applied a hydrodynamic code, which is a 1D upper atmosphere radiation absorption and hydrodynamic escape model that takes into account ionization, dissociation and recombination (Erkaev et al. 2016, 2017). The code provides mass-loss rates and atmospheric structure under different radiation conditions, and has been successfully applied multiple times to exoplanets (e.g., Erkaev et al. 2013, 2015; Lammer et al. 2014, 2016).

The code calculates the absorption of the stellar XUV flux by the thermosphere and solves the hydrodynamic equations for mass, momentum and energy conservation, and the continuity equations for neutrals and ions (both atoms and molecules). The model also accounts for dissociation, ionization, recombination and Ly$\alpha$ cooling. The quasi-neutrality condition determines the electron density. The model does not self-consistently calculate the ratio of the net local heating rate to the rate of the stellar radiative energy absorption. For this reason, we assumed a heating efficiency of 15 per cent, which is a mean value for a hydrogen-dominated atmosphere according to kinetic studies solving the Boltzmann equation (Shematovich et al. 2014). We run the model for 3, 7, 10, and 20 XUV and took the values for the number density, temperature, and gas outflow velocity at the exobase as an inner boundary input for the DSMC model.

Nitrogen atmospheres. For nitrogen-dominated atmospheres, we adopted the parameters from the literature, namely, from the paper by Lammer et al. (2008) based on the results by Tian et al. (2008), who studied the response of Earth’s atmosphere to enhanced levels of XUV radiation using a hydrodynamic atmosphere model. Tian et al. (2008) and Lammer et al. (2008) have shown that at high XUV fluxes, an Earth-like nitrogen-dominated atmosphere with a present atmospheric level of CO2 (1 PAL) expands to high altitudes and is dominated by atomic nitrogen. Venus, on the other hand, despite being exposed to a higher XUV flux than Earth, has a much less expanded atmosphere and a cooler exobase temperature due to the atmospheric composition dominated by CO2 (96% CO2; e.g., Lammer et al. 2008). High sensitivity of nitrogen-dominated atmospheres to stellar XUV flux is confirmed by other authors, who applied a hydrodynamic code combined with a sophisticated atmospheric chemistry model (Johnstone et al. 2018). Similar to the hydrogen-dominated atmospheres, we adopted the atmospheric parameters for 3, 7, 10, and 20 XUV (Lammer et al. 2008). While such radiation levels are typical mostly for young G dwarfs, M dwarfs often exhibit higher short wavelength radiation levels in their HZs even at older ages (e.g., West et al. 2008; France et al. 2016; Youngblood et al. 2016). At present, GJ 436 is a relatively quiet M dwarf with an XUV flux in the HZ of only 3 XUV (see Section 2.4), however, this star could have produced a higher level of XUV radiation in the past.

We note that the nitrogen atmosphere can possibly escape within geologically short times even at 10 and 20 XUV (Lichtenegger et al. 2010). This escape may be even stronger for non-magnetized planets such as those we consider in this article, especially taking into account that the stellar wind in the HZs of M dwarfs may be much denser than the solar wind at 1 AU, due to the proximity of the HZs of M dwarfs to their hosts. However, we still included these cases in our article to check if short-living nitrogen atmospheres may be detectable in the Ly$\alpha$ line, even though nitrogen-dominated atmospheres may not be stable for a geologically long time. An increased amount of CO2, which is a sufficient coolant in the thermosphere, may protect the atmosphere from expansion and enhanced non-thermal escape (Lichtenegger et al. 2010; Johnstone et al. 2018), but modeling of different compositions of atmospheres is beyond the scope of the present article.

We adopted the inner boundary for our DSMC simulations at the exobase, at a radius $R_{\rm ib}$, with a temperature, $T_{\rm exo}$. The values of $R_{\rm ib}$ and $T_{\rm exo}$ were adopted from Lammer et al. (2008). Since the densities were directly not available in this article, we estimated them from the temperature profiles as

[TABLE]

where $H$ is the local scale height and $\sigma$ is the collision cross section. The local scale height is determined as $H=k_{B}T_{\rm exo}/mg_{\rm exo}$, where $k_{B}$ is the Boltzmann constant, $m$ is the mass of an atom of the main species in the upper atmosphere (in this case, nitrogen), and $g_{\rm exo}$ is the local gravitational acceleration at the exobase. We assumed that nitrogen atoms were the only component of the upper atmosphere, which yielded a slight overestimate of their number density as Tian et al. (2008) have also taken into account other atmospheric constituents such as oxygen and others. However, as we showed below, even for a slightly overestimated density a nitrogen atmosphere does not produce any detectable Ly$\alpha$ signature, therefore, our results are not influenced by this factor.

To estimate collision cross sections between neutral species, we used a simplified approach following Atkins & de Paula (2000)

[TABLE]

where $\sigma_{a-b}$ is the elastic collision cross section between the two species, and $R_{a}$ and $R_{b}$ are the atomic radii of corresponding species. For the species included in the simulations, we obtained (in m2)

[TABLE]

We did not take into account dependence of collisional cross section on particles velocity. Elastic collisions with and between ions were not included in our model, however, we included charge exchange collisions between stellar wind protons and atmospheric neutrals.

Nitrogen-dominated atmosphere with an Earth-like hydrogen amount. As an additional test, we calculated also Ly$\alpha$ absorption of a nitrogen-dominated atmosphere at 3 XUV assuming the presence of an additional hydrogen content in the planetary exosphere, equal to the modern terrestrial level with the density at the exobase of $7.0\times 10^{10}$ m*-3* (Kislyakova et al. 2013). This hydrogen density is in agreement with the observations of the Earth exosphere by the IBEX satellite (Fuselier et al. 2010).

2.3 Stellar wind

Stellar wind interacts in various ways with planetary atmospheres. It compresses the magnetospheres, in case of planets with intrinsic magnetic fields, and compresses and penetrates planetary ionospheres of non-magnetized planets. Stellar wind is a powerful ionization source due to electron impact ionization and charge exchange between stellar wind particles and atmospheric neutrals. As a result of a charge exchange reaction, an Energetic Neutral Atom (ENA) and a slow ion of planetary origin are produced, because the electron is transferred from the planetary neutral particle to the stellar wind proton, which keeps its energy and velocity. Simulations show that stellar wind can penetrate deep into the atmospheres of non-magnetized planets (Cohen et al. 2015), and lead to intense atmospheric erosion and shape the hydrogen cloud surrounding the planet (Kislyakova et al. 2014b; Vidotto & Bourrier 2017; Shaikhislamov et al. 2016). A similar interaction pattern is also observed at Mars and Venus (Galli et al. 2008a, b; Lundin 2011). If charge exchange occurs between heavy ions of the stellar wind and neutral atoms from an atmosphere, a powerful secondary X-ray emission is generated (Kislyakova et al. 2015). However, in this article we disregarded the effects of heavy ions and focused on charge exchange between stellar wind protons and atmospheric neutrals, which generates ENAs. These atoms keep the velocity of a high energetic stellar wind proton and contribute to the absorption in the wings of the Ly$\alpha$ line. Due to the fact that the electron is transferred to ENAs from an atmospheric neutral atom, this process also contributes to overall ionization in the planetary exosphere. Since charge exchange takes place outside the planetary magnetosphere or ionosphere boundary, these newly produced ions can be easily picked up by the stellar wind and thus also contribute to the non-thermal losses. ENAs which impact the planetary atmosphere can also deposit their energy there and contribute to the atmosphere’s energy budget, but their influence on thermal escape has been estimated to be negligible (Lichtenegger et al. 2016).

In this article, we adopted the stellar wind data from Kislyakova et al. (2013). Kislyakova et al. (2013) applied a Versatile Advection Code (Tóth 1996) to model the plasma environment of a GJ436-like star. The model includes a self-consistent Parker-type corotating magnetic field and is based on the solution of the set of the ideal nonresistive nonrelativistic magnetohydrodynamic equations, and yields a self-consistent expanding stellar wind plasma flow. In the middle of the HZ (at 0.24 au), stellar wind parameters for a quiet case, which we adopted in the present study, were the following:

• •

Stellar wind density $n_{sw}=2.5\times 10^{8}$ m*-3*;

• •

Stellar wind velocity $v_{sw}=330$ km/s;

• •

Stellar wind temperature $T_{sw}=3.5\times 10^{5}$ K.

Here and below, the stellar wind densities are the proton densities. Assumed stellar wind parameters correspond to a total mass loss rate of $3.5\times 10^{-14}$ M⊙ year*-1*, or $\approx 2.5$ times the mass loss rate of the current Sun. This value is higher than the mass loss rate of GJ 436 estimated from the Ly$\alpha$ observations of GJ 436b by Vidotto & Bourrier (2017) of $(0.45-2.5)\times 10^{-15}$ M⊙ year*-1*. A higher mass loss rate leads to i) a higher electron impact ionization rate and ii) a higher amount of ENAs which are produced as a result of interaction between the atmospheric neutrals and incoming stellar wind. In Section 5.3, we investigated the influence of different stellar wind parameters on our results.

A lower electron impact ionization rate would decrease the ionization of neutral hydrogen already present in the atmosphere and, therefore, slightly increase the modeled absorption for hydrogen-dominated atmospheres. We used the database of energy-dependent absorption cross sections for electron impact ionization by the National Institute of Standards and Technology (NIST111https://physics.nist.gov/PhysRefData/Ionization/intro.html). We calculated the electron impact ionization rates used in the simulations as $\tau_{\rm ei}=\sigma(\epsilon)v_{\rm th}n_{\rm sw}$, where $\sigma(\epsilon)$ is the cross section for a given energy $\epsilon$ of a stellar wind, which is calculated according to its temperature as $\epsilon=3/2k_{B}T_{\rm sw}$, and $v_{\rm th}=\sqrt{3k_{B}T_{\rm sw}/m_{\rm el}}$ is the most probable Maxwellian velocity of a electron. Here, $T_{\rm sw}$ is the stellar wind temperature. Electron impact ionization affects only hydrogen atoms outside the magnetospheric-ionospheric boundary, while photoionization affects all atoms exposed to sunlight (not in the shadow of the planet).

2.4 Photoionization rates and the Ly$\alpha$ absorption rates

To calculate the photoionization rates we used the XUV spectrum of GJ 436 from the MUSCLES database222https://archive.stsci.edu/prepds/muscles (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016). Specifically, we used the panchromatic spectral energy distribution binned to a constant 1 Å resolution, version 1.2. The X-ray ($<$100 Å) part was derived from Chandra observations and the (unobservable) EUV (100-912 Å) part was reconstructed with the method of Linsky et al. (2014) from the star’s Ly$\alpha$ emission (Loyd et al. 2016). The spectra contain the apparent fluxes as observed at the Earth, so we adopted a distance of 10.14 pc (parallax of 98.61 mas; van Leeuwen 2007) and an orbital distance of 0.24 AU to scale the fluxes to the center of GJ 436’s HZ. For the photoionization cross sections $\sigma_{\mathrm{ph}}$ we used the analytic fits from Verner et al. (1996). We then calculated the photoionization rates $\beta_{\mathrm{ph}}=\int_{\lambda_{0}}^{\lambda_{1}}\sigma_{\mathrm{ph}}\phi_{\mathrm{XUV}}d\lambda$, where $\phi_{\mathrm{XUV}}$ is the XUV photon flux at 0.24 AU, $\lambda_{0}=5$ Å (minimum wavelength of the adopted spectrum), and $\lambda_{1}$ is the wavelength of the ionization threshold of the specific atom. The resulting photoionization rates were $1.434\times 10^{-7}$ s*-1* for H and $1.158\times 10^{-6}$ s*-1* for N. The adopted spectrum yielded a total XUV (5-912 Å) flux of 15.6 erg cm2 s*-1* at 0.24 AU from the star, which is about three times the present solar value at Earth (4.64 erg cm2 s*-1*; Ribas et al. 2005). The total XUV flux from MUSCLES is slightly lower than the value estimated by Ehrenreich et al. (2015), who derived a 35% higher XUV flux. More recently, King et al. (2018) estimated the XUV emission of GJ 436 to be 26% higher than from MUSCLES. Since the X-ray range does not contribute significantly to the photoionization rates because of the small cross sections, the rates would be raised by 20-28% if using these alternative XUV flux reconstructions. GJ 436 is variable in X-rays and likely also in EUV (King et al. 2018), so data taken at different epochs may yield varying XUV flux estimates and thus photoionization rates. The photoionization rates at 0.24 AU for a more active GJ 436-like star were calculated as follows. For a given Ly$\alpha$ flux, the XUV spectrum can be calculated following Linsky et al. (2013, 2014). Linsky et al. (2013) provide a scaling for the X-ray band ($<$100Å) and Linsky et al. (2014) for the EUV (100-912Å). For properly scaled Ly$\alpha$ fluxes, we computed the corresponding XUV fluxes by summing X-ray and EUV contributions to obtain 7, 10 and 20 XUV. The spectra matching these total XUV fluxes were then used to obtain the photoionization rates (Table 2) with cross-sections taken from Verner et al. (1996) and integrating over wavelength.

For obtaining the velocity-dependent Ly$\alpha$ absorption rate per atom $\beta_{\mathrm{abs}}(v)=\sigma_{\mathrm{abs}}\phi_{\mathrm{Ly\alpha}}(v)$ we used the reconstructed intrinsic line profile from the MUSCLES database (Youngblood et al. 2016). Here, $\phi_{\mathrm{Ly\alpha}}(v)$ is the Ly$\alpha$ photon flux at 0.24 AU, $\sigma_{\mathrm{abs}}=5.47\times 10^{-15}$ cm2 Å*-1* the total Ly$\alpha$ absorption cross section (Quémerais 2006), and $v=-c(\lambda/\lambda_{\mathrm{c}}-1)$ is the atom’s radial velocity (counted positive toward the star; $\lambda_{\mathrm{c}}=1215.67$ Å). The maximum absorption rate at the line peak is $9.2\times 10^{-3}$ s*-1*, falling off by a factor of ten within $\pm$120 km s*-1* from the line center. The total intrinsic Ly$\alpha$ flux from MUSCLES is in good agreement with Bourrier et al. (2015), but about 30-60% lower than the reconstructions shown in earlier studies (Ehrenreich et al. 2011; France et al. 2016). These differences can be attributed to slightly different reconstruction methods, as well as using data taken at different epochs. Calculated absorption rates for 3, 7, 10, and 20 XUV are shown in Fig. 1. To determine the UV absorption rates for XUV enhancement factors of 7, 10, and 20 we used the corresponding Ly$\alpha$ fluxes matching these enhancements from the method described above. We then scaled the present line profile of GJ 436 to these enhanced flux levels to obtain the corresponding UV absorption rates for the higher activity levels.

We did not consider acceleration by the radiation pressure for species other than H, because the stellar flux in the corresponding lines is much smaller and, therefore, species other than H are not accelerated efficiently. In this paper, we did not include self-shielding, meaning that we did not account for the optical depth of the hydrogen envelope in the Ly$\alpha$ line and assumed that all hydrogen atoms at all locations in the simulation domain (except in the optical shadow of the planet) have the same probability to scatter a stellar Ly$\alpha$ photon. This seems reasonable because in this paper we focused mostly on thin hydrogen envelopes (especially for nitrogen-dominated atmospheres). For hydrogen-dominated atmospheres, absence of self-shielding leads to excess acceleration of atmospheric particles away from the star. However, due to low radiation pressure in the habitable zone of GJ 436, this effect is rather minor and not so significant as for GJ 436b (see Fig. 1).

2.5 DSMC code

In this subsection, we describe the DSMC upper atmosphere-exosphere 3D particle model used to study the plasma interactions between the stellar wind of GJ 436 and the upper atmospheres of terrestrial planets. The software is based on the FLASH code written in Fortran 90 (Fryxell et al. 2000). The code follows particles in the simulation domain according to all forces acting on a hydrogen atom and includes up to three neutral species as well as their ions, with two of the species being neutral hydrogen and hydrogen ions (the latter also include stellar wind protons).

The comprehensive description of the model can be found in Holmström et al. (2008); Kislyakova et al. (2014a). The code has been successfully applied to modeling of hydrogen-dominated extrasolar giant (Kislyakova et al. 2014a) and terrestrial planets (Kislyakova et al. 2013, 2014b; Lichtenegger et al. 2016). The latest improvement of the code includes adding multiple species to it, which are allowed to charge exchange with stellar wind protons as well as being photo- and electron impact ionized and collide with each other. The main processes and forces included for an exospheric atom are:

For hydrogen atoms: collision with a UV Ly-$\alpha$ photon which defines the velocity-dependent radiation pressure. 2. 2.

Charge exchange with a stellar wind proton. 3. 3.

Elastic collision with another neutral atom. 4. 4.

Ionization by stellar photons or wind electrons. 5. 5.

Gravity of the star and planet (including the tidal forces), centrifugal, and Coriolis forces.

Summarizing this, our code includes three sources of ionization for atmospheric neutrals: photoionization, electron impact ionization, and charge exchange. Electron impact ionization and charge exchange can occur only outside the magnetospheric-ionospheric boundary, in the region, where the atmospheric neutrals come in contact and can collide with the stellar wind. While electron impact ionization is very important for hot Jupiters, in the HZ it is less significant due to a much less denser wind (Holzer 1977). If a planetary neutral is ionized, it is deleted from the simulation if it has been ionized inside the ionosphere and kept in the simulation, if the parental neutral particle has been ionized outside it. We assumed that all ions outside the ionosphere were picked up by the stellar wind and lost from the planetary atmosphere, thus contributing to ion pick up losses. The ions inside the ionosphere were deleted to shorten the simulation times. Our model does not include magnetic fields and collisions for the ions, therefore, the motion of the ions is governed only by the gravity forces of the planet and the star. If a simulation does not include atmospheric neutral hydrogen, it is produced only as a product of charge exchange between the incoming stellar wind protons and atmospheric neutral particles.

The magneto-ionopause of the planet was represented by a conic shaped obstacle:

[TABLE]

Here, $R_{\rm s}$ stands for the planetary obstacle stand-off distance and $R_{\rm t}$ for the width of the obstacle. We used a Cartesian coordinate system with the $x$-axis pointing toward the star (meaning that particles flying toward the star have positive velocities along the $x$-axis, and particles flying away from the star have negative velocities), $y$-axis pointing opposite to the planetary motion, and, finally, $z$-axis completing the right-hand coordinate system. Since the obstacle shape and location depend strongly on the planetary magnetic field strength, one can model the interaction of the stellar wind with magnetized as well as with non- or weakly magnetized planets by the appropriate choice of $R_{\rm s}$ and $R_{\rm t}$. In this work, we were interested to estimate the maximum possible absorption due to production of ENAs, tus we selected the values of $R_{\rm s}$ and $R_{\rm t}$ just above the exobase, which maximized the size of the region where the stellar wind protons could penetrate in the atmosphere and interact with it. The height of the exobase was determined based on results of MHD modeling. This choice is also justified by the fact that the pressure balance point between stellar wind ram pressure and the thermal pressure of the planetary atmosphere (not accounting for the RAM pressure of the planetary outflow and magnetic pressure, since we do not consider magnetic fields in this study) was located below the inner boundary.

The obstacle was rotated by an angle of $\arctan(v_{\rm pl}/v_{\rm sw})$, to account for the finite stellar wind speed relative to the planet’s orbital speed. This approach does not include polar openings (cusps), and as such can not reproduce all effects of a magnetosphere. However, it is a good choice for modeling of non-magnetized planets, on which we focused in this study. We assumed that the planets we consider lack intrinsic magnetic fields and experience a Venus-type interaction with the incoming stellar wind. In this case, the stellar wind plasma can directly impact the stellar atmosphere, leading to formation of a compressed and narrow ionospheric obstacle (Lundin 2011).

The code allows to track metaparticles of different weights. If charge exchange occurs between a stellar wind proton metaparticle and an atmospheric neutral particle of a different weight, then the bigger particle is split in two parts, one of which contains exactly the same amount of real atoms as the smaller metaparticle. The particle containing the excess atoms is excluded from the interaction and continues traveling with the same speed and trajectory. The particles with equal weights then charge exchange. In this paper, we include metaparticles of several different types (or species): neutral atomic hydrogen and nitrogen, ionized nitrogen of planetary origin, and two types of metaparticles for ionized hydrogen: of stellar and planetary origin.

After we modeled the hydrogen corona, we calculated the Ly$\alpha$ in-transit attenuation. Natural broadening was accounted for at this stage. Doppler broadening is included automatically by accurately taking into account the velocities of the atoms.

2.6 Calculation of the Ly$\alpha$ transmissivity

After a hydrogen cloud was simulated and by knowing the positions and velocities of all the hydrogen metaparticles at a certain time, we computed how these atoms attenuate the stellar Ly$\alpha$ radiation by using a post-processing software written in the Python programming language. To compute the transmissivity along the line-of-sight (LOS) we followed the approach of Semelin et al. (2007). The post-processing tools have been described in details in Kislyakova et al. (2014a). Here we repeat the main features of them.

To increase the particle statistics, especially for cases with a low number of metaparticles of a specific specie, particles were deposited in the cell according to the Cloud in cell (CIC) algorithm, where each particle is viewed as a cloud of particles with the same size as the cell. CIC was only applied for the spatial coordinates $(y,z)$. Each particle weight was distributed among four cells, which provided a self consistent smoothing of the solution.

We assumed that the test planet has zero inclination, therefore, it was located in the center of the stellar disk at mid transit. However, we took inclination into account for GJ 436b, which is a grazing planet and is transiting very close to the edge of the stellar disk. Only neutral hydrogen atoms absorb in the Ly$\alpha$ line. One has to take into account spectral line broadening. Real spectral lines are subject to several broadening mechanisms: i) natural broadening; ii) collisional broadening; iii) Doppler or thermal broadening. The “natural line width” is a result of quantum effects and arises due to the finite lifetime of an atom in a definite energy state. A photon emitted in a transition from this level to the ground state will have a range of possible frequencies: $\Delta f\sim\Delta E/\hbar\sim 1/\Delta t$. The distribution of frequencies can be approximated by a Lorentzian profile. Collisional broadening is caused by the collisions randomizing the phase of the emitted radiation. This effect can become very important in a dense environment, yet above the exobase it does not play a role and is important only in the lower parts of the atmosphere. The third type of broadening, which plays a significant role in the upper atmosphere of a hot exoplanet, is thermal broadening. If $f_{0}$ is the centroid frequency of the absorption line, the frequency will be shifted due to the Doppler effect. Combining Doppler shift with the Maxwellian distribution of $v_{x}$, one can obtain a Gaussian profile function, which is decreasing very rapidly away from the line center. The combination of thermal and natural (or collisional) broadening is described by the Voigt profile, which is the convolution of the Lorentz and Doppler profiles.

In the considered cases, an analytical solution for the absorption profile cannot be obtained, since it is not only thermal atoms that contribute to the broadening. The presence of a non-thermal population of hot atoms (ENAs and atoms accelerated by the radiation pressure) changes the picture. Mathematically it means that the line width cannot be described by the Voigt profile anymore. We calculated the natural broadening for all atoms and bin it by velocity, which automatically gives us the Doppler broadening for a particular velocity distribution.

The velocity spectrum of an atomic cloud was then converted to frequency via the relation $f=f_{0}+v_{x}/\lambda_{0}$ with $f_{0}=c/\lambda_{0}$, $\lambda_{0}=1215.65\times 10^{-10}$ m. To compute the transmissivity along the line-of-sight (LOS) we followed the approach of Semelin et al. (2007), where the relation between the observed intensity $I$ and the source intensity $I_{0}$ as a function of frequency $f$ could be written as

[TABLE]

Here, $\tau=\sigma(f)Q$ is the frequency dependent optical depth, $Q$ is the column number density of hydrogen atoms and $\sigma(f)$ is the frequency dependent crossection, which depends on the normalized velocity spectrum, the Ly$\alpha$ resonance wavelength and the natural absorption crossection in the rest frame of the scattered hydrogen atom (Peebles 1993). The quantity $T$ is called the transmissivity. For a hydrogen cloud in front of a star, the transmissivity of the stellar spectrum was computed on an $yz$-plane which was discretized into a grid with $N_{c}$ cells. For each cell in the grid along lines of sight in front of the star ($y^{2}+z^{2}<R_{\rm st}^{2}$), we calculated the velocity spectrum of all hydrogen atoms in the column along the $x$-axis. Then we calculated the transmissivity averaged over all columns in the $yz$-grid except those particles which fell outside the projected limb of the star or inside the planetary disk. The average transmissivity was then calculated as

[TABLE]

where $T_{i}$ is the transmissivity in a cell. For lines-of-sight in front of the planet $(y-y_{p})^{2}+(z-z_{p})^{2}<R_{\rm pl}^{2}$, where ($y_{p}$ , $z_{p}$) is the planet center position, we set $T_{i}=0$ (zero transmissivity). Then, we appied the average transmissivity to the observed out-of-transit spectrum $I_{0}$ yielding the modeled in-transit spectrum

[TABLE]

The frequency dependent crossection was defined by:

[TABLE]

where $\check{u}(v_{x})$ is the normalized velocity spectrum along the LOS, so that $\int_{-\infty}^{+\infty}\check{u}(v_{x})dv_{x}=1$, $c$ is the speed of light. We assumed the natural absorption crossection in the rest frame of the scattered hydrogen atom according to Peebles (1993)

[TABLE]

Here, $f_{\alpha}=c/\lambda_{\alpha}$ with $\lambda_{\alpha}$ being the Ly$\alpha$ resonant wavelength and $A_{21}=6.265\times 10^{8}$ s*-1* the rate of radiative decay from the $2p$ to the $1s$ energy level. Multiplying $\sigma(f)$ with Q, we computed the optical depth directly without normalizing the velocity spectrum:

[TABLE]

The expression for the natural absorption crossection given in Eq. 9 was approximated by the Lorentzian profile

[TABLE]

where $f_{12}=0.4162$ is the Ly$\alpha$ oscillator strength, $e$ is the elementary charge, $m_{e}$ is the electron mass, and $\Delta f_{L}=9.936\times 10^{7}$ [s*-1*] is the natural line width (Semelin et al. 2007). We used the crossection from Eq. 11 to include the contribution of broadening to the absorption. To account for the contribution of the lower atmosphere, for a case of a hydrogen-dominated atmosphere we added a Maxwellian velocity spectrum corresponding to a hydrogen gas with a specified column density and temperature to all pixels inside the inner boundary $R_{\rm ib}$.

3 Model validation. GJ 436b.

In this section, we validate our model by reproducing the Ly$\alpha$ transit observations of a “warm Neptune” exoplanet GJ 436b. Ly$\alpha$ transit observations by the Hubble Space Telescope revealed a strong absorption of about 8% at mid-transit and 33% post-transit (Kulow et al. 2014). A later result which used a corrected ephemeris revealed an even deeper transit depth of $\approx$56% in the Ly$\alpha$ line (Ehrenreich et al. 2015), with ultraviolet transits repeatedly starting about two hours before, and end more than three hours after the approximately one hour optical transit. For comparison, in the optical wavelengths the transit depth is only about 0.69 %. Later observations (Lavie et al. 2017) confirmed the presence of a deep in-transit absorption, and also showed that the Ly$\alpha$ transit started about two hours before the optical transit, and lasted for more than 10 hours after the end of the optical transit. Such deep absorption in the Ly$\alpha$ line is observed even despite GJ 436b grazing very close to the edge of the stellar disk, thus indicating a presence of a very extended hydrogen envelope enwrapping the planet. The time evolution of the observed Ly$\alpha$ emission line of GJ 436 before, during, and after the planetary transit is rather complex (Lavie et al. 2017); previously, significant efforts have been undertaken to model the observed phenomena (Bourrier et al. 2015, 2016). Here, we prove that our model was also capable of reproducing the observed mid-transit absorption of GJ 436b.

We applied the model described in the previous section to GJ 436b. First, we model the lower atmosphere using the hydrodynamic code described in Section 2.2. To calculate the lower atmosphere velocity, density, and temperature profiles, we assumed the lower atmosphere of GJ 436b to be comprised of atomic hydrogen only. Close to the planetary radius, the main atmospheric constituent is molecular H2, which we did not take into account. However, sophisticated 3D modeling shows that the H2 dissociation front is reached already at 1.35$R_{\rm pl}$, therefore, atomic hydrogen is the main atmospheric component thorough most of our simulation domain for the lower atmosphere (Shaikhislamov et al. 2018). Some works (Hu et al. 2015) have shown that planets like GJ 436b might accumulate a significant amount of helium due to faster atmospheric escape of hydrogen. The presence of helium can significantly influence the mass loss rate and the atmospheric structure (Shaikhislamov et al. 2018). However, although our hydrodynamic model for the lower atmosphere is 1D and not 3D and although we did not take into account the heavy species in this modeling, the model provided very similar values for hydrogen number density, temperature, and outflow velocity at 5.83 Rpl in comparison to the 3D multiple-species model (Shaikhislamov et al. 2018). We assumed the stellar wind parameters at GJ 436b’s orbit as estimated by Vidotto & Bourrier (2017) except of a slightly faster stellar wind, which, however, still fell well inside the range of the mass loss rates predicted for GJ 436 by (Vidotto & Bourrier 2017).

Our results confirmed that GJ 436b possesses an extended hydrogen envelope, the presence of which can be anticipated from the deep Ly$\alpha$ in-transit absorption signature and has been confirmed by previous modeling (Bourrier et al. 2016). Due to relatively low gravity of this planet, the exobase is located higher than the Roche lobe, above which the atmosphere is not bound to the planet gravitationally. This means that the atmosphere stays collisional up until the Roche lobe and beyond. According to our simulations, the Roche lobe and the exobase were located at 5.83 Rpl and 7.67 Rpl, respectively. We started our DSMC simulation at the Roche lobe of the planet, assuming the density, temperature, and outflow velocity of the atomic hydrogen gas obtained in the hydrodynamic simulation (described in Section 2.2) as an input at the spherical inner boundary. Table 1 summarizes simulation input parameters for the DSMC simulation. As in all other simulations, we took into account electron impact and photoionization, elastic collisions between the atmospheric particles, and charge exchange collisions with the stellar wind protons. We accounted for the radiation pressure at the orbital location of GJ 436b (Fig. 1; the data for 3 XUV in the habitable zone can be scaled according to the closer orbital distance of GJ 436b), but we did not take into account the self-shielding, which means the optical depth of the atmosphere for the Ly$\alpha$ radiation (Bourrier et al. 2016). Fig. 2 shows the modeled atomic corona around GJ 436b. Fig. 3 shows the calculated mid-transit Ly$\alpha$ absorption in comparison with observations. A test run without radiation pressure (not shown here) yielded similar result to the one shown in Fig. 3 with only slight overabsorption in the red wing of the Ly$\alpha$ line, which indicates that the atmosphere is escaping very efficiently due to gravitational forces. This seems reasonable, because due to the proximity to its host star, GJ 436b has a Roche Lobe located at only 5.83 Rpl, with a significant fraction of the atmosphere being above it. After the particles leave the planet’s gravitational well, they are no longer bound to the planet gravitationally and escape, forming a long tail curved by the Coriolis force (Fig. 2). We could reproduce the observations assuming a lower velocity of the planetary wind in comparison to earlier works (Bourrier et al. 2016), which assumed the outflow velocities of the planetary gas of 50–60 km s*-1*. In our model, we assumed an outflow velocity of 7.9 km s*-1* at the inner boundary, which was the gas velocity at this distance to the planet in accordance to our hydrodynamic model applied to the lower atmosphere. We could obtain a good fit to the observations also assuming the stellar wind speed of 85 km/s as estimated by Vidotto & Bourrier (2017), however, we got a slight overabsorption on the right part of the Ly$\alpha$ line in comparison to 110 km/s which is shown in Fig. 3.

We have calculated the Ly$\alpha$ in-transit signature of GJ 436b at mid-transit and compared it to the observed mid-transit signature by Lavie et al. (2017). Our model agrees well with the observations both on the left and right wings of the Ly$\alpha$ line with a overabsorption in only one observational point on the right wing, therefore confirming that out model can be used to estimate the in-transit absorption depth in the Ly$\alpha$ line of exoplanets. Detailed modeling of the Ly$\alpha$ observations of GJ 436b will be published in a following article (currently in preparation).

4 Simulation parameters for terrestrial planets

In this section, we summarize our simulation parameters which we have used as inputs for the DSMC code. In all cases, the planet has been assumed to have a mass and radius equal to the Earth, and to orbit in the middle of the HZ of GJ 436 at 0.24 au. The location of the habitable zone has been calculated as in Kislyakova et al. (2013).

Stellar wind parameters have been chosen according to Kislyakova et al. (2013) and have been the same in all simulations, as described in Section 2.3. Since we assumed a constant stellar wind, the ionization rate by stellar wind electrons is also constant in all simulations. Electron impact ionization rates have been caclulated according to Banks & Kockarts (1973). All simulation parameters which were kept constant for all XUV fluxes are summarized in Table 2.

In this paper, we were interested in modeling of the strongest possible interaction between the stellar wind and the planetary atmosphere. For this reason, we focused on non-magnetized planets. If the planet has no intrinsic magnetic field, we can assume that the magnetospheric-ionospheric obstacle lies very close to the inner atmosphere boundary and is very narrow (similar to the ionospheric boundary observed at Venus; Galli et al. 2008a; Lundin 2011). In our model, including an intrinsic magnetic field reduces the size of the region where the stellar wind can interact with the planetary atmosphere, thus reducing the number of produced ENAs and the ENAs absorption in the wings of the Ly$\alpha$ line. We did not model the polar cusp regions where the stellar wind can interact with neutral atmospheres of magnetized planets. Around the Earth, ENAs have been observed both at the magnetopause region (Fuselier et al. 2010) and in the cusps (Petrinec et al. 2011). Since the latter region of ENAs formation is typical only for a planet with an intrinsic magnetic field, we focus only at the interaction at the substellar region of the planetary magnetospheric-ionospheric boundary. Both the atmosphere and the ENAs contribute to the total in-transit absorption in the Ly$\alpha$ line. The atmosphere contribute mostly through the spectral line broadening (e.g., Ben-Jaffel 2007). Atmospheric particles can also be accelerated up to high velocities by the radiation pressure of the host star and increase the absorption in the left wing of the Ly$\alpha$ line. This mechanism is extremely effective for close-in planets such as HD 209458b or, to a lesser extent, GJ 436b (e.g., Bourrier & Lecavelier des Etangs 2013; Bourrier et al. 2016), but is of less importance for planets in the habitable zone, where the radiation pressure is naturally weaker due to the larger orbital distances (see Fig. 1). In this paper, we did account for the atmosphere contribution via line broadening and for the radiation pressure. For a given radiation pressure and an atmospheric density profile, the maximum amount of ENAs is produced at a non-magnetized planet.

Table 3 summarizes our atmospheric parameters for all simulations. As one can see, a nitrogen-dominated atmosphere expands significantly when the XUV flux is increased from three to seven times the modern Earth level. The exobase temperature increases from 3 to 7 XUV and then drops again due to adiabatic cooling caused by atmospheric expansion. Hydrogen-dominated atmospheres also expand if the XUV flux is increased. We always chose the inner boundary of our simulation, $R_{\rm ib}$, equal to the exobase height estimated in 1D simulations by Lammer et al. (2008) for nitrogen-dominated atmospheres and the 1D model described in section 2.2 adopted from Erkaev et al. (2017) for hydrogen-dominated atmospheres.

5 Results.

5.1 Modeling of atomic coronae around Earth-like planets.

Here, we present the results of our modeling of planetary exospheres of planets orbiting in the center of the habitable zone of GJ 436 (at 0.24 au). As stated earlier, we modeled on an Earth-like planet with nitrogen-dominated and hydrogen-dominated atmospheres at 3, 7, 10, and 20 XUV. At 3 XUV, we additionally model a nitrogen-dominated atmosphere with a terrestrial hydrogen content.

For every atmosphere, the velocity distribution of neutral hydrogen is modified by the stellar radiation pressure, which accelerates the atoms in the direction away from the star according to Fig. 1 (however, in the habitable zone this effect is less impotant than for hot Jupiters or for a close-in planet like GJ 436b). We assume that neutral nitrogen atoms are not accelerated by the stellar radiation, as for them this effect is negligible. The distribution of neutral atoms (both hydrogen and nitrogen) is also altered by photo- and electron impact ionizations, which predominantly remove particles with positive velocities along the x-axis (flying toward the star), as they are more probable to leave the shadow of the planet and the ionosphere boundary and undergo ionization.

Figs. 4, 5, and 6 show the simulated exospheric atomic clouds around the planet. All figures show all particles within a slice with coordinates $-10^{6}\leq z\leq 10^{6}$ m. The ionospheric boundary is not explicitly shown, but can be seen as an area close and behind the planet devoid of stellar wind protons.

Fig. 4 illustrates the appearance of the atomic clouds for 3 XUV for three different atmospheric compositions: a nitrogen-dominated atmosphere, a nitrogen-dominated atmosphere with a modern Earth hydrogen fraction, and, finally, a hydrogen-dominated atmosphere. As one can see, the hydrogen atmosphere is very expanded even if exposed to this comparably low XUV level. Fig. 5 shows three exospheric clouds simulated assuming a nitrogen-dominated composition exposed to 7, 10, and 20 XUV. As one can see, a pure nitrogen atmosphere exposed to XUV fluxes higher than $\sim$5–7 times the modern Earth value heats up efficiently, and, in absence of strong infrared coolants, the lower atmosphere expands to high altitudes, so that the exobase is located at very high altitudes in comparison to the 3 XUV case. In Fig. 5, the lower atmosphere is shown as a white empty area surrounding the black dot, the planet. High exobase altitudes lead to intense charge exchange between the stellar wind and the atmosphere due to the increased planetary cross section, and the production of a much higher amount of ENAs (see velocity distribution of neutral hydrogen atoms in section 5.2). Finally, Fig. 6 shows our modeling results for a hydrogen-dominated planet exposed to 7, 10 and 20 XUV. At a given XUV flux, a hydrogen-dominated atmosphere is even more expanded than a nitrogen-dominated one, which is due to the low mass of a hydrogen atom. Summarizing the results shown in Figs. 4, 5 and 6, we can conclude that both hydrogen and nitrogen-dominated atmospheres react with expansion to an enhanced level of stellar short wavelength radiation.

5.2 Modeling of in-transit absorption in the Ly$\alpha$ line.

In this section, we calculated the absorption in the Ly$\alpha$ line which would be produced if a planet surrounded by exospheric clouds shown in the previous section would transit in front of GJ 436. The method of calculating the Ly$\alpha$ transmissivity has been presented in section 2.6.

First, we calculated the velocity distribution of all neutral hydrogen atoms which were present in the simulation domain at the end of the simulation. For a hydrogen-dominated atmosphere, hydrogen atoms can be divided into two big groups: neutral exospheric particles of planetary origin with the velocity distribution slightly modified by the radiation pressure and ionization, and energetic neutral atoms produced as a result of the charge exchange between stellar wind protons and neutral particles. ENAs have a different velocity distribution in comparison with the planetary particles, since they initially keep the velocity of a stellar wind particle. For this reason, ENAs are much more energetic than the atmospheric particles. As we show below, ENAs contribute to a high velocity particle population. In case of nitrogen-dominated atmospheres, ENAs are the only source of exospheric hydrogen. For the case of 3 XUV we additionally modeled an atmosphere with an Earth-like hydrogen fraction, which in case of the Earth is produced by water photodissociation.

Figs. 7, 8, and 9 show velocity distributions for all H atoms in the simulations along the $x$-axis, which is directed toward the star. Fig. 7 shows velocity spectra for the case of 3 XUV corresponding to clouds shown in Fig. 4, from left to right: a nitrogen-dominated atmosphere, a nitrogen-dominated atmosphere with a terrestrial hydrogen fraction, and a hydrogen dominated atmosphere. Fig. 8 shows the H-atoms velocity distribution for a nitrogen-dominated atmosphere for 7, 10, and 20 XUV, and finally Fig. 9 shows the same for a hydrogen-dominated atmosphere.

Comparing the left panel of Fig. 7 to all cases shown in Fig. 8, one can see that the amount of hydrogen particles produced due to charge exchange increases significantly between 3 and 7 XUV, which indicates a more intense interaction. This coincides with an increase in the cloud size, therefore, it seems reasonable to suggest that the increase in the amount of produced H atoms is due to increase of the atmosphere cross section. Indeed, hydrodynamical models show that a nitrogen-dominated atmosphere with an Earth-like amount of CO2 is very sensitive to the flux in the XUV radiation domain (Tian et al. 2008; Johnstone et al. 2018). If one compares all three panels of Fig. 8 to each other, one can see that the amount of neutral hydrogen is gradually increasing with increasing XUV fluxes.

Comparing Figs. 8 and 9 to each other, one can see that the main central peak, which presents atmospheric H particles with small velocities and which is present for all hydrogen-dominated atmospheres is actually absent in case of nitrogen-dominated atmospheres, where one can see only the ENAs distribution. In Fig. 9 one can see that the central peak is actually slightly shifted toward the negative velocities. This is due to the combined effect of two factors, first, the radiation pressure, which accelerates hydrogen atoms away from the star, and second, due to enhanced ionization in the part of the cloud directed toward the star, where the particles are more likely to have positive velocities along the $x$-axis. In our model, ionization predominantly removes atoms with positive velocities, because electron impact ionization can occur to particles only outside the ionosphere, and photoionization is able to ionize only the particles outside the planetary shadow. This leads to some imbalanced removal of particles flying toward the star.

In all velocity spectra, one can see that the initial velocity of the ENAs is altered by collisions with other atmospheric neutrals. Initial ENAs velocities correspond to the stellar wind velocities, which present a Maxwellian distribution with a temperature of $2\times 10^{6}$ K and a maximum at 330 km/s. On the other hand, Figs. 7–9 show that there is some population of neutral H atoms with the positive velocities (flying toward the star), as well as a small maximum (in case of nitrogen-dominated atmospheres) at small velocities. This is due to the fact that collisions lead to i) isotropic scattering in a random direction and ii) to some energy exchange between ENAs and atmospheric particles.

Figures 10 and 11 show our final results, namely, Ly$\alpha$ in-transit spectra compared to out-of-transit observation of GJ 436 (Lavie et al. 2017). Simulated absorption along the line-of-sight shows similar features with the velocity spectra, as it is created by atoms with corresponding velocities. Fig. 10 shows in-transit absorption for a nitrogen-dominated atmosphere with the terrestrial hydrogen content at 3 XUV (left) and a nitrogen-dominated atmosphere at 20 XUV (right). As one can see, both cases reveal no significant in-transit absorption, even the case with a small amount of hydrogen present. One can see that it can produce a very weak absorption signature at low velocities, however, due to the low orbit of the Hubble Space Telescope this part of the spectrum is obscured by the geocorona and cannot be observed. The nitrogen-dominated atmosphere at 20 XUV produces some in-transit absorption due to ENAs in the blue and red wings of the Ly$\alpha$ line, however, the signal is almost negligible in comparison to the observational errors (shown as errorbars on the picture). We do not show simulated absorption for nitrogen-dominated atmospheres at 3, 7, and 10 XUV, because they produce even smaller signatures than the nitrogen-dominated atmosphere at 20 XUV, which are barely recognizable by eye.

Finally, Fig 11 summarizes our results for the in-transit absorption of hydrogen-dominated atmospheres. It shows calculated in-transit absorption signatures if a planet with a hydrogen-dominated atmosphere is exposed to 3 (upper left), 7 (upper right), 10 (lower left), and 20 (lower right) times the XUV flux of the modern Earth. As anticipated, absorption increases from 3 to 20 XUV. In the center of the Ly$\alpha$ line all atmospheres produce a deep absorption signature, however, this part of the line is contaminated by the geocorona. In the high velocity wings of the line, absorption increases with increasing XUV flux, being the strongest in the blue wing of the line at 20 XUV. One can see that in all cases the absorption exceeds the observational uncertainty and, therefore, can be considered detectable by the Hubble Space Telescope. However, the data by Lavie et al. (2017) have very low observational errors because they combine several observational visits. If one compares the modeled absorption to errorbars for just one visit (for instance, one of the visits by Ehrenreich et al. 2015), the in-transit signature only exceeds the observational errors for a hydrogen-dominated atmosphere exposed to 20 XUV.

5.3 Influence of stellar wind parameters

In this section, we investigated the influence of the assumed stellar wind parameters on our simulation results. The wind we assumed for the simulations corresponds to a mass-loss rate of $3.5\times 10^{-14}M_{\odot}$ year*-1*, which is approximately 2.5 times higher than the mass loss rate of the Sun, $\dot{M}_{\odot}$. This mass loss rate is approximately 30 times higher than the average value estimated for GJ 436 by Vidotto & Bourrier (2017). This discrepancy can indicate that additional studies are necessary on the stellar wind models, ideally coupled with the Ly$\alpha$ observations of astrospheres and planets (e.g, Wood et al. 2005). In this section, we performed two additional simulations to estimate the influence of the stellar wind on our results: first, we adopted the wind parameters from Vidotto & Bourrier (2017) and second, we tested the influence of a wind with lower temperature, lower velocity, and larger density than the ones assumed above. We assumed a hydrogen atmosphere exposed to 20 times the modern Earth’s XUV flux, with the deepest simulated in-transit Ly$\alpha$ absorption, as a proxy.

Stellar wind according to Vidotto & Bourrier (2017). For this simulation, we adopted the stellar wind proton density and velocity for the habitable zone from Fig. 3 of Vidotto & Bourrier (2017):

• •

Stellar wind density $n_{sw}=1.0\times 10^{7}$ m*-3*;

• •

Stellar wind velocity $v_{sw}=270$ km/s;

• •

Stellar wind temperature $T_{sw}=3.95\times 10^{5}$ K.

The temperature is adopted from Kislyakova et al. (2013), as before. For this simulation, we changed the stellar wind electron impact ionization rate according to the lower wind density to $\tau_{ei}=2.7\times 10^{-7}$ s*-1*. This stellar wind corresponds to a mass-loss rate of $1.2\times 10^{-15}M_{\odot}$ year*-1*. The left panel of Fig. 12 shows modeled in-transit absorption in the Ly$\alpha$ line. The absorption due to the ENAs decreased due to a lower flux of the stellar wind. The absorption due to atmosphere broadening almost did not change in comparison to the default 20 XUV case (Fig. 11, lower right panel), although it is marginally larger due to lower electron impact ionization.

Simulation for a denser and slower stellar wind. For this case, we assumed the following parameters:

• •

Stellar wind density $n_{sw}=8.25\times 10^{8}$ m*-3*;

• •

Stellar wind velocity $v_{sw}=100$ km/s;

• •

Stellar wind temperature $T_{sw}=4.0\times 10^{5}$ K.

These parameters still correspond to the same mass loss rate as in the standard stellar wind used in the simulations. Stellar wind electron impact ionization rate was assumed equal to $\tau_{ei}=2.3\times 10^{-5}$ s*-1*. Since a denser wind mowing with lower velocity can produce more ENAs in the observed velocity domain (-200 $\leq v_{x}\leq$ 200 km/s), one could expect a deeper absorption signature in comparison to the default case. However, right panel of Fig. 12 shows that this is not the case. This counterintuitive result is explained by the fact that the electron impact ionization rate is proportional to the wind density and is, therefore, higher for a denser wind. For this reason, although more ENAs are produced due to charge exchange with a denser wind, also more hydrogen atoms are ionized, which counterbalances the ENAs effect on the Ly$\alpha$ transit signature, and the absorption is only marginally increased.

Fig. 12 shows qualitatively the same absorption as Fig. 11 for the 20 XUV case. This indicates that in the considered case the Ly$\alpha$ absorption is to a large extent influenced by the planetary atmosphere, and to a lesser extent by the ENAs and the interaction with the stellar wind. A denser stellar wind can also influence the absorption of magnetized planets, because through its RAM pressure it is influencing the location of the substellar point and the general magnetosphere shape which, in turn, determines which portion of the atmospheric neutrals is exposed for the interaction with the wind (see, e.g., Kislyakova et al. 2014a for discussion). However, the study of magnetized planets is beyond the scope of the present article.

6 Discussion

As we have shown above, at the distance of $\sim$10 pc (the distance to GJ 436) only hydrogen-dominated atmospheres would be detectable around Earth-like planets in the HZ of M dwarfs using HST, and only in case of high XUV fluxes or increased observation times for stars with lower activities. Nitrogen-dominated atmospheres are non-detectable even at 20 XUV. According to our results, an Earth-like exosphere content also does not produce an observable Ly$\alpha$ signature. This result would likely be the same for other secondary atmospheres such as dominated by carbon dioxide. In fact, since CO2 is an efficient thermosphere coolant, it effectively decreases the temperature of the upper atmosphere thus protecting it from the expansion (Johnstone et al. 2018). A less expanded atmosphere would present a smaller cross section for the interaction with the stellar wind, therefore reducing the amount of hydrogen atoms produced by charge exchange reactions. Taking this into account, we can conclude that at 10 pc and beyond HST is able to detect only hydrogen-dominated atmospheres. This is in agreement with recent observations including non-detections of Ly$\alpha$ signatures by HST on small, presumably rocky, exoplanets (Ehrenreich et al. 2012; Bourrier et al. 2017b), possibly indicating that these planets either possess secondary (dominated by other species than hydrogen) atmospheres or are airless bodies due to efficient non-thermal escape. As we have shown in this paper, HST is not able to distinguish between these two possibilities for a target at 10 pc. A more advanced instrument such as LUVOIR or WSO UV (due to its higher orbit above the geocorona, which would allow to detect also the core of the Ly$\alpha$ line) might be able to do it.

One should keep in mind that we always compare the modeled absorption to the observed out-of-transit Ly$\alpha$ line of GJ 436, which corresponds to the XUV flux in the HZ of $\approx$3 XUV. On the other hand, our simulations have been performed for various XUV enhancement factors (3, 7, 10, and 20 XUV). If the star would be more active in XUV its Ly$\alpha$ line would also be stronger (see Fig. 1), in which case the errorbars of the observations would likely be smaller for the same observational setup.

Recently, Castro et al. (2018) also studied if exospheres of Earth-like planets in the HZs of M dwarfs can be observable by World Space Observatory UV (WSO UV). They did not perform a modeling of the upper atmosphere and used a more simplified approach, however, still arriving at similar conclusions. For an Earth-like exosphere (a nitrogen-dominated upper atmosphere with a small hydrogen content) and an Earth-like XUV flux they have assumed a modified barometric formula which accounts for a zenith angle. For higher XUV fluxes, they have adopted a hydrogen-dominated 1D atmosphere profile from Erkaev et al. (2013), calculated using a similar version of the code we used to obtain atmospheric profiles for hydrogen-dominated atmospheres in the current article. For both cases, Castro et al. (2018) did not take into account the interaction between the stellar wind and the neutral atmosphere. Even if a planet possesses a magnetic field, it may not protect the atmosphere from the dense and fast stellar wind of M dwarfs, which can sufficiently compress the atmosphere (e.g., Vidotto et al. 2013) and lead to an interaction pattern similar to the one of a non-magnetized planet considered in this article.

Despite different approaches, our results and the results of Castro et al. (2018) are in agreement. According to Castro et al. (2018), an Earth-like exosphere can be detected only as close as 1.35 pc, which is only slightly farther than the distance to the nearest star to the Solar system, Proxima Centauri (1.295 pc). This agrees with our conclusion that Ly$\alpha$ signatures of nitrogen-dominated atmospheres of a planet orbiting in the HZ of GJ 436 at 10.2 pc would be non-detectable by the HST. On the other hand, Castro et al. (2018) concluded that hydrogen-dominated atmospheres may be detectable by LUVOIR at further distances of about 10 pc, similar to our conclusion that HST may detect a hydrogen-dominated atmospheres of an Earth-like planet in the HZ of a star similar to GJ 436.

Therefore, there is a possibility that HST might be able to observe Earth-like exospheres at very close distances, such as that of Proxima Centauri. However, Proxima Centauri b is not a transiting planet (Anglada-Escudé et al. 2016), thus rendering these observations impossible. It is possible that the HST can detect hydrogen-dominated atmospheres around other nearby M dwarfs with low XUV fluxes within 5–10 pc. The TESS satellite, which is an all-sky survey which will detect transiting exoplanets orbiting the closest M dwarfs, may provide promising targets for observations by HST and, in the future, LUVOIR and WSO UV.

A recent study (dos Santos et al. 2018) estimated the observability of an Earth-like exosphere by LUVIOUR, coming to the conclusion that an Earth-like geocorona would be detectable by LUVOIR for planets orbiting stars not further than about 15 pc. They also arrived at the conclusion that at such distances the geocorona would not be detectable by the HST, which is in full agreement with our conclusions. In the current study, we have additionally considered different atmospheric structures and XUV enhancement factors, showing that the HST may detect a hydrogen atmosphere of an Earth-like planet orbiting an active nearby star, thus indicating that such observations are possible at present with the HST. However, our conclusions also point out that only future high-precision instruments such as LUVOIR could detect a hydrogen corona similar to Earth’s geocorona.

In this paper, we assumed quite large elastic collision cross sections and an isotropic scattering of atoms after a collision. In reality, cross sections decrease with increasing energy of the particles (e.g., Izmodenov et al. 2000) and the direction of the scattering is also non-isotropic (Izmodenov et al. 2000; Lindsay et al. 2005). This simplified assumption leads to an overestimate of ENAs in the red part of the Ly$\alpha$ line, because these particles with positive velocities (flying toward the star) are produced due to scattering on the neutral atmospheric particles. However, as we have shown in our results, for a terrestrial planet even in this case the ENAs do not produce a significant Ly$\alpha$ signature in the red part of the spectrum. Therefore, different collision cross sections would not change our conclusions. Implementing a velocity-dependent collision cross section and also taking into account the dependence of the scattering direction on the impact angle will be taken into account in future studies.

In our model, we did not take into account the influence of the interplanetary magnetic field (IMF) on the planetary ionosphere, therefore, our results are valid for a case of a small IMF, when the protons of the stellar wind can freely penetrate the planetary ionosphere. IMF can form a pile up region near the ionosphere (Erkaev et al. 2017), which will prevent the stellar wind from the interaction with the ionosphere, thus making the interaction region and the amount of generated ENAs smaller. However, since we were interested in estimating the Ly$\alpha$ signatures for the maximum interaction, but still did not find any significant absorption in the Ly$\alpha$ line for the nitrogen atmospheres, our main conclusions would not change for a case with a non-negligible IMF. For hydrogen-dominated atmospheres, the presence of an IMF may reduce the number of ENAs and the absorption in the left part of the spectrum, however, it would not affect the absorption caused by broadening of the spectral line by the atmosphere, or possibly even increase it due to decreased ionization by charge exchange and electron impacts.

We also do not take into account the effects of the formation of a bow shock and assume that stellar wind density and velocity do not change near the planetary ionospheric obstacle. Bow shock formation can lead to deceleration and temperature and density increase of the stellar wind near the planet. ENAs are mostly formed in the region between the ionopause and bowshock (Khodachenko et al. 2017). Taking into account these effects would slightly increase the amount of ENAs inside the observable velocity domain (from -200 to 200 km/s), but it is unlikely that even an increased ENA contribution would change our conclusion about observability of the nitrogen atmospheres in the Ly$\alpha$ line. As for hydrogen atmospheres, although ENAs can be a substantial contributor to the absorption signature, for the considered cases most of the absorption is produced due to line broadening, similar to hot Jupiters (Shaikhislamov et al. 2016). Therefore, it is unlikely that bow shock effects would alter our conclusions.

As a final note, we would like to repeat that nitrogen-dominated atmospheres may not be very common on planets orbiting M dwarfs. As we have stated before, nitrogen-dominated atmospheres with low CO2 ratios are subject to enhanced expansion and non-thermal losses if exposed to XUV levels higher than 5–7 (Tian et al. 2008; Lichtenegger et al. 2010). Even though GJ 436 does not exhibit a high level of XUV radiation in its habitable zone at present, it may have been much more active in the past. An increased level of CO2 may prevent the atmosphere from expansion and thus protect it (Kulikov et al. 2007; Lichtenegger et al. 2010; Johnstone et al. 2018), but this would require much higher CO2 levels in comparison to the present-day Earth, or a very inactive M dwarf thorough its history. Another argument which would suggest that nitrogen-dominated atmospheres may be rather rare is that Earth’s, its nitrogen atmosphere is sustained also due to plate tectonics and, indirectly, via influence of biologic activity (Lammer et al. 2018).

7 Conclusions

In this paper, we have performed modeling of the exospheres of an exoplanet orbiting in the HZ of GJ 436. We have assumed that the planet has a radius and a mass equal to those of the the Earth, and have considered hydrogen and nitrogen dominated atmospheres exposed to different XUV fluxes. We have applied sophisticated modeling which included 1D lower atmosphere models and a multi-species DSMC upper atmosphere code, and modeled atomic coronae surrounding the planets at 3, 7, 10, and 20 XUV. At 3 XUV, which is the current level of short-wavelength radiation in the middle of the HZ of GJ 436, we also additionally modeled a nitrogen-dominated atmosphere with an Earth-like hydrogen content. For all modeled atomic coronae, we have calculated the absorption in the Ly$\alpha$ line that a transit of such an atmosphere would produce, and compared it to the observed out-of-transit spectrum of the Ly$\alpha$ line of GJ 436 (Lavie et al. 2017). In case of hydrogen-dominated atmospheres, the signature in the Ly$\alpha$ line is produced due to the combined effects of spectral line broadening in the massive planetary atmosphere and the energetic neutral atoms generated due to charge exchange between stellar wind and the planetary neutrals. In case of a secondary atmosphere, which we assumed to be dominated by atomic nitrogen, the Ly$\alpha$ signature is produced due to ENAs alone.

According to our results, only hydrogen-dominated atmospheres of an Earth-like planet can be detected by the HST around GJ 436. At 3 XUV, a hydrogen-dominated atmosphere can be detected by the HST in case of multiple observations. Closer or more active M dwarfs (otherwise similar) will have stronger Ly$\alpha$ lines, therefore the SNR of the observations may also be better and merging multiple visits may not be necessary. Therefore, we conclude that hydrogen atmospheres may be detectable if the atmosphere of a planet is extended thus producing a very spacious exosphere. This is the case for close-in planets with sufficient hydrogen amount, such as hot Neptunes and Jupiters. We have reproduced the Ly$\alpha$ observations of GJ 436b assuming an escaping hydrogen-dominated atmosphere with a size significantly larger than the one of a terrestrial planet. Unlike terrestrial planets modeled in this article, GJ 436b produces a deep absorption signature in the Ly$\alpha$ line.

Our results indicate that further observations in the Ly$\alpha$ line can be more fruitful for expanded hydrogen-dominated planets orbiting active M dwarfs such as GJ 436b. Planets possessing secondary atmospheres do not produce a signature in the Ly$\alpha$ line strong enough to be detected by the Hubble Space Telescope at 10 pc or beyond, however, Ly$\alpha$ signatures of Earth-like planets with hydrogen-dominated atmospheres can be detected. It is possible that Ly$\alpha$ signatures of secondary atmospheres may be observable in the future using more advanced instrumentation, such as LUVOIR and WSO UV.

Acknowledgements.

We acknowledge the support by the Austria Science Fund (FWF) NFN project S116-N16 and the subprojects S11607-N16, S11606-N16 and S11604-N16. PO, HL, and NVE acknowledge support from the Austrian Science Fund (FWF) project P25256-N27 “Characterizing Stellar and Exoplanetary Environments via Modeling of Lyman-$\alpha$ Transit Observations of Hot Jupiters”. NVE also acknowledges support by the RFBR grant No 16-52-14006. MLK also acknowledges FWF projects I2939-N27 and the partial support by the Ministry of Education and Science of Russian federation (Grant No. RFMEFI61617X0084). IFS acknowleges support of Russian Science Foundation project 18-12-00080. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research was conducted using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North (HPC2N), Ume$\mathring{\rm a}$ University, Sweden. The authors are very thankful to Dr. David Ehrenreich for providing the Ly$\alpha$ spectra of GJ 436b, which were used in this article. We would also like to sincerely thank Dr. Vincent Bourrier and Baptiste Lavie for original processing of these spectra.