© ESO, 2016

1. Introduction

A substantial fraction of known exoplanets orbit extremely close to their stars, within a tenth of an astronomical unit. This population displays a wide diversity in nature, from gaseous giants to super-Earths and even disintegrating rocky cores, raising many questions about the formation and evolution of such objects and their possible relations. The intense X-ray and extreme ultraviolet irradiation from the parent star can lead a hydrogen-rich thermosphere to lose its stability, expand hydrodynamically, and evaporate (Vidal-Madjar et al. 2003).

An extended atmosphere produces a much deeper absorption than the lower atmospheric layers when observed in the spectral lines of elements that are abundant at high altitudes and/or associated with strong electronic transitions. Transit observations in the Lyman-α line of their host star first revealed the existence of extended atmospheres of neutral hydrogen around the hot Jupiters HD 209458 b (Vidal-Madjar et al. 2003, 2004; Ben-Jaffel 2007, 2008; Vidal-Madjar et al. 2008; Ehrenreich et al. 2008; Ben-Jaffel & Sona Hosseini 2010) and HD 189733 b (Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013). Observations of other lines in the far UV allowed the detection of heavy metals and ions carried to high altitudes by collisions with the hydrogen outflow (e.g., oxygen, carbon, and magnesium; Vidal-Madjar et al. 2004; Linsky et al. 2010; Vidal-Madjar et al. 2013; Ben-Jaffel & Ballester 2013; Fossati et al. 2010; Haswell et al. 2012), confirming that the upper atmosphere of these giant planets is in a state of blow-off.

Many questions remain to be answered about the effect of irradiation on the formation and the structure of extended atmospheres. The strong energy input into hot Jupiters is the source for their evaporation (e.g., Lammer et al. 2003; Lecavelier des Etangs et al. 2004; Yelle 2004, 2006; García Muñoz 2007; Owen & Jackson 2012; Koskinen et al. 2013a,b; Johnstone et al. 2015), but the detection of a Lyman-α absorption from the partially transiting upper atmosphere of the warm Jupiter 55 Cnc b first hinted that milder conditions of irradiation can lead to atmospheric expansion (Ehrenreich et al. 2012). Interestingly, 55 Cnc b is close to the orbital distance limit where the atmosphere of a giant planet is expected to lose its stability (Koskinen et al. 2007). The moderate depth of 55 Cnc b absorption signature (Ehrenreich et al. 2012), in comparison to HD 209458 b and HD 189733 b, can be explained if its exosphere fills about a third of the Roche lobe, consistent with the lower irradiation of the planet.

Recently, transit observations of the warm Neptune GJ 436 b in the Lyman-α line at three different epochs revealed deep, repeatable exospheric transits reaching up to 60%, and lasting for much longer than the optical transit (Kulow et al. 2014; Ehrenreich et al. 2015). “Radiative braking” (Bourrier et al. 2015b) and interactions with the wind of its M dwarf host star (Bourrier et al. 2016) were shown to shape the neutral hydrogen exosphere of GJ 436b into a giant coma that surrounds the planet and trails for millions of kilometers behind it. Counterintuitively, GJ 436 b observations revealed that much larger atmospheric signals could be retrieved from the upper atmosphere of moderately irradiated, low-mass planets than from more heavily irradiated massive planets. This opens thrilling perspectives for the characterization of the many small planets in the sub-Neptune and super-Earth mass regimes, which have no equivalent in the solar system and whose nature and evolution remain mysterious. Probing the upper extended atmosphere of such planets will be particularly important because their lower atmosphere, less extended and with a heavier composition, will be more difficult to characterize. The non-detection of a neutral hydrogen exosphere around the super-Earth 55 Cnc e (Ehrenreich et al. 2012) hinted at the presence of a high-weight atmosphere – or the absence of an atmosphere – recently supported by the study of its brightness map in the IR (Demory et al. 2016). Studying the atmospheric escape of small planets will not only bring insights into their nature, but will also help understand their formation and evolution.

The exoplanet HD 97658 b (M_p = 7.9 ± 0.7 M_⊕; R_p = 2.39 R_⊕) was discovered by a radial-velocity survey (Howard et al. 2011), and was later detected in transit across its 9.7 ± 2.8 Gy old parent star (Dragomir et al. 2013; Bonfanti et al. 2016). Measurement of the planetary radius at 4.5 μm with Spitzer have allowed the planet bulk density to be refined to  g cm^-3 (Van Grootel et al. 2014). While its exact nature remains to be unveiled, Van Grootel et al. (2014) have predicted that HD 97658 should have a large rocky core accounting for up to 60% of the planet mass, an envelope of water and ices accounting for up to 40% of the planet mass, and less than 2% of hydrogen and helium. If its composition is indeed water-rich, HD 97658 orbits close enough to its K1-dwarf star (P = 9.5 days, T_eq ~ 725 K; Knutson et al. 2014) that its temperature should allow for the formation of a thick envelope of steam. The brightness of HD 97658 b host star (V = 7.7, d = 21.1 pc) makes it one of the handful of small planets currently amenable to atmospheric characterization. However, given the shallowness of the planetary transit seen at optical or near-infrared (NIR) wavelengths (~0.09%) and the small expected scale height of a steam envelope, its detection would be extremely challenging through infrared absorption. As it is, transmission spectroscopy of HD 97658 b performed in the near IR with WFC3 (Knutson et al. 2014) is consistent with flat transmission spectrum models, indicating either a cloudy atmosphere or a cloudless, water-rich atmosphere. Consequently, HD 97658 b is a very good target for transit observations of its upper atmosphere. Photodissociation of a steam envelope could produce OH and H at high altitudes (Wu & Chen 1993; Jura 2004). Here, we present a search for neutral hydrogen escape through dedicated Lyman-α line observations. These observations are complemented by X-ray measurements to estimate the high-energy stellar irradiation. This study is comparable to the one performed for the hot Jupiter HD 189733 b, which transits a star with the same type as HD 97658 and was shown to be evaporating (Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013). Differences in the age and coronal activity of the host stars, in addition to the different orbital and bulk properties of their planets, may induce very different evaporation regimes.

