© ESO, 2016

1. Introduction

1.1. Atmospheric escape

Transit observations in the Lyman-α line of neutral hydrogen led to the detection of extended exospheres around the hot Jupiters HD 209458b (Vidal-Madjar et al. 2003, 2004) and HD 189733b (Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013), the warm Neptune GJ 436 b (Kulow et al. 2014; Ehrenreich et al. 2015), and the warm Jupiter 55 Cnc b (Ehrenreich et al. 2012). The intense stellar X-ray and extreme ultraviolet energy input at the base of a hydrogen-rich thermosphere has been shown to be responsible for the expansion of the upper atmospheric layers (e.g., Lammer et al. 2003; Lecavelier des Etangs et al. 2004; Koskinen et al. 2013a,b). Heavier species can be carried to high altitudes through collisions with the expanding flow of hydrogen, and several metals and ions were detected around these planets (Vidal-Madjar et al. 2004, 2013; Linsky et al. 2010; Ballester & Ben-Jaffel 2015; Ben-Jaffel & Ballester 2013; Fossati et al. 2010; Haswell et al. 2012), confirming that their atmospheres are in a state of hydrodynamic blow-off.

Transit observations at high resolution in the UV have also been used to probe the structure of these extended exospheres, revealing that they are shaped by interactions with the host star such as photoionization, radiation pressure, and stellar wind interactions (e.g., Holmström et al. 2008; Ekenbäck et al. 2010; Bourrier & Lecavelier des Etangs 2013; Ben-Jaffel & Ballester 2013; Bourrier et al. 2014; Kislyakova et al. 2014a; Guo & Ben-Jaffel 2016; Schneiter et al. 2016). Theoretical studies based on hydrodynamical simulations have also studied the processes that can affect the planetary outflow, such as charge-exchange reactions (e.g., Tremblin & Chiang 2013; Christie et al. 2016) or interactions with the planetary magnetic field (e.g., Khodachenko et al. 2015). The interpretation of absorption signatures with 3D numerical models of atmospheric escape allows studying not only the properties of the planetary outflow, but also obtaining direct constraints on the star, such as its X/EUV emission and wind properties. The detection of temporal variability in the exosphere of the hot Jupiter HD 189733b (Lecavelier des Etangs et al. 2012) additionally underlined the importance of multi-epoch observations to study the evolution of these properties.

While hot Jupiters have been the focus of most studies, they are subject to moderate escape rates that only weakly affect their long-term evolution. By contrast, theoretical studies (e.g., Lecavelier des Etangs et al. 2004; Lecavelier des Etangs 2007; Ehrenreich & Désert 2011; Owen & Jackson 2012; Lopez & Fortney 2013; Kislyakova et al. 2014b) and trends in the exoplanet population (e.g., Beaugé & Nesvorný 2013; Howard et al. 2012) show that lower-density planets like mini-Neptunes or super-Earths with a large volatile envelope may be most significantly affected by evaporation, leading in the more extreme cases to a massive erosion of the atmosphere and the formation of rocky remnant cores.

1.2. GJ 436 b

Located at the edge of the sub-Jupiter desert (Beaugé & Nesvorný 2013), the warm Neptune GJ 436 b (R_p = 4.2 R_Earth, P = 2.6 days, a = 0.0287 au; Butler et al. 2004; Gillon et al. 2007) is an ideal candidate to investigate the evaporation of low-mass gaseous planets. In contrast to other known evaporating planets orbiting G- and K-type stars, GJ 436 b is hosted by an M dwarf with moderate irradiation (M_⋆ = 0.45 M_⊙, R_⋆ = 0.44 R_⊙), which enables us to study a new regime of atmospheric escape and star-planet interactions. The brightness of the host star (V = 10.7) and its close proximity to Earth (d = 10.14 pc) makes GJ 436 a good target for transit observations in the Lyman-α line. Using HST/STIS, Kulow et al. (2014) identified a deep absorption signature from neutral hydrogen after the end of the optical transit, but their interpretation was misled by an inaccurate transit ephemeris and by the lack of an out-of-transit reference for the flux in the stellar Lyman-α line (Ehrenreich et al. 2015). Using two additional HST observations, Ehrenreich et al. (2015) revealed a deeper signature repeated over the three epochs of observations, which shows that GJ 436 b is surrounded by a giant coma of neutral hydrogen that is large enough to occult the stellar disk several hours before the optical transit, and that is trailed by a long cometary tail that could remain detectable for many hours after the optical transit (Fig. 1).

Fig. 1 Schematic representation of the neutral hydrogen cloud surrounding GJ 436 b, displaying the different regions of the exosphere. Colors distinguish between the two populations of hydrogen atoms with different origins that compose the cloud.

Following this detection, Bourrier et al. (2015a) studied the role played by stellar radiation pressure on the exosphere structure of GJ 436 b and its transmission spectrum. In contrast to hot evaporating planets, radiation pressure from its M-dwarf host star is too low to overcome stellar gravity and repel exospheric hydrogen atoms from the star. But it is still high enough to brake the gravitational deviation of the atoms toward the star, allowing their dispersion within a large volume around the planet. This effect is referred to as radiative braking. While radiative braking explains the size of the coma and the blueshifted velocity range of the observed absorption signatures up to about − 120 km s^-1 well, it does not account for the variations in the depth and duration of the absorption signal at the different phases of the transits. Furthermore, even though the overall signature of the GJ 436 b exosphere is very similar in the three different epochs, specific features to each visit cannot be explained by the radiation pressure, which is due to the Lyman-α line that was shown to be extremely stable over time (Ehrenreich et al. 2015; Bourrier et al. 2015a). Our goal in this paper is to investigate the coupled effects of stellar wind interactions and radiation pressure on the exosphere, using numerical simulations of the GJ 436 system with the EVaporating Exoplanet code (EVE). We also compare simulated spectra with the observations in each epoch to measure the values of planetary and stellar parameters that shape the exosphere.

This study is based on the three existing transit data sets of GJ 436 b observations (Table 2) taken in the H i Lyman-α line with the Space Telescope Imaging Spectrograph (STIS) instrument onboard the Hubble Space Telescope (HST). These datasets are described in Bourrier et al. (2015a). The code EVE was also described in this paper, and in Sect. 2 we summarize its main characteristics and detail the extension developed for stellar wind interactions. In Sect. 3 we investigate the differences between radiation pressure and stellar wind on the spatial and velocity structure of the exosphere and its spectral signature. EVE simulations are compared to the observed spectra in Sect. 4 to measure the properties of GJ 436 b environment at the different epochs of observations. These results are interpreted in Sect. 5. In Sect. 6 we discuss how the effects of radiation pressure and stellar wind interactions can be disentangled, and we conclude in Sect. 7.

2. Modeling the exosphere with the code EVE

EVE is a 3D numerical code developed to calculate the structure of an exoplanet upper atmosphere and its transmission spectra. The code was used in Ehrenreich et al. (2015) to perform a preliminary fit of the three combined observations of GJ 436 b neutral hydrogen exosphere. It was then used in Bourrier et al. (2015a) to investigate the influence of radiation pressure. In this paper, we implement a stellar wind extension to the code to study the effect of charge exchange on the exosphere and to measure the properties of GJ 436 b environment at the different observation epochs. The main physical parameters used for GJ 436 b and its host star are given in Table 1.

2.1. General description

