# Using Lyman-$\alpha$ transits to constrain models of atmospheric escape

Ethan Schreyer,1,3 James. E. Owen,1 R. O. Parke Loyd2 and Ruth Murray-Clay3
1Astrophysics Group, Imperial College London, Blackett Laboratory, Prince
Consort Road, London SW7 2AZ, UK
2Eureka Scientific, Inc., Oakland, CA 94602, USA
3Department of Astronomy and Astrophysics, University of California, Santa
Cruz, CA 95064, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Lyman-$\alpha$ transits provide an opportunity to test models of atmospheric
escape directly. However, translating observations into constraints on the
properties of the escaping atmosphere is challenging. The major reason for
this is that the observable parts of the outflow often comes from material
outside the planet’s Hill sphere, where the interaction between the planetary
outflow and circumstellar environment is important. As a result, 3D models are
required to match observations. Whilst 3D hydrodynamic simulations are able to
match observational features qualitatively, they are too computationally
expensive to perform a statistical retrieval of properties of the outflow.
Here, we develop a model that determines the trajectory, ionization state and
3D geometry of the outflow as a function of its properties and system
parameters. We then couple this model to a ray tracing routine in order to
produce synthetic transits. We demonstrate the validity of this approach,
reproducing the trajectory of the outflows seen in 3D simulations. We
illustrate the use of this model by performing a retrieval on the transit
spectrum of GJ 436 b. Our model constrains the sound speed of the outflow
$\gtrsim 10\text{ km s}^{-1}$, indicating that we can rule out core-powered
mass loss as the mechanism driving the outflow for this planet. The bound on
planetary outflow velocity and mass loss rates are consistent with a
photoevaporative wind.

###### keywords:

planets and satellites: atmospheres – planet-star interactions

††pubyear: 2022††pagerange: Using Lyman-$\alpha$ transits to constrain models
of atmospheric escape–B

##  1 Introduction

Many of the discovered exoplanets reside close to their host stars. These
planets are highly irradiated, causing their upper atmospheres to be extremely
hot. The resulting pressure gradients can drive the gas to escape
hydrodynamically. For small, low-mass planets ($1-4\,{\rm R}_{\oplus}$,
$\lesssim 20$ M⊕), the integrated effect of this mass loss can completely
strip them of their primordial hydrogen/helium envelopes (Valencia et al.,
2010; Owen & Wu, 2013; Lopez & Fortney, 2013; Chen & Rogers, 2016). This
process can explain the observed lack of Neptune-sized planets at short
orbital periods, the "hot Neptune desert" (Szabó & Kiss, 2011; Lundkvist et
al., 2016), and the observed sparsity of planets with radii $\sim
1.8R_{\oplus}$, referred to as "the radius gap" (Fulton et al., 2017; Van
Eylen et al., 2018). However, other interpretations characterize the
demographics of close-in, low-mass exoplanets as a result of gas-accretion in
a hot inner disc (e.g, Lee & Connors, 2021; Lee et al., 2022), or separate
populations of terrestrial planets and water worlds (e.g, Zeng et al., 2019;
Venturini et al., 2020; Luque & Pallé, 2022; Izidoro et al., 2022).

Atmospheric mass loss alters the observed properties of the close-in, low-mass
exoplanet population. Thus, models of hydrodynamic escape are needed to infer
the properties of the primordial population, which can be used to constrain
planet formation models. For low-mass planets, the two most prevalent
mechanisms of hydrodynamic escape are photoevaporation (e.g., Yelle, 2004;
Muñoz, 2007; Murray-Clay et al., 2009; Owen & Jackson, 2012; Kubyshkina et
al., 2018; Caldiroli et al., 2021) and core-powered mass loss (e.g., Ginzburg
et al., 2018; Gupta & Schlichting, 2019). In the photoevaporation model,
incident X-ray/EUV photons photoionize the upper atmosphere, heating it to up
to $\sim 10^{4}$ K. The resulting thermal pressure drives a transonic
hydrodynamic outflow, causing the gas to escape the planet at velocities of
the order of the sound speed $\sim$ 7-20 km/s, depending on the composition,
planet gravity, and radiative properties of the flow. In the core-powered mass
loss model, the internal luminosity of the planet and bolometric flux of the
star drive escape at velocities $\sim$ 1-2 km/s.

To test whether these mechanisms can reproduce the observed radius-period
distribution of planets, evolutionary models of mass loss were fit to the
current planet population for both photoevaporation (e.g., Owen & Wu, 2017;
Jin & Mordasini, 2018; Wu, 2019; Rogers & Owen, 2021; Rogers et al., 2023) and
core-powered mass loss (Gupta & Schlichting, 2019, 2020; Gupta et al., 2022).
Both mechanisms can reproduce the observed radius gap, and further analysis
showed the dominant mechanism in sculpting the radius cannot be
observationally distinguished from the current sample of planets (Loyd et al.,
2020; Rogers et al., 2021; Ho & Van Eylen, 2023). However, these mechanisms do
not infer the same primordial core-mass distribution; therefore, determining
the dominant mass-loss mechanism is of significant importance. Owen &
Schlichting (2023) analytically studied the transition between core-powered
and photoevaporative outflows and applied their results to the observed super-
Earth population. Although they found that each mechanism was responsible for
the removal of the gaseous envelopes of a significant fraction of the super-
Earth population, they argued photoevaporation is responsible for the final
carving of the radius valley.

Atmospheric escape also plays an important role in the evolution of the
atmosphere’s composition. Clearly, stripping the planet of its hydrogen
envelope greatly changes its composition. More subtle, though, is the effect
this process can have on the inventory of heavier elements in the atmosphere.
Generally, heavier elements are less susceptible to escape as the temperatures
required for them to reach the escape velocity are significantly higher than
that of hydrogen. However, a mechanism to remove these elements from the
atmosphere is to drag them out in a hydrogen-dominated outflow through
collisions. This process has been modelled in the case of early Earth and
Venus (e.g., Kasting & Pollack, 1983; Kasting et al., 1984) and close-in
exoplanets (e.g Luger & Barnes, 2015; Hu et al., 2015; Malsky & Rogers, 2020;
Malsky et al., 2023). Whether heavy elements can be dragged out depends on the
escape rate, and therefore, different escape rates can lead to varied final
atmospheric compositions.

An obvious way to test hydrodynamic escape models is via direct observations
of atmospheric escape using transmission spectroscopy. The two most commonly
used lines are Lyman-$\alpha$ (e.g., Vidal-Madjar et al., 2003; Lecavelier Des
Etangs et al., 2010; Ehrenreich et al., 2015; Bourrier et al., 2017; dos
Santos et al., 2020) and HeI 10830 Å (e.g., Spake et al., 2018; Allart et al.,
2018; Kirk et al., 2020; Zhang et al., 2023), although, H-$\alpha$ and some
metal lines have also been used to detect escaping gas (e.g., Vidal-Madjar et
al., 2004; Sing et al., 2019; Ben-Jaffel et al., 2022). However, translating
these observations into constraints on escape models is challenging due to a
lack of theoretical work.

Specifically for Lyman-$\alpha$ observations, the absorption by neutral
hydrogen in the interstellar medium obscures the line core $\sim-50\text{ to
}50\text{ km s}^{-1}$ (unless the star has a particularly large radial
velocity) (e.g., Landsman & Simon, 1993), such that in order interpret these
observations we rely on information contained in the wings of the spectrum. As
typical escape models predict outflow velocities to be < 50 km s-1, these
observations do not probe the launching regions of the outflow. The fact that
transits are indeed observed at such high velocities (generally in the blue
wing) indicates that interactions with the circumstellar environment have an
important role in making the planetary outflow observable. Furthermore,
observations often have a transit so deep it implies a size larger than the
planet’s Hill sphere and a duration much greater than the optical transit,
indicating that a significant absorption comes from parts of the outflow far
away (and unbound) from the planet. At these distances, interactions with the
circumstellar environment are important (e.g., Bourrier & Lecavelier des
Etangs, 2013; Bisikalo et al., 2013; Matsakos et al., 2015; McCann et al.,
2019).

Three mechanisms have been proposed to create high-velocity neutral hydrogen
atoms in the outflow: 1) ram-pressure from the stellar wind can accelerate the
outflowing gas away from the star; 2) Lyman-$\alpha$ radiation pressure can do
the same and 3) charge exchange between stellar wind protons and neutral
hydrogen in the outflow can create a population of hot, fast neutral atoms
known as energetic neutral atoms (ENAs) (e.g., Holmström et al., 2008;
Tremblin & Chiang, 2012). There have been a number of 3D particle (e.g.,
Bourrier & Lecavelier des Etangs, 2013; Bourrier et al., 2015, 2016), and
hydrodynamic/radiation-hydrodynamic (e.g., Bisikalo et al., 2013; Matsakos et
al., 2015; Carroll-Nellenback et al., 2016; Khodachenko et al., 2019; McCann
et al., 2019; Debrecht et al., 2020; Carolan et al., 2021a; Hazra et al.,
2022; MacLeod & Oklopčić, 2022) simulations focusing on the interaction
between the planetary outflow and the circumstellar environment. These works
show that the planetary outflow is shaped into a cometary tail by the stellar
tidal field, stellar wind ram pressure and radiation pressure. Typically,
these simulations are able to qualitatively reproduce observational features
such as the large blue-shifted transits and extended post-transit absorption,
indicating that they are capturing the basic physics driving these
observations. However, an unavoidable problem with 3D simulations is that they
are computationally expensive, so they are only able to be tested over a
limited parameter range. Consequently, it is challenging to match the transit
spectrum produced by these simulations to observations. Therefore, it is
difficult to use observations to test atmospheric models or retrieve outflow
properties rigorously. Furthermore, it is unclear whether the inability to
match certain observations arises from evaluating models in the wrong part of
parameter space or a more fundamental issue with the general framework.

In an attempt to gain a greater understanding of how to solve this problem,
Owen et al. (2023) developed a semi-analytic model of the outflowing tail that
could predict how Lyman-$\alpha$ transit observations would vary over a large
range of possible planetary and stellar parameters. With this model, they
highlighted the counter-intuitive nature of Lyman-$\alpha$ observations:
stronger EUV irradiation leads to larger mass-loss rates but often weaker
transits both in-depth and duration, something previously noted in simulations
(e.g., Bourrier et al., 2013, 2016; Villarreal D’Angelo et al., 2021). This
also provides a simple explanation of why several highly irradiated sub-
Neptunes do not have a detectable transit: HD97658b (Bourrier et al., 2017),
$\pi$ Men c (García Muñoz et al., 2020) and K2-25b (Rockcliffe et al., 2021).

Although this model gives useful insight into how to interpret Lyman-$\alpha$
transits, it simplified the problem to the extent that it cannot be used to
fit observations. One problem is that they assumed that the tail remains on
the same orbital separation of the planet, which is not the case for very
close-in planets like GJ 436 b (e.g. Khodachenko et al., 2019). Furthermore,
they assume the density of the gas in the tail is constant; however, by
conservation of mass, the density should decrease as the gas is accelerated.
Finally, the geometry of the tail (i.e. the vertical and horizontal extent of
the gas) was calculated by ignoring the stellar wind, which could compress the
tail.

