Effects of Radiation Pressure on the Evaporative Wind of HD 209458b

TL;DR

This study uses high-resolution 3-D hydrodynamic simulations to assess whether radiation pressure influences the photoevaporative wind of HD 209458b, concluding it is unlikely to be a significant factor within expected flux limits.

Contribution

The paper provides the first detailed 3-D hydrodynamic simulation analysis of radiation pressure effects on HD 209458b's wind, confirming its limited role in driving high-velocity atmospheric escape.

Findings

Radiation pressure is unlikely to significantly influence the wind at expected Lyman-$\alpha$ fluxes.

Simulations confirm analytic limits on radiation pressure effects.

High velocities in observed winds are not explained by radiation pressure alone.

Abstract

The role of radiation pressure in shaping exoplanet photoevaporation remains a topic of contention. Radiation pressure from the exoplanet's host star has been proposed as a mechanism to drive the escaping atmosphere into a "cometary" tail and explain the high velocities observed in systems where mass loss is occurring. In this paper we present results from high-resolution 3-D hydrodynamic simulations of a planet similar to HD 209458b. We self-consistently launch a wind flowing outward from the planet by calculating the ionization and heating resulting from incident high-energy radiation, and account for radiation pressure. We first present a simplified calculation, setting a limit on the Lyman-$\alpha$ flux required to drive the photo-evaporated planetary material to larger radii and line-of-sight velocities. We then present the results of our simulations, which confirm the limits determined by our analytic calculation. We thus demonstrate that, within the limits of our hydrodynamic simulation and for the Lyman-$\alpha$ fluxes expected for HD 209458, radiation pressure is unlikely to significantly affect photoevaporative winds or to explain the high velocities at which wind material is observed, though further possibilities remain to be investigated.

Tables

Table 1. Table 1: Run parameters

		No-Flux	Low-Flux	Int.-Flux	High-Flux
Planet Radius	$R_p\left(R_J\right)$	$1.529$
Planet Mass	$M_p\left(M_J\right)$	$0.73$
Planet Temperature	$T_p\left(\text{K}\right)$	$3\times {10}^3$
Planet Surface Density	${\rho }_p\left({\text{g cm}}^{-3}\right)$	$1.625\times {10}^{-15}$
Stellar Mass	$M_{\star }\left(M_{☉}\right)$	$1.23$
Stellar Radius	$R_{\star }\left(R_{☉}\right)$	$1.19$
Stellar Ionizing Flux	$F_{0,UV}\left({\text{phot cm}}^{-2}{\text{s}}^{-1}\right)$	$2\times {10}^{13}$
Stellar Lyman-$\alpha$ Flux	$F_{0,\alpha }\left({\text{phot cm}}^{-2}{\text{s}}^{-1}\right)$	$0$	$4.1\times {10}^{13}$	$4.1\times {10}^{14}$	$4.1\times {10}^{15}$
Orbital Separation	$a$ (AU)	$0.047$
Orbital Period	$P$ (days)	$3.525$
Orbital Velocity	$\Omega$ (rad/day)	$1.78$
Polytropic Index	$\gamma$	$\frac{5}{3}$

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathbf{v}=0,$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}\right)=-\nabla p-\rho \nabla \varphi +{\mathbf{f}}_{\mathbf{R}}+\mathcal{K},$$

$$\frac{\partial E}{\partial t}+\nabla \cdot ((E+p)\mathbf{v})=-\rho \mathbf{v}\cdot \nabla \varphi +{\mathbf{f}}_{\mathbf{R}}\cdot \mathbf{v}+\mathcal{G}-\mathcal{L}+\mathcal{J},$$

$$\frac{\partial n_{\text{mbox}H\text{sc}i}}{\partial t}+\nabla \cdot (n_{\text{mbox}H\text{sc}i}\mathbf{v})=\mathcal{R}-\mathcal{I},$$

$$\frac{\partial n_{\text{mbox}H\text{sc}ii}}{\partial t}+\nabla \cdot (n_{\text{mbox}H\text{sc}ii}\mathbf{v})=\mathcal{I}-\mathcal{R},$$

$$\mathcal{K}={\sigma }_{\alpha }{n}_{\text{mbox}H\text{sc}i}{F}_{\alpha }(x)p_{\gamma }\hat{x}$$

$$\mathcal{J}=\mathcal{K}\cdot \mathbf{v}$$

$${\rho }_{atm}(r)={\rho }_p{\left[\frac{R_0{R}_p}{R_0-R_p}\left(\frac{1}{r}-\frac{1}{R_0}\right)\right]}^{\frac{1}{\gamma -1}}$$

$$P_{atm}(r)={\rho }_p^{1-\gamma }{c}_{s,p}^2{\rho }^{\gamma },$$