Lyman-α observations of HD 97658 at different epochs are analyzed and interpreted in Sect. 2. They are used in Sect. 3 to estimate the intrinsic Lyman-α line of HD 97658 and the properties of the interstellar medium (ISM) in its direction. We use these measurements, along with that of the stellar X-ray emission, to study the XEUV irradiation of HD 97658 b in Sect. 4. In Sect. 5 we compare the Lyman-α observations with numerical simulations of the planet transit with the EVaporating Exoplanet (EVE) code to constrain the presence and properties of a putative extended thermosphere and exosphere of neutral hydrogen. We discuss our results and their importance for the nature of HD 97658 b in Sect. 6. The properties of the system relevant to our study are given in Table 1.

2. Analysis of the resolved Lyman-α line

2.1. Transit observations

We observed HD 97658 in the H i Lyman-α line (1215.6702 Å) with the Space Telescope Imaging Spectrograph (STIS) instrument on board the Hubble Space Telescope (HST). Three visits of five HST orbits each were obtained in the frame of the GO Program 13820 (PI: Ehrenreich), scheduled so that Visit 1 (April 2015) covered one transit of HD 97658 b, Visit 2 (December 2015) covered the transit and post-transit phase, and Visit 3 (March 2016) covered the pre-transit and transit phase (Table 2). The time-tagged data obtained with the G140M grating were reduced with the CALSTIS pipeline. Each HST orbit was divided into ~300 s exposures, as a compromise between signal-to-noise ratio (S/N) and temporal resolution. In all visits, this yielded six subexposures for the first orbit (the target acquisition reduces the duration of its scientific exposition) and seven for the other orbits. In the raw data (Fig. 1), the stellar Lyman-α emission line is superimposed with the geocoronal airglow emission from the upper atmosphere of the Earth (Vidal-Madjar et al. 2003). The airglow contamination was limited by the use of STIS narrow slit of 52′′ × 0.05′′, and can be estimated and removed from the final 1D spectra using CALSTIS. It is nonetheless recommended to treat with caution the regions where the airglow is stronger than the stellar flux. The strength and position of the airglow varies with the epoch of observation, and in a first step we excluded from our analysis the spectral ranges [−27; 12] km s^-1 (Visit 1), [−13; 40] km s^-1 (Visit 2), and [−21; 31] km s^-1 (Visit 3), defined in the stellar rest frame. We note that correcting manually for the airglow, using different areas of the 2D images to build its profile, did not change the final spectra significantly. For all visits, we also limited the analysis to the spectral range where the S/N is high enough, from −305 to 305 km s^-1.

2.2. Systematic variations

We checked for variations on short timescales within a given HST orbit caused by the telescope “breathing”. This effect is known to affect STIS observations in the UV and in the visible because of variations of the telescope throughput caused by the changes in temperature that the HST experiences during each orbit (e.g., Bourrier et al. 2013). The shape and amplitude of the breathing variations can change between visits of the same target, but the orbit-to-orbit variations within a single visit are both stable and highly repeatable (e.g. Brown et al. 2001; Sing et al. 2008; Huitson et al. 2012; Ehrenreich et al. 2015). We thus modeled the breathing effect using a Fourier series decomposition with a period equal to that of the HST around the Earth (P_HST = 96 min), and adjusted the following function to the 300 s exposure spectra in each visit, (1)with φ_HST(t) the HST orbital phase at a given absolute time t and f_breath the breathing model function. Its degree n and coefficients a_i are free parameters of the fit. The function F_nom corresponds to the nominal stellar flux unaffected by the breathing effect. Its value for each orbit, indexed by k_orb, is a free parameter of the model. This accounts for variations in the flux from one orbit to another arising for example from intrinsic stellar variability. We note that we tried defining F_nom as a polynomial function with its degree as a free parameter; there was no significant change in the results of the fit.

The breathing effect is achromatic, but it must be analyzed in wavelength bands where no spectral variations are found as a function of time that could be caused by the planetary atmosphere, which would risk being overcorrected. Typical exospheric signatures are located in the blue wing of the line, and the red wing can usually be used as a stable reference to characterize systematic variations. However, in the present case we found no spectro-temporal signatures that could clearly be attributed to the planet (this was validated a posteriori in Sect. 2.3), and for each visit a thorough search for the breathing effect in all parts of the line yielded similar results for f_breath. We thus fitted the F_corr function to the flux integrated in the entire Lyman-α line, using the Bayesian Information Criterion (BIC) as a merit function (Crossfield et al. 2012; Cowan et al. 2012).

Results are shown in Fig. 2. For Visits 1 and 3, we detected a significant breathing effect with order n = 1, yielding respective BICs of 57 and 63. By comparison, the same fits performed with no breathing component (n = 0) yielded BICs of 110 and 85. We detected no significant breathing variations in Visit 2: n = 0 and 1 yielded similar BICS of 66, and n> 2 yielded BICs > 72. Variations of the stellar flux from one orbit to the other were derived from F_nom in Visits 1 and 2, and are discussed in Sect. 2.3.2. We caution that the scientific exposure in the first HST orbit is shorter and begins at a later HST orbital phase than the other orbits (Fig. 2). Therefore, because the breathing-affected flux typically increases with HST phase (e.g., Bourrier et al. 2013; Ehrenreich et al. 2012; Ehrenreich et al. 2015), its average over the full duration of the first orbit can be overestimated compared to that of the other orbits. This underlines the importance of correcting for the breathing effect, which can only be done using data acquired in time-tag mode. Finally, we used the best fits for F_corr to correct each 300 s exposure spectrum by the value of the breathing function f_breath at the time of mid-exposure (Fig. 3).

2.3. Line profile variations

A preliminary analysis of the breathing-corrected data revealed variations between the observed spectra, but none that could be attributed to the transit of an exosphere. Because of these variations, no exposure could be clearly identified as a reference for the unocculted, average, stellar Lyman-α line. Therefore, this reference was defined as the error-weighted mean of all exposures in a given visit and used to perform a more thorough analysis of the observations.