In this work, we build a self-consistent model of the outflowing tail that
attempts to accurately reproduce the 3D hydrodynamics simulations that it is
based on. We prove the validity of this approach by comparing the trajectory
of the tails our model produces to multiple 3D hydrodynamics simulations. The
primary motivation of the model is to produce synthetic Lyman-$\alpha$
transits that can be computed quickly enough to perform “retrievals” of
atmospheric escape properties on the velocity-resolved transit spectrum. We
present a study of GJ 436 b, a test case and discuss interpretations of the
results.

##  2 Outflow Model

In this section, we build a model of the outflowing planetary gas in a steady
state. We split the outflow into two components: inside and outside the Hill
sphere of the planet. Inside the Hill sphere, the outflow is approximately
spherically symmetric (e.g. Murray-Clay et al., 2009); thus, we model the flow
as a 1D Parker wind (Parker, 1958). Outside the Hill sphere, 3D hydrodynamic
simulations show that planetary outflow is sculpted into a cometary tail by
the stellar tidal field and stellar wind ram pressure (e.g., Matsakos et al.,
2015; Carroll-Nellenback et al., 2016; Khodachenko et al., 2019; McCann et
al., 2019; Debrecht et al., 2020; Carolan et al., 2021a; Hazra et al., 2022;
MacLeod & Oklopčić, 2022). Inspired by methods used to study disc winds
(semi-) analytically (Fukue & Okada, 1990; Clarke & Alexander, 2016), we model
the cometary tail as a streambundle, solving for its dynamics. The transition
from a spherically symmetric outflow to a cometary tail cannot be modelled
semi-analytically. However, the primary goal of this work is to understand the
flow far outside the Hill sphere of the planet where a cometary tail is a good
approximation. This is because Lyman-$\alpha$ transits primarily probe
material outside the planet’s Hill sphere. Therefore, we choose to construct
the full planetary outflow by “glueing” these two outflow models to each other
at the Hill sphere of the planet, which we describe in Section 2.4. While the
transiting Hill sphere region often has a small contribution to the transit
profile, this approach allows us to have the outflow into the tail come from a
physically motivated (although parameterised) outflow model.

This model excludes parts of the outflow leading the planet, and therefore, we
cannot model the pre-transit absorption. We note that pre-transit absorption
may be highly variable because the gas leading the planet is unstable for
moderate stellar wind strengths (e.g., McCann et al., 2019), which makes it
challenging to model in this framework.

###  2.1 Model of Cometary Tail

In a frame co-rotating with the planet, we model the cometary tail as a
streamline bundle (e.g., Clarke & Alexander, 2016) with an elliptical cross-
section. This is composed of a curve representing the trajectory of the tail,
with ellipses stacked along it of variable height and depth. The variable
($s$) denotes the distance along the curve. The depth of the ellipse is small
compared to the radius of curvature of the tail, and therefore, the tail does
not intersect itself. At each point along this curve, we define a right-
handed, local Cartesian coordinate system $(l,a,\zeta)$ with the origin at the
ellipse’s centre. We orient this coordinate system such that
$\boldsymbol{\hat{l}},\boldsymbol{\hat{a}},\boldsymbol{\hat{\zeta}}$ point
tangent to the trajectory of the tail, normal to the trajectory of the tail
and parallel to the orbital angular momentum of the planet, respectively. As
the tail does not intersect with itself, any point can be specified by a
distance along the tail $s$, and coordinates $(a,\zeta)$ in the corresponding
ellipse.

We also define a global right-handed Cartesian coordinate system $(x,y,z)$ co-
rotating with the planet, with its origin at the star. The radial distance
from the star is denoted by $r$. We choose the $(x,y)$ coordinates to lie in
the planet’s orbital plane so that the trajectory of the tail also lies in the
orbital plane (i.e. $s=s(x,y)$). Given a point along the tail
$s_{0}=(x_{0},y_{0})$, a point in the corresponding cross-sectional plane is
given by:

$\displaystyle x=x_{0}-a\text{sin}\tau$ (1) $\displaystyle
y=y_{0}+a\text{cos}\tau$ (2) $\displaystyle z=\zeta$ (3)

where $\tau$ is the angle between $\boldsymbol{\hat{l}}$ and
$\hat{\textbf{x}}$ measured anticlockwise. Figure 1 shows a sketch of the
model and coordinate systems. For simplicity, we consider both the velocity
and the ionization fraction of the gas in the tail to be a function of the
distance $s$ along the tail only. These simplifications allow us to separately
model the bulk motion (i.e. the trajectory), ionisation state, geometry and
density structure of the outflow. The model, therefore, consists of three
parts:

(i) Self-consistently solving for the trajectory of the tail in the planet’s
orbital plane by assuming a generalised geostrophic balance.

(ii) Solving for the neutral fraction of hydrogen atoms along the tail using
the optically thin ionisation equation.

(iii) Solving for the height, depth and density structure of this
streambundle, assuming the gas is in hydrostatic equilibrium in the cross-
section of the streambundle.

Modelling the tail as a streambundle implicitly treats the planetary outflow
as a fluid as opposed to the particle approach used in some simulations of the
outflow (e.g., Ehrenreich et al., 2015; Bourrier et al., 2015). For many
planets we are studying, the density is high enough that the gas is
collisional throughout the outflow. This is illustrated by 3D simulations
(e.g., McCann et al., 2019; Debrecht et al., 2020) and therefore, we believe
this approach is valid.

Figure 1: A schematic diagram of the planetary outflow that is sculpted into a
cometary tail by the stellar tidal field and stellar wind ram pressure. The
light blue curve represents the trajectory of the tail, which is parameterised
by a coordinate ($s$), representing the distance along this curve. At each
point along this curve, we define a local cartesian coordinate system ($l$,
$a$, $\zeta$) that can be related to a global cartesian coordinate system
($x$, $y$, $z$) centred on the star. We solve for i) the trajectory, ii) the
ionization state, and iii) the geometry of the tail as a function of system
parameters shown in Table 2.

####  2.1.1 Trajectory of the Tail

In a reference frame co-rotating with the planet, in which we treat the flow
as steady, we solve for the trajectory of the tail as a 1D streamline subject
to momentum conservation.

$\displaystyle(\textbf{u}\cdot\boldsymbol{\nabla})\textbf{u}=\textbf{a}_{\text{sw}}-\boldsymbol{\nabla}\phi-2\boldsymbol{\Omega}\times\textbf{u}$
(5)

where $\textbf{a}_{\text{sw}}$ is the acceleration due to ram pressure of the
stellar wind, $\boldsymbol{\nabla}\phi$ is effective potential, and
$2\boldsymbol{\Omega}\times\textbf{u}$ is the Coriolis force. We assume the
planet has a circular orbit and ignore the Euler force (it’s not possible to
construct a steady-state model of an eccentric orbit). Although the stellar
wind ram pressure force is a surface force, we approximate that this force is
applied like a body force, as pressure gradients quickly distribute this force
across the tail. Furthermore, geostrophic balance keeps the tail stable. We
discuss the validity of this approximation later in this section. In this
framework, the force per unit length that the stellar wind applies to the tail
is:

$\displaystyle\frac{d\textbf{p}}{dt}=2H\rho_{*}(\textbf{u}_{*}-\textbf{u}\,\text{cos}\chi)^{2}\text{sin}\chi\>\hat{\textbf{u}}_{*}$
(6)

where $H$ is the height of the tail, $\rho_{*}$ and $\textbf{u}_{*}$ are the
stellar wind density and velocity, respectively, and $\chi$ is the angle
between the tangent of the tail and the direction stellar wind. In the case
where the stellar wind is strong, simulations show that all the mass lost from
the planet is directed down the cometary tail (e.g., McCann et al., 2019).
Therefore, we can define the mass per unit length in the tail:

$\displaystyle\lambda=\frac{\dot{M}_{p}}{u}$ (7)

where $\dot{M}_{p}$ is the mass loss rate of the planet and $u=|\textbf{u}|$.
Therefore, the acceleration of the tail due to stellar wind ram pressure is
given by:

$\displaystyle\textbf{a}_{\text{sw}}=\frac{2Hu\rho_{*}(\textbf{u}_{*}-\textbf{u}\,\text{cos}\chi)^{2}\text{sin}\chi}{\dot{M}_{p}}\hat{\textbf{u}}_{*}$
(8)

For this framework to be valid, the timescale for pressure waves to distribute
this force across the tail must be shorter than the timescale in which the
forces acting on the tail significantly change. We estimate this by comparing
the acceleration timescale to the sound crossing time as:

$\displaystyle\frac{\tau_{sc}}{\tau_{a}}=\frac{\text{dlog}u}{\text{dt}}\frac{2D}{c_{s}}$
(9)

where $c_{s}$ is the isothermal sound speed of the gas in the outflow and $D$
is the depth of the tail. Taking the conservative assumption that $u_{*}\gg u$
and $\text{sin}\chi=1$, so that
$a_{\text{sw}}\sim(2Hu_{*}^{2}\rho_{*})/\dot{M}_{p}$, this equals

$\displaystyle\frac{\tau_{sc}}{\tau_{a}}\lesssim\frac{1}{\pi}\left(\frac{\dot{M}_{*}}{\dot{M}_{p}}\right)\left(\frac{u_{*}}{c_{s}}\right)\left(\frac{HD}{r^{2}}\right)$
(10)

For a typical sub-Neptune around a solar-like star at 0.03 AU, this is smaller
than one.

$\displaystyle\begin{split}\frac{\tau_{sc}}{\tau_{a}}\lesssim\
&\left(\frac{\dot{M}_{*}}{10^{12}\text{ g
s}^{-1}}\right)\left(\frac{\dot{M}_{p}}{10^{10}\text{ g
s}^{-1}}\right)^{-1}\left(\frac{u_{*}}{100\text{ km s}^{-1}}\right)\\\
&\times\left(\frac{c_{s}}{10\text{ km
s}^{-1}}\right)^{-1}\left(\frac{HD}{3\times 10^{-3}r^{2}}\right)\end{split}$
(11)

where $H$ and $D$ are conservatively estimated using the distance that a
ballistic particle travels after being launched from the planet at 10 km s-1.

As this estimate is an upper limit, it is reasonable to approximate the total
effect of stellar wind ram pressure to apply a force equal to the momentum
flux throughout the tail. Additionally, in this framework, the height and
depth scales like the bow shock radius. Therefore, increasing the stellar wind
mass loss rate or velocity should reduce the height and depth equivalently so
that the argument still holds.

We also test the validity of this approach by comparing the tail trajectories
that we have calculated in our model to multiple 3D hydrodynamic simulations
(see Section 3 and Figure 2). Detailed inspection of publicly available
simulation outputs demonstrates that we are able to replicate the trajectories
of these simulations using the same system parameters, indicating that our
assumptions are reasonable.

####  2.1.2 Geometry of the Tail

A corollary of the previous argument is that the gas inside the tail should
adjust to its local conditions. Therefore, in the elliptical cross-section of
the tail, we can approximate the gas to be in hydrostatic equilibrium.
Perpendicular to $\boldsymbol{\hat{l}}$, the forces on a parcel of gas are:

$\displaystyle-\frac{1}{\rho(a,z)}\boldsymbol{\nabla}{P(a,z)}-\boldsymbol{\nabla}\phi(a,z)-2\boldsymbol{\Omega}\times\textbf{u}+\frac{\textbf{u}^{2}}{r_{c}(a)}\hat{\textbf{a}}=0$
(12)

where $\rho$ and $P$ are the density and pressure, respectively, and $r_{c}$
is the radius of curvature corresponding to streamlines in the outflow. We
note that this equation is also an implicit function of the distance along the
tail $s=(x,y)$. The effective potential is:

$\displaystyle\begin{split}\phi=-&\frac{GM_{*}}{\sqrt{(x-a\text{sin}\tau)^{2}+(y+a\text{cos}\tau)^{2}+z^{2}}}\\\
&-\frac{1}{2}\Omega^{2}\left((x-a\text{sin}\tau)^{2}+(y+a\text{cos}\tau)^{2}\right)\end{split}$
(13)

Since the height and depth of the tail are small compared to the distance from
the star, we expand the effective potential and radius of curvature to second-
order in terms of $\left(\frac{a}{r},\frac{z}{r}\right)$. Up to this order,
the horizontal ($a$) and vertical ($z$) coordinates decouple, therefore, the
density is separable, and each direction can be solved independently. For
simplicity, we solve for the density using an isothermal equation of state
such that $P=c_{s}^{2}\rho$. In the vertical direction, the density
distribution is given by

$\displaystyle-\frac{c_{s}^{2}}{\rho(z)}\frac{\partial\rho(z)}{\partial
z}-\frac{GMz}{r^{3}}=0$ (14)

which can be analytically solved as

$\displaystyle\rho(z)=\rho_{a,0}\text{exp}\left[-\frac{GMz^{2}}{2c_{s}^{2}r^{3}}\right]=\rho_{a,0}\text{exp}\left[-\frac{\Omega^{2}_{k}z^{2}}{2c_{s}^{2}}\right]$
(15)

where $\rho_{a,0}$ is the density in the orbital plane at horizontal
coordinate $a$ and $\Omega_{k}$ is the Keplerian angular velocity. This is
just the standard solution for vertical hydrostatic equilibrium of a rotating
fluid (e.g. in a disc).

Following the same procedure to calculate the horizontal density distribution
leads to unphysical solutions. This is a result of our assumption that the
velocity of the gas is constant across the cross-sectional slice of the tail.
We instead approximate the depth by assuming that the angular momentum
difference between the edge of the tail and the centre stays constant, where
the initial dispersion is approximated by launching particles at the sound
speed with varying angular momenta from the planet. In the limit where the
orbital velocity of the planet is greater than the velocity of the tail, the
depth
$\sim\frac{\Omega_{k}}{c_{s}}\sqrt{\frac{r}{a_{p}}}\sim\frac{\Omega_{k}}{c_{s}}$.
Therefore, a sensible choice for a density distribution is:

$\displaystyle\rho(a)\propto\text{exp}\left[-\frac{\Omega_{k}^{2}a^{2}}{c_{s}^{2}}\right]$
(16)

We note that in order to calculate the optical depth of the tail in transit,
we integrate this density profile (almost) horizontally through the tail, and
therefore, the transit depth will only be weakly dependent on this choice.
Given this choice, the total density structure is given by:

$\displaystyle\rho=\rho_{0}\text{exp}\left[-\frac{a^{2}}{\alpha^{2}}-\frac{z^{2}}{\beta^{2}}\right],\quad\alpha=\frac{\Omega_{k}}{c_{s}},\quad\beta=\frac{\Omega_{k}}{\sqrt{2}c_{s}}$
(17)

where $\rho_{0}$ is the density at the centre of the ellipse. In order to find
the height, depth and $\rho_{0}$ of the tail, three constraints must be
applied. The first of these comes from equating the mass-flux to the integral
of the density over the elliptical cross-section. The latter two are given by
specifying conditions at the boundary between the planetary and stellar wind.
As the stellar wind is incident on the inner edge of the tail, conditions vary
substantially on either side of the tail. We choose to evaluate this at the
inner edge and top of the ellipse (we could equivalently define it a the
bottom), where these boundaries are well-defined. These three constraints can
be written as:

$\displaystyle\text{i)}\
\iint_{\mathcal{A}}\rho(a,z)\,dadz=\frac{\dot{M}_{p}}{u}\quad\text{ (Mass per
unit length)}$ (18) $\displaystyle\text{ii)}\ \
\rho(-D,0)c_{s}^{2}=P_{\text{sw},D}\quad\text{ (Pressure at inner edge)}$ (19)
$\displaystyle\text{iii)}\ \rho(0,H)c_{s}^{2}=P_{\text{sw},H}\quad\text{
(Pressure at upper edge)}$ (20)

where $\mathcal{A}$ is the area of the ellipse and $P_{\text{sw},D}$ and
$P_{\text{sw},H}$ is the stellar wind pressure at the inner edge and top of
the tail respectively. We note that the obtained stellar wind pressure should
correspond to the shocked stellar wind. In the orbital plane, we approximate
the shock front to be parallel to the trajectory of the tail, and in the
vertical plane, it forms a bow shock. Since simulations indicate that the bow
shock is not thin, there are no general analytic ways to calculate the shape
of a bow shock or the downstream pressure. However, it is possible to solve
for the pressure at the nose of the shock using the oblique shock equations.
For simplicity, we assume that the pressure at the inner edge and top of the
tail are equal and are a linear scaling of the pressure at the nose.

$\displaystyle\kappa P_{\text{nose}}=P_{\text{sw},D}=P_{\text{sw},H}$ (21)

As we expect the pressure to decrease as the gas moves away from the nose,
$\kappa<1$. We set $\kappa=0.3$ using Schulreich & Breitschwerdt (2011) as a
guide. With this simplification, using Equations 17, 19, 20, the depth and
height are related as:

$\displaystyle\frac{D}{\alpha}=\frac{H}{\beta}$ (22)

The integral in Equation 18 cannot generally be solved analytically. However,
it can be written in terms of special functions (e.g Waugh, 1961).

$\displaystyle\frac{1}{2\pi}\iint_{\mathcal{A}}e^{-\frac{x^{2}+y^{2}}{2}}dxdy=p\left(\frac{A+B}{2},\frac{A-B}{2}\right)-p\left(\frac{A-B}{2},\frac{A+B}{2}\right)$
(23)

$\displaystyle
p(R,r)=e^{-\frac{r^{2}}{2}}\int_{0}^{R}e^{-\frac{t^{2}}{2}}I_{0}(rt)tdt$ (24)

where $A$, $B$ are the semi-major and semi-minor axes of the ellipse,
$\mathcal{A}$ is the area of the ellipse and $I_{0}$ is the modified Bessel
function of the first kind. In the special case given by Equation 22, this can
be solved exactly to give:

$\displaystyle D=\alpha\text{ln}\left(\frac{\dot{M}_{p}c_{s}^{2}}{u\pi\gamma
P_{\text{nose}}\alpha\beta}+1\right)$ (25) $\displaystyle
H=\frac{\beta}{\alpha}D$ (26)

####  2.1.3 Ionization of Tail

As absorption of Lyman-$\alpha$ radiation is from neutral hydrogen, it is
essential to calculate the ionization state of hydrogen in the tail. In the
photoevaporation model, where the absorption of EUV photons launches the
planetary outflow, the tail must be optically thin to EUV photons; otherwise,
the EUV flux absorbed at the base of the flow would be too small to drive a
large outflow (Owen et al., 2023). Owen et al. (2023) showed that this extends
to core-powered mass loss outflows, which are not driven by EUV radiation. As
a result, the photoionization rate of neutral hydrogen is approximately
constant in each cross-sectional slice of the tail (also noting that $D,H\ll
r$).

We also include hydrogen recombination in our model. For old small planets,
whilst observable, the gas in the tail is far way from photoionization-
recombination equilibrium, and therefore, the effect of recombination is
minimal. However, for giant planets or young planets with high mass loss
rates, the outflows may be dense enough for recombination to be important. As
recombination is density-dependent, it varies down the tail. For simplicity,
we approximate that the recombination rate of neutral hydrogen is constant in
each cross-sectional slice of the tail and equivalent to the mean density in
that slice. The number density, $n$, of protons and hydrogen atoms in a cross-
sectional slice of the tail is:

$\displaystyle n=\frac{\dot{M}_{p}}{\pi uHDm_{H}}$ (27)

where $m_{H}$ is the mass of a hydrogen atom. As such, we can simplify the
ionization fraction to a 1D problem, which depends on distance $(s)$ down the
tail. In this limit, the ionization equation as a parcel of gas moves down the
tail can be written as:

$\displaystyle\frac{DX}{Dt}=(1-X)\Gamma-nX^{2}\alpha_{A}$ (28)

where X is the ionisation fraction, $\Gamma$ is the photoionisation rate and
$\alpha_{A}$ is the case A recombination coefficient. We use the case A
recombination coefficient because the tail is optically thin to recombination
photons. This is because photons with energies $>13.6\text{ eV}$, that are
produced by a hydrogen atom recombining directly to the ground state generally
escape the tail rather than ionizing another hydrogen atom.

The optically thin photoionization rate is given by:

$\displaystyle\Gamma=\int_{13.6~{}{\rm eV}}^{\infty}\frac{L_{\nu}}{4\pi
r^{2}h\nu}\sigma_{\nu}d\nu$ (29)

where $\nu$ is the frequency of a photon, $\sigma_{\nu}$ is the
photoionization cross section of hydrogen as a function of frequency, and
$L_{\nu}$ is the luminosity of the star as a function of frequency. To map the
optically thin photoionization rate to an EUV luminosity, one needs to use a
stellar spectrum (see Appendix B). However, as the optically thin
photoionization rate is the physical quantity that determines the ionization
structure in the tail, we choose to use this quantity throughout this work.

Thus, in steady state, the governing equation for the ionization fraction
becomes:

$\displaystyle\frac{dN}{ds}=-\frac{\Gamma
N}{u}+\frac{n(1-N)^{2}\alpha_{A}}{u}$ (30)

where $N=1-X$ is the neutral fraction of hydrogen.

####  2.1.4 Stellar Wind

In order to calculate the acceleration of the tail due to the stellar wind, we
must specify the velocity and density of the stellar wind. For simplicity, we
treat the stellar wind as spherically symmetric, such that we can link the
density, velocity, and mass loss rate of the star through mass conservation.

$\displaystyle\dot{M}_{*}=4\pi r^{2}\rho_{*}u_{*}$ (31)

where $\dot{M}_{*}$ is the mass loss rate of the star. In this exploratory
work, we also choose the stellar wind to impact the tail radially. We note
that the outflow model, along with the simulations it is based on, is only
valid for planets located outside the Alfvén surface of the star where the ram
pressure of wind dominates over the magnetic pressure of the star. At several
stellar radii, the stellar wind is often supersonic, and therefore, its
velocity increases logarithmically as a function of radial distance r.
Therefore, we approximate that the stellar wind velocity is constant across
the radial extent of the tail. As a result, the ram pressure of the stellar
wind is:

$\displaystyle\rho_{*}u_{*}^{2}=\frac{\dot{M}_{*}u_{*}}{4\pi r^{2}}$ (32)

Although we have ignored many of the intricacies of the stellar wind, we
highlight that Lyman-$\alpha$ transits are probing the stellar wind conditions
near the planet. Therefore, a full 3D treatment of the stellar wind is
unnecessary. A possible complication with this arises when stellar wind
conditions can vary considerably along the orbit of the planet. However, more
complicated stellar wind conditions could be added to the model relatively
simply.

####  2.1.5 Method of Solution

It appears that we have to solve a system of coupled partial differential
equations for the trajectory, ionization fraction, and geometry of the tail.
However, as the domain of all the functions in these equations is the 1D
trajectory curve, they can be written as a function of distance along the tail
$(s)$. This system can equivalently be written as a set of coupled ordinary
differential equations:

$\displaystyle
u\frac{du_{x}}{ds}=(\textbf{a}_{\text{sw}})_{x}-(\boldsymbol{\nabla}\phi)_{x}+2\Omega
u_{y}$ (33) $\displaystyle
u\frac{du_{y}}{ds}=(\textbf{a}_{\text{sw}})_{y}-(\boldsymbol{\nabla}\phi)_{y}-2\Omega
u_{x}$ (34) $\displaystyle\frac{dx}{ds}=\frac{u_{x}}{u}$ (35)
$\displaystyle\frac{dy}{ds}=\frac{u_{y}}{u}$ (36)
$\displaystyle\frac{dN}{ds}=-\frac{\Gamma
N}{u}+\frac{n(1-N)^{2}\alpha_{A}}{u}$ (37) $\displaystyle
D=\alpha\text{ln}\left(\frac{\dot{M}_{p}c_{s}^{2}}{u\pi
P_{\text{sw},D}\alpha\beta}+1\right)$ (38) $\displaystyle
H=\frac{\beta}{\alpha}D$ (39)

where the notation $(\textbf{A})_{i}$ represents the component of the vector A
in the $i$ direction. We solve these implicitly using an implementation of a
backward differentiation formula method BDF (Byrne & Hindmarsh, 1975) provided
in the scipy library (Virtanen et al., 2020). The relative and absolute
tolerance used were $1\times 10^{-13}$ and $1\times 10^{-14}$, respectively.

###  2.2 Spherical Wind

Inside the Hill sphere, the outflow is modelled as a 1D constant sound speed
Parker wind, with the addition of tidal gravity, along a streamline
originating from the sub-stellar point on the planet. We do this in order to
stay agnostic about the mechanism driving the outflow. The outflow is
parameterised by a mass loss rate and a sound speed that is constant
throughout the flow. For a given mass loss rate and sound speed, the velocity
of the outflowing gas is analytically given by:

$\displaystyle u=\begin{cases}-W_{0}[-D(r^{\prime})]\quad\ r^{\prime}\leq
r_{cr}\\\ -W_{-1}[-D(r^{\prime})]\quad r^{\prime}>r_{cr}\end{cases}$ (40)

$\displaystyle
D(r^{\prime})=\left(\frac{r^{\prime}}{r_{cr}}\right)^{-4}\text{exp}\left[4r_{\alpha}\left(\frac{1}{r_{cr}}-\frac{1}{r^{\prime}}\right)+\frac{2r_{\alpha}}{R_{H}^{3}}(r_{cr}^{2}-r^{\prime
2})-1\right]$ (41)

where $r^{\prime}$ is the radial distance from the planet, $W_{k}$ is the
Lambert W function on branch $k$ (Cranmer, 2004), $r_{cr}$ is the sonic point
of the flow, $r_{\alpha}$ is the sonic point for the flow without tidal
gravity and $R_{H}$ is the Hill radius of the planet. We calculate the
fraction of ionized hydrogen by analogy with our method for the ionizing
fraction along the tail. Thus we use:

$\displaystyle\frac{\partial N}{\partial r^{\prime}}=-\frac{\Gamma
e^{-\tau_{\bar{\nu}}(r^{\prime})}N}{u}+\frac{\dot{M}_{p}(1-N)^{2}\alpha_{A}}{4\pi
r^{\prime 2}u^{2}m_{p}}$ (42)

We include the optical depth of the outflow to $20$ eV photons,
$\tau_{\bar{\nu}}$, in the photoionization part of the equation to account for
the fact that the outflow can be optically thick near the planet. We continue
to use the case A recombination coefficient as our primary interest is the
ionization structure in the optically thin region far from the planet, which
sets the initial ionization fraction of the gas in the tail.

Equation 42 has to be solved iteratively as $\tau_{\bar{\nu}}(r^{\prime})$
depends on the ionization fraction of the outflow. However, we can approximate
$\tau_{\bar{\nu}}(r^{\prime})$ by assuming the atmosphere is neutral and
solving the differential equation directly (e.g., Dos Santos et al., 2022). We
solve this using an implementation of the LSODA method (Petzold, 1983)
provided in the scipy library (Virtanen et al., 2020). The absolute and
relative tolerances used are $1\times 10^{-13}$ and $1\times 10^{-13}$
respectively.

Although modelling the outflow inside the Hill sphere as a Parker wind is not
as accurate as a self-consistent EUV heated (e.g., Murray-Clay et al., 2009)
or bolometrically heated (e.g., Schulik & Booth, 2023) model, we emphasize
that the primary motivation of modelling the Hill sphere is to find the
initial conditions of the gas entering the tail. As mentioned previously, this
also allows us to remain agnostic to the exact physics that is driving the
outflow. This agnostic nature then allows us to test the prediction of various
classes of escape models (e.g., photoevaporation vs core-powered mass-loss).

###  2.3 Energetic Neutral Atoms

The interaction between the stellar wind and the planetary outflow also
results in charge exchange between the stellar wind protons and neutral
hydrogen in the outflow. There is little momentum transfer between the atoms
in this process (e.g., Pinto & Galli, 2008), and therefore it generates a
high-velocity neutral atom. The production of these high-velocity neutral
atoms, known as energetic neutral atoms (ENAs), has been proposed as the
source of the large Lyman-$\alpha$ transits in some systems (e.g., Holmström
et al., 2008; Kislyakova et al., 2014; Bourrier et al., 2016; Khodachenko et
al., 2019).

In this work, we adopt a fluid description of the production of ENAs based on
single-fluid simulations (e.g., Tremblin & Chiang, 2012; Debrecht et al.,
2022). This differs from the works cited above, which either use a particle-
based or multi-fluid approach to calculate the production of ENAs. The
differences in these approaches can lead to different conclusions about the
contribution of ENAs to the Lyman-$\alpha$ optical depth.

In the single fluid framework, the stellar wind protons cannot penetrate into
the planetary outflow and are forced around the flow. The resulting velocity
shear at the boundary between the stellar and planetary wind creates Kelvin-
Helmholtz instabilities. Charge exchange between stellar wind protons and
planetary neutral hydrogen atoms in this turbulent mixing layer produces ENAs.

We estimate the number density of ENAs in this mixing layer using analytic
arguments from Tremblin & Chiang (2012); Owen et al. (2023). The rate of
production of ENAs is given by:

$\displaystyle\beta(n_{+,*}n_{0,p}-n_{0,*}n_{+,p})$ (43)

where $n_{+,*}$ and $n_{+,p}$ are protons of stellar wind and planetary wind
origin respectively, $n_{0,*}$ and $n_{0,p}$ are neutral hydrogen atoms of
stellar and planetary origin respectively and $\beta$ is the rate coefficient.
As the collisional timescale is much shorter than the timescale for ENAs to
get advected off the limb of the star, we can assume that the stellar wind
protons and planetary atoms are in chemical equilibrium. Therefore, the number
density of ENAs in the mixing layer can be approximated as:

$\displaystyle n_{\text{ENA}}=Nn_{*}$ (44)

where $n_{*}$ is the density of the post-shocked stellar wind and $N$ is the
neutral fraction of hydrogen in the tail.

The contribution of ENAs to the Lyman-$\alpha$ optical depth also depends on
the velocity and temperature of these atoms and the size of the mixing layer.
The collision and subsequent diversion of the stellar wind around the tail may
cause the bulk velocity of ENAs in the mixing layer to be different to the
initial stellar wind velocity. As we do not have an a priori way of
calculating the bulk velocity of ENAs, we treat it as a constant that is
allowed to be a free parameter of the model.

Collisions between atoms of stellar and planetary origin in the mixing layer
cool the atoms of the stellar origin. As the charge exchange collisional
timescale is short compared to the flow scale, such that the neutrals are
constantly being exchanged with protons, this process is dominated by the
cooling of the stellar wind protons via collisions with protons in the
planetary wind. For simplicity, we will ignore this process and set the
temperature of ENAs to be the same as the stellar wind temperature.

Raga et al. (1995) analytically studied the size of the mixing layer for the
interaction between a spherical outflow with a planet parallel flow. To do
this calculation, they had to assume that the stand-off distance between the
two shocks was small compared to the radius of curvature of the shock front.
This is not the case for the interaction between the stellar wind and the
planetary outflow, where the shocked regions can be thick. However, using this
as a guide and inspecting the simulations performed by Tremblin & Chiang
(2012); Debrecht et al. (2022) allows us to estimate the length of the mixing
layer to range from a few to tens of per cent of the radius of the tail. As
this range is fairly large, we choose to treat the size of the mixing layer as
a free parameter, specified by $L_{\text{mix}}$, which represents the
fractional size of the mixing layer compared to the radius of curvature of the
planetary outflow that the stellar wind impinges on.

For ENAs generated due to the interaction of the stellar wind with the Hill
sphere of the planet, the mixing layer is given by the area between spheres
with Hill sphere radius $r_{H}$ and $r_{H}(1+L_{\text{mix}})$. For ENAs
generated in the tail, the mixing layer is given by the area between an
ellipse with a height and depth given by $H,D$ and an ellipse with a height
and depth given by $H(1+L_{\text{mix}}),D(1+L_{\text{mix}})$ respectively.

###  2.4 Constructing the Full Outflow Model

In order to construct the full outflow model, we need a way to transition from
the spherically symmetric outflow inside the Hill sphere into the cometary
tail-like outflow. We do this by "glueing" the tail onto a point on the Hill
sphere of the planet. The initial speed and ionization fraction of the gas in
the tail are specified by the values given by the spherical wind model
evaluated at the Hill sphere. With this formulation, there is still freedom to
choose exactly at what position to glue the tail onto the Hill sphere and the
initial direction of the outflow. To simplify this, we assume that the initial
position and velocity vector of the tail, in a frame centred on the planet,
are parallel. We can then encode this information in a free parameter $h$,
which we interpret as the relative specific angular momentum between the
outflow and the planet. We shall refer to this quantity as the angular
momentum deficit of the outflow, and it can be written as:

$\displaystyle
h=\Omega(a_{p}^{2}+2a_{p}R_{H}\text{sin}\theta+R_{H}^{2})+u\text{cos}\theta\sqrt{a_{p}^{2}+2a_{p}R_{H}\text{sin}\theta+R_{H}^{2}}$
(45)

where $\theta$ is the angle at which the tail is glued onto the Hill sphere of
the planet, and $a_{p}$ is the semimajor axis of the planet. We estimate the
projected area of this transition region blocking the star relative to the
rest of the tail is $\sim\frac{R_{H}^{2}}{R_{*}H}$, which is generally small,
and therefore this uncertainty should not have a large effect on our results.
Furthermore, the primary aim of this work is to model the transit on scales
much larger than the Hill sphere of the planet where the geometry of the model
matches 3D simulations.

The full outflow model requires 14 input parameters, shown in Table 1. These
can be roughly divided into "tightly constrained" and "weakly constrained"
parameters. The "tightly constrained" are ones, like the planet’s or the
star’s radius, that are generally well constrained by independent
observations. In general, when comparing to observations, we would allow the
"weakly constrained" parameters to be free and the "tightly constrained"
parameters to be fixed or use these independent observations as tight priors.

