© ESO, 2016

1. Introduction

Gas giants are ubiquitous among the known extrasolar planets, because they produce rather strong transit signals that can be easily detected. Today, hot gas planets with semimajor axes smaller than 0.1 au constitute more than 15% of all verified planets (e.g., exoplanets.org, Wright et al. 2011). Such hot planets have been found as close as two stellar radii above the photosphere of their host stars (Hebb et al. 2009) and some of them also orbit highly active host stars in close proximity (Alonso et al. 2008; Schröter et al. 2011). This places their high-energy irradiation level as much as 10⁵ times higher than the irradiation of Earth. The detection of gas planets in such extreme environments immediately raises questions about the stability of their atmospheres, about possible evaporation, and about the resulting lifetime of these planets (e.g., Schneider et al. 1998; Lammer et al. 2003; Lecavelier des Etangs 2007).

In the case of close-in planets, the absorption of high-energy radiation causes temperatures of several thousand degrees in the planetary thermosphere¹. The onset of the temperature rise at the base of the thermosphere was confirmed in several observational studies of HD 209458 b and HD 189733 b (Ballester et al. 2007; Sing et al. 2008; Vidal-Madjar et al. 2011b; Huitson et al. 2012). Even at 1 au the thermosphere of Earth is heated to up to 2000 K (e.g., Banks & Kockarts 1973). Under these conditions the upper atmospheric layers of hot gas planets are probably unstable and expand continuously, unless there is a confining outer pressure (Watson et al. 1981; Chassefière 1996). The formation of such a hydrodynamic planetary wind is similar to the early solar wind models by Parker (1958), with the major difference that the energy source for the wind is the external irradiation. Watson et al. (1981) also showed that under certain circumstances irradiated atmospheres undergo an energy-limited escape, which simply means that the radiative energy input is completely used to drive the planetary wind. Naturally one has to account for the fact that the heating efficiency of the absorption is less than unity. Values ranging from 0.1 to 1.0 have been used for the heating efficiency in the literature. Recently, Shematovich et al. (2014) performed a detailed analysis of the atmosphere of HD 209458 b and found a height-dependent value between 0.1 and 0.25.

In the case of energy-limited escape the mass-loss rate is proportional to the XUV flux² impinging on the atmosphere and inversely proportional to the planetary density (Sanz-Forcada et al. 2011). Since the mass-loss rate only depends on the planetary density but not on its mass, small planets with low densities are most heavily affected by the fractional mass-loss rate. While for most detected hot gas planets the predicted mass-loss rates do not result in the loss of large portions of their total mass over their lifetime (Ehrenreich & Désert 2011), a complete evaporation of volatile elements is possible for smaller gas planets. Indeed this might be the reason for the lack of detected hot mini-Neptunes (Lecavelier des Etangs 2007; Carter et al. 2012).

The existence of expanded thermospheres in hot gas planets was first confirmed in the exoplanet HD 209458 b, where Vidal-Madjar et al. (2003) found absorption of up to 15% in the line wings of the stellar hydrogen Lyα line during the transit of the hot Jupiter. Since the optical transit depth in this system is only 1.5% (Henry et al. 2000; Charbonneau et al. 2000), the excess absorption is most likely caused by a neutral hydrogen cloud around the planet. As Vidal-Madjar et al. also noted, the variability of solar-like stars is not sufficient to mimic this absorption signal. Although the Lyα emission of the Sun is strongly inhomogeneously distributed over the solar disk (e.g., see images from The Multi-Spectral Solar Telescope Array, MSSTA), this is only true for the line core, which is produced at temperatures of around 20 000 K in the solar chromosphere (Vernazza et al. 1973). The line wings originate deeper within the chromosphere at around 6000 to 7000 K. Emission at these temperatures is more homogeneously distributed over the solar disk, as seen for example by images in Hα or Ca ii K filters – lines that are produced at similar temperatures in the chromosphere. In analogy to the Sun, it is reasonable to assume that for solar-like stars the variability in the calcium line cores is also indicative of the variability in the wings of the Lyα line, where the planetary absorption is detected.

The variability of stars in the Ca ii H&K emission lines has been studied in detail for many stars (e.g. Vaughan et al. 1981). In particular, the inactive star HD 209458 (Henry et al. 2000) shows a fractional standard deviation of only 1.0% over several years compared to the solar variation of 14.4% (Hall et al. 2007), and the short term variability is on the same order (using data from the California Planet Search program, Isaacson & Fischer 2010). For this star the variability in the Lyα line wings has been studied directly by Ben-Jaffel (2007), who also found a variability of about 1% (with an uncertainty of 3%) in their three out-of-transit time bins. Compared to the Lyα absorption signal this is small, therefore, variability or spatial inhomogeneity of the stellar emission is unlikely to be the dominating source of the absorption signal seen during the planetary transit.

In the case of HD 209458 b the confidence in the detection of an expanded atmosphere was further increased by several transit observations, which revealed excess absorption in spectral lines of H i, O i, C ii, Si iii, and Mg i (Vidal-Madjar et al. 2004, 2013; Ballester et al. 2007; Ehrenreich et al. 2008; Linsky et al. 2010; Jensen et al. 2012), although the detection of silicon was recently questioned (Ballester & Ben-Jaffel 2015). There is also an indication of possible absorption in a Si iv line during the transit of this planet (Schlawin et al. 2010). Today, three additional systems exist with a comparable coverage of observations. In the more active host star HD 189733 excess absorption was detected in H i, O i, and possibly in C ii (Lecavelier des Etangs et al. 2010, 2012; Jensen et al. 2012; Ben-Jaffel & Ballester 2013) during the transit of the hot Jupiter. The superposition of several X-ray observations also indicates an excess transit depth in this system (Poppenhaeger et al. 2013). In WASP-12 three observations confirmed the existence of an expanded atmosphere (Fossati et al. 2010; Haswell et al. 2012; Nichols et al. 2015). In GJ 436 b a large cometary tail like hydrogen cloud was detected, absorbing as much as 50% of the Lyα emission of the host star in three individual transit observations (Kulow et al. 2014; Ehrenreich et al. 2015). In addition, a tentative detection of excess absorption was obtained in the grazing transit of the hot Jupiter 55 Cancri b, while in the same system the transit of the super-Earth 55 Cancri e yielded a non-detection of atmospheric excess absorption in the Lyα line (Ehrenreich et al. 2012).

2. Approach

Although observations show that an extended atmosphere exists in at least the four planets HD 209458 b, HD 189733 b, WASP-12 b, and GJ 436 b there is no generally accepted theory for the observed absorption signals. The planetary atmospheres are affected by photoevaporation, by interactions with the stellar wind, by radiation pressure, and by planetary magnetic fields. All four processes also affect the absorption signals, for example, the observed signal of HD 209458 b has been explained by several models based on different physical assumptions (Koskinen et al. 2013b; Bourrier & Lecavelier des Etangs 2013; Tremblin & Chiang 2013; Trammell et al. 2014; Kislyakova et al. 2014). In our view, one of the most promising ways to disentangle the different influences is to observe the planetary atmospheres in different systems with a wide parameter range (Fossati et al. 2015). Since all the processes depend on different parameters (e.g., Lyα emission line strength, stellar wind density and velocity, magnetic field strength) the comparison of absorption signals from different systems can reveal the dominating processes in the atmospheres.

In this context, we identified the most promising targets in the search for yet undetected expanded atmospheres. Here, we present coupled radiative-hydrodynamical simulations of the planetary winds in 18 of these systems. Our simulations were performed with the PLUTO-CLOUDY interface (Salz et al. 2015a). In contrast to simple energy-limited estimates, our simulations identify the best targets for future transit observations with high certainty. The spherically symmetric simulations are comparable to several wind simulations in previous studies (e.g., Yelle 2004; Tian et al. 2005; García Muñoz 2007; Penz et al. 2008; Murray-Clay et al. 2009; Koskinen et al. 2013a; Shaikhislamov et al. 2014), but for the first time we applied them to possible future targets. Furthermore, our use of a detailed photoionization and plasma simulation code is the most important improvement compared to previous simulations. Additionally, we used a detailed reconstruction of the planetary irradiation levels in four spectral ranges for our simulations.

We first introduce our code and verify our assumptions (Sect. 3) and then compare the results for the different systems, explaining in detail the atmospheric structures (Sect. 4). We explain the transition from strong planetary winds in small planets to hydrodynamically stable thermospheres of compact, massive planets. Finally, we compute the Lyα absorption and emission signals from the simulated atmospheres and investigate their detectability (Sect. 5).

2.1. Target selection

Our target selection has been described in Salz et al. (2015b). Basically we select targets where expanded atmospheres can be detected by Lyα transit spectroscopy. We predict the Lyα luminosity of all host stars based on their X-ray luminosity, and if not available, based on spectral type and stellar rotation period using the relations presented by Linsky et al. (2013). We further only select transiting systems with an orbital separation <0.1 au to ensure high levels of irradiation. Note that this is not necessarily a boundary for hydrodynamic escape. Earth-sized planets can harbor escaping atmospheres at even smaller irradiation levels, because their atmospheres are relatively weakly bound as a result of the low planetary mass; this was already studied by Watson et al. (1981) for Earth and Venus. While we selected planets with a transit depth >0.5% in Salz et al. (2015b), we drop this criterion here, because recent observations indicate that also small planets can produce large absorption signals (Kulow et al. 2014). This extends our sample further toward planets with smaller masses and increases the parameter space.

Finally we limit the sample by selecting planets with a Lyα flux stronger than 1/5 of the reconstructed flux of HD 209458 b. Considering the remaining uncertainties of the scaling relations and the interstellar absorption, these targets are potentially bright enough for transit spectroscopy.

In total we find 18 suitable targets (see Tables 1 and 3), all located within 120 pc, since the stellar distance has a strong impact on the detectability. This probably places the targets within the Local Bubble of hot ionized interstellar material (Redfield & Linsky 2008), limiting the uncertainty introduced by interstellar absorption. As expected, systems with detected absorption signals (HD 189733, HD 209458, GJ 436) also lead the ranking according to our estimates. 55 Cnc e is also among the top targets, but for this planet an expanded atmosphere was not detected (Ehrenreich et al. 2012). We add three further targets to the list, which do not strictly fulfill our selection criteria: WASP-12 is too distant for Lyα transit spectroscopy, but excess absorption was measured in metal lines (Fossati et al. 2010). The distance of CoRoT-2 is too large, but the host star is extremely active (Schröter et al. 2011), which makes the system interesting. 55 Cnc b does not transit its host star, but a grazing transit of an expanded atmosphere has been proposed (Ehrenreich et al. 2012).

This finally leaves us with 21 planets in 20 systems. We present simulations for 18 of these targets. We did not simulate the atmospheres of WASP-38 b and WASP-18 b, because they host stable thermospheres (see Sect. 4.6). 55 Cnc b was not simulated, because the radius is unknown in this non-transiting planet. Planets similar to 55 Cnc b show densities ranging from 0.1 to 2.4 g cm^-3 (WASP-59 b, Kepler-435 b; Hébrard et al. 2013; Almenara et al. 2015). Reasonable values for the radius of 55 Cnc b range from 0.8 to 2.0 Jupiter-radii, which has a major impact on the planets atmosphere. Therefore, this system should be studied with a simulation grid in a dedicated study.

2.2. Reconstruction of the XUV flux and the SED

One of the main input parameters for our simulations is the irradiation strength and the spectral energy distribution (SED) of the host stars. Therefore, we first explain our reconstruction of the stellar SED, before going into details about the simulations.

