Free Access
Issue
A&A
Volume 638, June 2020
Article Number A49
Number of page(s) 21
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201834814
Published online 10 June 2020

© ESO 2020

1 Introduction

The first observational evidence of an expanding and escaping exoplanetary atmosphere was the Lyα transit absorption detected for the hot Jupiter HD 209458b (Vidal-Madjar et al. 2003). Subsequently, more observations of excess absorption in Lyα and UV lines of other elements succeeded, the most prominent ones being the hot Jupiters HD 209458b (Vidal-Madjar et al. 2004, 2008, 2013; Ben-Jaffel 2007, 2008; Ballester et al. 2007; Ehrenreich et al. 2008; Linsky et al. 2010; Ben-Jaffel & Sona Hosseini 2010; Ballester & Ben-Jaffel 2015), HD 189733b (Lecavelier des Etangs et al. 2010, 2012; Bourrier et al. 2013; Ben-Jaffel & Ballester 2013), WASP-12b (Fossati et al. 2010; Haswell et al. 2012), 55 Cnc b (Ehrenreich et al. 2012), and the warm Neptune GJ 436b (Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017). Numerous hydrodynamic models for hot Jupiter upper atmospheres were developed, aiming to explain the observations (Yelle 2004; Tian et al. 2005; García Muñoz 2007; Penz et al. 2008; Stone & Proga 2009; Murray-Clay et al. 2009; Guo 2011, 2013; Koskinen et al. 2013; Shaikhislamov et al. 2014; Salz et al. 2016a; Guo & Ben-Jaffel 2016; Erkaev et al. 2016; Khodachenko et al. 2017; Debrecht et al. 2018). Other studies employed Direct Simulation Monte Carlo (DSMC) models, which include the generation of energetic neutral atoms (ENAs) due to interaction with the stellar wind and take into account the effects of radiation pressure on the shaping of the large hydrogen clouds (Holmström et al. 2008; Ekenbäck et al. 2010; Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014; Bourrier et al. 2015, 2016).

It is therefore important to study the interaction of the expanding planetary atmosphere with the magnetized stellar wind. The escaping atmospheric particles mixed with stellar wind plasma can have a strong influence on the wind plasma parameters in the vicinity of the planet. The main effect of an intrinsic planetary magnetic field on atmospheric escape is to suppress the outflow and to makeit highly anisotropic (Adams 2011; Trammell et al. 2011, 2014; Owen & Adams 2014; Khodachenko et al. 2015). Many studies investigated the interaction between a close-in planet and its host star’s wind, but some of them neglected magnetic fields and just considered a purely hydrodynamic interaction (Stone & Proga 2009; Bisikalo et al. 2013; Tremblin & Chiang 2013; Christie et al. 2016). Other studies applied magnetohydrodynamic (MHD) models instead (Cohen et al. 2011; Matsakos et al. 2015; Tilley et al. 2016), but employed mostly simplified descriptions of the planetary wind. Shaikhislamov et al. (2016) used a 2D multi-fluid code to study the interaction of a non-magnetized HD 209458b-like hot Jupiter with the stellar wind, taking into account heating by the stellar X-ray and extreme ultraviolet (XUV) flux and hydrogen photochemistry to self-consistently model the planetary outflow. However, they did not include the interplanetary magnetic field (IMF) and its effect on the formation of the planetary obstacle. This was addressed by Erkaev et al. (2017), who used a 1D hydrodynamic planetary upper atmosphere model in combination with a 3D MHD stellar wind flow model to investigate the build-up of a planetary obstacle by the interaction of the partially ionized planetary wind with the plasma flow of a magnetized stellar wind.

The hot Jupiter HD 189733b was discovered by Bouchy et al. (2005) based on both radial velocity and photometric transit observations. The most recent determinations of its mass and radius are 1.13 MJup and 1.13 RJup, based on parallaxes from the first data release of Gaia (Stassun et al. 2017). The host star HD 189733 has a spectral type of K2V (Gray et al. 2003) and is the primary of a double system with a mid-M dwarf companion located at a projected separation of ~216 AU (Bakos et al. 2006). The activity level and rotation period (11.95 d; Henry & Winn 2008) of HD 189733, comparable to the similar K star ɛ Eri, indicate an age of 1-2 Gyr (Poppenhaeger & Wolk 2014). However, the accompanying M dwarf is rather inactive and has an estimated age of > 5 Gyr. It was therefore suggested that HD 189733’s high rotation rate and associated strong magnetic activity may be caused by interactions with its hot Jupiter (Poppenhaeger & Wolk 2014). Independent of the reasons for the star’s high activity level, HD 189733b is exposed to the intense stellar XUV emission, and possibly also to a dense and fast stellar wind. This may lead to enhanced atmospheric losses, which is also indicated by UV transit observations. Moreover, Pillitteri et al. (2015) suggested that the observed FUV variability could stem from accretion of matter from the planet onto the star. Route (2019) analyzed multiwavelength data of HD 189733 to find possible indications of star-planet interaction, but could not identify a clear relation between the stellar activity characteristics and the planetary orbital phases.

Lecavelier des Etangs et al. (2010) observed three planetary transits and detected an absorption of 5.05 ± 0.75% in the unresolved Lyα line, significantly higher than the absorption by the planetary disk at optical wavelengths. Later, Lecavelier des Etangs et al. (2012) obtained Lyα observations with higher spectral resolution at two different epochs (April 2010, September 2011). While during the first epoch no excess absorption could be detected (just absorption of the total flux comparable to the 2.4% absorption by the planetary disk at optical wavelengths, but no spectrally resolved absorption), they found absorption of 14.4 ± 3.6% (i.e., 12.3 ± 3.6% excess absorption) in the blue wing of the line in a velocity range of − 230 to − 140 km s−1, indicating absorbing material moving away from the star. They also found absorption in the red wing from 60 to 110 km s−1 of 7.7 ± 2.7%, meaning 5.5 ± 2.7% excess, indicating absorbing material moving towards the star, but this detection was not statistically significant. Using Monte Carlo simulations, Bourrier & Lecavelier des Etangs (2013) were able to reproduce the excess absorption observed in 2011 with a neutral hydrogen outflow rate of 4 × 108−1011 g s−1, an ionizing flux of 6–23 times the solar value, a stellar wind with a temperature of 3 × 104 K, a velocity of 200 ± 20 km s−1, and a density between 103 − 107 cm−3. Lecavelier des Etangs et al. (2012) suggested that the absence of excess absorption in 2010 could be due to a much lower escape rate or a less dense stellar wind. They explained the discrepancy between the observations in 2010 and 2011 by the influence of a flare that occurred ~8 h prior to the transit in 2011. Bourrier et al. (2013) presented a more detailed analysis of the 2011 observations, reaching similar conclusions.

Further attempts were made to explain the temporal variability of the Lyα transit absorption of HD 189733b. Guo & Ben-Jaffel (2016) investigated the effect of the stellar XUV spectral energy distribution (SED) on atmospheric profiles and mass-loss rates. They found that the mass-loss rate is mainly determined by the total absorbed energy, whereas the ionization is strongly affected by the SED. For SEDs dominated by the low-energy part of the spectrum, the H/H+ transition moves closer to the planet, and the amount of H atoms at a certain altitude can differ by 1–2 orders of magnitude, in comparison to SEDs dominated by the high-energy part. They used the method of Ben-Jaffel (2008) to investigate the differences in absorption signature depending on the stellar XUV SED and found that they can explain the differences in observationsbetween 2010 and 2011 by assuming a harder stellar spectrum in 2011. The model of Ben-Jaffel (2008) assumes an extended thermosphere with an absorption profile broadened by natural broadening and does not include nonthermal H populations like ENAs. Therefore, the resulting transmission spectrum is symmetric around the Lyα line center and does not reproduce the localized absorption at higher velocities in the blue wing. Recently, Chadney et al. (2017) studied the influence of flares on the upper atmospheres and escape rates of hot Jupiters. They found a maximum mass-loss enhancement of about a factor of two, and much less if the limited duration of the radiation enhancement during a flare is taken into account. However, they suggested that an extreme proton event associated with the flare could have led to sufficiently enhanced escape.

Ben-Jaffel & Ballester (2013) reported a 6.4 ± 1.8% absorption in O I and a marginal detection of early-ingress C II absorption. However, they could not exclude that the latter had a stellar or instrumental origin. Transit absorption in the Hα line was also detected in several observations, as well as a pre-transit signature (Jensen et al. 2012; Cauley et al. 2015, 2016, 2017b). However, both in- and pre-transit absorption signatures were highly variable in time and strongly affected by stellar activity, making their interpretation difficult (Barnes et al. 2016; Cauley et al. 2017a; Kohl et al. 2018). Furthermore, a 6–8% transit absorption in X-rays was reported by Poppenhaeger et al. (2013), but more observations with higher sensitivity are needed to confirm this result (Marin & Grosso 2017).

In this paper, we investigate the possible causes for the variations of the Lyα transit absorption of HD 189733b. To achieve this, we modeled the structure of the upper atmosphere and the associated planetary mass loss, taking into account stellar XUV heating at both quiescent and flaring conditions. Then, we modeled the interaction between the stellar wind and the planetary atmosphere, including the related production of ENAs, and its effect on the UV transit signature. In Sect. 2, we describe the hydrodynamic model and show the resulting upper atmosphere profiles and planetary mass-loss rates in Sect. 3. In Sect. 4, we describe the MHD flow model and show the stellar wind flow maps. In Sect. 5, we compute the Lyα absorption for the various scenarios and compare them with the observations. In Sect. 6, we discuss our findings and compare them with other studies. Conclusions are presented in Sect. 7.

2 Hydrodynamic upper atmosphere modeling

2.1 Model description

The hydrodynamic model solves the time-dependent system of the equations of mass, momentum and energy conservation in 1D spherical geometry along the star-planet line,

The gas parameters ρ, u, T, Eth = p∕(γ − 1) are the mass density, velocity, temperature, and thermal energy of the upper atmosphere. The gas pressure is p = ρRTμ, with the mean molecular weight μ and the gas constant R. The distance from the planet’s center is r, and t is the time. For the specific heat ratio γ, we adopted 5/3 for monatomic gas. The gravitational force g = −∂Φ∂r was derived fromthe Roche potential Φ along the star-planet line: (4)

(e.g., García Muñoz 2007). Here, G is the gravitational constant, Mp the planet’s mass, M* the stellar mass, and a the star-planet separation. The right-hand side of Eq. (3) includes the net volume heating rate Q and thermal conduction, which are both described in Sect. 2.2.

The simulations assumed an atmosphere composed entirely of hydrogen, meaning we neglected helium and other heavier elements, which are minor compared to hydrogen. For HD 189733b, we also neglected molecular species present in the lower atmosphere, but the strong ionization of the upper atmosphere is expected to destroy them efficiently (cf., Guo & Ben-Jaffel 2016, and Appendix C). Therefore, we only considered hydrogen atoms (H) and protons (H+). The production of H+ is calculated as (5)

where αion is the photoionization rate, αrec the radiative recombination rate, and αcol the collisional ionization rate (all quantities in cgs units). The parameters nH, nH+, ne are the number densities of neutral hydrogen atoms, protons, and electrons, respectively. We assumed quasi-neutrality, nH+ = ne. The photoionization rate is given by (6)

where σion,λ is the photoionization cross-section, FXUV,λ the spectral XUV flux at the planet’s orbit outside of the atmosphere, τλ (r) the optical depth at distance r from the planet’s center, and Eλ the photon energy (e.g., Murray-Clay et al. 2009). We computed the average photoionization cross-sections per spectral bin from the fits of Verner et al. (1996) and adopted photon energies corresponding to the central wavelengths of the bins. The recombination and collisional ionization rates were taken from Glover & Jappsen (2007). For the former, we adopted the case B1 coefficient (7)

while the latter is given by

where θ = ln Te with the electron temperature Te in eV.

The system of equations is normalized as described in Erkaev et al. (2016) and solved using the MacCormack scheme (MacCormack 1969). We slightly modified the code in comparison with the version described in Erkaev et al. (2016). Mainly, we extended the scheme with the total variation diminishing (TVD) property. The code extensions are described in detail in Appendix A. Another modification is introduced here. Our initial modeling attempts of HD 189733b resulted in unphysical behavior close to the lower boundary of the computational domain. Specifically, the velocity was increasing towards r0, which violated the mass conservation. This happened mainly because the outflow of HD 189733b’s upper atmosphere remains subsonic up to the L1 point, and is therefore highly subsonic close to the lower boundary. Schemes like the MacCormack method are not well suited to describe such flows. As the kinetic energy is very small for highly subsonic flows, we modified the code to solve for the thermal energy instead of the total energy. We thereby discretized the third term on the right-hand side of Eq. (3) equally to the flux terms in both MacCormack steps.

The code evolves the system of equations until a steady-state solution is achieved. We assumed that this condition is fulfilled when the mass flux throughout the simulation domain is approximately constant, meaning mass conservation is fulfilled to 1%. Moreover, we compared the total heating and cooling rates to confirm that energy conservation is fulfilled as well (cf. Appendix A.2).

2.2 Heating and cooling processes

The main source of heating is the stellar XUV radiation. The XUV volume heating rate can be written as (13)

where the term 1 + ατλ(r) takes into account 2D effects of the energy absorption in an approximate way (Sekiya et al. 1980). The optical depth is calculated as (14)