Table 1: The input parameters for our tail model. The fixed parameters are well constrained from other observations. The free parameters are the ones we would like to infer from observations, particularly the mass loss rate and sound speed of the escaping gas.  Tightly Constrained |  Weakly Constrained
---|---
Planet radius, $R_{p}$ |  Sound speed, $c_{s}$
Planet mass, $M_{p}$ |  Planet mass loss rate, $\dot{M}_{p}$
Stellar mass, $M_{*}$ |  Stellar wind velocity (at planet), $u_{*}$
Semimajoraxis of orbit, $a_{p}$ |  Star mass loss rate, $\dot{M}_{*}$
Orbital inclination, $i$ |  Optically thin photoionization rate (at planet), $\Gamma_{p}$
|  Stellar wind temperature, $T_{\text{sw}}$
|  Angular momentum deficit of outflow, $h$
|  Bulk velocity of ENAs, $u_{\text{ENA}}$
|  Length of turbulent mixing layer, $L_{\text{mix}}$

##  3 Comparison to 3D Hydrodynamics Simulations

To demonstrate that this model is accurate in predicting the geometry of the
outflowing gas, we compare the trajectory that our model predicts to those
obtained from 3D hydrodynamics simulations. We do this for four simulation
setups described in McCann et al. (2019); Khodachenko et al. (2019); MacLeod &
Oklopčić (2022); and Hazra et al. (2022). These simulations examine the
outflow for a range of planetary and stellar properties to highlight the
model’s applicability over a wide range of parameters. Specifically, the
masses of the planets we compare to range from $0.07-1.14\ M_{J}$, with mass
loss rates ranging from $0.2-10\times 10^{10}\text{ g s}^{-1}$ and stellar
type ranging from M to F.

For the comparison to be meaningful, we fix the free parameters in our model
to be the same as those used in the respective simulation. In each case, the
planetary mass loss rate, stellar wind mass loss rate, and stellar wind
temperature are provided in the simulation paper. They also provide a
temperature plot of the outflow, which we can use to estimate the sound speed
of the gas. These are estimated to be 11, 11, 9, and 6 km s-1 respectively. We
fix the value of the stellar wind velocity to match the value of stellar wind
velocity at the planet’s position in the simulations. Finally, the angular
momentum deficit, $h$, is not fixed from the initial simulation parameters, so
we approximate it by looking at the direction of the gas flow at the beginning
of the tail. Normalized to the specific angular momentum of the planet, the
angular momentum deficit is given by -0.05, 0.07, 0.02, and -0.06,
respectively.

Figure 2 shows the trajectory of the tail calculated using our model
superimposed on the outflow produced from 3D hydrodynamics simulations. Our
model is successful in reproducing the trajectory of the simulated tails in
all the simulations, and this gives us confidence that it can be used to test
atmospheric escape models quantitatively.

Figure 2: Each panel shows the geometry of the planetary outflow produced from
3D hydrodynamics simulations, covering a variety of planet and system
parameters. These are: a) Hot Jupiter based on HD209458b (McCann et al.,
2019), b) HD189733b, a Hot Jupiter (Hazra et al., 2022) , c) GJ436b, a Warm
Neptune (Khodachenko et al., 2019) and d) Sub-Jovian based on WASP-107
(MacLeod & Oklopčić, 2022). Using the same input parameters as were used in
the simulations, we calculate the trajectory of the tail using our model,
shown by the thick black line. The model shows good agreement with the
simulations in all cases.

##  4 Computing Synthetic Transit Profiles

The primary motivation of this model is to produce synthetic transmission
spectra that can be compared to observations. We use a ray-tracing scheme to
calculate the attenuation of Lyman-$\alpha$ due to the transit of the outflow
over the stellar disc. We produce observations at different orbital phases by
rotating the outflow around the system’s barycenter, which we approximate as
the star’s centre. To compute the attenuation of Lyman-$\alpha$ at a
particular orbital phase, we create a 3D cylindrical grid from the product of
a 2D grid on the stellar disc and a 1D grid along the line connecting the
centre of the stellar disc to the observer. We then define a new cartesian
coordinate system with the origin at the centre of the stellar disc. The
$x_{t}$ and $y_{t}$ axes are parallel to the stellar disc, and the $z_{t}$
axis is given by the line connecting the centre of the stellar disc to the
observer. The 2D grid consists of 705 cells. For ray tracing through the tail,
the 1D grid consists of a minimum of 10 cells across the width of the tail.
For ray tracing through the Hill sphere, the 1D grid consists of 15 cells. We
did resolution tests in order to check that this number of cells was
sufficient.

The optical depth through the outflow for a corresponding point on the stellar
disc is given by

$\displaystyle\tau_{\nu}(x_{t},y_{t})=\int_{0}^{\infty}n_{\text{HI}}(x_{t},y_{t},z_{t})\sigma_{\nu}(u_{z_{t}},T)dz_{t}$
(46)

where $n_{\text{HI}}$ is the number density of hydrogen atoms, $\sigma_{\nu}$
is the absorption cross-section of hydrogen to photons of frequency $\nu$,
which is a function of the line of sight velocity of the gas and the
temperature. This accounts for Lyman-$\alpha$ optical depth due to both
neutral hydrogen of planetary origin and ENAs of stellar origin.

The absorption cross-section is given by

$\displaystyle\sigma_{\nu}=\frac{\pi e^{2}}{m_{e}c}f\Phi(u_{z_{t}},T)$ (47)

where $e$ and $m_{e}$ is the charge and mass of the electron, $c$ is the speed
of light, $f=4.1641\times 10^{-1}$ is the oscillator strength of the
transition and $\Phi$ is the Voigt line profile (Kramida et al., 2022). The
Gaussian part of the Voigt profile has a standard deviation of $\sigma_{N}$,
the Lorentzian part has a half-width at half-maximum of $\gamma$, and the line
centre has been Doppler shifted to $\nu$. These are given below:

$\displaystyle\sigma_{N}=\sqrt{\frac{k_{b}T}{m_{H}c^{2}}}\nu_{0},\quad\gamma=\frac{A}{4\pi},\quad\nu=\nu_{0}\left(1-\frac{u_{z_{t}}}{c}\right)$
(48)

where $k_{b}$ is the Boltzmann constant. Here $\nu_{0}$ is the rest frequency
of Lyman-$\alpha$ and $A=6.2649\times 10^{8}\text{ s}^{-1}$ is the Einstein
$A$ coefficient of the transition (Kramida et al., 2022).

The obscuration fraction averaged over the stellar disc is:

$\displaystyle\mathcal{O}_{\nu}=\frac{\iint_{D}e^{-\tau_{\nu}(x_{t},y_{t})}dx_{t}dy_{t}}{\pi
R_{*}^{2}}$ (49)

Here, we do not account for spatial variations in the Lyman-$\alpha$ profile
across the stellar surface, as seen on the Sun (Gordino et al., 2022). To
numerically solve for the obscuration fraction, the integral in Equation 46 is
calculated by performing a middle Riemann sum using the 1D grid in the
$z_{t}$-direction. The integral in Equations 48 is calculated by performing a
double Riemann sum using the 2D grid over the stellar disc (e.g., Schreyer et
al., 2023).

##  5 Comparison to Real Systems

In order to assess the biases in our model, we check the ability to recover
known parameters from mock observations. For a parameter set
$\boldsymbol{\Theta}$, we calculate the obscuration due to the planetary
outflow as a function of time and wavelength. We bin the resulting spectrum
into 30-minute time bins and three velocity bins between [-150, -50] km s-1,
and add error bars $\sim 10\%$ to the resulting points. We make twenty
realisations of the synthetic data by re-sampling the points randomly within
its error. We perform a retrieval on each of these realizations. For each, we
are able to recover the true parameters within 1$\sigma$, and the combination
of the twenty realizations was not biased.

###  5.1 Fitting the transit of GJ 436 b

We use the outflow model to perform a retrieval of the warm Neptune GJ 436 b
in a Bayesian framework. We choose this planet because the Lyman-$\alpha$
transit of this planet has been extensively observed with the Hubble Space
Telescope (Kulow et al., 2014; Ehrenreich et al., 2015; Lavie et al., 2017;
dos Santos et al., 2019). The fixed and free parameters of the model are
listed in Table 2. Here, we list the values used for the fixed parameters and
the prior bounds for the free parameters. In order to reduce the number of
free parameters, we set the temperature of the stellar wind to 0.5 MK. We
checked that this choice does not affect our conclusions by running models
with $T_{\text{sw}}=1$ MK, in which we found similar results. In the fitting
procedure, we convert the angular momentum deficit, $h$, into the "launch
angle" of the outflow, given by $\theta$ in 45. We do this so that we can put
a uniform prior on the launch angle. For all other free parameters, apart from
the inclination, we use log uniform priors, as the scale of these parameters
is not constrained. The prior bounds enclose realistic parameter bounds as
determined from theory or other observations. We use the energy-limited mass
loss rate,

$\displaystyle\dot{M}_{\text{EL}}=\eta\frac{\pi
F_{\text{XUV}}R_{P}^{3}}{GM_{p}}$ (50)

with the efficiency $\eta=1$ to set an upper bound on the planetary mass loss
rate. Here, the EUV flux is approximated from the optically thin
photoionization rate by assuming all ionizing photons deliver $20\text{ eV}$
(see Murray-Clay et al. (2009), Appendix B). We note that the appropriate
radius over which the planet absorbs EUV flux may be larger than the optical
radius, which can enhance this mass loss rate limit. However, this is often
balanced out by the real efficiency for escaping planets being significantly
less than one. Therefore, in our case, it is reasonable to simply use the
optical radius.

Table 2: The Fixed Parameters  Parameter |  Value/Prior
---|---
Fixed |
Planet Mass, $M_{p}$ ($M_{J}$) |  0.07
Planet Radius, $R_{p}$ ($R_{J}$) |  0.35
Semimajoraxis, $a_{p}$ (AU) |  0.029
Stellar Mass, $M_{*}$ ($M_{\odot}$) |  0.45
Stellar Radius, $R_{*}$ ($R_{\odot})$ |  0.425
Stellar Wind Temperature, $T_{\text{sw}}$ (MK) |  0.5
Free |
Log Sound speed (cm s-1) |  U(5.2, 6.5)
Log Planet mass loss rate (g s-1) |  U(8, $\dot{M}_{\text{EL}}$)
Log Stellar wind velocity (at planet) (cm s-1) |  U(6.5, 8)
Log Star mass loss rate (g s-1) |  U(10.3, 13)
Log Optically thin photoionization rate (at planet) (s-1) |  U(-5.6, -2.6)
Launch Angle (radians) |  U($\frac{\pi}{2}$, $\pi$)
Inclination (radians) |  N(1.51, 0.02)
Log Bulk velocity of ENAs (cm s-1) |  U(6.4, $u_{*}$)
Log Length of turbulent mixing layer |  U(-2, -0.4)

We sample the posterior using the Markov chain Monte Carlo (MCMC) ensemble
sampler emcee (Foreman-Mackey et al., 2013). This is a Python implementation
of Goodman & Weare’s affine invariant MCMC ensemble sampler (Goodman & Weare,
2010). We calculate the log-likelihood by comparing the observed flux
$\mathcal{F}_{o}$ to the modelled flux $\mathcal{F}$ as functions of both time
and wavelength within a collection of observed time series spectra as:

$\displaystyle\text{ln}\mathcal{L}(\boldsymbol{\Theta})=-\frac{1}{2}\sum_{n}\left[\frac{(\mathcal{F}_{o}-\mathcal{F})^{2}}{\sigma^{2}_{n}}+\text{ln}(2\pi\sigma^{2}_{n})\right]$
(51)