Clearly, the strength of a planetary wind crucially depends on the available energy, which is mostly supplied by the EUV emission of the host stars, but X-rays also contribute a large percentage of the high-energy emission of active host stars like CoRoT-2 (Schröter et al. 2011). Interstellar absorption extinguishes EUV emission even from close-by targets so that observations cannot be used to determine the irradiation level in this spectral range. The upper atmosphere, which is simulated here, is mostly translucent to FUV, NUV, and optical emission, except for line absorption as shown by the observations (see Sect. 1). The simulations presented here only include hydrogen and helium (H+He), and any absorption beyond the Lyα line is negligible in terms of radiative heating. At the bottom of Fig. 1, we anticipate the spectrum of radiation that has passed through the simulated atmosphere of HD 209458 b to justify this statement. This radiation would be absorbed in lower atmospheric layers or close to the planetary photosphere. The atmospheric absorption from 10 to 912 Å is clear as well as the transmission of the optical emission, and furthermore the complete extinction of the Lyα line.

The reconstruction of the SED must be sufficiently accurate in the spectral range up to and including the Lyα line. For this purpose, we developed a piecewise reconstruction, which is mainly based on the prediction of the EUV luminosity of late type dwarfs by Linsky et al. (2014). Note that several methods exist to predict the EUV luminosity of stars (Salz et al. 2015b). The results of these methods differ by up to one order of magnitude for active stars, which directly translates into the uncertainty of the mass-loss rates in the simulations (see Sect. 3.8).

The reconstruction of the SEDs is split into four spectral regions and for each range we need a spectral shape and a luminosity to assemble the full SED:

• 1.

  For the X-ray spectrum (0 to100 Å) we use a 2 or4 × 10⁶ K plasma emission model from CHIANTI (Dere et al. 1997, 2009) for inactive or active host stars respectively (active: L_X> 10²⁸ erg s^-1). The SED is normalized to the observed X-ray luminosities of the host stars. Only if observations are unavailable, we use an estimate based on the stellar rotation period and the stellar mass (Pizzolato et al. 2003). The models are computed with a low resolution (bin width =10 Å) adapted to the resolution in the EUV range. Since the absorption of X-rays causes ionizations in the planetary atmospheres, resolving emission lines is irrelevant in the X-ray regime and the low resolution is sufficient.

• 2.

  For the hydrogen Lyα emission line of the host stars a Gaussian with a full width at half maximum (FWHM) of 9.4 Å is used. The irradiation strength is normalized according to the host star’s total Lyα luminosity. Only for four host stars has the luminosity been reconstructed based on HST observations (see references in Table 3). For most of the targets the Lyα luminosity is predicted based on the X-ray luminosity following Linsky et al. (2013). The line width is adapted to the resolution of our photoionization solver, which is 6.1 Å at the Lyα line. Our approach places about 90% of the line flux in the spectral bin at the Lyα line center. The line width is about a factor of 10 broader than usual Lyα stellar emission lines (compare Wood et al. 2005). Based on test simulations, we confirm that this only affects the strength of the radiation pressure significantly (see Sect. 4.8).

• 3.

  For the EUV range (100 Å to ~912 Å) the luminosity of the host stars is predicted based on the Lyα luminosity given by step 2 (Linsky et al. 2014). For most targets this procedure reverts back to an X-ray based estimate, because the Lyα luminosity is based on the X-ray luminosity. If a measurement of the Lyα luminosity is available, we obtain two predictions for the EUV luminosities (X-ray and Lyα based) and use the average. The shape of the SED is taken from an active or inactive solar spectrum (Woods & Rottman 2002), depending on the activity of the host star defined in step 1. This shape is continued to a connecting point with the photospheric blackbody, where we use a visual best fit for the connection point between 1500 to 4000 Å. Certainly the SED of the host stars will differ from the solar type emission in the EUV range, however, only the relative strength of broader ranges is relevant, because the radiation causes ionizations and the presence of spectral lines is unimportant (except for the He ii Lyα see below).

• 4.

  For the remaining part of the spectrum up to 5 × 10⁴ Å we choose a blackbody according to the host star’s effective temperature. This range does not affect the simulations, but can be used to check absorption in the hydrogen Balmer lines for example.

Our approach does not include a prediction of the He ii Lyα line (304 Å), which accounts for 20% of the EUV emission in the solar spectrum. The use of the solar SED shape introduces another error at this point (Linsky, priv. comm.). While a more detailed reconstruction of the SEDs is conceivable, it certainly has a smaller effect on the results than the uncertainty on the overall EUV luminosity of the host stars.

Figure 1 shows the reconstructed SEDs for the host stars and the derived luminosities for each spectral range are given in Table 3 alongside the results from our simulations. The total XUV irradiation is a factor of 200 stronger for CoRoT-2 b than in GJ 436 b. In these two systems X-rays contribute 60% and 30% of the total energy of hydrogen ionizing radiation. Also the relative Lyα emission line strength varies strongly for the individual targets. HD 209458 shows powerful Lyα emission given the inactive state of the host star. In contrast, the upper limit for GJ 1214 derived by France et al. (2013) is extremely low.

3. Simulations

We simulate the hydrodynamically expanding atmospheres of 18 hot gas planets in spherically symmetric, 1D simulations. Absorption and emission, and the corresponding heating and cooling is solved self-consistently using a photoionization and plasma simulation code, which incorporates most elementary processes from first principles. The computational effort of the presented simulations corresponds approximately to 300 000 h on a standard 1 GHz CPU.

The structure of this section is as follows. We introduce the code (Sect. 3.1) used for the simulations and explain the numerical setup (Sect. 3.2). We verify the convergence of the simulations (Sect. 3.3). Resolution and boundary effects are investigated in Sects. 3.4 to 3.6. The impact of metals and molecules on the results is checked in Sect. 3.7. We conclude this section with an overview of the uncertainties that generally affect 1D simulations such as ours both from numerical and observational points of view (Sect. 3.8).

3.1. Code

We use our photoionization-hydrodynamics simulation code the PLUTO-CLOUDY interface (TPCI) for the simulations of the escaping planetary atmospheres (Salz et al. 2015a). The interface couples the photoionization and plasma simulation code CLOUDY to the magnetohydrodynamics code PLUTO. This allows simulations of steady-state photoevaporative flows. PLUTO is a versatile MHD code that includes different physical phenomena, among others, gravitational acceleration and thermal conduction (Mignone et al. 2007, 2012). Our interface introduces radiative heating and cooling, as well as acceleration caused by radiation pressure from CLOUDY in PLUTO (Salz et al. 2015a). CLOUDY solves the microphysical state of a gas under a given irradiation in a static density structure (Ferland et al. 1998, 2013). The user can choose to include all elements from hydrogen to zinc, and CLOUDY then solves the equilibrium state regarding the degree of ionization, the chemical state, and the level populations in model atoms. Radiative transfer is approximated by the escape probability mechanism (Castor 1970; Elitzur 1982). The impact of the advection of species on the equilibrium state is included via a CLOUDY internal steady state solver as described in Salz et al. (2015a).

Recently Shematovich et al. (2014) have found that photoelectrons and suprathermal electron populations contribute significantly to excitation and ionization rates in the atmospheres of hot Jupiters. CLOUDY includes these processes assuming local energy deposition (Ferland et al. 1998, 2013), a transport of photoelectrons is not included. In partially ionized gases, this assumption is valid as long as the column densities of the individual cells exceed 10¹⁷ cm^-2 (Dalgarno et al. 1999), which is approximately fulfilled in neutral parts of the presented atmospheres. For highly ionized plasmas distant collision of the ionized species reduce the mean free path of photoelectrons by more than one order of magnitude (Spitzer 1962). The mean free path of photoelectrons becomes larger than one cell width only in the upper thermosphere of CoRoT-2 b and the planets which have stable thermospheres at a height where the atmospheres also become collisionless for neutral hydrogen (see Sect. 4.6, using equations from Goedbloed & Poedts 2004).

The use of CLOUDY is one of the main improvements compared with previous models of escaping planetary atmospheres because it solves not only the absorption of ionizing radiation in detail but also the emission of the gas. Most previous studies have either chosen a fixed heating efficiency for the absorbed radiation to account for a general radiative cooling effect or they specifically included individual radiative cooling agents like Lyα cooling (e.g., Murray-Clay et al. 2009; Koskinen et al. 2013a; Shaikhislamov et al. 2014). With such an approach it is impossible to find radiative equilibrium, i.e., a state where the complete radiative energy input is re-emitted, because not all relevant microphysical processes are included. Therefore, these simulations can only be used for strongly escaping atmospheres, where radiative cooling has a minor impact, but not for almost stable thermospheres, which are close to radiative equilibrium. TPCI solves this issue and can be used to simulate the transition from rapidly escaping atmospheres, where the radiative energy input is virtually completely used to drive the hydrodynamic wind, all the way to a situation with a hydrostatic atmospheres, where the absorbed energy is re-emitted.

3.2. Simulation setup

We simulate the expanding planetary atmospheres on a spherical, 1D grid. A similar setup is chosen for all simulations. The planetary atmosphere is irradiated from the top. We simulate the substellar point, which results in a maximal outflow rate because of the high irradiation level at this point. Forces due to the effective gravitation potential in the rotating two-body system are included. The resulting force is a superposition of the planetary and stellar gravitational forces with the centrifugal force caused by the planetary orbital motion. Furthermore, radiation pressure is included and thermal conduction is considered in a simplified manner by adding the thermal conductivity coefficients of an electron plasma and of a neutral hydrogen gas. A detailed explanation of the numerical implementation of these terms is given in Salz et al. (2015a).

Our simulations have several input parameters, which are fixed by observations, i.e., the planetary mass and radius, the semimajor axis, and the stellar mass. Additionally, we need to know the irradiating SED (see Sect. 2.2). The last two required parameters are not well defined: first, the temperature in the lower atmosphere, which is close to the equilibrium temperature but with a considerable uncertainty (see Sect. 3.5); and second, the boundary density, which reflects how deep our simulations proceed into the lower atmosphere (see Sect. 3.6). Fortunately, these two parameters have no strong impact on the outcome, therefore, the resulting atmosphere is unique for every planet and we do not need to sample a parameter space.

The simulated atmospheres consist of hydrogen and helium only; metals or molecules are neglected in our current setup (see Sect. 3.7). CLOUDY solves a plane parallel atmosphere, allowing radiation to escape in both directions toward the host star and toward the lower atmosphere. We double the optical depth at the bottom of the simulated atmosphere to ensure that radiation does not escape freely into a vacuum in this direction. We further use a constant microturbulence of 1 km s^-1, which is adapted to the conditions close to our lower boundary; the microturbulence has little impact on the overall atmospheric structure.

A stretched grid with 500 grid cells spans the range from 1 to 12/15 planetary radii. The resolution is increased from 2 × 10^-4R_pl to 0.22R_pl at the top boundary. We emphasize the high resolution at the lower boundary, which is not always visible in figures showing the complete atmospheres. The photoionization solver runs on an independent grid, which is autonomously chosen and varies throughout the simulation progress. It uses 600 to 1200 grid cells. The resolution in the dense atmosphere at the lower boundary can be up to a factor ten higher than in the hydrodynamic solver.