2.3.1. Long-term evolution

First, we compared the reference spectra of the different visits to assess long-term variability in the shape of the Lyman-α line (Fig. 4). All visits show remarkably similar line profiles at high velocities in the wings of the line, with no variations at more than 2σ beyond about ±187 km s^-1 (see also Fig. 5). The only exception is a localized flux increase of 27 ± 11%, from about 170 to 215 km s^-1, in Visit 2. The reference spectra are also similar in the core of the line (−36 to 36 km s^-1). Although the Visit 2 spectrum displays spurious pixels that are likely biased by its higher airglow emission (Fig. 1), this gave us confidence that the background was properly corrected in the core of the Visits 1 and 3 spectra, and hereafter this range was included in our analysis of these visits.

Except for the core and the red wing signature, the reference spectra of Visits 2 and 3 are nearly identical. In contrast, the Visit 1 spectrum shows a significant flux increase (~11.3 ± 2.6% at 4σ) in two wide bands roughly symmetrical in the stellar rest frame (36 < | v | < 187 km s^-1; Figs. 4, 5). These regions of higher flux at the peaks of the observed line likely trace the parts of the intrinsic line emitted by the lower transition region between the upper chromosphere of the star and its corona, that is at temperatures near 20 000 K in the case of the quiet Sun (Vernazza et al. 1981). The flux increase in Visit 1 may thus indicate a higher level of chromospheric activity at this epoch. The similarity between the spectra of Visits 2 and 3, and their difference with the Visit 1 spectrum, is also consistent with long-term variability of the star since Visit 2 was secured ~8.1 month after Visit 1, while only ~2.8 months separated Visits 2 and 3. We note that ISM absorption is highest in the core of HD 97658 Lyman-α line, and its intrinsic flux is thus lost from about −30 to 50 km s^-1 (Sect. 3; Fig. 6). The flux in the core of the observed line comes from the convolution of the ISM-absorbed spectrum farther from the core by STIS instrumental response, and while it still contains information that can be used to reconstruct the intrinsic line (Sect. 3), we found that it is too blended to reflect directly the line variability between the epochs. Finally, the line wings beyond 187 km s^-1 are little affected by ISM absorption and directly trace the intrinsic line. The stability that they show over the three visits is consistent with their formation in the colder regions of the chromosphere.

2.3.2. Short-term variations

Our goal in this section is to characterize the variability between the exposures of each visit, and to search for transit signatures of an extended atmosphere around HD 97658 b. This was done by calculating the ratio of each exposure spectrum to the reference spectrum of the corresponding visit. The S/N per pixel is too low in a 300 s exposure to study spectral features, and we used the spectra averaged over the full duration of each HST orbit. We identified features characterized by relative flux variations with S/N larger than 2, and extending over more than three pixels (i.e., we excluded intervals narrower than STIS spectral resolution, about 0.1 Å). Using the bootstrap-like method detailed in Bourrier et al. (2013), we also calculated the false-positive probability (p-value) to have features as significant as the ones we detected but caused by noise alone. The results of this analysis are given in Table 3.

Most known evaporating giants are characterized by an extended upper atmosphere of neutral hydrogen partially filling or overflowing the planetary Roche lobe. Gravity and interactions with the host star (in particular radiation pressure and charge-exchange) shape the escaping gas from a dense, low-velocity coma around the planet into a diluted, fast-moving cometary tail trailing behind the planet. Observationally, this atmospheric structure results in a deep absorption in the core of the Lyman-α line caused by the transit of the coma, and the absorption decreases and shifts to negative velocities in the blue wing of the Lyman-α line when the exosphere later transits the star. We found none of these features in our analysis of HD 97658 Lyman-α line transit observations. While deviations with respect to the reference spectrum were detected in each visit, most of them have S/N lower than 3 and p-values up to 70% (Table 3). Even if an exosphere-absorbed spectrum had been included in the weighted mean of the spectra used to build the reference spectrum, residual variations should still occur in the spectral range of this absorption signature. Instead, the measured variations do not show any correlation from one spectrum to another: they arise in all parts of the line, at no specific time relative to the planet transit, and with no consistency between the different epochs.

Interestingly, as can be seen in Table 3, Visit 1 is more variable than the other visits and Visit 2 displays a significant variation (S/N> 3 σ, p-value < 5%) in the red wing. To further investigate this short-term variability we compare in Fig. 5 the light curves obtained by integrating the flux in the five bands identified in Sect. 2.3.1: the blue and red “far” wings, the core of the observed line, and its peaks. Within the statistical variations, the flux is stable in the far wings and the core of the line, not only between the different epochs, but also over a timescale of a few hours. In contrast, the regions of higher flux at the peaks of the observed line display significant systematic variations. Given that these spectral bands are the same that show long-term variability between Visit 1 and the other epochs (Sect. 2.3.1), this short-term variability could also be caused by activity in the hot upper chromosphere, stronger in Visits 1 and 2.

We thus conclude that there are no detectable transit signatures of an extended atmosphere around HD 97658 b. The changes observed between the different spectra in each visit arise mainly from statistical variations around the reference spectrum, possibly from intrinsic variability of the stellar chromosphere, and in any case with no variations consistent with a transit light curve.

3. Intrinsic Lyman-α line

The observed Lyman-α line is affected by absorption from the ISM, which needs to be corrected for in order to estimate the intrinsic stellar Lyman-α line profile. Knowledge of the intrinsic line is useful for several reasons: it is proportional to the stellar radiation pressure, which is essential in understanding the structure of a planetary exosphere (Sect. 3.2); intrinsic lines in different epochs can be directly compared, even when they are derived from observations obtained with different instrumental settings (Sect. 3.2); and empirical relations between the intrinsic Lyman-α flux and the stellar EUV spectrum can be used to estimate the irradiation of the planet (Sect. 4.2.1).

3.1. Reconstruction method