A detailed description of the code EVE and the numerical settings used for GJ 436 b can be found in Bourrier et al. (2015a), and we summarize its main features here. The upper planetary atmosphere is divided into two different regimes that are joined at the mean altitude of the Roche lobe. The bottom layers of the atmosphere are described analytically, while Monte Carlo particle simulations are used to compute the dynamics of neutral hydrogen metaparticles in the upper atmospheric layers. Particles are subjected to the stellar and planetary gravities, the stellar radiation pressure, and the inertial force linked to the non-Galilean stellar reference frame. They are also affected by stellar photoionization and charge exchange with the stellar wind (Sect. 2.2). We use the term projected velocity to refer to the projection of a particle velocity on the star-Earth line of sight, and refer to its projection on the star-particle axis as radial velocity (i.e., the radial coordinate in the star reference frame). While we measure spectra as a function of projected velocity, radiation pressure is proportional to the flux in the intrinsic Lyman-α line and therefore it varies with the radial velocity of hydrogen atoms. The intrinsic stellar Lyman-α line can only be reconstructed by accounting for interstellar medium (ISM) absorption, and we used the lines reconstructed for Visits 2 and 3 in Bourrier et al. (2015a). As shown by these authors, the lack of out-of-transit observations in Visit 1 prevents reconstructing the intrinsic line, and we used the line obtained for Visit 2 as a proxy. The density and velocity structures of the gas in the exosphere depend on three free model parameters: the escape rate of neutral hydrogen (Ṁ_H⁰ in g s^-1, at the distance of the semi-major axis), the photoionization rate per atom (Γ_ion in s^-1, at the distance of the semi-major axis), and the velocity of the planetary wind at the Roche lobe (${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$, in km s^-1). We accounted for the effect of the orbital eccentricity (e = 0.16, Lanotte et al. 2014) that causes the escape rate and photoionization rate to vary in time with the inverse distance to the star squared. The effect of radiation pressure and photoionization on neutral hydrogen particles was calculated taking self-shielding within the exospheric cloud into account.

Constraints on the model parameters come from the direct comparison between the STIS observations and theoretical Lyman-α line spectra calculated at each time step with EVE at a resolution Δλ = 0.04 Å  corresponding to Δv = 10 km s^-1 (about half the resolution of the STIS spectra at 1215.67 Å). These theoretical spectra are affected by the planetary occultation, the exospheric absorption (taking the bulk motion, thermal, and natural broadening of the hydrogen gas into account), the interstellar medium (ISM) absorption, and STIS line spread function (LSF). The merit function for a given visit is the sum of the χ² yielded by the comparison of the observed spectra with the theoretical spectra averaged during the time window of each observation (see Table 2). Data obtained in time-tag mode was sliced into two exposures per HST orbit (Ehrenreich et al. 2015), yielding a good compromise between signal quality and time sampling of the simulations. The fits were calculated in the velocity ranges [− 200; − 25] and [40; 200]  km s^-1 for Visit 1 and [− 200; − 40] and [20; 200]  km s^-1 for Visits 2 and 3. The line core was excluded because of airglow contamination and ISM absorption, while high velocities in the wings of the line did not have enough signal. To summarize, the fits were performed on eight exposures per visit for a total of about 210 data points.

2.2. Extension: stellar wind interactions

A stellar wind proton may gain an electron from the interaction with a neutral hydrogen atom in the planetary exospheric outflow. With the usual assumption (Lindsay & Stebbings 2005) that charge transfer collisions result in little deflection of the interacting proton and no significant change in its kinetic energy, we consider that the population of neutralized protons keeps the velocity distribution of the stellar wind. With regard to observations in the Lyman-α line, charge exchange replaces the contribution of an exospheric neutral atom by that of a neutralized proton at the same position but with different velocity properties. We therefore handle this process in EVE as an impulsion given to the neutral atom undergoing charge exchange, so that afterward it moves with the velocity of the interacting proton. In that way, stellar wind protons need not be treated as an independent particle population because we only need to know the probability dP that a given neutral hydrogen atom is accelerated by a proton during a simulation time step dt: $\mathrm{d}\mathit{P}\mathrm{=}\mathrm{1}\mathrm{-}\exp \mathrm{[}\mathrm{-}{\mathit{\sigma }}_{\mathrm{H}{\mathrm{H}}^{\mathrm{+}}}\mathrm{(}\mathrm{\Delta }\mathit{V}\mathrm{)}\hspace{0.17em}\mathrm{\Delta }\mathit{V}\hspace{0.17em}{\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}}\hspace{0.17em}\mathrm{d}\mathit{t}\mathrm{]}\mathit{,}$(1)with n_H⁺ the stellar wind proton density in the vicinity of the hydrogen atom and ΔV = | | V_H − V_H⁺ ∥ the relative velocity between the neutral atom and the interacting proton. Compared to Bourrier et al. (2015a), we reduced the time step dt to 2.6 min to better account for the fast dynamics of the stellar wind. The cross section of the interaction σ_HH⁺ is energy dependent. For relative velocities in our simulations lower than 1000 km s^-1 (energy lower than 5.2 keV), the formula from Lindsay & Stebbings (2005) can be approximated to ${\mathit{\sigma }}_{{\mathrm{HH}}^{\mathrm{+}}}\mathrm{=}{\mathrm{10}}^{-20}\mathrm{(}\mathrm{10.61}\mathrm{-}\mathrm{1.062}\hspace{0.17em}\ln \mathrm{(}\mathrm{\Delta }\mathit{V}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(2)with σ_HH⁺ in m², ΔV in km s^-1. The velocity of the interacting proton is taken from the Maxwellian speed distribution of the stellar wind, defined analytically through its radial bulk velocity ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ and kinetic temperature ${\mathit{T}}_{\mathrm{wind}}^{\mathrm{st}}$, which are assumed to vary little over the spatial extension of the exosphere. We note that this kinetic temperature corresponds to the Maxwellian velocity dispersion of the proton population ${\mathit{v}}_{\mathrm{therm}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$=$\sqrt{\mathit{k}\hspace{0.17em}{\mathit{T}}_{\mathrm{wind}}^{\mathrm{st}}\mathit{/}{\mathit{m}}_{{\mathrm{H}}^{\mathrm{+}}}}$ around the bulk motion of the stellar wind. This bulk motion is associated with a different temperature, generally in the order of millions of kelvins, in the frame of the Parker theory (e.g., Vidotto et al. 2010).

Outside of the exosphere, we assumed that the stellar wind density decreases as a function of the distance from the star r according to a quadratic law. Within the exosphere, hydrogen atoms shield each other from the stellar protons in the same way as from stellar photons. The proton density n_H⁺ at the distance r from the star decreases with the penetration depth Λ into the atmosphere as $\begin{array}{ccc}& & {\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}\hspace{0.17em}{\left(\frac{{\mathit{a}}_{\mathrm{p}}}{\mathit{r}}\right)}^{\mathrm{2}}\hspace{0.17em}\exp \mathrm{[}\mathrm{-}{\mathit{\tau }}_{\mathrm{H}{\mathrm{H}}^{\mathrm{+}}}\mathrm{(}\mathrm{\Lambda }\mathrm{\left)}\mathrm{\right]}\\ & & \end{array}$(3)with ${\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}$ the value of the proton density at the distance of the planet semi-major axis a_p, and n_H the neutral hydrogen density at an intermediate penetration depth μ. In the upper layers of the exosphere, each metaparticle contributes to the optical depth τ_HH⁺ at its own velocity. To calculate the contribution of the lower atmosphere, we used the equations from Fahr & Bzowski (2004) to integrate ΔV(μ) over the velocity distribution of the neutral hydrogen gas, assumed to be a Maxwellian centered on the planet orbital velocity. We simplified these calculations by using the average velocity of the stellar wind ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$.

Our description for the stellar wind adds three more free parameters to the model: the bulk velocity (${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ in km s^-1) and kinetic dispersion (${\mathit{v}}_{\mathrm{therm}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ in km s^-1) of the proton distribution, which are representative of the wind conditions at the location of the planet, and the proton density (${\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}$ in cm^-3, given at the semi-major axis of the planet).

3. Structure of the exosphere

3.1. Radiative braking and stellar wind interactions

The radiation pressure from the M dwarf GJ 436 has a strong influence on the exosphere of its warm-Neptune companion and reproduces the velocity range of its absorption signature well (Bourrier et al. 2015a). However, radiative braking alone does not explain the variations of the absorption depth observed at the different phases of the transit well (Sect. 1.2). While it allows the formation of a coma that is large enough to occult about half of the stellar disk at the center of the optical transit, this coma also extends too far ahead of the planet and produces a deeper and earlier ingress than observed (Bourrier et al. 2015a). Furthermore, each observation epoch shows specific features that cannot be explained by a stable radiation pressure (see Fig. 2). The flux in the blue wing of the Lyman-α line varies more smoothly over time during Visit 1, with a longer transit duration than other visits, but a lower absorption depth at the center of the transit. In contrast, Visit 2 shows sharper flux variations at the ingress and egress and a shorter, deeper transit. Finally, ingress starts the latest in Visit 3 with a dramatic increase in absorption depth, which surprisingly remains at about the same level during the transit and post-transit phases. This feature is at odds with the gradual decrease in absorption depth caused by stellar photoionization and the dilution of the gas subjected to radiative braking (Bourrier et al. 2015a).

An additional mechanism is therefore needed 1) to reduce the size of the coma ahead of the planet in all epochs and 2) to explain variations in the dynamics and the geometry of the exosphere between the different epochs. In this section, we show that by abrading the hydrogen cloud formed by planetary escape and by creating a secondary population of neutral hydrogen atoms, stellar wind interactions might be able to explain these two points. In Sect. 4 we compare observations and EVE simulations with a stable radiation pressure and a variable stellar wind and derive the corresponding properties of the stellar wind and planetary environment at the different epochs.

Fig. 2 GJ 436 spectra during Visits 1 (left column), 2 (central column), and 3 (right column). Spectra in the same row were measured at about the same phase of the planet orbital position. We show the spectra gathered over the full HST orbits. The shaded gray area corresponds to the range affected by ISM absorption and geocoronal emission, which was excluded from the fits. Black solid spectra are the reconstructed out-of-transit stellar line profiles. Dashed black spectra correspond to the best-fit theoretical spectra for Visits 2 (reported for comparison in Visit 1 epoch) and 3. Note that the transit of the exosphere has little effect on the red wing, where the best fits nearly overlap with the out-of-transit spectra. In the bottom part of the plot, light curves for each epoch show the observed and theoretical flux integrated between − 120 and − 40 km s-1. The solid black line shows the optical transit.

3.1.1. Abrasion

Through charge exchange, stellar wind protons ionize the neutral hydrogen exosphere arising from the planetary escape. Hereafter we refer to this phenomenon as abrasion to distinguish from the effect of stellar photoionization. The regions of the coma ahead of the planet and facing the star are strongly affected by abrasion, while self-shielding from the protons can protect the farthest regions opposite the star (Fig. 1). Nonetheless, our best-fit simulations (see Sect. 4.2) show that hydrogen densities are low enough around GJ 436 b for the stellar wind to reach far within the core of the coma, interacting with many hydrogen atoms even before they can move away from the planet into the outer regions of the exosphere. Stellar wind abrasion thus also indirectly effects the cometary tail and the front of the coma, which are fueled by the expanding inner regions of the exosphere (see Figs. 1 and 3). This leads to a steeper decrease in absorption depth before and after the optical transit (Fig. 4).

Fig. 3 Views of GJ 436 b exosphere (dots) within the orbital plane, seen from the top, at the time of the optical transit. The two panels correspond to two simulations performed with the same escape rate, planetary wind velocity, and photoionization rate, but without stellar wind interactions (left panel) or including them (right panel, corresponding to the best fit for Visit 2). The planet is represented by a small black disk, to scale, between the dashed black lines limiting the LOS toward the stellar disk. Arrows display the average velocity field of the neutral hydrogen atoms in the stellar rest frame, with particles colored as a function of the time they escaped the atmosphere (for planetary neutrals) or were created through charge exchange (for neutralized protons). With no stellar wind interactions the dynamics of the escaping gas is initially dominated by the planetary orbital velocity and later constrained by radiative braking. With charge exchange most protons are neutralized in the dense inner coma, abrading at the source planetary neutrals that would have fed the outer regions of the exosphere, but also creating a compact population of neutrals dominated by the radial bulk velocity of the stellar wind. It takes much more time for stellar gravity to eventually disperse this population of neutralized protons.

3.1.2. Neutralized protons

In addition to abrading neutral hydrogen atoms in the planetary exosphere, another effect of charge-exchange interactions is to neutralize protons in the stellar wind. This leads to the formation of a secondary population of neutral hydrogen in the exosphere (hereafter, neutralized protons), which differs from the population of neutral hydrogen atoms that escaped from the upper atmosphere of the planet (hereafter, planetary neutrals). This population of neutralized protons has two important characteristics. First, it arises primarily from the inner coma of GJ 436 b where most protons are neutralized because of the higher density of planetary neutrals close to the planet (Fig. 3, right panel). Then, planetary neutrals escaping the atmosphere initially have a low radial velocity that naturally increases as they move away from the planet orbit¹. In contrast, neutralized protons continue to move with the velocity distribution of the stellar wind, which is dominated by a high radial bulk velocity. Consequently, while it takes several hours for planetary neutrals to move into the outer regions of the exosphere, neutralized protons swiftly move away from the planet and can go farther before they are photoionized. Because of these two features, the population of neutralized protons is shaped into a compact cloud originating from the core of the coma and extending into a long comet-like tail that moves with the persistent velocity distribution of the stellar wind (Fig. 3, right panel). This tail spreads more slowly than the planetary one, both dynamically and geometrically (Fig. 3, see left and right panel). For several hours after the optical transit the population of neutralized protons thus produces a secondary absorption profile with a fairly stable spectral range and slowly decreasing depth (Fig. 4). Eventually, this profile spreads and shifts toward more positive radial velocities as stellar gravity slows down atoms in the tail and creates a strong velocity gradient between its outer C-shaped region and the inner C-shaped region closer to the star (Fig. 3, right panel). We caution that for the absorption signature of neutralized protons to become directly observable, their bulk velocity must be lower than about 200 km s^-1. With higher absolute velocity there is not enough flux in the Lyman-α line of GJ 436 (see Fig. 2), and only the abrading effect of the wind on the exosphere can be detected.

Fig. 4 Views of GJ 436 b exosphere (lower panels) and its theoretical absorption profile (upper panels). We caution that absorption in the shaded range ([− 40; 20] km s-1) cannot be observed from Earth and that the STIS LSF would spread the absorption signal from the exosphere over more data points. Time increases from left to right and is indicated in each plot. The dashed black lines limit the LOS toward the stellar disk, while the dotted line indicates the star-planet axis. Black arrows show the velocity field of the gas in the stellar rest frame, with H0 atoms in the orbital plane colored as a function of their LOS-projected velocity. Red lines in the upper panels show the absorption profiles from the exosphere simulated with Visit 3 best-fit parameters. The population of neutralized protons arising from the inner coma can be identified by its radial velocities and stable high-velocity absorption peak, while planetary neutrals remain present at the front of the coma and in the outer regions of the exosphere opposite the star with dynamics dominated by radiative braking. The blue line profiles in the upper panels correspond to the same simulation with no stellar wind, illustrating how charge exchange reduces the size of the exosphere before and after the optical transit but only weakly affects the absorption profile at the center of the transit. It also shows that neutralized protons produce a stable high-velocity absorption peak.

3.2. Influence of the physical parameters

We describe here the first-order influence of the escape rate and velocity of the planetary outflow and of the stellar photoionization rate on the structure of the exosphere and its transmission spectrum when it is subjected to stellar wind interactions. We describe the role of each parameter of the modelled stellar wind in more detail (previous descriptions in the literature can be found in Holmström et al. 2008; Ekenbäck et al. 2010; Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014a).

• Planetary wind velocity ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$: this parameter corresponds to the initialupward bulk velocity of the atoms escaping the upper atmosphereand influences the spatial and spectral dispersion of the gas in theexosphere. With higher escape velocities, the wider distributionof Doppler velocities for the planetary neutrals spreads their ab-sorption signature over a broader spectral range, but the velocitydistribution of the neutralized protons remains that of the stellarwind (Sect. 2.2) and their absorption profile cov-ers the same spectral range. Nonetheless, the effective variationin absorption depth with time and wavelength is complex and de-pends on the column densities and stellar disk area that is occultedby neutral hydrogen coming from both populations. For exam-ple, an increase in ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$ expands the coma of planetary neutrals and increases its dilution. Stellar wind protons are thus neutralized within a larger volume, but with lower local densities.

• Photoionization rate Γ_ion: photoionization leads to a general decrease in density and absorption from all regions of the exosphere, albeit with differences caused by self-shielding and the influence of stellar wind interactions. The structure of the exosphere is less affected by the photoionization of neutralized protons than planetary neutrals. The latter move through the exosphere on longer timescales, making their population more vulnerable to photoionization, which additionally removes the neutral planetary material required for charge exchange with the stellar wind (Sect. 3.1.2).

• Escape rate Ṁ_H⁰: higher values for the escape rate increase the densities of planetary neutrals and neutralized protons in the exosphere, which deepens the absorption profiles of both populations. Different regions of the exosphere have different velocity distributions and contribute to the absorption in different spectral ranges (Sect. 3.1.2). However, self-shielding remains low enough for GJ 436 b that variations in the escape rate propagate roughly uniformly throughout the whole exosphere and do not affect the velocity distribution of the gas and the spectral range of its absorption profile. We caution that variations related to Ṁ_H⁰ would be different for saturation regimes of the stellar wind, in which either all escaping planetary neutrals would interact with the incoming protons (as for HD 189733b, Bourrier & Lecavelier des Etangs 2013), or all protons in the stellar wind would be neutralized by the planetary exosphere. Simulations show, however, that GJ 436 b is not observed in such regimes (Sect. 4.2).

• Stellar proton density ${\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}$: this parameter influences the number of planetary neutrals abraded by charge exchange and consequently the number of neutralized stellar wind protons. A higher stellar wind density increases the abrasion of the exosphere, which reduces the size of the exosphere envelop and the absorption depth before and after the optical transit (Sect. 3.1.1). This can be partly compensated for by the neutralized protons that contribute to the observed Lyman-α absorption in the spectral range that corresponds to the stellar wind velocity distribution. This additional absorption is mainly visible during and after the optical transit because protons neutralized ahead of the planet move away into the tail much more quickly than the planetary neutrals (Sect. 3.1.2).

• Stellar proton kinetic temperature ${\mathit{T}}_{\mathrm{wind}}^{\mathrm{st}}$: the kinetic temperature of the stellar wind controls the spread of the protons velocity distribution and thus the breadth of the neutralized protons absorption profile population. A higher temperature decreases the absorption depth in the core of the profile in favor of its wings. The overall effect on the exosphere is stronger after the optical transit because a higher dispersion of the protons neutralized in the coma feeds back the spatial and spectral dispersion of the gas in the cometary tail. The value of ${\mathit{T}}_{\mathrm{wind}}^{\mathrm{st}}$ has very little influence on the abrasion of the exosphere because the cross section for charge exchange varies slowly with the relative velocities at play between protons and planetary neutrals (Sect. 2.2), and their probability of interaction therefore does not depend strongly on the thermal dispersion of the proton velocity.

• Stellar proton velocity ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$: the velocity distribution of the stellar wind is dominated by its radial bulk velocity. Therefore the projected velocity of the neutralized protons is close to their absolute velocity in the stellar rest frame when they are transiting the star, and ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ determines the center of their absorption profile in the spectra. By contrast, the velocity distribution of the planetary neutrals is strongly dependent on the planet orbital velocity, and their absorption profile is closer to the core of the Lyman-α line (Bourrier et al. 2015a). Regarding the abrasion of the exosphere, variations in the bulk velocity affect both the relative velocities between protons and planetary neutrals and their interaction cross-section (Sect. 2.2). An increase in ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ tends to slightly increase the probability for a planetary neutral to undergo charge exchange.

4. Measuring the properties of GJ 436 b environment

4.1. Stellar XEUV emission and photoionization rate

We searched for prior constraints to place on the EVE parameters used to interpret the Lyman-α line observations of GJ 436 b. It is not possible to use the planetary mass-loss properties derived by Kulow et al. (2014), which are biased by their interpretation of Visit 1 (Sect. 1.2). Ehrenreich et al. (2015) did not account for possible temporal variability between each visit and made a strong assumption on the presence of a sharp density transition in the outer regions of the exosphere facing the star. Bourrier et al. (2015a) obtained rough estimates for the planetary outflow velocity, but with the assumption that radiation pressure alone acts on the exosphere. Regarding the host star, no measurements of the stellar wind are available. However, Bourrier et al. (2015a) reconstructed the intrinsic stellar Lyman-α line for Visits 2 and 3, and this can be used to estimate the EUV spectrum of GJ 436 and the resulting photoionization rate per neutral hydrogen atom at 1 au from the star: $\begin{array}{c}\mathrm{􏽥}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}\mathrm{=}{\mathrm{\int }}^{\mathrm{911.8}\hspace{0.17em}\AA }\frac{{\mathit{F}}_{\mathrm{EUV}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{ion}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}}{\mathit{hc}}\hspace{0.17em}\mathit{\lambda }\hspace{0.17em}\mathrm{d}\mathit{\lambda ,}$(4)with $\begin{array}{c}{\mathrm{􏽥}}^{}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}$ in s^-1, F_EUV(λ) the stellar flux at 1 au (in erg s^-1 cm^-2 Å^-1) and σ_ion the cross section for photoionization (in cm²). The integration was performed up to the ionization threshold at 911.8 Å. The cross section for hydrogen photoionization is wavelength dependent, and we used the expression from Verner et al. (1996) and Bzowski et al. (2013): $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{ion}}& \mathrm{=}& \mathrm{6.538}\mathrm{\times }{\mathrm{10}}^{-32}\hspace{0.17em}{\left(\frac{\mathrm{29.62}}{\sqrt{\mathit{\lambda }}}\mathrm{+}\mathrm{1}\right)}^{-2.963}\\ & & \mathrm{\times }\mathrm{(}\mathit{\lambda }\mathrm{-}\mathrm{28846.9}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\lambda }}^{\mathrm{2.0185}}\mathit{,}\end{array}$(5)with σ_ion in cm² and wavelengths in Å. The X-ray emission of GJ 436 from 5 to 100 Å  was measured with Chandra at the epoch of Visit 2 (Ehrenreich et al. 2015). The EUV emission from 100 to 912 Å  is mostly absorbed by the ISM and must be estimated through indirect methods. We used the scaling relations of Linsky et al. (2014)², which are based on the integrated intrinsic Lyman-α flux at 1 au from GJ 436 that we measured to be 0.90$\stackrel{\mathrm{+}\mathrm{0.20}}{{-0.14}_{}}$ erg s^-1 cm^-2 for Visit 2 and 0.88$\stackrel{\mathrm{+}\mathrm{0.30}}{{-0.21}_{}}$ erg s^-1 cm^-2 for Visit 3. The uncertainties on these values were obtained by varying the peak flux of the Lyman-α line constrained by Bourrier et al. (2015a) within its 1σ error bars. We note that France et al. (2013) performed a similar estimation using the intrinsic Lyman-α line of GJ 436 observed in 2010 (Ehrenreich et al. 2011). With a maximum uncertainty of 30%, the total flux they obtained ranges from 1.1 to 2.0 erg s^-1 cm^-2 at 1 au from the star. This is marginally higher (~1.2σ) than our estimates, possibly because of long-term variations of the Lyman-α line flux between the 2−3.5 yr that separate this earlier observation from Visits 2 and 3, or because of short-term temporal variability through impulsive flares in chromospheric and transition region emission lines of the host star at this epoch (France et al. 2013).

