<figure><img src="image_1.png" /></figure>

HAL
open science

# Colliding planetary and stellar winds: charge exchange and transit spectroscopy in neutral hydrogen

Pascal Tremblin, Eugene Chiang

► To cite this version:

Pascal Tremblin, Eugene Chiang. Colliding planetary and stellar winds: charge exchange and transit spectroscopy in neutral hydrogen. Monthly Notices of the Royal Astronomical Society, 2013, 428, pp.2565-2576. 10.1093/mnras/sts212. insu-03621069

HAL Id: insu-03621069

https://insu.hal.science/insu-03621069

Submitted on 27 Mar 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Colliding planetary and stellar winds: charge exchange and transit spectroscopy in neutral hydrogen

Pascal Tremblin<fnref n="1" />* and Eugene Chiang<fnref n="2" />

Accepted 2012 October 15. Received 2012 October 9; in original form 2012 June 19

## ABSTRACT

When transiting their host stars, hot Jupiters absorb about 10 per cent of the light in the wings of the stellar Lyman α emission line. The absorption occurs at wavelengths Doppler-shifted from line centre by ±100 km s⁻<fnref n="1" /> – larger than the thermal speeds with which partially neutral, ~10⁴ K hydrogen escapes from hot Jupiter atmospheres. It has been proposed that the absorption arises from ~10⁶ K hydrogen from the host stellar wind, made momentarily neutral by charge exchange with planetary H I. The ±100 km s⁻<fnref n="1" /> velocities would then be attributed to the typical velocity dispersions of protons in the stellar wind – as inferred from spacecraft measurements of the solar wind. To test this proposal, we perform 2D hydrodynamic simulations of colliding hot Jupiter and stellar winds, augmented by a chemistry module to compute the amount of hot neutral hydrogen produced by charge exchange. We observe the contact discontinuity where the two winds meet to be Kelvin–Helmholtz unstable. The Kelvin–Helmholtz instability mixes the two winds; in the mixing layer, charge exchange reactions establish, within tens of seconds, a chemical equilibrium in which the neutral fraction of hot stellar hydrogen equals the neutral fraction of cold planetary hydrogen (about 20 per cent). In our simulations, enough hot neutral hydrogen is generated to reproduce the transit observations, and the amount of absorption converges with both spatial resolution and time. Our calculations support the idea that charge transfer between colliding winds correctly explains the Lyman α transit observations – modulo the effects of magnetic fields, which we do not model but which may suppress mixing. Other neglected effects include, in order of decreasing importance, rotational forces related to orbital motion, gravity and stellar radiation pressure; we discuss quantitatively the errors introduced by our approximations. How hot stellar hydrogen cools when it collides with cold planetary hydrogen is also considered; a more careful treatment of how the mixing layer thermally equilibrates might explain the recent detection of Balmer Hα absorption in transiting hot Jupiters.

**Key words:** line: formation – methods: numerical – planets and satellites: atmospheres – stars: winds, outflows – ultraviolet: planetary systems.

## 1 INTRODUCTION

Gas-laden planets lose mass to space when their upper atmospheres are heated by stellar ultraviolet (UV) radiation. Ubiquitous in the Solar system, thermally driven outflows modify the compositions of their underlying atmospheres over geologic time (e.g. Weissman, McFadden & Johnson 1999). Thanks to the *Hubble Space Telescope* (*HST*), escaping winds are now observed from extrasolar hot Jupiters: Jovian-sized planets orbiting at distances ≲0.05 au from

their host stars and bathed in intense ionizing fields. Spectroscopy with *HST* reveals absorption depths of ~2–10 per cent in various resonance transitions (H I, O I, C II, Si III and Mg II) when the planet transits the star, implying gas outflows that extend for at least several planetary radii (e.g. Vidal-Madjar et al. 2003, 2004, 2008, VM03; Ben-Jaffel 2007, 2008; Fossati et al. 2010; Lecavelier Des Etangs et al. 2010; Linsky et al. 2010). Recent observations of HD 189733b also indicate temporal variations in H I Lyman α absorption, possibly correlated with stellar X-ray activity (Lecavelier des Etangs et al. 2012). These data promise to constrain the compositions of hot Jupiter atmospheres and the degrees to which they are vertically mixed (Liang et al. 2003; Moses et al. 2011).

<footnote n="1">¹Laboratoire AIM Paris-Saclay (CEA/Irfu - Uni. Paris Diderot - CNRS/INSU), Centre d'études de Saclay, F-91191 Gif-Sur-Yvette, France</footnote>

<footnote n="2">²Departments of Astronomy and of Earth and Planetary Science, University of California at Berkeley, Hearst Field Annex B-20, Berkeley, CA 94720-3411, USA</footnote>

<footnote>* E-mail: pascal.tremblin@cea.fr</footnote>
The HST observations of hot Jupiter winds are accompanied by theoretical studies that model planetary outflows starting from first principles (e.g. Yelle 2004, 2006; Tian et al. 2005; García Muñoz 2007; Murray-Clay, Chiang & Murray 2009, M09). These 1D hydrodynamic models generally agree that hot Jupiters like HD 209458b and HD 189733b are emitting $\dot{M} \sim 10^{10}-10^{11}$ g s$^{-1}$ in mostly hydrogen gas. Three-dimensional (3D) models include Lecavelier des Etangs et al. (2004) and Jaritz et al. (2005), who emphasize the importance of tidal forces.

Do the models agree with the observations? Linsky et al. (2010) find that their observations of C II absorption in HD 209458b can be made consistent with modelled mass-loss rates, assuming the carbon abundance of the wind is not too different from solar. More comparisons between observation and theory would be welcome – particularly for hydrogen, the dominant component of the wind. But the observations of H I absorption have proven surprisingly difficult to interpret. On the one hand, the original measurements by VM03 indicate substantial (~10 per cent) absorption at Doppler shifts of $\pm 100$ km s$^{-1}$ from the centre of the H I Lyman $\alpha$ line. On the other hand, theory (e.g. M09) indicates that planetary outflows, heated by photoionization to temperatures $T \lesssim 10^4$ K, blow only at ~10 km s$^{-1}$. How can such slow planetary winds produce significant absorption at $\pm 100$ km s$^{-1}$? Measurements of blueshifted velocities as large as -230 km s$^{-1}$ in the case of HD 189733b only accentuate this problem (Lecavelier des Etangs et al. 2012).

Holmström et al. (2008, hereafter H08) propose that the observed energetic neutral H atoms arise from charge exchange between planetary H I and protons from the incident stellar wind. In this interpretation, the $\pm 100$ km s$^{-1}$ velocities correspond to the thermal velocities of $10^6$ K hydrogen from the star–hydrogen which is made neutral by electron exchange with planetary H I. The situation is analogous to that of the colliding winds of O star binaries (Stevens, Blondin & Pollock 1992; Lamberts, Fromang & Dubus 2011, and references therein). The H I Lyman $\alpha$ absorption arises from the contact discontinuity where the two winds meet, mix and charge exchange to produce hot neutral hydrogen.

The calculations of H08, and those of the follow-up study by Ekenbäck et al. (2010, hereafter E10), are based on a Monte Carlo algorithm that tracks individual ‘meta-particles’ of neutral hydrogen launched from the planet. The meta-particles collide and charge exchange with stellar wind protons outside a presumed planetary magnetosphere, which is modelled as an ‘obstacle’ in the shape of a bow shock. Good agreement with the Ly$\alpha$ observations is obtained for a range of stellar and planetary wind parameters, and for a range of assumed obstacle sizes.

In this work we further test the hypothesis of charge exchange first explored by H08 and E10. Our methods are complementary: instead of adopting their kinetic approach, we solve the hydrodynamic equations. We do not prescribe any obstacle to deflect the stellar wind, but instead allow the planetary and stellar winds to meet and shape each other self-consistently via their ram and thermal pressures. Some aspects of our solution are not realistic – we ignore the Coriolis force, the centrifugal force, stellar tidal gravity and magnetic fields.<fnref n="1" /> Our goal is to develop a first-cut hydrodynamic-chemical model of the contact discontinuity between the two winds where material mixes and charge exchanges. Simple and physically motivated scaling relations will be developed between the amount

of H I absorption and the properties of the stellar and planetary winds.

The plan of this paper is as follows. In Section 2 we describe our numerical methods, which involve augmenting our grid-based hydrodynamics code to solve the chemical reactions of charge exchange, and specifying special boundary conditions to launch the two winds. In Section 3 we present our results, including a direct comparison with the H I Ly$\alpha$ transit spectra of VM03, and a parameter study to elucidate how the absorption depth varies with stellar and planetary wind properties. A summary is given in Section 4, together with an assessment of the shortcomings of our study and pointers towards future work.

## 2 NUMERICAL METHODS

