Modeling the Ly$\alpha$ transit absorption of the hot Jupiter HD 189733b

TL;DR

This study models the Lyα transit absorption of hot Jupiter HD 189733b, showing that atmospheric broadening dominates the absorption signal and that typical stellar activity conditions can explain observed variations.

Contribution

The paper combines 1D hydrodynamic and 3D MHD models to analyze how stellar XUV and wind variations affect Lyα absorption, highlighting the dominance of atmospheric broadening over ENA contributions.

Findings

Lyα absorption can be explained by typical stellar activity conditions.

Atmospheric broadening dominates the Lyα transit signature.

ENAs have negligible impact on the absorption in this model.

Abstract

Hydrogen-dominated atmospheres of hot exoplanets expand and escape due to the intense heating by the X-ray and extreme ultraviolet (XUV) irradiation of their host stars. Excess absorption of neutral hydrogen has been observed in the Ly$\alpha$ line during transits of several close-in exoplanets, indicating such extended atmospheres. For the hot Jupiter HD 189733b, this absorption shows temporal variability. Variations in stellar XUV emission and/or stellar wind conditions have been invoked to explain this effect. We apply a 1D hydrodynamic upper atmosphere model and a 3D MHD stellar wind flow model to study the effect of variations of the stellar XUV and wind conditions on the neutral hydrogen distribution, including the production of energetic neutral atoms (ENAs), and the related Ly$\alpha$ transit signature. We obtain comparable, albeit slightly higher Ly$\alpha$ absorption as observed in 2011 with a stellar XUV flux of $1.8\times10^4$ erg cm$^{-2}$ s$^{-1}$, rather typical activity conditions for this star. Flares similar to the one observed 8 h before the transit are unlikely to have caused a significant modulation of the transit signature. The resulting Ly$\alpha$ absorption is dominated by atmospheric broadening, whereas the contribution of ENAs is negligible, as they are formed inside the bow shock from decelerated wind ions that are heated to high temperatures. Thus, within our modeling framework and assumptions, we find an insignificant dependence on the stellar wind parameters. Since the transit absorption can be modeled with typical stellar XUV and wind conditions, it is possible that the non-detection of the absorption in 2010 was affected by less typical stellar activity conditions, such as a very different magnitude and/or shape of the star's spectral XUV emission, or temporal/spatial variations in Ly$\alpha$ affecting the determination of the transit absorption.

Figures

Tables

Table 1. Table 1: Stellar and planetary parameters from Stassun et al. ( 2017 ) . The orbital distance a 𝑎 a was calculated from their quoted value of a / R ∗ 𝑎 subscript 𝑅 a/R_{*} and R ∗ subscript 𝑅 R_{*} .

Parameter	Value
$d$ (pc)	$19.84\pm 0.09$
$M_{\ast }$ (M☉)	$0.79\pm 0.08$
$R_{\ast }$ (R☉)	$0.75\pm 0.01$
$F_{\mathrm{bol}}$ (erg cm-2 s-1)	$2.68\times {10}^{-8}\pm 3.94\times {10}^{-10}$
$M_p$ (MJup)	$1.13\pm 0.08$
$R_p$ (RJup)	$1.13\pm 0.01$
$a$ (AU)	$0.031\pm 0.001$
$i$ (°)	85.71

Table 2. Table 2: Modeled mass-loss rates M ˙ ˙ 𝑀 \dot{M} for different stellar XUV fluxes/spectra and lower boundary parameters (number density n 0 subscript 𝑛 0 n_{0} , temperature T 0 subscript 𝑇 0 T_{0} ). Rows 3 and 4 give the results for the X-ray and XUV flares, respectively (see Section  3.2 ).

$F_{\mathrm{XUV}}$	XUV	$n_0$	$T_0$	$\dot{M}$
(erg cm-2 s-1)		(cm-3)	(K)	(g s-1)
$1.80\times {10}^4$	SF11	${10}^{15}$	1200	$5.4\times {10}^{10}$
$1.27\times {10}^4$	L14	${10}^{15}$	1200	$2.5\times {10}^{10}$
$3.65\times {10}^4$	SF11+X	${10}^{15}$	1200	$8.7\times {10}^{10}$
$7.18\times {10}^4$	SF11+XUV	${10}^{15}$	1200	$1.2\times {10}^{11}$
$1.80\times {10}^4$	SF11	${10}^{15}$	800	$4.8\times {10}^{10}$
$1.80\times {10}^4$	SF11	${10}^{15}$	1600	$6.1\times {10}^{10}$
$1.80\times {10}^4$	SF11	$5\times {10}^{15}$	1200	$1.7\times {10}^{11}$

Table 3. Table 3: Adopted stellar wind parameters at the planet’s orbit.

Parameter	mean	max	CME
$N_{\mathrm{sw}}$ (cm-3)	$4.4\times {10}^5$	$4.9\times {10}^5$	$4.4\times {10}^6$
$V_{\mathrm{sw}}$ (km s-1)	326	472	1000
$T_{\mathrm{sw}}$ (K)	$1.3\times {10}^6$	$2\times {10}^6$	$2\times {10}^6$
$B_{\mathrm{sw},p}$ (mG)	44	23	0
$B_{\mathrm{sw},n}$ (mG)	12	9.5	100
$B_{\mathrm{sw},\mathrm{tot}}$ (mG)	46	25	100
${\theta }_B$ (°)	15	22	90
$P_{\mathrm{dsw}}$ (dyn cm-2)	$7.82\times {10}^{-4}$	$1.83\times {10}^{-3}$	$7.36\times {10}^{-2}$

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