Our results for the stellar flux and luminosity in Visits 2 and 3 are shown in Table 3 for complementary wavelength bands of the XEUV domain. Uncertainties on the integrated Lyman-α flux were propagated on the flux in each domain to estimate 1σ error bars on $\begin{array}{c}{\mathrm{􏽥}}^{}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}$. We neglected uncertainties on the X-ray emission because it is about four times lower than the EUV emission and σ_ion steeply decreases with $\mathit{\lambda }{}^{\begin{array}{c}\begin{array}{c}\mathit{<}\end{array}\\ \mathrm{~}\end{array}}\mathrm{100}$Å, making the contribution of the X-ray flux to the photoionization rate negligible. We found photoionization rates at 1 au from the star of 1.9$\stackrel{\mathrm{+}\mathrm{0.7}}{{}_{-0.4}}\mathrm{\times }{\mathrm{10}}^{-8}$ s^-1 for Visit 2 and 1.8$\stackrel{\mathrm{+}\mathrm{1.1}}{{}_{-0.6}}\mathrm{\times }{\mathrm{10}}^{-8}$ s^-1 for Visit 3. This corresponds to 2.3$\stackrel{\mathrm{+}\mathrm{0.8}}{{}_{-0.5}}\mathrm{\times }{\mathrm{10}}^{-5}$ s^-1 and 2.2$\stackrel{\mathrm{+}\mathrm{1.3}}{{}_{-0.7}}\mathrm{\times }{\mathrm{10}}^{-5}$ s^-1 at the distance of the semi-major axis. The values for both visits are remarkably similar, as expected from the very stable Lyman-α line profiles between the two epochs (Ehrenreich et al. 2015; Bourrier et al. 2015a). Because of this stability and the short six-month interval between Visits 1 and 2, we assumed these photoionization rates to be a good estimate for Visit 1.