In Section 2.1, we describe the hydrodynamics code used to simulate the colliding planetary and stellar winds. In Section 2.2 we detail the charge exchange reactions that were added to the code. In Section 2.3, we outline our post-processing procedure for computing the Lyman $\alpha$ transmission spectrum. As a convenience to readers, in Section 2.4 we re-cap the differences between our treatment of colliding winds and that of E10/H08.

### 2.1 Hydrodynamics: code and initial conditions

Our simulations are performed with HERACLES (González, Audit & Huynh 2007),<fnref n="2" /> a grid-based code using a second-order Godunov scheme to solve the Euler equations:

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{V}) &= 0 \\
\frac{\partial \rho \mathbf{V}}{\partial t} + \nabla \cdot [\rho \mathbf{V} \otimes \mathbf{V} + p \mathbf{I}] &= 0 \\
\frac{\partial E}{\partial t} + \nabla \cdot [(E + p)\mathbf{V}] &= 0 \\
\frac{\partial \rho x_i}{\partial t} + \nabla \cdot (\rho x_i \mathbf{V}) &= 0.
\end{aligned}
\quad (1) $$

Here $\rho$, $\mathbf{V}$, $p$ and $E$ are the mass density, velocity, pressure and total energy density, respectively (e.g. Clarke & Carswell 2003). The code tracks abundances of individual species: $x_i$ is the mass fraction of the $i$th species of hydrogen, where $i \in \{1, 2, 3, 4\}$ to cover four possible combinations of ionization state (either neutral or ionized) and temperature (either ‘hot’ because it arises from the star or ‘cold’ because it arises from the planet). The outer product is denoted as $\otimes$, and $\mathbf{I}$ is the identity matrix.

All our simulations are 2D Cartesian in the dimensions $x$ (stellar radius) and $y$ (height above the planet’s orbital plane). Equivalently the simulations may be regarded as 3D, but with no rotation and with a star and a planet that are infinite cylinders oriented parallel to the $z$-axis. At fixed computational cost, 2D simulations enjoy better spatial resolution than 3D simulations and thus better resolve the fluid instabilities at the interface of the two winds. The standard box size is $(L_x, L_y) = (40R_p, 60R_p)$, where $R_p = 10^{10}$ cm is the planet radius. The number of grid points ranges up to $(N_x, N_y) = (6400, 9600)$ (see Table 1).

The star and its wind are modelled after the Sun and the solar wind. The stellar wind is injected through the left edge of the simulation box; the densities, velocities and temperatures in the vertical column of cells at the box’s left edge are fixed in time. Stellar wind

<footnote n="1">¹ For recent explorations of star–planet interactions including magnetic forces, see Cohen et al. (2011a,b). These simulations do not resolve the mixing layer interface between the stellar and planetary winds.</footnote>

<footnote n="2">² http://irfu.cea.fr/Projets/Site_heracles/index.html</footnote>
<table><caption>**Table 1.** Parameters of the winds at launch, and of the simulation box.</caption><thead><tr><th>Stellar wind</th><th>Planetary wind</th></tr></thead><tbody><tr><td>r<sub>launch,*</sub> = 5R<sub>☉</sub></td><td>d<sub>launch,p</sub> = 4R<sub>p</sub></td></tr><tr><td>n<sub>*</sub> = 2.9E4 cm<sup>-3</sup></td><td>n<sub>p</sub> = 3.9E6/cm<sup>3</sup></td></tr><tr><td>T<sub>*</sub> = 1E6 K</td><td>T<sub>p</sub> = 7000 K</td></tr><tr><td>v<sub>*</sub> = 130 km s<sup>-1</sup></td><td>v<sub>p</sub> = 12 km s<sup>-1</sup></td></tr><tr><td>f<sub>*</sub><sup>+</sup> = 1</td><td>f<sub>p</sub><sup>+</sup> = 0.8</td></tr><tr><td>c<sub>*</sub> = 129 km s<sup>-1</sup></td><td>c<sub>p</sub> = 10 km s<sup>-1</sup></td></tr><tr><td>M<sub>*</sub> = 1.01</td><td>M<sub>p</sub> = 1.2</td></tr><tr><td>Radial (x) direction</td><td>Vertical (y) direction</td></tr><tr><td>L<sub>x</sub>/R<sub>p</sub> = 40</td><td>L<sub>y</sub>/R<sub>p</sub> = 60</td></tr><tr><td>N<sub>x</sub> = 50, 100, 200,400,</td><td>N<sub>y</sub> = 75, 150, 300, 600,</td></tr><tr><td>800, 1600, 3200, 6400</td><td>1200, 2400, 4800, 9600</td></tr></tbody></table>

properties as listed in Table 1 are given for a stellocentric distance $r = r_{\text{launch},*} = 5R_{\odot}$, near the box's left edge. Here the stellar wind density, temperature, sound speed and flow speed are set to $n_* = 2.9 \times 10^4$ cm$^{-3}$, $T_* = 10^6$ K, $c_* = 129$ km s$^{-1}$ (computed for a mean molecular weight equal to half the proton mass, appropriate for an $f_*^+ = 100$ per cent ionized hydrogen plasma) and $v_* = 130$ km s$^{-1}$ (Sheeley et al. 1997; Quémérais et al. 2007; see also Lemaire 2011), respectively. Our stellar wind parameters are such that the implied spherically symmetric (3D) mass-loss rate is $1 \times 10^{12}$ g s$^{-1}$ or $2 \times 10^{-14}$ M$_{\odot}$ yr$^{-1}$.

Our stellar wind parameters are similar to those of the 'slow' solar wind in the Sun's equatorial plane. Compare our choices with those of Ekenbäck et al. (2010), who adopt a stellar wind speed of 450 km s$^{-1}$. Their speed is closer to that of the 'fast' solar wind which emerges from coronal holes. At solar minimum, the fast wind tends to be confined to large heliographic latitudes (polar regions), but at solar maximum, the coronal holes migrate to lower latitudes and the fast wind can more readily penetrate to the ecliptic (Kohl et al. 1998; McComas et al. 2003; S. Bale, private communication). Evaporating hot Jupiters like HD 209458b and HD 189733b have orbit normals that are nearly aligned with the spin axes of their host stars (Winn et al. 2005, 2006). Because such planets reside near their stellar equatorial planes, the slow equatorial solar wind seems a better guide than the fast, more polar wind; nevertheless, as noted above, the fast wind is known to extend to low latitudes, and the speeds and densities of both winds vary by factors of order-unity or more with time.

The stellar wind velocity at the left boundary is not plane-parallel but points radially away from the central star (located outside the box). The density, velocity and temperature in each cell at the boundary are computed by assuming that the central star emits a spherical isothermal wind whose velocity grows linearly with stellocentric distance $r$ and whose density decreases as $1/r^3$. These scalings, which are modelled after empirical solar wind measurements (e.g. Sheeley et al. 1997) and which maintain a constant mass-loss rate with $r$, are used only to define the left-edge boundary conditions and are not used in the simulation domain. Outflow boundary conditions are applied at the top, bottom and right edges of the box.

As a final comment about our choice of stellar wind parameters, we note that they are valid for the left-edge boundary at $r = r_{\text{launch},*} = 5R_{\odot}$ – not for the planet's orbital radius of $r = 10R_{\odot}$. The left-edge boundary must be far enough away from the planet that the stellar wind properties at the boundary are well-approximated by their 'free-stream' values in the absence of any planetary obstacle. We will see in Section 3 that the stellar wind slows considerably

between $r = 5R_{\odot}$ and $r = 10R_{\odot}$ as a consequence of the oncoming planetary wind. This region of deceleration is absent from the models of H08 and E10.

A circle of radius $d_{\text{launch},p} = 4R_p$, centred at position $(l_x, l_y) = (30R_p, 30R_p)$ (where the origin is located at the bottom left corner of the domain), defines the boundary where the assumed isotropic and radial planetary wind is launched. The properties of our simulated planetary wind, which are similar to those of the standard supersonic models of HD 209458b by García-Muñoz (2007) and Murray-Clay et al. (2009), are listed in Table 1, and are constant in time along the circular boundary. The density and velocity of the planetary wind at this boundary are such that if the wind were spherically symmetric, the mass-loss rate would be $1.6 \times 10^{11}$ g s$^{-1}$. This value lies within the range estimated from observations by Linsky et al. (2010) and from energetic considerations (e.g. Lecavelier des Etangs 2007; Ehrenreich & Désert 2011). Note that $1 - f_p^+ = 20$ per cent of the planetary wind at launch is neutral (Murray-Clay et al. 2009) and available for charge exchange. This neutral fraction represents a balance between photoionizations by extreme UV radiation and gas advection of neutrals at a planetary altitude of 4–5 $R_p$ (Murray-Clay et al. 2009). The planetary and stellar winds are barely supersonic at launch (Mach numbers $M_p = 1.2$ and $M_* = 1.01$).