The number of boundary conditions that are fixed are two at the bottom (density, pressure) and none at the top boundary. This is the correct choice for a subsonic inflow and a supersonic outflow (e.g. Aluru et al. 1995). The bottom boundary condition is anchored in the lower atmosphere at a number density of 10¹⁴ cm^-3. This is sufficient to resolve the absorption of the EUV emission, which is mostly absorbed in atmospheric layers with a density of around 10¹¹ cm^-3. The density at the planetary photosphere is usually argued to be around 10¹⁷ to 10¹⁸ cm^-3, which we only state for a comparison (Lecavelier des Etangs et al. 2008a,b). Furthermore, we fix the pressure according to the equilibrium temperature of the planets (see Table 1). This results in pressures around 14 dyn cm^-2 equivalent to 14 μbar. The velocity in the boundary cell is adapted at each time step according to the velocity in the first grid cell. However, we only allow inflow, setting the boundary velocity to zero otherwise. This choice dampens oscillations that appear in simulations with steep density gradients at the lower boundary. The impact of the lower boundary conditions on the outcome of the simulations is discussed in Sects. 3.5 and 3.6.

The chosen numerical setup of the PLUTO code is a third order Runge-Kutta integrator (Gottlieb & Shu 1996), with the weighted essentially non-oscillatory finite difference scheme (WENO3, Jiang & Shu 1996) for interpolation and the Harten, Lax and van Leer approximate Riemann solver with contact discontinuity (HLLC, Toro et al. 1994). The radiative source term is not included in the Runge-Kutta integrator, but a simple forward Euler time stepping is used.

3.3. Convergence and initial conditions

All simulations were started from the same initial conditions, which are adjusted to an atmosphere with a medium mass-loss rate. Every atmosphere undergoes a series of shock waves in the first few hours of real time. This is depicted in Fig. 2 for the density structure in the atmosphere of WASP-80 b. Two shock waves propagate from the lower atmosphere (lefthand side) into the upper atmosphere (righthand side), increasing the average density in the thermosphere by a factor of ten. The depicted 7 h evolution in the atmosphere corresponds to less than 10^-3 of the final converging time. The final steady-state does not depend on the chosen initial conditions.

The convergence of the simulations was analyzed by checking the mass-loss rate and the temperature structure. We define a system specific time scale as (1)and choose v_pl = 1 km s^-1 as a typical velocity scale in the middle of the planetary thermospheres. 1000 of such time steps correspond to between half a year and four years in the systems. This time is sufficient for atmospheric material to propagate at least once through the computational domain. We ensured that the temperature structure and the total mass-loss rate have converged to a level of 0.1 on this time scale (ΔṀ/Ṁ< 0.1). Individual simulations were advanced a factor of 10 longer and did not show any significant further changes.

3.4. Resolution effects

We have checked that the chosen resolution does not affect the results of our simulations. The resolution of the stretched grid at the lower boundary has been increased in three steps by factors of two in the simulations of HD 209458 b and HD 189733 b. HD 189733 b has a steeper density gradient in the lower atmosphere and is more challenging for the numerical solver. Decreasing the resolution leads to larger oscillations in the atmospheres, which remain after the convergence of the simulations as defined above. The remaining oscillations are larger in the atmosphere of HD 189733 b than in HD 209458 b, hence, atmospheres with lower mass-loss rates are affected more strongly by such oscillations. They do not affect the time averaged density, temperature, or velocity structure and, therefore, also not the mass-loss rate (fractional change of the mass-loss rate ΔṀ/Ṁ< 10^-2).

The oscillations are numerical artifacts with two sources. First, the small acceleration of the planetary wind in the lower atmospheres results from a slight imbalance of the gravitational force and the force due to the hydrostatic pressure gradient. Subtraction of these large forces results in a small value, which is prone to numerical inaccuracy. Increasing the resolution results in better estimates of the pressure gradient, which in turn reduces the errors and thus the oscillations. The second source of oscillations in the lower atmosphere of planets with a small mass-loss rate is a consequence of these atmospheric layers being almost in radiative equilibrium. When the microphysics solver is called within TPCI the best level of equilibrium between radiative heating and cooling is about 10^-2. This results in small transient heating and cooling events in regions of radiative equilibrium, which then induce oscillations.

In some cases the oscillations are large enough to be seen in our figures. If this is the case, we present time averaged structures. The convergence to the steady-state as defined above is unaffected by these oscillations, but the accuracy of the final atmospheric structure is worse than in simulations without oscillations. This is one reason for our conservative convergence level of 10%. Furthermore, the computational effort of the individual simulation is increased by the oscillations.

The independence of the mass-loss rates on the resolution distinguishes our simulations from those of Tian et al. (2005). The mass-loss rates in their simulations depended on the chosen number of grid points. However, the authors used a Lax-Friedrich solver with a higher numerical diffusion compared with the more advanced Godunov-type scheme with the approximate HLLC Riemann solver. Tian et al. identified the numerical diffusion to be responsible for the dependency of their mass-loss rates on the resolution, which is consistent with our findings.

3.5. Boundary temperature

The temperature at the lower boundary is one of the free parameters in our simulations. It corresponds to a temperature in the lower atmosphere of the planet at a certain height. Basically, the equilibrium temperature T_eq of the planet provides a simple approximation for this temperature (with zero albedo, e.g., Charbonneau et al. 2005): (2)Here, T_eff is the effective temperature of the host star, R_st is the stellar radius, and a is the semimajor axis of planetary orbit. In some cases measurements of the average dayside brightness temperatures are available from infrared transit observations (see references of Table 1). However, neither the equilibrium temperature nor the brightness temperature necessarily represents the temperature at the atmospheric height, which corresponds to the lower boundary of our simulations. Detailed radiative transfer models are another method to obtain an accurate temperature structure for the lower atmosphere, but such models only exist for well studied systems like HD 209458 b. For HD 209458 b and HD 189733 b pressure-temperature profiles were reconstructed for the atmosphere above the terminator based on sodium absorption measurements (Vidal-Madjar et al. 2011b; Huitson et al. 2012; Wyttenbach et al. 2015), but these probably do not apply to the substellar point. Even for the probably best studied case, HD 209458 b, Koskinen et al. (2013a) conclude that despite several different observations and models the temperature structure in the lower atmosphere is uncertain.

Fortunately, we do not need to know the temperature at the lower boundary with high precision. To verify this, we tested the impact of this temperature on our results. In the simulations of the atmospheres of HD 209458 b and HD 189733 b, which have boundary temperatures of 1320 K and 1200 K respectively, we increased and decreased the temperature in steps of 100 K/200 K. Figure 3 shows that increasing the boundary temperature slightly expands the planetary atmosphere, but the large scale structure remains unaffected as does the temperature structure in the upper planetary thermosphere. The change in the mass-loss rate scales linearly with the change of the boundary temperature over the tested range. The overall impact is small: A change of 40% in the boundary temperature causes the mass-loss rate to change by less than 5%. Our results are comparable to those of Murray-Clay et al. (2009), who tested a wider range for the lower boundary temperature. This result is not surprising, because the gain in advected energy due to a 1000 K temperature change is about a factor ten smaller than the total radiative energy input. The change in the mass-loss rate is about twice as strong in the simulation of HD 189733 b. A temperature minimum in the atmosphere of HD 209458 b decouples the simulation better from the boundary conditions, because radiative processes have more time to dispose of additional energy.

Even considerable changes in the lower boundary temperature affect the total mass-loss rates in our simulations by less than 10%. With respect to the uncertainty on other parameters, especially the irradiation level, we can neglect the effect at this point.

3.6. Boundary density

The density at the lower boundary is another free parameter in our simulations. Changing the boundary density has a larger effect on the mass-loss rate than changing the boundary temperature. This is a result of a slight inconsistency in the setup of our simulations. The density at the lower boundary is a density that is reached at a certain height above the planetary photosphere³. However, our simulation setup always starts at 1R_pl. If we start with a density two orders of magnitude smaller, this layer is located higher up in the planetary atmosphere; the change in height of the maximal volume heating rate is about 0.025R_pl in the case of HD 209458 b. In the spherical simulation the irradiated area of the grid cells increases with , thus in this example the similar volume heating rate results in a 5% increase of the total heating rate, which in turn affects the mass-loss rate. This effect is depicted in Fig. 4, which shows that the density structures of the individual simulations almost merge in the upper thermosphere of HD 209458 b, but larger differences occur in the lower atmosphere. These would be mitigated if we had started the simulations with a lower boundary density at the correct heights as depicted in the insert.

As a secondary effect, the temperature minimum in the lower atmosphere becomes deeper as the mass loss increases. In this region of the atmosphere, radiative heating is small and the increased adiabatic cooling caused by a higher mass-loss rate reduces the local temperature (Watson et al. 1981; García Muñoz 2007). At even lower boundary densities (10¹² cm^-3) a significant fraction of the XUV radiation passes through our atmosphere without being absorbed in the computational domain, which further decreases the mass-loss rate.

For HD 189733 b the impact of the boundary density on the mass-loss rate is larger, because the steeper atmosphere shifts stronger in response to a change of the boundary density. Decreasing the boundary density to 10¹² cm^-3 causes a discontinuity in the density structure from the boundary cell to the first grid point, and the total mass-loss rate is reduced by 60%. We confirmed that none of the presented simulations is affected by such a density discontinuity.

Our standard simulation has a density of 10¹⁴ cm^-3 at the lower boundary. The uncertainty in the correct height of our boundary layer above the planetary photosphere introduces an error on the mass-loss rates of less than 50%, consistent with the results of Murray-Clay et al. (2009). We further note that the observationally determined planetary radii are usually not accurate to the percent level, so that this error is inevitable.

3.7. Metals and molecules

Our current simulations include the chemistry of hydrogen and helium. Adding metals substantially increases the computational load and will be studied for individual systems in the future. To justify that we can neglect metals in this study, we anticipate results from a test simulation of the atmosphere of HD 209458 b including all 30 lightest metals with solar abundances. We find that the temperature structure of the lower atmosphere (< 1.05R_pl) is dominated by line heating and cooling of metals, but the temperature in this region only has a small effect on the mass-loss rate (see Sect. 3.5). Cooling of Ca ii and Fe ii dominates from 1.1 to 1.6R_pl and reduces the maximum temperature in the atmosphere by about 1000 K. Additionally, Mg ii and Si ii are important cooling agents. While such heavy metals can be advected by the strong planetary wind (Yelle et al. 2008; García Muñoz 2007), their occurrence in the upper thermosphere depends on the formation of condensates in the lower atmosphere (e.g., Sudarsky et al. 2000). Including metals reduces the mass-loss rate in our test simulation by a factor of two. While this value may be used for guidance, the result certainly depends on the metal abundance, because increasing the abundances to supersolar values will cool the planetary atmosphere more efficiently. Hence, metals have an effect on the results, but as long as they remain minor constituents their impact likely does not exceed the uncertainty in the EUV irradiation level (Salz et al. 2015b).

Metals also efficiently absorb X-rays, which increases the available energy. Strongly expanded thermospheres like that of HD 209458 b are opaque to X-rays longward of 25 Å, but more compact atmosphere only are opaque longward of 60 Å in our simulations (defined as 50% absorption). Including metals shifts these cutoff values to approximately 10 and 20 Å. In the case of our most active host star CoRoT-2 metals would increase the absorbed energy in the planetary atmosphere by 13%.

We also disable the formation of molecules in CLOUDY. This choice is unphysical, because photodissociation of H₂ would occur close to the lower boundary (Yelle 2004; García Muñoz 2007; Koskinen et al. 2013a). Molecules strongly affect the temperature structure of the lower atmospheres, but the impact of this temperature structure on the mass-loss rate is quite small (see Sect. 3.5). Furthermore, planets with a moderate irradiation level can host a molecular outflow. A detailed study of this phenomenon is beyond the scope of this paper, but we performed a test simulation of one of the planets with the smallest irradiation level, GJ 1214 b, including molecules in the planetary atmosphere. The simulation shows an H₂ fraction of 10 to 25% throughout the thermosphere of the planet. H⁻ is a significant cooling agent, as for example also seen in the solar chromosphere at similar conditions (Vernazza et al. 1981). The mass-loss rate is reduced by 15%, thus, our simulations without molecules remain valid as long as hydrogen is the major constituent of the atmospheres.