4.2. Constraints from the observed spectra

Using the settings described in Sect. 2, we compared spectra from EVE simulations with Lyman-α observations to measure the properties of GJ 436 b environment that best explain the observations. The spectra in each visit were fit independently to strengthen the reliability of the estimations and test the scenario with different data sets. Because of the runtime for a single simulation (from ~5 h to several days), it is not practical to explore the parameter space using MCMC algorithms. Instead, we computed the χ² of the fits on a grid scanning all possible values for the six model parameters (see Sect. 2): the escape rate of neutral hydrogen Ṁ_H⁰, the planetary outflow velocity ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$, the photoionization rate Γ_ion, and the stellar wind properties (bulk velocity ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$, kinetic dispersion ${\mathit{v}}_{\mathrm{therm}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$, and density ${\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}$ of the protons distribution). When the absolute minimum χ² and corresponding best values for the parameters were obtained, we calculated their error bars from an analysis of χ² variations. A given parameter was pegged at various trial values, and for each trial value we searched for the minimum χ² with the five other parameters that were allowed to vary freely. The 1σ error bar for the pegged parameter was obtained when its value yields a χ² increase of 1 from the absolute minimum (see, e.g., Hébrard et al. 2002). We used the independent estimates of the photoionization rate $\begin{array}{c}{}^{\mathrm{􏽥}}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}\mathrm{\pm }{\mathit{\sigma }}_{\begin{array}{c}{}^{\mathrm{􏽥}}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}}$ derived in Sect. 4.1 as constraints on the model value Γ_ion, adding to the χ² the term $\mathrm{\left(}\mathrm{\right(}{\mathrm{\Gamma }}_{\mathrm{ion}}\mathrm{-}\begin{array}{c}{}^{\mathrm{􏽥}}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}\mathrm{)}\mathit{/}{\mathit{\sigma }}_{\begin{array}{c}{}^{\mathrm{􏽥}}\\ {\mathrm{\Gamma }}_{\mathrm{ion}}\end{array}}{\mathrm{)}}^{\mathrm{2}}$.