The intrinsic Lyman-α line can be reconstructed directly from the STIS spectra, using the method detailed in Bourrier et al. (2015a; see, e.g., Wood et al. 2005; Ehrenreich et al. 2011; Youngblood et al. 2016). A model profile of the stellar Lyman-α line is absorbed by the interstellar hydrogen and deuterium, and convolved by the STIS line spread function (LSF). A Metropolis-Hasting Markov chain Monte Carlo algorithm is used to find the best fit between the model spectrum and the observed spectrum, which are compared over the entire spectral line (see the spectral ranges in Table 4). The medians of the posterior probability distributions are chosen as the final best-fit values for the model parameters, and their 1σ uncertainties are obtained by finding the intervals on both sides of the median that contain 34.15% of the accepted steps.

In addition to the reference spectra calculated in Sect. 2.3 for Visits 1, 2, and 3, we used the measurement of HD 97658 Lyman-α line secured by France et al. (2016) in February 2015, about a quarter orbital period after the transit (hereafter, Visit 0). Because it was obtained with STIS E140M mode this spectrum could not be directly compared with our G140M measurements in the previous sections. Nonetheless, it can be used to reconstruct the intrinsic line at another epoch, and the higher resolution of the echelle mode combined with the use of the 0.2′′ × 0.06′′ slit significantly limits the airglow contamination and enables the measurement of the ISM thermal broadening, constrained by the ISM deuterium absorption (Fig. 6).

The Lyman-α line can be double-peaked (e.g., the K dwarf HD 189733; Bourrier et al. 2013) or single-peaked (e.g., the M dwarf GJ 436; Ehrenreich et al. 2015), depending on the presence of a self-reversal arising from the absorption by the upper chromospheric layers of the Lyman-α emission from the lower chromosphere. We thus tried different combinations of Gaussian and Voigt profiles to model the possible presence of a self-reversal (Bourrier et al. 2015a; Youngblood et al. 2016). The fit of the STIS LSF is included in the reconstruction process because there is no published tabulated LSF profile for the 50′′ × 0.05′′ slit used in Visits 1 to 3, and because it can be contaminated by systematic noise sources (e.g., stray light). Past studies (e.g., Ehrenreich et al. 2012; Bourrier et al. 2013; Bourrier et al. 2015b) have shown that the fitted LSF can be wider than the 50′′ × 0.1′′ slit LSF, and can possess broad, shallow wings. We thus modeled the LSF as a Gaussian core, and tried including an additional Gaussian representing these wings. We used the Bayesian Information Criterion (BIC; de Wit et al. 2012) to prevent over-fitting in our search for the best model for the Lyman-α and LSF profiles.

3.2. ISM and line properties

In a first step we performed the reconstruction independently for each visit, and found that the lowest χ² and BIC was always obtained when the intrinsic line is modeled as a single-peaked Voigt profile. Surprisingly, Bourrier et al. (2013) found that the intrinsic line of HD 189733, which has the same stellar type as HD 97658, was best reconstructed by including a self-reversal. Furthermore, the Mg II line of HD 97658 (France et al. 2016) does present a self-reversed emission. The absence of a self-absorption at the peak of HD 97658 Lyman-α line would be related to the temperature and density structure of neutral hydrogen in the upper chromosphere of the K dwarf (e.g., Curdt et al. 2001).

In Visit 0, a single Gaussian best fits the echelle mode LSF, while in the three visits obtained with the G140M the LSF is best modeled with a double-Gaussian profile. The free parameters of the fit are thus, for each epoch, the heliocentric velocity of the star γ_∗, the Lyman-α line profile parameters (maximum stellar flux at 1 au f_peak(1 au), Doppler width Δv_D and damping parameter a_damp), and the LSF profile parameters (the width of the core , and when required the width of the wings and ratio between the peaks of the core and the wings ). The ISM parameters for neutral hydrogen are the column density log₁₀ N(H i), Doppler broadening parameter b(H i), and velocity relative to the star γ(H i)_/∗. The LISM Kinematic Calculator¹ (Redfield & Linsky 2008) predicts that the line of sight (LOS) toward HD 97658 passes within 20° of the LIC, Leo, and NGP clouds. We tried including two ISM components in the reconstruction process, but found that the second component could not be well constrained. We thus modeled the ISM opacity from a single component as the combination of two Voigt profiles for the atomic hydrogen and deuterium, separated by about 0.33 Å. We used a D i/H i ratio of 1.5 × 10^-5 (e.g., Hébrard & Moos 2003; Linsky et al. 2006) and assumed there is no turbulent broadening (implying that b(D i) = b(H i)/, Wood et al. 2005). Assuming that the ISM kept the same properties over the year of observations of HD 97658, we used the weighted mean of the spectra of Visits 1, 2, and 3 as a reference to estimate the ISM column density and radial velocity (Visit 0 was not included because of its different instrumental response). Shifted to the heliocentric rest frame, the derived ISM velocity (7.2 ± 0.3 km s^-1) is in between the predicted velocities for the Leo and NGP clouds (4.7 ± 0.8 km s^-1 and 9.5 ± 0.8 km s^-1). Based on the Visit 0 spectrum alone, Youngblood et al. (2016) found a degeneracy between a low b(H i) value (6 km s^-1), high column density (18.45) solution and a high b(H i) value (12 km s^-1), low column density (18.2) solution. We thus confirm their expectation that the second solution is more realistic, and that the ISM component along HD 97658 LOS has a high Doppler broadening (13.1 km s^-1) and a low column density (18.18). The parameters log₁₀ N(H i) and γ(H i)_/∗ were set to these best-fit values in the final reconstructions. ISM Doppler broadening b(H i) cannot be constrained by the low-resolution G140M spectra and was estimated through the reconstruction of the Visit 0 echelle spectrum, along with the LSF and intrinsic line properties at this epoch. We set b(H i) to its derived best-fit value and reconstructed the intrinsic line for Visits 1, 2, and 3 independently.