Gravity is neglected, as are all rotational forces. The pressure $p$ is related to the internal energy density $e = E - \rho V^2/2$ via $p = (\gamma - 1)e$, where $\gamma = 1.01$. That is, gas is assumed to behave nearly isothermally. This isothermal assumption should not be taken to mean that the temperature is the same across the simulation domain; the temperature of the stellar wind at injection is $T_* = 10^6$ K, while that of the planetary wind is $T_p = 7000$ K.<fnref n="3" /> Rather, the two winds, as long as they remain unmixed, tend to maintain their respective temperatures as they rarefy and compress. In reality, the stellar wind can keep, in and of itself, a near-isothermal profile on length scales of interest to us because thermal conduction times (estimated, e.g. using the Spitzer conductivity) are short compared to dynamical times. Treating the planetary wind as an isothermal flow is less well justified, as cooling by adiabatic expansion can be a significant portion of the energy budget (García-Muñoz 2007; M09). Nevertheless, the error incurred by assuming the planetary wind is isothermal is small for our standard model because the planetary wind hardly travels beyond its launch radius of $4R_p$ before it encounters a shock; thus rarefaction factors are small. Furthermore, as noted above, shock compression factors are modest because the speed of the planetary wind is only marginally supersonic. Where the stellar and planetary winds meet and mix, the code ascribes an intermediate temperature $10^4$ K $< T < 10^6$ K. This temperature, as computed by HERACLES, is used only for the hydrodynamic evolution; it is not used for computing either the charge exchange reactions (Section 2.2) or the transmission spectrum (Section 2.3).

Each simulation is performed in two steps. First only the planetary wind is launched from its boundary and allowed to fill the entire domain for $2 \times 10^5$ s. Secondly, the stellar wind is injected through the left side of the box, by suitable assignment of ghost cells. This two-step procedure was found to minimize transients. The simulations typically run for $2 \times 10^6$ s, which corresponds to ~60 box-crossing times for the stellar wind in the horizontal direction.

<footnote n="3">³ In the standard model of M09, the temperature starts at ~$10^4$ K at a planetocentric radius of $1.1R_p$ – consistent with the observations by Ballester, Sing & Herbert (2007) – and cools to ~3000 K at $4R_p$. The temperatures calculated by García Muñoz (2007) at $4R_p$ are 6000–7000 K.</footnote>
## 2.2 Charge exchange

Charge exchange consists of the following forward and reverse reactions:

$$H_h^+ + H_c^0 \rightleftharpoons H_h^0 + H_c^+. \quad (2)$$

Hot (subscript h) ionized (superscript +) hydrogen emitted by the star can collide with cold (subscript c) neutral (superscript 0) hydrogen emitted by the planet, neutralizing the former and ionizing the latter while preserving their kinetic energies. The reverse reaction occurs with an identical rate coefficient β (units of cm³ s⁻¹; β is the cross-section multiplied by the relative velocity).

We have added reaction (2) to HERACLES by integrating the following equations in every grid cell (we refer to this portion of the calculation as the 'chemistry step'.):

$$\begin{aligned}
\frac{d(n_{\text{H}}x_{\text{h}}^{+})}{dt} &= \beta n_{\text{H}}^{2} (x_{\text{h}}^{0}x_{\text{c}}^{+} - x_{\text{h}}^{+}x_{\text{c}}^{0}) \\
\frac{d(n_{\text{H}}x_{\text{c}}^{0})}{dt} &= + \frac{d(n_{\text{H}}x_{\text{h}}^{+})}{dt} \\
\frac{d(n_{\text{H}}x_{\text{c}}^{+})}{dt} &= - \frac{d(n_{\text{H}}x_{\text{h}}^{+})}{dt} \\
\frac{d(n_{\text{H}}x_{\text{h}}^{0})}{dt} &= - \frac{d(n_{\text{H}}x_{\text{h}}^{+})}{dt} \\
x_{\text{h}}^{+} + x_{\text{h}}^{0} + x_{\text{c}}^{+} + x_{\text{c}}^{0} = 1.
\end{aligned} \quad (3)$$

Here n_H is the total hydrogen number density (regardless of ionization state or temperature), and $x_{(c,h)}^{(0,+)}$ is a number fraction (equivalently a mass fraction because the only element treated in the simulation is hydrogen). The rate coefficient $\beta = 4 \times 10^{-8}$ cm³ s⁻¹ is calculated by combining the energy-dependent cross-section of Lindsay & Stebbings (2005) with a Maxwellian distribution for the relative velocity between hydrogen atoms at the two (constant) temperatures $T_*$ and $T_p$. The finite-difference forms of equations (3) are

$$\begin{aligned}
x_h^{+(n+1)} - x_h^{+(n)} &= b \left( x_h^{0(n+1)} x_c^{+(n+1)} - x_h^{+(n+1)} x_c^{0(n+1)} \right) \\
&\quad x_c^{0(n+1)} - x_c^{0(n)} = x_h^{+(n+1)} - x_h^{+(n)} \\
&\quad x_c^{+(n+1)} - x_c^{+(n)} = -x_h^{+(n+1)} + x_h^{+(n)} \\
&\quad x_h^{0(n+1)} - x_h^{0(n)} = -x_h^{+(n+1)} + x_h^{+(n)} \\
&\quad x_h^{+(n)} + x_h^{0(n)} + x_c^{+(n)} + x_c^{0(n)} = 1
\end{aligned} \quad (4)$$

where the superscript (n) refers to the nth time-step, $b \equiv \beta n_{\text{H}} \Delta t$ and $\Delta t$ is the integration time-step of HERACLES. Because the right-hand side of the first of these equations is evaluated at step (n + 1) instead of step (n), our scheme is implicit. The first equation combines with the others to yield

$$x_h^{+(n+1)} = \frac{[x_h^{+(n)} + b(x_h^{+(n)} + x_c^{0(n)})(x_h^{+(n)} + x_c^{+(n)})]}{1 + b} \quad (5)$$

from which the remaining number fractions at time-step (n + 1) are derived. Because our solution is implicit, the dimensionless time-step b can exceed unity (as it does for our runs at low spatial resolution), and the system will still relax to its correct equilibrium. This chemical equilibrium is discussed further in Section 3.2.3.

Note that in contrast to H08 and E10, our calculations account for the reverse reaction $H_h^0 + H_c^+ \rightarrow H_h^+ + H_c^0$. Accounting for the reverse reaction helps us to avoid overestimating the amount of hot neutral hydrogen. Our calculations of $n_h^0$ are still overestimated, however, because we neglect thermal equilibration, i.e. cooling of

hot hydrogen by collisions with cold hydrogen. In Section 4.1 we estimate the error incurred to be of the order of unity.

Our calculations of the neutral fraction in the mixing layer do not explicitly account for photoionizations by Lyman continuum photons, radiative recombinations or advection of neutral hydrogen from the planetary wind – but these effects are already included by Murray-Clay et al. (2009) whose planetary wind parameters we use (see Section 2.1).

## 2.3 Lyman α absorption

The transmission spectrum in the Lyman α line is post-processed, i.e. calculated after HERACLES has finished running. Both hot and cold neutral hydrogen ($n_h^0$ and $n_c^0$) contribute to the Lyman α optical depth. It is assumed that the hot and cold neutral hydrogen do not thermally equilibrate (see Section 4.1 where we question this assumption). Thus, in computing the opacity due to hot hydrogen, we adopt a kinetic temperature of $T_* = 10^6$ K, and in computing the opacity due to cold hydrogen we take $T_p = 7000$ K. In each grid cell, the wavelength at line centre is Doppler shifted according to the horizontal component of the bulk velocity (the observer is to the far right of the simulation box). Voigt line profiles are used with a damping constant (Einstein A coefficient) equal to $\Gamma = 6.365 \times 10^8$ s⁻¹ (e.g. Verhamme, Schaeerer & Maselli 2006).

For each wavelength λ, the line-of-sight optical depth $\tau_\lambda(y)$ is evaluated along each horizontal row of cells pointing to the star (lying between the white dashed lines in Figs 1 and 2). The total absorption is then computed as

$$A(\lambda) = \langle 1 - \exp(-\tau_{\lambda}) \rangle \quad (6)$$

where $\langle \rangle$ denotes a 1D spatial average over y. Of course, the star actually presents a circular disc, but because the simulation is only 2D, our simple 1D average seems fair. The absorption profile $A(\lambda)$ can be computed for every snapshot (time-step) of the simulation.

## 2.4 Differences between this work and E10/H08

The main difference between our methods and those of E10/H08 is that we numerically solve the equations of hydrodynamics in a 2D geometry, whereas E10/H08 simulate collisions of hydrogen ‘meta-particles’ in a more kinetic, 3D treatment. Neither we nor they compute magnetic forces explicitly.

