© ESO, 2015

1. Introduction

1.1. Atmospheric escape

More than 40% of the known exoplanets orbit extremely close to their star (0.1 au), raising many questions about the state of their atmosphere and its interactions with the host star. Theoretical studies indicate that the intense stellar X-ray and extreme ultraviolet energy input into the lower thermosphere of hydrogen-rich planets can lead to a hydrodynamic expansion of the upper gas layers. If the flow of expanding hydrogen reaches high velocities or high altitudes, it can escape the gravitational pull of the planet in large amounts, leading to the “evaporation” or “blow-off” of the atmosphere (Lammer et al. 2003; Lecavelier des Etangs et al. 2004; Yelle 2006; García Muñoz 2007; Koskinen et al. 2007; Murray-Clay et al. 2009; Koskinen et al. 2013a,b).

Transit spectroscopy in the ultraviolet has proved an invaluable tool for studying the outer regions of an evaporating planet. Because of its expansion, the atmosphere produces a much deeper absorption than the planetary disk when observed in the lines of elements that are abundant at high altitudes and/or associated with strong electronic transitions. Observations in the bright stellar Lyman-α line of neutral hydrogen have led to the detections of evaporation from the hot-Jupiter HD 209458b (Vidal-Madjar et al. 2003), HD 189733b (Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013), the hot-Neptune GJ 436 b (Kulow et al. 2014; Ehrenreich et al. 2015), and the warm-Jupiter 55 Cnc b (Ehrenreich et al. 2012). Heavier metals and ions can be carried upward through collisions with the expanding flow of hydrogen, which was confirmed by detection of several such species at high altitudes around HD 209458 b (O⁰, C⁺, Vidal-Madjar et al. 2004; Linsky et al. 2010; Mg⁰, Vidal-Madjar et al. 2013; and tentatively Si²⁺, Linsky et al. 2010; Ballester & Ben-Jaffel 2015) and HD 189733 b (O⁰ and possibly C⁺, Ben-Jaffel & Ballester 2013), and allowed detection of an extended exosphere around the hot-Jupiter WASP-12 b through transit observations in the near-UV (Mg⁺, Fe⁺; Fossati et al. 2010; Haswell et al. 2012).

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 radiation pressure, stellar wind, and magnetic interactions (e.g., Llama et al. 2011; Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014). Observations at different epochs in the Lyman-α line of HD 189733 b led further to the first detection of temporal variations in the physical properties of an extended exosphere, which probably come from the energetic and/or dynamical influence of the host star (Lecavelier des Etangs et al. 2012; Bourrier et al. 2013).

Many close-in exoplanets are expected to show signs of evaporation. The impact of this phenomenon on the formation and the evolution of planets has been studied through energy diagrams (Lecavelier des Etangs 2007; Ehrenreich & Désert 2011; Lammer et al. 2009) and theoretical modeling applied to specific systems (e.g., Lopez et al. 2012; Lopez & Fortney 2013) or populations of exoplanets (e.g., Jin et al. 2014; Kurokawa & Nakamoto 2014). While massive gaseous giants like HD 209458b or HD 189733 b are subjected to moderate escape rate, leading to the loss of a few percentage points of the planet mass over the lifetime of the system, low-density planets like mini-Neptunes or super-Earths with a large volatile envelope may be significantly eroded by evaporation. This scenario is supported by trends in the exoplanet population: the lack of sub-Jupiter size planets (Lecavelier des Etangs 2007; Davis & Wheatley 2009; Beaugé & Nesvorný 2013) and the large number of close-in super-earths (Howard et al. 2012) are consistent with the formation of rocky planetary remnants through the evaporation of their hydrogen/helium envelope. It is thus crucial to observe the atmosphere of close-in Neptune-mass planets to understand in which conditions their evaporation can lead to the formation of atmosphereless remnant cores like CoRoT-7b and Kepler-10 b.

1.2. The case of GJ 436 b

In this framework, the warm Neptune GJ 436 b (Butler et al. 2004; Gillon et al. 2007) is a good target for transit observations in the Lyman-α line, as it orbits a bright star (V = 10.7) whose close proximity to the Earth (d = 10.3 pc) limits the effects of interstellar medium absorption in the core of this line (Ehrenreich et al. 2011). The low mass of GJ 436 b (M_p = 23.1M_Earth), the moderate irradiation from its M-dwarf host star (a = 0.029 au; Ehrenreich et al. 2015), and the puzzling eccentricity of its orbit (e = 0.16, Beust et al. 2012 and references therein) offer the opportunity to study a new evaporation regime. Kulow et al. (2014) reported the detection of a deep absorption signature in the stellar Lyman-α line starting near the end of the optical transit, thought to be caused by the occultation from a large cometary tail of neutral hydrogen H⁰. However, interpretation of these observations was misled by the use of 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 much deeper H i absorption that starts several hours before the optical transit (~2–4 h) and that may remain detectable for more than 20 h afterward. GJ 436 b is thus not only trailed by a cometary tail, but also surrounded by a giant coma of neutral hydrogen extending tens of planetary radii in front of the planet.