3.8. Overview of the uncertainty budget

One-dimensional simulations of the planetary winds are subject to several uncertainties. These are introduced by observational limitations and by the simplifications of the 1D model. We provide a short overview and estimate the impact of probably the most important effects that are neglected (see Table 2).

A large uncertainty is introduced by neglecting the 3D structure of the system geometry, however, Stone & Proga (2009) have shown that 1D models produce valid estimates for the planetary mass-loss rates, if the irradiation is averaged correctly over the planetary surface. In Salz et al. (2015b) we indicated that the irradiation level of hot gas planets is uncertain, because the EUV spectral range is completely absorbed by interstellar hydrogen and different reconstruction methods differ by up to a factor of 10. Most recently, Chadney et al. (2015) revised one of the reconstruction methods, reducing the differences for active stars. Nevertheless, we adopt a conservative error estimate of a factor of 3 for the EUV reconstruction. The impact of metals and molecules on the wind has been explained in the previous section. Related to our simulations, the boundary conditions introduce a small uncertainty on the resulting mass-loss rates.

1D-HD simulations of the planetary winds ignore any interactions with possible stellar winds. This approximation is valid if the interaction with the stellar wind occurs after the planetary wind becomes supersonic. However, Murray-Clay et al. (2009) showed that the stellar wind can suppress the formation of a planetary wind at least on the dayside of the planet. There are likely two regimes of confinement: first, the planetary wind becomes too weak to sustain a supersonic expansion and a subsonic planetary wind persists with a lower mass-loss rate. Second, if the stellar wind is strong enough, the planetary wind is completely suppressed.

Furthermore, planetary winds can interact with planetary magnetic fields. This has been studied for example by Trammell et al. (2014), who stated that, if the planet’s magnetic field is strong enough, the outflow can be suppressed over larger fractions of the planetary surface. The authors also found that there are always open field lines along which a planetary wind can develop freely.

The formation of planetary winds is a complex physical system, which currently can only be studied by neglecting certain aspects. We focus on accurately solving the energy balance throughout the thermospheres regarding the absorption and emission of radiation by the atmospheric gas in detail. In the following sections we also show atmospheric structures above the planetary Roche lobe height, which is only for guidance. Our spherical approximation is invalid above this height and interactions with stellar winds can deform and accelerate the planetary wind material.

4. Results and discussion of the simulations

We now present the simulated photoevaporative winds of 18 exoplanets. The atmospheres of HD 209458 b and HD 189733 b are compared in Sect. 4.2, where we explain in detail the individual processes affecting the planetary atmospheres. The complete sample of planets is presented and compared in Sect. 4.3, where the impact of individual system parameters on the planetary atmospheres is investigated. Here, we also explain the approach to radiative equilibrium in the thermospheres of massive and compact planets. We then compute the mass-loss rates (Sect. 4.4), explain the case of WASP-12 b (Sect. 4.5), identify stable thermospheres (Sect. 4.6), consider the ram pressure of the planetary winds (Sect. 4.7), and estimate the strength of the radiation pressure (Sect. 4.8).

4.1. Heating fraction and heating efficiency

For the analysis of our simulations we define the height-dependent heating fraction as (3)Note that the heating fraction differs from the often used heating efficiency, which is defined as: (4)Here, the denominator holds the total energy contained in the absorbed radiative flux. The difference of the two definitions is mainly the energy that is used to ionize hydrogen and helium; this energy fraction is not available for heating the planetary atmosphere. Equation (3) excludes this fraction of the radiative energy input, hence, the heating fraction approaches one as radiative cooling becomes insignificant, and it is zero in radiative equilibrium. The heating efficiency is also zero in radiative equilibrium, but it never reaches a value of one. The heating efficiency is not directly available in our simulations.

4.2. Structure of the escaping atmospheres

To understand the general behavior of escaping planetary atmospheres, we compare two examples: HD 209458 b and HD 189733 b. For HD 209458 b we also juxtapose a simulation with only hydrogen as constituent (H-only). In Fig. 5 we show the density, temperature, velocity, and specific heating rate throughout the atmospheres in the three simulations. Panel e further shows individual heating and cooling terms in the atmosphere of HD 209458 b with hydrogen and helium (H+He). Finally, Fig. 6 splits the radiative heating and cooling terms into the most important agents. In all our figures the righthand side of the plot is the top of the atmosphere, while the lefthand side is the lower boundary, which is located in the denser atmosphere close to the planetary photosphere.

The general structure of all planetary atmospheres is similar. The density in the atmospheres decreases with height; the density gradient in the thermosphere is shallow due to the high temperature, which results from the absorption of stellar XUV emission. The high temperature further leads to a persistent expansion of the atmospheric gas, which accelerates the atmosphere against the gravitational force of the planet. The initial velocity is small (~ 10^-5 km s^-1) at the lower boundary, but reaches the sonic point within the simulation box, and the atmospheric material exits the simulated region as a supersonic wind. Thus, we simulate the transonic hydrodynamic planetary winds. At the Roche lobe height, all simulated atmospheres reach velocities between 10 and 20 km s^-1.

4.2.1. Atmosphere of HD 209458 b

First, we focus on the atmosphere of HD 209458 b, which is depicted by the green solid lines in Fig. 5. Panel a shows how the density decreases from the boundary density of 10¹⁴ to 10⁵ cm^-3. The extent of the Roche lobe for HD 209458 b is 4.2R_pl above the substellar point and the sonic point is reached at 3.5R_pl. The optical depth for EUV photons (60 eV) reaches one at a height of 1.03R_pl. This height corresponds well with the atmospheric layer that experiences the maximum volume heating rate.

Panel b in Fig. 5 shows that the absorption of XUV emission around the height of 1.03R_pl increases the temperature. This raises the atmospheric scale height and the density gradient becomes more shallow above this height. The temperature reaches a maximum of 9100 K. The bulk velocity strongly increases along with the temperature in the lower atmosphere (see panel c).

Panel d shows the net specific radiative heating rate throughout the atmosphere. We plot the specific heating rate, which is the heating rate per gram of material, because it clearly shows where the atmospheric material gains most of the energy for the escape. The net heating rate has a maximum close to the temperature maximum, but declines above 1.6R_pl where the atmosphere is more strongly ionized. There is a distinct dip in the net heating rate at 1.15R_pl, which is caused by hydrogen line cooling. This can be seen in panel e, where we show the radiative heating and cooling rates separately as thin, solid lines. Subtracting the radiative cooling rate from the heating rate gives the net radiative-heating rate (solid thick line), which is displayed in panels d and e. Between 1.05 and 1.8R_pl about 45% of the radiative heating produced by the absorption of XUV radiation is counterbalanced by hydrogen line emission, which is mainly Lyα cooling.

In panel e we further compare the radiative heating rate to the hydrodynamic sources and sinks of internal energy. In the atmosphere of HD 209458 b the advection of internal energy is a heat sink below 1.7R_pl, but acts as a heat source above this level, providing 50% of the energy in the upper thermosphere. The adiabatic expansion of the gas is a heat sink throughout the atmosphere. This term stands for the energy, which is converted into gravitational potential and kinetic energy. It now becomes clear why the temperature decreases in the upper thermosphere: Radiative heating becomes inefficient due to a high degree of ionization, but the adiabatic cooling term remains large. Therefore, the advected thermal energy is used to drive the expansion of the gas. Panel e also shows cooling due to thermal conduction, which is largest where the temperature increases most strongly, but it is not relevant in any of the simulations.

In Fig. 6 the H+He atmosphere of HD 209458 b is depicted in the middle column; panel a shows the fractional abundances. Ionized hydrogen is the major constituent above 1.6R_pl (H/H⁺ transition layer). This correlates with the height above which the radiative heating declines. The He/He⁺ transition occurs slightly below, but charge exchange between H and He⁺ couples the ionization heights of the two elements closely. Panel b shows the individual radiative heating agents. Ionization of hydrogen is the most important heating agent, but also ionization of helium contributes throughout the atmosphere. Ionization of He⁺ to He²⁺ contributes less than 10% of the heating rate in this atmosphere and is most important above the Roche lobe (not depicted in Fig. 6). Charge exchange between H and He⁺ is a heat source when hydrogen is weakly ionized.

Panel c of Fig. 6 shows the cooling agents and the heating fraction. Where the total radiative cooling rate is small (< 10% of the heating rate) the plot is shaded. The most important cooling agent is hydrogen line cooling between 1.06 and 1.7R_pl (see the discussion of Fig. 5e above). Beyond 1.8R_pl radiative cooling is small; the main cooling agents in this region are hydrogen recombination and free-free emission of electrons (thermal bremsstrahlung). Only at the bottom of the atmosphere (< 1.02R_pl) free-free emission is stronger than the heating rate. At this level, the absorption of Lyα radiation with following collisional de-excitation is the main local energy source, because ionizing radiation has been absorbed in higher atmospheric layers. The specific cooling at this height is negligible. Helium contributes less than 10% to the radiative cooling rate in the atmosphere of HD 209458 b.

While our simulation shows that HD 209458 b produces a strong, transonic wind as a result of the absorption of XUV radiation, it also demonstrates that the hot atmosphere re-emits a significant fraction of the radiative energy input mostly by Lyα line emission in intermediate thermospheric layers.

4.2.2. The H-only atmosphere and the H/H⁺ transition layer

If the atmosphere of HD 209458 b only consisted of hydrogen, the general structure would be similar to an atmosphere including hydrogen and helium (see Fig. 5). However, an H-only atmosphere of HD 209458 b produces a slightly stronger planetary wind. On average, including helium in our simulations decreases the mass-loss rate by a factor of two, which contradicts results from García Muñoz (2007). Since the involved Lyα cooling is important for all atmospheres, we will explain the difference here in more detail.

The change in density throughout the atmosphere is larger than the change in velocity when including helium, so let’s assume the velocity would remain constant. Including helium increases the mean molecular weight of the atmosphere, therefore, a higher temperature is needed to sustain the atmospheric density structure and the planetary mass-loss rate. A higher temperature, however, increases the cooling rate by increasing the collisional excitation of hydrogen. Lyα cooling is proportional to the occupation number of the second level; the relative occupation number can be approximated using the Boltzmann factor, exp(−hν/kT). Lyα radiation mainly escapes from a depth around 1.15R_pl, where the temperatures are about 7900 and 8400 K in the H-only and H+He simulations respectively. The higher temperature in the H+He simulation raises the occupation number of the second level by a factor of two, which results in the increased cooling rate (see Fig. 6c). The simulations of García Muñoz (2007) did not include Lyα cooling, hence, he could not find this effect.

The H/H⁺ transition in the H-only simulation occurs at a height of 2.3R_pl in contrast to 1.6R_pl in the H+He simulation. This higher transition layer results from the stronger mass-loss rate, because more neutral hydrogen is advected into the thermosphere. Additionally, the shift of the transition layer further reduces the Lyα cooling rate by increasing the neutral hydrogen column density above the Lyα emission layer.

Koskinen et al. (2013a) found the H/H⁺ transition to occur at 3.4R_pl in the atmosphere of HD 209458 b. The height of this layer depends on the cooling in the lower thermosphere as seen by the comparison of the H-only and H+He simulations. This was also noted by Koskinen et al.; they compared their results with those of Yelle (2004), who found strong H cooling and a transition layer height of 1.7R_pl. Koskinen et al. argued that H is probably not formed if metals are present in the atmosphere, but this argument does not hold for Lyα cooling. Unfortunately, the considerable uncertainties in this type of simulations (see Table 2) do not allow a final statement regarding the transition layer height in the atmosphere of HD 209458 b.