We adopted a constant photoelectron heating efficiency ηph of 50%. This is intermediate between the values calculated for Jupiter (63%; Waite et al. 1983) and HD 209458b (20–40%; Shematovich et al. 2014), and was adopted by studies similar to ours (Shaikhislamov et al. 2014, 2016). The heating rate from Eq. (13) yields very similar results to the 2D method described in previous papers (Erkaev et al. 2013, 2016) for α = 4. We reverted to a 1D calculation of QXUV, because we introduced the usage of XUV spectra here, motivated by the previously demonstrated effects of the assumed SED (Guo & Ben-Jaffel 2016). This also allowed inclusion of X-ray heating, which is important for hot Jupiters. The calculation of a non-gray 2D heating function would be computationally more expensive, but the resulting QXUV is very similar to the 1D method shown in Eq. (13) (see Appendix B for a comparison). The adopted XUV spectra are described in Sect. 2.4.

We included several cooling processes, namely Lyα cooling (H excitation), radiative recombination, free–free emission (Bremsstrahlung), and collisional ionization. All cooling rates were taken from Glover & Jappsen (2007) and are given in cgs units. The Lyα cooling rate (15)

was multiplied by a factor of 0.1, as suggested by Koskinen et al. (2013) based on detailed calculations of the photon escape probability in the upper atmosphere of HD 209458b (Menager et al. 2013), to account for the optical thickness. The cooling rate by radiative recombination is (16)

by collisional ionization (17)

and by free–free emission (18)

with gff = 0.79464 + 0.1243log10(T). The net heating rate is then computed via (19)

Thermal conduction is also included in the model. The conductivity coefficient is calculated as κ = n−1jnjκj (e.g. García Muñoz 2007), where nj are the number densities of the constituents and κj = AjTsj, where Aj and sj are the fitting parameters for the individual conductivities. We adopted AH = 379, AH+ = 7.37 × 10−8, Ae = 1.2 × 10−6, sH = 0.69, and sH+ = se = 2.5 (García Muñoz 2007). We note that the contribution of H+ is negligible compared to that of the electrons.

We computed the total heating and cooling rates to confirm the energy conservation in our code (Appendix A.2). The total cooling rate consists of the four explicitly included cooling processes described above, in addition to adiabatic cooling, which is implicitly included in Eq. (3). Adiabatic cooling includes the contributions of advection (20)

and expansion (21)

(e.g. Salz et al. 2016a). The total heating rate is mostly dominated by XUV heating (Eq. (13)), but locally, advection can also contribute significantly, especially at larger heights. We note that conduction can also both cool and heat the gas, but we found it to be negligible, because its contribution to the energy balance is significantly smaller than the other processes onthis specific planet.

2.3 Initial and boundary conditions

The simulation domain extends from the planet’s optical transit radius r0 = Rp to r1 = 4.5 Rp. The upper boundary lies near the L1 point RL1 = (μμ2∕3)a ~ 4.3 Rp, where (e.g., Erkaev et al. 2007). We note that due to the high gravity of HD 189733b, the distance between Rp and typical mesopause pressure levels is very small. We used a nonuniform grid with typically N = 5000 grid points. This high number is necessary because of the steep density gradient and the small velocities near the inner boundary that are present on this planet.

At the lower, subsonic inflow boundary, we fixed the values for temperature T0 and number density n0. For the temperature, we adopted the equilibrium temperature of 1200 K (cf., Sect. 2.4) and for the number density a value of 1015 cm−3. This corresponds to a pressure of 0.16 mbar. Changing T0 to 800 or 1600 K decreases or increases, respectively, the mass-loss rate by only ~10%, and affects the atmospheric profiles negligibly (cf., Sect. 3.1). We checked that the assumed density n0 is chosen largeenough so that the optical depth in all XUV spectral bins exceeds 10. Only then is the stellar XUV flux completely absorbed in the simulation domain, and the resulting mass-loss rates do not strongly depend on the choice of n0. Smaller values of n0 would lead to significant ionization below the optical transit radius, which would be inconsistent with its measurement (see Sect. 6.3). For runs that resulted in very subsonic flows, we set u = 0 at the lower boundary. This does not affect the results, but improves the speed of convergence and the quality of the atmospheric profiles (e.g., Shaikhislamov et al. 2014). At the upper boundary, we adopted free outflow (i.e., zero gradient) conditions on all parameters.

As initial conditions, we used a constant temperature and a monotonically increasing velocity profile u(r) = 0.1(r − 1). For the number density, we assumed hydrostatic conditions n(r) = exp(Φ(r) − Φ(r0)), but modified it in the upper region to n(r) ∝ 1∕r2. This was found to be necessary because of the steep density gradient from the hydrostatic solution for this high-gravity planet. Initially, we assumed that the atmosphere consists only of atomic H, meaning nH+ (r) = 0. Practically, most runs were started from previous solutions with similar parameter sets to reduce the computing time. We considered a run as properly converged when the mass flux ρur2 is spatially constant within 1% throughout the computational domain (excluding a small region close to the lower boundary where the mass flux may drop to zero for numerical reasons; cf., Fig. A.1). We checked that the adopted initial conditions and the chosen resolution of the grid do not affect the results. Moreover, we verified that the chosen location of the upper boundary does not affect the results.

Table 1

Stellar and planetary parameters from Stassun et al. (2017). The orbital distance a was calculated from their quoted value of aR* and R*.

2.4 Adopted physical parameters

The adopted planetary and stellar parameters (Table 1) were taken from Stassun et al. (2017) and are based on the accurate parallaxes from the first data release of Gaia (Gaia Collaboration 2016). The distance d, stellar mass M* and radius R*, stellar bolometric flux at Earth Fbol, planetary mass Mp and radius Rp, as well as the orbital inclination i were taken from their study. The orbital separation a was calculated from their quoted valueof the parameter aR* = 8.84 ± 0.27 and the stellar radius. The average orbital velocity is therefore 151 km s−1, consistent with the maximum orbital radial velocity of km s−1 measured by de Kok et al. (2013). Using the measured bolometric stellar flux, we find an equilibrium temperature of 1200 K, assuming full redistribution (f = 4) and zero albedo (A = 0). Here, σ is the Stefan-Boltzmann constant. This is similar to, but slightly higher than, the apparent effective dayside temperature of 1163 ± 37 K (Schwartz et al. 2017). Such small differences in the adopted Teq do not affect our results.

The XUV flux of HD 189733 is strongly variable due to the high activity of this star. Observed X-ray luminosities are in a range of 1.1−2.8 × 1028 erg s−1, which were partly obtained from different instruments with slightly different bandpasses, but were also measured at different epochs (Hünsch et al. 1999; Pillitteri et al. 2010, 2011, 2014; Lecavelier des Etangs et al. 2012; Poppenhaeger et al. 2013). Since the extreme ultraviolet (EUV) part of the stellar spectrum is largely unobservable due to absorption by the interstellar medium (ISM), it has to be inferred indirectly. We compared two approaches to estimate the unobservable EUV spectrum. Firstly, we used the scaling relations from Linsky et al. (2014, hereafter L14) based on the intrinsic Lyα flux, as derived from observations of the short-wavelength part of the EUV range (<400 Å) in nearby stars, and solar models for >400 Å. Since Lyα fluxes are also affected by the ISM absorption, the intrinsic stellar line profile needs to be reconstructed. From HST observations of HD 189733 in 2010, the intrinsic Lyα flux at Earth amounts to 7.5 × 10−13 erg cm−2 s−1 with quoted uncertainties of 15–30% (France et al. 2013). The reconstructed line profile from Bourrier et al. (2013) for the 2011 observations corresponds to a flux of 6.8 × 10−13 erg cm−2 s−1, slightly lower, but consistent with France et al. (2013), considering the typical uncertainties of the reconstruction. We adopted the latter value, because we are focusing on modeling the 2011 observations. This flux value, together with the scalings from L14, gives a total EUV (100–912 Å) flux at the planet’s orbit of 6.5 × 103 erg cm−2 s−1 with a spectral energy distribution shown in Fig. 1.

Secondly, we examined the synthetic spectrum from the X-exoplanets2 database (Sanz-Forcada et al. 2010, 2011, hereafter SF11). These authors used observed X-ray and UV spectra, reconstructed the emissionmeasure distribution, and used it as an input for coronal models to infer the unknown EUV part of the spectrum. We show their synthetic spectrum with the same binning as for L14 in Fig. 1. The errorbars in x-direction give the width of the bins (100 Å), those along the y-axis the estimated flux errors based on the uncertainties quoted in the respective studies. It is apparent that the SF11 fluxes are systematically higher than the L14 ones (except for the 800–912 Å bin, which is likely due to the peak of the Lyman continuum, which is included in the solar models used by L14, but not in SF11). The λ > 300Å range is consistent within the estimated uncertainties. The λ < 300Å range differs significantly. This could be due to shortcomings in the models of either study, and/or the fact that the X-ray spectrum used by SF11 was taken at a different epoch (2007) than the Lyα flux we adopted for the L14 method (2011). Some intrinsic stellar variability could therefore also be a cause of the differences. Pillitteri et al. (2014) gave a comparison of the X-ray fluxes obtained with XMM-Newton at different epochs (2007 = SF11, 2009, 2011, 2012) and the associated temperatures and emission measures. In 2011, the temperatures were comparable, but the emission measures and X-ray fluxes were higher than in 2007. However, the X-ray observations in 2011 were not taken at the same time as the Lyα observations. Simultaneous observations with Swift/XRT during the 2011 transit yielded an X-ray (0.3–3 keV) flux of 3.6 × 10−13 erg cm−2 s−1 (Lecavelier des Etangs et al. 2012), whereas the XMM-Newton (0.12-2.48 keV) flux in 2007 was 3.4 × 10−13 erg cm−2 s−1 (SF11), and 3.2−3.9 × 10−13 erg cm−2 s−1 in 2011 (excluding the flare; Pillitteri et al. 2014). Hereafter, we adopted the average XMM-Newton value (excluding the flare) from 2011, 3.55× 10−13 erg cm−2 s−1 (Pillitteri et al. 2014), which amounts to 6.2 × 103 erg cm−2 s−1 at the orbit of HD 189733b.

The total EUV fluxes (100–912 Å) at the planet’s orbit obtained with the two different methods are 6.5 × 103 erg cm−2 s−1 for the L14 method and 1.2 × 104 erg cm−2 s−1 for SF11. As a comparison, we also calculated the total EUV flux using the method of Chadney et al. (2015), which uses a scaling with X-ray flux. With our adopted X-ray flux, we obtain an orbital EUV flux of 1.1 × 104 erg cm−2 s−1, comparable to SF11. We note, however, that the Chadney et al. (2015) scaling is based on the SF11 coronal models. On the basis of the discrepancies between the different reconstruction methods, we estimate that the uncertainty in EUV flux is at least a factor of two. The total XUV (5–912 Å) fluxes at the planet’s orbit from each method amount to 1.3 × 104 and 1.8 × 104 erg cm−2 s−1 for L14 and SF11, respectively (cf., Table 2). We explore the influence of XUV flux variations further in Sect. 3.2, where we model the effects of a flare.

thumbnail Fig. 1

XUV (5–912 Å) spectral flux at the orbit of HD 189733b in bins of 100 Å, obtained with the methods of L14 and SF11. The X-ray flux is an average value outside of flares obtained with XMM-Newton in 2011 (Pillitteri et al. 2014).

Open with DEXTER
Table 2

Modeled mass-loss rates for different stellar XUV fluxes/spectra and lower boundary parameters (number density n0, temperature T0).

3 Hydrodynamic modeling results

3.1 Atmospheric profiles and mass-loss rate

Here, we show the upper atmosphere profiles obtained using both the SF11 and L14 XUV spectra (Fig. 2). The isotropic (i.e., maximum) mass-loss rates = 4πr2ρu (assuming that the star–planet line value is representative for entire atmosphere) are 5.4 × 1010 and 2.5 × 1010 g s−1 for SF11 and L14, respectively (Table 2).

In Fig. 2, we show the atmospheric profiles for number density, velocity, temperature, and ionization fraction (nH+n) for both the SF11 and L14 reconstructions of the stellar XUV spectra. The results are rather similar for both spectra, although all parameters are slightly lower for the L14 spectra because of the lower total XUV flux. The number density is dominated by neutral H below about 1.85 Rp and by H+ above. The transition from a H-dominated atmosphere to one dominated by H+ occurs where the ionization fraction exceeds 0.5. The outflow reaches the sonic speed just slightly below the L1 point at 4.3 Rp in both cases, whereas it already exceeds the escape velocity at points that are ~ 0.3 Rp closer to the planet. The temperature reaches a maximum of ~1.1 × 104 K at ~2 Rp for SF11, whereas the temperature maximum is slightly lower for L14 and located closer to the planet.