Our goal is to investigate the effects of radiation pressure and the basic physical parameters of an evaporating exoplanet on the structure of GJ 436 b hydrogen exosphere and its transmission spectrum. To this aim we reconstructed the intrinsic Lyman-α line of the host star and compared the results of 3D numerical simulations of extended atmospheres based on the EVaporating Exoplanets (EVE) code with the three sets of transit observations. Observations are presented in Sect. 2, and the reconstruction of the intrinsic stellar Lyman-α line is explained in Sect. 3. The EVE code is described in Sect. 4 and used in Sect. 5 to study the properties of GJ 436 b exosphere in the low radiation pressure field of its host star. In Sect. 6 we discuss the part played by radiation pressure in the observed features and investigate several mechanisms that could also intervene. We summarize our conclusions in Sect. 7.

2. HST/STIS observations in the Lyman-α line

We used the three existing transit data sets of GJ 436 b observations in the H i Lyman-α line (1215.6702 Å), all taken with similar settings with the Space Telescope Imaging Spectrograph (STIS) instrument onboard the Hubble Space Telescope (HST). The HST data consist of time-tagged spectra obtained with the G140M grating, with a spectral resolution of ~20 km s^-1 at 1215.67 Å. Data was reduced as described in Ehrenreich et al. (2015). Throughout the paper we use the term transit to refer to the transit of the extended exosphere of GJ 436 b, while we explicitely speak of optical transit to refer to the occultation caused by the planetary disk alone. Each data set was obtained through four HST orbits performed at different phases near the optical transit. During Visit 1 (December 2012; Kulow et al. 2014), one exposure was taken before the optical transit, one during and two after. Visits 2 and 3 (June 2013 and 2014; Ehrenreich et al. 2015) have two exposures before the optical transit, one during and one after. The times of mid-exposures were calculated with the ephemeris in Knutson et al. 2011, and are reported in Table 1.

3. Reconstruction of the Lyman-α line

3.1. Method

The observed Lyman-α line is affected by absorption from the interstellar medium (ISM), which needs to be corrected for in order to estimate the intrinsic stellar Lyman-α line profile as seen from the planet. It is used to calculate the radiation pressure on the escaping exosphere (which is dependent on the stellar flux), and as a reference to calculate the theoretical transmission spectra of the exosphere as seen with HST/STIS (Sect. 4.3). The intrinsic line can be reconstructed directly using the STIS observations, and we applied the same method to all visits independently (Wood et al. 2005; Ehrenreich et al. 2011; France et al. 2012; Bourrier et al. 2013). We calculate a model profile of the stellar Lyman-α line, which is absorbed by the interstellar hydrogen and deuterium, convolved by STIS line spread function (LSF), and then compared with an observed reference spectrum. For each visit we used as reference the spectrum obtained during the first HST orbit, assuming it is unabsorbed by the exosphere of GJ 436 b. We used the Bayesian Information Criterion (BIC; de Wit et al. 2012) to prevent over-fitting in our search for the best model for the Lyman-α and LSF profiles. Trying out different models for the stellar line and the LSF, we found that in all visits the lowest BICs were obtained with a single Voigt profile for the stellar emission line, and a single Gaussian profile for the LSF. The free parameters of the fit are thus the heliocentric velocity of the star γ_∗, the Lyman-α line profile parameters (maximum stellar flux at 1 au f_peak (1 au), Doppler width Δ v_D and damping parameter a_damp), the width of the LSF profile σ_LSF, and the ISM parameters for neutral hydrogen (column density log ₁₀ N(H i), Doppler broadening parameter b(H i), and velocity relative to the star γ(H i)_/∗). Although there may be multiple ISM components along the line of sight (LOS) toward GJ 436, we assumed their absorption of the stellar Lyman-α profile can be interpreted as that of a single component with average properties. We thus modeled the ISM opacity as the combination of two Voigt profiles for the atomic hydrogen and deuterium, separated by about 0.33 Å. We used a D i/H i ratio of 1.5 × 10^-5 (e.g., Hébrard & Moos 2003; Linsky et al. 2006) and assumed there is no turbulent broadening with b(H i) = b(D i)/$\sqrt{\mathrm{2}}$ (Wood et al. 2005).

A preliminary analysis of χ² variations was used to locate the best model, and revealed no multiple local minima in the parameter space. We then refined the best-fit parameter values and calculated their uncertainties using a Metropolis-Hasting Markov chain Monte Carlo (MCMC) algorithm. To better sample the posterior distribution in the case of non-linear correlations between parameters, we applied an adaptive principal component analysis to the chains and jumped the parameters in an uncorrelated space (Díaz et al. 2014). The step size was adjusted to ensure an acceptation rate of about 25%, and the system was analyzed with chains of 5 × 10⁵ accepted steps. The comparison between the model spectra and the observed spectra was performed over the entire spectral line, excluding the core where ISM absorption is saturated and noise from the airglow correction is high (the spectral range of the fits are given in Table 2). The medians of the posterior probability distributions were chosen as the final best-fit values for the model parameters, and their 1σ uncertainties were obtained by finding the intervals on both sides of the median that contain 34.15% of the accepted steps.