4.2.3. Impact of the gravitational potential on the mass loss: HD 189733 b versus HD 209458 b

Comparing the atmospheres of HD 189733 b and HD 209458 b, one finds three major differences: The atmosphere of HD 189733 b is hotter, the expansion velocity is smaller below the Roche lobe, and the net radiative-heating rate is also smaller (see Fig. 5), although the XUV irradiation level is actually 16 times higher. The maximum temperature in the atmosphere of HD 189733 b is 2700 K higher than in HD 209458 b, which is also a result of less adiabatic cooling of the weaker planetary wind. The atmosphere reaches the sonic point just below the Roche lobe height, above which the material is accelerated freely toward the stellar surface in the 1D simulation. The temperature decreases more strongly in the upper thermosphere, because the degree of ionization is high and only little radiative heating remains.

In the case of HD 189733 b it is remarkable that the simulated wind is weaker in terms of the total mass-loss rate although the irradiation level is higher than in HD 209458 b. The weaker wind is a result of increased radiative cooling, which actually depends on the depth of the gravitational potential well of the planet. This can be demonstrated with a thought experiment: let us start with a planet that hosts a strong wind and only weak radiative cooling. What happens if we increase the mass of the planet?

• The first reaction will be a contraction of the atmosphere, because the gravitational attraction was increased. Let us only slightly increase the planetary mass, so that the upward flux of the planetary wind does not collapse but only becomes weaker.

• The smaller mass flux has two effects: first, the adiabatic cooling decreases, because a smaller velocity means that the atmospheric expansion proceeds slower. This causes a higher atmospheric temperature. Second, the advection of neutral hydrogen into the upper atmosphere is weaker, which shifts the ionization front closer to the planetary surface.

• The increased temperature partially counterbalances the higher gravity and prevents the atmosphere to contract further, but it also increases radiative cooling (Lyα and free-free emission). Lyα cooling is further increased by the now smaller neutral hydrogen column density above the emission layer, which increases the escape probability for the radiation.

• The radiative cooling reduces the available energy for driving the planetary wind.

• The outflow settles at a lower mass-loss rate and the resulting atmosphere is hotter and more strongly ionized.

If the radiative cooling did not increase with the temperature, the down settling of the atmosphere in the first step would be a transient effect and the mass-loss rate would be restored to the original value. Of course, we can drive this thought experiment further by again increasing the mass of the hypothetical planet. At a certain mass the radiative energy input will be completely re-emitted by a hot and more compact thermosphere. This explains a transition in the atmospheres of hot exoplanets. Smaller planets host strong and cool planetary winds, and massive and compact planets host stable, hot, hydrostatic thermospheres.

The intermediate case is precisely what occurs in the atmosphere of HD 189733 b. Most of the EUV photons are absorbed at a height of 1.08R_pl. While the heating fraction is close to unity at this height in the atmosphere of HD 209458 b, in HD 189733 b it ranges from 5% to 20%, which results in only weak acceleration of the planetary wind. Hydrogen and helium are mainly ionized above 1.2R_pl and helium is doubly ionized above 2.1R_pl (see Fig. 6). The overall structure of the heating and cooling agents is similar to that in the atmosphere of HD 209458 b, but shifted toward smaller heights. Photoionization of He⁺ to He²⁺ is the main radiative heating agent above 2.0R_pl. However, the total radiative heating in the upper thermosphere (> 2.4R_pl) is small compared to the advected energy. The lack of radiative heating in the upper thermosphere of HD 189733 b makes it cooler than that of HD 209458 b.

4.2.4. Atmospheric structure at the thermobase

For HD 209458 b and HD 189733 b we can compare the temperature structure of the lower atmosphere with those derived from atmospheric absorption measurements in the sodium D lines (Vidal-Madjar et al. 2011b,a; Huitson et al. 2012). In contrast, the inferred inverted or non-inverted temperature profiles obtained from Spitzer data or recently based on high-dispersion spectroscopy in the infrared (e.g., CHRIRES on the Very Large Telescope) are sensitive to atmospheric layers below the 10^-5 bar level and, thus, below our simulated region (Line et al. 2014; Schwarz et al. 2015).

Our atmospheric temperature structures show a reasonable agreement with the Na D observations (see Fig. 7a). For the lower atmosphere of HD 209458 b, Vidal-Madjar et al. (2011b,a) found a temperature of 365 ± 200 K, which is below our adopted temperature of 1320 K derived from Spitzer observations (Crossfield et al. 2012). Decreasing the boundary temperature in our simulations slightly shifts the atmospheric profile toward smaller heights (see Fig. 3), so that we could certainly produce a better fit by adopting a lower boundary temperature. However, we simulated the substellar point at which the atmosphere is expected to be hotter than at the terminator, which is probed by the observations. The temperature minimum of 500 K in our simulated atmosphere of HD 209458 b agrees remarkably well with the low temperatures deduced from the observations.

While the observed temperature-pressure profile of HD 189733 b is reproduced by our simulations, the atmosphere of HD 209458 b shows more pronounced differences (see Fig. 7b). In our simulations the temperature rise at the thermobase occurs at lower pressure levels, thus higher in the planetary atmosphere than in the observations. This disagreement has been noted before for various other 1D escape models of this planetary atmosphere (Vidal-Madjar et al. 2011b,a; Huitson et al. 2012). Our test simulations including molecules show a potential source for this discrepancy, namely H⁻ dominates the atmospheric heating and cooling at the relevant heights (see Sect. 3.7). While these test simulations are not fully converged and require a more detailed analysis, preliminary results reproduce the observed temperature-pressure profiles more closely in both HD 209458 b and HD 189733 b.

4.3. The atmospheres of hot gas planets in comparison

We now compare the atmospheres of our sample of hot gas planets. HAT-P-2 b, HAT-P-20 b, WASP-8 b, and WASP-10 b are excluded, because these compact planets host stable thermospheres (see Sect. 4.6). In the following, we discuss how the atmospheres depend on system parameters such as irradiation level, planetary density, or the planetary gravitational potential. We either compare two planets, which differ mostly in a single system parameter, or we compare two subsamples, for which we divide our systems loosely at a gravitational potential of log − Φ_G = log GM_pl/R_pl = 13.0 in units of erg g^-1. The cut allows us to explain differences in the atmospheres of high-potential and low-potential planets, which we use as a short term for the strength of the gravitational potential per unit mass.

The atmospheres of the complete sample of planets are shown in Fig. 8, which is split into two columns: low-potential planets are depicted in the lefthand column and high-potential planets in the righthand column. In the figure each subsample is ordered with respect to the mass-loss rate. GJ 3470 b has the strongest mass-loss rate among the low-potential planets and WASP-12 b is its counterpart among the high-potential planets. HD 97658 b and CoRoT-2 b are the planets with the smallest mass-loss rates among the subsamples. The general behavior of all atmospheres follows our explanations in Sect. 4.2. The density decreases with height, the temperature reaches several thousand degrees maximally and decreases in the upper thermosphere, and all atmospheres reach a sonic point between 3 and 5R_pl, except for the atmosphere of WASP-12 b (see Sect. 4.5).

4.3.1. Atmospheric temperature and the H/H⁺ transition layer

A comparison of both columns in Fig. 8 shows that high-potential planets have a lower average density in the thermosphere. The temperature maximum always occurs between 1.3 and 3R_pl in the atmospheres, and rises from low-potential (2700 K in GJ 1214 b) to high-potential planets (16 800 K in CoRoT-2 b). The level of irradiation only weakly affects the maximal temperature in the atmospheres, e.g., the irradiation level of WASP-12 b is 16 times higher than in HD 209458 b, but the maximal temperature is only about 10% higher. Strong adiabatic cooling lowers the thermospheric temperature of WASP-12 b.

Super-Earth-sized planets have atmospheres that are almost completely neutral throughout our simulation box. As the planets become more massive and compact, the H/H⁺ transition layer moves closer to the planetary surface (see Fig. 8d), which is mainly a result of the higher temperature in the thermospheres of high-potential planets. To some degree the irradiation strength also affects the height of the transition layer, which can be demonstrated by comparing GJ 436 b with HAT-P-11 b. Both planets have similar densities and gravitational potentials. Despite the stronger planetary wind of HAT-P-11 b, the H/H⁺ transition occurs at a height of 3.7R_pl in its thermosphere compared to 7.8R_pl in GJ 436 b. Thus, the higher XUV irradiation level of HAT-P-11 b shifts the transition layer to smaller heights. The almost neutral atmospheres of GJ 1214 b and HD 97658 b result from the small gravitational potentials of the exoplanets in combination with moderate XUV irradiation levels. These two thermospheres would also contain a significant fraction of H₂, which is neglected here (see the discussion about molecules in Sect. 3.7).

The specific heating rate is more homogeneously distributed throughout the thermospheres of low-potential planets (see Fig. 8e), which is a result of a higher neutral hydrogen fraction in their atmospheres. In contrast, the upper thermospheres of WASP-77 b, WASP-43 b, and CoRoT-2 b experience negligible amounts of radiative heating, because both hydrogen and helium are almost completely ionized. For example, the fractional abundance of He²⁺ is > 0.99 above 2R_pl in CoRoT-2 b.

4.3.2. Approaching radiative equilibrium: the effect of the gravitational potential

The τ(60 eV) = 1 layer always occurs between 1.03 and 1.15R_pl for all planets⁴. The planetary wind is strongly accelerated at this height in low-potential planets, but planets with a higher potential (log − Φ_G> 13.2) show almost no acceleration in this atmospheric layer (see Fig. 8c). Their atmospheres are more gradually accelerated between 1.1 and 2R_pl. The small amount of acceleration in the lower thermospheres of high-potential planets is a result of strong radiative cooling. Consider, for example, CoRoT-2 b with the strongest irradiation among our sample. This planet experiences the smallest net radiative-heating rate. Below 1.3R_pl the atmosphere is close to radiative equilibrium (see Fig. 8f), so that Lyα and free-free emission re-emit almost the complete radiative energy input of 1.55 × 10⁵ erg cm^-2 s^-1. Radiative cooling of these two cooling agents becomes more and more important going from planets with a low gravitational potential to high-potential planets.

Lyα cooling starts to re-emit larger fractions of the energy input at a height around 1.2R_pl in the atmospheres of WASP-80 b and HD 149026 b (see panels f in Fig. 8). This cooling is small in the smallest planets, and only becomes important once the atmospheres are heated to more than 8000 K. For example, at a temperature of 5000 K the collisional excitation of the second level in hydrogen is reduced by a factor of 7000, thus, Lyα cooling is small. In high-potential planets the τ(60 eV) = 1 layer occurs close to the bottom of the Lyα emission layer. As a result, in HD 209458 b 80% of the absorbed energy delivered by these photons is re-emitted, in HD 189733 b 99% is re-emitted and in WASP-43 b this fraction increases to about 99.8%.

Radiation with higher energies penetrates deeper into the atmospheres and heats the atmosphere of HD 209458 b below the Lyα emission layer. Free-free emission, the second major cooling agent, is weak in the atmosphere of HD 209458 b, which is reflected by a heating fraction of almost one around 1.04R_pl. This leads to heating and acceleration of the planetary wind by stellar irradiation with photon energies higher than 60 eV. The heating below the Lyα emission layer becomes less important in the atmosphere of HD 189733 b and the heating fraction is reduced to 0.001 for WASP-43 b, which is a result of increased free-free emission in hot, lower atmospheric regions. This temperature increase of the dense atmospheric layers close to the lower boundary in massive and compact planets can be seen in panel b of Fig. 8.