The resulting atmospheric profiles are only considered to be valid well within the Roche lobe, and less reliable close to or even above the Roche lobe where 3D effects become significant (e.g., Bisikalo et al. 2013). Moreover, it is important to check if the outflow remains collisional within our computational domain so that the hydrodynamic treatment can be justified. The transition level to the collisionless regime is commonly taken as Knudsen number Kn = Λ∕X = 1, where Λ =1∕(col) is the mean free path and X is an appropriate system scale, for example,the scale height H = kT∕(mpg) for regions close to the planet, or the planetary radius for more distant regions (Shaikhislamov et al. 2014). For the collision cross-section σcol, we adopted the H–H+ charge exchange cross-section (~ 2 × 10−15 cm2; e.g., Lindsay & Stebbings 2005), which is the appropriate choice for partially ionized hydrogen atmospheres (Guo 2011; Salz et al. 2016a; Shaikhislamov et al. 2016). Taking for X either H or Rp, we find that Kn ≪ 1 inside our computational domain, justifying the hydrodynamic approach.

Figure 3 details the individual heating and cooling processes in our model for the SF11 run. One can see that the main heating source is the stellar XUV radiation, although advection contributes to heating in the uppermost regions. Advection cools the gas at lower heights, although its contribution is much smaller than the cooling from expansion. Expansion is even the dominant cooling mechanism above about 2.5 Rp. Radiative cooling is also very important for HD 189733b and is dominated by Lyα emission close to the planet, whereas, in the upper layers, recombination radiation and free-free emission are more important. Both collisional ionization and conduction (not shown) are negligible cooling mechanisms.

We calculated the heating efficiency using two different definitions commonly found in the literature. Firstly, we calculated the heating efficiency ηXUV(r) defined as the ratio of the XUV volume heating rate QXUV (Eq. (13)) to the locally absorbed XUV radiation (22)

Secondly, wecalculated ηnet using the net localheating rate Qnet(r) (Eq. (19)) divided by Qabs(r) (e.g., Salz et al. 2016a). A comparison is shown in Fig. 4. The corresponding mean heating efficiencies (, ) in the atmosphere, obtained by integrating QXUV(r) and Qnet(r), respectively,over r and dividing by the stellar XUV flux at the planet’s orbit (Table 2) amount to 12 and 3%. The former is in good agreement with detailed studies of the hot Jupiter HD 209458b (Shematovich et al. 2014; Ionov & Shematovich 2015). We note that these heating efficiencies are different quantities than the photoelectron heating efficiency ηph described in Sect. 2.2.

The effective XUV absorption radius (23)

(Erkaev et al. 2007, 2015) is about 2.6 Rp for SF11 and 2.1 Rp for L14. Using Qnet instead of QXUV in Eq. (23) yields 1.15 (SF11) and 0.78 Rp (L14), respectively. We note that this radius can be smaller than 1 if Qnet is used in Eq. (23) and radiative cooling is very efficient, and/or a too large value is chosen for the mean heating efficiency in the denominator. Evaluating the energy-limited mass-loss rate (24)

(Erkaev et al. 2007) with these parameters yields 1.5−3.1 × 1011 g s−1 if adopting a typical mean heating efficiency for hot Jupiters of % (e.g., Shematovich et al. 2014) and calculating the tidal enhancement factor K as described in Erkaev et al. (2007). This overestimates the mass-loss rate of HD 189733b by almost an order of magnitude. Adopting the RXUV values obtained with Qnet yields 2.1−6.2 × 1010 g s−1, in much better agreement with the hydrodynamic rates. We note that the simple assumption of RXUV~Rp would coincidentally yield, for this specific planet, mass-loss rates of 3.4−4.7 × 1010 g s−1, which are similar to the hydrodynamic results.

thumbnail Fig. 2

Planetary atmospheric profiles adopting the XUV spectra of SF11 (solid) and L14 (dashed). Shown are the number densities (total, atomic H, and protons), outflow velocity, ionization fraction, and temperature. In the velocity plot, crosses indicate the sonic points and plus signs denote the locations where the outflow velocity exceeds the escape velocity. The vertical dotted line indicates the location of the L1 point.

Open with DEXTER
thumbnail Fig. 3

All heating (red) and cooling (blue) processes included in our model, shown for the SF11 run. The cyan lines detail the individual contributions to the radiative cooling rate (solid blue line).

Open with DEXTER
thumbnail Fig. 4

Heating efficiencies ηXUV(r) (solid) and ηnet(r) (dashed).

Open with DEXTER

3.2 The effect of a flare

The star HD 189733 is magnetically active, and flares have been detected frequently (Pillitteri et al. 2010, 2011, 2014). An excess absorption of ~14% in Lyα was reported to have occurred during a planetary transit about eight hours after a strong X-ray flare was detected (Lecavelier des Etangs et al. 2012). Here, we study if the enhanced XUV flux emitted by such a flare could have increased the planetary mass-loss rate enough to have caused this absorption.

The X-ray flux in the 0.3–3 keV band of Swift/XRT was increased by an average factor of four for a duration of ~27 min, which corresponds to one bin of the temporal resolution of the observations. Together with the average pre-flare flux of 3.6 × 10−13 erg cm−2 s−1, this yields an estimated energy of ~8 × 1031 erg in this bandpass, and is thus a lower limit to the total radiated energy of this flare.

We tested the effect of such a flare on the planetary mass-loss rate and the atmospheric profiles. Two cases were considered: firstly, we increased the total XUV flux by the factor of four found in X-rays; secondly, we assumed that only the X-ray flux increased by a factor of four and the EUV flux remained unchanged, since the hot flare plasma could have radiated predominantly at shorter wavelengths, changing the shape of the XUV SED (like in solar flares which radiate mostly in X-rays; Chadney et al. 2017). We did not take into account the duration of the flux enhancement, but instead calculated the steady-state solution with these enhanced fluxes. If this solution does not show a significant increase of the neutral mass-loss rate, a shorter pulse of enhanced radiation would not either. The same conclusions were reached by Chadney et al. (2017), who investigated both the time evolution and steady-state enhancement effects of flares. We note that the true peak flux of the flare is unknown due to the low time resolution of the observations. The duration of the observed flare was approximately 27 min, since there was a significant enhancement in only one temporal bin, and the flux in the following bin was comparable to (even slightly lower than) pre-flare levels (see Fig. 4 in Lecavelier des Etangs et al. 2012). Thus, even though the true peak flux could have been higher, its duration cannot have exceeded more than a few minutes, limiting its effect on the planet.

Figure 5 shows the atmospheric profiles under flaring conditions compared to those for the average XUV (SF11) spectrum. The profiles for all cases are rather similar. The most apparent differences are the stronger ionization for the XUV flare case and the slightly lower ionization for the X-ray flare case. Moreover, the temperatures in the uppermost regions are slightly higher for the flare cases than for the average XUV flux, but the peak is only slightly increased for the XUV flare. Other than that, the profiles are not strongly altered by exposure to the flare radiation, similar to what was found in the study of Chadney et al. (2017).

The total mass-loss rates are 8.7 × 1010 (X-ray flare) and 1.2 × 1011 g s−1 (XUV flare), higher by factors of 1.6 and 2.2, respectively, compared to the non-flaring state. Our flare-related mass-loss rate enhancements are very similar to the range of 1.8–2.2 found by Chadney et al. (2017). To compare with the modeling of Bourrier & Lecavelier des Etangs (2013), the neutral mass-loss rates at 2.95 Rp (the particle launch radius in their model) are 2.2 × 109 (X-ray flare) and 8.9 × 108 g s−1 (XUV flare), meaning a factor of two higher (X-ray flare) and about 10% lower (XUV flare) than in the non-flaring state (109 g s−1). This indicates that the enhanced radiation from the flare does not increase the neutral H densities and outflow fluxes sufficiently to account for the observed drastically different atmospheric absorption. More interestingly, although the total mass-loss rate increases for the XUV flare compared to both the non-flaring state and the X-ray flare because of the larger energy input, the neutral loss rate is actually smaller due to increased ionization. We note that the neutral outflow rates we obtain are in agreement with those required for the detected absorption (5 × 108−1.5 × 109 g s−1) in the model of Bourrier & Lecavelier des Etangs (2013), but for both flaring and non-flaring states. This indicates that the occurrence of a flare is likely insufficient to produce such strong variability in the planetary atmospheric absorption, especially if recalling that the actual flare duration was neglected. The nondetection in 2010 would require a reduction of the neutral escape rate by a factor of 5–20 compared to 2011 (Bourrier & Lecavelier des Etangs 2013). This could only happen if the star were much less active in 2010 compared to our adopted XUV fluxes (reducing the total escape rate), or if it were much more active to increase the ionization sufficiently to reduce the neutral escape rate, despite increasing the total one. From existing observations, the variability in X-rays is about a factor of three (cf., Sect. 2.4), which is likely not sufficient to account for the variability of atmospheric absorption. However, this contradicts the necessity of having similar ionizing fluxes at both epochs in the Bourrier & Lecavelier des Etangs (2013) model. Therefore, other processes may be responsible for, or contribute to, the variable planetary absorption, like variations of the stellar wind. Since strong X-ray flares on the Sun are frequently accompanied by coronal mass ejections (CMEs), we also consider the effect of a possible CME impact in Sect. 4.3. We note that other upper atmosphere models of this planet find higher ionization (and therefore lower H densities) at 2.95 Rp, yielding much lower neutral loss rates (Guo 2011; Salz et al. 2016a; Chadney et al. 2017). This is mainly related to the smaller number densities assumed at the lower boundary in these studies.

4 MHD flow modeling

4.1 3D MHD flow model

We used the 3D MHD flow model described in Erkaev et al. (2017) to calculate the plasma flow around HD 189733b (assuming no intrinsic planetary magnetic field) and to obtain the plasma parameters and magnetic field in the region around the planet. This model takes into account radiation and charge-exchange processes acting on the hydrodynamically expanding upper planetary atmosphere penetrating into the stellar wind plasma.

Magnetic field and plasma parameters are determined by the following system of equations for mass, momentum and energy conservation, which are completed by the magnetic induction equation

where ρ, V, P, and B are the mass density, velocity, plasma pressure, and magnetic field of the stellar wind, respectively. The parameter γ is the polytropic index (assumed to be equal to 5/3) and I is the identity matrix. The parameters VH and TH are the velocity and temperature of the atmospheric neutral hydrogen atoms escaping from the planet.

The mass conservation equation for ions includes an interaction source term, which is related to photoionization (31)

and charge exchange ionization (32)

of the hydrogen atoms. Here, NH is the number density of the neutral planetary hydrogen atoms, mp the particle mass, σex (~10−15 cm2) the charge exchange cross-section, ⟨Vrel⟩ the average relative speed of the stellar wind and atmospheric particles, and αi = 5.9 × 10−8FXUV is the ionization rate proportional to the XUV flux FXUV = 1.8 × 104 erg cm−2 s−1 (cf., Table 2). The numerical scheme for the solution of this system of equations is described in Erkaev et al. (2017).

The neutral hydrogen atoms can be ionized via photoionization or charge exchange processes. However, the latter is a dominating ionization mechanism for the atmospheric atoms in the stellar wind region, which leads to the generation of low-energy planetary ions and high-energy neutral hydrogen atoms (ENAs) originating from the stellar wind. These newly born ions are immediately accelerated by the local electric field and start to move together with the stellar wind plasma around the planetary obstacle. An acceleration of the picked up ions is accompanied by deceleration of the stellar wind plasma, with conservation of the total momentum.

The streamlined obstacle is considered to be a semi-sphere. The position of the stellar wind stagnation point (Rs) is determined by the pressure balance condition, which means that the total pressure of the external stellar wind has to be equal to the momentum flux of the internal ionized atmospheric particles at the boundary. The ratio of the curvature radius of the obstacle to the distance between the stagnation point and the planetary center was taken as 1.3, similar to the value used by Erkaev et al. (2017).

The calculation domain for the MHD stellar wind flow is bounded by the external semi-sphere related to the undisturbedstellar wind region and the internal semi-sphere corresponding to the planetary obstacle. At the outer boundary, we applied the undisturbed stellar wind parameter values of density, velocity, temperature, and magnetic field. At the obstacle boundary, we assumed that the normal components of the stellar wind velocity and magnetic field vanish. Finally, we obtained a stationary solution for the stellar wind flow as a result of time relaxation of the nonsteady MHD solution. As initial conditions, we set the undisturbed stellar wind parameters in the whole computational domain. The final stationary solution is unique and does not depend on the particular initial conditions.

thumbnail Fig. 5

Atmospheric profiles for average XUV (solid), X-ray flare (dashed), and XUV flare (dash-dotted) cases based on the SF11 spectrum.

Open with DEXTER
Table 3

Adopted stellar wind parameters at the planet’s orbit.

4.2 Stellar wind interaction

Table 3 summarizes the adopted stellar wind parameters used for the flow model. The mean and maximum values were taken from a 3D stellar wind model (Llama et al. 2013) based on the measured magnetic field map of the star (Fares et al. 2010). Since all wind parameters vary along the orbit depending on longitude, we considered two cases. For the first case (mean wind), we computed the mean values of number density Nsw, velocity Vsw, temperature Tsw, and magnetic field strength along the orbit. We computed both the parallel Bsw,p (parallel to Vsw) and normal Bsw,n (normal to Vsw) components of its magnitude, as well as the total field strength Bsw,tot. The velocity is dominated by its radial component, and the magnetic field at the orbit is also dominated largely by its radial (i.e., parallel) field component. The angle θB is the angle of the total magnetic field vector and the star-planet line. For the second case (maximum wind), we used the parameter values where the ram pressure of the wind along the orbit is maximal. This corresponds to higher values compared to the mean wind (except for B), but not to the maximum values of the individual parameters. We chose this approach because the parameter maxima usually lie at different longitudes (e.g., the maximum velocity corresponds to the minimum density; cf., Fig. D.1). We ignored the velocity of the planetary orbital motion in the following calculations, as the ram pressure of the orbital motion is a factor of five smaller than that of the radial wind for the mean wind conditions, and even correspondingly smaller for the other scenarios.