3.2. Results

The best-fit reconstructed Lyman-α stellar profiles are displayed in Fig. 1, with the corresponding model parameters given in Table 2. The reconstructed Lyman-α lines for Visits 2 and 3 are remarkably similar, with all fit parameters consistent within their 1σ uncertainties. On the other hand we found spurious values for the best-fit parameters in Visit 1, with the ISM absorption profile anomalously deep and strongly blue-shifted. This is consistent with a lower flux in the blue wing of the Lyman-α line caused by absorption from hydrogen escaping the planet, and confirms that the spectrum obtained during the first exposure in Visit 1 is already absorbed by the giant exosphere of GJ 436 b (Ehrenreich et al. 2015). Without a good knowledge of both wings of the line, it is not possible to reconstruct the intrinsic stellar line for Visit 1. However the atmospheric absorption occurs in a localized spectral range in the blue wing of the line, and Ehrenreich et al. (2015) showed that the Lyman-α flux measured at larger velocity blue shifts, and over the whole red wing of the line, is very stable between the three visits, which is in agreement with our line reconstructions in Visits 2 and 3. Hereafter we thus use the Lyman-α line profile obtained for Visit 2 as a proxy for Visit 1.

The best-fit parameters we obtained for Visits 2 and 3 correspond to heliocentric radial velocity for the ISM hydrogen of –1.7 ± 2.1 km s^-1 and –1.0 ± 2.7 km s^-1, respectively. Using the LISM Kinematic Calculator¹ (Redfield & Linsky 2008), we identified three ISM clouds near (less than 20°) the line-of-sight toward GJ 436 with radial velocities of 0.37 ± 1.38 km s^-1 (LIC cloud), 2.33 ± 0.76 km s^-1 (Leo cloud), and 5.81 ± 0.74 km s^-1 (NGP cloud). The ISM absorption derived from our reconstruction of the Lyman-α line thus corresponds more likely to the LIC cloud (<1σ), and may be consistent with the Leo cloud (<2σ).

Fig. 1 Lyman-α line profiles of GJ 436 for Visits 2 (top) and 3 (bottom). The yellow line shows the theoretical intrinsic stellar emission line profile as seen by hydrogen atoms escaping the planetary atmosphere, scaled to the Earth distance. It also corresponds to the profile of the ratio β between radiation pressure and stellar gravity (values reported on the right axis). The solid black line shows the Lyman-α line profile after absorption by the interstellar hydrogen (1215.6 Å) and deuterium (1215.25 Å), whose cumulated profile is plotted as a dotted black line (ISM absorption in the range 0–1 has been scaled to the entire vertical axis range). The line profile convolved with the LSF (blue line) is compared to the observations shown as a red histogram. Above the horizontal dashed line (β = 1), radiation pressure would overcome stellar gravity and accelerate particles away from the star. However with β values below unity, radiation pressure from GJ 436 can only brake the gravitational deviation of particles toward the star.

4. EVE: the EVaporating Exoplanets code

EVE is a 3D numerical code that we developed to calculate the structure of an exoplanet upper atmosphere and the transmission spectrum of the species it contains. The direct comparison of these spectra with observational data allows the inferrance of several physical properties of both the planetary atmosphere and the star, and the code has been used in this way to study the neutral hydrogen exospheres of the hot Jupiters HD 209458b and HD 189733b (Bourrier & Lecavelier des Etangs 2013), the warm Jupiter 55 Cnc b, and the super-Earth 55 Cnc e (Ehrenreich et al. 2012), as well as the neutral and ionized magnesium populations in the thermosphere/exosphere of HD 209458b (Bourrier et al. 2014). The EVE 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. In the present paper, we further investigate the effects of stellar radiation pressure, stellar photo-ionization, and the planetary escape properties on the structure of the exosphere and its transmission spectrum. The following subsections offer a synthetic description of the model used in this frame.

4.1. Bottom layers of the atmosphere

The upper atmosphere is divided in two different regimes, joined at a transition altitude set at R_trans = 4 R_planet, approximately the mean altitude of the Roche lobe. We assume that the atmosphere within the Roche lobe, which is by definition under the gravitational influence of the planet, keeps a global cohesion and can be described analytically with a radial density profile corresponding to a hydrostatic equilibrium. In the numerical simulation, the density profile of neutral hydrogen in this part of the atmosphere is reassessed every time step by adjusting its value at R_trans to match the hydrogen density obtained in the upper atmospheric layers (see Sect. 4.2). Because of the eccentricity of the planetary orbit (e = 0.16), the density profile thus varies over time. The gas filling the Roche Lobe mainly contributes to the absorption in the core of the Lyman-α line, which cannot be observed because of ISM absorption and geocoronal emission. Therefore, this part of the atmosphere has little influence on the observed spectra, except through the Lorentzian wings of the absorption profile (the so-called damping wings, e.g., Ben-Jaffel 2008). This is however excluded because no significant absorption signal was observed in the red wing of the Lyman-α line during the optical transit of GJ 436 b (Ehrenreich et al. 2015).