E10 include forces arising from the orbit of the planet about the star, including the Coriolis force, the centrifugal force and stellar tidal gravity. We do not. Our focus is on resolving mixing and charge exchange in the interface between the two winds. To that end, we solve for both the forward and reverse reactions of charge exchange (equations 2–3), whereas E10/H08 solve only for the forward reaction. Our equations permit a chemical equilibrium to be established in the mixing layer (see Section 3.2.3). Furthermore, the structure and geometry of the interaction region between the two winds are direct outcomes of our simulations, whereas the shape of the interface layer is imposed as a fixed ‘obstacle’ in the simulations of E10.

Other differences include our treatments of the planetary and stellar winds. We account for both the neutral and ionized components of the planetary wind; E10 assume the planetary outflow is purely neutral. We draw our parameters of the stellar wind from those of the slow equatorial solar wind, which blows at ~130 km s⁻¹ at a stellar centric distance of $r = 5 R_\odot$ (Sheeley et al. 1997; Quémerais et al. 2007). E10 take the stellar wind to blow at 450 km s⁻¹, while H08 take the stellar wind to blow at 50 km s⁻¹. Neither work accounts
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Snapshots of total density and velocity (left-hand panel), density of hot neutral hydrogen ($n_h^0$, middle panel) and temperature (right-hand panel) of the 50×75 simulation, with parameters listed in Table 1. Snapshots are taken at $t = 2 \times 10^6$ s. The temperature map shown in the right-hand panel is computed by HERACLES and used only to compute the hydrodynamic evolution; it is not used to compute the charge exchange reactions or the Lyman $\alpha$ spectrum (see Sections 2.2–2.3). The two dashed white lines represent sightlines to the stellar limbs.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Same as Fig. 1 and for the same simulation parameters but at a grid resolution of 3200 × 4800. Kelvin–Helmholtz rolls at the contact discontinuity first appear at this resolution. An even higher resolution of 6400×9600 yields the same star-averaged absorption (see Fig. 3).</figcaption></figure>

for how the stellar wind decelerates due to its interaction with the planetary wind, whereas in our simulations the deceleration zones are well-resolved.

We will review again our simulation methods, and assess the severity of our approximations, in Section 4.1.

# 3 RESULTS

Results for Lyman $\alpha$ absorption by the mixing layer, including numerical convergence tests and a direct comparison with observations, are given in Section 3.1. A parameter study is described in Section 3.2.

## 3.1 Absorption versus spatial resolution and time

In Figs 1 and 2, we present results at our lowest (50 × 75) and near-highest (3200 × 4800) spatial resolutions, respectively. The simulations agree on the basic properties of the flow. The planetary wind is launched from the red circle and encounters a bow shock,

visible in the left panels as a curved boundary separating orange (unshocked planetary wind) from red (shocked planetary wind). The radius of curvature of the planetary wind shock is roughly ~6$R_p$. Outside, the red region of thickness ~5$R_p$ contains shocked planetary wind.

The stellar wind encounters a weak shock – visible as a near-vertical line separating dark blue from lighter blue in the left-hand panels of Figs 1 and 2 – at a distance of ~5$R_p$ from the left edge of the box. The shocked stellar wind is diverted around the planet by the pressure at the stagnation point where the two winds collide head on.

We observe that both winds accelerate somewhat before they encounter shocks. For our standard model, the Mach numbers are $M_* \lesssim 1.3$ and $M_p \lesssim 1.5$ (for the parameter study simulations of Section 3.2, $M_p$ can grow up to 2–3). Density enhancements are thus modest – less than a factor of 2.

The contact discontinuity between the stellar and planetary winds separates light blue from dark red in the left panels. It is laminar at low resolution but breaks up into turbulent Kelvin–Helmholtz rolls at high resolution (cf. Stone & Proga 2009 whose spatial resolution
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Lyman $\alpha$ absorption $A$ (equation 6), evaluated at a Doppler-shift velocity of +100 km s$^{-1}$ from line centre, versus time and spatial resolution. The absorption converges in time for all simulations, but only for grid resolutions of 3200 × 4800 or greater does a unique value for the absorption emerge. The 3200 × 4800 simulation is also the lowest resolution run to resolve Kelvin–Helmholtz billows (see Fig. 2).</figcaption></figure>

was probably too low to detect the Kelvin–Helmholtz instability). The middle panels plot the density of hot neutral hydrogen produced by charge exchange in the mixing layer. The ‘head’ of the mixing layer, located near the stagnation point, spans only one or two grid cells in the low-resolution simulation. The high-resolution simulation resolves much better the head of the mixing layer. Zoomed-in snapshots of the head will be presented in Section 3.2.

In Fig. 3, the star-averaged absorption $A$ at an equivalent Doppler velocity of +100 km s$^{-1}$ (redshifted away from the observer) is plotted against time for a range of spatial resolutions. From $t = 0$ to $2 \times 10^5$ s, the planetary wind fills the simulation domain; the absorption quickly settles down to a value of ~2 per cent. At these early times, only cold ($T_p = 7000$ K) neutral hydrogen from the planet is available to absorb in Ly$\alpha$, and it is clearly insufficient to explain the absorption observed with *HST*.

Starting at $t = 2 \times 10^5$ s, the stellar wind is injected into the box. The absorption attains a first peak when the planetary and stellar winds reach a rough momentum balance and a mixing layer containing hot ($T_* = 10^6$ K) neutral hydrogen is established. The height of the first peak decreases with each factor of 2 improvement in grid resolution until a resolution of 3200 × 4800 is reached. Encouragingly, all of the absorption values calculated in the various simulations converge at late times.

The 3200 × 4800 run is the best behaved, with the absorption holding steady at $A \approx 9$ per cent for $10^6$ s. Compared to all other simulations at lower resolution, the 3200 × 4800 run is the only one in which Kelvin–Helmholtz rolls appear (more on the Kelvin–Helmholtz instability in Section 3.2.2).

We further tested the convergence of the 3200 × 4800 run by performing an even higher resolution simulation with 6400 × 9600

grid cells. Because of the expense of such a simulation, the initial conditions of the 6400 × 9600 run were taken from the 3200 × 4800 run at $t = 10^6$ s, and integrated forward for only $3 \times 10^5$ s (approximately nine box crossing times for the stellar wind in the horizontal direction). The absorption values versus time for the 6400 × 9600 run are overlaid in Fig. 3 and are practically indistinguishable from those of the 3200 × 4800 run. Having thus satisfied ourselves that the 3200 × 4800 run yields numerically convergent results, we will utilize this grid resolution (0.0125$R_p$ per grid cell length) for further experiments to understand the dependence of the absorption on input parameters, as described in Section 3.2.

Fig. 4 plots the absorption spectrum for our standard 3200 × 4800 simulation at $t = 2 \times 10^6$ s. The absorption $A$ is evaluated at wavelengths offset from the central rest-frame wavelength of the Lyman $\alpha$ transition by nine Doppler-shift velocities $\Delta v$. Absorption at $-50$ km s$^{-1}$ is stronger than at $+50$ km s$^{-1}$, a consequence of neutral, charge-exchanged hydrogen from the star accelerating from the stagnation point towards the observer. At larger velocities $|\Delta v| > 100$ km s$^{-1}$, the spectrum is more nearly reflection-symmetric about $\Delta v = 0$, because the broadening is purely thermal at $T_* = 10^6$ K.

Fig. 5 displays the same information as in Fig. 4 but in the full context of the *HST* observations. The agreement between the modelled and observed in-transit spectra is encouraging.

## 3.2 Scaling relations for absorption in the mixing layer

To understand how absorption in the mixing layer depends on input parameters, we performed three additional simulations varying the launch properties $n_p$, $n_*$, $v_p$ and $T_p$. The altered parameters are listed in Table 2. For all three simulations, the box size was maintained
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Lyman α absorption *A* versus Doppler-shift velocity *Δv* from line centre, evaluated for our standard 3200 × 4800 simulation at *t* = 2 × 10⁶ s. Absorption at −50 km s⁻¹ is stronger than at +50 km s⁻¹ because of the bulk motion of charge-exchanged neutral hydrogen streaming from the star towards the observer. The line wings at larger Doppler shifts are primarily thermally broadened at *T*<sub>*</sub> = 10⁶ K. The absorption *A* ≈ 9 per cent at *Δv* = ±100 km s⁻¹, in accordance with *HST* observations (see Fig. 5).</figcaption></figure>

at (L<sub>x</sub>, L<sub>y</sub>) = (40R<sub>p</sub>, 60R<sub>p</sub>) and the grid resolution was (N<sub>x</sub>, N<sub>y</sub>) = (3200, 4800).

Note that our parameter study is not exhaustive. For example, none of the simulations listed in Table 2 varies the Mach number at launch of either the planetary or stellar wind. Actually we