$$\rho (r)=\{\begin{array}{ll}{\rho }_{atm}(R_{mask}),& r<R_{mask}\\ {\rho }_{atm}(r),& R_{mask}\le r\le R_{ob}\\ {\rho }_{atm}(R_{ob})\times 10^{-4},& r>R_{ob}\end{array},$$

$$P(r)=\{\begin{array}{ll}{P}_{atm}(R_{mask}),& r<R_{mask}\\ P_{atm}(r),& R_{mask}\le r\le R_{ob}\\ P_{atm}(R_{ob})\times 10^{-2},& r>R_{ob}\end{array}.$$

$$\Delta v=\sqrt{\frac{2G{M}_{\star }}{a}}=2.16\times 10^7\text{mbox}cm/s.$$

$$\frac{d{p}_x}{dt}=\dot{M}\Delta v=\dot{M}\sqrt{\frac{2G{M}_{\star }}{a}}.$$

$$\frac{d{p}_x}{dt}=FA{p}_{\gamma },$$

$$F_{th}=\sqrt{\frac{2G{M}_{\star }}{a}}\frac{\dot{M}}{A{p}_{\gamma }}=\sqrt{\frac{2G{M}_{\star }}{a}}\frac{\dot{M}}{2\pi aH{p}_{\gamma }}.$$

$$\frac{M_p}{M_{\star }}{a}^2=\frac{r_H^2}{a}[r_H+a(\beta -1)+2r_H\beta ],$$

Full text

††pagerange: Effects of Radiation Pressure on the Evaporative Wind of HD 209458b–References

1 Introduction

Characterization of exoplanetary atmospheres is now a dynamic field of research. Much of the literature has focused on atmospheric structure, composition, and the detection of biosignatures, but the interaction of a planetary atmosphere with the circumplanetary environment remains an important area of inquiry, as it plays a key role in determining the long-term evolution of planetary atmospheres. In particular, atmospheric escape, caused by the absorption of stellar X-ray and extreme ultraviolet radiation by the atmosphere of close-in planets, is of particular interest, both theoretically (e.g. Lammer et al., 2003; Yelle, 2004; Lecavelier des Etangs et al., 2004; Tian et al., 2005; García Muñoz, 2007; Murray-Clay et al., 2009; Koskinen et al., 2013; Khodachenko et al., 2017) and observationally (e.g. Vidal-Madjar et al., 2003; Fossati et al., 2010; Lecavelier des Etangs et al., 2012; Ehrenreich et al., 2015; Bourrier et al., 2018).

From the theoretical side, many physical processes can affect the evolution of evaporative planetary winds, including stellar and planetary magnetic fields (Matsakos et al., 2015; Villarreal D’Angelo et al., 2018; Owen & Adams, 2014; Khodachenko et al., 2015) and time-dependent phenomena such as flares (Bisikalo et al., 2018), coronal mass ejections (Cherenkov et al., 2017), stellar variability (Lecavelier des Etangs et al., 2012; Bisikalo & Cherenkov, 2016), and atmospheric circulation (Teyssandier et al., 2015).

Observationally, photoevaporative outflows have been directly studied in some detail for a handful of planets. These include HD 209458b (Vidal-Madjar et al., 2003), HD 189733b (Lecavelier des Etangs et al., 2010, 2012; Bourrier et al., 2013), GJ 436b (Kulow et al., 2014; Ehrenreich et al., 2015), and GJ 3470b (Bourrier et al., 2018). For many of these planets, asymmetric Hi absorption of the stellar Lyman-$\alpha$ emission line by material of planetary origin has been observed, with enhanced absorption seen in the blue wing when compared to the red wing. The high velocities of the absorbing material of up to $\pm 150$ km/s observed in these systems challenge models of exoplanet evaporation, which predict lower planetary wind speeds of order the planetary escape speed ($\sim 30$ km/s).

Several processes have been suggested as the cause of this asymmetric high-velocity absorption, including charge exchange between the stellar and planetary wind (Holmström et al., 2008; Ekenbäck et al., 2010; Tremblin & Chiang, 2013; Kislyakova et al., 2014; Bourrier et al., 2016; Christie et al., 2016), confinement by the stellar wind (Schneiter et al., 2007; Schneiter et al., 2016; McCann et al., 2019), and acceleration of neutrals by radiation pressure from stellar Lyman-$\alpha$ emission (Lecavelier Des Etangs et al., 2008; Schneiter et al., 2016; Bourrier & Lecavelier des Etangs, 2013; Khodachenko et al., 2017; Cherenkov et al., 2018). The differing conclusions resulting from these studies show that it is not yet clear which if any of these processes can fully account for the observations. Here we address this question by focusing on one process: radiation pressure from the host star.

A number of studies have supported the idea that radiation pressure can change the photoevaporative flows. Schneiter et al. (2016), for example, performed detailed simulations including both the star and the planet, finding that the planetary wind was driven back into a cometary tail. However, in these simulations both planet and star were modeled as isotropic sources, with the planetary wind launched from outside the planet’s atmosphere at 3 planetary radii ($R_{p}$). In addition, radiation pressure from the star was modeled as a reduction in the stellar gravity throughout the grid. Bourrier & Lecavelier des Etangs (2013) also evolved the wind from between the Roche lobe and the planet’s surface isotropically. However, while they also reduced the stellar gravity to model radiation pressure, they included self-shielding from absorption of the Lyman-$\alpha$ photons. In addition, Bourrier & Lecavelier des Etangs (2013) removed hydrogen particles from the simulation once they were ionized. This study also found that radiation pressure produces a cometary tail-like structure in the planetary wind, with much of the gas escaping the system.

On the other hand, Khodachenko et al. (2017) have recently presented 2-dimensional hydrodynamic simulations of a four-fluid model of HD 209458b, using a simplified plateau-like Lyman-$\alpha$ profile and a velocity-dependent absorption of Lyman-$\alpha$ photons by neutral hydrogen. The planetary wind is launched by heating from EUV radiation and interacts with the stellar wind through charge exchange. They find that radiation pressure does not significantly affect the planetary wind. Cherenkov et al. (2018) calculated the ionization and radiation pressure self-consistently in a 3D simulation at moderate resolution, using a hydrogen envelope to simulate HD 209458b. They find that the planetary wind would be affected by Lyman-$\alpha$ radiation only if it were two orders of magnitude greater than that predicted for HD 209458b.

Thus the question of the role of radiation pressure on exoplanetary photoevaporative flows remains quite open. In this paper we present results from high-resolution 3-D hydrodynamic simulations of a planet similar to HD 209458b, launching the wind self-consistently by calculating the ionizing radiation and accounting for radiation pressure. In section 2 we present the computational method and parameters used in our simulations. In section 3 we show the results of our simulations. In section 4, we compare our results to the analytic calculation and the results of previous studies and suggest possible reasons for our lack of a significant cometary tail and high-velocity absorption signature. We conclude in section 5.

2 Methods and Model

Our simulations were conducted with AstroBEAR111https://astrobear.pas.rochester.edu/ (Cunningham et al., 2009; Carroll-Nellenback et al., 2013), a massively parallelized adaptive mesh refinement (AMR) code that includes a variety of multiphysics solvers, such as self-gravity, heat conduction, magnetic resistivity, radiative transport, and ionization dynamics. The equations solved for these simulations are those of fluid dynamics in a rotating reference frame, including gravitational effects of both planet and star:

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the mass density, $\boldsymbol{v}$ is the fluid velocity, $p$ is the thermal pressure, $\phi$ is the gravitational potential, $\boldsymbol{f_{R}}$ combines the Coriolis and centrifugal forces so that $\boldsymbol{f_{R}}=\rho\left(-2\boldsymbol{\Omega}\times\boldsymbol{v}-\boldsymbol{\Omega}\times\left(\boldsymbol{\Omega}\times\boldsymbol{r}\right)\right)$ (where $\boldsymbol{\Omega}$ is the orbital velocity), $E=p/(\gamma-1)+\rho v^{2}/2$ is the combined internal and kinetic energies, $\mathcal{G}$ and $\mathcal{L}$ are respectively the heating and cooling rates, $\boldsymbol{\mathcal{K}}$ is the force due to Lyman-$\alpha$ radiation, and $\mathcal{J}$ is the kinetic energy added by the radiation pressure from Lyman-$\alpha$ radiation.

The simulation also tracked the advection, photoionization, and recombination of neutral and ionized hydrogen. We use the photon-conserving update scheme from Krumholz et al. (2007) to solve the following equations:

[TABLE]

[TABLE]

where $n_{\mbox{\small{H}\sc{i}}}$ is the number density of neutral hydrogen, $n_{{\mbox{\small{H}\sc{ii}}}}$ is the number density of ionized hydrogen, and $\mathcal{R}$ and $\mathcal{I}$ are the recombination and ionization rates.

2.1 Radiation transfer

We perform the radiation transfer as in Debrecht et al. (2019). In addition to the ionizing flux, we apply a monochromatic Lyman-$\alpha$ flux $F_{0,\alpha}$. We calculate its radiative transfer using the same method used for the ionizing flux, but with an absorption cross-section $\sigma_{\alpha}=5.9\times 10^{-14}\left(10^{4}\mbox{ K}/T\right)^{1/2}$. Each absorbed Lyman-$\alpha$ photon deposits its entire momentum into the gas, so the (volumetric) force due to radiation pressure is

[TABLE]

and the rate of kinetic energy deposition is

[TABLE]

where $p_{\gamma}=10.2\mbox{ eV}/c$ is the momentum per Lyman-$\alpha$ photon, $p_{0,x}$ is the initial momentum of the cell in the $x$ direction, $dt$ is the time step of the radiation calculation, and $\rho$ is the density of the cell. Note that the gas density does not change during the radiation calculation, and that our radiation is represented by a single frequency bin, thereby ignoring the Doppler shift of the absorption line. Consequences of this assumption are discussed further in section 4.1.

2.2 Planet atmosphere model

In our simulations, we have modeled the planet as a sphere of hydrogen in hydrostatic equilibrium as in Debrecht et al. (2019), with a density profile of

[TABLE]

and pressure profile of

[TABLE]

where $R_{p}$ is the radius of the planet and $R_{0}$ is the radius where a physical atmosphere would end. For completeness, we now describe the internal planetary boundaries used in the simulation to model the planet. $R_{ob}$, three cells inside $R_{0}$, is the radius at which the planetary profile is initially cut off and the ambient conditions begin. This is done to reduce numerical noise at the interface between the edge of the planet and the ambient material. $R_{ob}=1.04R_{p}$. $R_{ib}=R_{0}/(1+(R_{0}-R_{p})/(R_{p}e^{3(\gamma-1)}))$ is the radius at which the planetary profile is reset at every time step in order to replenish the supply of atmospheric hydrogen flowing into the outer layers of the planet from its core. $R_{ib}=0.68R_{p}$. $R_{mask}$, set five cells within $R_{ib}$, is the radius at which the planetary profile is cut off to prevent the singularity at $r=0$. $R_{mask}=0.62R_{p}$. As long as the radius from which the wind is launched (the $\tau=1$ radius) does not reach $R_{ib}$, the structure of the wind will not be significantly affected by any of these choices (see Appendix A of Murray-Clay et al., 2009). In addition, the mass lost from the planet during the simulation, approximately $10^{15}$ g at $5\times 10^{9}$ g/s, is only about 1% of the total initial mass, suggesting that the structure of the layers does not have significant bearing on the evolution of the wind. $\rho_{p}$ is the density and $c_{s,p}^{2}$ is the sound speed at the planetary surface $R_{p}$.

In addition, we now set up a static ambient pressure two orders of magnitude lower than the final value of the planet profile, rather than matching the pressure at the boundary. The complete specification of the initial conditions of the simulation is thus

[TABLE]

with corresponding pressures of

[TABLE]

A sketch of the different layers is shown in figure 1.

2.3 Description of simulation

The input parameters of the simulations were chosen to model a planet similar to HD 209458b. Our planet has a radius of $1.529R_{J}$, chosen to be 1.1 times the measured broadband optical transit radius $R_{tr}=1.39R_{J}$ (based on the discussion in section 2 of Murray-Clay et al. (2009) of the predicted wind launch radius), a mass of $0.73M_{J}$, and a temperature of 3000 K. While this is higher than the predicted equilibrium temperature of 1500 K, it has been shown that the planet’s surface temperature has a negligible effect on the wind structure for any temperature below the wind temperature at the launching surface (Murray-Clay et al., 2009, their Appendix A), which is generally $\sim 10^{4}\mbox{ K}$. The planet orbits the host star of mass $1.23M_{\astrosun}$ and radius $1.19R_{\astrosun}$ at a separation of $a=0.047\mbox{ AU}$. Planetary and stellar parameters were taken from Stassun et al. (2017).

The Cartesian simulation domain ranges from $[-10,-10,-10]R_{p}$ to $[10,10,10]R_{p}$, with the planet centered at $[0,0,0]$. We apply outflow-only extrapolating boundary conditions at all boundaries, with the initial ambient conditions applied if the extrapolated conditions would result in inflow. The simulation has a base resolution of $105^{3}$ and 4 levels of additional refinement, giving an effective resolution of $1680^{3}$, which resolves the entire grid to $0.012R_{p}$. The maximum resolution is forced in a sphere out to the termination radius $R_{0}$ of the hydrostatic atmosphere (defined further in Section 2.2). The planetary radius is therefore resolved by 84 cells. We allow the mesh to evolve outside of the planet based on the density gradient. The stellar location is not included in these simulations; only its radiation and gravitational effects are considered. In order to isolate the effect of radiation pressure, we have not included a stellar wind in these simulations.

We have run three simulations, varying the Lyman-$\alpha$ flux by an order of magnitude between each. Our fiducial value of $4.1\times 10^{14}\mbox{ phot cm}^{-2}\mbox{s}^{-1}$ was calculated from the value of the Lyman-$\alpha$ flux of HD 209458 measured by Wood et al. (2005). The hydrogen density $\rho_{p}=1.625\times 10^{-15}\mbox{ g cm}^{-}3$ at the planetary radius $R_{p}$ was adjusted so that the surface of unit optical depth to ionizing radiation at 16 eV for the un-irradiated planet, along the substellar line, was near the planet radius $R_{p}$.

For our base simulation, we evolved the planet wind for 4.58 days (1.3 planetary orbits) to a steady state without Lyman-$\alpha$ flux. After introducing the Lyman-$\alpha$ flux, we ran the simulations for 2.5 days (0.72 planetary orbits) for the low- and intermediate-flux runs and 3.4 days (0.97 planetary orbits) for the high-flux run, after which the low- and intermediate-flux simulations (runs 1 and 2 in table 1) had reached a steady state, by which we mean the flow had achieved a stable ionization front and wind morphology. The high-flux simulation (Run 3 in table 1), however, enters a cycle whereby the wind drives away from the planet, is disrupted by radiation pressure, and drives away from the planet again, repeating the cycle. This simulation was run for a time sufficient to see two repetitions of this cycle.

2.4 Assumption of collisionality

As AstroBEAR is a hydrodynamics code, we implicitly assume that the fluid is collisional throughout the simulation domain, and can therefore be treated hydrodynamically rather than kinetically. We check this assumption by plotting the Knudsen number, $Kn=\lambda/L$, for the final states of the no-flux, intermediate-flux, and high-flux simulations, as shown in figure 2. A gas is collisional for $Kn\lesssim 1$. Here the mean free path $\lambda=(\sigma_{col}n)^{-1}$ and the characteristic length scale $L=0.1R_{p}$. We plot both the Knudsen number due to proton-proton (Coulomb) collisions (top row), where $\sigma_{col}=10^{-13}(10^{4}/T)^{2}$ and the number density $n=n_{{\mbox{\small{H}\sc{ii}}}}$, and the Knudsen number due to neutral-proton, neutral-neutral, and proton-proton collisions (bottom row), all assumed to be hard-body collisions where $\sigma_{col}=3.53\times 10^{-16}$ and the number density $n=n_{{\mbox{\small{H}}}}$.

These plots show that the entire planetary wind is collisional, as also found in e.g. Salz et al. (2016), save for a small portion of the blown-off material in the high-flux simulation. In particular, the arms of the wind are collisional primarily thanks to proton-proton interactions, while the base of the wind is collisional due to the high neutral density. The lack of collisionality in the ambient material can be disregarded, as the ambient does not affect the evolution of our simulations. The average mean free path at $2R_{p}$ using the hard-body assumption is $2.46\times 10^{8}\mbox{ cm}$, resulting in an average Knudsen number at $2R_{p}$ of $0.225$.

3 Results

3.1 Analytic treatment

To understand the influence of radiation pressure on the escaping planetary atmosphere, we first discuss an analytic approximation for the stellar Lyman-$\alpha$ flux required to completely unbind a torus of material, hereafter referred to as the ”threshold flux”. It was found in previous studies that hydrodynamic planetary winds will often form an extended torus-like structure filling a significant portion of the orbit (Matsakos et al., 2015; Carroll-Nellenback et al., 2017; Debrecht et al., 2018). We therefore take a simplification of such a torus as the basis for our analytic treatment of the effects of radiation pressure. We assume that the planet is injecting material into the torus with a mass-loss rate of $\dot{M}=5.7\times 10^{9}\mbox{ g/s}$, which is that found in our initial-condition simulation with no Lyman-$\alpha$ flux. We remark that this value agrees with that of most other theoretical estimates (e.g., see Table 1 of Kubyshkina et al., 2018), as well as with observational estimates (e.g. Vidal-Madjar et al., 2003). In addition, we assume that the torus is optically thick to Lyman-$\alpha$ radiation, so that all stellar Lyman-$\alpha$ flux incident upon it is absorbed, maximizing the absorbed momentum.

To liberate material from the torus and form a cometary tail, it must be supplied with a $\Delta v$ of approximately the escape speed from the central star at the orbital distance, that is,

[TABLE]

To liberate material as it is introduced into the torus, we require a force of

[TABLE]

The force due to the stellar Lyman-$\alpha$ flux is

[TABLE]

where the area on which the flux acts is the rectangular cross-section of the torus, $A=2\pi aH$, with $H$ being the full height of the torus and $F$ the flux. Equating these forces, we find that the threshold flux $F_{th}$ is given by

[TABLE]

With a torus height of $H\approx 6R_{p}$, we find $F_{th}\approx 8\times 10^{14}\mbox{ phot cm}^{-2}\mbox{ s}^{-1}$. This is a minimum threshold, as it applies when all particles of the torus acquire only up to the escape speed. Were most of the acquired momentum instead carried off by a small fraction of the torus, a significantly higher flux would be needed. It should be noted that this analysis does not treat the acceleration of small portions of the neutral torus material to the observed velocities of $\sim 10-100$ km/s.

3.2 Simulation results

All of our simulations are centered on the planet and carried out in the planet’s orbiting frame of reference, with the orbital velocity vector $\boldsymbol{\Omega}$ in the $+z$ direction. Therefore, ”up-orbit” is approximately in the $+y$ direction and ”down-orbit” is approximately in the $-y$ direction. We first discuss the case without Lyman-$\alpha$ radiation pressure on the torus, which is the steady state from which the other cases are begun. We then summarize the intermediate-flux (Run 2 of table 1) and high-flux (Run 3 of table 1) cases. As the low-flux and intermediate-flux cases are essentially identical, we do not present the results of the low-flux case. Movies for each of these cases are available on the AstroBEAR YouTube channel222https://www.youtube.com/user/URAstroBEAR . Much of the following discussion is similar to that from our previous study on the parameter space of planetary winds (Debrecht et al., 2019).

3.2.1 Initial steady-state: no Lyman-$\alpha$ flux

The no-flux case is shown in figure 3. The leftmost panels of this figure show the density (hue, logarithmically scaled), in $\mbox{g cm}^{-3}$, and velocity field (quivers, normalized to the maximum velocity in the plane) for this simulation. In these panels we also show the $\tau=1$ surface for our ionizing radiation (black contour), Mach surface (magenta contour), and nominal planetary radius $R_{p}$ (green contour). The simulation is carried out in the co-rotating frame. In the top row we show a 2-D slice looking down on the orbital plane (hereafter referred to as the top view). In the bottom row we show a slice in the orbital plane looking up-orbit (hereafter referred to as the side view).

The $\tau=1$ surface denotes the position from which the wind is launched. At this point, most of the ionizing photons have been absorbed. Below the $\tau=1$ surface the planetary hydrogen density and the corresponding optical depth quickly increase such that by $R=0.94R_{p}$, 99.9% of the incident radiation has been absorbed. Murray-Clay et al. (2009) found that the wind solution is insensitive to conditions below the $\tau=1$ surface. Although the details of the flow below $\tau=1$ are not expected to accurately model conditions of a real giant planet, this region still plays a role in the simulation by providing a flow of neutral material to larger radii that is subsequently ionized and continually supplies the wind.

The $\tau=1$ surface also denotes the extent of the ionization shadow of the planet. The center panels of figure 3 show the neutral hydrogen fraction (hue, logarithmically scaled) for the same slices through the planet. Here we see that the planet’s ionization shadow results in material leaving the night side and remaining neutral over a significant distance. This distance is defined by the ionization timescale and the wind velocity. Also seen in this panel is that the rest of the planetary wind carries a small fraction of neutrals with it. The neutral fraction $n_{\mbox{\small{H}\sc{i}}}/n_{{\mbox{\small{H}}}}$ is of order $10^{-2}$ in the bulk of the wind.

The right panels in figure 3 show the temperature, in K, on a logarithmic scale. Because the wind is mostly transparent to ionizing radiation outside of the planet’s atmosphere, the temperature of the wind is determined hydrodynamically, beginning at $T=8200\mbox{ K}$ at the base of the wind and cooling primarily by expansion (though radiative and recombination cooling are still present). At the surface of the planet, slight asymmetry can be seen, highlighting the effect of the ionizing radiation on the equilibrium structure of the wind.

As in Carroll-Nellenback et al. (2017), Matsakos et al. (2015), and Debrecht et al. (2019), we see the formation of up- and down-orbit wind trajectories. Material otherwise bound by the planet’s Hill sphere is driven into a wind by pressure forces (from ionization heating) and tidal forces. Rotational forces then turn this material into the two arms of the wind. Because of the rotational forces, the wind does not fill the whole computational domain. Rather, it is confined to a torus with a quasi-cylindrical cross section, as shown in the bottom row.

We now consider the flow pattern in the model via figure 4, created by convolving random noise integrated along the streamlines of the velocity field with a color plot of the density. In this figure, hue shows density (as in the left panels of figure 3) and the texture of the plot represents flow streamlines. The left panel shows the top view, while the right panel shows the side view. In the top view, we can see the origin of the material in each of the arms of the wind. The up-orbit arm originates from gas in the Hill sphere to the left of the red line, while the down-orbit arm originates from gas to the right of the red line. In addition, vortices formed at the points of divergence of the arms are apparent from this plot.

The side view again shows vortices at the corners of the wind. In addition, material flowing laterally around the planet after being heated can be seen joining the bulk of the wind.

3.2.2 Comparison with previous studies

The steady state of this wind differs significantly from that seen in Debrecht et al. (2019), with much of the difference due to the significantly denser planet. Movies of the no-flux simulation (see footnote 2) show that the up-orbit and down-orbit arms form slowly, with material that is peeled out of the Hill sphere of the planet by tidal forces. In contrast with the previous study, these arms are not significantly redirected by the Coriolis force after formation, so that there is little complex structure to the planetary wind. The neutral tail, composed of material pulled from the night side of the planet, still forms as it did in the simulations of Debrecht et al. (2019). But as the planetary wind is not turned completely around the planet by Coriolis forces the tail is not redirected by any shocks and merely flows along with the down-orbit arm. Thus it continues to exist farther away (in the $-y$ direction) from the planet than previously seen. Finally, the Mach surfaces of the wind form at a significant distance from the ionization front in the simulation of HD 209458b, while in the previous study the ionization front and Mach surface were nearly identical.

It is also interesting to compare our simulation with that of Cherenkov et al. (2018), particularly their figure 1. In broad strokes, their flow appears similar to ours, with up-orbit and down-orbit arms formed by material escaping the Roche lobe of the planet. There are some subtle differences, perhaps the most significant of which is that their mass loss is much more focused through the L1 and L2 points than ours. In the simulations shown here the mass loss appears much broader in terms of the planetary surface area (though still averaging to the L1 and L2 points). The remaining differences are likely due to their inclusion of the stellar wind, which appears to force the up-orbit arm back and results in turbulence that somewhat disrupts the down-orbit arm.

3.2.3 Intermediate Lyman-$\alpha$ flux

Figure 5 shows the steady state of the intermediate-flux run, which uses the predicted actual Lyman-$\alpha$ flux of the host star HD 209458. It can be seen that the intermediate-flux case differs only slightly from the initial steady state. The most apparent difference is that instabilities in the trailing edges of the arms of the wind are suppressed by the increase in external pressure, leading to arms that are effectively wider by approximately one planetary radius. On the other hand, instabilities along the leading edge of the wind are enhanced by the additional pressure coming from the star. However, no cometary tail forms, with both the up-orbit and down-orbit arms forming as before.

The similarities between the no-flux and intermediate-flux case can also be seen in the velocity streamlines plotted in Figure 6. The origins of the up- and down-orbit arms are nearly identical to the no-flux case, with the exception that some material from near the substellar point flows completely around the planet to join the down-orbit material. In addition, the side view shows that the high-velocity portions of the wind are closer to the planet, outlining the Hill sphere more clearly in the z direction. The lateral flow around to the night side of the planet is also somewhat enhanced.

We have plotted the acceleration of wind material due to radiation pressure in figure 7, in units of km/s/$t_{c}$, where $t_{c}=20\mbox{ hr}$ is the crossing time for half the box (10 $R_{p}$). Here $0.72\mbox{ km/s/}t_{c}=1\mbox{ cm/s}^{2}$. The left panel shows the top view, while the right panel shows the side view. Note that the maximum acceleration, located along the edge of the up-orbit arm around $5R_{p}$, is around 100 cm/s2, while the acceleration from the stellar gravity at the orbit of HD 209458b is approximately 330 cm/s2. Because the wind is collisional (see section 2.4), this acceleration is acting not only on the material directly subject to the radiation pressure, but also on the material deeper in the wind.

3.2.4 High Lyman-$\alpha$ flux

The high-flux case does not reach a steady state. Instead, the up-orbit arm is pushed back periodically to join into the down-orbit arm, forming a cometary tail. We note, however, that the pressure from Lyman-$\alpha$ radiation is insufficient to completely suppress the up-orbit arm. Eventually, ionized material from the planet’s Hill sphere begins to form an up-orbit arm again, which is then eventually blown back into a cometary tail. The intermittency of the cometary tail distinguishes our result from the permanent cometary tail seen in Bourrier & Lecavelier des Etangs (2013) and Schneiter et al. (2016).

The timescale for this cycle is the crossing time for the Hill sphere. The Hill sphere is enlarged due to the effects of radiation pressure. Modifying the standard Hill sphere derivation, we find:

[TABLE]

where $\beta$ is the reduction in stellar gravity due to radiation pressure ($\beta=(g_{\star}-a_{rad})/g_{\star}$) and $r_{H}$ is the radius of the Hill sphere. The average size of the Hill sphere for this simulation varies, but is approximately this size of the simulation domain for much of the simulation period. The timescale for the cycling of this simulation is therefore approximately the crossing time of the simulation domain.

Figure 8 shows the cycle using plots of the density (hue), mach number (magenta contour), velocity field (quivers), and surface of optical depth equal to one (black contour). The top row shows the top view, while the bottom row shows the side view, with time progressing from left to right. First we examine the top view. We take the left panel to be the initial state ($t=0$ hours) for the purpose of this description. Here the radiation pressure has suppressed the up-orbit arm. Material on the up-orbit side of the planet has essentially been confined to the Hill sphere. Material that has managed to expand outwards is almost entirely subsonic.

The second panel, at $t\approx 18$ hours, shows the expansion of material from the Hill sphere to re-form a supersonic up-orbit arm. However, this arm does not extend off the grid as in the previous simulations. Instead, as we see in the third panel ($t\approx 30$ hours), as the arm expands its inertia and thermal pressure drop, allowing radiation pressure to drive it away from the star. Finally, the right panel shows the arm being blown completely backward at $t\approx 40$ hours, with most of the material flowing around the planet to join the down-orbit arm. The material on the up-orbit side is once again subsonic and inside the Hill sphere.

We now examine the side view. In the left panel, we see a relatively symmetric core of the wind with a quasi-cylindrical cross-section as in the previous simulations. We also see a large flow of diffuse material off the grid on the far side of the wind. In the second panel, the substellar portion of the wind is compressed such that its cross-section is no longer symmetric. Most of the diffuse material has blown off the grid or been recaptured by the central flow. The third panel shows material on the upper and lower edges of the wind essentially ablating due to radiation pressure. The final panel shows this material once again being driven off the grid in a diffuse stream.

As was done for the intermediate flux case, we have plotted the acceleration of the wind material due to radiation pressure in figure 9. The panels are arranged as described for figure 8. Note that in this case the maximum acceleration is greater than 1000 cm/s2 along the starward edge of the wind, significantly greater than the stellar acceleration of 330 cm/s2. As in the intermediate-flux case, the wind is collisional (see section 2.4), and the acceleration is therefore acting on material deeper in the wind in addition to the material acted upon directly by radiation pressure.

3.3 Synthetic observations

We also compute synthetic observations for the no-flux, intermediate-flux (Run 2), and high-flux (Run 3) cases, as described in Carroll-Nellenback et al. (2017). Because they are done in post-processing, where computational expediency is not an issue (unlike the radiation pressure calculations, which are done on the fly), we include Doppler shifting of Lyman-$\alpha$ flux. The stellar Lyman-$\alpha$ profile used for the synthetic observations is always that taken from Bourrier & Lecavelier des Etangs (2013), and is not scaled with the incident flux used in the calculation of the radiation pressure. This does not affect the width of the line absorption (i.e. observed velocity), which is the principle metric of interest in the study.

In the top panels of figure 10 we show the synthetic Lyman-$\alpha$ observations for these three cases, including their variability over time. Thick lines represent observations within one standard deviation of the mean time-variable absorption, while the dashed line gives the observed average flux. Standard deviations were taken over 24.4, 18.3, and 39.66 hours for the no-flux, intermediate-flux, and high-flux cases, respectively. We see that the variability of the observations is more noticeable in the high-flux case, due to the cyclic nature of the outflow. We note that the no-flux case has a slightly larger variability than the medium-flux case. This is likely due to the previously-mentioned suppression of the edge instabilities by the external pressure.

The no-flux and intermediate-flux cases both show enhanced absorption in the red wind pre-transit and the blue wing post-transit, due to the presence of the up-orbit and down-orbit arms of the wind, respectively. Note that the line-of-sight velocity of the up-orbit arm is away from the observer, while the line-of-sight velocity of the down-orbit arm is toward the observer, while the line-of-sight component of the orbital velocity is opposite these directions. Also note the absence of enhanced absorption in the pre-transit red wing and a greater spread in velocity in the post-transit blue wing for the high-flux case, due to the periodic absence of the up-orbit arm and the acceleration of material in the down-orbit arm.

It should be noted that the out-of-transit absorptions are artificially lowered by the size of our simulation box, as there are some rays from the camera to the star that either do not or only partially pass through the simulation. 100% of the rays pass completely through the simulation only at $\pm 17$ minutes. The out-of-transit observations should therefore be interpreted cautiously. Note however that this restriction does not alter the principle conclusion that there is no high velocity component to the outflow.

The bottom panels of figure 10 show the fractional absorption of the Lyman-$\alpha$ line for each simulation. It can be seen that the absorption is confined almost entirely to the center of the line, with only the planetary transit having an effect on the wings in the lower-flux cases. The absorption differs slightly in the high-flux case, where a small proportion of excess absorption can be seen in the blue wing out to almost 50 km/s. We remark that nearly all of the absorption lies in the region covered by the interstellar medium and geocoronal emissions. This is a consequence of radiation pressure not driving sufficient material to high enough velocities.

Figure 11 is similar to Figure 10, but shows the transmission fraction as a function of time after transit (orbital angle) and velocity (wavelength). Again we see that the bulk of the absorption is confined to the center of the Lyman-$\alpha$ line, with deeper absorption extending to approximately -40 km/s in the high-flux case. Finally, we note that the Lyman-$\alpha$ absorption found in these simulations is lower than that found in Carroll-Nellenback et al. (2017), which is due to a significantly higher ionization fraction for the bulk of the wind in the current simulations, which in turn is due to an ionization timescale in the optically-thin wind of a factor of $\sim 4$ shorter in the current simulations.

4 Discussion

4.1 Comparison with the analytic treatment

In our analytic calculation, we show that we expect a flux of $\sim 8\times 10^{14}\mbox{ phot/s/cm}^{2}$ to be required to drive the planetary outflow away (i.e. into a cometary tail. See section 3.1). Note that this is for a full torus of material, while in our simulations the ”torus” has an angular extent of $\approx\pi/10$, with the ends of the torus outside the simulation boundaries. Estimating an actual extent of $\pi/2$, we would have a threshold flux $F_{th}=3.2\times 10^{15}\mbox{ phot/s/cm}^{2}$. In our simulations, we find that a flux greater than our maximum of $4.1\times 10^{15}\mbox{phot/s/cm}^{2}$ is required to create a permanent cometary tail; however, we find that this flux, which is only $1.3\times$ the calculated threshold flux for the truncated torus, is sufficient to periodically drive the wind into a cometary tail. Thus our analytic model correctly suggests that the observed flux from HD 209458 is too low to drive cometary tails but underestimates the flux that would be needed to do so.

This underestimation may be due to two possibilities: i) the threshold flux calculated above is sufficient to blow the planetary wind away, but insufficient to completely suppress the wind; ii) the starward layer which absorbs most of the Lyman-$\alpha$ flux in our simulations is accelerated to significantly greater than escape velocity, representing an extra source of lost energy unaccounted for in the analytical calculation.