4.2. Upper layers of the atmosphere

We use Monte-Carlo particle simulations to compute the dynamics of the escaping gas above R_trans. Neutral hydrogen atoms are represented by meta-particles, and the total number of metaparticles launched every time step dt, set to 5 mn, depends on the escape rate of neutral hydrogen along the eccentric orbit. To take the eccentricity into account, we assumed that GJ 436 b is in the energy-limited regime of evaporation, and the escape rate is proportional to the stellar energy input into the atmosphere. We thus varied the escape rate as the inverse square of the distance to the star, setting the reference Ṁ_H⁰ at the semi-major axis. The particles are released from the entire atmosphere and their initial velocity distribution, relatively to the planet, is the combination of the upward bulk velocity of the planetary wind v_pl−wind and an additional thermal speed component from a Maxwell-Boltzman velocity distribution (see Bourrier & Lecavelier des Etangs 2013).

The dynamics of the particles is calculated in the stellar reference frame and results from the stellar and planetary gravities, the stellar radiation pressure, and the inertial force linked to the non-Galilean reference frame. Radiation pressure on an atom is the combination of a radial² impulsion in the direction opposite to the star, caused by the absortion of a photon, followed by an isotropic impulsion after the photon reemission. We process the cumulated radial impulses from all photons absorbed by a hydrogen atom during dt as a force proportional to stellar gravity (Lagrange et al. 1998), and apply an equivalent number of velocity changes with constant speed 3.26 m s^-1 (this value comes from the conservation of momentum) and random direction to reflect the isotropic impulsions. A crucial aspect of radiation pressure is its dependence on a hydrogen atom radial velocity, as the number of photons it absorbs and reemits is proportional to the flux in the Lyman-α line Doppler-shifted at this velocity. The effect of radiation pressure thus strongly depends on the shape of the line profile, and the single peak profile obtained for GJ 436 b (Fig. 1) implies that the strength of radiation pressure monotonously decreases from the core to the wings of the line, ie with increasing absolute radial velocities.

The hydrogen atoms lifetime depends on the stellar photo-ionization rate Γ_ion which varies as the inverse square of the distance to the star. Once ionized, hydrogen atoms are removed from the simulation.

We do not take collisions into account in the upper atmospheric layers, as neutral hydrogen densities reach about 4 × 10⁵ m^-3 at R_trans for our best-fit simulations. With our simulation settings this value corresponds to a Knudsen number (the ratio between the mean free path of the gas and the atmospheric scale height) of about 40, which is above the level of the exobase defined by a Knudsen of one, even if the neutral fraction of hydrogen in the upper atmosphere of GJ 436 b is about 0.1 (Kulow et al. 2014; Koskinen et al. 2013a).

4.3. Theoretical spectra

Self-shielding for radiation pressure and photo-ionization is taken into account by calculating their respective opacities in the bottom and upper atmospheric layers. Similarly, EVE calculates the opacity along the star-Earth line of sight to obtain at each time step the theoretical absorption spectra in the Lyman-α line (for details regarding the calculation of these opacities, see Bourrier & Lecavelier des Etangs 2013). During the optical transit, the stellar flux is absorbed at every wavelength by the planetary disk. In contrast, the neutral hydrogen atmosphere may contribute to the absorption on longer time scales and at specific wavelengths that depends on the dynamics of the gas and the effects of thermal and natural broadening. At a given time step in the simulation, the theoretical intrinsic stellar line profile is multiplied by the absorption profile of the planet and its atmosphere, calculated with a resolution Δv = 10 km s^-1 (more than twice the resolution of the STIS spectra). The resulting spectrum is then attenuated by the ISM absorption and convolved with the LSF (using the fits obtained in Sect. 3.2), and interpolated to the observed velocities. To compare the simulation results with the STIS observations, we average the theoretical spectra calculated during the time window of each observation.

Fig. 2 Views of GJ 436 b exosphere (blue dots) from the perpendicular to the orbital plane. The planet is the small disk at the intersection of the star/Earth LOS (dashed black line) and the planetary orbit (green line). All system properties are similar in the three panels, except for radiation pressure. Upper panel: without radiation pressure, stellar gravity shears the extended cloud of hydrogen, and gas at shorter orbital distance than the planet falls toward the star. Middle panel: with even a moderate radiation pressure (~60% of stellar gravity at maximum), all the escaping gas decelerates and moves to larger orbits, forming a longer comet-like tail trailing behind the planet. This is the real radiation pressure that corresponds to the observed Lyman-α flux. Lower panel: with a high radiation pressure (~400% of stellar gravity at maximum), the escaping gas is swiftly blown away. The planet is surrounded by a much smaller coma but trailed by a narrower, more radial, cometary-tail.

5. Structure of the exosphere

While mechanisms such as stellar wind or bow-shocks may play a role in shaping the exosphere of an exoplanet, their description requires strong assumptions and the introduction of several free parameters in models of atmospheric escape. By contrast radiation pressure is known to be present with a well defined behavior, and we showed in Sect. 3 that it is very well constrained from direct measurements of the stellar Lyman-α line. In this paper we thus focus on the effects of radiation pressure alone on the structure and the spectral signature of the exosphere.