<table><caption>**Table 2** Launch parameters for the three additional 40R<sub>p</sub> × 60R<sub>p</sub> simulations at 3200 × 4800 resolution. The stellar parameters v<sub>*</sub> and T<sub>*</sub> are kept at their nominal values from Table 1. Note that none of the Mach numbers changes.</caption><thead><tr><th></th><th>Nominal</th><th>n<sub>p</sub>↑</th><th>n<sub>*</sub>↓</th><th>v<sub>p</sub>, T<sub>p</sub>↑</th></tr></thead><tbody><tr><td>n<sub>*</sub></td><td>2.9E4 cm<sup>-3</sup></td><td>2.9E4 cm<sup>-3</sup></td><td>9.7E3 cm<sup>-3</sup></td><td>2.9E4 cm<sup>-3</sup></td></tr><tr><td>n<sub>p</sub></td><td>3.9E6 cm<sup>-3</sup></td><td>1.2E7 cm<sup>-3</sup></td><td>3.9E7 cm<sup>-3</sup></td><td>3.9E6 cm<sup>-3</sup></td></tr><tr><td>v<sub>p</sub></td><td>12 km s<sup>-1</sup></td><td>12 km s<sup>-1</sup></td><td>12 km s<sup>-1</sup></td><td>√3×12 km s<sup>-1</sup></td></tr><tr><td>T<sub>p</sub></td><td>7000 K</td><td>7000 K</td><td>7000 K</td><td>21000 K</td></tr><tr><td>R</td><td>R<sub>0</sub> = 0.11</td><td>3 × R<sub>0</sub></td><td>3 × R<sub>0</sub></td><td>3 × R<sub>0</sub></td></tr></tbody></table>

have performed simulations varying the Mach number of the stellar wind. These behave as we would expect – in particular, increasing M<sub>*</sub> increases the amount of absorption because of the increased compression in the stellar shock. Nevertheless, we elect not to include these extra simulations in our parameter study below. Magnetic fields, neglected by our simulations but certainly present in the stellar wind if not also the planetary wind, would stiffen the gas and prevent the kind of compression that we see when we raise M<sub>*</sub>.

In the following subsections, we explain our numerical results on the properties of the mixing layer with order-of-magnitude scaling relations. The mixing layer’s location is analysed in Section 3.2.1; its thickness in Section 3.2.2; the densities of its constituent species in Section 3.2.3; and the column density and absorptivity of hot neutral hydrogen in Section 3.2.4.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Observed out-of-transit (highest blue curve) and observed in-transit (green curve) Lyman α spectra, reproduced from fig. 2 of Vidal-Madjar et al. (2003). In the line ‘core’ from −42 to +32 km s⁻¹, where interstellar absorption is too strong to extract a planetary transit signal, the flux is set to zero. Our theoretical in-transit spectrum (red curve) is computed by multiplying the observed out-of-transit spectrum by 1 − A, where A is plotted in Fig. 4. The agreement between the theoretical and observed in-transit spectra is good, supporting the idea that charge exchange between the stellar and planetary winds correctly explains the observed absorption at Doppler-shift velocities around ±100 km s⁻¹.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Zoomed-in snapshots of hot neutral hydrogen in the four simulations used to study how the properties of the mixing layer depend on input parameters (see Table 2). Snapshots are taken near $t = 2 \times 10^6$ for the standard model, and near $t = 1 \times 10^6$ s for the others. The bottom white dashed line is the line of sight to the lower stellar limb. The upper two red dashed lines bracket the ‘sampling interval’ over which the hot neutral hydrogen density is vertically averaged to produce the density profiles shown in Fig. 7 (the uppermost red dashed line is also the sightline to the upper stellar limb). In cases (b), (c) and (d), the mixing layer is located farther from the planet (centred at $l_x = 30R_p$) than is the case in (a). In cases (b) and (c), the horizontal thickness of the mixing layer is greater than in cases (a) and (d). In case (c), the density of hot neutral hydrogen is lowest. See Section 3.2 for explanations.</figcaption></figure>

### 3.2.1 Location of the mixing layer

Along the line joining the planet to the star, the mixing layer – equivalently, the contact discontinuity – is located approximately where the two winds reach pressure balance:

$$ \rho_{*} (v_{*}^{2} + c_{*}^{2}) = \rho_{p} (v_{p}^{2} + c_{p}^{2}). \qquad (7) $$

In equation (7), quantities are evaluated near the mixing layer, not at launch. Note further that in equation (7) and in equations to follow, we ignore the distinction between shocked and unshocked gas, as wind Mach numbers are near unity. Idealizing each wind velocity as constant, we substitute $\rho_p = j\dot{M}_p/(2\pi v_p d_p)$ and $\rho_* = M_*/(2\pi v_* d_*)$ into equation (7), as appropriate for the 2D circular winds in our simulations. Here $d_p$ measures distance from the planet, and $d_*$ measures distance from the star. Then the distance from the planet to the mixing layer – i.e. the approximate radius of curvature of the mixing layer – is given by

$$ d_p = d_* \mathcal{R} \qquad (8) $$

where

$$ \mathcal{R} = \frac{\dot{M}_p(v_p^2 + c_p^2)/v_p}{\dot{M}_*(v_*^2 + c_*^2)/v_*} \qquad (9) $$

For our standard model, $\mathcal{R} = \mathcal{R}_0 \approx 0.11$. Note that for 3D spherical winds, $d_p = d_* \sqrt{\mathcal{R}}$ (Stevens et al. 1992), but this relation is not relevant for our 2D Cartesian simulations – we will use (8) instead.

The parameters in Table 2 were chosen to increase $\mathcal{R}$ by a factor of 3 compared to its value in our fiducial model. By equation (8), when $\mathcal{R} = 3\mathcal{R}_0$, the mixing layer should be displaced three times farther away from the planet compared to its location in our standard model, assuming the star is far enough away that $d_*$ is essentially fixed at the star–planet separation. Fig. 6 displays zoomed-in snapshots of the mixing layers for all simulations in Table 2. Looking at the $l_x$-positions of the mixing layers, and recalling that the planet sits at $l_x = 30R_p$, we find that the layer is displaced $(30 - 8)/(30 - 21) \approx 2.4$ times farther away in the three new simulations as compared to the standard model. We consider this close enough to our expected

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Density profiles of charge-exchanged hot neutral hydrogen in the mixing layer, averaged vertically (between the dashed red lines in Fig. 6) and plotted against horizontal position. The planet is located to the right at $l_x = 30R_p$. The mixing layers of the three non-standard simulations are all displaced farther from the planet than in the standard model, a consequence of increasing the ratio $\mathcal{R}$ of the momentum carried by the planetary wind to that of the stellar wind (Section 3.2.1). The thicknesses of the mixing layers as shown by the red and green curves are larger than those shown by the blue and cyan curves, a consequence of changing the growth rate for the Kelvin–Helmholtz instability (Section 3.2.2). The two dashed lines are predictions of equation (16) based on considerations of chemical equilibrium; the blue, green and cyan curves correctly intersect the upper dashed line, while the red curve correctly intersects the lower dashed line (Section 3.2.3).</figcaption></figure>

factor of 3, given our neglect of the rather thick layers of shocked gas surrounding the mixing layer.

Fig. 7 shows density profiles for hot neutral hydrogen in the mixing layer for the three simulations plus our standard model. Densities are averaged over $l_y$ and plotted against $l_x$. The fact that the mixing layers in the simulations having $\mathcal{R} = 3\mathcal{R}_0$ align in position confirms that $\mathcal{R}$ is the dimensionless parameter controlling the location of the mixing layer.
### 3.2.2 Thickness of the mixing layer

Fig. 7 also indicates that the thickness of the mixing layer, $L_{\text{mix}}$, varies when we change input parameters. Empirically, we find that the variations are consistent with the relation

$$L_{\text{mix}} \sim 0.1 R_p \left( \frac{n_p}{n_*} \right)^{0.5} \quad (10)$$

where, as before, the distinction between shocked and unshocked gas densities is ignored.

We can rationalize (10) as follows. The time-scale for a mode of wavelength $\lambda_{KH}$ to grow exponentially by the linear Kelvin-Helmholtz instability (KHI) is given by

$$t_{KH} \sim \frac{\lambda_{KH}}{v_* - v_p} \frac{\rho_p + \rho_*}{2\pi(\rho_p \rho_*)^{1/2}} \sim \frac{\lambda_{KH}}{2\pi v_*} \sqrt{\frac{\rho_p}{\rho_*}} \quad (11)$$

(e.g. Chandrasekhar 1961). We assume that the thickness of the mixing layer saturates when a certain mode first becomes non-linear. Near saturation, the velocity perpendicular to the background shear flow becomes comparable to the shear flow velocity: $v_{\perp} \sim v_* - v_p \sim v_*$. Thus when the mode becomes non-linear, the mixing layer has thickness $L_{\text{mix}} \sim v_{\perp}t_{KH} \sim (\lambda_{KH}/2\pi)\sqrt{\rho_p/\rho_*}$. This result matches (10), if we assume the initial disturbance that develops into the mixing layer has a characteristic length scale that is fixed at $\lambda_{KH} = R_p$. Our description of mode saturation can only apply to locations not too far downstream from the stagnation point; far away, the flows are too strongly perturbed to be described by the linear growth time-scale (11).