In a first scenario (Sect. 4.2.1), we explored moderate values for the stellar wind bulk velocities in the range of the Lyman-α line Doppler width (≲200 km s^-1) to allow for both abrasion and proton neutralization in the exosphere. Abrasion alone from a high-velocity stellar wind (≳200 km s^-1) is discussed as a second scenario in Sect. 4.2.2.

4.2.1. Stellar wind abrasion and proton neutralization

Best-fit values for Visits 2 and 3 are given in Table 4, along with their 1σ uncertainties. They provide a good fit to the data, with χ² of 203 and 182 for 200 and 198 degrees of freedom (d.o.f.), respectively. Observations in both epochs require significant interactions between the exosphere and the stellar wind, with proton densities in the order of 10³ cm^-3 at the location of the planet. These interactions abrade the coma and produce a transit ingress starting ~3 h before the optical transit center, consistent with the observations (Fig. 3). Without these interactions, the ingress would start too early, about 5 h before the optical transit center (Bourrier et al. 2015a). The proton density is higher in Visit 3 than in Visit 2, which is consistent with a stronger abrasion causing the sharper ingress observed during Visit 3 (Sect. 3.1.1) and with a more abundant population of neutralized protons in the exosphere producing a stable absorption signature at later orbital phases (Sect. 3.1.2). Except for variations in the proton density and a marginally lower atmospheric escape velocity in Visit 2, the planetary outflow and stellar properties inferred for the GJ 436 system are very similar in both epochs, with all parameters consistent at the 1σ level. The consistency between these parameters, derived from independent datasets obtained at two different epochs, gives credence to the present scenario and shows the stability of the planetary mass loss over time. We note that removing the prior constraint on the photoionization rate has little influence on the quality of the fits over a broad range of values because it leads to changes in the local density of planetary neutrals that can be compensated for by variations in the escape rate and/or in the proton density.

We found that we were unable to perform a χ² analysis for Visit 1, with no clear minimum and spurious χ² variations that did not allow us to constrain the parameter values at this epoch. Using the lines reconstructed for Visits 2 and 3 as reference for Visit 1 may have biased the evaluation of radiation pressure and the calculation of the theoretical transmission spectra for this epoch. We also investigated whether these problems might have been caused by the flux variations in the red wing of the Lyman-α line (Fig. 2), since absorption at positive velocities cannot not be explained by stellar wind interactions. But limiting the fit to the blue wing of the line did not improve the χ² analysis. It is possible that these variations in the shape of the line are caused by active Lyman-α regions at the surface of the stellar disk (Llama & Shkolnik 2016), although this is made unlikely by the very localized spectral ranges of these variations and the stability of the Lyman-α line over time (Fig. 2). Finally, we compared the best fits to Visits 2 and 3 (Table 4) with Visit 1 observations and obtained χ² values of 204 and 212 for 104 points in the blue wing (445 and 447 for 208 points in both wings). Although Visit 1 shares similarities with the other epochs (see Ehrenreich et al. 2015), this hints at more drastic differences in the physical conditions of GJ 436 b exosphere during the first epoch.

4.2.2. High-velocity stellar wind

Because the Lyman-α line extends between ~±200 km s^-1, interactions between a stellar wind moving faster than this limit and the planetary exosphere will produce a different observational signature than in the first scenario (Sect. 4.2.1). In that case, neutralized protons created by charge exchange have radial velocities that are too high to be visible in transmission in the Lyman-α line, and only their abrasion of planetary neutrals affects the observed signature (Sects. 3.1.2 and 6). We investigated whether the observations of GJ 436 b might be explained by interactions between the exosphere and a fast stellar wind.

To explore this scenario, ${\mathit{V}}_{\mathrm{bulk}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}$ was arbitrarily fixed to a value higher than 200 km s^-1. With the planetary neutral velocities below ~120 km s^-1 (Bourrier et al. 2015a), any value of the stellar wind bulk velocity beyond 200 km s^-1 will produce the same results (because it will only influence the probability for a planetary neutral to be abraded; see Eq. (1)), and it was fixed to 350 km s^-1. The kinetic dispersion of the proton distribution has little influence on the charge-exchange probability as well and was fixed to ${\mathit{v}}_{\mathrm{therm}\mathrm{-}\mathrm{wind}}^{\mathrm{st}}\mathrm{=}\mathrm{20}$ km s^-1. We then searched for the best fit by varying other free parameters. We were unable to find a good fit for Visit 3, as expected from the similar absorption signatures observed during both the transit and post-transit phases, which can only be produced by low-velocity neutralized protons in the exospheric tail (thus visible in transmission in the stellar line; Sect. 3.1.2). The best fit for Visit 2 was obtained for the following values: Ṁ_H⁰ = 1.6× 10⁹ g s^-1, Γ_ion = 1.2× 10^-5 s^-1, ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}\mathrm{=}\mathrm{45}$ km s^-1, and ${\mathit{n}}_{\mathrm{wind}}^{\mathrm{st}}\mathrm{=}\mathrm{3.6}$× 10³ cm^-3, yielding a χ² of 209.

These values are similar to those obtained in the low wind-speed scenario (Sect. 4.2.1), but the quality of the fit is lower for Visit 2 data, and a high-velocity wind scenario does not match Visit 3 data. This is in stark contrast with the low-velocity scenario that could explain both Visits 2 and 3 with very similar stellar wind and mass-loss properties, consistent with the stability of GJ 436 in the Lyman-α line and in the X-rays (Ehrenreich et al. 2015; Bourrier et al. 2015a). Hereafter, we therefore assume that the most likely scenario for the observations of GJ 436 b is a low-velocity stellar wind allowing for both abrasion and proton neutralization in the planetary exosphere.

5. Interpretation of the results

5.1. Stellar winds and planetary exospheres