Wind	stand-off	ENA layer	ENA density	ion density
	distance	thickness	maximum	maximum
	$\left(R_p\right)$	$\left(R_p\right)$	(cm-3)	(cm-3)
mean	2.6	1.3	$2.3\times {10}^7$	$1.0\times {10}^7$
max	2.4	1.2	$5.6\times {10}^7$	$1.13\times {10}^7$
CME	2	0.25	$1.5\times {10}^8$	$0.67\times {10}^7$

Equations

$$\frac{\partial \rho r^2}{\partial t}+\frac{\partial \rho u{r}^2}{\partial r}$$

$$\frac{\partial \rho u{r}^2}{\partial t}+\frac{\partial \left(\rho u^2+p\right)r^2}{\partial r}$$

$$\frac{\partial E_{\mathrm{th}}{r}^2}{\partial t}+\frac{\partial E_{\mathrm{th}}u{r}^2}{\partial r}$$

$$\Phi =-\frac{G{M}_{\mathrm{p}}}{r}-\frac{G{M}_{\ast }}{a-r}-\frac{G\left(M_{\mathrm{p}}+M_{\ast }\right)}{2a^3}{\left(a\frac{M_{\ast }}{M_{\mathrm{p}}+M_{\ast }}-r\right)}^2$$

$$\frac{\partial n_{\mathrm{H}+}}{\partial t}+\frac{1}{r^2}\frac{\partial n_{\mathrm{H}+}u{r}^2}{\partial r}={\alpha }_{\mathrm{ion}}{n}_{\mathrm{H}}-{\alpha }_{\mathrm{rec}}{n}_{\mathrm{H}+}{n}_{\mathrm{e}}+{\alpha }_{\mathrm{col}}{n}_{\mathrm{H}}{n}_{\mathrm{e}},$$

$${\alpha }_{\mathrm{ion}}(r)=\sum _{\lambda }{\sigma }_{\mathrm{ion},\lambda }\frac{F_{\mathrm{XUV},\lambda }{\mathrm{e}}^{-{\tau }_{\lambda }(r)}}{E_{\lambda }},$$

$${\alpha }_{\mathrm{rec}}=2.753\times 10^{-14}{\left(\frac{315614}{T}\right)}^{1.5}{\left[1+{\left(\frac{115188}{T}\right)}^{0.407}\right]}^{-2.242},$$

$${\alpha }_{\mathrm{col}}$$

$$-5.73932875\times 10^0{\theta }^2+1.56315498\times 10^0{\theta }^3$$

$$-2.87705600\times 10^{-1}{\theta }^4+3.48255977\times 10^{-2}{\theta }^5$$

$$-2.63197617\times 10^{-3}{\theta }^6+1.11954395\times 10^{-4}{\theta }^7$$

$$-2.03914985\times 10^{-6}{\theta }^8),$$

$$Q_{\mathrm{XUV}}(r)=\sum _{\lambda }{\eta }_{\mathrm{ph}}{\sigma }_{\mathrm{ion},\lambda }{n}_{\mathrm{H}}(r)\frac{(E_{\lambda }-13.6\mathrm{eV})}{E_{\lambda }}{F}_{\mathrm{XUV},\lambda }\frac{\exp (-{\tau }_{\lambda }(r))}{1+\alpha {\tau }_{\lambda }(r)},$$

$${\tau }_{\lambda }(r)={\int }_r^{\infty }{\sigma }_{\mathrm{ion},\lambda }{n}_{\mathrm{H}}(r)dr.$$

$${\Lambda }_{\mathrm{Ly}\alpha }=0.1\times 7.5\times 10^{-19}{\left(1+\sqrt{\frac{T}{10^5}}\right)}^{-1}\exp \left(-\frac{118348}{T}\right)n_{\mathrm{e}}{n}_{\mathrm{H}}$$

$${\Lambda }_{\mathrm{rec}}=1.38\times 10^{-16}T{\alpha }_{\mathrm{rec}}{n}_{\mathrm{e}}{n}_{\mathrm{H}+},$$

$${\Lambda }_{\mathrm{col}}=2.179\times 10^{-11}{\alpha }_{\mathrm{col}}{n}_{\mathrm{e}}{n}_{\mathrm{H}},$$

$${\Lambda }_{\mathrm{ff}}=1.426\times 10^{-27}{T}^{1/2}{g}_{\mathrm{ff}}{n}_{\mathrm{e}}{n}_{\mathrm{H}+},$$

$$Q_{\mathrm{net}}=Q_{\mathrm{XUV}}-({\Lambda }_{\mathrm{Ly}\alpha }+{\Lambda }_{\mathrm{rec}}+{\Lambda }_{\mathrm{col}}+{\Lambda }_{\mathrm{ff}}).$$

$${\Lambda }_{\mathrm{ad}}=\frac{1}{r^2}\frac{\partial E_{\mathrm{th}}u{r}^2}{\partial r},$$

$${\Lambda }_{\mathrm{ex}}=\frac{p}{r^2}\frac{\partial (u{r}^2)}{\partial r}$$

$$Q_{\mathrm{abs}}(r)=\sum _{\lambda }{\sigma }_{\mathrm{ion},\lambda }{n}_{\mathrm{H}}(r)F_{\mathrm{XUV},\lambda }\frac{\exp (-{\tau }_{\lambda }(r))}{1+\alpha {\tau }_{\lambda }(r)}.$$

$$R_{\mathrm{XUV}}={\left(\frac{4\pi {\int }_{r_0}^{r_1}{Q}_{\mathrm{XUV}}{r}^{2}dr}{\pi \bar{\eta }{F}_{\mathrm{XUV}}}\right)}^{1/2}$$