The best-fit reconstructed Lyman-α stellar profiles are displayed in Fig. 6 and the corresponding model parameters are given in Table 4. The best-fit LSF is stable in the three G140M visits. The quality of the fits is good in all epochs, with reduced χ² lower than 1.1. For Visits 2 and 3 the reduced χ² were found to be lower than 1, and the uncertainties on the parameters derived from their fits can be considered as conservative. These results confirm that the mean spectra used as reference (Sect. 2.3) trace the intrinsic stellar line well. We do not detect significant variations between the four different epochs². Nonetheless, it is noteworthy that the median value for the peak amplitude of the line is higher in Visits 0 and 1, separated by only 2 months, and that the line-integrated flux is about 2σ higher in Visit 1 than in later visits distant by more than 8 months (Table 4; Fig. 6). This is consistent with the higher flux measured at the peaks of the observed line in Visit 1 (Sect. 2.3.1). Its origin is likely an increased emission in the line core from a more active upper stellar chromosphere.

Finally, we note that the intrinsic Lyman-α line of HD 97658 is strong enough that hydrogen atoms escaping from the planet would be accelerated away from the star by radiation pressure. Similarly to HD 189733, this force can indeed overcome stellar gravity up to a factor of ~3 (Fig. 6), although the different shapes of the intrinsic lines and the larger orbital distance of HD 97658 b suggests that the dynamics of a putative exosphere would be different from that of HD 189733 b.

4. XUV irradiation of the upper atmosphere

4.1. Measurement of the X-ray stellar emission

We observed HD 97658 (PI: P.J Wheatley) in June 2015 with XMM-Newton/EPIC-PN (Obs ID: 0764100601; duration 33.9 ks), and contemporaneously with the Lyman-α observations in December 2015 (Obs ID: 18724; duration 19.82 ks) and March 2016 (Obs ID: 18725; duration 19.82 ks) with Chandra/ACIS-S (Table 2). The XMM EPIC-PN spectrum shows HD 97658 to be about ten times fainter and much softer than HD 189733, even though both stars are about at the same distance (~20 pc) and have the same stellar type. This indicates a significantly lower coronal activity for HD 97658, possibly due to its older age (9.7 ± 2.8 Gyr, compared to 5.3 ± 3.8 Gyr for HD 189733; Bonfanti et al. 2016). We use the flux measured in 0.15−2.4 keV with XMM EPIC-PN, and 0.234−2.4 keV with Chandra to fit the X-ray spectrum of the star at the different epochs of observations. We used C-statistics to fit spectral models to the data, a more reliable method than χ² statistics for low count rates. The XMM EPIC-PN spectrum was thus binned to a minimum of 10 counts per bin, while the Chandra spectra were set to a minimum of 3 counts per bin. In a first step we tested the merits of three different spectral models by fitting only the XMM EPIC-PN spectrum. The first model (APEC; Smith et al. 2001), generally used for X-ray observations of coronal emission from late-type stars, is a single-temperature model of optically thin thermal plasma. Here we use two single-temperature components in order to approximate a multi-temperature plasma. The second (CEMEKL; Schmitt et al. 1990; Singh et al. 1996) and third (C6PMEKL; Lemen et al. 1989; Singh et al. 1996) models both generate a spectrum from a continuous range of temperatures, with the emission measure distribution set respectively by a power law and a Chebyshev polynomial. The C6PMEKL model differs from the CEMEKL model by including the possibility for higher order terms. Fluxes derived from the three models showed good agreement in the observed energy band 0.2−2.4 keV in which the PN camera is well calibrated. However, we need to extrapolate the X-ray spectrum down to 0.1 keV (124 Å) to derive the stellar EUV flux with the relation from Chadney et al. (2015). For our two-temperature APEC model fit there is no constraint on the relative temperatures and emission measures of the two components, and this seems to allow too much freedom when extrapolating the fitted spectrum, with a very large flux implied between 0.1 and 0.2 keV. We regard this extrapolation as unreliable. The emission measure distribution in the other two models is more tightly constrained and probably represents more closely the true emission measure distribution of the star, and yields more consistent and plausible results down to 0.1 keV. Since the goodness of fit is marginally better for the simpler CEMEKL model, it was used to fit the entire X-ray dataset.

We carried out joint fits of the XMM and Chandra spectra, forcing the spectral profile to be the same but allowing the flux level to vary between observations. We note that the Chandra ACIS instrument is accumulating a contaminant on the filter, which is reducing the soft-band sensitivity. Nonetheless, we used the most recent update of the Chandra calibration pipeline (CALDB 4.6.2; July 2014), which provides a realistic model for the contamination, and we checked that including soft energies in the simultaneous fit to the XMM and Chandra spectra did not change significantly the fit to the high energies alone. This gives us confidence in the derived Chandra fluxes, which were found to be consistent with each other but are up to twice as low as the XMM flux (at 2−2.5σ). Indeed, in the common Chandra band (0.243–2.4 keV) we found fluxes at Earth distance of 2.1810^-14 erg cm^-2 s^-1 in June 2015 (XMM), 1.3 ± 0.4 × 10^-14 erg cm^-2 s^-1 in December 2015 (Chandra), and 1.1 ± 0.4 × 10^-14 erg cm^-2 s^-1 in March 2016 (Chandra). The stability of the Chandra measurements is consistent with that of the contemporaneous Lyman-α measurements in Visits 2 and 3. The spectra measured in each epoch are shown in Fig. 7 along with their best fits that yield extrapolated fluxes between 0.1−2.4 keV (5.17−124 Å ) of 14.6 erg cm^-2 s^-1 in June 2015, 8.9 erg cm^-2 s^-1 in December 2015, and 7.2 erg cm^-2 s^-1 in March 2016.

4.2. XUV irradiation

4.2.1. XUV spectrum

We calculated the XEUV spectrum of HD 97658, from 5 Å  to 912 Å , to infer the photoionization rate of neutral hydrogen in the upper atmosphere of the planet and to constrain its energy-limited mass-loss rate. The EUV emission is mostly absorbed by the ISM and must be estimated through indirect methods.