The time-series spectra consist of all Lyman-$\alpha$ observations of GJ 436
to date taken with the Hubble Space Telescope’s (HST’s) Space Telescope
Imaging Spectrograph (STIS) configured with the G140M grating. This consists
of 26 individual observations obtained between 2010 and 2016 (program PIs
Ehrenreich, Green, and France). We uniformly re-extract the spectra using the
calSTIS pipeline, manually identifying the position of the spectral trace on
the detector. We use an extraction ribbon of 11 pixels and background ribbons
of 20 pixels at $\pm$20 pixels from the trace location to minimize effects
from a spatially-variable background systematic known as FUV glow. The
observations are divided into subexposures of roughly 500 s each, with minor
adjustments to ensure the full use of each exposure.

For comparison to the data, we forward-model the observed Lyman-$\alpha$
spectrum during each subexposure. We do this by applying the absorption
predicted by the model ray tracing to a de-broadened version of the stellar
Lyman-$\alpha$ profile as observed at earth reconstructed using the technique
described in Wilson et al. (2022). The profile is then broadened by the
instrument line spread function and binned to match the wavelength grid of the
data. The observed and modelled spectra are then normalized to an identical
flux in the range of 90 and 250 km s-1. This accounts for known variability in
instrument throughput and stellar variability in which the full spectrum is
scaled by a constant. We note that this method cannot take into account
variability in which the ratio of flux in the red wing to the blue wing
changes. After normalization, we use the predicted spectrum and flux scale of
the instrument to estimate the expected number of counts in each wavelength
bin, then use this to compute the expected variance due to Poisson statistics.
This mitigates problems due to zero-flux, zero-error points in the calSTIS-
reduced data and the disproportional influence of low-flux, low-error points
on fits and is more appropriate for our Bayesian framework. We use the
resulting normalized spectra and estimated errors to compute the data
likelihood with Equation 51.

We ran the MCMC for a minimum of 15000 steps with 100 walkers. In order to
monitor convergence, we ran two independent chains for each model so that we
could calculate the split chains potential scale reduction factor,
split-$\hat{R}$ (Gelman et al., 2014). The maximum value across all was 1.02
across all parameters. We also calculated the effective sample size for each
of the parameters, which was a minimum of 15. Based on these results, we
conclude that the posterior we have calculated from the MCMC is representative
of the true posterior.

##  6 Results and Discussion

To explore how including different parts of the outflow affects inferences, we
run four models of increasing complexity. In the first, we only include the
tail outflow; in the second, we include both the tail and the hill sphere
outflow; in the third, we include the tail outflow and ENAs; and in the
fourth, we include all model components (i.e. the tail, hill sphere outflows,
and ENAs). In fitting the data, we must only compare the parts of the transit
in which the model is applicable. In the models where we do not include the
outflow inside the Hill sphere, we only compare to data after the egress of
the Hill sphere ($>1.3\text{ hours}$). In the cases where we include the
outflow inside the Hill sphere, we compare to data after the egress of the
planet ($>0.7\text{ hours}$). This is in order to avoid comparing to data
where parts of the outflow in front of the planet may be blocking the star,
where our model is not applicable. The retrieved parameters from all four
models are consistent, indicating the tail part of the outflow drives our
inferences. This is expected for GJ 436 b as the tail dominates the transit
duration. One reason for this is that we find that ENAs contribute very little
to the Lyman-$\alpha$ optical depth. We verify this by calculating the
spectrum for our best-fit model without the contribution of ENAs and find it
virtually indistinguishable from the spectrum including ENAs. We do find a
slight inconsistency in retrieved inclinations between the models with and
without the inclusion of the Hill sphere. In the cases with the Hill sphere,
the inclination has a median of 1.55 radians, whilst in the cases without the
Hill sphere, it has a median of 1.52 radians, more in line with observational
constraints. We think this discrepancy comes from the fact that cutting the
gas off at the Hill sphere underpredicts the extent of the gas in this region,
and therefore, the inclination has to change to compensate. However, the fact
the rest of the retrieved parameters are the same gives us confidence that
most of the information is coming from the tail. As a result, we focus on the
full model from now on.

Figure 3: Left: Posterior distribution of the sound speed, planet mass loss
rate, and photoionization rate for GJ 436 b. On the joint sound speed and
photoionization rate posterior, we have plotted lines of constant ratio
corresponding to $\frac{c_{s}}{\Gamma_{p}}=10^{9.7-10.1}$ cm. On the joint
mass loss rate and photoionization rate posterior, we have plotted lines of
constant efficiency from the energy-limited mass loss formula (Equation 50).
Right: Posterior distribution of the planet mass loss rate, stellar wind ram
pressure, stellar mass loss rate. On the joint planet mass loss rate and
stellar wind ram pressure posterior, we plot lines of constant ratio
corresponding to $\frac{\dot{M}_{*}u_{*}}{\dot{M}_{p}}=10^{9.1-9.7}\text{ cm
s}^{-1}$. Although planetary mass loss rate and stellar wind ram pressure
cannot be tightly constrained individually, the tight correlation between
these parameters means that this ratio is tightly constrained.

In Figure 3, we present the posterior distributions for the sound speed,
planetary mass loss rate, optically thin photoionization rate, mass loss rate
of the star and ram pressure of the stellar wind. We provide the full
posterior distribution in Appendix A. In Figure 4, we plot the model
transmission spectra at different times in the transit for the maximum a
posteriori parameter estimates (i.e. posterior mode). It is clear that the
model is a good match for the spectrum, providing evidence that the model
accurately represents the planetary outflow. This is further illustrated by
Figure 5, where we show the transit integrated between -50 and -150 km s-1.

Figure 4: Lyman-$\alpha$ transit spectrum of GJ 436 b at different times
(relative to optical mid-transit) measured using Hubble Space Telescope (black
points) and the model (solid blue) using the MAP estimates from the MCMC
sampler. The dashed green line is the out-of-transit model that is used to
calculate the transmission spectrum. Figure 5: Lyman-$\alpha$ transit of GJ
436 b integrated over the blue wing interval of [-150, 50] km s-1. The shaded
region is before the egress of the planet, where the model is not applicable.

###  6.1 Mass Loss Rate

We find a total hydrogen mass loss rate of $2.6^{+3.5}_{-1.6}\times
10^{9}\text{ g s}^{-1}$. This is similar to the mass loss rate of $2\times
10^{9}$ g s-1 found in the 1D simulations by Schulik & Booth (2023), and
$5.5\times 10^{9}\text{ g s}^{-1}$ found in 3D hydrodynamics simulations by
Villarreal D’Angelo et al. (2021), which were able to reproduce the depth and
duration of the transit. Considering the joint mass loss rate and
photoionization, we can also put constraints on the efficiency of the outflow,
$\eta$, given in the energy-limited formula (Equation 50). On the posterior,
we plot lines of constant mass loss efficiency and are able to constrain this
to range from 0.03-0.5. This is consistent with simulations of atmospheric
escape from GJ 436 b that accounted for both EUV and bolometric heating and
found an efficiency of $\sim 0.05$ (Schulik & Booth, 2023). More generally, it
is consistent with the theoretical estimates for a photoevaporative flow (e.g.
Owen & Jackson, 2012; Owen & Alvarez, 2016). We highlight that as a real
stellar spectrum has significant flux upwards of 20 eV, these efficiencies are
likely underestimated by a factor of a few (see Appendix B).

The slope of the correlation between mass loss rate and flux is much weaker
than the energy-limited slope. Even though higher EUV fluxes lead to a higher
mass loss rate, the outflow is more ionized to the extent that the density of
neutral hydrogen actually decreases. Therefore, in order to match the observed
optical depth to Lyman-$\alpha$ photons from neutral hydrogen, the model has
to compensate for a higher EUV flux with more mass loss relative to the
energy-limited relationship.

Interestingly, we put a tight constraint on the ratio of the ram pressure of
the stellar wind at the planet and the mass loss rate of the planet of
$\frac{\dot{M}_{*}u_{*}}{\dot{M}_{p}}=10^{9.4\pm 0.3}\text{ cm s}^{-1}$. The
reason for this is that this ratio determines the radial acceleration of the
tail and, hence the radial velocity of the tail as a function of time. We can
see the effect acceleration of the tail by considering the ratio of observed
flux in the (-150,-100) km s-1 band and (-100,-50) km s-1 band, shown in
Figure 6. As the tail is accelerated, more flux is absorbed at higher
velocities, so the ratio decreases. We note that this ratio also changes as
the column density of neutral hydrogen changes, and therefore, to extract the
effect of the acceleration of the tail, we have shown how this ratio evolves
if the gas in the tail was fully neutral.

We highlight that acceleration of the tail (Equation 8) is not exactly
proportional to the ram pressure of the stellar wind but depends on the
difference between the radial velocity of the tail and the stellar wind. The
maximum velocity that the tail can be accelerated to is the stellar wind
velocity. If the tail reaches this terminal velocity whilst still attenuating
Lyman-$\alpha$ photons, then very strong constraints can be placed on the
stellar wind velocity. In the case of GJ 436 b, the tail is still
accelerating, and therefore, we are only able to put a lower limit of 300 km
s-1 on the stellar wind velocity at the planet.

We emphasize that we do not dismiss the case where the stellar wind velocity
is very high (> 800 km s-1), and the mass loss rate is low as GJ 436 b has a
polar orbit (Bourrier et al., 2022). For the Sun, the polar wind is
significantly faster and less dense than the equatorial wind, and therefore,
it is possible that this is the case for GJ 436 as well.

####  6.1.1 Sound Speed of the Outflow

The $2\sigma$ lower bound for the retrieved sound speed is $13\text{ km
s}^{-1}$, corresponding to a temperature of $\sim$ 10, 000 K for a fully
ionized purely hydrogen gas, which is consistent with a photoevaporative
outflow. The bulk of the posterior encompasses higher sound speeds, with a MAP
estimate of $20.5\text{ km s}^{-1}$, which corresponds to outflow temperature
$\sim$ 40,000 K. This is higher than expected for photoevaporative outflows
which are usually thermostat at temperatures $\sim$ 10,000 K, due to strong
radiative cooling of hydrogen above those temperatures.

It is possible that another physical mechanism is accelerating the outflow to
high velocities by the time it reaches the Hill sphere, raising the apparent
sound speed. This interpretation is compatible with the results of Bourrier et
al. (2016), which used a collisionless outflow model of the exosphere, and
required outflow velocities of $\sim 50-60\text{ km s}^{-1}$ in order to fit
the data. A possible mechanism to produce these high velocities is via the
dissipation of MHD waves in the upper atmosphere (e.g Tanaka et al., 2014);
however, consistent multi-dimensional models are needed to prove that this is
a viable mechanism to drive high-velocity outflows.

It is also feasible that charge exchange is responsible for more
Lyman-$\alpha$ optical depth than predicted by our model at early times, and
the sound speed of the outflow needed to be raised in order to compensate for
this. Our single-fluid model for charge exchange predicts a significantly
lower contribution to ENAs compared to 3D multi-fluid hydrodynamic simulations
performed by Khodachenko et al. (2019), which are able to reproduce the depths
of the large transit at early times with 3D multi-fluid hydrodynamic
simulations. This is predominantly because they find a large number of ENAs
generated throughout a thick shocked layer between the stellar wind bow shock
and the planetary outflow. More simulations are needed to figure out which of
these models are correct.