The first possibility can be thought of in terms of the ram and thermal pressure of the wind interacting with the radiation pressure from the star in two regimes: one where $P_{Ly-\alpha}\lesssim P_{th}+P_{ram}$ and one where $P_{Ly-\alpha}>P_{th}+P_{ram}$. During the period of bubbling out (right panels in figures 8 and 9), the combination of the ram and thermal pressures of the planetary wind near the sonic surface is sufficient to just counter the radiation pressure from the stellar Lyman-$\alpha$ line, leading to the expansion of material and the formation of the up-orbit arm. During the period of blowback (second panels of figures 8 and 9), the expansion of the wind means radiation pressure becomes greater than the sum of the planetary wind pressures and so is able to drive the up-orbit arm to larger radii.

Figures 7 and 9 show that the radiation pressure acts on a thin shell before the optical depth is sufficient to completely absorb the Lyman-$\alpha$ flux. This shell is blown around the edges of the wind while still blocking significant portions of the Lyman-$\alpha$ flux from being absorbed by the core of the outflow. This phenomenon is particularly apparent in the high-flux case (figure 9). How significant this shielding is to the formation of a cometary tail can be quantified by comparing the average speed of a particle being significantly affected by radiation pressure to the escape speed $v_{esc}=2.16\times 10^{7}\mbox{cm/s}$ of the system. In the intermediate-flux case, this average velocity is $1.6\times 10^{6}\mbox{ cm/s}$, while in the high-flux case, it is $3.32\times 10^{6}\mbox{ cm/s}$. It is therefore likely that material driven along the edges of the wind is a moderate, though not necessarily significant, source of energy loss not accounted for in the analytic treatment.