In a first step, we estimated the flux in 124–912 Å  for each epoch of observation using two independent methods. The EUV spectrum was estimated with the empirical relation from Chadney et al. (2015), which is based on the high-energy X-ray flux measured in 5.17−124 Å  (Sect. 4.1), and with the semi-empirical relations from Linsky et al. (2014), which are based on the lower energy integrated intrinsic Lyman-α flux (Table 4). We obtained the following fluxes at 1 au from the star (in erg cm^-2 s^-1): 4.26 (February 2015, from HST Visit 0), 4.87 (April 2015, from HST Visit 1), 6.57 (June 2015, from XMM), 4.44 and 4.93 (December 2015, from HST Visit 2 and Chandra, respectively), 4.37 and 4.36 (March 2016, from HST Visit 3 and Chandra, respectively). Estimating uncertainties on these EUV fluxes is not straightforward, as they depend on uncertainties on the X-ray and the Lyman-α measurements, and on the empirical relations themselves. The median of the seven derived values yields a flux of 4.44 erg cm^-2 s^-1, and a rms dispersion of 19%. Within this uncertainty we did not detect any trend in the EUV flux as a function of time, and although the stellar emission is higher in June 2015 it remains within 2.5 rms of the median flux. The EUV emission of HD 97658 can thus be considered stable within 20% variations, which shows not only that the Chadney et al. (2015) and Linsky et al. (2014) formulae can provide consistent first-order predictions of the EUV emission, but also that a flux of 4.44 ± 0.84 erg cm^-2 s^-1 is a good estimate of HD 97658 average emission in 124−912 Å.

Finally, we used the Linsky et al. (2014) relations to calculate the EUV spectrum in complementary wavelengths bands between 100 and 1170 Å , adjusting the global level of the fluxes so that they yield the value derived in 124–912 Å . For each band, we set an uncertainty of 19%. The stellar emission in 5–100 Å  was derived from the spectra obtained in Sect. 4.1, and fixed to the median value of the three epochs. Flux values for the entire XUV spectrum are given in Table 5.

4.2.2. Photoionization rate

The photoionization rate per neutral hydrogen atom at the semi-major axis of the planet is defined as (2)with Γ_ion in s^-1, F_XUV(λ) the stellar flux at 0.08 au in erg s^-1 cm^-2 Å^-1 and, σ_ion the cross-section for photoionization in cm². The integration is performed up to the ionization threshold at 911.8 Å. The cross-section for hydrogen photoionization is wavelength-dependent, and we used the expression from Verner et al. (1996) and Bzowski et al. (2013), (3)with σ_ion in cm² and wavelengths in Å. Using the XUV spectrum given in Table 5 and propagating its uncertainties, we obtained a photoionization rate at the semi-major axis of HD 97658 b of 4.8 s^-1 (3.1 s^-1 at 1 au from the star). This corresponds to a photionization lifetime of about 5.8 ± 1.8 h at the semi-major axis.

Bourrier & Lecavelier des Etangs (2013) found ionizing fluxes between 6 and 23 times the solar value in their study of HD 189733 b. Using their conventions, this converts into photoionization rates at 1 au in [4.6−17.7] × 10^-7 s^-1, higher than the photoionization rates from HD 97658 at the same distance and consistent with our aforementioned proposition that HD 189733 has a more active upper chromosphere and corona. Since the hot Jupiter HD 189733 b orbits about 2.5 times closer than HD 97658 b from its host star, the photoionization rate at HD 189733 b is at least six times higher than at HD 97658 b. In contrast, Bourrier et al. (2016) found that the evaporating warm Neptune GJ 436 b, in orbit around a colder M dwarf, is subjected to low photoionization rates in [1.4−3.4] × 10^-5 s^-1 at the planet, explaining in part the formation of the giant exosphere of neutral hydrogen surrounding the planet (Ehrenreich et al. 2015). HD 97658 b is thus closer to GJ 436b in terms of photoionization, while being subjected to more than four times its radiation pressure (see Bourrier et al. 2015a, and Fig. 6 in the present paper). This strong radiation pressure would swiftly blow away neutral hydrogen atoms escaping HD 97658 b, which would then remain in this state for many hours owing to the low photoionization rates, leading eventually to the formation of a comet-like tail with a very long extension. We investigate in the next sections reasons that could explain our non-detection of such an exosphere.

5. Observational constraints on the upper atmosphere

5.1. EVE simulations of HD 97658 b

We use the non-detection of neutral hydrogen absorption from HD 97658 b (Sect. 2.3.2) and the measurement of its radiation environment (Sect. 4.2) to constrain the structure of its upper atmosphere. We interpret the transit observations of HD 97658 b in the Lyman-α line using the EVaporating Exoplanets (EVE) code (for a detailed description see, e.g., Bourrier et al. 2015a). The upper atmosphere is described in 3D with two different regions: the bottom layers are parameterized with a fixed set of atmospheric profiles, while Monte Carlo particle simulations are used to compute the dynamics of neutral hydrogen meta-particles in the exosphere. The transition between these two regions is set at R_trans, which represents the altitude up to which the atmosphere keeps a global cohesion. For the atmosphere below the transition, hereafter called the isotropic atmosphere, we assume an upper atmosphere of pure neutral hydrogen with an isotropic, hydrostatic density profile. This profile is constrained by a mean temperature T_iso and scaled so that it matches the exospheric hydrogen density at R_trans. To represent the expansion of the hydrogen outflow we set a radial velocity profile derived from Fig. 3 in Koskinen et al. (2013a), which increases from a few meters per second at the planet surface to a velocity v_trans at R_trans. We note that the choice of this velocity profile does not have much influence on our results compared to the value of v_trans. At each time step, meta-particles representing groups of neutral hydrogen atoms with the same properties are released from the entire atmosphere at R_trans with a rate Ṁ_H⁰. Because the orbit is slightly eccentric (e = 0.078), the escape rate varies as the inverse square of the distance to the star, with a reference set at the semi-major axis. The initial velocity distribution of the particles relative to the planet is the combination of the outflow bulk velocity v_trans, and an additional thermal speed component from a Maxwell-Boltzman velocity distribution defined by T_iso. The dynamics of the particles is calculated in the stellar reference frame and results from the stellar and planetary gravities, the velocity-dependent stellar radiation pressure, and the inertial force linked to the non-Galilean stellar reference frame. Neutral hydrogen particles are also subjected to stellar photoionization, set by the rate 4.8 × 10^-5 s^-1 at the planet semi-major axis (as derived in Sect. 4.2.2) but varying with the particles’ distance to the star. The impact of radiation pressure and photoionization onto neutral hydrogen particles is calculated taking into account self-shielding within the exospheric cloud.