We adopted the hydrodynamic solution of the planetary upper atmosphere obtained with the SF11 XUV spectrum as an input for the MHD flow modeling. Considering the stellar wind parameters given in Table 3, we applied the 3D MHD flow model to calculate the spatial distribution of the magnetic field and plasma parameters in the vicinity of the planet. By solving the nonsteady MHD equations, a stationary solution is established as a result of time relaxation. As an initial condition, we assumed the uniform undisturbed stellar wind flow, which suddenly stops at the planetary obstacle. This leads to the appearance of a shock-like wave front propagating outwards from the obstacle. Since we have a super-Alfvenic stellar wind flow, this shock approaches its stationary position located at some distance from the obstacle.

In the case of the mean wind, the radial distance to the magnetopause is about 2.6 Rp. Here, the ion and neutral densities of the atmospheric particles are 3.9 × 107 cm−3 and 1.25× 106 cm−3, respectively, and the temperature is 8.8 × 103 K. In the maximum wind case, the magnetopause position is slightly closer, at 2.4 Rp, since the stellar wind dynamic pressure is larger. The corresponding atmospheric ion and neutral densities are 6.1 × 107 cm−3 and 2.9 × 106 cm−3, respectively, and the temperature is 1 × 104 K (cf. Fig. 2).

Figure 6 shows the spatial distribution of the total pressure Ptot (sum of the magnetic and gas pressures) obtained from the MHD model employing the mean stellar wind parameters (see Table 3). Here, the total pressure is normalized to the stellar wind dynamic pressure Pdsw = Pa. The origin of the coordinate system is placed at the planet’s center, and the star is located along the positive X-axis. The direction of the Z-axis is chosen to have coplanarity between the XZ plane and the interplanetary magnetic field vector Bsw (the arrow in Fig. 6). The white area around the center of the coordinate system indicates the region filled by atmosphericparticles exclusively, while the embedded dark blue semicircle indicates the optical planetary radius Rp. Figure 7 is similar to Fig. 6, but corresponds to the maximum stellar wind parameters. The total pressure in this case is also normalized to the corresponding unperturbed stellar wind dynamic pressure of Pdsw = 1.8 × 10−4 Pa.

Figure 8 shows the profiles of the ENA and ion densities along the X-axis between the magnetopause and the bow shock (top panel), as well as the profiles of the total pressure and temperature (bottom panel), bothcorresponding to the mean stellar wind. The total ion number density has a maximum of about 4 × 104 cm−3 due to charge exchange interaction between the stellar wind plasma and atmospheric neutral atoms, which are flowing through the magnetopause. In front of the magnetopause, the ion density and ENA density maxima are about 1.0 × 107 cm−3 and 2.3 × 107 cm−3, respectively.The ENAs form a layer around the magnetopause with a thickness of about 1.3 Rp. Figure 9 is similar to Fig. 8, but corresponds to the maximum stellar wind parameters. In this case, the ion density and ENA density maxima in front of the magnetopause are about 1.13× 107 cm−3 and 5.6 × 107 cm−3, respectively. The ENA layer thickness around the magnetopause is about 1.2 Rp.

thumbnail Fig. 6

Cut at Y = 0 of the simulation showing the distribution of the total pressure (sum of the magnetic and gas pressures) around HD 189733b normalized to the stellar wind dynamic pressure for the mean stellar wind parameters. The white area close to the origin indicates the atmospheric region around the planet, while the semicircle indicates the planetary optical radius. The star is located along the X-axis. The arrow shows the direction of the interplanetary magnetic field.

Open with DEXTER
thumbnail Fig. 7

Same as Fig. 6, but for the maximum stellar wind.

Open with DEXTER
thumbnail Fig. 8

Profiles of the parameters along the stagnation line from the magnetopause to the bow shock for the mean stellar wind conditions. Top panel: ENA and ion densities. Bottom panel: total pressure and temperature.

Open with DEXTER
thumbnail Fig. 9

Same as Fig. 8, but for the maximum stellar wind conditions.

Open with DEXTER

4.3 Impact of a coronal mass ejection

We also consider the case of a stellar CME to study the effect of a potential CME impact. Interaction with a CME could have provided the elevated plasma density levels necessary to produce the observed variability of the planetary Lyα transit signature. On the Sun, large flares are frequently accompanied by CMEs. Empirical models based on flare-CME relationships from the Sun predict high CME occurrence rates for active stars (e.g., Odert et al. 2017, and references therein), because of their high flare rates. The estimated X-ray flare energy of ~ 8 × 1031 erg (cf., Sect. 3.2) corresponds to a flare that is (almost) always accompanied by a CME on the Sun (Yashiro & Gopalswamy 2009). If the flare was indeed accompanied by a CME, its estimated mass would be in the order of 1016 g, and its velocity about 1300 km s−1, based on solar scalings (Drake et al. 2013; Odert et al. 2017). We applied an empirical CME prediction model to HD 189733, which calculates CME occurrence rates for Sun-like and cooler main-sequence stars based on their X-ray luminosities LX (Odert et al. 2017). Adopting LX = 1.67 × 1028 erg s−1 (cf., Sect. 2.4), we obtain about 10−2000 CMEs per day, depending on the power law index of the stellar flare energy distribution, d N∕dEEα. We assumed α = 1.5−2.5, which corresponds to an observationally determined range typical for the Sun and other stars (Güdel et al. 2003). The mass-loss rate from these CMEs are consistent with observations of stellar mass-loss rates (e.g., Wood 2004, and references therein). However, recent modeling results suggest that CME rates may be lower than estimated from solar scalings due to the stronger magnetic fields on active stars (Alvarado-Gómez et al. 2018), so the obtained numbers are possibly overestimates.

Knowledge of stellar CMEs is still sparse, and therefore the CME activity of HD 189733 is not constrained by observations either. One method to detect stellar CMEs is from transient blue-wing asymmetries or blue-shifted extra-emissions or absorptions (depending on viewing geometry) in Balmer lines, which probe the (partly) neutral prominence material embedded in the CME core (e.g., Leitzinger et al. 2014, and references therein). We analyzed archival optical spectroscopic observations of HD 189733 (430 Stokes I spectra from PolarBase; Petit et al. 2014) as described in Leitzinger et al. (2020). The total on-source time of the spectra is about 122 h, and the observations span the years 2006–2015. We find no signatures of CMEs in the Balmer lines. With this nondetection, we estimate <0.6 observable CMEs per day (95% confidence). It is important to correct for the expected Hα emission and geometrical constraints, since not all occurring CMEs are necessarily observable with a given method. We used a semiempirical model that takes into account the expected maximum possible CME core emission in the Balmer lines and geometrical constraints, considering the total on-source time and average signal-to-noise ratio (S/N) of the spectra (Odert et al. 2020). It predicts that ≲1% of the intrinsically occurring CMEs would be observable in Hα on HD 189733 for the given observational parameters. The estimated maximum observable CME rates are comparable to the upper limit from the analyzed observations, indicating that much longer observing times and/or higher quality of the spectra would be required to determine the CME rate of this star with the given method. We note that none of the analyzed Hα observationswere taken during the time of the studied flare and planetary transit in 2011. Due to the nondetection of CMEs on this star, we have to rely on reasonable extrapolations from the Sun to estimate plausible CME parameters for HD 189733.

We assumed a velocity of 1000 km s−1, which is commonly found for energetic CMEs already close to the Sun (Yashiro et al. 2004), and which is also similar to the CME velocity derived above from the estimated flare energy. For the density, we used an enhancement factor of ten compared to the mean wind, which is also common for solar CMEs at separations of a few solar radii (Schwenn et al. 2006). We note that we cannot use solar CME density profiles as adopted, for instance, in Khodachenko et al. (2007), because the modeled stellar wind densities at HD 189733b’s orbit are already higher than the solar CME densities at the same distance (Khodachenko et al. 2007). Due to the lack of knowledge regarding stellar CME temperatures, we simply adopted the maximum wind value. For the magnetic field, we used 100 mG, similar to values of solar CMEs at around 10 R from the Sun (Patsourakos & Georgoulis 2016). The angle θB is assumed to be 90°, corresponding to the geometry of a centrally impacting ejected flux rope, as commonly observed in solar CMEs. With the assumed velocity, a CME would need about 1 h to reach the planet’s orbit. The typical duration of solar CMEs at a few solar radii is about 8 h (Lara et al. 2004). This means that if a CME occurred simultaneously with the flare, the planet may have still been exposed to a plasma environment dominated by the CME. The adopted CME parameters are summarized in Table 3.

Figure 10 shows the spatial distribution of the total pressure obtained with the MHD flow model for the adopted CME parameters (Table 3). In this case, the magnetopause is located closer to the planet compared to both wind cases because of the higher dynamic pressure, namely at a distance of ~ 2 Rp. The corresponding atmospheric ion and neutral densities are about 1.8 × 108 cm−3 and 4.5 × 107 cm−3, respectively, and the temperature is 1.16 × 104 K (cf. Fig. 2). The total pressure is again normalized to the CME dynamic pressure, Pdsw = 7.35 × 10−3 Pa. The ion density and ENA density maxima in front of the magnetopause (Fig. 11) are about 0.67× 107 cm−3 and 1.5 × 108 cm−3, respectively.The ENA layer is much thinner in the CME case, with a thickness of about 0.25 Rp.

Table 4 summarizes the results of the stellar wind interaction modeling. It compares the stand-off distances, the thickness of the ENA layers, and the maximum number densities of ENAs and ions in front of the magnetopause. Stronger winds or CMEs confine the planetary atmosphere to smaller regions because of the higher ram pressures. Higher stellar wind ram pressures also lead to more compressed, thinner ENA layers, but with higher maximum densities.

thumbnail Fig. 10

Same as Fig. 6, but for the CME conditions.

Open with DEXTER

5 Modeling the Lyα absorption signature

To calculate the transmissivity in the Lyα line, we assumed that a hydrogen cloud surrounding the planet consists of two parts: a spherically symmetric lower atmosphere corresponding to the 1D atmospheric profile, and an upper exospheric part consisting of the ENAs calculated by the 3D MHD model. We considered one atmospheric profile (Fig. 2) and the three different stellar wind cases (mean, maximum, and CME; Table 3). After a hydrogen cloud was simulated (Figs. 6, 7, and 10) and through knowledge of the positions and velocities of all hydrogen particles, we computed how these atoms attenuate the stellar Lyα radiation byusing a post-processing software written in the Python programming language. To compute the transmissivity along the line of sight, we followed the approach of Semelin et al. (2007). The post-processing tool is described in detail in Kislyakova et al. (2014). Here, we summarize its main features.

Only neutral hydrogen atoms absorb in the Lyα line. One has to take into account spectral line broadening. Real spectral lines are subject to several broadening mechanisms: (i) natural broadening; (ii) collisional broadening; (iii) Doppler or thermal broadening. The “natural line width” is a result of quantum effects and arises due to the finite lifetime of an atom in a definite energy state. A photon emitted in a transition from this level to the ground state will have a range of possible frequencies Δf ~ ΔE ~ 1∕Δt, which can be approximated by a Lorentzian profile. Collisional broadening is caused by the collisions randomizing the phase of the emitted radiation. This effect can become very significant in a dense environment, yet above the exobase, it does not play a role and is important only in the lower parts of the atmosphere, therefore, we did not take it into account.

The third type of broadening, which plays a significant role in the upper atmosphere of a hot exoplanet, is thermal broadening, which arises because the frequency of the absorption is shifted due to the Doppler effect. In the considered cases, an analytical solution for the absorption profile cannot be obtained, since it is not only thermal atoms that contribute to the broadening, but also ENAs. For this reason, we cannot use a Voigt profile (which can only be used for a pure Maxwellian distribution). We calculated the natural broadening for all atoms and bin it by velocity, which automatically gives us the Doppler broadening for a particular velocity distribution.

To compute the transmissivity along the line of sight, we followed the approach of Semelin et al. (2007) and Kislyakova et al. (2014). We used a Cartesian coordinate system with the x-axis pointing towards the star, the y-axis directed antiparallel to the planetary motion, and the z-axis completing the right-handed coordinate system. We calculated the relation between the observed intensity I and the source intensity I0 as a function of frequency f of the stellar spectrum in the yz-plane by dividing the computational domain into a grid with Nc cells. For each cell in the grid along lines of sight in front of the star (), the velocity spectrum of all hydrogen atoms in the column along the x-axis can be calculated. We accounted for the planetary inclination by shifting the cloud by z = acosi = 3.67 × 1010 cm relative to the center of the stellar disk at mid transit. Then, the transmissivity could be averaged over all columns in the yz-grid except for those particles that fell outside the projected limb of the star or inside the planetary disk.

We used the frequency-dependent cross-section, which depends on the normalized velocity spectrum, the Lyα resonance wavelength and the natural absorption cross-section in the rest frame of the scattered hydrogen atom (Peebles 1993). For lines of sight in front of the planet , where (yp, zp) is the planet center position, we set the transmissivity to zero.