5.1. Radiation pressure: shaping the cloud

Fig. 3 View of GJ 436 b exosphere (gray dots) from the perpendicular to the orbital plane, 3.5 h after the optical transit. The white dashed line shows the LOS toward Earth. Arrows display the velocity field of the hydrogen atoms with respect to the star, colored as a function of the time they escaped the atmosphere (counted from the center of the optical transit). In the reduced gravity field from the star, particles move beyond the orbit of the planet. However, radiation pressure does not overcome stellar gravity, and particles are still deviated from the tangent to the orbit at the time of their escape (solid colored lines). Because radiation pressure is velocity-dependent, the particles do not stay on Keplerian orbits. To better illustrate the effect of radiative braking, this simulation was performed with a low ionization rate artificially lengthening the exospheric tail.

If the atmosphere of an exoplanet is significantly expanded, various regions of the exosphere move on very different orbits and revolve around the star in differential rotation. In the absence of any radiation pressure from the star, gas escaping beyond the Roche Lobe would follow ballistic trajectories subjected to stellar gravity alone. The exosphere would suffer shear, the gas within the orbit of the planet accelerating and falling toward the star while the gas beyond would decelerate and move toward larger distances (top panel in Fig. 2). However, this is an unrealistic case for neutral hydrogen escape, as radiation pressure from the Lyman-α line will always compensate for some fraction of the stellar gravity. In the case of close-in planets around K-type stars and earlier, such as HD 189733b and HD 209458b, the Lyman-α line is actually strong enough for radiation pressure to overcome up to 3–5 times the gravity, thus swiftly accelerating the H⁰ atoms away from the star at radial velocities exceeding 100 km s^-1 (bottom panel in Fig. 2). The exosphere is compressed on its dayside, while gas escaping at the limbs of the atmosphere is blown away and shaped into a narrow cometary tail, bent toward the radial direction of the stellar photons (e.g., Schneider et al. 1998; Vidal-Madjar et al. 2003; Bourrier & Lecavelier des Etangs 2013). Orbiting a M dwarf at moderate orbital distance, the warm Neptune GJ 436 b comparatively receives a lower flux from the stellar Lyman-α line. With radiation pressure compensating for up to ~60% of the gravitational attraction at low radial velocities (Fig. 1), the escaping hydrogen atoms cannot be accelerated away from the planet but are nonetheless subjected to a strong radiative braking. Indeed, even a small reduction of stellar gravity is enough for all the exospheric gas to decelerate and move to larger orbits, preventing the formation of a stream of matter infalling toward the star. In contrast to highly irradiated planets, the gas can diffuse within a larger volume around the planet and forms a wider cometary tail closer to the planet orbit (middle panel in Fig. 2).

5.2. Effect of radiation pressure on the velocity field

Radiation pressure has also a major effect on the velocities of the escaping hydrogen atoms that correspond to the wavelength range of the observed absorption profile. In the stellar reference frame, the initial velocity of an escaping particle is the combination of the planet orbital velocity and the velocity of the planetary wind at the Roche Lobe. For close-in exoplanets (a< 0.1 au), orbital velocity is generally higher than 100 km s^-1, significantly higher than planetary wind velocity in the order of 1–25 km s^-1 with respect to the planet (Yelle 2004; Koskinen et al. 2013a; Bourrier et al. 2014). To a first approximation the average dynamics of the gas in the exosphere is thus set by the planet orbital velocity v_pl. With even a low radiation pressure, hydrogen atoms decelerate with respect to the planet. Because radiation pressure does not fully overcome the star gravity, hydrogen atoms nevertheless remain attracted toward the star and will not move beyond the plane tangent to the planet orbit at the time of their escape, as illustrated in Fig. 3. As a result of both the orbital velocity of the planet and the effect of radiation pressure, the projected velocity of the gas thus ranges approximately between the projection of v_pl (for atoms close to the planet) and –135 km s^-1, which corresponds to the highest planet orbital velocity. Therefore, this dynamics is consistent with the detection of blue-shifted absorption up to about –120 km s^-1 during the transit of GJ 436 b exosphere (Ehrenreich et al. 2015).

Fig. 4 Views of GJ 436 b exosphere (lower panels) and its theoretical absorption profile, as it would be observed close to the star (upper panels). We caution that absorption in the range [–40; 20] km s-1 cannot be observed from Earth, and that STIS LSF would spread this absorption signal over several pixels. Time increases from left to right and is indicated on each plot. The dashed black lines limit the LOS toward the stellar surface, while the dotted line indicates the star-planet axis. Black arrows show the velocity field of the H0 atoms relative to the star. Gas is seen in the orbital plane and colored as a function of its velocity projected on the LOS. Left panel: as they circumvent the planet to move away from the star, atoms in the front of the coma absorb the stellar flux at blueshifted velocities before the optical transit. Middle panel: at the time of the optical transit the strongest absorption comes from the core of the coma, and is centered on the projection of the planet orbital velocity at –16 km s-1. The outer regions of the exospheric tail contribute to absorption at higher velocities in the blue wing. Right panel: after the optical transit, absorption in the blue wing of the line comes exclusively from the tail.