### 3.2.3 Density of hot neutral hydrogen in the mixing layer

The density of hot neutral hydrogen in the mixing layer is set by chemical equilibrium. Suppose that within the layer, the total density $n_{\text{H, mix}}$ is approximately the average of the planetary wind density and the stellar wind density:

$$n_{\text{H,mix}} \sim \frac{n_p + n_*}{2} \sim \frac{n_p}{2}. \quad (12)$$

The densities in (12) are those of shocked gas, but as is the case for all of Section 3.2, we ignore for simplicity the difference in density between pre-shock and post-shock gas (see Section 3.2.1). Because $n_p \gg n_*$, charge exchange hardly alters the ionization state of the shocked – and still cold – planetary wind. That is, the values of $x_c^0$ and $x_c^+$ do not change as the dense planetary wind mixes with the dilute stellar wind. In particular, the ratio $x_c^0/x_c^+$ is fixed at its initial value of $(1 - f_p^+)/f_p^+ = 1/4$.

The time-scale for charge exchange is $(n_{\text{H, mix}}\beta)^{-1} \sim 10$ s, much shorter than the hours required for stellar-occulting gas to travel from the stagnation point to regions off the projected stellar limb. Thus nearly all of the gas seen in transit is driven quickly into chemical equilibrium, which from equation (3) demands that

$$\frac{x_h^0}{x_h^+} = \frac{x_c^0}{x_c^+} \quad (13)$$

$$= \frac{1 - f_p^+}{f_p^+} = \frac{1}{4}. \quad (14)$$

In other words, in the mixing layer, the ionization fraction of stellar wind material quickly slaves itself to the ionization fraction of planetary wind material. Now all of the hot hydrogen (both neutral

and ionized) in the mixing layer originates from the stellar wind; from (12), we have

$$(x_h^0 + x_h^+) n_{\text{H,mix}} \sim \frac{n_*}{2}. \quad (15)$$

Combining (14) with (15) yields

$$n_h^0 = x_h^0 n_{\text{H,mix}} \sim \frac{1}{2} (1 - f_p^+) n_* . \quad (16)$$

Equation (16) is approximately confirmed by our numerical results in Fig. 7; the horizontal dashed lines predicted by (16) roughly match the densities from our numerical simulations.

### 3.2.4 Column density and absorptivity of hot neutral hydrogen

Combining (10) with (16) gives the total column density of hot neutral hydrogen:

$$N_h^0 \sim n_h^0 L_{\text{mix}} \sim 0.05 (1 - f_p^+) (n_p n_*)^{1/2} R_p. \quad (17)$$

For our standard model, $N_h^0 \sim 3 \times 10^{13} \text{ cm}^{-2}$.

During planetary transit, the hot absorbing gas that covers the face of the star is located near the stagnation point. As such, the bulk line-of-sight velocity of transiting gas is much less than its thermal velocity, which is of the order of 100 km s$^{-1}$. Assuming that the gas is only thermally broadened, and that the gas is optically thin at wavelengths Doppler-shifted from line centre by velocities $\Delta v$, we construct an approximate, semi-empirical formula for the absorption:

$$A(\Delta v) \sim N_h^0 \sigma_{\text{line-ctr}} \exp[-m_H(\Delta v)^2/2kT_*] \quad (18)$$

$$\sim 0.1 \left( \frac{1 - f_p^+}{0.2} \right) \left( \frac{n_p}{4 \times 10^6 \text{ cm}^{-3}} \right)^{1/2} \left( \frac{n_*}{3 \times 10^4 \text{ cm}^{-3}} \right)^{1/2} \times \left( \frac{10^6 \text{ K}}{T_*} \right)^{1/2} \left( \frac{\exp[-m_H(\Delta v)^2/2kT_*]}{0.5} \right) \quad (19)$$

where $\sigma_{\text{line-ctr}} = 6 \times 10^{-15} (10^6 \text{ K}/T_*)^{1/2} \text{ cm}^2$ is the line-centre cross-section for the Lyman $\alpha$ transition, $m_H$ is the mass of the hydrogen atom and $k$ is Boltzmann's constant. Strictly speaking, the quantities in equation (19) should be evaluated in the vicinity of the contact discontinuity, but we have instead normalized equation (19) to the wind properties at launch (evaluated at $d_{\text{launch, p}} = 4R_p$ and $r_{\text{launch, *}} = 5 R_{\odot}$). We have verified in our simulations that the launch properties differ only by factors of the order of unity from the values at the contact discontinuity, and so equation (19) may be used to predict the absorption by inserting only the launch properties. The exponential in equation (19) is evaluated for nominal parameters $\Delta v = 100 \text{ km s}^{-1}$ and $T_* = 10^6 \text{ K}$.

As a further check, we show in Fig. 8 the absorption values $A$ to which the four simulations converge. They compare well with the values predicted by (19) using only the launch properties.

Had we kept the dependence of the mixing layer properties on the stellar wind Mach number $M_*$, equation (16) would be modified such that $n_h^0 \propto M_*^2 n_*$ – where $n_*$ is the pre-shock (launch) density – in accordance with the usual Rankine–Hugoniot jump condition that states that the density increases by the square of the Mach number across a plane-parallel isothermal shock. Equations (17)–(19) would be modified such that $A \propto N_h^0 \propto M_*$. Indeed, our numerical simulations (not shown) confirm this linear dependence of $A$ on $M_*$. We mention this result only in passing because it is not likely to remain true once we account for the real-life magnetization of the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Lyman $\alpha$ absorption A versus time, evaluated at a Doppler-shift velocity of +100 km s$^{-1}$ from line centre, for our standard model plus three additional models with different input parameters as indicated in the legend (see also Table 2). The coloured jagged lines are the results from our numerical simulations. The dashed black lines are the predictions from our physically motivated scaling relation (19); the simulations converge fairly well to the predicted values.</figcaption></figure>

stellar wind. Magnetic fields stiffen gas and reduce the dependence of A on $M_*$.

## 4 SUMMARY AND DISCUSSION

Using a 2D numerical hydrodynamics code, we simulated the collisional interaction between two winds, one emanating from a hot Jupiter and the other from its host star. The winds were assumed for simplicity to be unmagnetized. Properties of the stellar wind were drawn directly from observations of the equatorial slow solar wind (Sheeley et al. 1997; Quémerais et al. 2007; Lemaire 2011), while those of the planetary wind were taken from hydrodynamic models of outflows powered by photoionization heating (García-Muñoz 2007; Murray-Clay et al. 2009). For our standard parameters, the mass-loss rate of the star is $\dot{M}_* = 2 \times 10^{-14} \text{ M}_\odot \text{yr}^{-1} = 10^{12} \text{ g s}^{-1}$ and the mass-loss rate of the planet is $\dot{M}_p = 1.6 \times 10^{11} \text{ g s}^{-1} = 2.7 \times 10^{-3} \text{ M}_J \text{Gyr}^{-1}$. At the relevant distances, each wind is marginally supersonic – the stellar wind blows at $\sim 130-170 \text{ km s}^{-1}$ (sonic Mach number $M_* \lesssim 1.3$) and the planetary wind blows at $\sim 12-15 \text{ km s}^{-1}$ (Mach number $M_p \lesssim 1.5$). Thus, shock compression is modest, even without additional stiffening of the gas by magnetic fields.

A strong shear flow exists at the contact discontinuity between the two winds. At sufficiently high spatial resolution, we observed the interfacial flow to be disrupted by the Kelvin-Helmholtz instability. The Kelvin-Helmholtz rolls mix cold, partially neutral planetary gas with hot, completely ionized stellar gas. Charge exchange in the mixing layer produces observable amounts of hot ($10^6$ K) neutral hydrogen. Upon impacting the planetary wind, the hot stellar wind acquires, within tens of seconds, a neutral component whose fractional density equals the neutral fraction of the planetary wind (about $1 - f_p^+ = 20$ per cent). Seen transiting against the star, hot neutral hydrogen in the mixing layer absorbs $\sim 10$ per cent of the light in the thermally broadened wings of the stellar Lyman $\alpha$ emission line, at Doppler shifts of $\sim 100 \text{ km s}^{-1}$ from line centre. Just such a transit signal has been observed with the HST (Vidal-Madjar et al. 2003). The $\pm 100 \text{ km s}^{-1}$ velocities reflect the characteristic velocity dispersions of protons in the stellar wind – as inferred from in situ spacecraft observations of the solar wind (e.g. fig. 3 of Marsch 2006).