In the models matching the observations, GJ 436 b is not in saturation regimes where either all escaping planetary neutrals interact with the stellar wind, or all stellar wind protons crossing the exosphere are neutralized by charge exchange. The observations of the warm Neptune constrain the exosphere to be composed of both planetary neutrals and neutralized protons. Planetary neutrals are mainly present in the coma close to the planet, while neutralized protons dominate the cometary tail (see Fig. 3 and Sect. 3). This situation is different for the hot Jupiter HD 189733 b, where the isolated absorption signature observed at very high velocity in the blue wing of the Lyman-α line was representative of a neutralized proton population alone (Lecavelier des Etangs et al. 2012). Interactions of the K-type star HD 189733 with its hot-Jupiter companion (a = 0.031 au or 8.8 R_∗) are stronger than for the GJ436 system, with stellar wind velocities in excess of ~200 km s^-1 and a proton density in the range 4 × 10³−5 × 10⁷ cm^-3 (Bourrier & Lecavelier des Etangs 2013; Ben-Jaffel & Ballester 2013). By comparison, the warm Neptune GJ 436 b (a = 0.0287 au or 14.1 R_∗) is subject to wind velocities in the order of 70−90 km s^-1 and proton densities in the range 1 × 10³−5 × 10³ cm^-3 from its M dwarf host star. The properties of HD 189733 stellar wind at the orbit the planet are most likely variable over time (see Bourrier & Lecavelier des Etangs 2013; Llama et al. 2013; Cauley et al. 2015; Fares et al., in prep.), and this also seems to be the case for GJ 436 with a proton density in Visit 2 lower than in Visit 3. This may be linked to temporal variability or non-homogeneities in the stellar magnetic field at the orbital distance of the planet, which could be further investigated using spectropolarimetric observations of the M dwarf.

5.2. Planetary mass loss

5.2.1. Energy-limited escape rate and ionization fraction

In the energy-limited regime, a fraction of the stellar X/EUV energy received by the upper planetary atmosphere is converted into mass loss and compensates for the gravitational potential energy required by the gas to escape. Based on results of Owen & Alvarez (2016) (see their Figs. 1 and 4), we infer that GJ 436 b is located in this regime, since the warm Neptune has a mass of 1.4 × 10²⁹ g (Butler et al. 2004) and radius of 2.5 × 10⁹ cm (Knutson et al. 2011), and orbits a star with EUV emission in the order of 10² erg s^-1 cm^-2 at the semi-major axis³. The escape rate of neutral hydrogen can therefore be expressed as (e.g., Bourrier et al. 2015b) $\mathit{Ṁ}\mathrm{tot}\mathrm{=}\mathit{\eta }\hspace{0.17em}\frac{\mathrm{3}\hspace{0.17em}{\mathit{F}}_{\mathrm{X}\mathit{/}\mathrm{EUV}}\mathrm{\left(}\mathrm{1}\mathrm{au}\mathrm{\right)}}{\mathrm{4}\hspace{0.17em}\mathit{G}\hspace{0.17em}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}\hspace{0.17em}\overline{\mathit{\rho }}\hspace{0.17em}{\mathit{K}}_{\mathrm{tide}}}\mathit{,}$(6)with $\overline{\mathit{\rho }}$ the mean density of the planet, K_tide a correction factor accounting for the contribution of tidal forces to the potential energy (Erkaev et al. 2007), and η the heating efficiency, which the most recent theoretical estimations estimate at between 10 and 20% (e.g., Lammer et al. 2013; Shematovich et al. 2014; Owen & Alvarez 2016). The total X/EUV flux per unit area at 1 au from GJ 436, measured from 0.5 nm to the Lyman-α line (included), is F_X/EUV(1 au) = 2.3 ± 0.5 erg s^-1 cm^-2 (Table 3). The corresponding total mass-loss rate from GJ 436 b atmosphere is Ṁ^tot = ηṀ^100%, with Ṁ^100% = 2.2 ± 0.6 × 10¹⁰ g/s. Defining Ṁ_H⁰ = f_H⁰Ṁ^tot, where f_H⁰ is the neutral fraction of hydrogen in the upper thermosphere of the planet and Ṁ_H⁰ the escape rate of neutral hydrogen obtained with EVE, we found that ηf_H⁰ should be 1.2 ± 0.5 × 10^-2.

5.2.2. Magnetically driven outflow

Our simulations show that the atmospheric escape velocity must be in the range 40−70 km s^-1 to explain the size of the exosphere (Table 4). This outflow velocity is faster than predicted for hot Jupiters (~1−10 km s^-1; see, e.g., Murray-Clay et al. 2009; Koskinen et al. 2013a). To solve this discrepancy, we investigated the possibility that the outflow might arise from MHD waves in the upper atmosphere of GJ 436 b. In the presence of a magnetic field, the turbulence of the gas can excite MHD waves that dissipate in the upper atmosphere and drive the outflow (Suzuki & Inutsuka 2005, 2006). We used the model described in Tanaka et al. (2014) to calculate the 1D atmospheric structure of GJ 436 b for different values of the velocity dispersion at the surface. The velocity that we observed at the Roche lobe can be obtained with dispersions of about [6−17]% of the sound speed (2.4 km s^-1 at the planet surface, assuming a temperature of 800 K). Under these conditions, the surface velocity dispersion is high enough to drive a fast planetary wind, but low enough that the increasing density of the outflow does not hinder its acceleration (Fig. 5). Because of the nonlinear dissipation of the MHD wave energy in the upper atmosphere, the structure of the outflow is time variable (Tanaka et al. 2015), and the values measured here for ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$ are reached when the outflow is close to its maximum speeds. Theoretical studies have shown that Alfvenic waves can be severely damped at low altitudes in the atmosphere, confining the heating and enhanced mass loss in specific regions such as the magnetic poles (Trammell et al. 2011; Tanaka et al. 2015; Khodachenko et al. 2015). For GJ 436 b the large measured transit depth requires a giant coma that surrounds the planet on all sides, which in turn requires gas to escape from equatorial as well as polar regions. We conclude that the planetary outflow is unlikely to be confined to a specific region of the upper atmosphere but can be magnetically driven by MHD waves that dissipate in the upper atmosphere, which would explain the high velocity observed at the base of GJ 436 b exosphere. A self-consistent 3D model including both a MHD-described thermosphere and particle-described exosphere is nonetheless required to fully explore this scenario.

Moreover, when the time-variable outflow is close to its maximum speed, simulations of the magnetically driven outflow yield a total atmospheric mass-loss rate Ṁ^mag in the range [8 × 10⁷, 5 × 10¹⁰] g/s, which agrees remarkably well with the observations. The comparison between this calculated total escape rate Ṁ^mag, the total escape rate in the energy limited regime Ṁ^tot, and the observed escape rate of neutral hydrogen Ṁ_H⁰ allows us to further constrain the heating efficiency η and neutral fraction of hydrogen f_H⁰. Using the relations Ṁ^mag = Ṁ^tot = ηṀ^100% = Ṁ_H⁰/f_H⁰, the estimated Ṁ^100% (Sect. 5.2.1) and the measured H⁰ (Table 4), we find lower limits of ~0.5% for both η and f_H⁰. This is consistent with the value of about 1% for η estimated by Ehrenreich et al. (2015).

In conclusion, a magnetically driven outflow provides a consistent scenario for the observed escaping planetary wind, in which the measured values for the wind velocity and the mass-loss rate (and the related heating efficiency and neutral fraction) are well explained.

Fig. 5 Dependence of the maximum wind velocity (dashed line) and total planetary mass-loss rate (solid line) at the Roche radius as a function of the velocity dispersion at the planet surface (note that the two vertical axes do not correspond). Results come from 1D simulations of atmospheric outflows induced by MHD waves. Ranges of ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$ values derived from the observations are highlighted as blue (Visit 2) and red (Visit 3) shaded areas, allowing an estimation of the required velocity dispersions and corresponding total mass loss.

5.3. Planetary magnetic moment