One potential way to break the degeneracy between the contribution of the
planetary outflow and ENAs is to consider the distinction between the
velocity-resolved transit from each of these factors. In the case where ENAs
generated throughout a thick layer are responsible for most of the
Lyman-$\alpha$ optical depth, the velocity of peak absorption should stay
fairly constant. This is because the stellar wind protons must be able to
penetrate far into the planetary outflow and, therefore, should not have lost
their original velocity; hence, the velocity of ENAs in the stellar frame
should be fairly constant as a function of transit time. In contrast, if the
planetary wind is responsible for the majority of Lyman-$\alpha$ optical
depth, then the velocity of peak absorption should increase over the duration
of the transit as the tail is accelerated away from the star. In this case, we
cannot definitively tell that the absorbing material is being accelerated as a
function of time (model independently) because the gas becomes ionized too
quickly for the velocity of peak absorption to change dramatically (see Figure
6). However, for planets with lower ionization rates and therefore longer
tails or when the gas in the tail is in ionization-recombination equilibrium,
it may be possible to "see" the acceleration of the gas in the lightcurve.

Figure 6: The observed flux ratio between (-150, -100) and (-100, -50) km s-1
as a function of time after the transit (black dots). The blue line is the
best-fit model, the orange line is the best-fit model, setting the
photoionization rate to zero, and the dashed green line is the out-of-transit
ratio. There is a smooth decrease in this ratio from $\sim$ 1 hour to 5 hours
as the tail is accelerated to higher velocities and becomes more ionized,
absorbing more photons in the (-150, -100) km s-1 band compared to the (-100,
-50) km s-1 band. The model with ionization rate set to zero indicates the
length of the transit required in order to unambiguously detect the
acceleration of the tail (as ionization also changes the flux ratio).

####  6.1.2 Limitations and Future Improvements

Since we only model the gas that is funnelled behind the planet into a
cometary tail, we do not attempt to match pre-transit absorption. As a result,
we are not able to account for the information about the outflow in the pre-
transit data. The challenge with modelling the outflow leading the planet is
that it is unstable, resulting in the variable occultation of the star (e.g.,
McCann et al., 2019) as seen in the case of AU Mic b (Rockcliffe et al.,
2023). The reason for this instability is that the Coriolis force and the
stellar wind pressure both act in the same direction; therefore, the outflow
cannot be in generalised geostrophic balance and reach a stable density
structure. This problem is only pertinent when the planetary wind and the
stellar wind have similar strengths. However, Lyman-$\alpha$ observations,
which show an extended outflow leading the planet (e.g., Rockcliffe et al.,
2023) make understanding this case important. We leave modelling of the
outflow leading the planet to future work. In the case of GJ 436 b, where
there is minimal pre-transit absorption, the stellar wind is most likely
strong enough to stop a leading outflow. In this case, the time of ingress
provides an approximate location for the boundary between the stellar and
planetary wind, which can be used to infer their relative pressures at this
location.

Our model does not include the effect of magnetic fields, as we use purely
hydrodynamic simulations as the basis for our work. Both the planetary
magnetic field and the stellar wind magnetic field can affect the geometry of
the outflow (e.g., Carolan et al., 2021b). However, it is not yet clear
whether close-in exoplanets host strong magnetic fields (e.g., Schreyer et
al., 2023) or whether the stellar wind magnetic field plays an important role
in the interaction with the planetary outflow as the ram pressure of the
stellar wind dominates over the magnetic pressure outside the Alfvén surface.

Similarly, we do not include the effects of radiation pressure. Like the
stellar wind, this can impart momentum to the planetary outflow, accelerating
it radially. For GJ 436 b, radiation pressure is unlikely to be important
(e.g., Khodachenko et al., 2019). However, this may not be the case for all
systems. In the future, we will extend the model to include a term accounting
for this force. Another possible improvement to this model would be to include
a more realistic stellar wind accounting for the orbital velocity of the
planet and variations along the orbit.

We also imposed a circular orbit on the planet. Since highly irradiated
planets are close to their host star, they are generally assumed to have
undergone orbital circularization. However, deviations from a circular orbit
can cause a line of sight velocity shift of $v_{k}e\text{cos}\omega$, where
$e$ is the eccentricity of the planet and $\omega$ is the argument of
periastron. In the case of GJ 436 b, the eccentricity and argument of
periastron are 0.156 and 325.8 degrees, respectively (Bourrier et al., 2022),
which equates to a line of sight velocity shift of $\sim-15\text{ km s}^{-1}$
at mid-transit. This is small compared to the velocities probed by
Lyman-$\alpha$ transits and, therefore cannot be fully responsible for the
high radial velocities observed. However, it may partially be responsible for
the higher-than-expected sound speed retrieved. In the future, it may be
worthwhile folding constraints on the planetary orbital dynamics into the
model.

##  7 Conclusions

We have created a model of the interaction of a planetary outflow with the
circumstellar environment, which can be used to perform synthetic
Lyman-$\alpha$ transits. We demonstrate that this model is able to reproduce
the geometry of outflows produced by 3D hydrodynamics simulations.
Furthermore, we illustrate the usage of this model to perform the first
retrieval of atmospheric escape properties using velocity-resolved
Lyman-$\alpha$ observations.

For GJ 436 b, we are able to constrain the sound speed of the escaping gas to
be $\gtrsim 10\text{ km s}^{-1}$, indicating that we can rule out core-powered
mass loss as the mechanism of escape for this planet. This is consistent with
the predictions from the combined photoevaporation and core-powered mass-loss
model of Owen & Schlichting (2023). Although the retrieved planetary mass loss
rate has wide bounds, we find that the ratio of planetary mass loss rate to
stellar wind ram pressure is tightly constrained. Independent constraints on
the space environment around this planet (e.g Bellotti et al., 2023) would
mean that we are able to put a stronger constraint on the mass loss rate.
Equally, as we expect this ratio to be tightly constrained if the tail is
observable, we can also use Lyman-$\alpha$ transits as a way to measure the
stellar winds of stars.

The usefulness of the simple model is that we are able to efficiently explore
the large parameter space to find solutions that match the data well. In doing
so, we are able to put constraints on parameters and find correlations between
parameters, which can both inform whether current models can fit the
observations and guide future 3D simulations. We do caution that making a
simplified model requires sacrificing some accuracy as compared to 3D
simulations, and therefore, further simulations are needed to confirm the
conclusions drawn from this study. At the same time, although much progress
has been made with large-scale simulations, we also note the physics included,
and hence, the validity to describe Lyman-$\alpha$ transits is still debated.
Therefore, we believe simple models that can provide reasonably accurate tests
of these are essential to progress in this field.

## Acknowledgements

JEO is also supported by a Royal Society University Research Fellowship. ES &
JEO. has also received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 853022, PEVAP). The authors acknowledge support from NASA XRP
grant 80NSSC23K0282. This research is based on observations made with the
NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are
associated with program(s) 11817, 12034, 12965, 13650, and 14222. This work
benefited from the 2022 Exoplanet Summer Program in the Other Worlds
Laboratory (OWL) at the University of California, Santa Cruz, a program funded
by the Heising-Simons Foundation. For the purpose of open access, the authors
have applied a Creative Commons Attribution (CC-BY) license to any Author
Accepted Manuscript version arising.

## Data Availability

The data used in this work are publicly available through the Barbara A.
Mikulski Archive for Space Telescopes and have been aggregated under DOI
10.17909/pjvt-h370111http://dx.doi.org/10.17909/pjvt-h370. The code used to
generate these results will be made public upon acceptance of this manuscript.