To account for the contribution of the lower atmosphere, a Maxwellian velocity spectrum corresponding to a hydrogen gas with a specified column density and temperature is added to all pixels up to the outer extent of the atmosphere, according to the atmosphericprofile calculated with the 1D hydrodynamic code. Similarly, to account for the ENAs, we added a Maxwellian spectrum corresponding to their temperatures with the central velocity located at their bulk velocity, which differs depending on their position relative to the magnetopause boundary. Due to the formation of a bow shock, ENAs are generated from decelerated stellar wind ions and are thus much slower than the stellar wind, but at the same time obtain a very high temperature, manifesting in a very broad spectrum spanning from very low up to very high velocities along the x-axis.

Figure 12 shows the velocity spectra of all neutral H atoms along the line of sight, including atmospheric particles (the main central peak in the plot) and the ENA population (seen as the wide “wings” on both sides of the main peak) for allthree considered stellar wind cases. Positive velocities and negative velocities correspond to atoms flying toward and away from the star, respectively. The velocity spectrum of the atomic cloud can then be converted to frequencies via the relation f =f0 + vxλ0 with f0 = cλ0, λ0 = 1215.65Å, and where c is the speed of light. As one can see, ENAs form very wide wings around the central atmospheric peak. Different initial values for temperature, density, and velocity of the stellar wind influence the height and the flatness of the wings. However, one can see that all three cases produce very flat wings, which is due to a sharp temperature increase and wind deceleration near the planetary boundary.

Figure 13 shows the calculated average transmissivity and absorption for the mean wind spectrum shown in Fig. 12. Despite differences in the velocity spectra in the ENA part, the transmissivity spectra are indistinguishable, so we do not show the other cases. This is due to the fact that the contribution from the ENAs produces very flat and low spectra, with the amount of particles in a given velocity bin not high enough to produce any significant contribution. Therefore, the absorption is mostly determined by the atmospheric broadening from the dense central peak, which is identical in all three cases, because we used the same atmospheric profile.

Finally, Fig. 14 compares the calculated absorption to the Lyα in-transit observations of HD 189733b from 2011. It was obtained by multiplying the out-of-transit spectrum (blue) by the transmissivity (Fig. 13). This simple method may generate a bias if applied to the analyzed observations because of the combination of the broad line-spread function of STIS and the ISM absorption, yet this is limited to the region close to the line core where the ISM absorption is saturated, but should not affect the line wings where the planetary absorption is detected. In general, there is a good agreement between the calculated and observed spectra, although the absorption predicted by our model is slightly larger. Outside of the geocoronal emission region, the modeled blue-wing (− 400… − 40 km s−1) absorption is 10.9%, 1.7σ larger than the observed value (6.8 ± 2.4%); for the red wing (+ 40… + 400 km s−1), we obtain, with 12.9%, too much absorption by 5.5σ (observed: 4.5 ± 1.5%). This means that the model with our default parameters predicts neutral densities that are slightly too large. If we scaled the modeled absorption profile to match the blue-wing absorption, the red-wing absorption would be within 3σ of the observations. The overabsorption could be due to limitations of our 1D hydrodynamic atmosphere model, as we used the atmosphericprofiles along the star-planet line as the global input to the 3D MHD model. Recent 3D hydrodynamic models obtain moreradially extended atmospheres along this direction compared to the terminator region (Shaikhislamov et al. 2018).

As one could expect from Fig. 13, we find insignificant differences of the absorption for the different wind scenarios due to the negligible contribution of ENAs in our model. This conclusion contradicts the one made by Kislyakova et al. (2014), who were able to reproduce the Lyα observations of HD 209458b only assuming a specific wind configuration. This contradiction can be easily explained by the fact that the Direct Simulation Monte Carlo model by Kislyakova et al. (2014) did not account for the deceleration and temperature increase of the ENAs near the planetary obstacle. For this reason, their results only accounted for a Maxwellian spectrum according to the initial density, temperature, and velocity distribution. On the contrary, our results represent a better approximation of the wind properties near the planet, and show that different stellar wind conditions can produce similar Lyα signatures. Differences to previous studies are discussed in more detail in Sect. 6.3.

One should keep in mind that we did not account for the compression and additional ionization of the planetary atmosphere by the stellar wind. Therefore, our results still represent an approximation, even though they present a significant improvement in comparison to earlier works by Holmström et al. (2008) and Kislyakova et al. (2014). However, these neglected processes could also be a reason for the slight overabsorption we obtained with our models.

thumbnail Fig. 11

Same as Fig. 8, but for the CME conditions.

Open with DEXTER
Table 4

Stand-off distances, ENA layer thickness, and maximum number densities of ENAs and ions in front of the magnetopause.

thumbnail Fig. 12

Velocity spectra of neutral H atoms for the atmosphere of HD 189733b for mean (blue), maximum (orange), and CME conditions (green). The main peak represents the atmospheric neutrals, the low and wide wings are the ENAs with different bulk velocities and high temperatures.

Open with DEXTER
thumbnail Fig. 13

Calculated transmissivity (blue) and corresponding absorption (orange) in the Lyα line as a function of velocity. The results for the mean wind parameters are shown, but the other two cases produce indistinguishable results due to the dominating contribution of atmospheric broadening.

Open with DEXTER
thumbnail Fig. 14

Modeled absorption of HD 189733b compared to observations from 2011 (Lecavelier des Etangs et al. 2012). The blue line shows the observed out-of-transit Lyα line profile of HD 189733. The green line shows the modeled in-transit absorption. The orange line shows the in-transit absorption observed in September 2011. The shaded region surrounding the orange line shows the observational errors. Our modeled spectrum reproduces the observed profile well (also in its red part) and is mostly within the error boundaries. Mean wind conditions were adopted, but the other cases produce identical results. The region between − 40 and + 40 km s−1 was affectedby geocoronal emission in 2011 and should be ignored (Bourrier & Lecavelier des Etangs 2013).

Open with DEXTER

6 Discussion

6.1 Comparison with other hydrodynamic models

As described above, our simulations yield a total mass-loss rate of = 2.5−5.4 × 1010 g s−1, depending on the adopted XUV spectrum. Previous studies of this planet found = 4.8 × 1010−2 × 1011 g s−1 for FXUV = 2 × 104−105 erg cm−2 s−1 (Guo 2011), = 4.5 × 1011−9 × 1011 g s−1 for FXUV = 24 778 erg cm−2 s−1 (but different spectral energy distributions; Guo & Ben-Jaffel 2016), and = 1.64 × 1010 g s−1 for FXUV = 20 893 erg cm−2 s−1 (Salz et al. 2016a)3. Despite using similar XUV fluxes, the model of Salz et al. (2016a) yielded a mass-loss rate lower by more than an order of magnitude compared tothe other studies (Guo 2011; Guo & Ben-Jaffel 2016), which is likely because not all relevant radiative cooling processes were included in the latter models, leading to an overestimate of the escape rate. Our mass-loss rates are up to a factor of three higher than those of Salz et al. (2016a). Estimating the spectral shape with a ratio of fluxes like in Guo & Ben-Jaffel (2016) we find β =F50−400ÅF50−900Å~F0−400ÅF0−912Å~0.8 for both spectra, corresponding to a mass-loss rate of ~ 9 × 1011 g s−1 in their model (see their Fig. 13), which is a factor of 15 higher than our results.

The atmospheric profiles are compared in Fig. 15. For Guo (2011) we adopted their results for FXUV = 2 × 104 erg cm−2 s−1 (their Fig. 3). One can see that our model yields lower velocities, but higher densities throughout the computational domain. Moreover, ionization occurs at a larger height in our model compared to their results. This could be partly related to the much lower densities n0 adopted in the other studies, as well as different XUV spectra. The temperature maximum is comparable to Salz et al. (2016a), but ours is located further out at about 2 Rp instead of 1.5 Rp.

The main differences between the discussed models are: gray atmosphere approximation, no molecules, only H, only Lyα cooling in Guo (2011); XUV spectra, molecules, H+He, but only H cooling in Guo & Ben-Jaffel (2016); XUV spectra, no molecules, H+He, all atomic radiative cooling processes in Salz et al. (2016a); XUV spectra, no molecules, only H, all atomic cooling processes in the present study. In addition, there are some differences in the adopted stellar and planetary parameters and lower boundary conditions, as well as usage of different solution methods of the hydrodynamic equations. Therefore, some differences in the mass-loss rates and atmospheric profiles are to be expected. We find the best qualitative agreement with Salz et al. (2016b) in that HD 189733b is an intermediate case between planets with high escape rates and stable planets in radiative equilibrium.

6.2 Neglected processes

There are several processes that are not considered in our present model that may have an impact on the planetary escape rates, and thus the modeled transit absorption. Firstly, usage of a 1D hydrodynamic code for the upper atmosphere limits the validity of the results to a region well within the Roche lobe. The average Roche lobe radius of HD 189733b is ~ 3.1 Rp (Eggleton 1983), smaller than the distance to the L1 point (cf., Sect. 2.3). Close to and beyond the Roche lobe, the gas dynamics are dominated by 3D effects (e.g., Bisikalo et al. 2013). However, we find that even for mean stellar wind conditions, the pressure balance distance is at about 2.6 Rp (well within the Roche lobe),which justifies using the 1D results for the atmosphere up to this point. However, we note that the sonic point is reached close to the L1 point above 4 Rp, which means that the dynamic pressure of the stellar wind at this point is larger than that of the planetary outflow. This means we could be overestimating the escape rates, because the stellar wind may confine the planetary mass-loss. However, it is possible that the confinement on the dayside is (partly) compensated by a stronger outflow on the nightside (Murray-Clay et al. 2009; Shaikhislamov et al. 2016). On the other hand, the absorbed XUV energy could also be radiated away by Lyα and free-free cooling due to the enhanced temperatures, which would lead to a reduction of the escape rates (Salz et al. 2016a). Self-consistent modeling of planetary-stellar wind interaction would require a multidimensional multi-fluid model.

Due to the strong ionization of hot Jupiter atmospheres, an intrinsic planetary magnetic field may lead to a reduction of the escape rate (Owen & Adams 2014; Trammell et al. 2014; Khodachenko et al. 2015). It is not known if HD 189733b possesses an intrinsic magnetic field strong enough to affect its mass loss. However, we do consider the generation of an induced magnetic field at the obstacle (magnetic barrier) in our MHD model (Erkaev et al. 2017). An intrinsic planetary magnetic field may push the pressure balance distance between the planetary and stellar winds further out.

We include only hydrogen atoms and protons in our simulations. Neglecting H2 and its related ions is justified for this planet, because their presence is confined to a very small region close to the optical radius (see Guo & Ben-Jaffel 2016, and Appendix C). This also means that cooling by H is not important here due to its small number density. However, helium is likely to be present in the atmosphere with cosmic abundances. Salz et al. (2016a) compared hydrodynamic simulations with and without including He and found mass-loss rates lower by a factor of two in the former case.

thumbnail Fig. 15

Atmospheric profiles for the SF11 run (solid) compared with previous model results of (Guo 2011; dash-dotted) and (Salz et al. 2016a; dashed).

Open with DEXTER

6.3 Comparison with previous transit signature modeling

Bourrier & Lecavelier des Etangs (2013) described the details of the 3D particle model, which they used to explain the Lyα observations of HD 189733b (Lecavelier des Etangs et al. 2012). It depends on the planetary mass-loss rate of neutral H, ionizing flux (limiting the lifetime of H atoms), and stellar wind parameters at the planet’s orbit (Vsw, Nsw, Tsw). It includes radiation pressure, planetary and stellar gravities, charge exchange with stellar wind protons, and self-shielding of stellar photons and stellar wind protons by the H cloud. They found their best fits to the 2011 observations with Vsw = 200 km s−1, Tsw = 3 × 104 K, Nsw = 103−3 × 105 cm−3, an ionizing flux ten times the solar one, and neutral H escape rates of 5 × 108−1.5 × 109 g s−1. To reproduce the nondetection in 2010, for fixed wind parameters the H escape rates would need to be 5–20 times lower, whereas for a fixed ionizing flux and escape rate, the stellar wind proton densities would need to be about a factor of ten lower. This indicates that the 2011 observations can be explained with a higher H escape rate (but not much higher ionizing flux) or a denser wind compared to the conditions in 2010. Lecavelier des Etangs et al. (2012) suggested that these discrepant observations could be due to a flare that occurred ~8 h prior to the transit in 2011.