In addition, figure 12 plots contours of density over the speed of the gas in the x direction for the no-flux, intermediate-flux, and high-flux cases. The contours are as follows: green, $1.00\times 10^{-16}$ g/cm3; red, $3.16\times 10^{-18}$ g/cm3; blue, $1.00\times 10^{-19}$ g/cm3; magenta, $3.16\times 10^{-21}$ g/cm3; white, $1.00\times 10^{-22}$ g/cm3. No gas exceeds a velocity of 100 km/s. The left panels show that the bulk of the up-orbit and down-orbit arms are between 10 and 50 km/s with no incident Lyman-$\alpha$ flux. The center panels indicate that the inclusion of the intermediate flux has almost no effect on the gas speed throughout the simulation. On the other hand, the right panels show that during some phases of the cycle, small amounts of neutral material are accelerated to velocities greater than 50 km/s. However, as we have seen in figure 10, the density is too low for absorption at these velocities to be significant, resulting in no absorption at the observed velocities of $\sim\pm 100$ km/s. Also note that the orbital escape speed is $\sim 200$ km/s; therefore, as discussed above, we confirm that even the high flux is insufficient to accelerate any of the wind to stellar escape speed within the space of our simulation.

Importantly, we assume in our radiation transfer that the Lyman-$\alpha$ photons are absorbed regardless of gas velocity, and thus ignore the Doppler shifts. Our calculations therefore provide an upper limit on the efficacy of radiation pressure in the formation of a cometary tail. As has been seen in figures 7 and 9, all of the Lyman-$\alpha$ radiation is absorbed in a thin layer along the starward edge of the planetary wind. Most material along the leading edge has a speed between 10 and 50 km/s, with the majority of that in the $y$ or $z$, rather than $x$, directions. If we were to include the Doppler broadening in our absorption, this material would absorb less of the total Lyman-$\alpha$ radiation, allowing material deeper in the wind to absorb some momentum. At best, however, the same total momentum would be absorbed e.g. in the up-orbit arm, with the possibility that some significant frequency bins would have little absorption. As the creation of a cometary tail is a bulk phenomenon, only the total momentum absorbed by the up-orbit arm will have a significant influence on its formation.