Our work supports the proposal by Holmström et al. (2008) and Ekenbäck et al. (2010) that charge exchange between the stellar and planetary winds is responsible for the Ly$\alpha$ absorption observed by HST. This same conclusion is reached by Lecavelier des Etangs et al. (2012) in the specific case of HD 189733b. Our ability to reproduce the observations corroborates the first-principles calculations of hot Jupiter mass loss on which we have relied (e.g. Yelle 2004; García-Muñoz 2007; Murray-Clay et al. 2009, M09). Time variations in Ly$\alpha$ absorption are expected both from the variable stellar wind – the solar wind is notoriously gusty – and from the variable planetary wind, whose mass-loss rate tracks the time-variable UV and X-ray stellar luminosity.

## 4.1 Neglected effects and directions for future research

Although the general idea of photoionization-powered planetary outflows exchanging charge with their host stellar winds seems correct, details remain uncertain. We list below some unresolved issues, and review the effects that our simulations have neglected, in order of decreasing concern.

(i) **Thermal equilibration in the mixing layer.** Our calculations overestimate the amount of hot neutral hydrogen produced by charge exchange because they neglect thermal equilibration. A hot neutral hydrogen atom cools by colliding with cold gas, both ionized and neutral, from the planetary wind. The concern is that hot neutral gas cools before it transits off the face of the star. Starting from where the mixing layer is well-developed (say the lower red dashed line in Fig. 6), hot neutral gas is advected off the projected stellar limb in a time

$$t_{\text{adv}} \sim 2R_p/v_* \sim 2 \times 10^3 \text{s}. \quad (20)$$

By comparison, the cooling time is of the order of

$$t_{\text{cool}} \sim \frac{1}{n_c^+ \sigma v_{\text{rel}}} \\ \sim 500 \left( \frac{2 \times 10^6 \text{ cm}^{-3}}{n_c^+} \right) \left( \frac{10^{-16} \text{ cm}^2}{\sigma} \right) \left( \frac{100 \text{ km s}^{-1}}{v_{\text{rel}}} \right) \text{s}, \quad (21)$$

where $n_c^+$ is the density of cold ionized hydrogen in the mixing layer, $v_{\text{rel}}$ is the relative speed between hot and cold hydrogen, and $\sigma$ is the H-H$^+$ cross-section for slowing down fast hydrogen, here taken to be the ‘viscosity’ cross-section calculated by Schultz et al. (2008).<fnref n="4" /> Our estimate of $t_{\text{cool}}$ in (21) neglects cooling by neutral-neutral collisions, but we estimate the correction to be small, as $n_c^0$ is lower than $n_c^+$ by a factor of $1/(1 - f_p^+) \sim 5$, and the cross-section for H-H collisions is generally not greater than for H-H$^+$ collisions (A. Glassgold, private communication; see also Swenson, Tupa & Anderson 1985; note that E10 take the relevant H-H cross-section to be $10^{-17} \text{ cm}^2$ but do not provide a reference).

That $t_{\text{cool}} \sim t_{\text{adv}}$ indicates our simulated column densities of hot neutral hydrogen may be too large, but hopefully not by factors of more than a few. Keeping more careful track of the velocity distributions – and excitation states – of neutral hydrogen in the mixing layer would not only improve upon our calculations of Lyman $\alpha$ absorption, but would also bear upon the recent detection

<footnote n="4">⁴For slowing down fast H in a sea of cold H$^+$, there may also be a contribution to $\sigma$ from ‘momentum transfer’ in ‘elastic’ (non-charge-exchange) collisions. This contribution increases $\sigma$ over the viscosity cross-section by only $\sim 30$ per cent; compare figs 6 and 7 of Schultz et al. (2008).</footnote>
of Balmer Hα absorption in the hot Jupiters HD 209458b and HD 189733b (Jensen et al. 2012).

(ii) Magnetic fields. Insofar as our results depend on Kelvin–Helmholtz mixing, our neglect of magnetic fields is worrisome because magnetic tension can suppress the Kelvin–Helmholtz instability (Frank et al. 1996). For numerical simulations of magnetized planetary winds interacting with magnetized stellar winds, see Cohen et al. (2011a,b). These magnetohydrodynamic simulations can track how planetary plasma is shaped by Lorentz forces, but as yet do not resolve how the planetary wind mixes and exchanges charge with the stellar wind.

(iii) Dependence of Lyα absorption A on the planetary wind density $n_p$. In the same vein as item (ii), we found empirically that $A \propto n_p^{1/2}$, and argued that this result arose from the Kelvin–Helmholtz growth time-scale. Ekenbäck et al. (2010) found a much weaker dependence: increasing $n_p$ by a factor of 100 only increases A in their models by a factor of ~2 at −100 km s⁻¹ and even less at positive velocities – see their figs 8 and 9. The true dependence of A on $n_p$ remains unclear.

(iv) Rotational effects and gravity. There are a few order-unity geometrical corrections that our study is missing. Our standard stellar wind velocity of $v_*$ = 130 km s⁻¹ is comparable to the planet's orbital velocity of $v_{\text{orb}}$ = 150 km s⁻¹, so that in reality the stellar wind strikes the planet at an angle of roughly 45°. The Coriolis force will also deflect the planetary wind by an order-unity angle after a dynamical time of $r/v_{\text{orb}} \sim 5 \times 10^4$ s, by which time the wind will have travelled ~5R<sub>p</sub> from the planet. These geometrical effects are potentially observable [see e.g. Schneiter et al. (2007) and Ehrenreich et al. (2008) for modelling of HD 209458b, and Rappaport et al. (2012) for a real-life example of a transit light curve that reflects the trailing comet-tail-like shape of the occulting cloud]. However, these geometrical effects seem unlikely to change the basic order of magnitude of the absorption $A \sim 10$ per cent that we have calculated.

We have also neglected planetary gravity, stellar tidal gravity and the centrifugal force, all of which can change the planetary wind velocity. But this omission seems minor, since we have drawn our input planetary wind velocities from calculations that do account for such forces (M09), at least along the substellar ray. According to fig. 9 of M09, the planetary wind accelerates from $v_p \approx 10$ km s⁻¹ at a planetocentric distance $d = 4R_p$, to $v_p \approx 30$ km s⁻¹ at $d = 10R_p$. This range of velocities and corresponding distances overlap reasonably well with the range of velocities and distances characterizing our simulations.

(v) Hydrodynamic approximations for the stellar and planetary winds. We have not formally justified our use of the hydrodynamic equations to describe the wind–wind interaction. The problem is that the collisional mean-free path in the stellar wind is much longer than the length scales of the flow: $\lambda_{\text{Coulomb, *}} = 1/(n_*\sigma_{\text{Coulomb}}) \sim 10^{13}(10^4 \text{ cm}^{-3}/n_*)(10^{-17} \text{ cm}^2/\sigma_{\text{Coulomb}}) \text{ cm}$, where $\sigma_{\text{Coulomb}} \sim 10^{-17}(T_*/10^6 \text{ K})^{-2} \text{ cm}^2$ is the cross-section for protons scattering off protons. The fact that the solar wind is collisionless and does not necessarily admit a one-fluid treatment is well-known.

Nevertheless, it is perhaps just as well-known that Parker's (1958, 1963) use of the fluid equations to describe the collisionless solar wind is surprisingly accurate, capturing the leading-order features of the actual solar wind. The role of Coulomb collisions in relaxing the velocity distribution functions of protons and electrons is fulfilled instead by plasma instabilities and wave–particle interactions – see e.g. reviews of solar wind physics by Marsch, Axford & McKenzie (2003) and Marsch (2006). The gross properties of collisionless shocks can still be modelled with the hydrodynamic

equations insofar as those properties depend only on the macroscopic physics of mass, momentum and energy conservation, and not on microphysics (e.g. Shu 1992).

Note that the planetary wind is fully collisional because of its higher density and lower temperature, and modelling it as a single fluid appears justified: $\lambda_{\text{Coulomb, p}} \sim 10^7(10^6 \text{ cm}^{-3}/n_p)(T_p/10^4 \text{ K})^2 \text{ cm}$, which is smaller than any other length scale in the problem.

(vi) Non-Maxwellian behaviour of the stellar proton velocity distribution. Lyman α absorption at the redshifted velocity of +100 km s⁻¹ arises from charge-exchanged neutral hydrogen at the assumed stellar wind temperature of 10⁶ K. We have assumed a Maxwellian distribution function for hydrogen in the stellar wind, and have ignored non-Maxwellian features that have been observed in the actual solar wind, including high-energy tails and temperature anisotropies. Accounting for non-Maxwellian behaviour may introduce order-unity corrections to our results for the absorption. For the more polar fast solar wind, proton temperatures parallel to and perpendicular to the solar wind magnetic field differ by factors of a few at heliocentric distances of 5–10 solar radii (McKenzie, Axford & Banaszcziewicz 1997). For the more equatorial slow solar wind – which our simulations are modelled after – temperature anisotropies are more muted (Marsch et al. 2003, p. 391).