Magnetic interactions between close-in planets and their host star may play a part in shaping planetary exospheres (e.g., Vidotto et al. 2011a; Matsakos et al. 2015). Analysis of radio emission from an exoplanet would constrain the planetary magnetic field strength, but there has been no confirmed detection to date (e.g., Bastian et al. 2000; Lazio et al. 2004; Smith et al. 2009; Hallinan et al. 2013, but see also Lecavelier des Etangs et al. 2013; Sirothia et al. 2014). Because the charged stellar wind protons can interact with the exosphere of an exoplanet and be deflected by its magnetosphere, observations in the Lyman-α line can potentially be used to estimate the planetary magnetic field in addition to the stellar wind properties. In their study of the hot Jupiter HD 209458b, this approach was chosen by Holmström et al. (2008) and Ekenbäck et al. (2010), who prescribed a magnetic obstacle surrounding the planet that prevented penetration by protons, and by Kislyakova et al. (2014a), who estimated the magnetic moment of the planet. In that case, however, the uncertainties and phase coverage of the Lyman-α observations lead to strong degeneracies between the relative contribution of the planetary escape, stellar wind, and magnetic obstacle properties, and a scenario with radiation pressure alone was even found to explain the observations well (Vidal-Madjar et al. 2003; Lecavelier Des Etangs et al. 2008; Bourrier & Lecavelier des Etangs 2013). Similar conclusions were reached by Ben-Jaffel & Ballester (2013) when adjusting the observations of ionized carbon around HD 189733b, with many possible solutions for the stellar wind and planetary magnetosphere properties.

We showed for GJ 436 b that the Lyman-α transit observations are well explained by the combination of radiative braking and stellar wind interactions. Therefore we did not include a magnetic obstacle in the EVE simulations because it would bring unnecessary additional free parameters to the interpretation of the data. It is nonetheless possible to use our simulation results to set an upper limit on the magnetic moment of GJ 436 b. Consider that the planet is surrounded by a magnetosphere, which would prevent stellar wind protons from interacting with the planetary exosphere. Since our best-fit simulations show significant stellar wind interactions down to about 20 R_p from the planet at the time of the transit (see Figs. 3 and 4), this corresponds to the largest possible stand-off distance R_s of the magnetosphere. At this altitude, a temperature of several 10⁶ K would be needed for the thermal pressure from the planetary wind to overcome the stellar wind forcing. Therefore we consider that the contribution from thermal pressure is negligible and that the magnetospheric size is set by the balance between the planetary magnetic pressure and the stellar wind ram pressure, yielding a planetary magnetic moment (e.g., Grießmeier et al. 2004; See et al. 2014) $\mathrm{ℳ}\mathrm{=}{\left(\frac{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{s}}^{\mathrm{6}}{\mathit{\rho }}_{\mathrm{wind}}{\mathit{v}}_{\mathrm{rel}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{0}}{\mathit{f}}_{\mathrm{0}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(7)where we assumed the planet to have a dipolar field, v_rel is the velocity of the stellar wind relative to the planet, μ₀ the permeability of vacuum, and f₀ ≈ 1.22 is a form factor of the magnetosphere. We set the mass density of the stellar wind at the planetary orbit, ρ_wind ≃ (2 − 5) × 10^-18 kg m^-3, to its average value at the center of the transit for Visits 2 and 3 (Table 4). At this time, the planet is close to its semi-minor axis, and its orbital velocity is nearly perpendicular to the radial bulk velocity of the stellar wind (Vidotto et al. 2011b), with v_rel = $\sqrt{{\mathit{v}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{bulk}}^{\mathrm{2}}}\mathrm{~}$ 140 km s^-1. From Eq. (7), we find $\mathrm{ℳ}\lesssim \mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{26}}{\hspace{0.17em}\mathrm{A}\hspace{0.17em}\mathrm{m}}^{\mathrm{2}}\mathrm{~}\mathrm{0.16}\hspace{0.17em}{\mathrm{ℳ}}_{\mathrm{Jup}}\mathit{,}$(8)where ℳ_Jup = 1.56 × 10²⁷ A m² is the magnetic moment of Jupiter (Grießmeier et al. 2004). For a planetary radius R_p = 0.64 R_Jup, this implies an upper limit for the equatorial magnetic field strength of GJ 436 b of $\frac{{\mathit{B}}_{\mathrm{p}}}{{\mathit{B}}_{\mathrm{Jup}}}\mathrm{=}\frac{\mathrm{ℳ}}{{\mathrm{ℳ}}_{\mathrm{Jup}}}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{Jup}}}\right)}^{-3}\lesssim \hspace{0.17em}\mathrm{0.6}\mathit{,}$(9)or B_p ≲ 2.7 G, assuming an equatorial magnetic field strength of B_Jup ~ 4.3G for Jupiter (Bagenal 2013). This is consistent with the magnetic field strength of ~1 G required to accelerate a magnetically driven outflow to the observed velocities (Sect. 5.2.2). We caution that the magnetic field strength we derive here is a conservative upper limit, and, at the strength of 1 G, a self-consistent 3D model would be required to explore the possible escape anisotropies. The magnetosphere of the planet is likely to lie well inside the 20 R_p estimate for the stand-off distance, as our analysis shows that stellar wind particles interact with the (unprotected) planetary exosphere at such distances.

6. Discussion: separating radiation pressure and stellar wind contributions

We discuss in this section the general possibility of distinguishing the contributions of radiation pressure and stellar wind interactions when analyzing Lyman-α line observations of an evaporating exoplanet exosphere.

• The most direct case to identify the contribution of the stellar windis when its projected velocity is lower than the Doppler width ofthe Lyman-α line, but beyond the highest velocity that can be reached under radiation pressure acceleration. An absorption signature beyond this limit can then be attributed to a population of neutralized protons independent of the planetary neutrals, as was the case for the hot Jupiter HD 189733b (Lecavelier des Etangs et al. 2012). Interestingly, no other Lyman-α absorption signature was detected for this planet at radiation-pressure-induced velocities, presumably because of high photo-ionization rates or massive stellar wind abrasion of the escaping planetary neutrals (Bourrier & Lecavelier des Etangs 2013).

• If interacting stellar wind protons have projected velocities higher than the Doppler width of the Lyman-α line, they can still be detected through their abrading effect. However, the phase coverage of the transit must then allow very different regions of the exosphere to be probed. For the hot Jupiter HD 209458b, observations after the optical transit have been obtained only once at low signal-to-noise ratio (Ehrenreich et al. 2008), and therefore available data can be interpreted either by radiative blow-out alone (Bourrier & Lecavelier des Etangs 2013) or by the addition of a ~400 km s^-1 stellar wind with higher escape rates (Holmström et al. 2008; Ekenbäck et al. 2010; Kislyakova et al. 2014a). In that case, the decrease in absorption depth caused by the fast wind abrasion can indeed be compensated for by variations in other parameters such as the planetary escape rate or stellar photoionization rate, leading to degeneracies in the stellar wind and planetary escape properties.

• When stellar wind and radiatively induced projected velocities overlap, the observed absorption results from a balance between planetary neutral hydrogen atoms and neutralized stellar wind protons. As in the previous case, different combinations of parameters such as proton density and planetary escape rate can yield similar absorption profiles at a given time, but Bourrier & Lecavelier des Etangs (2013) proposed that the two mechanisms may be distinguished by analyzing the spectro-temporal variations of the absorption profile. For radiation-pressure-driven mechanisms, a radiative blow-out creates a narrow cometary tail with a strong velocity gradient (e.g., for HD 209458b; Bourrier & Lecavelier des Etangs 2013), whereas radiative braking leads to a massive expansion of the exosphere and its dilution within a broad cometary tail (Bourrier et al. 2015a; see also Fig. 2). In both cases, the depth of the absorption at a given wavelength varies strongly over time during the transit. In contrast, for the stellar-wind-driven mechanism, we showed in Sect. 3.1.2 that a population of neutralized stellar wind protons is characterized by a time-stable repartition of the absorption depth with wavelength. These spectro-temporal variations of the absorption profile could not be analyzed for HD 209458b because of the reduced phase coverage and were only tentatively studied for HD 189733 b (Bourrier et al. 2013).

With a good phase coverage of the exospheric transit and very large absorption depths that magnify the spectro-temporal variations of the absorption, the observations of GJ 436 b allow us, for the first time, to clearly separate the contributions of radiation pressure and stellar wind and to probe the regions of the exosphere that are shaped by each mechanism.

7. Conclusion

We investigated the physical conditions in the exosphere of the warm Neptune GJ 436 b observed at three different epochs in the Lyman-α line. We independently interpreted the spectra at each epoch through direct comparison with 3D numerical simulations performed with the code EVE. While radiative braking was shown to play a major role in shaping the exosphere of GJ 436 b, specific features remained to be explained in the observations, such as sharp egresses in Visits 2 and 3.