5.3. Influence of the physical parameters

We used simulations performed with EVE to study the effects of the most influential star and planet properties on the structure of the exosphere and its absorption signature. The planetary escape rate and the stellar photo-ionization rate are two basic properties of neutral hydrogen escape, but in contrast to hot-Jupiters subjected to strong radiation pressure the planetary wind velocity also has a significant influence on the exosphere of a moderately irradiated, lower-mass planet like GJ 436 b.

• Planetary wind velocity v_pl−wind: this parameter influences the spatial and spectral dispersion of the gas in the exosphere. Higher values for v_pl−wind increase the width of the coma and the exospheric tail, extending in particular the duration of the pre-transit and post-transit absorptions. The same amount of gas escaping faster from the atmosphere is more diffused and has a wider distribution of Doppler velocities, yielding lower depths in the core of the absorption profile but spreading the absorption over a larger spectral range.

• Photo-ionization rate Γ_ion: this parameter acts on the lifetime of neutral hydrogen atoms, and yields local variations in the size of the exosphere and its density. Because of self-shielding, photo-ionization has indeed more impact on the regions facing the star, which leads to a more significant erosion of the tail and the front of the exosphere than its denser coma. Higher values for Γ_ion thus decrease more strongly the pre-transit and post-transit absorption depths.

• Escape rate Ṁ_H⁰: this parameter changes the local density roughly uniformly within the exosphere. An increase in the escape rate can thus partly balance a stronger photo-ionization as long as self-shielding is not significant. Higher values for Ṁ_H⁰ increase the depth of the absorption profile at all wavelengths, and lengthen the exospheric tail.

When hydrogen densities are high enough, self-shielding can protect the interior of the coma from radiation pressure. Gas in these regions and escaping from the planet dayside is subjected to the nominal stellar gravity, and can move far ahead of the planet toward the star. We note however that the very high escape rate/very low photo-ionization rate that would lead to the formation of this extended leading tail are excluded, as it would occult the star too early compared to the observations.

6. Interpretation of the observations

6.1. Qualitative description

The main features of the Lyman-α line observations of GJ 436 b were described by Ehrenreich et al. (2015). They showed that the flux in the entire red wing, and at high velocities in the blue wing, is stable over the transit duration and between the different epochs. Visit 1 displays a puzzling absorption in the red wing just after the optical transit, which is not reproduced in subsequent visits. Absorption attributed to the exosphere occurs for all visits in a similar range of velocities up to about –120 km s^-1 in the blue wing of the line. Overall the light curves measured in this range (Fig. 5) present similar variations over time for the three visits: 1) the exospheric transit begins about 2 h before the optical transit; 2) the maximum absorption depth is found at about the time of the optical transit; 3) although it decreases after the optical transit, absorption may still be observed for several hours afterwards. There are nevertheless some differences between the visits in the way the absorption varies over time. Visit 1 displays smooth variations during the atmospheric ingress and egress. The decrease in flux is sharper at the ingress of the other visits, in particular in Visit 3. Regarding the post-transit phases, absorption decreases very quickly in Visit 2 while Visit 3 displays absorption depths similar to those of the optical transit.

Fig. 5 Transit light curves of the Lyman-α line in Visit 1 (upper panel), Visit 2 (middle panel), and Visit 3 (lower panel), integrated between –120 and –40 km s-1. Black points correspond to the observations, while the black line shows the optical transit. The blue lines corresponds to best-fit simulations when the exosphere is subjected to stellar gravity and radiation pressure alone.

6.2. Possible scenarios

In this section we investigate different physical scenarios that could explain the features described in Sect. 6.1.

6.2.1. Radiative braking

As shown in Sect. 5.2, the velocities expected from the dynamics of the radiatively braked exosphere are consistent with the range of the absorption signature up to –120 km s^-1. Using the model described in Sect. 4, we further found that the combination of low radiation pressure with high planetary wind velocities of about 30 km s^-1 allows the diffusion of the escaping gas within a coma large enough to explain the deep absorption observed in all visits. This wind velocity, higher than those obtained from the modeling of hot-Jupiters expanding thermosphere (~1–10 km s^-1; see, e.g., Murray-Clay et al. 2009; Koskinen et al. 2013a), may be linked to the lower density of the Neptune-mass planet GJ 436 b, or caused by a magnetically driven planetary wind (Tanaka et al. 2014).

The extended coma also produces a deeper and earlier ingress than observed, in particular for Visit 3 (Fig. 5). We also found that the photo-ionization of the gas and its natural dilution as it moves away from the planet leads to a stronger decrease of the absorption depth than observed during the post-transit phases in Visit 3. Therefore, while the exosphere is always subjected to radiation pressure and this force naturally reproduces the observed velocities of the escaping gas, radiative braking alone does not explain the time variations of the measured absorption depth. In other words, radiation pressure is the main driver of the velocity field in the exosphere, and secondary mechanisms must reshape the geometry of the gas cloud.