The model of Bourrier & Lecavelier des Etangs (2013) can explain the observed Lyα absorption in 2011 if the neutral H mass-loss rate at 2.95 Rp (the lower boundary of their model) is in the order of 109 g s−1. Our H mass-loss rate at this altitude is comparable to this requirement, but similar for average XUV and flaring cases (see Sect. 3.2), even if ignoring the limited duration of the enhanced XUV flux during a flare. Hence, our neutral loss rates are consistent with their model for the 2011 observations, but they cannot explain the nondetection in 2010, because the XUV fluxes must have been very different. Either the fluxes were much lower, reducing the total mass-loss rate, or much higher, reducing the neutral H density due to enhanced ionization. This contradicts the necessity of similar ionizing fluxes at both epochs (Bourrier & Lecavelier des Etangs 2013). Test simulations with significantly enhanced XUV fluxes show that even an increase by a factor of 1000 only decreases the neutral H mass-loss rate by a factor of three. For lower XUV fluxes, the total mass-loss rate decreases, but the neutral H mass-loss rate increases relative to our results (e.g., about a factor of two larger for a factor of five lower XUV). However, we stress that other hydrodynamic models of this planet find higher ionization (and therefore lower H densities) at 2.95 Rp, yielding much lower neutral loss rates (Guo 2011; Salz et al. 2016a; Chadney et al. 2017). We find that this is mainly related to the lowernumber densities assumed at the lower boundary in these studies. If we run our model with lower n0 values, we also obtain neutral densities at 2.95 Rp that are too low compared to those required in Bourrier & Lecavelier des Etangs (2013). This indicates that the neutral density at a given altitude does not just depend on the irradiating stellar XUV flux, but strongly depends on n0, stronger than the total mass-loss rate. Higher values of n0 shift the neutral-ion transition to larger heights. However, we find that lowering n0 results in incomplete absorption of the XUV radiation in the computational domain, which makes the resulting mass-loss rates strongly dependent on the assumedn0.

One may argue that if the true value of n0 were indeed much lower than our assumptions, the absorption signal from the planetary atmosphere could be much weaker thanwhat our results show and not compatible with the 2011 observations. However, it is unlikely that n0 could be much smaller, because too much ionizing radiation would then reach layers close to the optical radius or even below, which is inconsistent with the observed value of the optical radius. The only possibility that the upper atmosphere could absorb most of the ionizing radiation above Rp despite a lower H density could be efficient absorption by non-hydrogen species (such as the absorption of X-rays by metals). Future studies should investigate whether this effect could be relevant for reasonable atmospheric metallicities.

The other possible explanation for the transit variability in the Bourrier & Lecavelier des Etangs (2013) model is a higher stellar wind density in 2011 by a factor of about ten compared to 2010. Our results are not consistent with this picture, as we find a negligible influence of the stellar wind on the planetary absorption, even in the case of a putative CME that may have accompanied the flare. The absorption in our model is completely dominated by the contribution from the parts of the planetary atmosphere below the pressure balance distance (Sect. 5). We also tested the effect of adopting stellar wind parameters more similar to the best fitting model of Bourrier & Lecavelier des Etangs (2013) in our MHD model, namely Vsw = 200 km s−1, Tsw = 3 × 104 K, Nsw = 106 cm−3, and keeping the magnetic field parameters of our mean wind case. We find that the maximum ENA density is much smaller than in the cases from Table 3, and that their bulk speed is low compared to the local thermal speed of the stellar wind. Therefore, such a wind does not produce a tail and a localized blue-shifted absorption signature in our model.

Previous modeling approaches employed Monte Carlo particle codes to model the generation and distribution of ENAs around exoplanets (Holmström et al. 2008; Bourrier & Lecavelier des Etangs 2013; Kislyakova et al. 2014). ENAs generated in such models by charge exchange of neutral planetary atoms with the stellar wind have a velocity distribution peaked at the stellar wind speed. After formation, ENAs are generally not coupled to the plasma flow due to lack of collisions; their motion is controlled by gravitational, centrifugal and Coriolis forces, as well as radiation pressure (e.g., Kislyakova et al. 2014). For HD 189733b, the situation is different. The high densities in the circumplanetary environment lead to an efficient coupling of ENAs and ions by collisions4. This justifies modeling ENAs in the HD 189733 system with our MHD code. In our model, the ENAs are generated inside the bow shock. These ENAs thus have velocities much smaller than that of the stellar wind because they are formed from decelerated wind ions. This leads to ENA velocity spectra that peak at small velocities (Fig. 12). Due to strong ionization by the intense XUV flux, the upper planetary atmosphere is almost completely ionized, so the generation of additional ENAs outside the bow shock (which would have velocities comparable to the stellar wind) can likely be neglected here due to the very small amount of available planetary neutrals.

As described in Sect. 6.1, Guo & Ben-Jaffel (2016) used a hydrodynamic upper atmosphere model similar to ours in combination with different spectral energy distributions of the stellar XUVradiation. They are able to reproduce the 2011 observations with a spectral shape β =0.76, which is very similar to our adopted SF11 spectrum (β = 0.82) that also reproduces the transit absorption. The nondetection in 2010 would require a spectral shape of β =0.38 (Guo & Ben-Jaffel 2016). Our results are consistent with theirs in that the 2011 observations can be reproduced by solely considering the planetary atmosphere, ENAs being negligible. However, due to the difficulties in observing the EUV emission of stars, it is unknown if such strong temporal changes in spectral shape would be possible while leaving the total XUV emission rather similar.

6.4 Effects of flares

As described before, with our adopted mean XUV spectrum we can reproduce the transit observations from 2011. Thus, no transient stellar variability, like a flare, is needed to produce the planetary absorption. As discussed in Sect. 3.2, the X-ray flare observed 8 h before the planetary transit in 2011 was unlikely able to induce any measurable difference in planetary absorption. Stellar emissions higher by a factor of four, either in X-rays only or XUV, raise the total mass-loss rates by 60 or 120%, respectively. In the latter case, however, the neutral H mass-loss rate drops because of increased ionization. Using a realistic time evolution of a flare and not just a constant elevation of the irradiation would make the effects even smaller. Our findings are in agreement with Chadney et al. (2017) who studied the effect of flares on the upper atmospheres of close-in exoplanets. They also found that a flare is likely not capable of producing the strong modulation of the transit absorption as seen in HD 189733b, but suggested that a stellar proton event could cause such differences.

6.5 Stellar activity effects on the Lyα transit absorption

Stellar activity can also affect the study of Lyα transit absorption in other ways than a modulation of the planet’s atmospheric properties. To obtain the amount of absorption during the transit, the in-transit spectrum is usually compared to an out-of-transit spectrum. This can cause difficulties in the case of an active star for two main reasons. Firstly, the surface of a star may not be homogeneous due to the presence of active regions. Secondly, the star may be variable on time scales comparable to the planetary transit, so the stellar background emission during the transit could be different from the stellar spectrum recorded at some time before the transit.

The first effect was studied in Llama & Shkolnik (2015), who simulated hot Jupiter transits over solar X-ray, EUV, FUV, and optical images. Both occulted and unocculted active regions can modify the obtained absorption. It also depends if the active regions are brighter or darker than the inactive stellar disk in a given wavelength range. For bright active regions, unocculted ones lead to shallower transits, and occulted ones to deeper transits. If a star has an activity belt like the Sun, it depends on where a planet transits the disk if the observed absorption is deeper or shallower than if assuming a homogeneous disk. For the case of HD 189733b, Llama & Shkolnik (2016) stated that other studies found no bumps in lightcurves, meaning no spot crossings, so the planet is unlikely to transit an active belt (in the case that the star actually has one). Since unocculted regions make a transit shallower, it is possible that the 2010 observations were affected by an increased number of unocculted active regions. However, we remind the reader that the planetary absorption in Lyα can only be studied in the line wings due to geocoronal emission and ISM absorption. The wings of the Lyα line show a much more homogeneous activity pattern than the core, since they are formed at lower temperatures deeper in the chromosphere (Salz et al. 2016a).

For the second effect, Llama & Shkolnik (2016) simulated transits in light curves of the disk-integrated solar Lyα line flux. Although the recovered planetary radii can be up to 50% larger for solar-like variability, they found that the effect is unlikely to be responsible for the strong transit absorptions detected on stars like HD 189733b, even if the solar variability is scaled up in amplitude.

7 Conclusions

We modeled the Lyα transit absorption of the hot Jupiter HD 189733b using a 1D hydrodynamic code for the upper atmosphere and a 3D MHD code for the planetary-stellar wind interaction and related production of ENAs. We find that the transit absorption observed in 2011 can be reproduced reasonably well with typical stellar XUV emissions (FXUV ~ 1.8 × 104 erg cm−2 s−1). We note that with our adopted default parameters, the model produces some extra-absorption across the Lyα line, which could be due to limitations of the 1D approach and/or the unaccounted compression of the atmosphere by the stellar wind. The influence of enhanced stellar irradiation during a flare similar to the one observed about 8 h before the transit in 2011 is too small to significantly affect the neutral hydrogen outflow rates. Moreover, we find with our modeling approach that the absorption signature is dominated by the atmospheric neutral hydrogen and that the ENAs produced by charge-exchange with the stellar wind have a negligibly small effect, as they are produced inside the bow shock and have small velocities and high temperatures due to the stellar wind deceleration at the planetary obstacle. Moreover, we find no differences in planetary absorption for all tested stellar wind cases, including a CME event. This raises the question if the nondetection in 2010 was actually “anomalous”, and not the absorption signature detected in 2011 as previously suggested. The variations between the two observations may be related to significant differences in stellar activity, by magnitude or spectral shape (Guo & Ben-Jaffel 2016), or its variability (Llama & Shkolnik 2016). Observations of this system in Lyα at a further epoch, preferentially with simultaneous X-ray monitoring to constrain the stellar XUV emission, would help to clarify this issue.

Acknowledgements

We thank the referee for valuable comments and suggestions. We thank A. Vidotto for providing the modeled stellar wind parameters. This work was supported by the Austrian Science Fund (FWF) project P27256-N27 “Characterizing Stellar and Exoplanetary Environments via Modeling of Lyman-α Transit Observations of Hot Jupiters”. The authors further acknowledge the support by the FWF NFN project S11601-N16 “Pathways to Habitability: From Disks to Active Stars, Planets and Life”, and the related FWF NFN subproject S11607-N16 “Particle/Radiative Interactions with Upper Atmospheres of Planetary Bodies Under Extreme Stellar Conditions” (K.K., H.L., N.V.E.), as well as FWF projects P30949-N36 (P.O., M.L.) and I2939-N27. L.F., D.K. and N.V.E. acknowledge also the Austrian Forschungsförderungsgesellschaft FFG project “TAPAS4CHEOPS” P853993. N.V.E. and V.A.I. acknowledge support by the Russian Science Foundation grant No 18-12-00080.

Appendix A MacCormack TVD scheme

The hydrodynamic model, as described in Erkaev et al. (2016), solves the set of 1D hydrodynamic equations using the MacCormack scheme (MacCormack 1969), which is second-order accurate in space and time. This scheme is well suited for supersonic outflows, but may experience instabilities and oscillations if the outflows are largely subsonic, as in the case of the high-gravityplanet HD 189733b. Several studies developed corrections to the original MacCormack scheme to give it total variation diminishing (TVD) properties (e.g., Davis 1984). This suppresses spurious numerical oscillations, which may occur at steep gradients and destabilize the code. Specifically, we adopted the method of Glaister (1991b) designed for nonuniform grids (Glaister 1991a).

The system of hydrodynamic equations (Eqs. (1)–(3)) was first normalized as described in Erkaev et al. (2016) and then written in vector form (A.1)

where U is the vector of the variables, F the fluxes, and S the source terms. Here, (A.2)

The MacCormack method consists of a predictor step: (A.3)

and a corrector step: (A.4)

where the indices i and n refer to spatial and temporal steps, respectively. The superscript n + 1∕2 refers to quantities evaluated after the predictor step. The time step was updated according to the Courant-Friedrichs-Lewy condition , where Δxmin is the smallest grid size, c is the adiabatic sound speed, is the maximum wave speed on thegrid, and CCFL ≤ 1 is the Courant-Friedrichs-Lewy number, which we took to be 0.8.

Here, we newly implemented the improved scheme of MacCormack (1971) to solve the momentum equation. This is related to the destabilizing effect of cases where the velocities ui < 0 and ui+1 > 0. Since the momentum flux ρu2 loses the information on the sign of u, MacCormack (1971) suggested, in such cases, replacing the momentum flux differences with in Eq. (A.3), and with correspondingly adjusted indices in Eq. (A.4). This preserves the sign of the flux and removes instabilities.

The TVD corrections from Glaister (1991b) are implemented as follows. After every corrector step, TVD correction terms are added to the solution, (A.5)

where , , and Δri = (ri+1ri−1)∕2. The coefficients are determined by (A.6)

where |λi±1∕2|max = |ui±1∕2| + ci±1∕2 is the maximum local wave speed and ψ(M) = max(0, min(1, M)) is the minmod limiter. The coefficients M± are calculated via (A.7)

where s denotes the sign function of the superscript5 and the expression (* ⋅ *) refers to the inner product of the difference vectors ΔU. The parameter (A.8)

is the ratio of successive grid spacings.

A.1 Modification of the momentum equation

Due to the highly subsonic nature of the outflow for planets with high gravity, the solution of the velocity profile is affected by unphysical behavior (e.g., negative values, u increasing towards the planet) close to the lower boundary of the computational domain. The situation improves if the code is run for very long times (because of the small time steps), but this is inconvenient. The problem is (partly) due to the solution method. We tried to solve the momentum equation in a nonconservative form, treating both the pressure and gravity terms as source terms, meaning the vector components (Eq. (A.2)) are modified to F2 = ρu2r2 and S2 = −ρr2∂Φ∂rr2∂p∂r. In the limit of u → 0, hydrostatic equilibrium should be reached, and the two terms in S2 should cancel each other out. Schemes that maintain hydrostatic equilibrium in this limit are called well-balanced schemes. We modified the discretization of the source terms according to the suggestions of Käppeli & Mishra (2016). Firstly, the derivatives of Φ and p in S2 are discretized the same way as the respective advection terms in Eqs. (A.3) and (A.4); secondly, the prefactors of these derivatives (ρr2, r2) are calculated as averages of their values at the grid points involved in the derivatives. For Eq. (A.3), this yields (A.9)