To this aim, we studied the additional effect of the stellar wind and found that significant interactions with the exosphere are required to explain the observations. With low velocities in the order of ~85 km s^-1, the stellar wind both abrades the planetary exosphere and leads to the formation of a secondary population of neutralized protons that contribute to the observed absorption in the blue wing of the Lyman-α line. The combination of radiation pressure and stellar wind abrasion allows for the formation of a giant coma with a reduced front ahead of the planet that reproduces the observed early ingresses. Furthermore, the different dynamics of the neutralized proton and planetary neutral populations allow for the formation of a cometary tail that either produces a sharp (Visit 2) or a flat egress (Visit 3), depending on the relative balance between the two populations. This balance seems to be determined mainly by changes in the density of the stellar wind between the two epochs (from ~10³ cm^-3 in Visit 2 to 3 × 10³ cm^-3 in Visit 3) and a marginally lower planetary outflow velocity in the second epoch (~50 km s^-1 against 60 km s^-1), while the other physical conditions in the upper atmosphere of GJ 436 b are otherwise very stable for the planetary mass loss (~2.5 × 10⁸ g s^-1), hydrogen photoionization rate (~2 × 10^-5 s^-1), and stellar wind bulk (~85 km s^-1) and thermal velocities (~10 km s^-1).

Using EVE simulations, we detailed how these properties influence the structure of the exosphere and compared the best-fit results with independent theoretical estimations. Comparisons with energy-limited escape rates place constraints on the heating efficiency η of the upper atmosphere and its neutral hydrogen content f_H⁰, with ηf_H⁰ ~ 10^-2 and a lower limit on both η and f_H⁰ of about 0.5%. The properties of the stellar wind for the low-mass star GJ 436 are not unexpected, with lower density and bulk velocity than for earlier host stars, but it is the first time that these properties are directly measured from observations for a M dwarf. Future observations of extended atmospheres in similar systems therefore have a high potential for stellar wind characterization. On the other hand, the velocity of the escaping gas is faster than the outflows of evaporating hot Jupiters. The shallower gravity well of the Neptune-mass GJ 436 b may play some part, although its thermosphere is less irradiated. We showed that a possible mechanism for the fast ouflows might be turbulence-driven MHD waves at the surface of the planet. While the planetary magnetic field cannot be strongly constrained by Lyman-α observations, we used a geometric argument to set an upper limit on GJ 436 b of about a tenth of Jupiter magnetic moment.

We note that Visit 2 observations might also be explained through pure abrasion of the exosphere from a high-velocity stellar wind. However, this scenario involves escape rates expected from more strongly irradiated hot Jupiters and is not consistent with Visit 1 or Visit 3 observations. Given the stability of GJ 436 stellar Lyman-α line in four epochs covering four years of observations (see Ehrenreich et al. 2015) and the fact that Visits 2 and 3 can be explained with very similar conditions for the stellar irradiation, the stellar wind, and planetary mass-loss properties, this second scenario seems very unlikely. We also note that we were unable to adjust Visit 1 spectra with the model used in this paper, most probably because of the lack of reference at this epoch for the intrinsic stellar Lyman-α line, but also because the observations hint at more significant variations in the physical conditions of the exosphere at this epoch.

We emphasize that future models of the GJ 436 b exosphere should make use of the entire spectral content of the Lyman-α observations, distributed over 24 exposures and more than 600 data points. The exploration of the six-parameter space for the three epochs of observations and the different scenarii investigated required about 18 000 simulations, running full-time for nearly a year on 15 dual-processors compute nodes totalling 276 cores. EVE simulations show that the neutralized protons populating the exosphere would be photoionized much farther from GJ 436 b than planetary neutrals, making the transit of the cometary tail visible at Lyman-α for more than half the revolution period of the planet. New observations covering later phases than previously observed will allow for a full caracterization of the shape and properties of the exosphere, refining the measurements of the stellar wind properties.

Acknowledgments

This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work has been carried out in the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). The authors wish to thank the referee for their useful report. V.B., D.E., and A.A.V. acknowledge the financial support of the SNSF. A.L.E. acknowledges financial support from the Centre National d‚Etudes Spatiales (CNES). The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under program ANR-12-BS05-0012 “Exo-Atmos”. Y.A.T. thanks T. K. Suzuki for developing a simulation code for magnetically driven outflow and fruitful discussions.

References

All Tables

All Figures

	Fig. 1 Schematic representation of the neutral hydrogen cloud surrounding GJ 436 b, displaying the different regions of the exosphere. Colors distinguish between the two populations of hydrogen atoms with different origins that compose the cloud.
In the text

Fig. 2 GJ 436 spectra during Visits 1 (left column), 2 (central column), and 3 (right column). Spectra in the same row were measured at about the same phase of the planet orbital position. We show the spectra gathered over the full HST orbits. The shaded gray area corresponds to the range affected by ISM absorption and geocoronal emission, which was excluded from the fits. Black solid spectra are the reconstructed out-of-transit stellar line profiles. Dashed black spectra correspond to the best-fit theoretical spectra for Visits 2 (reported for comparison in Visit 1 epoch) and 3. Note that the transit of the exosphere has little effect on the red wing, where the best fits nearly overlap with the out-of-transit spectra. In the bottom part of the plot, light curves for each epoch show the observed and theoretical flux integrated between − 120 and − 40 km s-1. The solid black line shows the optical transit.

In the text

Fig. 3 Views of GJ 436 b exosphere (dots) within the orbital plane, seen from the top, at the time of the optical transit. The two panels correspond to two simulations performed with the same escape rate, planetary wind velocity, and photoionization rate, but without stellar wind interactions (left panel) or including them (right panel, corresponding to the best fit for Visit 2). The planet is represented by a small black disk, to scale, between the dashed black lines limiting the LOS toward the stellar disk. Arrows display the average velocity field of the neutral hydrogen atoms in the stellar rest frame, with particles colored as a function of the time they escaped the atmosphere (for planetary neutrals) or were created through charge exchange (for neutralized protons). With no stellar wind interactions the dynamics of the escaping gas is initially dominated by the planetary orbital velocity and later constrained by radiative braking. With charge exchange most protons are neutralized in the dense inner coma, abrading at the source planetary neutrals that would have fed the outer regions of the exosphere, but also creating a compact population of neutrals dominated by the radial bulk velocity of the stellar wind. It takes much more time for stellar gravity to eventually disperse this population of neutralized protons.

In the text

Fig. 4 Views of GJ 436 b exosphere (lower panels) and its theoretical absorption profile (upper panels). We caution that absorption in the shaded range ([− 40; 20] km s-1) cannot be observed from Earth and that the STIS LSF would spread the absorption signal from the exosphere over more data points. Time increases from left to right and is indicated in each plot. The dashed black lines limit the LOS toward the stellar disk, while the dotted line indicates the star-planet axis. Black arrows show the velocity field of the gas in the stellar rest frame, with H0 atoms in the orbital plane colored as a function of their LOS-projected velocity. Red lines in the upper panels show the absorption profiles from the exosphere simulated with Visit 3 best-fit parameters. The population of neutralized protons arising from the inner coma can be identified by its radial velocities and stable high-velocity absorption peak, while planetary neutrals remain present at the front of the coma and in the outer regions of the exosphere opposite the star with dynamics dominated by radiative braking. The blue line profiles in the upper panels correspond to the same simulation with no stellar wind, illustrating how charge exchange reduces the size of the exosphere before and after the optical transit but only weakly affects the absorption profile at the center of the transit. It also shows that neutralized protons produce a stable high-velocity absorption peak.

In the text

Fig. 5 Dependence of the maximum wind velocity (dashed line) and total planetary mass-loss rate (solid line) at the Roche radius as a function of the velocity dispersion at the planet surface (note that the two vertical axes do not correspond). Results come from 1D simulations of atmospheric outflows induced by MHD waves. Ranges of ${\mathit{v}}_{\mathrm{wind}}^{\mathrm{p}}$ values derived from the observations are highlighted as blue (Visit 2) and red (Visit 3) shaded areas, allowing an estimation of the required velocity dispersions and corresponding total mass loss.

In the text