Note that this does not necessarily mean that smaller amounts of material cannot be accelerated to the observed velocities of $\sim 10-100$ km/s. Indeed, if we were to include the effects of Doppler shifting in the Lyman-$\alpha$ absorption, one can imagine that material on the far edge of the wind initially moving with some small velocity in the $x$ direction could be continually accelerated by photons of increasing Doppler velocity, allowing some material to reach the observed velocities. This requires that the accelerated material not have significant density in the direction of acceleration (as we have shown that the wind is collisional) and that material between it and the star is transparent to the necessary frequencies of Lyman-$\alpha$ radiation.

4.2 Comparison to Previous Work

The periodic behavior we observe in the high-flux case is qualitatively similar to the behavior found in the intermediate stellar wind regime of McCann et al. (2019) (see their Figure 12). They describe the regime as being characterized by a stellar wind pressure not capable of completely confining the planetary outflow, but strong enough to disrupt the outflow once the cross-sectional area between the outflow and stellar wind becomes large enough. We note that the radiation pressure from our high-flux case ($P_{Ly-\alpha}\approx 2.2\times 10^{-6}\mbox{ dyn/cm}^{2}$) is comparable to the pressure of McCann et al. (2019)’s intermediate stellar wind ($P_{\star}\approx 1.3\times 10^{-6}\mbox{ dyn/cm}^{2}$). Therefore, it is not surprising that the same qualitative behavior is found (though the planetary outflows do differ as discussed in 3.2.2). This lends credence to their suggestion that stellar radiation pressure and stellar winds should behave similarly.