4.4. Mass loss rates

In the simulations we evaluate the mass-loss rates using the formula Ṁ^sim = 1/4 × 4πR²ρv. The additional factor of 1/4 results from averaging the irradiation over the planetary surface. In our spherical simulation the planet is irradiated from all sides, but in reality only the dayside is illuminated. The difference of the spherical surface and the surface of the planetary dayside disk is 1/4. The factor is also supported by two-dimensional simulations (Stone & Proga 2009). A sorted list of the calculated mass-loss rates is provided in Table 3. We also provide two more mass loss columns: Ṁ is the fractional mass loss per Ga, and is an estimate for the total, fractional mass loss during the first 100 Ma after formation of the planetary system. This value is derived by linearly scaling the mass loss according to the expected irradiation from the young host stars. The XUV emission is derived in the same fashion as before, but assuming average X-ray luminosities for young stars between the T Tauri phase and the age of the Pleiades (100 Ma) from Stelzer & Neuhäuser (2001), i.e., G-type log L_X ≈ 29.7, K-type log L_X ≈ 29.5, and M-type log L_X ≈ 29.0.

In general, the current mass-loss rates of the planets in our sample are most strongly affected by the gravitational potential, the irradiation level and by the planetary density. The effect of the gravitational potential is directly visible in Fig. 8, looking at the low thermospheric densities of high-potential planets. The impact of the irradiation level can be demonstrated by coming back to our comparison of GJ 436 b with HAT-P-11 b. HAT-P-11 b experiences a five times higher irradiation level and its wind is 4.4 times stronger. The effect of the planetary density on the mass-loss rates can be seen by comparing HD 149026 b and WASP-80 b, which both have almost the same gravitational potential. The irradiation of HD 149026 b is almost twice as strong as in WASP-80 b, but the mass-loss rate is lower. This is a result of the higher planetary density.

Our simulations show that none of the planets is currently losing mass at a rate that would evaporate the complete planet over its lifetime. For WASP-12 b this contradicts earlier mass loss estimates (see Sect. 4.5). Among our sample 55 Cnc e is highly interesting, because the planetary density of 4 g cm^-3 suggests that the planet must consist of a significant amount of volatiles (~20%, Gillon et al. 2012). Our calculated mass-loss rate of 1.4 × 10¹⁰ g s^-1 is higher than the upper limit set by Ehrenreich et al. (2012) and would erode all volatiles in 23 Ga. Nevertheless, the simulations do not contradict the observations of Ehrenreich et al. (2012) as we show in Sect. 5.1.2. Furthermore, the planet may have lost 5% of its initial mass at young ages. We caution that these results are based on hydrogen-dominated atmospheres. The planetary composition has most likely been altered over the lifetime of 55 Cnc e and although the planet can maintain a hydrogen-dominated atmosphere on long time scales, large amounts of volatiles are lost in the planetary wind.

4.5. The extreme case of WASP-12 b

WASP-12 b is an exception among all simulated planets. Its extreme proximity of only 2 stellar radii above the stellar photosphere reduces the size of the Roche lobe to only 1.8R_pl above the substellar point. According to the derivation of Erkaev et al. (2007), the gravitational potential of the planet is reduced by a factor of K = 0.28. This causes the strongest mass-loss rate among the planets in our sample, which is, however, 10 000 times weaker than the proposed mass loss of Li et al. (2010) and still 400 times weaker than value of Lai et al. (2010). In this case a multidimensional simulation can result in a higher mass-loss rate, because of the elongated shape of the planetary Roche lobe. The Roche lobe height perpendicular to the star-planet axis is 2/3 × 1.8R_pl = 1.2R_pl (Li et al. 2010), which results in a gravitational reduction factor of K = 0.04. However, this is not sufficient to explain the differences of our result and the published values.

For their mass loss estimates, both Li et al. and Lai et al. derived a density of 3 × 10^-12 g cm^-3 at 1.2R_pl and assume a bulk flow with the speed of sound. The density in our simulation decreases more strongly with height, due to a temperature minimum of 1000 K at 1.014R_pl, resulting in a density that is a factor of 600 lower. Thus, the authors assumption of an isothermal atmosphere with a temperature of 3000 K is incorrect. Furthermore the sound speed is reached at the L1 point and at 1.2R_pl it is a factor 10 lower than the sound speed. These two effects cause the large differences between our simulated mass-loss rate and the estimates of Li et al. (2010) and Lai et al. (2010). While according to Li et al. the planet would have a lifetime of only 10 Ma, our simulation results in about 0.28/0.04 × 0.48% = 3.4% fractional mass loss per Ga. In respect to the planetary mass loss, this resolves the problem of the small chance to discover the planet in a brief period of its lifetime.

4.6. Exobase and stable thermospheres

Hydrodynamic simulations are only valid as long as collisions are frequent enough and a Maxwellian velocity distribution is maintained. This is not always true in planetary atmospheres because of the exponential density decrease with height. At a certain height collisions become too infrequent and atmospheric particles move without collisions on parabolic or hyperbolic paths. This height is called the exobase. Energetic atoms from the tail of the Maxwellian velocity distribution can escape from the exobase, a process called Jeans escape. Jeans escape appears on all planetary bodies, but the escape rates are quite small. For expanded thermospheres of hot gas planets Jeans escape is typically one order of magnitude smaller than hydrodynamic escape (Tian et al. 2005; García Muñoz 2007). Thus, hydrodynamic escape is the dominant escape process, if the hydrodynamic approximation is valid. We verify the validity of the hydrodynamic simulations in retrospect by determining the height of the exobase, which is defined as the height above which the collision probability for atmospheric constituents is unity in the outward direction (e.g., Tian et al. 2005).

The exobase height crucially depends on the collision cross-section. For planetary atmospheres a typical collision cross-section often used is σ = 3.3 × 10^-15 cm² (Chamberlain 1963; Hunten 1973; Tian et al. 2005; García Muñoz 2007). The calculations of Dalgarno (1960) confirm that this value also holds for collisions of neutral hydrogen, σ = 4.3 × 10^-15 cm². However, some authors argue for a collisionless regime closer to the planetary photosphere. They use an H-H collision cross-section of σ = 3 × 10^-17 cm² (Ekenbäck et al. 2010; Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014), which places the exobase deeper in the atmosphere. However, the atmospheres of these close-in planets are always highly ionized at the exobase level. Thus, charge exchange between H-H⁺ with a cross-section of σ ~ 2 × 10^-15 cm² (Ekenbäck et al. 2010; Kislyakova et al. 2014) becomes the dominant process for determining the height at which neutral hydrogen enters a collisionless regime (Guo 2011). Note that the cross-section for proton-proton scattering is yet a factor of 100 larger (e.g. Spitzer 1978; Murray-Clay et al. 2009), so that the ionized wind stays collisional to even larger heights above the planetary photosphere. This result is in agreement with considerations based on collisional relaxation times of protons (using the equations in Goedbloed & Poedts 2004).

The exobase level as defined by Tian et al. (2005) is marked in all plots of the simulated atmospheres for neutral hydrogen and protons, if it is inside the simulated region. For most of our planets the hydrodynamic simulations are valid within the simulation box. In the atmospheres of CoRoT-2 b, WASP-43 b, and WASP-77 b neutral hydrogen decouples from the ionized species between 2.6 and 5.5R_pl (see Fig. 8). The atmosphere of HD 189733 b becomes collisionless for neutral hydrogen at a height of 12R_pl. Neutral hydrogen above this height can be accelerated out of the ionized planetary wind by radiation pressure in the Lyα line (see Sect. 4.8). Nevertheless, the thermospheres are collisional for ionized particles. Above 2R_pl ionization of He⁺ is the main heat source for the wind of CoRoT-2 b and this ion is coupled to H⁺ and electrons, thus, an ionized planetary wind can form.

More importantly, we find that the atmospheres of the massive and compact planets HAT-P-2 b, HAT-P-20 b, WASP-10 b, and WASP-8 b become collisionless before the atmospheric gas is significantly accelerated. Figure 9 shows this for WASP-8 b. The EUV radiation is absorbed deep in the thermosphere, but does not lead to strong acceleration because of efficient radiative cooling. The density gradient is steep and the atmosphere becomes collisionless for neutral hydrogen at 1.7R_pl and for protons at 3.0R_pl. The strongest acceleration takes place above these layers, where the hydrodynamic simulation is invalid and internal energy cannot be converted into a bulk flow. The simulations of these four planets result in very low mass-loss rates (< 3 × 10⁵ g s^-1), but these are unphysical and cannot be interpreted as valid predictions. Even if we consider the uncertainties in our simulations, we conclude that these planets do not support a hydrodynamic wind, but rather compact and stable thermospheres. WASP-38 b and WASP-18 b are similarly compact and massive as WASP-10 b or WASP-8 b, thus, these two hot Jupiters also host stable thermospheres. Naturally, these thermospheres are again subject to Jeans escape and non-thermal escape processes like ion pick-up.

4.7. Ram pressure of the winds

In general, the atmospheres of the compact planets WASP-77 b, WASP-43 b, and CoRoT-2 b are collisional so that planetary winds can form. However, the ram pressures of the planetary winds are only between 0.2 and 0.01 pbar at the sonic point. If we compare these values to a stellar wind pressure of about 10 pbar for a hot Jupiter at 0.05 au (Murray-Clay et al. 2009), it is most likely that these planetary thermospheres are confined over a large portion of the planetary surface. Murray-Clay et al. (2009) argued that a suppressed dayside wind could be compensated for by a stronger wind from the planetary nightside, because the radiative energy input must be spent in some way. Our simulations allow a different scenario. Confined planetary atmospheres exhibit increased thermospheric temperatures, which can lead at least to a partial re-emission of the radiative energy input via Lyα and free-free emission (see Sect. 4.3.2). This would reduce the overall planetary mass-loss rates and not only the dayside wind.

Because of stellar wind confinement, the mass-loss rates of WASP-77 b, WASP-43 b, and CoRoT-2 b are likely even smaller than the results from our simulations. The same argument holds for HD 189733 b with a ram pressure of only 1.0 pbar at the sonic point. In contrast, the wind of HD 209458 b has a ram pressure of 4.3 pbar and that of WASP-80 b reaches 16 pbar. The wind of WASP-80 b is probably strong enough to blow against the stellar wind pressure over the full dayside of the planet. The ram pressures of the other planets range from 6.5 pbar (GJ 436 b) to 480 pbar (WASP-12 b) at the sonic point.

4.8. Radiation pressure

The absorption of stellar Lyα emission by neutral hydrogen in a planetary atmosphere leads to a net acceleration that can exceed the stellar gravitational attraction. If neutral hydrogen is transported beyond the planetary Roche lobe and if it is in a collisionless regime, it can be accelerated to velocities of about 100 km s^-1 (Bourrier & Lecavelier des Etangs 2013). However, our hydrodynamic simulations show that escaping atmospheres stay collisional to large distances above the planetary photosphere (see Fig. 8). Hence, when considering radiation pressure, we must take into account that neutral hydrogen is coupled to the ionized planetary wind through charge exchange (Guo 2011). This reduces the radiation pressure according to the ionization fraction of the atmosphere. The coupling of the ionized and neutral planetary wind constituents should be accounted for in simulation that model the escape with direct simulation Monte Carlo codes (e.g., Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014).