## References

  * Allart et al. (2018) Allart R., et al., 2018, Science, 362, 1384
  * Bellotti et al. (2023) Bellotti S., et al., 2023, arXiv e-prints, p. arXiv:2306.15391
  * Ben-Jaffel et al. (2022) Ben-Jaffel L., et al., 2022, Nature Astronomy, 6, 141
  * Bisikalo et al. (2013) Bisikalo D., Kaygorodov P., Ionov D., Shematovich V., Lammer H., Fossati L., 2013, ApJ, 764, 19
  * Bourrier & Lecavelier des Etangs (2013) Bourrier V., Lecavelier des Etangs A., 2013, A&A, 557, A124
  * Bourrier et al. (2013) Bourrier V., et al., 2013, A&A, 551, A63
  * Bourrier et al. (2015) Bourrier V., Ehrenreich D., Lecavelier des Etangs A., 2015, A&A, 582, A65
  * Bourrier et al. (2016) Bourrier V., Lecavelier des Etangs A., Ehrenreich D., Tanaka Y. A., Vidotto A. A., 2016, A&A, 591, A121
  * Bourrier et al. (2017) Bourrier V., Ehrenreich D., King G., Lecavelier des Etangs A., Wheatley P. J., Vidal-Madjar A., Pepe F., Udry S., 2017, A&A, 597, A26
  * Bourrier et al. (2022) Bourrier V., et al., 2022, A&A, 663, A160
  * Byrne & Hindmarsh (1975) Byrne G. D., Hindmarsh A. C., 1975, ACM Trans. Math. Softw., 1, 71
  * Caldiroli et al. (2021) Caldiroli A., Haardt F., Gallo E., Spinelli R., Malsky I., Rauscher E., 2021, A&A, 655, A30
  * Carolan et al. (2021a) Carolan S., Vidotto A. A., Villarreal D’Angelo C., Hazra G., 2021a, MNRAS, 500, 3382
  * Carolan et al. (2021b) Carolan S., Vidotto A. A., Hazra G., Villarreal D’Angelo C., Kubyshkina D., 2021b, MNRAS, 508, 6001
  * Carroll-Nellenback et al. (2016) Carroll-Nellenback J., Frank A., Liu B., Quillen A. C., Blackman E. G., Dobbs-Dixon I., 2016, Monthly Notices of the Royal Astronomical Society, 466, 2458
  * Chen & Rogers (2016) Chen H., Rogers L. A., 2016, ApJ, 831, 180
  * Clarke & Alexander (2016) Clarke C. J., Alexander R. D., 2016, Monthly Notices of the Royal Astronomical Society, 460, 3044
  * Cranmer (2004) Cranmer S. R., 2004, American Journal of Physics, 72, 1397
  * Debrecht et al. (2020) Debrecht A., Carroll-Nellenback J., Frank A., Blackman E. G., Fossati L., McCann J., Murray-Clay R., 2020, MNRAS, 493, 1292
  * Debrecht et al. (2022) Debrecht A., Carroll-Nellenback J., Frank A., Blackman E. G., Fossati L., Murray-Clay R., McCann J., 2022, MNRAS,
  * Dos Santos et al. (2022) Dos Santos L. A., et al., 2022, A&A, 659, A62
  * Ehrenreich et al. (2015) Ehrenreich D., et al., 2015, Nature, 522, 459
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fukue & Okada (1990) Fukue J., Okada R., 1990, PASJ, 42, 249
  * Fulton et al. (2017) Fulton B. J., et al., 2017, The Astronomical Journal, 154, 109
  * García Muñoz et al. (2020) García Muñoz A., Youngblood A., Fossati L., Gandolfi D., Cabrera J., Rauer H., 2020, ApJ, 888, L21
  * Gelman et al. (2014) Gelman A., Carlin J. B., Stern H. S., Dunson D. B., Vehtari A., Rubin D. B., 2014, Bayesian Data Analysis
  * Ginzburg et al. (2018) Ginzburg S., Schlichting H. E., Sari R., 2018, Monthly Notices of the Royal Astronomical Society, 476, 759
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Gordino et al. (2022) Gordino M., et al., 2022, A&A, 657, A86
  * Gupta & Schlichting (2019) Gupta A., Schlichting H. E., 2019, Monthly Notices of the Royal Astronomical Society, 487, 24
  * Gupta & Schlichting (2020) Gupta A., Schlichting H. E., 2020, MNRAS, 493, 792
  * Gupta et al. (2022) Gupta A., Nicholson L., Schlichting H. E., 2022, MNRAS, 516, 4585
  * Hazra et al. (2022) Hazra G., Vidotto A. A., Carolan S., Villarreal D’Angelo C., Manchester W., 2022, MNRAS, 509, 5858
  * Ho & Van Eylen (2023) Ho C. S. K., Van Eylen V., 2023, MNRAS, 519, 4056
  * Holmström et al. (2008) Holmström M., Ekenbäck A., Selsis F., Penz T., Lammer H., Wurz P., 2008, Nature, 451, 970
  * Hu et al. (2015) Hu R., Seager S., Yung Y. L., 2015, ApJ, 807, 8
  * Izidoro et al. (2022) Izidoro A., Schlichting H. E., Isella A., Dasgupta R., Zimmermann C., Bitsch B., 2022, ApJ, 939, L19
  * Jin & Mordasini (2018) Jin S., Mordasini C., 2018, The Astrophysical Journal, 853, 163
  * Kasting & Pollack (1983) Kasting J. F., Pollack J. B., 1983, Icarus, 53, 479
  * Kasting et al. (1984) Kasting J. F., Pollack J. B., Ackerman T. P., 1984, Icarus, 57, 335
  * Khodachenko et al. (2019) Khodachenko M. L., Shaikhislamov I. F., Lammer H., Berezutsky A. G., Miroshnichenko I. B., Rumenskikh M. S., Kislyakova K. G., Dwivedi N. K., 2019, ApJ, 885, 67
  * Kirk et al. (2020) Kirk J., Alam M. K., López-Morales M., Zeng L., 2020, AJ, 159, 115
  * Kislyakova et al. (2014) Kislyakova K. G., Holmström M., Lammer H., Odert P., Khodachenko M. L., 2014, Science, 346, 981
  * Kramida et al. (2022) Kramida A., Yu. Ralchenko Reader J., and NIST ASD Team 2022, NIST Atomic Spectra Database (ver. 5.10), [Online]. Available: https://physics.nist.gov/asd [2017, April 9]. National Institute of Standards and Technology, Gaithersburg, MD.
  * Kubyshkina et al. (2018) Kubyshkina D., et al., 2018, A&A, 619, A151
  * Kulow et al. (2014) Kulow J. R., France K., Linsky J., Loyd R. O. P., 2014, ApJ, 786, 132
  * Landsman & Simon (1993) Landsman W., Simon T., 1993, ApJ, 408, 305
  * Lavie et al. (2017) Lavie B., et al., 2017, A&A, 605, L7
  * Lecavelier Des Etangs et al. (2010) Lecavelier Des Etangs A., et al., 2010, A&A, 514, A72
  * Lee & Connors (2021) Lee E. J., Connors N. J., 2021, ApJ, 908, 32
  * Lee et al. (2022) Lee E. J., Karalis A., Thorngren D. P., 2022, arXiv e-prints, p. arXiv:2201.09898
  * Lopez & Fortney (2013) Lopez E. D., Fortney J. J., 2013, ApJ, 776, 2
  * Loyd et al. (2020) Loyd R. O. P., Shkolnik E. L., Schneider A. C., Richey-Yowell T., Barman T. S., Peacock S., Pagano I., 2020, ApJ, 890, 23
  * Luger & Barnes (2015) Luger R., Barnes R., 2015, Astrobiology, 15, 119
  * Lundkvist et al. (2016) Lundkvist M. S., et al., 2016, Nature Communications, 7, 11201
  * Luque & Pallé (2022) Luque R., Pallé E., 2022, Science, 377, 1211
  * MacLeod & Oklopčić (2022) MacLeod M., Oklopčić A., 2022, ApJ, 926, 226
  * Malsky & Rogers (2020) Malsky I., Rogers L. A., 2020, ApJ, 896, 48
  * Malsky et al. (2023) Malsky I., Rogers L., Kempton E. M. R., Marounina N., 2023, Nature Astronomy, 7, 57
  * Matsakos et al. (2015) Matsakos T., Uribe A., Königl A., 2015, A&A, 578, A6
  * McCann et al. (2019) McCann J., Murray-Clay R. A., Kratter K., Krumholz M. R., 2019, ApJ, 873, 89
  * Muñoz (2007) Muñoz A., 2007, Planetary and Space Science, 55, 1426
  * Murray-Clay et al. (2009) Murray-Clay R. A., Chiang E. I., Murray N., 2009, ApJ, 693, 23
  * Owen & Alvarez (2016) Owen J. E., Alvarez M. A., 2016, ApJ, 816, 34
  * Owen & Jackson (2012) Owen J. E., Jackson A. P., 2012, MNRAS, 425, 2931
  * Owen & Schlichting (2023) Owen J. E., Schlichting H. E., 2023, arXiv e-prints, p. arXiv:2308.00020
  * Owen & Wu (2013) Owen J. E., Wu Y., 2013, ApJ, 775, 105
  * Owen & Wu (2017) Owen J. E., Wu Y., 2017, The Astrophysical Journal, 847, 29
  * Owen et al. (2023) Owen J. E., et al., 2023, MNRAS, 518, 4357
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Petzold (1983) Petzold L., 1983, SIAM Journal on Scientific and Statistical Computing, 4, 136
  * Pinto & Galli (2008) Pinto C., Galli D., 2008, A&A, 484, 17
  * Raga et al. (1995) Raga A. C., Cabrit S., Canto J., 1995, MNRAS, 273, 422
  * Rockcliffe et al. (2021) Rockcliffe K. E., et al., 2021, The Astronomical Journal, 162, 116
  * Rockcliffe et al. (2023) Rockcliffe K. E., Newton E. R., Youngblood A., Duvvuri G. M., Plavchan P., Gao P., Mann A. W., Lowrance P. J., 2023, AJ, 166, 77
  * Rogers & Owen (2021) Rogers J. G., Owen J. E., 2021, MNRAS, 503, 1526
  * Rogers et al. (2021) Rogers J. G., Gupta A., Owen J. E., Schlichting H. E., 2021, MNRAS, 508, 5886
  * Rogers et al. (2023) Rogers J. G., Janó Muñoz C., Owen J. E., Makinen T. L., 2023, MNRAS, 519, 6028
  * Schreyer et al. (2023) Schreyer E., Owen J. E., Spake J. J., Bahroloom Z., Di Giampasquale S., 2023, arXiv e-prints, p. arXiv:2302.10947
  * Schulik & Booth (2023) Schulik M., Booth R. A., 2023, MNRAS, 523, 286
  * Schulreich & Breitschwerdt (2011) Schulreich M. M., Breitschwerdt D., 2011, A&A, 531, A13
  * Sing et al. (2019) Sing D. K., et al., 2019, AJ, 158, 91
  * Spake et al. (2018) Spake J. J., et al., 2018, Nature, 557, 68
  * Spitzer (1978) Spitzer L., 1978, Physical processes in the interstellar medium, doi:10.1002/9783527617722.
  * Szabó & Kiss (2011) Szabó G. M., Kiss L. L., 2011, The Astrophysical Journal, 727, L44
  * Tanaka et al. (2014) Tanaka Y. A., Suzuki T. K., Inutsuka S.-i., 2014, ApJ, 792, 18
  * Tremblin & Chiang (2012) Tremblin P., Chiang E., 2012, Monthly Notices of the Royal Astronomical Society, 428, 2565
  * Valencia et al. (2010) Valencia D., Ikoma M., Guillot T., Nettelmann N., 2010, A&A, 516, A20
  * Van Eylen et al. (2018) Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., Fulton B. J., Petigura E., Snellen I., 2018, Monthly Notices of the Royal Astronomical Society, 479, 4786
  * Venturini et al. (2020) Venturini J., Guilera O. M., Haldemann J., Ronco M. P., Mordasini C., 2020, A&A, 643, L1
  * Vidal-Madjar et al. (2003) Vidal-Madjar A., Lecavelier des Etangs A., Désert J. M., Ballester G. E., Ferlet R., Hébrard G., Mayor M., 2003, Nature, 422, 143
  * Vidal-Madjar et al. (2004) Vidal-Madjar A., et al., 2004, ApJ, 604, L69
  * Villarreal D’Angelo et al. (2021) Villarreal D’Angelo C., Vidotto A. A., Esquivel A., Hazra G., Youngblood A., 2021, MNRAS, 501, 4383
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Waugh (1961) Waugh 1961, Technical report, Evaluation Of Integral Of Elliptic Gaussian Distribution Over A Centred Ellipse, https://apps.dtic.mil/sti/pdfs/AD0271430.pdf. U.S Navy, https://apps.dtic.mil/sti/pdfs/AD0271430.pdf
  * Wilson et al. (2022) Wilson D. J., et al., 2022, ApJ, 936, 189
  * Wu (2019) Wu Y., 2019, The Astrophysical Journal, 874, 91
  * Yelle (2004) Yelle R. V., 2004, Icarus, 170, 167
  * Zeng et al. (2019) Zeng L., et al., 2019, Proceedings of the National Academy of Science, 116, 9723
  * Zhang et al. (2023) Zhang M., Knutson H. A., Dai F., Wang L., Ricker G. R., Schwarz R. P., Mann C., Collins K., 2023, AJ, 165, 62
  * dos Santos et al. (2019) dos Santos L. A., et al., 2019, A&A, 629, A47
  * dos Santos et al. (2020) dos Santos L. A., et al., 2020, A&A, 634, L4

##  Appendix A Posterior Plots

We provide the full posterior plots for the model outlined in Section 6 in
Figure 7.

Figure 7: The full posterior plot for our model shown in Section 6.

##  Appendix B Mapping the optically thin ionization rate to an EUV luminosity

In order to map the optically thin ionization rate at the planet to an
integrated EUV luminosity of the star, one needs a stellar spectrum. As the
interstellar medium extinguishes much EUV radiation from stars, both the
integrated EUV luminosity of planet-hosting stars and the shape of the
spectrum have significant uncertainties. This means that the range of
optically thin photoionization rates accounts for variations in both these
factors.

Assuming we have a "perfect" stellar spectrum, $f_{\nu}$, that has been
normalised to unity in the EUV range, we can map the integrated EUV luminosity
to the optically thin ionization rate as follows:

$\displaystyle\Gamma=L_{\text{EUV}}\int_{13.6~{}{\rm eV}}^{100~{}{\rm
eV}}\frac{f_{\nu}}{4\pi r^{2}h\nu}\sigma_{\nu}d\nu$ (52)

where $L_{\text{EUV}}$ is the luminosity integrated over the EUV range. As
$\sigma_{\nu}$ peaks at $13.6$ eV and decreases strongly as a function as a
function of frequency, $\propto\left(\frac{\bar{\nu}}{\nu}\right)^{3}$ (e.g.,
Spitzer, 1978), harder spectra lead to a smaller photoionization rate for the
same integrated EUV luminosity.