6.2.2. Radiative blow-out

As is the case for hot Jupiters, a higher radiation pressure could repel the escaping hydrogen and reduce the front of the coma responsible for an extended early ingress. However, because radiation pressure is proportional to the stellar flux, this would require a higher flux in the entire Lyman-α line, implying that the spectra measured during the first orbits in Visits 2 and 3 are already absorbed by the exosphere. This is unlikely, because these first orbit spectra are remarkably similar (as would be the case if they corresponded to the stable intrinsic stellar Lyman-α line), in contrast to other spectra measured at similar phases of Visits 2 and 3 (Fig. 5). Furthermore to explain absorption as early as the first orbits in Visits 2 and 3, the coma should be even more extended than is the case with a low radiation pressure, in contradiction with its narrowing from the stronger radiation pressure.

6.2.3. Transition region

Ehrenreich et al. (2015) used the EVE code with similar settings as those described in Sect. 4 to fit the three combined visits of GJ 436 b in the frame of a radiative braking scenario. They reproduced the observed sharp early ingress by introducing an additional free parameter to the model, which is a multiplicative factor to artificially increase the ionization cross-section of H⁰ atoms. The corresponding increase in self-shielding coupled with high ionization rates allows for the thin outer envelop of the exosphere facing the star to be eroded, while preserving the denser inner regions. The sharp density transition region thus created is akin to the boundary of a H ii region, or to a bow-shock environment. However this model requires extremely high ionization rates in the order of 10^-6 s^-1 at 1 au from the star, and makes strong assumptions on self-shielding that affect the overall structure of the exosphere. Detailed models of magnetic interactions between the star and planet will be required to investigate this scenario further.

6.2.4. Stellar wind interactions

Charge exchange between the escaping H⁰ atoms and stellar protons could also ionize the front layers of the coma (e.g., Holmström et al. 2008) and reduce the depth of the early ingress. These interactions would also ionize planetary neutral atoms in every region of the exosphere, decreasing the overall absorption, but this might be compensated for by the concomitant creation of an independent population of neutralized protons repopulating the exosphere. Moving with the high radial velocity of the stellar wind, this population would be less sensitive to photo-ionization and spatial dilution than the planetary atoms and would keep the same distribution over longer time scales, possibly explaining the post-transit absorption depths in Visit 3 (Sect. 6.1). A combination of the stable radiation pressure with a variable stellar wind may explain the different variations in absorption depths observed between the three visits.

7. Discussion

New HST observations of the Neptune-mass exoplanet GJ 436 b revealed a spectacular atmospheric escape (Kulow et al. 2014; Ehrenreich et al. 2015), repeatable over the three epochs of observations. Where the exospheres of known evaporating hot Jupiters cover ~15% of the stellar disk, the cloud of hydrogen surrounding GJ 436 occults up to ~60% of the star with significant absorption observed several hours before and after the optical transit. The amplitude of the absorption and the corresponding size of the exosphere may appear surprising, considering that GJ 436 b is only moderately irradiated by its M-dwarf host star.

To understand this puzzle, we have adapted and run EVE, our 3D particle simulation of atmospheric escape, with the specific properties of the GJ 436 system. In this paper we focused on a scenario where escaping atoms are submitted to radiation pressure, stellar and planet gravitational forces, and photo-ionization to assess the impact of these phenomenon on the observed Lyman-α transit signature.

Radiation pressure is a velocity-dependent force that can be calculated directly from observations of the Lyman-α line, and we provide detailed reconstructions of the intrinsic stellar line at the different epochs. We found that radiation pressure from GJ 436 is too low to exceed stellar gravity and to repel the escaping hydrogen atoms away from the star, as is the case for hot Jupiters. Even this moderate radiation pressure is nevertheless enough to brake the gravitational deviation of the gas toward the star, and it allows the formation of a large coma comoving with the planet (although not gravitationally bound to it) that extends in a broad cometary tail trailing the planet and progressively eroded by photo-ionization.

Our simulations show that the low radiation pressure from GJ 436 explains the velocity range of the observed absorption signatures well up to about –120 km s^-1 and suggest that a high planetary wind velocity of about 30 km s^-1 is required to create a sizeable coma reproducing the observed absorption depths. However, we found that the variations in the depth and duration of the absorption at the different phases of the transit, some specific to the visit, cannot be explained by radiation pressure alone. In the low radiation pressure field of the star, the coma also extends far ahead of the planet, producing a deeper and earlier ingress than observed. An additional mechanism, such as interactions with the stellar wind, is thus needed to reduce the size of the coma front. By creating a new population of neutral hydrogen atoms with the properties of the stellar wind, charge exchange may also explain peculiar features observed in the last visit.

We emphasize that future studies intent on reproducing the observations of GJ 436 b should take the entire spectral information into account and not only the photometric light curves. Given the size of the exosphere, the regions that were observed transiting the star only sample a small part of the atmospheric structure, and more observations will be required to fully constrain the dynamics of the escaping gas and the properties of GJ 436 b environment.