Radiation pressure is included in our simulations as described in Salz et al. (2015a), but in the standard simulations the resolution of the SED does not resolve the stellar Lyα line. While the line flux is correct, the peak is about a factor 10 too small, which also decreases the maximal acceleration by a factor of 10. To evaluate the relevance of radiation pressure, we performed a simplified CLOUDY simulation of a hydrogen cloud with similar conditions to those at the top of the atmosphere of HD 209458 b. We ran the test with a fixed ionization fraction of zero and 90%, using the same SED as in the standard simulation, but with an increased resolution at the Lyα line to fully resolve the central emission line strength.

Figure 10 shows the resulting radiative acceleration. The neutral atmosphere experiences 2.5 times the gravitational acceleration of the host star (a_st = GM_st/a, a: semimajor axis) at the upper most layers of the atmosphere. Shadowing effects lead to a rapid decrease of the acceleration below the irradiated surface. In principle neutral hydrogen could be accelerated away from the planetary thermosphere, which is consistent with the results from Bourrier & Lecavelier des Etangs (2013). However, the upper atmosphere of HD 209458 b is about 90% ionized and collisional. This reduces the radiative acceleration by a factor of ten as depicted in Fig. 10. In the escaping atmosphere of HD 209458 b the radiative acceleration reduces the gravitational attraction of the host star only slightly. The radiation pressure in the Lyα line is more than a factor of 30 smaller than the ram pressure of the planetary wind at the Roche lobe of HD 209458 b. The wind from the dayside of the planet cannot be stopped by radiation pressure, it will rather collide with the stellar wind as described in the previous section.

Among all simulated planets, HAT-P-11 b and HD 97658 b experience the strongest radiative acceleration due to a favorable combination of the distance from the host star, the host star’s Lyα luminosity, and the neutral hydrogen fraction in the upper thermosphere. In a thin layer (width < 0.5R_pl) at the top of the simulated atmosphere these two planets experience a radiative acceleration of half of the gravitational acceleration of the host star. The reduced stellar attraction close to our upper boundary does not affect the structure of the thermospheres below 14R_pl. In all other systems the radiative acceleration is less than a tenth of the gravitational acceleration of the host star.

In summary, acceleration due to radiation pressure in the Lyα line can be neglected in 1D simulations of the planetary wind formation. It only reduces the effective stellar gravitational attraction in the uppermost thermospheric layers. 3D Models simulating the interaction of the evaporating atmosphere with the stellar wind and the stellar radiation pressure should consider the coupling of neutral hydrogen to the ionized planetary wind instead of removing ionized planetary wind particles from their simulations (e.g., Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014).

5. Planetary Lyα absorption and emission signals

The systems simulated in this paper were selected based on a prediction of the host stars’ Lyα fluxes at Earth, so that all host stars are potentially bright enough to detect the escaping planetary atmospheres in absorption. We first calculate Lyα absorption signals from the simulated atmospheres, which are used to identify the most promising targets for future transit spectroscopy. For systems with observed signals, we compare our results with the measurements. Second, we estimate the planetary Lyα emission line strength and consider the detectability of these emission signals.

5.1. Lyα absorption signals

Our simulations provide information on the possible strength of Lyα absorption during transit observations, because the more neutral hydrogen is transported into the planetary thermosphere and beyond, the more absorption occurs. However, the observed planetary Lyα absorption signals show relatively large radial velocity offsets and an asymmetric blue shift (e.g., Lecavelier des Etangs et al. 2012; Ehrenreich et al. 2015). Among others, direct simulation Monte Carlo codes have been successfully used to model the observed signals by including interactions with the stellar wind and the Lyα radiation pressure (Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014). It is clear that the complex interaction region between the stellar and planetary winds cannot be approximated with our simplified 1D models. However, also in more complex models the neutral hydrogen mass-loss rate is adjusted to fit the measured planetary absorption signal and larger mass-loss rates result in stronger Lyα absorption (Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014).

In addition to the planetary neutral hydrogen mass-loss rate, the measured absorption signals are affected by several other system parameters. The stellar high energy irradiation level affects the lifetime of neutral hydrogen and the planet-to-star size ratio influences the contrast of the signals. Stellar wind parameters and Lyα radiation pressure affect the production of energetic neutrals and the shape of the escaping atmosphere. Both processes cause enhanced absorption in the Lyα line wings, where the stellar flux is less affected by interstellar absorption. Furthermore, strong planetary magnetic fields can shield the planetary atmosphere and affect the Lyα absorption signals (Trammell et al. 2014; Kislyakova et al. 2014).

To estimate the strength of the planetary Lyα absorption, we compute the absorption signals from our spherically symmetric planetary atmospheres (see Appendix A). The resulting values reflect the combined impact of the strength of the planetary mass-loss rate, the lifetime of neutral hydrogen in the upper atmosphere, and the planet-to-star size ratio. The absorption is first calculated below the Roche lobe height and second for the complete simulated atmospheres (up to a height of 12/15R_pl). This twofold approach was chosen, because the signal from below the Roche lobe is caused by the rather unperturbed planetary wind, which can be approximated by our simulations. The signal of the full atmospheres is a measure of the neutral hydrogen transported beyond the Roche lobe. In Table 3 the equivalent width of both signals is given; the values in brackets are the additional equivalent width of unbound hydrogen. In Appendix A we further show absorption profiles and simulated transmission images. In the following, we do not use the term exosphere, because according to our estimates the escaping atmospheres are collisional up to large distances above the planetary limb (see Sect. 4.6), rather we distinguish between bound and unbound hydrogen.

5.1.1. Absorption signals from nearby systems

Our results show that in general low-potential planets produce deeper absorption signals, because they host strong and neutral winds (see Figs. A.2 and A.3). All of the low-potential planets qualify for Lyα transit measurements in terms of the expected absorption signal. WASP-80 b is one of the most promising systems in the search for yet undetected expanded atmospheres. The planet’s size is large compared to its host star and the Roche lobe has an extent of 0.85 R_st. Thus, the planetary atmosphere can cover almost the complete host star. This results in the strongest predicted absorption signal with an equivalent width of 320 mÅ not including excess absorption of unbound hydrogen.

The simulated absorption signal of GJ 3470 b is even stronger than that of WASP-80 b if unbound hydrogen is included, thus, strong absorption in the post-transit phase could be present. Although GJ 1214 b also produces a large absorption signal, the host star was not detected in a Lyα observation with the Hubble Space Telescope (France et al. 2013), so that the absorption probably cannot be measured today.

HD 97658 b produces the smallest total absorption signal of the low-potential planets due to the unfavorable planet to star size ratio. We do not predict absorption of unbound hydrogen for this planet, because the Roche lobe is larger than our simulation box. Our computed signal is slightly stronger than that of 55 Cnc e, where evaporation was not detected during a dedicated observation campaign (see below). HD 97658 b is highly interesting, because among those systems with detectable evaporation, it is probably the closest proxy to the young Venus, which is believed to have been enshrouded in a dense protoatmospheric hydrogen envelope (Lammer et al. 2008). The mass of this planet is only 8 Earth-masses and the irradiation level is similar to that of Venus at a solar age of 100 Ma (~1000 erg cm^-2 s^-1, Ribas et al. 2005).

The large mass-loss rate of WASP-12 b results in strong absorption of unbound hydrogen. In principle, the strong absorption agrees with the detection of excess absorption in metal lines (Fossati et al. 2010), but Lyα emission cannot be measured in this system due to its large distance. In the atmospheres of WASP-43 b, WASP-77 b, or CoRoT-2 b neutral hydrogen remains confined within the Roche lobe. The equivalent width of the absorption signals for these three planets remains lower than that of 55 Cnc e. The hot and stable thermospheres of HAT-P-2 b, HAT-P-20 b, WASP-10 b, and WASP-8 b produces about 1% of excess absorption in the Lyα line, which is challenging for a detection. The same applies to the planets WASP-18 b and WASP-38 b, which were not simulated.

5.1.2. Comparison with measured absorption signals

In HD 209458 b, HD 189733 b, GJ 436 b, and for the non-detection in 55 Cnc e, we can compare the computed absorption signals with observational results. The equivalent widths of the measured and simulated absorption signals of the four systems are provided in Table 4. Planets with strong neutral winds produce large absorption signals in our simulations, which is qualitatively consistent with the observations. Figure 11 shows the spectral shape of the simulated absorption signals in comparison with the measurements. Note that the strong absorption in the line center cannot be detected because it would be extinguished by interstellar absorption. The figure also contains calculated transmission images, which correlate with the absorption signals as follows: the strong absorption signal in the line center is produced by the thin upper atmospheric layers that cover a large percentage of the stellar disk, and the absorption in the line wings is caused by the dense, lower thermosphere. The results for 55 Cnc e are given in Fig. 12.

The top panel in the left column of Fig. 11 shows the absorption signal of HD 209458 b. The system exhibits the most symmetric absorption signal measured so far. The absorption was discussed in detail by Murray-Clay et al. (2009) and by Koskinen et al. (2013a). HD 209458 b is the only system, where the absorption signal can be solely explained by a spherical atmosphere, depending on the height of the H/H⁺ transition layer. In the model of Koskinen et al. the transition occurs at a height of 3.4R_pl, so that their predicted absorption fits the transit observations. Our simulation results in an H/H⁺ transition layer at a height of 1.6R_pl, which does not produce sufficient absorption in the line wings. A twice as strong planetary wind already shifts the ionization front to 2.3R_pl as seen in our H-only simulation (see Sect. 4.2), but the absorption remains insufficient (see Fig. 11). The fact that changes within the uncertainty of the current simulations decisively affect the absorption signal renders it difficult to determine its origin. For example, Trammell et al. (2014) reproduce the observed signal with an atmosphere that is largely stabilized through a strong planetary magnetic field; in contrast, Kislyakova et al. (2014) reproduce the signal with an exospheric model including interactions with the stellar wind and radiation pressure. The current data is likely insufficient to distinguish between the different models.

In the case of HD 189733 b it is impossible to explain the spectral shape of the measured absorption signal with a spherically expanding atmosphere. The absorption of HD 189733 b is found at radial velocity shifts of up to −250 km s^-1 and around + 100 km s^-1 (see Fig. 11). Even considering all uncertainties, a spherical wind will not produce such a signal.

GJ 436 b shows similar absorption signals during and after the planetary transit, which indicates a cometary tail like structure of the escaping atmosphere (Kulow et al. 2014; Ehrenreich et al. 2015). Again, the signal cannot be explained by a spherical atmosphere, but interaction with the stellar wind or radiation pressure is required to reproduce the detected radial velocity shift and the temporal evolution of the absorption (Ehrenreich et al. 2015). In our simulations, GJ 436 b has one of the largest amounts of unbound neutral hydrogen, which agrees with the detection of a cometary tail. Although both HD 209458 b and GJ 436 b produce strong neutral winds, the detected Lyα absorption of HD 209458 b is much weaker and this planet shows no post-transit absorption. There are several possibilities to explain these differences. If our estimate of the XUV irradiation level is too low, HD 209458 b would host a stronger but also more ionized planetary wind, which reduces the amount of unbound neutral hydrogen. Also a partially confined atmosphere or the interaction with the stellar wind could reduce the amount of unbound neutral hydrogen in HD 209458 b. Finally, Bourrier et al. (2015) explain such differences based on their radiation pressure driven exosphere model.

Excess absorption was not detected during the transit of 55 Cnc e. There is a weak signal in the blue line wing, but the false alarm probability is 89% (Ehrenreich et al. 2012). The absorption produced by the simulated escaping atmosphere including unbound hydrogen is no larger than 8% in the line center (see Fig. 12). Although the planet can host a strong hydrogen-dominated wind, the simulated absorption signal remains small, mainly because the planet is small compared to its host star, but also because of the short lifetime of neutral hydrogen in the highly irradiated atmosphere. The simulated signal is close to the detection limit and, considering additional interstellar absorption, the simulation results agree with the non-detection. Therefore, the observations by Ehrenreich et al. (2012) do not exclude the existence of a hydrogen-dominated thermosphere that currently undergoes strong evaporation.