$${\dot{M}}_{\mathrm{el}}=\frac{\pi R_{\mathrm{p}}{R}_{\mathrm{XUV}}^2\bar{\eta }{F}_{\mathrm{XUV}}}{G{M}_{\mathrm{p}}K}$$

$$\frac{\partial (\rho \text{@vec}V)}{\partial t}$$

$$=Q_{\mathrm{i}}\text{@vec}{V}_{\mathrm{H}}-Q_{\mathrm{ex}}(\text{@vec}V-\text{@vec}{V}_{\mathrm{H}}),$$

$$\nabla \cdot \text{@vec}B$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial \text{@vec}B}{\partial t}$$

$$\frac{\partial W}{\partial t}$$

$$=\left(Q_{\mathrm{i}}+Q_{\mathrm{ex}}\right)\left(\frac{1}{2}{V}_{\mathrm{H}}^2+\frac{3k{T}_{\mathrm{H}}}{2m_{\mathrm{p}}}\right)-Q_{\mathrm{ex}}\left(\frac{1}{2}{V}^2+\frac{3kT}{2m_{\mathrm{p}}}\right),$$

$$W$$

$$Q_{\mathrm{i}}={\alpha }_{\mathrm{i}}{N}_{\mathrm{H}}{m}_{\mathrm{p}}$$

$$Q_{\mathrm{ex}}=\rho ⟨V_{\mathrm{rel}}⟩N_{\mathrm{H}}{\sigma }_{\mathrm{ex}}$$

$$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial r}=S,$$

$$\begin{gathered} U=\text{pmatrix}\rho r^2 \\ \rho u{r}^2 \\ {E}_{\mathrm{th}}{r}^2,F=\text{pmatrix}\rho u{r}^2 \\ \left(\rho u^2+p\right)r^2 \\ {E}_{\mathrm{th}}u{r}^2,S=\text{pmatrix}0 \\ g\rho r^2+2pr \\ {Q}_{\mathrm{net}}{r}^2-p\frac{\partial u{r}^2}{\partial r}. \end{gathered}$$

$$U_i^{n+1/2}=U_i-\Delta t\frac{F_{i+1}^n-F_i^n}{r_{i+1}-r_i}+\Delta t{S}_i^n,$$

Full text

*Aims. *For the hot Jupiter HD 189733b, this absorption shows temporal variability. We aim to study if variations in stellar XUV emission and/or variable stellar wind conditions may explain this effect.

*Methods. *We applied a 1D hydrodynamic planetary upper atmosphere model and a 3D magnetohydrodynamic stellar wind flow model to study the effect of variations of the stellar XUV irradiation and wind conditions at the planet’s orbit on the neutral hydrogen distribution. This includes the production of energetic neutral atoms (ENAs) and the related Ly$\alpha$ transit signature.

Results. We obtain comparable, albeit slightly higher Ly$\alpha$ absorption than that observed in 2011 with a stellar XUV flux of $1.8\times 10^{4}$ erg cm-2 s*-1*, rather typical activity conditions for this star. Flares with parameters similar to that observed eight hours before the transit are unlikely to have caused a significant modulation of the transit signature. We find that the resulting Ly$\alpha$ absorption is dominated by atmospheric broadening, whereas the contribution of ENAs is negligible, as they are formed inside the bow shock from decelerated wind ions that are heated to high temperatures. Thus, within our modeling framework and assumptions, we find an insignificant dependence of the absorption on the stellar wind parameters.

*Conclusions. *Since the transit absorption can be modeled with typical stellar XUV and wind conditions, it is possible that the nondetection of the absorption in 2010 was affected by less typical stellar activity conditions, such as a very different magnitude and/or shape of the star’s spectral XUV emission, or temporal and/or spatial variations in Ly$\alpha$ affecting the determination of the transit absorption.

1 Introduction