HST observations of a partially transiting exosphere around the warm giant 55 Cancri b (Ehrenreich et al. 2012) previously hinted that moderately irradiated planets could have significantly extended exospheres. The detection of GJ 436 b exosphere offers the opportunity to probe the atmospheric escape of a “cool” exoplanet in detail (~800 K), whose atmospheric photochemistry may also be significantly affected by Lyman-α radiation (Miguel et al. 2015). Generalizing these results to other systems, in particular those that will be detected around M dwarfs by CHEOPS, TESS, and PLATO, will open tantalizing new perspectives for characterizating warm and temperate exoplanets.

Acknowledgments

We thank the referee, Jeffrey Linsky, for useful and fruitful comments. 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 bas been carried out in the framework of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). V.B. and D.E. acknowledge the financial support of the SNSF. A.L.E acknowledges financial support from the Centre National d’Études Spatiales (CNES). The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under program ANR-12-BS05-0012 “Exo-Atmos”. We warmly thank J. B. Delisle and A. Wyttenbach for fruitful discussions.

References

All Tables

All Figures

Fig. 1 Lyman-α line profiles of GJ 436 for Visits 2 (top) and 3 (bottom). The yellow line shows the theoretical intrinsic stellar emission line profile as seen by hydrogen atoms escaping the planetary atmosphere, scaled to the Earth distance. It also corresponds to the profile of the ratio β between radiation pressure and stellar gravity (values reported on the right axis). The solid black line shows the Lyman-α line profile after absorption by the interstellar hydrogen (1215.6 Å) and deuterium (1215.25 Å), whose cumulated profile is plotted as a dotted black line (ISM absorption in the range 0–1 has been scaled to the entire vertical axis range). The line profile convolved with the LSF (blue line) is compared to the observations shown as a red histogram. Above the horizontal dashed line (β = 1), radiation pressure would overcome stellar gravity and accelerate particles away from the star. However with β values below unity, radiation pressure from GJ 436 can only brake the gravitational deviation of particles toward the star.

In the text

Fig. 2 Views of GJ 436 b exosphere (blue dots) from the perpendicular to the orbital plane. The planet is the small disk at the intersection of the star/Earth LOS (dashed black line) and the planetary orbit (green line). All system properties are similar in the three panels, except for radiation pressure. Upper panel: without radiation pressure, stellar gravity shears the extended cloud of hydrogen, and gas at shorter orbital distance than the planet falls toward the star. Middle panel: with even a moderate radiation pressure (~60% of stellar gravity at maximum), all the escaping gas decelerates and moves to larger orbits, forming a longer comet-like tail trailing behind the planet. This is the real radiation pressure that corresponds to the observed Lyman-α flux. Lower panel: with a high radiation pressure (~400% of stellar gravity at maximum), the escaping gas is swiftly blown away. The planet is surrounded by a much smaller coma but trailed by a narrower, more radial, cometary-tail.

In the text

Fig. 3 View of GJ 436 b exosphere (gray dots) from the perpendicular to the orbital plane, 3.5 h after the optical transit. The white dashed line shows the LOS toward Earth. Arrows display the velocity field of the hydrogen atoms with respect to the star, colored as a function of the time they escaped the atmosphere (counted from the center of the optical transit). In the reduced gravity field from the star, particles move beyond the orbit of the planet. However, radiation pressure does not overcome stellar gravity, and particles are still deviated from the tangent to the orbit at the time of their escape (solid colored lines). Because radiation pressure is velocity-dependent, the particles do not stay on Keplerian orbits. To better illustrate the effect of radiative braking, this simulation was performed with a low ionization rate artificially lengthening the exospheric tail.

In the text

Fig. 4 Views of GJ 436 b exosphere (lower panels) and its theoretical absorption profile, as it would be observed close to the star (upper panels). We caution that absorption in the range [–40; 20] km s-1 cannot be observed from Earth, and that STIS LSF would spread this absorption signal over several pixels. Time increases from left to right and is indicated on each plot. The dashed black lines limit the LOS toward the stellar surface, while the dotted line indicates the star-planet axis. Black arrows show the velocity field of the H0 atoms relative to the star. Gas is seen in the orbital plane and colored as a function of its velocity projected on the LOS. Left panel: as they circumvent the planet to move away from the star, atoms in the front of the coma absorb the stellar flux at blueshifted velocities before the optical transit. Middle panel: at the time of the optical transit the strongest absorption comes from the core of the coma, and is centered on the projection of the planet orbital velocity at –16 km s-1. The outer regions of the exospheric tail contribute to absorption at higher velocities in the blue wing. Right panel: after the optical transit, absorption in the blue wing of the line comes exclusively from the tail.

In the text

	Fig. 5 Transit light curves of the Lyman-α line in Visit 1 (upper panel), Visit 2 (middle panel), and Visit 3 (lower panel), integrated between –120 and –40 km s-1. Black points correspond to the observations, while the black line shows the optical transit. The blue lines corresponds to best-fit simulations when the exosphere is subjected to stellar gravity and radiation pressure alone.
In the text