5.2. Lyα emission signals

Our simulations show that planetary atmospheres are strongly cooled by Lyα emission. Especially high-potential planets re-emit large amounts of the XUV energy input through Lyα emission. This poses the question whether the planetary Lyα emission is detectable from Earth. Based on the simulations we can predict the planetary Lyα emission line strength, which can be compared with the stellar line fluxes (see Sect. 2.2).

We compute the ratio between the planetary and stellar Lyα line fluxes φ_pl/φ_st, similar to the approach of Menager et al. (2013). Our photoionization solver does not spectrally resolve the stellar and planetary Lyα lines, hence, in a first step we scale our values to derive the flux in the line center. For the stellar emission lines we used a ten times broader line width than usually measured for stellar line profiles (see Sect. 2.2), thus, the central stellar Lyα line flux will be a factor 10 higher than in our reconstructions. The planetary emission is given by the photoionization solver with a bin width of Δλ = 6.1 Å. We derive the central emission line strength by assuming a Gaussian shape. In the models of Menager et al. (2013) the planetary emission has a FWHM of about 0.13 Å. Assuming an atmospheric temperature of 10 000 K we derive a similar value of (5)Here, b is the Doppler parameter, v_th the thermal Doppler velocity, and for the microturbulence v_micro we use the sound speed. Comparing the two results, we regard a FWHM of 0.15 Å a reasonable value for the planetary line profile.

The central emission line flux from the planet is given by conserving the integrated line flux: (6)Here φ_pl,CLOUDY Δλ is the integrated planetary flux given by the photoionization solver and the righthand side of the equation is the integrated flux from the assumed Gaussian line profile. Finally, we scale the flux ratios φ_pl/φ_st with the relative sizes of planetary and stellar disks, where we assume that only the hot dayside of the planet emits Lyα radiation. We adopt a height of 1.2R_pl for the Lyα emission layer in the planetary atmospheres, which is indicated by the dip in the heating fractions (see panels f in Fig. 8).

The derived values for HD 209458 b and HD 189733 b can be compared to those of Menager et al. (2013), who find flux ratios of 2.1 × 10^-3 and 1.5 × 10^-3 respectively. Our value of 1.8 × 10^-4 for HD 209458 b is one order of magnitude smaller. However, the two results are based on different atmospheric models. Menager et al. (2013) used the atmospheric models from Koskinen et al. (2013a) with a 1200 K hotter thermosphere and a larger H/H⁺ transition layer height, which is a possible reason for the differences. Our flux ratio of 1.3 × 10^-3 for HD 189733 b agrees with the results from Menager et al. (2013), however, these authors used a solar irradiation level, which is about 20 times weaker than our value. Since a higher XUV irradiation increases the flux ratios, it seems that our simulations generally produce lower planetary Lyα luminosities than those obtained by Menager et al. (2013). Yet, the reasons for the differences can only be assessed with a detailed analysis based on an identical atmospheric model.

We now focus on the maximum flux ratios (φ_pl/φ_st)_max, which approximately occur when the planetary emission experiences the maximal radial velocity shift through the orbital motion of the planet at phases 0.25 and 0.75. This causes radial velocity shifts between 0.3 and 0.9 Å in our sample of planets. The reconstructed stellar Lyα line profiles from Wood et al. (2005) show that the line flux drops by a factor of 10 compared to the central emission line strength if the line shift exceeds 0.5 Å. Only for stars with unusually wide line profiles this assumption breaks down. Therefore, we apply a factor of 10 to the flux ratios, assuming that the observations are placed at the corresponding phases.

The calculated flux ratios depend on several parameters. Naturally, a larger size of the planet increases the ratio and a smaller semimajor axis also increases the ratio through higher irradiation levels. Figure 13 shows the flux ratios for the simulated planets versus the planetary gravitational potential. Massive and compact planets show larger flux ratios, because their thermospheres are close to radiative equilibrium and convert the XUV energy input efficiently into Lyα radiation. The total planetary Lyα line fluxes at the top of the planetary atmospheres reach maximally 30 to 40% of the combined irradiating stellar Lyα and XUV flux. The maximal flux ratios at orbital phases 0.25 and 0.75 are smaller than 10%. The strong emission of WASP-12 b is caused by the extreme proximity of this planet to its host star. GJ 1214 b shows a relatively large flux ratio, because the host star emits little Lyα radiation, but the planet converts part of the XUV energy input into Lyα emission. The Lyα emission of WASP-43 b and WASP-77 b reach 3.5% and 2.6% of the stellar flux in the line wing. These signals may be detectable, depending on the stellar line flux at Earth and the specific stellar and planetary line profiles.

The presented results are somewhat rough estimates because of the described limitations. They remain valid if metals are included according to our test simulation of HD 209458 b with solar metal abundances (see Sect. 3.7). In principle, one could assume that the Lyα emission of planetary thermospheres containing metals is significantly reduced, because metals are efficient cooling agents, but in HD 209458 b the emission was only reduced by a factor of two. This general factor should not impact the comparison between the individual systems.

A detection of such Lyα emission from planets would substantiate our finding that massive and compact hot gas planets host stable thermospheres. Only the conversion of stellar XUV flux into Lyα emission in the hot and stable thermospheres boosts the flux ratios to detectable values. Furthermore, the mechanism allows a distinction between two classes of thermospheres in hot gas planets. Figure 13 shows that small planets with strong and cool winds create large absorption signals and only weak Lyα emission; massive and compact planets with close to stable thermosphere create little absorption but re-emit the XUV energy input partially via Lyα emission. Therefore, planetary Lyα absorption and emission are to some extent mutually exclusive and the onset of strong planetary winds around log − Φ_G = 13.2 erg g^-1 provides an observational test of the hydrodynamic escape model. According to our simulations, WASP-43 and HD 189733 are equally suited for an observational campaign aimed at the detection of planetary Lyα emission. While the expected signal is a factor 2.7 larger in the system WASP-43, the predicted Lyα flux is a factor 10 smaller than that in HD 189733, hence, the signal-to-noise ratio should be similar.

6. Summary and conclusion

We present coupled photoionization-hydrodynamics simulations of the planetary winds from hot gas planets in the solar neighborhood. The most promising candidates in terms of Lyα transit spectroscopy were selected for the simulations, which were performed with our TPCI code. We have described the behavior of the escaping atmospheres in detail with a focus on the differences between winds from small planets and massive, compact planets. The use of a photoionization and plasma simulation code, which gives a precise solution for the radiative cooling, has proven to be crucial for the results presented here. In particular, we find that compact and massive hot Jupiters host stable thermospheres, because of strong radiative cooling in their hot thermospheres. Radiative cooling is weak in smaller planets, because adiabatic cooling in the rapidly expanding atmospheres efficiently reduces the thermospheric temperatures. In our simulations of hydrogen and helium atmospheres, Lyα and free-free emission are the most important cooling agents. Specifically, we found stable thermospheres in the planets HAT-P-2 b, HAT-P-20 b, WASP-8 b, WASP-10 b, WASP-38 b, and WASP-18 b.

The simulated thermosphere of HD 189733 b reproduces the temperature-pressure profile reconstructed from measured sodium absorption signals, but in HD 209458 b our simulations and the observations show significant differences that currently remain unexplained.

Lyα absorption signals of the spherically symmetric atmospheres are presented to assess the systems’ potential for transit spectroscopy. The calculations reproduce large absorption signals, where strong Lyα absorption has been detected. Our simulations show that smaller planets tend to cause larger Lyα absorption signals, because they produce strong and cool winds that transport large amounts of neutral hydrogen into the upper thermospheres. In contrast, the stable thermospheres of massive and compact planets are hot and highly ionized and, thus, only create small Lyα absorption signals. Instead, such hot thermospheres are strong Lyα emitters, and the planetary emission could be detectable with current instrumentation. The distinction between planets with strong Lyα absorption but no emission and planets with little Lyα absorption but strong emission can be used to distinguish between these two classes of thermospheres. Based on the predicted strong absorption signals, WASP-80 b and GJ 3470 b are the most promising targets among our sample in the search for yet undetected escaping atmospheres.

---

Acknowledgments

M.S. acknowledges support through Verbundforschung (50OR 1105) and the German National Science Foundation (DFG) within the Research Training College 1351. P.C.S. acknowledges support by the DLR under 50 OR 1307 and by an ESA Research Fellowship. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.

References

Appendix A: Calculation of the Lyα absorption signals

Here, we describe the method used to calculate the Lyα absorption in the simulated atmospheres. Although we use the term Lyα absorption, we note that the actual physical process is scattering of photons out of our line of sight. Absorption of Lyα with collisional de-excitation only occurs deep in the simulated atmospheres (see Fig. 6). We do not account for any emission or line transfer in the computation of the absorption signals, or for scattering of photons into our line of sight. We further assume a homogeneous brightness of the stellar disk in Lyα emission.

To compute the absorption signals we shoot rays through the spherical atmospheres (see Fig. A.1). The wavelength dependent transmission T^λ in a single ray with index j is given by the product of the absorption in each cell: (A.1)Here, n_i is the neutral hydrogen density in cell i with the cell width Δs_i. is the standard Lyα absorption cross-section with a Voigt-profile, shifted by the local tangential velocity of the gas. The Doppler width is given by the sum of the local thermal and turbulent widths (see Eq. (5)); the latter is given by the speed of sound. Deuterium absorption is included with a deuterium fraction of (D/H) = 2.2 × 10^-5 according to measurements from Jupiter (Lellouch et al. 2001).

The total transmission is given by: (A.2)Here, n ≈ 100 is the number of rays and Δh is half of the vertical resolution. In this equation, the first term accounts for the fraction of radiation that bypasses the simulated atmospheres (transmission equals one). In most of the systems this term is zero, because the simulated atmospheres cover the full stellar disk. The sum adds the transmission of each ray weighted with the fractional area that a ring with the given height R_j above the planetary photosphere covers of the stellar disk. The transmission of the planetary disk is zero and does not add to the total transmission (R₁ = R_pl). For each planet we compute the equivalent width of the absorption signal over the full line (± 500 km s^-1) (A.3)These values are given in Table 3 for the atmosphere below the Roche lobe and the values for unbound atmospheric material above the Roche lobe is given by the value in brackets.

We further fold transmission signal with the line spread function of the 52′′ × 0.1′′ aperture of the Space Telescope Imaging Spectrograph on board the Hubble Space Telescope. This reduces the sharpness of the predicted absorption signals and provides results comparable with measurements. The spectral shapes are shown in Figs. A.2 and A.3. Our procedure is similar to the description in Koskinen et al. (2010), but we do not include interstellar absorption in the calcuations, because the main result is the potential Lyα absorption of the atmsosphere. For example, interactions with the stellar wind can cause large radial velocity shifts of the absorbing neutral hydrogen, which will result in absorption in the Lyα line wings rather than in the line center. Such signals are less affected by interstellar absorption. For visualization we also compute transmission images, which show how the planets would be observed during transit locally in the systems (without interstellar absorption). To create these images we apply the wavelength dependent transmission to an average stellar Lyα profile with a double Gaussian shape (Salz et al. 2015b) and then integrate the total absorption of stellar Lyα flux for each ray. Plotting the individual atmospheric rings over the stellar disk creates the transit images.

All Tables

All Figures