It is also instructive to compare our simulations with those performed by Khodachenko et al. (2017), Bourrier & Lecavelier des Etangs (2013), Schneiter et al. (2016), and Cherenkov et al. (2018). Our results are in good agreement with the simulations of Khodachenko et al. (2017), who used a self-consistent 2-dimensional model of the wind launching due to incident EUV radiation. Those models included the effects of radiation pressure and charge exchange with the stellar wind, and found that the effects of radiation pressure on the Lyman-$\alpha$ observations were smaller than 1% in the most favorable cases. Their EUV flux of $1.75\times 10^{13}$ phot/s/cm2 and Lyman-$\alpha$ flux of $4.9\times 10^{14}\mbox{ phot/s/cm}^{2}$ are close to those used in our simulations of the intermediate-flux case (Run 2), making them a good point of comparison. Our results are also in good agreement with the similar simulations of Cherenkov et al. (2018), who resolved the planet at a moderately lower resolution than ours and the outer portions of the wind at a significantly lower resolution. They propagated their ionizing and Lyman-$\alpha$ radiation using a similar ray-tracing method as ours, and included the effects of the Doppler broadening in the Lyman-$\alpha$ line, and found that the wind envelope is essentially unaffected by Lyman-$\alpha$ flux below approximately $4\times 10^{16}\mbox{ phot/s/cm}^{2}$, which supports our assertion that our flux values represent lower bounds due to our ignoring the Doppler effects.