The EVE code generates theoretical Lyman-α spectra as seen with HST/STIS after absorption by the planet, its upper atmosphere, and the interstellar medium using a limited number of parameters. The direct comparison of these spectra with observational data allows us to study the influence of R_trans (in R_pl), Ṁ_H⁰ (in g s^-1), v_trans (in km s^-1), and T_iso (in K). We create a single observational dataset by combining the spectra from Visits 2 and 3³, which were shown to be very similar (Sect. 2.3.1). The exposures from these two visits do not overlap in terms of planetary orbital phase, except for the first orbit of Visit 3 and the last orbit of Visit 2, which were merged into a single exposure. These observations were compared to theoretical Lyman-α line spectra calculated every 5 mn with EVE at a resolution Δλ = 0.02 Å  corresponding to Δv = 5 km s^-1 (about a quarter of the resolution of the STIS spectra at 1215.67 Å). The intrinsic line flux required to calculate the theoretical spectra and radiation pressure on hydrogen atoms was set to the average line between Visits 2 and 3 (Fig. 6). The merit function is the sum of the χ² yielded by the comparison of the observed spectra with the theoretical spectra averaged during the time window of each observation, convolved with STIS LSF (set to the average reconstruction from Visits 2 and 3), and limited to the wavelength range corresponding to the velocity range [−300; 300] km s^-1.

5.2. Structure of the upper atmosphere

The model atmosphere can be constrained in two ways by the observations. Neutral hydrogen in the isotropic atmosphere contributes mainly to absorption in the core of the Lyman-α line, which cannot be observed because of ISM absorption (Sect. 2.3.1). If the column density is high enough, the isotropic atmosphere can influence the observed spectra through the Lorentzian wings of its absorption profile (the so-called damping wings). However, this also requires a significant atmospheric expansion because damping wings can only absorb as much as the opaque surface of the atmosphere projected onto the stellar disk. It is highly unlikely that the atmosphere of a super-Earth like HD 97658 b would be dense within a large enough volume to produce detectable damping wings. Even if that were the case and the upper atmosphere extended up to the Roche lobe (~11.8 R_p), its transit would still last for about the same duration as the optical transit (~3 h). In contrast, the strong radiation pressure from the host star (Sect. 3.2) shapes the exosphere into an extended cometary tail trailing behind the planet (Fig. 9). We estimate from our simulations that the escaping neutral hydrogen can be accelerated by radiation pressure up to about 95 km s^-1 away from the star, which implies that the exosphere can produce absorption mainly in the blue wing of the Lyman-α line after the optical transit. While the bulk motion of the gas in the exosphere is set by the combination of radiation pressure and stellar gravity, Fig. 9 illustrates how its spatial extension in the transverse direction increases with increasing outflow velocity (v_trans) and thermal dispersion (T_iso).

Figure 8 shows upper limits on the escape rate of neutral hydrogen as a function of the mean temperature of the isotropic atmosphere and the outflow velocity. The transit observations of HD 97658 b are most consistent with the absence of an extended upper atmosphere, but they do not completely exclude a moderate atmospheric escape. With high values for the outflow temperature and bulk velocity, the escaping neutral hydrogen is diluted enough for the exosphere to be optically thin and its absorption roughly proportional to the local density. Since exospheric density scales at first order with the escape rate, the non-detection of an absorption signature from HD 97658 b limits the escape rate roughly independently of the other model parameters for temperatures above ~4000 K (Fig. 8). Conversely, the entire exosphere is more compact when neutral hydrogen escapes from low altitudes and/or with a lower outflow velocity and temperature (Fig. 9). In that case regions close to the planet become optically thick and their absorption increases more slowly with increasing escape rates, allowing for larger mass loss to be consistent with the observations (Fig. 8). Finally, we note that the upper limit on the escape rate drops quickly when the temperature decreases below ~1000 K. This occurs because the density profile in the isotropic atmosphere is scaled by its value at the transition R_trans, which depends on the escape rate. Since the density profile then increases steeply closer to the planet with low temperatures, even moderate escape rates can lead to high column density and the formation of strong damping wings at such temperatures, which are excluded by the observations. This effect is stronger when the size of the isotropic atmosphere increases as its larger projected area allows for larger damping wings.

Salz et al. (2016a) performed 1D, spherically symmetric simulations of HD 97658 b and predicted that the upper atmosphere of the planet should expand hydrodynamically, reaching velocities between ~1 and 10 km s^-1 and a mean temperature of about 3500 K above 2 R_p. In these conditions, our simulations set an upper limit of ~10⁸ g s^-1 on the neutral hydrogen mass loss, which remains unchanged for larger temperature or stronger atmospheric expansion (Fig. 8). Since it is unlikely that the upper atmosphere of HD 97658 b could expand with lower temperatures Salz et al. 2016a, we consider that the upper limit on the neutral hydrogen mass loss from HD 97658 b is 10⁸ g s^-1 at the 3σ level.

5.3. Mass loss, heating efficiency, and neutral fraction