The first observational evidence of an expanding and escaping exoplanetary atmosphere was the Ly$\alpha$ transit absorption detected for the hot Jupiter HD 209458b (Vidal-Madjar et al., 2003). Subsequently, more observations of excess absorption in Ly$\alpha$ 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 make it 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\,M_{\mathrm{Jup}}$ and $1.13\,R_{\mathrm{Jup}}$, 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 $\sim$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 $\epsilon$ 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\pm 0.75$% in the unresolved Ly$\alpha$ line, significantly higher than the absorption by the planetary disk at optical wavelengths. Later, Lecavelier des Etangs et al. (2012) obtained Ly$\alpha$ 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\pm 3.6$% (i.e., $12.3\pm 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\pm 2.7$%, meaning $5.5\pm 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\times 10^{8}-10^{11}$ g s*-1*, an ionizing flux of 6-23 times the solar value, a stellar wind with a temperature of $3\times 10^{4}$ K, a velocity of $200\pm 20$ km s*-1*, and a density between $10^{3}-10^{7}$ 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 $\sim$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$\alpha$ 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 observations between 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$\alpha$ 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\pm 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$\alpha$ 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$\alpha$ 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 Section 2, we describe the hydrodynamic model and show the resulting upper atmosphere profiles and planetary mass-loss rates in Section 3. In Section 4, we describe the MHD flow model and show the stellar wind flow maps. In Section 5, we compute the Ly$\alpha$ absorption for the various scenarios and compare them with the observations. In Section 6, we discuss our findings and compare them with other studies. Conclusions are presented in Section 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,

[TABLE]

The gas parameters $\rho$, $u$, $T$, $E_{\mathrm{th}}=p/(\gamma-1)$ are the mass density, velocity, temperature, and thermal energy of the upper atmosphere. The gas pressure is $p=\rho RT/\mu$, with the mean molecular weight $\mu$ and the gas constant $R$. The distance from the planet’s center is $r,$ and $t$ is the time. For the specific heat ratio $\gamma,$ we adopted 5/3 for monatomic gas. The gravitational force $g=-\partial\Phi/\partial r$ was derived from the Roche potential $\Phi$ along the star-planet line:

[TABLE]

(e.g., García Muñoz, 2007). Here, $G$ is the gravitational constant, $M_{\mathrm{p}}$ 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 Section 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

[TABLE]

where $\alpha_{\mathrm{ion}}$ is the photoionization rate, $\alpha_{\mathrm{rec}}$ the radiative recombination rate, and $\alpha_{\mathrm{col}}$ the collisional ionization rate (all quantities in cgs units). The parameters $n_{\mathrm{H}},n_{\mathrm{H+}},n_{\mathrm{e}}$ are the number densities of neutral hydrogen atoms, protons, and electrons, respectively. We assumed quasi-neutrality, $n_{\mathrm{H+}}=n_{\mathrm{e}}$. The photoionization rate is given by

[TABLE]

where $\sigma_{\mathrm{ion,\lambda}}$ is the photoionization cross-section, $F_{\mathrm{XUV,\lambda}}$ the spectral XUV flux at the planet’s orbit outside of the atmosphere, $\tau_{\lambda}(r)$ the optical depth at distance $r$ from the planet’s center, and $E_{\lambda}$ 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 B111The 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. coefficient

[TABLE]

while the latter is given by

[TABLE]

where $\theta=\ln T_{\mathrm{e}}$ with the electron temperature $T_{\mathrm{e}}$ 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 $r_{0}$, 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

[TABLE]

where the term $1+\alpha\tau_{\lambda}(r)$ takes into account 2D effects of the energy absorption in an approximate way (Sekiya et al., 1980). The optical depth is calculated as

[TABLE]

We adopted a constant photoelectron heating efficiency $\eta_{\mathrm{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 $\alpha=4$. We reverted to a 1D calculation of $Q_{\mathrm{XUV}}$, 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 $Q_{\mathrm{XUV}}$ is very similar to the 1D method shown in Eq. 13 (see Appendix B for a comparison). The adopted XUV spectra are described in Section 2.4.

We included several cooling processes, namely Ly$\alpha$ 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$\alpha$ cooling rate

[TABLE]

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

[TABLE]

by collisional ionization

[TABLE]

and by free-free emission

[TABLE]

with $g_{\mathrm{ff}}=0.79464+0.1243\log_{10}(T)$. The net heating rate is then computed via

[TABLE]

Thermal conduction is also included in the model. The conductivity coefficient is calculated as $\kappa=n^{-1}\sum_{j}n_{j}\kappa_{j}$ (e.g. García Muñoz, 2007), where $n_{j}$ are the number densities of the constituents and $\kappa_{j}=A_{j}T^{s_{j}}$, where $A_{j}$ and $s_{j}$ are the fitting parameters for the individual conductivities. We adopted $A_{\mathrm{H}}=379$, $A_{\mathrm{H+}}=7.37\times 10^{-8}$, $A_{\mathrm{e}}=1.2\times 10^{-6}$, $s_{\mathrm{H}}=0.69$, and $s_{\mathrm{H+}}=s_{\mathrm{e}}=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

[TABLE]

and expansion

[TABLE]

(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 on this specific planet.

2.3 Initial and boundary conditions

The simulation domain extends from the planet’s optical transit radius $r_{0}=R_{\mathrm{p}}$ to $r_{1}=4.5R_{\mathrm{p}}$. The upper boundary lies near the L1 point $R_{\mathrm{L1}}=(\mu-\mu^{2}/3)a\sim 4.3R_{\mathrm{p}}$, where $\mu=(M_{\mathrm{p}}/(3M_{*}))^{1/3}$ (e.g., Erkaev et al., 2007). We note that due to the high gravity of HD 189733b, the distance between $R_{\mathrm{p}}$ and typical mesopause pressure levels is very small. We used a nonuniform grid $r_{i}=r_{1}^{i/(N-1)}$ 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 $T_{0}$ and number density $n_{0}$. For the temperature, we adopted the equilibrium temperature of 1200 K (cf., Section 2.4) and for the number density a value of $10^{15}\,\mathrm{cm}^{-3}$. This corresponds to a pressure of 0.16 mbar. Changing $T_{0}$ to 800 K or 1600 K decreases or increases, respectively, the mass-loss rate by only $\sim$10%, and affects the atmospheric profiles negligibly (cf., Section 3.1). We checked that the assumed density $n_{0}$ is chosen large enough 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 $n_{0}$. Smaller values of $n_{0}$ would lead to significant ionization below the optical transit radius, which would be inconsistent with its measurement (see Section 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(\Phi(r)-\Phi(r_{0}))$, but modified it in the upper region to $n(r)\propto 1/r^{2}$. 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 $n_{\mathrm{H+}}(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 $\rho ur^{2}$ 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. 16). 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.

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 et al., 2016). The distance $d$, stellar mass $M_{*}$ and radius $R_{*}$, stellar bolometric flux at Earth $F_{\mathrm{bol}}$, planetary mass $M_{\mathrm{p}}$ and radius $R_{\mathrm{p}}$, as well as the orbital inclination $i$ were taken from their study. The orbital separation $a$ was calculated from their quoted value of the parameter $a/R_{*}=8.84\pm 0.27$ and the stellar radius. The average orbital velocity is therefore 151 km s*-1*, consistent with the maximum orbital radial velocity of $154^{+4}_{-3}$ km s*-1* measured by de Kok et al. (2013). Using the measured bolometric stellar flux, we find an equilibrium temperature $T_{\mathrm{eq}}=\left(F_{\mathrm{bol}}(1-A)/(f\sigma)\right)^{1/4}$ of 1200 K, assuming full redistribution ($f=4$) and zero albedo ($A=0$). Here, $\sigma$ is the Stefan-Boltzmann constant. This is similar to, but slightly higher than, the apparent effective dayside temperature of $1163\pm 37$ K (Schwartz et al., 2017). Such small differences in the adopted $T_{\mathrm{eq}}$ 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\times 10^{28}$ 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$\alpha$ 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$\alpha$ 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$\alpha$ flux at Earth amounts to $7.5\times 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\times 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\times 10^{3}$ erg cm*-2* s*-1* with a spectral energy distribution shown in Fig. 1.

Secondly, we examined the synthetic spectrum from the X-exoplanets222http://sdc.cab.inta-csic.es/xexoplanets/jsp/homepage.jsp database (Sanz-Forcada et al., 2010, 2011, hereafter SF11). These authors used observed X-ray and UV spectra, reconstructed the emission measure 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 $\lambda{>}300$ Å range is consistent within the estimated uncertainties. The $\lambda{<}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$\alpha$ 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$\alpha$ observations. Simultaneous observations with Swift/XRT during the 2011 transit yielded an X-ray (0.3-3 keV) flux of $3.6\times 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\times 10^{-13}$ erg cm*-2* s*-1* (SF11), and $3.2{-}3.9\times 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\times 10^{-13}$ erg cm*-2* s*-1* (Pillitteri et al., 2014), which amounts to $6.2\times 10^{3}$ 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\times 10^{3}$ erg cm*-2* s*-1* for the L14 method and $1.2\times 10^{4}$ 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\times 10^{4}$ 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\times 10^{4}$ and $1.8\times 10^{4}$ erg cm*-2* s*-1* for L14 and SF11, respectively (cf., Table 2). We explore the influence of XUV flux variations further in Section 3.2, where we model the effects of a flare.

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 $\dot{M}=4\pi r^{2}\rho u$ (assuming that the star-planet line value is representative for entire atmosphere) are $5.4\times 10^{10}$ and $2.5\times 10^{10}$ 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 ($n_{\mathrm{H+}}/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.85R_{\mathrm{p}}$ 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.3R_{\mathrm{p}}$ in both cases, whereas it already exceeds the escape velocity at points that are ${\sim}0.3R_{\mathrm{p}}$ closer to the planet. The temperature reaches a maximum of ${\sim}1.1\times 10^{4}$ K at ${\sim}2R_{\mathrm{p}}$ 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=\Lambda/X=1$, where $\Lambda=1/(n\sigma_{\mathrm{col}})$ is the mean free path and $X$ is an appropriate system scale, for example, the scale height $H=kT/(m_{\mathrm{p}}g)$ for regions close to the planet, or the planetary radius for more distant regions (Shaikhislamov et al., 2014). For the collision cross-section $\sigma_{\mathrm{col}}$, we adopted the H-H+ charge exchange cross-section (${\sim}2\times 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 $R_{\mathrm{p}}$, we find that $Kn{\ll}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.5R_{\mathrm{p}}$. Radiative cooling is also very important for HD 189733b and is dominated by Ly$\alpha$ 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 $\eta_{\mathrm{XUV}}(r)$ defined as the ratio of the XUV volume heating rate $Q_{\mathrm{XUV}}$ (Eq. 13) to the locally absorbed XUV radiation

[TABLE]

Secondly, we calculated $\eta_{\mathrm{net}}$ using the net local heating rate $Q_{\mathrm{net}}(r)$ (Eq. 19) divided by $Q_{\mathrm{abs}}(r)$ (e.g., Salz et al., 2016a). A comparison is shown in Fig. 4. The corresponding mean heating efficiencies ($\bar{\eta}_{\mathrm{XUV}}$, $\bar{\eta}_{\mathrm{net}}$) in the atmosphere, obtained by integrating $Q_{\mathrm{XUV}}(r)$ and $Q_{\mathrm{net}}(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 $\eta_{\mathrm{ph}}$ described in Section 2.2.

The effective XUV absorption radius

[TABLE]

(Erkaev et al., 2007, 2015) is about $2.6R_{\mathrm{p}}$ for SF11 and $2.1R_{\mathrm{p}}$ for L14. Using $Q_{\mathrm{net}}$ instead of $Q_{\mathrm{XUV}}$ in Eq. 23 yields 1.15 (SF11) and $0.78R_{\mathrm{p}}$ (L14), respectively. We note that this radius can be smaller than 1 if $Q_{\mathrm{net}}$ is used in Eq. 23 and radiative cooling is very efficient, and/or a too large value is chosen for the mean heating efficiency $\bar{\eta}$ in the denominator. Evaluating the energy-limited mass-loss rate

[TABLE]

(Erkaev et al., 2007) with these parameters yields $1.5{-}3.1\times 10^{11}$ g s*-1* if adopting a typical mean heating efficiency for hot Jupiters of $\bar{\eta}=15$% (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 $R_{\mathrm{XUV}}$ values obtained with $Q_{\mathrm{net}}$ yields $2.1{-}6.2\times 10^{10}$ g s*-1*, in much better agreement with the hydrodynamic rates. We note that the simple assumption of $R_{\mathrm{XUV}}{\sim}R_{\mathrm{p}}$ would coincidentally yield, for this specific planet, mass-loss rates of $3.4{-}4.7\times 10^{10}$ g s*-1*, which are similar to the hydrodynamic results.

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 $\sim$14% in Ly$\alpha$ 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 $\sim$27 min, which corresponds to one bin of the temporal resolution of the observations. Together with the average pre-flare flux of $3.6\times 10^{-13}$ erg cm*-2* s*-1*, this yields an estimated energy of ${\sim}8\times 10^{31}$ 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\times 10^{10}$ (X-ray flare) and $1.2\times 10^{11}$ 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.95R_{\mathrm{p}}$ (the particle launch radius in their model) are $2.2\times 10^{9}$ (X-ray flare) and $8.9\times 10^{8}$ 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 ($10^{9}$ 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\times 10^{8}{-}1.5\times 10^{9}$ 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., Section 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 Section 4.3. We note that other upper atmosphere models of this planet find higher ionization (and therefore lower H densities) at $2.95R_{\mathrm{p}}$, 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

[TABLE]

where $\rho$, $\@vec{V}$, $P$, and $\@vec{B}$ are the mass density, velocity, plasma pressure, and magnetic field of the stellar wind, respectively. The parameter $\gamma$ is the polytropic index (assumed to be equal to 5/3) and $\@vec{I}$ is the identity matrix. The parameters $\@vec{V}_{\mathrm{H}}$ and $T_{\mathrm{H}}$ 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

[TABLE]

and charge exchange ionization

[TABLE]

of the hydrogen atoms. Here, $N_{\mathrm{H}}$ is the number density of the neutral planetary hydrogen atoms, $m_{\mathrm{p}}$ the particle mass, $\sigma_{\mathrm{ex}}$ ($\sim$10*-15* cm2) the charge exchange cross-section, $\langle V_{\mathrm{rel}}\rangle$ the average relative speed of the stellar wind and atmospheric particles, and $\alpha_{\mathrm{i}}=5.9\times 10^{-8}F_{\mathrm{XUV}}$ is the ionization rate proportional to the XUV flux $F_{\mathrm{XUV}}=1.8\times 10^{4}$ 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 ($R_{\mathrm{s}}$) 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 undisturbed stellar 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.

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 $N_{\mathrm{sw}}$, velocity $V_{\mathrm{sw}}$, temperature $T_{\mathrm{sw}}$, and magnetic field strength along the orbit. We computed both the parallel $B_{\mathrm{sw,p}}$ (parallel to $V_{\mathrm{sw}}$) and normal $B_{\mathrm{sw,n}}$ (normal to $V_{\mathrm{sw}}$) components of its magnitude, as well as the total field strength $B_{\mathrm{sw,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 $\theta_{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. 19). 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.6R_{\mathrm{p}}$. Here, the ion and neutral densities of the atmospheric particles are $3.9\times 10^{7}$ cm*-3* and $1.25\times 10^{6}$ cm*-3*, respectively, and the temperature is $8.8\times 10^{3}$ K. In the maximum wind case, the magnetopause position is slightly closer, at $2.4R_{\mathrm{p}}$, since the stellar wind dynamic pressure is larger. The corresponding atmospheric ion and neutral densities are $6.1\times 10^{7}$ cm*-3* and $2.9\times 10^{6}$ cm*-3*, respectively, and the temperature is $1\times 10^{4}$ K (cf. Fig. 2).

Figure 6 shows the spatial distribution of the total pressure $P_{\mathrm{tot}}$ (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 $P_{\mathrm{dsw}}$ = $N_{\mathrm{sw}}m_{\mathrm{p}}V_{\mathrm{sw}}^{2}=7.8\times 10^{-5}$ 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 $\@vec{B}_{\mathrm{sw}}$ (the arrow in Fig. 6). The white area around the center of the coordinate system indicates the region filled by atmospheric particles exclusively, while the embedded dark blue semicircle indicates the optical planetary radius $R_{\mathrm{p}}$. 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 $P_{\mathrm{dsw}}$ = $1.8\times 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), both corresponding to the mean stellar wind. The total ion number density has a maximum of about $4\times 10^{4}$ 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\times 10^{7}$ cm*-3* and $2.3\times 10^{7}$ cm*-3*, respectively. The ENAs form a layer around the magnetopause with a thickness of about $1.3R_{\mathrm{p}}$. 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\times 10^{7}$ cm*-3* and $5.6\times 10^{7}$ cm*-3*, respectively. The ENA layer thickness around the magnetopause is about $1.2R_{\mathrm{p}}$.

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$\alpha$ 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 ${\sim}8\times 10^{31}$ erg (cf., Section 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 $10^{16}$ 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 $L_{\mathrm{X}}$ (Odert et al., 2017). Adopting $L_{\mathrm{X}}=1.67\times 10^{28}$ erg s*-1* (cf., Section 2.4), we obtain about $10{-}2000$ CMEs per day, depending on the power law index of the stellar flare energy distribution, $dN/dE{\propto}E^{-\alpha}$. We assumed $\alpha=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$\alpha$ 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 of the spectra (Odert et al., 2020). It predicts that $\lesssim$1% of the intrinsically occurring CMEs would be observable in H$\alpha$ 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$\alpha$ observations were 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 $10R_{\sun}$ from the Sun (Patsourakos & Georgoulis, 2016). The angle $\theta_{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 ${\sim}2R_{\mathrm{p}}$. The corresponding atmospheric ion and neutral densities are about $1.8\times 10^{8}$ cm*-3* and $4.5\times 10^{7}$ cm*-3*, respectively, and the temperature is $1.16\times 10^{4}$ K (cf. Fig. 2). The total pressure is again normalized to the CME dynamic pressure, $P_{\mathrm{dsw}}=7.35\times 10^{-3}$ Pa. The ion density and ENA density maxima in front of the magnetopause (Fig. 11) are about $0.67\times 10^{7}$ cm*-3* and $1.5\times 10^{8}$ cm*-3*, respectively. The ENA layer is much thinner in the CME case, with a thickness of about $0.25R_{\mathrm{p}}$.

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.

5 Modeling the Ly$\alpha$ absorption signature

To calculate the transmissivity in the Ly$\alpha$ 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, 10) and through knowledge of the positions and velocities of all hydrogen particles, we computed how these atoms attenuate the stellar Ly$\alpha$ radiation by using 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$\alpha$ 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 $\Delta f{\sim}\Delta E/\hbar{\sim}1/\Delta 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 $I_{0}$ as a function of frequency $f$ of the stellar spectrum in the $yz$-plane by dividing the computational domain into a grid with $N_{\mathrm{c}}$ cells. For each cell in the grid along lines of sight in front of the star ($y^{2}+z^{2}<R_{*}^{2}$), 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=a\cos{i}=3.67\times 10^{10}$ 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$\alpha$ 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 $(y-y_{\mathrm{p}})^{2}+(z-z_{\mathrm{p}})^{2}<R_{\mathrm{p}}^{2}$, where ($y_{\mathrm{p}}$,$z_{\mathrm{p}}$) 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 atmospheric profile 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 all three 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=f_{0}+v_{x}/\lambda_{0}$ with $f_{0}=c/\lambda_{0}$, $\lambda_{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$\alpha$ 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\ldots-40$ km s*-1*) absorption is 10.9%, 1.7$\sigma$ larger than the observed value ($6.8\pm 2.4$%); for the red wing ($+40\ldots+400$ km s*-1*), we obtain, with 12.9%, too much absorption by 5.5$\sigma$ (observed: $4.5\pm 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$\sigma$ of the observations. The overabsorption could be due to limitations of our 1D hydrodynamic atmosphere model, as we used the atmospheric profiles along the star-planet line as the global input to the 3D MHD model. Recent 3D hydrodynamic models obtain more radially 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$\alpha$ 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$\alpha$ signatures. Differences to previous studies are discussed in more detail in Section 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.

6 Discussion

6.1 Comparison with other hydrodynamic models

As described above, our simulations yield a total mass-loss rate of $\dot{M}=2.5{-}5.4\times 10^{10}$ g s*-1*, depending on the adopted XUV spectrum. Previous studies of this planet found $\dot{M}=4.8\times 10^{10}{-}2\times 10^{11}$ g s*-1* for $F_{\mathrm{XUV}}=2\times 10^{4}{-}10^{5}$ erg cm*-2* s*-1* (Guo, 2011), $\dot{M}=4.5\times 10^{11}{-}9\times 10^{11}$ g s*-1* for $F_{\mathrm{XUV}}=24778$ erg cm*-2* s*-1* (but different spectral energy distributions; Guo & Ben-Jaffel, 2016), and $\dot{M}=1.64\times 10^{10}$ g s*-1* for $F_{\mathrm{XUV}}=20893$ erg cm*-2* s*-1* (Salz et al., 2016a)333Their given value of $4.1\times 10^{9}$ g s*-1* corresponds to 1/4 of the isotropic mass-loss rate.. 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 to the 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 $\beta=F_{50-400\AA}/F_{50-900\AA}{\sim}F_{0-400\AA}/F_{0-912\AA}{\sim}0.8$ for both spectra, corresponding to a mass-loss rate of ${\sim}9\times 10^{11}$ 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 $F_{\mathrm{XUV}}=2\times 10^{4}$ 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 $n_{0}$ 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 $2R_{\mathrm{p}}$ instead of $1.5R_{\mathrm{p}}$.

The main differences between the discussed models are: gray atmosphere approximation, no molecules, only H, only Ly$\alpha$ cooling in Guo (2011); XUV spectra, molecules, H+He, but only H${}_{3}^{+}$ 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 ${\sim}3.1R_{\mathrm{p}}$ (Eggleton, 1983), smaller than the distance to the L1 point (cf., Section 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\,R_{\mathrm{p}}$ (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\,R_{\mathrm{p}}$, 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$\alpha$ 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${}_{3}^{+}$ 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.

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$\alpha$ 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 ($V_{\mathrm{sw}}$, $N_{\mathrm{sw}}$, $T_{\mathrm{sw}}$). 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 $V_{\mathrm{sw}}=200$ km s*-1*, $T_{\mathrm{sw}}=3\times 10^{4}$ K, $N_{\mathrm{sw}}=10^{3}{-}3\times 10^{5}$ cm*-3*, an ionizing flux ten times the solar one, and neutral H escape rates of $5\times 10^{8}{-}1.5\times 10^{9}$ 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 $\sim$8 h prior to the transit in 2011.

The model of Bourrier & Lecavelier des Etangs (2013) can explain the observed Ly$\alpha$ absorption in 2011 if the neutral H mass-loss rate at 2.95 $R_{\mathrm{p}}$ (the lower boundary of their model) is in the order of $10^{9}$ 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 Section 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 $R_{\mathrm{p}}$, 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 lower number densities assumed at the lower boundary in these studies. If we run our model with lower $n_{0}$ values, we also obtain neutral densities at 2.95 $R_{\mathrm{p}}$ 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 $n_{0}$, stronger than the total mass-loss rate. Higher values of $n_{0}$ shift the neutral-ion transition to larger heights. However, we find that lowering $n_{0}$ results in incomplete absorption of the XUV radiation in the computational domain, which makes the resulting mass-loss rates strongly dependent on the assumed $n_{0}$.

One may argue that if the true value of $n_{0}$ were indeed much lower than our assumptions, the absorption signal from the planetary atmosphere could be much weaker than what our results show and not compatible with the 2011 observations. However, it is unlikely that $n_{0}$ 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 $R_{\mathrm{p}}$ 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 (Section 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 $V_{\mathrm{sw}}=200$ km s*-1*, $T_{\mathrm{sw}}=3\times 10^{4}$ K, $N_{\mathrm{sw}}=10^{6}$ 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 collisions444Taking $n{\gtrsim}N_{\mathrm{sw}}$, $\sigma_{\mathrm{col}}{\sim}10^{-15}$ cm2 and $X{\sim}R_{\mathrm{p}}$ yields a Knudsen number $Kn{<}1$ in the modeling domain.. 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 Section 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 XUV radiation. They are able to reproduce the 2011 observations with a spectral shape $\beta=0.76$, which is very similar to our adopted SF11 spectrum ($\beta=0.82$) that also reproduces the transit absorption. The nondetection in 2010 would require a spectral shape of $\beta=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 Section 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$\alpha$ transit absorption

Stellar activity can also affect the study of Ly$\alpha$ 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$\alpha$ can only be studied in the line wings due to geocoronal emission and ISM absorption. The wings of the Ly$\alpha$ 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$\alpha$ 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$\alpha$ 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 ($F_{\mathrm{XUV}}\sim 1.8\times 10^{4}$ erg cm*-2* s*-1*). We note that with our adopted default parameters, the model produces some extra-absorption across the Ly$\alpha$ 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$\alpha$ 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-$\alpha$ 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” (KK, HL, NVE), as well as FWF projects P30949-N36 (PO, ML) and I2939-N27. LF, DK and NVE acknowledge also the Austrian Forschungsförderungsgesellschaft FFG project “TAPAS4CHEOPS” P853993. NVE and VAI 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-gravity planet 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

[TABLE]

where $U$ is the vector of the variables, $F$ the fluxes, and $S$ the source terms. Here,

[TABLE]

The MacCormack method consists of a predictor step:

[TABLE]

and a corrector step:

[TABLE]

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 $\Delta t=C_{\mathrm{CFL}}\Delta x_{\mathrm{min}}/(|u|+c)_{\mathrm{max}}$, where $\Delta x_{\mathrm{min}}$ is the smallest grid size, $c$ is the adiabatic sound speed, $(|u|+c)_{\mathrm{max}}$ is the maximum wave speed on the grid, and $C_{\mathrm{CFL}}\leq 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 $u_{i}{<}0$ and $u_{i+1}{>}0$. Since the momentum flux $\rho u^{2}$ loses the information on the sign of $u$, MacCormack (1971) suggested, in such cases, replacing the momentum flux differences $(\rho u^{2})_{i+1}-(\rho u^{2})_{i}$ with $[(\rho u)_{i+1}-(\rho u)_{i}](u_{i+1}+u_{i})/2$ in Eq. 35, and with correspondingly adjusted indices in Eq. 36. 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,

[TABLE]

where $\Delta U_{i+1/2}^{n}=U_{i+1}^{n}-U_{i}^{n}$, $\Delta U_{i-1/2}^{n}=U_{i}^{n}-U_{i-1}^{n}$, and $\Delta r_{i}=(r_{i+1}-r_{i-1})/2$. The coefficients $\hat{K}$ are determined by

[TABLE]

where $|\lambda_{i\pm 1/2}|_{\mathrm{max}}=|u_{i\pm 1/2}|+c_{i\pm 1/2}$ is the maximum local wave speed and $\psi(M)=\max(0,\min(1,M))$ is the minmod limiter. The coefficients $M^{\pm}$ are calculated via

[TABLE]

where $s$ denotes the sign function of the superscript555Meaning $s=1$ for superscript “$+$” and $s=-1$ for superscript “$-$”. and the expression $(*\cdot*)$ refers to the inner product of the difference vectors $\Delta U$. The parameter

[TABLE]

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. 34) are modified to $F_{2}=\rho u^{2}r^{2}$ and $S_{2}=-\rho r^{2}\partial\Phi/\partial r-r^{2}\partial p/\partial r$. In the limit of $u\rightarrow 0$, hydrostatic equilibrium should be reached, and the two terms in $S_{2}$ 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 $\Phi$ and $p$ in $S_{2}$ are discretized the same way as the respective advection terms in Eqs. 35 and 36; secondly, the prefactors of these derivatives ($\rho r^{2}$, $r^{2}$) are calculated as averages of their values at the grid points involved in the derivatives. For Eq. 35, this yields

[TABLE]

and similar for Eq. 36 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 16 shows the mass flux $\rho ur^{2}$ 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.

Appendix B Comparison of 1D and 2D heating

Figure 17 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 is distributed 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, $\alpha=4$ represents the best match and was consequently used in all runs.

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 $R_{\mathrm{p}}$ where velocities are very small. The H2 density can be written as

[TABLE]

where $\gamma_{\mathrm{H}}=8\times 10^{-33}(300/T)^{0.6}$ is the rate of the reaction producing H2 (H+H$\rightarrow$H2) and $\nu_{\mathrm{diss}}=1.5\times 10^{-9}\exp(-49000/T)$ the thermal dissociation rate (Yelle 2004; Erkaev et al. 2016). The photoionization rate $\nu_{\mathrm{H_{2}}}$ is calculated after Eq. 6, but replacing $\sigma_{\mathrm{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. 18. One can see that H2 would only be the dominating species at distances ${\lesssim}1.04R_{\mathrm{p}}$, 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.

Appendix D Stellar wind parameters

In Fig. 19, 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).