On the other hand, both Schneiter et al. (2016) and Bourrier & Lecavelier des Etangs (2013) find a cometary tail forming in their simulations. Schneiter et al. (2016) apply isotropic fixed boundary conditions at $3R_{p}$, and use a decreased value for the stellar gravity to approximate the effects of radiation pressure. Bourrier & Lecavelier des Etangs (2013) also use isotropic fixed boundary conditions rather than self-consistently calculated winds, with the wind launched at $2.8R_{p}$ for HD 209458b. Their radiation pressure is implemented as a reduction in the gravity of the star, but is calculated with the self-shielding one would expect due to the absorption of Lyman-$\alpha$ photons in the outer layers of the wind. Perhaps most significantly, Bourrier & Lecavelier des Etangs (2013) use a kinetic model of the planetary wind, with no pressure forces present, and remove hydrogen particles from their simulation once they’ve become ionized. We have seen that ionized hydrogen is a significant indirect sink for radiation pressure. Both Schneiter et al. (2016) and Bourrier & Lecavelier des Etangs (2013) include stellar winds, and Bourrier & Lecavelier des Etangs (2013) also include the effect of charge exchange. Thus it is difficult to perform direct comparisons to our simulations. Their results are similar to the semi-analytic models presented in Carroll-Nellenback et al. (2017), where a small Coriolis radius and significant radiation pressure (at least 10 percent of the stellar gravity) resulted in a cometary tail. The results from Schneiter et al. (2016) and Carroll-Nellenback et al. (2017) highlight the importance of considering self-shielding when simulating the effects of radiation pressure on the evaporating winds of exoplanets. A uniform reduction in gravity is evidently insufficient to capture the full behavior of the envelope. In addition, the continued presence of the stellar wind and, in Bourrier & Lecavelier des Etangs (2013), charge exchange makes it difficult to disentangle which effect is primarily responsible for the cometary tail and consequently the excess absorption in the blue wing of the Lyman-$\alpha$ line.