Assuming that the upper atmosphere of HD 97658 b is in the energy-limited regime (see, e.g., Owen & Alvarez 2016), we calculated its total mass-loss rate as (4)with η the fraction of the energy input that is available for atmospheric heating (e.g., Ehrenreich & Désert 2011; Lammer et al. 2013; Shematovich et al. 2014), R_XUV the mean altitude at which the XUV energy is absorbed in the upper atmosphere (e.g., Lammer et al. 2003), F_XUV(sma) the total XUV flux per unit area at the semi-major axis (sma) of HD 97658 b, ρ_p the density of the planet, and K_tide a correction factor accounting for the contribution of tidal forces to the potential energy (Erkaev et al. 2007)⁴. Using our measurement for the XUV flux (835 erg s^-1 cm^-2, Table 5), we found  g s^-1. Defining Ṁ_H⁰ = f_H⁰Ṁ^tot, where f_H⁰ is the fraction of neutral hydrogen in the escaping gas and Ṁ_H⁰ the upper limit on the escape rate of neutral hydrogen obtained with EVE, we can derive that . This constraint hints at a low heating efficiency and/or low neutral hydrogen fraction in the upper atmosphere. Assuming R_XUV ~ R_p and using theoretical estimations of η between 10 and 30% (e.g., Lammer et al. 2013; Shematovich et al. 2014; Owen & Alvarez 2016; Salz et al. 2016b), the upper limit on the neutral hydrogen fraction would range between about 12 and 36%.

An independent constraint can be derived from the total mass-loss rate Ṁ^tot = 3 × 10⁹ g s^-1 predicted by Salz et al. (2016a). They used an XUV irradiation of 955 erg s^-1 cm^-2, which is consistent with the value we measure. However, the comparison between their Ṁ^tot and our 3σ upper limit on the neutral hydrogen escape rate Ṁ_H⁰ yields %. This value is surprising considering the low irradiation and photoionization rate derived from our observations (Sect. 4.2.2), and is inconsistent with the high neutral hydrogen fraction (>50%) obtained by Salz et al. (2016a). This suggests either that the upper atmosphere of HD 97658 b does not expand hydrodynamically, or that its thermospheric and mass-loss properties are much different from the ones predicted by Salz et al. (2016a).

6. Discussion

Despite appearing quiet in the optical and NIR (Van Grootel et al. 2014), HD 97658 has detectable high-energy emission. It is a weak and soft X-ray source with a variable upper chromosphere and corona. Our HST and XMM-Newton/Chandra observations of HD 97658 at five different epochs reveal signs of both short-term and long-term variability in the Lyman-α and X-ray emission of this star. We were nonetheless able to determine a stable, average reference for the intrinsic Lyman-α line and XUV spectrum of the star, refining the properties derived by Youngblood et al. (2016). These estimates will be useful in determining the impact of the incident high-energy radiation on the atmospheric heating and chemistry for this super-Earth, and more generally in understanding the conditions that lead to evaporation. The HD 97658 system might be crucial in understanding the stability of the atmospheres of small planets in terms of atmospheric composition, planetary density, and irradiation. Indeed, HD 9758 b is a super-Earth more than twice as dense (3.9 g cm^-3) as the evaporating warm Neptune GJ 436 b (1.6 g cm^-3). Furthermore, HD 97658 b orbits at 0.08 au from its K dwarf host star and we measure an XUV energy input of 835 ± 160 erg s^-1 cm^-2 (Sect. 4.2), which is about three times lower than the irradiation of GJ 436b (at 0.027 au from its M dwarf host; Ehrenreich et al. 2015; Bourrier et al. 2016). The upper atmosphere of HD 97658 b is also subjected to a photoionization rate of 4.8 s^-1, which is more than six times lower than for the evaporating hot Jupiter HD 189733 b (at 0.031 au from its K dwarf host; Bourrier & Lecavelier des Etangs 2013).

Our analysis of three Lyman-α transits of the super-Earth HD 97658 b at independent epochs yields no signature arising from an extended atmosphere of neutral hydrogen. This suggests that the thermosphere of HD 97658 b is not hydrodynamically expanding, and in any case that the escape of neutral hydrogen is below 10⁸ g s^-1 at 3σ. This might be linked to the low XUV stellar irradiation, to an inefficient conversion of the stellar energy input, and/or to a low hydrogen content in the upper atmosphere (possibly linked to a helium enrichment, Hu et al. 2015, or to an enhanced mass loss from the irradiation and stellar wind of the young host star). A more exotic possibility would be that HD 97658 b is in a saturation regime with the wind from its host star, as defined by Bourrier & Lecavelier des Etangs (2013), where the stellar wind is so dense that it interacts with all neutral hydrogen atoms escaping the upper atmosphere. Provided that the stellar wind protons move faster than about 300 km s^-1, all hydrogen atoms escaping HD 97658 b would then be ionized through charge-exchange, while the high-velocity neutralized protons would contribute to Lyman-α absorption too far from the line core to be detectable (Bourrier et al. 2016).

The methodology used in this study is the same as led to the non-detection of Lyman-α signatures around the super-Earth 55 Cnc e and to the detection of hydrogen exospheres around the giant planets HD 209458 b (Vidal-Madjar et al. 2003), HD 189733 b (Lecavelier des Etangs et al. 2010, 2012), 55 Cnc b (Ehrenreich et al. 2012), and GJ 436b (Ehrenreich et al. 2015). There is a clear distinction between the stellar activity of HD 97658 b, revealed as short- and long-term variations in the core of the Lyman-α line, and the transit of a neutral hydrogen exosphere. The latter occurs around and after the time of the optical transit, with the depth and spectral range of the absorption evolving coherently over time. Confirmed signatures of hydrogen escape have been detected within specific wavelength ranges in the blue wing of the Lyman-α line, which can be explained well by the physical mechanisms acting on the upper planetary atmosphere. In this frame our study of HD 97658 b can be seen as a control experiment; the results we obtained confirm that detections of atmospheric escape are not biased by our methodology, and underline its strength and reliability. Future observations of this system in other FUV lines will be needed to search for escape signatures of heavier species than hydrogen and to further constrain the nature of this super-Earth.

---

Acknowledgments

We thank the referee, Jeffrey Linsky, for his comments that helped improve the science and clarity of this paper. We warmly thank Diana Dragomir for her validation of the HD 97658 b transit ephemeris. This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The authors thank the STScI for extending the proprietary period of program 13820, allowing a careful analysis of the full dataset. This work has been carried out in the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). V.B. and D.E. acknowledge the financial support of the SNSF. A.L.E. acknowledges financial support from the Centre National d’Études Spatiales (CNES). The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under program ANR-12-BS05-0012 “Exo-Atmos”.

References

All Tables

All Figures