(vii) Stellar radiation pressure. Stellar Lyman α photons can radiatively accelerate neutral hydrogen away from the star (e.g. Vidal-Madjar et al. 2003; M09). Both the planetary wind and the charge-exchanged stellar wind in the mixing layer are subject to a radiation pressure force that exceeds the force of stellar gravity by a factor $\beta$ on the order of unity.

Radiative repulsion of the charge-exchanged stellar wind in the mixing layer may not be observable, because once hot neutral hydrogen is created in the mixing layer, it is advected off the projected limb of the star before radiation pressure can produce a significant velocity: $\delta v_{\text{rad}} \sim GM_*/r^2 \times \beta \times t_{\text{adv}} \sim 6\beta \text{ km s}^{-1}$, which does not exceed the hot neutral hydrogen's thermal velocity of $\sim 100 \text{ km s}^{-1}$.

What about radiative acceleration of the planetary wind? The travel time of the planetary wind from the planet to the mixing layer is $\sim 10R_p/v_p \sim 10^5$ s, long enough for neutral hydrogen to attain radiative blow-out velocities in excess of 100 km s⁻¹. However, the amount of hydrogen that suffers radiative blow-out is limited to the column that presents optical depth unity to Lyman α photons. This column is $1/\sigma_{\text{line-ctr}} \sim 2 \times 10^{13}(T_p/10^4 \text{ K})^{1/2} \text{ cm}^{-2}$, and is much smaller than the typical column in the planetary wind, which is $(1 - f_p^+) n_p R_p \sim 10^{16} \text{ cm}^{-2}$. Thus the bulk of the planetary wind is shielded from radiative blow-out, and our neglect of radiation pressure appears safe.

Note that Lecavelier des Etangs et al. (2012) find that radiation pressure cannot explain the largest blueshifted velocities observed for HD 189733b; like us, they favour charge exchange.

## ACKNOWLEDGMENTS

This project was initiated during the 2011 International Summer Institute for Modelling in Astrophysics (ISIMA) program hosted by the Kavli Institute for Astronomy and Astrophysics (KIAA) at Beijing University. We thank Pascale Garaud and Doug Lin for organizing this stimulating workshop, and Edouard Audit, Stuart Bale, Andrew Cumming, Sebastien Fromang, Pascale Garaud, Astrid Lamberts, Eliot Quataert and Jim Stone for helpful exchanges. We are indebted to Al Glassgold for his guidance in helping us understand H–H and H–H⁺ collisions, and an anonymous referee for a report that improved the presentation of this paper. ISIMA is funded
by the American and Chinese National Science Foundations, the Center for the Origin, Dynamics and Evolution of Planets and the Center for Theoretical Astrophysics at the University of California at Santa Cruz, the Silk Road Project and the Excellence Cluster Universe at TU Munich. Support for EC was provided by NASA through a Hubble Space Telescope Theory grant. This work was given access to the HPC resources of [CCRT/CINES/IDRIS] under the allocations c2011042204 and c2012042204 made by GENCI (Grand Equipement National de Calcul Intensif).

REFERENCES

Ballester G. E., Sing D. K., Herbert F., 2007, Nat, 445, 511

Ben-Jaffel L., 2007, ApJ, 671, L61

Ben-Jaffel L., 2008, ApJ, 688, 1352

Chandrasekhar S., 1961, in Marshall W., Wilkinson D. H., eds, Hydrodynamical and Hydromagnetic Stability. Dover Publications, New York

Clarke C. J., Carswell R. F., 2003, Principles of Astrophysical Fluid Dynamics. Cambridge Univ. Press, Cambridge

Cohen O., Kashyap V. L., Drake J. J., Sokolov I. V., Gombosi T. I., 2011a, ApJ, 733, 67

Cohen O., Kashyap V. L., Drake J. J., Sokolov I. V., Gombosi T. I., 2011b, ApJ, 738, 166

Ehrenreich D., Désert J.-M., 2011, A&A, 529, 136

Ehrenreich D., Lecavelier des Etangs A., Hébrard G., Désert J.-M., Vidal-Madjar A., McConnell J. C., Parkinson C. D., Ballester G. E., 2008, A&A, 483, 933

Ekenbäck A., Holmström M., Wurz P., Grießmeier J.-M., Lammer H., Selsis F., Penz T., 2010, ApJ, 709, 670 (E10)

Fossati L. et al., 2010, ApJ, 714, L222

Frank A., Jones T. W., Ryu D., Gaalaas J. B., 1996, ApJ, 460, 777

García Muñoz A., 2007, Planet. Space Sci., 55, 1426

González M., Audit E., Huynh P., 2007, A&A, 464, 429

Holmström M., Ekenbäck A., Selsis F., Penz T., Lammer H., Wurz P., 2008, Nat, 451, 970 (H08)

Jaritz G. F., Endler S., Langmayr D., Lammer H., Grießmeier J.-M., Erkaev N. V., Biernat H. K., 2005, A&A, 439, 771

Jensen A. G., Redfield S., Endl M., Cochran W. D., Koesterke L., Barman T., 2012, ApJ, 751, 86

Kohl J. L. et al., 1998, ApJ, 501, L127

Lamberts A., Fromang S., Dubus G., 2011, MNRAS, 418, 2618

Lecavelier des Etangs A., 2007, A&A, 461, 1185

Lecavelier des Etangs A., Vidal-Madjar A., McConnell J. C., Hébrard G., 2004, A&A, 418, L1

Lecavelier Des Etangs A. et al., 2010, A&A, 514, A72

Lecavelier des Etangs A. et al., 2012, A&A, 543, L4

Lemaire J. F., 2011, preprint (arXiv e-prints)

Liang M.-C., Parkinson C. D., Lee A. Y.-T., Yung Y. L., Seager S., 2003, ApJ, 596, L247

Lindsay B. G., Stebbings R. F., 2005, J. Geophys. Res., 110, A12213

Linsky J. L., Yang H., France K., Froning C. S., Green J. C., Stocke J. T., Osterman S. N., 2010, ApJ, 717, 1291

McComas D. J., Elliott H. A., Schwadron N. A., Gosling J. T., Skoug R. M., Goldstein B. E., 2003, Geophys. Res. Lett., 30, 24

McKenzie J. F., Axford W. I., Banaszkiewicz M., 1997, Geophys. Res. Lett., 24, 2877

Marsch E., 2006, Liv. Rev. Solar Phys., 3, 1

Marsch E., Axford W. I., McKenzie J. F., 2003, in Dwivedi B. N., ed., Dynamic Sun. Cambridge Univ. Press, Cambridge, p. 374

Moses J. I. et al., 2011, ApJ, 737, 15

Murray-Clay R. A., Chiang E. I., Murray N., 2009, ApJ, 693, 23

Parker E. N., 1958, ApJ, 128, 664

Parker E. N., 1963, in Marshak R. E., ed., Interplanetary Dynamical Processes. Interscience Publishers, New York

Quémerais E., Lallement R., Koutroumpa D., Lamy P., 2007, ApJ, 667, 1229

Rappaport S. et al., 2012, ApJ, 752, 1

Schneiter E. M., Velázquez P. F., Esquivel A., Raga A. C., Blanco-Cano X., 2007, ApJ, 671, L57

Schultz D. R., Krstic P. S., Lee T. G., Raymond J. C., 2008, ApJ, 678, 950

Sheeley J. et al., 1997, ApJ, 484, 472

Shu F. H., 1992, in Kelly A., ed., Physics of Astrophysics, Vol. II. University Science Books, Sausalito

Stone J. M., Proga D., 2009, ApJ, 694, 205

Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386, 265

Swenson D., Tupa D., Anderson L., 1985, J. Phys. B: Atomic Mol. Phys., 18, 4433

Tian F., Toon O. B., Pavlov A. A., De Sterck H., 2005, ApJ, 621, 1049

Verhamme A., Schaerer D., Maselli A., 2006, A&A, 460, 397

Vidal-Madjar A., Lecavelier des Etangs A., Désert J.-M., Ballester G. E., Ferlet R., Hébrard G., Mayor M., 2003, Nat, 422, 143

Vidal-Madjar A. et al., 2004, ApJ, 604, L69

Vidal-Madjar A., Lecavelier des Etangs A., Désert J.-M., Ballester G. E., Ferlet R., Hébrard G., Mayor M., 2008, ApJ, 676, L57

Weissman P. R., McFadden L.-A., Johnson T. V., 1999, in Helé J., ed., Encyclopedia of the Solar System. Academic Press, New York

Winn J. N. et al., 2005, ApJ, 631, 1215

Winn J. N. et al., 2006, ApJ, 653, L69

Yelle R. V., 2004, Icarus, 170, 167

Yelle R. V., 2006, Icarus, 183, 508

This paper has been typeset from a *TEX/LATEX* file prepared by the author.