and similar for Eq. (A.4) with correspondingly adjusted indices. This procedure improved the resulting velocity profiles. We checked that this implementation did not affect the results by comparing runs for a lower gravity planet (where this modification is not necessary) using both the conservative and the modified formulation of the momentum equation.

A.2 Mass and energy conservation

Figure A.1 shows the mass flux ρur2 normalized to its value at the upper boundary as a function of r, as well as the total heating and cooling rates for the SF11 run. One can see that both mass and energy are well conserved in our model. We adopted the convergence criterion that the momentum flux should be conserved to 1% throughout the computational domain. This is only violated close to the lower boundary where the highly subsonic velocities introduce numerical difficulties for this planet. Also, the total heating and cooling rates agree well, indicating good energy conservation properties of the code, except close to the planet where the numerical errors lead to slightly lower adiabatic cooling than expected.

thumbnail Fig. A.1

Demonstration of the conservation properties our code in a converged solution (SF11 run). Left: normalized mass flux (blue), the dashed and dotted black lines indicate unity and ±1%, respectively. Right: total heating (red) and cooling rates (dashed blue).

Open with DEXTER

Appendix B Comparison of 1D and 2D heating

Figure B.1 shows a comparison between Eq. (13) and the 2D heating function (e.g., Erkaev et al. 2016). For the latter, it was assumed that the density profile along the star-planet line is valid for all other directions. This is approximately valid if the Roche lobe is far from the planet, and the planet is not tidally locked so that the stellar XUV flux isdistributed efficiently over the whole atmosphere. However, for a hot Jupiter, these assumptions are likely not valid, and a 3D model would be needed. Moreover, since we introduced a wavelength-dependent heating function, which raises the computational demand, we used a more simple 1D heating function. Adopting the neutral hydrogen density profile from the SF11 run, we calculated the 2D volume heating rate a posteriori and compared it with Eq. (13). One can see that for the adopted planetary and stellar parameters, α =4 represents the best match and was consequently used in all runs.

thumbnail Fig. B.1

Comparison of 2D and 1D XUV volume heating rates.

Open with DEXTER

Appendix C Molecular hydrogen

Here, we test if neglecting molecular hydrogen could have an impact on the results. We calculated the equilibrium abundance of H2, which can be justified considering that it will only be abundant close to Rp where velocities are very small. The H2 density can be written as (C.1)

where is the rate of the reaction producing H2 (H+H →H2) and νdiss = 1.5 × 10−9 exp(−49 000∕T) the thermaldissociation rate (Yelle 2004; Erkaev et al. 2016). The photoionization rate is calculated after Eq. (6), but replacing σion with that for H2 (Huebner & Mukherjee 2015). The resulting number density profile of H2, along with the densities of H atoms and protons for the SF11 run are shown in Fig. C.1. One can see that H2 would only be the dominating species at distances ≲1.04 Rp, very close to the optical radius. The main part of the upper atmosphere is composed of H and H+. This is consistent with Guo & Ben-Jaffel (2016), who found the same H2 /H transition radius with their hydrodynamic model (which does include hydrogen chemistry) and that the planetary wind is dominated by H and H+. Therefore, we do not expect that inclusion of H2 would noticeably change our results.

thumbnail Fig. C.1

Number densities of H2, H, and H+, where the latter two are taken from the SF11 run and the H2 density is calculated with Eq. (C.1).

Open with DEXTER

Appendix D Stellar wind parameters

In Fig. D.1, we show the stellar wind parameters along the orbit based on the measured magnetic field map (Fares et al. 2010) and a 3D wind model (Llama et al. 2013).

thumbnail Fig. D.1

Variation of the stellar wind parameters along the orbit of HD 189733b. Shown are the number density, velocity, temperature, strength of the total magnetic field and its components (parallel and normal to Vsw), the angle between the magnetic field and the star-planet line, and the total pressure with its components (thermal, magnetic, and ram pressures).

Open with DEXTER

References

  1. Adams, F. C. 2011, ApJ, 730, 27 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alvarado-Gómez, J. D., Drake, J. J., Cohen, O., Moschou, S. P., & Garraffo, C. 2018, ApJ, 862, 93 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bakos, G. Á., Pál, A., Latham, D. W., Noyes, R. W., & Stefanik, R. P. 2006, ApJ, 641, L57 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ballester, G. E., & Ben-Jaffel, L. 2015, ApJ, 804, 116 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ballester, G. E., Sing, D. K., & Herbert, F. 2007, Nature, 445, 511 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  6. Barnes, J. R., Haswell, C. A., Staab, D., & Anglada-Escudé, G. 2016, MNRAS, 462, 1012 [NASA ADS] [CrossRef] [Google Scholar]
  7. Ben-Jaffel, L. 2007, ApJ, 671, L61 [NASA ADS] [CrossRef] [Google Scholar]
  8. Ben-Jaffel, L. 2008, ApJ, 688, 1352 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ben-Jaffel, L., & Ballester, G. E. 2013, A&A, 553, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Ben-Jaffel, L., & Sona Hosseini, S. 2010, ApJ, 709, 1284 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bisikalo, D., Kaygorodov, P., Ionov, D., et al. 2013, ApJ, 764, 19 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bouchy, F., Udry, S., Mayor, M., et al. 2005, A&A, 444, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bourrier, V., & Lecavelier des Etangs, A. 2013, A&A, 557, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bourrier, V., Lecavelier des Etangs, A., Dupuy, H., et al. 2013, A&A, 551, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bourrier, V., Ehrenreich, D., & Lecavelier des Etangs, A. 2015, A&A, 582, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Cauley, P. W., Redfield, S., Jensen, A. G., et al. 2015, ApJ, 810, 13 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cauley, P. W., Redfield, S., Jensen, A. G., & Barman, T. 2016, AJ, 152, 20 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cauley, P. W., Redfield, S., & Jensen, A. G. 2017a, AJ, 153, 217 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cauley, P. W., Redfield, S., & Jensen, A. G. 2017b, AJ, 153, 185 [NASA ADS] [CrossRef] [Google Scholar]
  21. Chadney, J. M., Galand, M., Unruh, Y. C., Koskinen, T. T., & Sanz-Forcada, J. 2015, Icarus, 250, 357 [NASA ADS] [CrossRef] [Google Scholar]
  22. Chadney, J. M., Koskinen, T. T., Galand, M., Unruh, Y. C., & Sanz-Forcada, J. 2017, A&A, 608, A75 [CrossRef] [EDP Sciences] [Google Scholar]
  23. Christie, D., Arras, P., & Li, Z.-Y. 2016, ApJ, 820, 3 [NASA ADS] [CrossRef] [Google Scholar]
  24. Cohen, O., Kashyap, V. L., Drake, J. J., et al. 2011, ApJ, 733, 67 [NASA ADS] [CrossRef] [Google Scholar]
  25. Davis, S. F. 1984, ICASE Report, 84-20, 1 [Google Scholar]
  26. de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Debrecht, A., Carroll-Nellenback, J., Frank, A., et al. 2018, MNRAS, 478, 2592 [CrossRef] [Google Scholar]
  28. Drake, J. J., Cohen, O., Yashiro, S., & Gopalswamy, N. 2013, ApJ, 764, 170 [NASA ADS] [CrossRef] [Google Scholar]
  29. Eggleton, P. P. 1983, ApJ, 268, 368 [NASA ADS] [CrossRef] [Google Scholar]
  30. Ehrenreich, D., Lecavelier Des Etangs, A., Hébrard, G., et al. 2008, A&A, 483, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Ehrenreich, D., Bourrier, V., Bonfils, X., et al. 2012, A&A, 547, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459 [NASA ADS] [CrossRef] [Google Scholar]
  33. Ekenbäck, A., Holmström, M., Wurz, P., et al. 2010, ApJ, 709, 670 [NASA ADS] [CrossRef] [Google Scholar]
  34. Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Erkaev, N. V., Lammer, H., Odert, P., et al. 2013, Astrobiology, 13, 1011 [NASA ADS] [CrossRef] [Google Scholar]
  36. Erkaev, N. V., Lammer, H., Odert, P., Kulikov, Y. N., & Kislyakova, K. G. 2015, MNRAS, 448, 1916 [NASA ADS] [CrossRef] [Google Scholar]
  37. Erkaev, N. V., Lammer, H., Odert, P., et al. 2016, MNRAS, 460, 1300 [NASA ADS] [CrossRef] [Google Scholar]
  38. Erkaev, N. V., Odert, P., Lammer, H., et al. 2017, MNRAS, 470, 4330 [NASA ADS] [CrossRef] [Google Scholar]
  39. Fares, R., Donati, J.-F., Moutou, C., et al. 2010, MNRAS, 406, 409 [NASA ADS] [CrossRef] [Google Scholar]
  40. Fossati, L., Haswell, C. A., Froning, C. S., et al. 2010, ApJ, 714, L222 [NASA ADS] [CrossRef] [Google Scholar]
  41. France, K., Froning, C. S., Linsky, J. L., et al. 2013, ApJ, 763, 149 [NASA ADS] [CrossRef] [Google Scholar]
  42. Gaia Collaboration (Brown, A. G. A., et al. ) 2016, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. García Muñoz, A. 2007, Planet. Space Sci., 55, 1426 [NASA ADS] [CrossRef] [Google Scholar]
  44. Glaister, P. 1991a, Comput. Math. Appl., 21, 39 [CrossRef] [Google Scholar]
  45. Glaister, P. 1991b, Comput. Math. Appl., 22, 45 [CrossRef] [Google Scholar]
  46. Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1 [NASA ADS] [CrossRef] [Google Scholar]
  47. Gray, R. O., Corbally, C. J., Garrison, R. F., McFadden, M. T., & Robinson, P. E. 2003, AJ, 126, 2048 [NASA ADS] [CrossRef] [Google Scholar]
  48. Güdel, M., Audard, M., Kashyap, V. L., Drake, J. J., & Guinan, E. F. 2003, ApJ, 582, 423 [NASA ADS] [CrossRef] [Google Scholar]
  49. Guo, J. H. 2011, ApJ, 733, 98 [NASA ADS] [CrossRef] [Google Scholar]
  50. Guo, J. H. 2013, ApJ, 766, 102 [NASA ADS] [CrossRef] [Google Scholar]
  51. Guo, J. H., & Ben-Jaffel, L. 2016, ApJ, 818, 107 [NASA ADS] [CrossRef] [Google Scholar]
  52. Haswell, C. A., Fossati, L., Ayres, T., et al. 2012, ApJ, 760, 79 [NASA ADS] [CrossRef] [Google Scholar]
  53. Henry, G. W., & Winn, J. N. 2008, AJ, 135, 68 [NASA ADS] [CrossRef] [Google Scholar]
  54. Holmström, M., Ekenbäck, A., Selsis, F., et al. 2008, Nature, 451, 970 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  55. Huebner, W. F., & Mukherjee, J. 2015, Planet. Space Sci., 106, 11 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hünsch, M., Schmitt, J. H. M. M., Sterzik, M. F., & Voges, W. 1999, A&AS, 135, 319 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Ionov, D., & Shematovich, V. 2015, Sol. Sys. Res., 49, 339 [CrossRef] [Google Scholar]
  58. Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86 [NASA ADS] [CrossRef] [Google Scholar]
  59. Käppeli, R., & Mishra, S. 2016, A&A, 587, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Khodachenko, M. L., Lammer, H., Lichtenegger, H. I. M., et al. 2007, Planet. Space Sci., 55, 631 [NASA ADS] [CrossRef] [Google Scholar]
  61. Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., & Prokopov, P. A. 2015, ApJ, 813, 50 [NASA ADS] [CrossRef] [Google Scholar]
  62. Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2017, ApJ, 847, 126 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kislyakova, K. G., Holmström, M., Lammer, H., Odert, P., & Khodachenko, M. L. 2014, Science, 346, 981 [NASA ADS] [CrossRef] [Google Scholar]
  64. Kohl, S., Salz, M., Czesla, S., & Schmitt, J. H. M. M. 2018, A&A, 619, A96 [CrossRef] [EDP Sciences] [Google Scholar]
  65. Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013, Icarus, 226, 1678 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kulow, J. R., France, K., Linsky, J., & Loyd, R. O. P. 2014, ApJ, 786, 132 [NASA ADS] [CrossRef] [Google Scholar]
  67. Lara, A., Gonzalez-Esparza, J. A., & Gopalswamy, N. 2004, Geofis. Int., 43, 75 [Google Scholar]
  68. Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Lecavelier des Etangs, A., Ehrenreich, D., Vidal-Madjar, A., et al. 2010, A&A, 514, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Lecavelier des Etangs, A., Bourrier, V., Wheatley, P. J., et al. 2012, A&A, 543, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Leitzinger, M., Odert, P., Greimel, R., et al. 2014, MNRAS, 443, 898 [NASA ADS] [CrossRef] [Google Scholar]
  72. Leitzinger, M., Odert, P., Greimel, R., et al. 2020, MNRAS, 493, 4570 [CrossRef] [Google Scholar]
  73. Lindsay, B. G., & Stebbings, R. F. 2005, J. Geophys. Res. Space Phys., 110, A12213 [NASA ADS] [CrossRef] [Google Scholar]
  74. Linsky, J. L., Yang, H., France, K., et al. 2010, ApJ, 717, 1291 [NASA ADS] [CrossRef] [Google Scholar]
  75. Linsky, J. L., Fontenla, J., & France, K. 2014, ApJ, 780, 61 [NASA ADS] [CrossRef] [Google Scholar]
  76. Llama, J., & Shkolnik, E. L. 2015, ApJ, 802, 41 [NASA ADS] [CrossRef] [Google Scholar]
  77. Llama, J., & Shkolnik, E. L. 2016, ApJ, 817, 81 [NASA ADS] [CrossRef] [Google Scholar]
  78. Llama, J., Vidotto, A. A., Jardine, M., et al. 2013, MNRAS, 436, 2179 [NASA ADS] [CrossRef] [Google Scholar]
  79. MacCormack, R. W. 1969, AIAA, Paper 69–354 [Google Scholar]
  80. MacCormack, R. W. 1971, Lecture Notes in Physics (Berlin: Springer Verlag), 8, 151 [CrossRef] [Google Scholar]
  81. Marin, F., & Grosso, N. 2017, ApJ, 835, 283 [NASA ADS] [CrossRef] [Google Scholar]
  82. Matsakos, T., Uribe, A., & Königl, A. 2015, A&A, 578, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Menager, H., Barthélemy, M., Koskinen, T., et al. 2013, Icarus, 226, 1709 [NASA ADS] [CrossRef] [Google Scholar]
  84. Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23 [NASA ADS] [CrossRef] [Google Scholar]
  85. Odert, P., Leitzinger, M., Hanslmeier, A., & Lammer, H. 2017, MNRAS, 472, 876 [NASA ADS] [CrossRef] [Google Scholar]
  86. Odert, P., Leitzinger, M., Guenther, E. W., & Heinzel, P. 2020, MNRAS, 494, 3766 [CrossRef] [Google Scholar]
  87. Owen, J. E., & Adams, F. C. 2014, MNRAS, 444, 3761 [NASA ADS] [CrossRef] [Google Scholar]
  88. Patsourakos, S., & Georgoulis, M. K. 2016, A&A, 595, A121 [CrossRef] [EDP Sciences] [Google Scholar]
  89. Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton: Princeton University Press) [Google Scholar]
  90. Penz, T., Erkaev, N. V., Kulikov, Y. N., et al. 2008, Planet. Space Sci., 56, 1260 [NASA ADS] [CrossRef] [Google Scholar]
  91. Petit, P., Louge, T., Théado, S., et al. 2014, PASP, 126, 469 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pillitteri, I., Wolk, S. J., Cohen, O., et al. 2010, ApJ, 722, 1216 [NASA ADS] [CrossRef] [Google Scholar]
  93. Pillitteri, I., Günther, H. M., Wolk, S. J., Kashyap, V. L., & Cohen, O. 2011, ApJ, 741, L18 [NASA ADS] [CrossRef] [Google Scholar]
  94. Pillitteri, I., Wolk, S. J., Lopez-Santiago, J., et al. 2014, ApJ, 785, 145 [NASA ADS] [CrossRef] [Google Scholar]
  95. Pillitteri, I., Maggio, A., Micela, G., et al. 2015, ApJ, 805, 52 [NASA ADS] [CrossRef] [Google Scholar]
  96. Poppenhaeger,K., & Wolk, S. J. 2014, A&A, 565, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Poppenhaeger, K., Schmitt, J. H. M. M., & Wolk, S. J. 2013, ApJ, 773, 62 [NASA ADS] [CrossRef] [Google Scholar]
  98. Route, M. 2019, ApJ, 872, 79 [CrossRef] [Google Scholar]
  99. Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016a, A&A, 586, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2016b, A&A, 585, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Sanz-Forcada, J., García-Álvarez, D., Velasco, A., et al. 2010, IAU Symp., 264, 478 [Google Scholar]
  102. Sanz-Forcada, J., Micela, G., Ribas, I., et al. 2011, A&A, 532, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Schwartz, J. C., Kashner, Z., Jovmir, D., & Cowan, N. B. 2017, ApJ, 850, 154 [NASA ADS] [CrossRef] [Google Scholar]
  104. Schwenn, R., Raymond, J. C., Alexander, D., et al. 2006, Space Sci. Rev., 123, 127 [CrossRef] [Google Scholar]
  105. Sekiya, M., Nakazawa, K., & Hayashi, C. 1980, Prog. Theor. Phys., 64, 1968 [NASA ADS] [CrossRef] [Google Scholar]
  106. Semelin, B., Combes, F., & Baek, S. 2007, A&A, 474, 365 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Shaikhislamov, I. F., Khodachenko, M. L., Sasunov, Y. L., et al. 2014, ApJ, 795, 132 [NASA ADS] [CrossRef] [Google Scholar]
  108. Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2016, ApJ, 832, 173 [NASA ADS] [CrossRef] [Google Scholar]
  109. Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2018, MNRAS, 481, 5315 [NASA ADS] [CrossRef] [Google Scholar]
  110. Shematovich, V. I., Ionov, D. E., & Lammer, H. 2014, A&A, 571, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136 [NASA ADS] [CrossRef] [Google Scholar]
  112. Stone, J. M., & Proga, D. 2009, ApJ, 694, 205 [NASA ADS] [CrossRef] [Google Scholar]
  113. Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, ApJ, 621, 1049 [NASA ADS] [CrossRef] [Google Scholar]
  114. Tilley, M. A., Harnett, E. M., & Winglee, R. M. 2016, ApJ, 827, 77 [CrossRef] [Google Scholar]
  115. Trammell, G. B., Arras, P., & Li, Z.-Y. 2011, ApJ, 728, 152 [NASA ADS] [CrossRef] [Google Scholar]
  116. Trammell, G. B., Li, Z.-Y., & Arras, P. 2014, ApJ, 788, 161 [NASA ADS] [CrossRef] [Google Scholar]
  117. Tremblin, P., & Chiang, E. 2013, MNRAS, 428, 2565 [NASA ADS] [CrossRef] [Google Scholar]
  118. Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487 [NASA ADS] [CrossRef] [Google Scholar]
  119. Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J.-M., et al. 2003, Nature, 422, 143 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  120. Vidal-Madjar, A., Désert, J., Lecavelier des Etangs, A., et al. 2004, ApJ, 604, L69 [NASA ADS] [CrossRef] [Google Scholar]
  121. Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J., et al. 2008, ApJ, 676, L57 [NASA ADS] [CrossRef] [Google Scholar]
  122. Vidal-Madjar, A., Huitson, C. M., Bourrier, V., et al. 2013, A&A, 560, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Waite, J. H., Cravens, T. E., Kozyra, J., et al. 1983, J. Geophys. Res., 88, 6143 [NASA ADS] [CrossRef] [Google Scholar]
  124. Wood, B. E. 2004, Liv. Rev. Sol. Phys., 1, 2 [Google Scholar]
  125. Yashiro, S., & Gopalswamy, N. 2009, IAU Symp., 257, 233 [NASA ADS] [Google Scholar]
  126. Yashiro, S., Gopalswamy, N., Michalek, G., et al. 2004, J. Geophys. Res., 109, A07105 [NASA ADS] [CrossRef] [Google Scholar]
  127. Yelle, R. V. 2004, Icarus, 170, 167 [NASA ADS] [CrossRef] [Google Scholar]