5 Conclusion

Based on these comparisons, we suggest that radiation pressure alone is insufficient to create a cometary tail in the outflow from HD 209458b. In addition, radiation pressure appears to be insufficient to create the observed absorption spectra of HD 209458b. To reproduce the observations, substantial amounts of gas must be accelerated to significantly greater than 50 km/s, which we do not see in our simulations. While the Lyman-$\alpha$ flux of HD 209458 is low in comparison to similar stars (Wood et al., 2005), it is not an order of magnitude lower. Thus even hot Jupiter host stars with more activity are likely to have insufficient radiation pressure to significantly affect their planets’ winds.

Both confinement by a stellar wind and charge exchange with the stellar wind are still potential sources of the fast neutral hydrogen observed in planetary winds, and a more thorough treatment of the Doppler shifting of the Lyman-$\alpha$ radiation may result in higher velocities due to radiation pressure. In fact, charge exchange has been shown in some cases to be effective in producing fast neutral hydrogen (e.g. Ekenbäck et al., 2010; Bourrier & Lecavelier des Etangs, 2013; Bourrier et al., 2016), though Bourrier & Lecavelier des Etangs (2013) suggest that charge exchange is not a significant contributor to those results. We will further investigate the effects of these phenomena in future studies.

6 Acknowledgments

We thank the Other Worlds Laboratory (OWL) at University of California, Santa Cruz for facilitating this collaboration by way of the OWL Exoplanets Summer Program, funded by the Heising-Simons Foundation. This work used the computational and visualization resources in the Center for Integrated Research Computing (CIRC) at the University of Rochester and the computational resources of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, provided through allocation TG-AST120060 from the Extreme Science and Engineering Discovery Environment (XSEDE) (Towns et al., 2014), which is supported by National Science Foundation grant number ACI-1548562. Financial support for this project was provided by the Department of Energy grant DE-SC0001063, the National Science Foundation grants AST-1515648, AST-1813298, and AST-1411536, and the Space Telescope Science Institute grant HST-AR-12832.01-A. EB acknowledges additional support from KITP UC Santa Barbara, funded by NSF Grant PHY-1748958, and the Aspen Center for Physics, funded by NSF Grant PHY-1607611.