1

The case B recombination coefficient takes into account that photons emitted by recombinations to the ground level immediately ionize a nearby atom and do not effectively contribute to recombination.

3

Their given value of 4.1 × 109 g s−1 corresponds to 1/4 of the isotropic mass-loss rate.

4

Taking nNsw, σcol ~ 10−15 cm2 and X ~ Rp yields a Knudsen number Kn < 1 in the modeling domain.

5

Meaning s = 1 for superscript “+” and s = −1 for superscript “−”.

All Tables

Table 1

Stellar and planetary parameters from Stassun et al. (2017). The orbital distance a was calculated from their quoted value of aR* and R*.

Table 2

Modeled mass-loss rates for different stellar XUV fluxes/spectra and lower boundary parameters (number density n0, temperature T0).

Table 3

Adopted stellar wind parameters at the planet’s orbit.

Table 4

Stand-off distances, ENA layer thickness, and maximum number densities of ENAs and ions in front of the magnetopause.

All Figures

thumbnail Fig. 1

XUV (5–912 Å) spectral flux at the orbit of HD 189733b in bins of 100 Å, obtained with the methods of L14 and SF11. The X-ray flux is an average value outside of flares obtained with XMM-Newton in 2011 (Pillitteri et al. 2014).

Open with DEXTER
In the text
thumbnail Fig. 2

Planetary atmospheric profiles adopting the XUV spectra of SF11 (solid) and L14 (dashed). Shown are the number densities (total, atomic H, and protons), outflow velocity, ionization fraction, and temperature. In the velocity plot, crosses indicate the sonic points and plus signs denote the locations where the outflow velocity exceeds the escape velocity. The vertical dotted line indicates the location of the L1 point.

Open with DEXTER
In the text
thumbnail Fig. 3

All heating (red) and cooling (blue) processes included in our model, shown for the SF11 run. The cyan lines detail the individual contributions to the radiative cooling rate (solid blue line).

Open with DEXTER
In the text
thumbnail Fig. 4

Heating efficiencies ηXUV(r) (solid) and ηnet(r) (dashed).

Open with DEXTER
In the text
thumbnail Fig. 5

Atmospheric profiles for average XUV (solid), X-ray flare (dashed), and XUV flare (dash-dotted) cases based on the SF11 spectrum.

Open with DEXTER
In the text
thumbnail Fig. 6

Cut at Y = 0 of the simulation showing the distribution of the total pressure (sum of the magnetic and gas pressures) around HD 189733b normalized to the stellar wind dynamic pressure for the mean stellar wind parameters. The white area close to the origin indicates the atmospheric region around the planet, while the semicircle indicates the planetary optical radius. The star is located along the X-axis. The arrow shows the direction of the interplanetary magnetic field.

Open with DEXTER
In the text
thumbnail Fig. 7

Same as Fig. 6, but for the maximum stellar wind.

Open with DEXTER
In the text
thumbnail Fig. 8

Profiles of the parameters along the stagnation line from the magnetopause to the bow shock for the mean stellar wind conditions. Top panel: ENA and ion densities. Bottom panel: total pressure and temperature.

Open with DEXTER
In the text
thumbnail Fig. 9

Same as Fig. 8, but for the maximum stellar wind conditions.

Open with DEXTER
In the text
thumbnail Fig. 10

Same as Fig. 6, but for the CME conditions.

Open with DEXTER
In the text
thumbnail Fig. 11

Same as Fig. 8, but for the CME conditions.

Open with DEXTER
In the text
thumbnail Fig. 12

Velocity spectra of neutral H atoms for the atmosphere of HD 189733b for mean (blue), maximum (orange), and CME conditions (green). The main peak represents the atmospheric neutrals, the low and wide wings are the ENAs with different bulk velocities and high temperatures.

Open with DEXTER
In the text
thumbnail Fig. 13

Calculated transmissivity (blue) and corresponding absorption (orange) in the Lyα line as a function of velocity. The results for the mean wind parameters are shown, but the other two cases produce indistinguishable results due to the dominating contribution of atmospheric broadening.

Open with DEXTER
In the text
thumbnail Fig. 14

Modeled absorption of HD 189733b compared to observations from 2011 (Lecavelier des Etangs et al. 2012). The blue line shows the observed out-of-transit Lyα line profile of HD 189733. The green line shows the modeled in-transit absorption. The orange line shows the in-transit absorption observed in September 2011. The shaded region surrounding the orange line shows the observational errors. Our modeled spectrum reproduces the observed profile well (also in its red part) and is mostly within the error boundaries. Mean wind conditions were adopted, but the other cases produce identical results. The region between − 40 and + 40 km s−1 was affectedby geocoronal emission in 2011 and should be ignored (Bourrier & Lecavelier des Etangs 2013).

Open with DEXTER
In the text
thumbnail Fig. 15

Atmospheric profiles for the SF11 run (solid) compared with previous model results of (Guo 2011; dash-dotted) and (Salz et al. 2016a; dashed).

Open with DEXTER
In the text
thumbnail Fig. A.1

Demonstration of the conservation properties our code in a converged solution (SF11 run). Left: normalized mass flux (blue), the dashed and dotted black lines indicate unity and ±1%, respectively. Right: total heating (red) and cooling rates (dashed blue).

Open with DEXTER
In the text
thumbnail Fig. B.1

Comparison of 2D and 1D XUV volume heating rates.

Open with DEXTER
In the text
thumbnail Fig. C.1

Number densities of H2, H, and H+, where the latter two are taken from the SF11 run and the H2 density is calculated with Eq. (C.1).

Open with DEXTER
In the text
thumbnail Fig. D.1

Variation of the stellar wind parameters along the orbit of HD 189733b. Shown are the number density, velocity, temperature, strength of the total magnetic field and its components (parallel and normal to Vsw), the angle between the magnetic field and the star-planet line, and the total pressure with its components (thermal, magnetic, and ram pressures).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.