11institutetext: Anton Pannekoek Institute for Astronomy, University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
22institutetext: Center for Astrophysics — Harvard & Smithsonian, 60 Garden
Street, MS-16, Cambridge, MA 02138, USA

# Constraining planetary mass-loss rates by simulating Parker wind profiles
with Cloudy

D.C. Linssen E-mail<EMAIL_ADDRESS>A. Oklopčić 11 M. MacLeod 22

(Received XX-XX-XXXX; accepted YY-YY-YYYY)

Models of exoplanet atmospheres based on Parker-wind density and velocity
profiles are a common choice in fitting spectroscopic observations tracing
planetary atmospheric escape. Inferring atmospheric properties using these
models often results in a degeneracy between the temperature and the mass-loss
rate, and thus provides weak constraints on either parameter. We present a
framework that can partially resolve this degeneracy by placing more stringent
constraints on the expected thermospheric temperature. We use the
photoionization code Cloudy within an iterative scheme to compute the
temperature structure of a grid of 1D Parker wind models, including the
effects of radiative heating/cooling, as well as the hydrodynamic effects
(expansion cooling and heat advection). We constrain the parameter space by
identifying models that are not self-consistent through a comparison of the
simulated temperature in the He 10830 Å line-forming region to the temperature
assumed in creating the models. We demonstrate this procedure on models based
on HD 209458 b. By investigating the Parker wind models with an assumed
temperature between $4000$ and $12000$ K, and a mass-loss rate between
$10^{8}$ and $10^{11}$ g s-1, we are able to rule out a large portion of this
parameter space. Furthermore, we fit the models to previous observational data
and combine both constraints to find a preferred thermospheric temperature of
$T=8200^{+1200}_{-1100}$ K and a mass-loss rate of
$\dot{M}=10^{9.84^{+0.24}_{-0.27}}$ g s-1 assuming a fixed atmospheric
composition and no gas pressure confinement by the stellar wind. Using the
same procedure, we constrain the temperatures and mass-loss rates of WASP-69
b, WASP-52 b, HAT-P-11 b, HAT-P-18 b and WASP-107 b.

##  1 Introduction

Atmospheric escape is thought to significantly influence the evolution of
planets at short orbital distances, in particular for lower mass planets which
can lose a sizeable fraction of their mass (Owen, 2019). It has been proposed
that this process shapes planetary demographics, resulting in the observed hot
Neptune desert and evaporation valley (Owen & Wu, 2013; Fulton et al., 2017;
Ginzburg et al., 2018). However, planet formation (Lopez & Rice, 2018; Lee et
al., 2022), migration (Matsakos & Königl, 2016; Owen & Lai, 2018) and the
local environment (Kruijssen et al., 2020) may play a role in sculpting these
features in the population as well. Direct observational evidence of escaping
atmospheres has come in the form of strong transit absorption in the
Lyman-$\alpha$ line (e.g. Vidal-Madjar et al., 2003; Lecavelier des Etangs et
al., 2010; Lavie et al., 2017), the Balmer-$\alpha$ line (e.g. Yan & Henning,
2018; Yan et al., 2021) the helium line at 10830 Å (e.g. Spake et al., 2018;
Allart et al., 2018; Vissapragada et al., 2020) and (ionized) metal lines in
the UV (e.g. Linsky et al., 2010; Sing et al., 2019; Cubillos et al., 2020).

Isothermal Parker wind models (Parker 1958, see also Lamers & Cassinelli 1999)
are a common choice for modelling the upper atmospheric structure when fitting
transmission spectra, such as those obtained in observations of the helium
10830 Å line. Parker wind models parametrize the outflow in terms of the
temperature and mass-loss rate (Oklopčić & Hirata, 2018; Lampón et al., 2020;
dos Santos et al., 2022), enabling the retrieval of these parameters from
observations, while staying agnostic about the mechanism(s) responsible for
controlling the outflow.111An often-used complementary approach is based on
self-consistent forward models, which assume the photoevaporation scenario of
atmospheric escape to obtain a unique solution for a given set of parameters
(such as Murray-Clay et al., 2009; Salz et al., 2016; Wang & Dai, 2021a;
Caldiroli et al., 2021a).  However, studies that have tried to constrain
planetary mass-loss rates in this way have often found large degeneracies with
temperature, resulting in only loose constraints on both parameters (e.g.
Mansfield et al., 2018; Vissapragada et al., 2020; Palle et al., 2020; Lampón
et al., 2020, 2021; Paragas et al., 2021). In some cases, such a degeneracy
can be somewhat reduced by fitting the full line profile observed at high
resolution (Lampón et al., 2021; dos Santos et al., 2022). Since this is not
possible for low-resolution observations such as those with the Hubble Space
Telescope Wide Field Camera 3 (Spake et al., 2018; Mansfield et al., 2018), an
ultra-narrowband photometric filter (Vissapragada et al., 2020; Paragas et
al., 2021), or with the James Webb Space Telescope, a different approach to
reducing the degeneracy is required for these cases.

One promising way of reducing this degeneracy is observing multiple spectral
lines for the same planet and performing a joint fit. For example, Lampón et
al. (2020, 2021) constrained mass-loss rates and H/He ratios by combining
information from He 10830 Å and Ly$\alpha$ observations, but these analyses
may be complicated by the fact that Ly$\alpha$ transits do not primarily probe
mass-loss rates, as pointed out by Owen et al. (2021). Interpreting these
transmission signatures of the upper atmosphere will require detailed models
that produce accurate excited and ionized level populations, taking into
account non-local thermodynamic equilibrium (NLTE) processes (e.g. Young et
al., 2020). For the metal species absorbing in the UV, the atmospheric
temperature plays an important role, as collisional effects may contribute
significantly to their excitation and ionization state. With these goals in
mind, we develop a framework around the NLTE plasma simulation code
Cloudy222https://gitlab.nublado.org/cloudy/cloudy/-/wikis/home (Ferland et
al., 1998, 2017) to model the temperature structure of escaping atmospheres of
exoplanets.

In this work, we use Cloudy to simulate the temperature structure of 1D Parker
wind models. We thus combine the flexibility of parametrized density and
velocity profiles with a detailed treatment of the chemical and thermal state
of the outflow, instead of fixing an isothermal temperature. We show that the
resulting temperatures can be used to identify Parker wind models that are not
self-consistent. Restricting the ‘allowed’ parameter space in this way allows
for tighter observational constraints on the mass-loss rate, even when only a
single spectral line is observed. We demonstrate our method here on a sample
of planets with helium observations at 10830 Å, but it is not limited to this
specific line. Vissapragada et al. (2022a) recently pursued a similar goal
using a different approach. They constrained the self-consistent Parker wind
parameter space by analytically calculating the maximum mass-loss efficiency
of an energy-limited photoevaporative outflow.

The outline of this paper is as follows. In Sec. 2 we review the basic
assumptions underlying isothermal Parker wind models and how we can test them
to assess the model self-consistency. In Sec. 3 we present an algorithm to
model the temperature structure of upper exoplanet atmospheres with Cloudy. In
Sec. 4 we present the resulting temperature structure constraints for
different planets and combine these with fits to the observed line strengths
to better constrain their mass-loss rates. In Sec. 5 we compare our results
and methods to similar studies and discuss the advantages and potential
drawbacks of our modeling approach. We summarize in Sec. 6.

##  2 Assumptions of the Parker wind model

The three equations that describe a steady-state planetary outflow in 1D are
those of mass conservation (Eq. 1), momentum conservation (Eq. 2) and energy
conservation (Eq. 3):

$\dot{M}=4\pi r^{2}\rho v,$ (1)
$v\frac{dv}{dr}+\frac{1}{\rho}\frac{dP}{dr}+\frac{GM_{p}}{r^{2}}=0,$ (2)
$-\rho
v\frac{d}{dr}\Bigg{[}\frac{kT}{(\gamma-1)\mu}\Bigg{]}+\frac{kTv}{\mu}\frac{d\rho}{dr}+\Gamma+\Lambda=0.$
(3)

In these equations, $\dot{M}$ is the mass-loss rate, $r$ the distance from the
center of the planet, $\rho$ the density, $v$ the outflow velocity, $P$ the
gas pressure, $G$ the graviational constant, $M_{p}$ the planet mass, $M_{*}$
the stellar mass, $k$ the Boltzmann constant, $T$ the temperature, $\gamma$
the adiabatic index (5/3 for an ideal gas), $\mu$ the mean molecular weight
and $\Gamma$ and $\Lambda$ the radiative volumetric heating and cooling rates,
respectively. The gravitational influence of the star and the Coriolis force
are neglected in Eq. 2. Solving this set of equations gives the three main
parameters describing the outflow, $\rho(r)$, $v(r)$ and $T(r)$, but is a non-
trivial task as it involves finding the radiative heating and cooling rates,
which depend on $\rho$, $v$ and $T$. An example of a solution to this problem
was provided by Murray-Clay et al. (2009), under the assumption that hydrogen
photo-ionization and Lyman-$\alpha$ emission are the main contributors to
$\Gamma$ and $\Lambda$, respectively. Recently, Caldiroli et al. (2021b)
presented a more advanced solution that includes additional contributors such
as heating by helium ionization and cooling by bremsstrahlung and
recombination.

Isothermal Parker wind profiles simplify the problem by ignoring Eq. 3, taking
$T/\mu$ to be constant, and solving only Eqs. 1 and 2 for the density and
velocity structure (see e.g. Lamers & Cassinelli 1999 for the calculation).
Maintaining a constant $T/\mu$ as $r\rightarrow\infty$ requires energy input
as can be seen from the second term of Eq. 3, which is not physical at very
large altitudes. This is not a problem however at altitudes that are relevant
for launching the outflow (until at least the Roche radius), since Parker wind
structures are typically quite similar to self-consistent simulations such as
those of Salz et al. (2016), as shown in Oklopčić & Hirata (2018). Even though
the constant $T(r)/\mu(r)$ does not immediately imply an isothermal structure,
a fixed temperature $T_{0}$ can be assigned if a fixed value for $\mu(r)$ is
calculated. This $\bar{\mu}$ is usually taken as a weighted mean and requires
assuming a gas composition (see e.g. Oklopčić & Hirata, 2018; Lampón et al.,
2020). The simplification of the Parker wind model allows the freedom to
choose $T_{0}$ and $\dot{M}$, which helps in retrieving these parameters for a
given planet, as they can be adjusted to fit observational results. This is in
contrast to self-consistent forward models that provide one outflow solution,
and would have to alter their physical assumptions in order to match
observations that are different from the model prediction (e.g. Zhang et al.,
2022).

In principle, a Parker wind solution can be calculated for any value of
$T_{0}$ and $\dot{M}$ (as well as for different gas compositions), but not all
of these models are physical, which can manifest itself in various ways.
Vissapragada et al. (2022a) presented one such assessment of the self-
consistency of Parker wind profiles. Under the assumption of a maximally
efficient energy-limited outflow heated by photoionization, they checked
whether the total amount of heat absorbed by the atmosphere provided enough
energy to launch the wind, which allowed them to rule out a part of the
$T_{0}-\dot{M}$ parameter space for which this was not the case.

Here, we assess the self-consistency of Parker wind solutions through a
different route, focusing on the assumed temperature $T_{0}$. In particular,
we investigate whether the assumed temperature of a model is reasonable, given
its density and velocity structure and host star spectrum. We do so by
reintroducing Eq. 3 and solving it for $T(r)$, using $\rho$ and $v$ of the
Parker wind profile. The comparison between $T(r)$ and $T_{0}$ then acts as a
‘sanity check’ on the Parker wind assumption that Eq. 3 is well approximated
by the parametrized $T_{0}$. By quantifying the discrepancy between $T(r)$ and
$T_{0}$, we construct a $T_{0}-\dot{M}$ Parker wind parameter space of self-
consistent models. We note that a slightly more sophisticated approach would
be to compare the parametrized $T/\mu$ to the $T(r)/\mu(r)$ structure we get
from solving Eq. 3. In literature however, Parker wind fits to observations
are usually presented in terms of $T_{0}$ rather than $T/\mu$, and thus
benefit more from constraints on $T_{0}$. Furthermore, for a few test models
we found that the relative differences were similar for $T$ and $T/\mu$.

##  3 Modelling exoplanet atmospheres with Cloudy

We used Cloudy v17.02, a NLTE plasma-simulation code suited for a wide variety
of astrophysical conditions, to model Parker wind models. Cloudy assumes
hydrostatic equilibrium and simulates a 1D slab of material under irradiation
of a light source. In the context of our simulations, this means that given an
atmospheric density profile and incident stellar spectrum, Cloudy calculates
the (electron) temperature, ionization and excitation states of all of the
most relevant atomic species, as well as the line and continuum radiative
transfer and the associated heating and cooling rates, at each depth into the
cloud. Alternative to letting Cloudy solve for the electron temperature, the
code can also be run for a prescribed temperature profile, in which case it
may return potentially unequal heating and cooling rates. The chemical
composition of the plasma includes all elements up to atomic number 30 and can
be changed as desired. We used the default solar composition throughout our
simulations, but were forced to exclude the element calcium due to an error in
Cloudy. For a more detailed description of Cloudy’s main features, we refer to
Ferland et al. (2017) or Fossati et al. (2021), who provide a concise summary
in their Section 2.2.1.

###  3.1 Previous uses of Cloudy in exoplanet context

When using Cloudy to simulate escaping exoplanet atmospheres, the assumption
of hydrostatic equilibrium is not valid, and attention must be given to
treating hydrodynamic effects properly. Perhaps the most notable example of
this was provided by Salz et al. (2016), who coupled Cloudy to magneto-
hydrodynamics code PLUTO (Mignone et al., 2012) to create The PLUTO-CLOUDY
Interface (Salz et al., 2015). In their framework, PLUTO and Cloudy are
alternately executed. Cloudy provides accurate net radiative heating or
cooling rates, which are then passed on to PLUTO to calculate the next time
step in the hydrodynamic simulation. By coupling to PLUTO, hydrodynamic
effects on the temperature structure, such as heat advection and expansion
cooling, are included. Using their code, Salz et al. (2016) simulated
photoevaporative winds of 18 exoplanets in 1D. Their simulations are self-
consistent and physically elaborate, but computationally more expensive than
our approach, and therefore less suitable for exploring a large parameter
space when fitting models to observations.

Alternatively to using Cloudy’s heating and cooling rates to solve for the
temperature structure self-consistently, Turner et al. (2020) and Young et al.
(2020) used Cloudy with an assumed temperature-pressure profile to model the
transmission spectra of KELT-9 b and HD 209458 b, respectively. They projected
the one-dimensional T-P profile onto a spherically symmetric grid and ran
Cloudy simulations for stellar light rays at different impact parameters from
the center of the planet. Similarly, Fossati et al. (2020) used this
computational scheme with a family of T-P profiles to retrieve the atmospheric
structure of KELT-9 b from published observations of the H$\alpha$ and
H$\beta$ lines. Fossati et al. (2021) expanded on this work by forward
modelling the T-P profile of KELT-9 b, using Cloudy to model the upper
atmosphere ($P\lesssim 10^{-4}$ bar) self-consistently. They ran Cloudy
iteratively to solve for the temperature structure, which was then used to
update the pressure scale heights and density structure. Their calculations
were performed under the assumption of hydrostatic equilibrium, thus not
considering an outflowing atmosphere nor determining mass-loss rates.

###  3.2 Including hydrodynamic effects on the temperature structure

The temperature structure of the atmosphere is governed by the equation for
energy conservation (Eq. 3), which contains heating and cooling rates due to
various physical processes. The first term of this equation represents
advection: heating or cooling due to the bulk transport of the fluid. The
second term represents cooling due to adiabatic expansion. The third and
fourth term represent radiative heating and cooling, respectively, and are the
only sources included in Cloudy. We developed an iterative algorithm centered
around Cloudy in order to solve for the temperature structure while including
advection and expansion heating and cooling rates. For a given density and
velocity structure, the procedure is as follows.

We start by running a Cloudy simulation at a constant temperature or
previously solved temperature structure of a similar Parker wind profile (the
initially assumed temperature does not affect the final temperature
structure). In this fixed temperature mode, Cloudy returns the radiative
heating and cooling rates as a function of altitude ($\Gamma$ and $\Lambda$ in
Eq. 3), as well as the ionization structure. We use the latter to calculate
the mean molecular weight structure, considering the ionization fractions of
hydrogen and helium and neglecting all other elements, which are included in
the Cloudy simulations at solar metallicity. Together with the given density
and velocity structures, we calculate the advection and expansion terms of Eq.
3.

Adding together all heating and cooling rates separately results in a ratio of
total heating to cooling rate as a function of altitude, which we use to
construct a new temperature structure. In regions where this ratio is $>1$, we
construct a new temperature structure according to
${T_{new}=T_{old}(1+0.3\mathrm{log}_{10}(\mathcal{H}/\mathcal{C}))}$, where
$\mathcal{H}$ and $\mathcal{C}$ are the total heating and cooling rates,
respectively. In regions where the ratio is $<1$, we construct a new
temperature structure according to
${T_{new}=T_{old}/(1-0.3\mathrm{log}_{10}(\mathcal{H}/\mathcal{C}))}$. We then
pass the new temperature structure on to Cloudy to get $\Gamma$ and $\Lambda$
for the next iteration, and this process is repeated until we reach a
converged temperature structure for which the total heating and cooling rate
are in balance. The specific functional dependence of $T_{new}$ on the
$\mathcal{H}/\mathcal{C}$ ratio does not have a physical motivation, rather,
we found that this expression ensures relatively fast convergence of $T$. We
set a convergence threshold of
$\nicefrac{{1}}{{1.1}}<\mathcal{H}/\mathcal{C}<1.1$, which for a typical
temperature structure means $\lesssim 50$ K convergence. Higher precision is
usually limited by round-off errors in Cloudy’s interpolation of the given
density and temperature structures onto its internal grid, combined with the
fact that this internal grid changes throughout successive iterations of our
algorithm, resulting in slight variations of the advection rate. Simulating a
temperature structure in this way usually takes between 2 and 10 iterations.

In regions where advection is the dominant heating process, we are able to
speed up the convergence. Instead of constructing a new temperature structure
based on the $\mathcal{H}/\mathcal{C}$ ratio as described above, we construct
the new temperature structure by solving Eq. 3 for $T$, using the $\mu$,
$\Gamma$ and $\Lambda$ of the previous iteration. Since $\Gamma$ and $\Lambda$
depend on the temperature, this approach only works well when $\Gamma$ and
$\Lambda$ of the newly constructed temperature structure are comparable to
those of the last iteration, or if they are small compared to the advection
and expansion rates. In an advection-dominated regime, the latter is indeed
the case and we use this alternative way of solving for $T$, usually requiring
on the order of $\sim$2 iterations to converge. In radiation-dominated
regimes, it is exactly the dependence of $\Gamma$ and $\Lambda$ on $T$ that
determines the temperature, and using $\Gamma$ and $\Lambda$ of the last
iteration to construct the new temperature structure with Eq. 3 does not work,
which is why we use the $\mathcal{H}/\mathcal{C}$ method described above.

We tested our algorithm on the model of HD 209458 b provided by Salz et al.
(2016). They ran simulations with The PLUTO-CLOUDY Interface (Salz et al.,
2015, see also Sec. 3.1) of a photoevaporation driven escaping atmosphere.
Their use of a hydrodynamics code ensures proper treatment of hydrodynamic
thermal effects, and these simulations thus act as a useful validation of our
code. We extracted the spectral energy distribution (SED) of HD 209458 from
their Fig. 1 using the
WebPlotDigitizer333https://automeris.io/WebPlotDigitizer/index.html software
package (Ankit, 2020). Using this SED in combination with the density and
velocity structure shown in their Fig. 5, we converged on the temperature
structure shown in Fig. 1444A reproduction package for all figures of this
manuscript is available at doi:10.5281/zenodo.6798207.

Figure 1: Comparison of simulated temperature structures for the density
profile of HD 209458 b from Salz et al. (2016). The orange dashed line shows
the density structure (read from the right y-axis). Although not shown here,
the qualitative behavior of the different heating and cooling rates as a
function of altitude (such as shown for a different model in the bottom panel
of Fig. 2) is similar to the results of Salz et al. (2016), shown in their
Fig. 5.

Our results agree quite well, but there is a constant $\sim 1000$ K offset
between the simulations. Part of this difference can be explained by Cloudy’s
wind advection command, which was used by Salz et al. (2016). This command
includes the advection of ion species, such that neutral atoms are transported
upwards with the outflow velocity, leading to higher ionization heating in the
upper atmosphere. The command does not include the advection of heat (the
first term of Eq. 3), however. In our simulations, we did not use this
functionality as it slows down Cloudy very significantly. We enabled it
temporarily and observed an increase in temperature but were still unable to
perfectly match the reported temperature structure of Salz et al. (2016).
Changing to the older Cloudy version v13.01 that is used by The PLUTO-CLOUDY
Interface also made a minimal difference. We could therefore not find the
cause for the remaining discrepancy between the temperature structures.

##  4 Results

###  4.1 Two types of temperature structure

We ran our code on a suite of Parker wind models for the well-studied
exoplanet HD 209458 b, a hot Jupiter orbiting a G0 star. We used a gas
composition of 90% hydrogen and 10% helium by number, a radius of
$R_{p}=1.35R_{J}$ and a mass of $M_{p}=0.71M_{J}$ to create the models. We
explored mass-loss rates between $10^{8}<\dot{M}<10^{11}$ g s-1 and
temperatures between $4000<T_{0}<12000$ K. We used the SED of Salz et al.
(2016), shown in their Fig. 1, with a solar gas composition to run the Cloudy
simulations. Fig. 2 shows the result for the Parker wind profile with
$T_{0}=8000$ K, $\dot{M}=10^{10}$ g s-1. The bottom panel shows the heating
and cooling rates due to different physical processes. The qualitative
behavior of these rates is the same across all calculated Parker wind models
for this planet: radiative heating and cooling determine the temperature
structure at lower altitudes, while advection heating and expansion cooling
set the temperature at higher altitudes. Therefore, the temperature at lower
altitudes does not depend sensitively on the outflow velocity.

Figure 2: Top: Velocity (black) and density (dashed orange, read from the
right y-axis) structure for the Parker wind profile of HD 209458 b with
$T_{0}=8000$ K, $\dot{M}=10^{10}$ g s-1. Middle: Simulated temperature
structure (black) and optical depth (dashed orange, read from the right
y-axis) in the major component of the metastable helium line (dashed line).
Bottom: Relevant heating and cooling rates, shown in red and blue,
respectively. We observe that at lower altitudes ($\lesssim 2R_{p}$),
radiative heating and cooling dominate, while at higher altitudes ($\gtrsim
3R_{p}$), advection heating and expansion cooling are the dominant thermal
processes.

This behavior changes for planets with lower gravitational potential
($|\phi|=GM_{p}/R_{p}$). We ran a grid of Parker wind models for WASP-69 b, a
Saturn-mass planet orbiting a K5 star. We adopted a mass of $M_{p}=0.26M_{J}$,
a radius of $R_{p}=1.057R_{J}$ and used the version 2.2 SED of the K6 star
HD85512 from the MUSCLES survey (France et al., 2016; Youngblood et al., 2016;
P. Loyd et al., 2016). We ran models with mass-loss rates between
$10^{9}<\dot{M}<10^{11}$ g s-1 and temperatures between $5000<T_{0}<9000$ K.
Fig. 3 shows a typical temperature structure, in this case for $T_{0}=6000$ K,
$\dot{M}=10^{10}$ g s-1. WASP-69 b’s lower gravitational potential compared to
HD 209458 b (see Table 1) results in a qualitatively different atmospheric
structure. For this type of planet, expansion cooling is the dominant coolant
throughout the atmosphere while the temperature decreases monotonically with
altitude for most models. A consequence of this is that the velocity
influences the thermal structure even at the base of the atmosphere, as the
expansion cooling rate depends on velocity (see Eq. 3).

Figure 3: Similar to Fig. 2, but for the Parker wind profile of WASP-69 b with
$T_{0}=6000$ K, $\dot{M}=10^{10}$ g s-1. Here, we observe that at nearly all
altitudes, expansion cooling is the dominant cooling agent.

When we turn to other planets (Table 1 lists our sample, although we do not
show their temperature structures), we find a similar behavior: planets with
lower gravitational potential like HAT-P-18 b, HAT-P-11 b and WASP-107 b are
dominated by expansion cooling throughout the atmosphere and showcase a
monotonically decreasing temperature structure. WASP-52 b falls in the
transition regime, with some Parker wind models dominated by expansion cooling
and others by radiative cooling at lower altitudes. We thus identify two types
of upper atmospheres. For high-gravity gaseous planets (loosely taken as
$|\phi|\gtrsim 10^{12.82}$ erg g-1 here), radiative heating is balanced by
radiative cooling and the thermospheric temperature increases with altitude
until expansion cooling takes over and the exospheric temperature decreases
with altitude. For lower gravity gaseous planets, expansion cooling always
dominates over radiative cooling, resulting in generally cooler temperatures
that decrease with altitude. A similar difference between the dominant cooling
process in planets with high and low gravitational potential was identified by
Salz et al. (2016) and Caldiroli et al. (2021b).

###  4.2 Self-consistency of Parker wind profile temperatures

For many models, the temperature structure we converged on with Cloudy does
not agree with the constant temperature assumed to create the model. In
principle, complete agreement is never expected, as the temperature structures
found with Cloudy are not isothermal. Still, some models show reasonable
overlap between the range of temperatures returned by Cloudy and the assumed
constant value in the atmospheric regions that we are interested in. If this
is not the case, such a discrepancy indicates that the Parker wind profile is
not self-consistent and therefore not a physical model solution. While fitting
observations, models that are not self-consistent can be excluded from
consideration, which can help to put narrower constraints on the mass-loss
rate and temperature.

The criterion to assess and ‘quantify’ the self-consistency of a Parker wind
profile depends on the way in which Cloudy’s temperature structure is compared
to the assumed isothermal profile. In principle, this choice is arbitrary and
can be adjusted to fit the observations that we are aiming to interpret. In
other words, one has the freedom to define a characteristic temperature of the
converged structure that is suitable to the science case at hand, and compare
it to the assumed temperature.

In this work, we aim to aid the interpretation of metastable helium
observations and we defined the characteristic temperature as the average in
the He 10830 Å line-forming region, since it is this temperature that is
probed by the observations. One way to calculate this (weighted) average
temperature is to weigh the temperature in each bin by the contribution of
that bin to the optical depth at 10830 Å. The problem with this definition is
that higher altitude bins are relatively down-weighed, as their higher radial
velocity shifts the line away from the rest-frame wavelength such that the
optical depth contribution at 10830 Å (line center) is lower. To avoid this,
we instead opted to weigh the temperature of each bin by that bin’s metastable
helium column density, which is independent of velocity but still directly
related to the optical depth contribution. Thus using the metastable helium
column densities $N_{He*}$ as given by our Cloudy simulations, we calculated
the characteristic temperature

$T_{He}=\frac{\sum T(r)\cdot N_{He*}(r)}{\sum N_{He*}(r)}.$ (4)

The difference between the temperature of the He 10830 Å line-forming region
$T_{He}$ and the temperature assumed to create the Parker wind profile $T_{0}$
then provides a quantifier for the self-consistency of each Parker wind
profile. Following this definition, the particular Parker wind profiles
presented in Figs 2 and 3 are considered self-consistent, as $T_{He}\approx
T_{0}$. As a measure for the spread of $T(r)$ around the characteristic
temperature $T_{He}$, we calculated the weighted standard deviation

$\sigma_{T}^{2}=\frac{\sum\big{(}T(r)-T_{He}\big{)}^{2}\cdot N_{He*}(r)}{\sum
N_{He*}(r)}.$ (5)

This quantity indicates the range of atmospheric temperatures that would be
probed by the metastable helium line, and thus can be used to judge the
significance of the difference between $T_{He}$ and $T_{0}$.

###  4.3 Mass-loss rate constraints

For each of our simulated Parker wind profiles for HD 209458 b, we calculated
$T_{He}-T_{0}$ and $\sigma_{T}$. The results are shown in the left panels of
Fig. 4. We find a rather narrow range of self-consistent Parker wind profiles
for which $T_{He}\approx T_{0}$, at an almost constant temperature of $T\sim
8000$ K for the range of mass-loss rates explored. Our code always points us
towards this temperature: models with $T_{0}>8000$ K result in helium line
formation temperatures lower than $T_{0}$, as visible by the blue part in the
bottom left panel of Fig. 4, indicating that Cloudy prefers a cooler model.
Vice versa, the red region in the figure for $T_{0}<8000$ K reveals that
Cloudy prefers hotter models.

Figure 4: Top left: helium line-forming temperatures $T_{He}$ (Eq. 4) for
simulated temperature structures of HD 209458 b Parker wind models. The
colormap linearly interpolates on the values of the models that were run
(indicated by grey dots). The grey region in the upper left corner could not
be simulated as Cloudy could not handle the high density of these models well.
Bottom left: Difference between $T_{He}$ (top left panel) and the assumed
isothermal value $T_{0}$. The dashed and dotted lines indicate the $1\sigma$
and $2\sigma$ standard deviations of the self-consistent region based on Eq.
5, respectively. Top right: Compatibility between the He 10830 Å EWs of the
Parker wind models and the observed value from Alonso-Floriano et al. (2019).
The dark blue region indicates good agreement. Bottom right: The orange line
and dark and light areas are the self-consistent models, 1$\sigma$ and
2$\sigma$ contours from the bottom left panel, respectively, which we use as
our prior. The blue line and dark and light areas are the best-fit models,
1$\sigma$ and 2$\sigma$ contours from the top right panel, respectively, which
are our likelihood values. The 1$\sigma$ contour of the posterior is shown in
black. The 1$\sigma$ posteriors of the 3x and 1/3x EUV flux simulations are
shown in grey.

We next fit the Parker wind models to the He 10830 Å observations of HD 209458
b by Alonso-Floriano et al. (2019). They report an average absorption of 0.44%
in a 0.915 Å bandpass centered on the helium line. In order to estimate the
uncertainty on this value, we extracted the light curve in this bandpass from
their Fig. 6 using WebPlotDigitizer (Ankit, 2020) and followed a bootstrap
analysis similar to that described in Sec. 3.4 of Salz et al. (2018). We
randomly drew half of the points in the light curve between second and third
contact and resampled those points from normal distributions with a standard
deviation of 0.167%, which is approximately the average error on the light
curve values. We calculated the mean of the new light curve and repeated this
procedure 5000 times, resulting in a sample of average absorption values. The
standard deviation of this sample was 0.0596%, which we used as the error.
Thus, the average absorption of $0.44\pm 0.06$% in a 0.915 Å bandpass could be
translated to an equivalent width of EW=$4.03\pm 0.55$ mÅ.

For every simulated Parker wind profile, we calculate the transmission
spectrum following the procedure described in Sec. 3.4 of Oklopčić & Hirata
(2018), using the metastable helium number density and temperature profile as
given by the Cloudy simulations. We made mid-transit spectra at the planet’s
impact parameter while including stellar limb darkening using the parameters
obtained from
EXOFAST555https://astroutils.astronomy.osu.edu/exofast/limbdark.shtml (Eastman
et al., 2013). We integrated these helium transmission spectra to obtain the
EW and compared it to the observed value. This approach means we did not fit
the resolved spectrum because we want to develop a framework that works on
both low- and high-resolution data. The fact that we fit the average in-
transit absorption with a synthetic mid-transit spectrum can be justified in
this case because of the relatively large scatter in the light curve, visible
in Fig. 6 of Alonso-Floriano et al. (2019). Using the full in-transit light
curve gives a higher S/N absorption value than using the few points around
mid-transit. Furthermore, the in-transit light curve appears quite flat
despite the scatter and does not clearly reveal asymmetric absorption
indicative of a planetary tail or center-to-limb variations, which would
otherwise prompt us to consider fitting only the mid-transit observations. The
results of the fits are shown in the top right panel of Fig. 4. A degeneracy
between the mass-loss rate and temperature of good-fitting models is visible,
which was also found by Lampón et al. (2020) who performed similar
calculations.

We then combined the observational and theoretical constraints into one joint
constraint on the mass-loss rate and temperature. We do this in a Bayesian
framework, treating the measure for self-consistency of the Parker wind model
(bottom left panel of Fig. 4) as our prior and the goodness-of-fit to the
observations (top right panel of Fig. 4) as the likelihood. For the prior, we
convert the temperature difference $(T_{He}-T_{0})$ to a p-value assuming it
is normally distributed with standard deviation $\sigma_{T}$. For example, a
temperature difference of $T_{He}-T_{0}=2\sigma_{T}$ would be considered a
2$\sigma$ discrepancy and thus translate to a prior p-value of $\sim$0.05. For
the likelihood, we convert the offset between modeled and observed EW to a
p-value using the normally distributed observed error bar. Applying Bayes’
theorem then gives the posterior distribution, which is shown in the bottom
right panel of Fig. 4. As visible in the figure, the constraints on self-
consistent models complement the fits to the observations extremely well.

To estimate our modeling uncertainties, we repeated our calculations for host
star spectra with different levels of EUV flux. The EUV part of the spectrum
strongly influences the temperature because it is the main driver of radiative
heating of the atmosphere through photoionization. Furthermore, the metastable
helium number density depends on the EUV flux, since these photons populate
the metastable helium state through ionization of ground-state atoms
(Oklopčić, 2019). Despite its importance, accurate stellar EUV spectra are
notably hard to constrain through observations and thus not readily available.
EUV spectra reconstructed through scaling relations (such as those used in
this work) have typical uncertainties of one order-of-magnitude (Youngblood et
al., 2019; Drake et al., 2020), which makes them a large source of uncertainty
in our modeling efforts. In order to explore the effect of the EUV uncertainty
on our constraints, we performed simulations throughout the $T_{0}-\dot{M}$
Parker wind parameter space with 3x and 1/3x flux at $\lambda<1000$ Å. We
followed the same steps as in the ‘fiducial’ flux case to calculate the
posterior distribution for both the high and low EUV flux case. The $1\sigma$
contours of these two posteriors are shown by grey dashed lines in the bottom
right panel of Fig. 4. In total then, we find a temperature of
$T=8200^{+1200}_{-1100}$ K and a mass-loss rate of
$\dot{M}=10^{9.84^{+0.24}_{-0.27}}$ g s-1, where we quote the fiducial
simulations for our central value and the $1\sigma$ contours of the scaled EUV
flux simulations for the upper and lower errors.

We repeated our analysis for five other planets. Table 1 lists the input
parameters, observed helium signal strengths and inferred temperature and
mass-loss rate constraints for these planets, while Fig. 5 shows the posterior
distributions. We chose our sample based on the availability of published
metastable helium observations. The WASP-69 b observations were performed by
Vissapragada et al. (2020) and reanalyzed by Vissapragada et al. (2022b), who
report an excess absorption of $0.512\pm 0.049$% in their narrowband filter.
We convolved our helium model spectra with the transmission profile of this
filter (S. Vissapragada, private communication) to obtain the excess
absorption and compared these against the observed value. The WASP-52 b
observations were performed by Kirk et al. (2022). We extracted the EW from
their Fig. 9, which shows a mean absorption of $1.65\pm 0.25$% in a 2.44 Å
bandpass, indicating an EW of $40.3\pm 6.1$ mÅ. The HAT-P-11 b observations
were performed by Mansfield et al. (2018). They report a white-light transit
depth of $3371\pm 15$ ppm and a $3560\pm 34$ ppm transit depth in the
10809-10858 Å channel. Adding the uncertainties in quadrature, we obtain
EW=$9.26\pm 1.82$ mÅ. The HAT-P-18 b observations were performed by Paragas et
al. (2021) and reanalyzed by Vissapragada et al. (2022b), who report $0.70\pm
0.16$% excess absorption in the narrowband filter. The WASP-107 b observations
were performed by Spake et al. (2018), who report EW=$48.02\pm 10.78$ mÅ. We
note that for WASP-69 b, HAT-P-11 b and WASP-107 b, additional spectrally
resolved observations exist (Nortmann et al., 2018; Allart et al., 2018, 2019;
Kirk et al., 2020; Spake et al., 2021), but they are generally consistent with
the unresolved observations, so we took this opportunity to demonstrate the
ability to constrain mass-loss rates from unresolved observations. With the
exception of HD 209458, the SEDs of host stars in our sample are not available
in the literature. Therefore, for each host star we used the closest spectral
type SED from the MUSCLES survey (France et al., 2016; Youngblood et al.,
2016; P. Loyd et al., 2016), listed in Table 1. For HAT-P-11, which is a K4
star, we used an averaged MUSCLES spectrum of the K2 star HD 85512 and K6 star
$\epsilon$ Eridani. We do not know how well all of these SEDs describe the
true SEDs of our planets’ host stars, especially in the EUV part of the
spectrum.

Table 1: Model parameters, helium line strengths extracted from observational
papers, and the constrained temperatures and mass-loss rates for investigated
exoplanets. The latter were inferred from Fig. 5.

Planet | $R$ | $M$ | $\mathrm{log_{10}}(|\phi|)$ | $a$ | SED | He $\lambda$10830 signal | $T$ |  $\mathrm{log_{10}}(\dot{M}$)
---|---|---|---|---|---|---|---|---
| $[R_{J}]$ | $[M_{J}]$ | [cgs] | [AU] |  |  | [K] | [$\mathrm{log_{10}}$(g s-1)]
HD 209458 b | 1.39 | 0.73 | 12.98 | 0.04747 | Salz et al. (2016) |  $4.03\pm 0.55$ mÅ 11$1$EW from Alonso-Floriano et al. 2019 | $8200^{+1200}_{-1100}$ | $9.84^{+0.24}_{-0.27}$
WASP-69 b | 1.057 | 0.26 | 12.65 | 0.04525 | HD 85512 (K6)mm$m$From the MUSCLES survey, version 2.2 (France et al., 2016; Youngblood et al., 2016; P. Loyd et al., 2016) |  $0.512\pm 0.049$ % 22$2$Excess absorption from Vissapragada et al. 2022b | $5300^{+1300}_{-300\;u}$ | $9.91^{+0.19}_{-0.12\;u}$
WASP-52 b | 1.27 | 0.46 | 12.82 | 0.0272 |  $\epsilon$ Eridani (K2)mm$m$From the MUSCLES survey, version 2.2 (France et al., 2016; Youngblood et al., 2016; P. Loyd et al., 2016) |  $40.3\pm 6.1$ mÅ 33$3$EW from Kirk et al. 2022 | $9000^{+1800}_{-1500}$ | $11.06^{+0.26}_{-0.25}$
HAT-P-11 b | 0.389 | 0.0736 | 12.53 | 0.053 | Avg. K2 & K6aa$a$Averaged MUSCLES spectra of $\epsilon$ Eridani (K2) and HD 85512 (K6) |  $9.26\pm 1.82$ mÅ 44$4$EW from Mansfield et al. 2018 | $4500^{+1300}_{-1100}$ | $9.70^{+0.29}_{-0.15}$
HAT-P-18 b | 0.995 | 0.197 | 12.55 | 0.0559 |  $\epsilon$ Eridani (K2)mm$m$From the MUSCLES survey, version 2.2 (France et al., 2016; Youngblood et al., 2016; P. Loyd et al., 2016) |  $0.70\pm 0.16$ % 55$5$Excess absorption from Vissapragada et al. 2022b | $5400^{+1300}_{-1200}$ | $10.33^{+0.29}_{-0.24}$
WASP-107 b | 0.94 | 0.096 | 12.27 | 0.055 | HD 85512 (K6)mm$m$From the MUSCLES survey, version 2.2 (France et al., 2016; Youngblood et al., 2016; P. Loyd et al., 2016) |  $48.02\pm 10.78$ mÅ 66$6$EW from Spake et al. 2018 | $3400^{+1400}_{-400\;u}$ | $9.96^{+0.23}_{-0.19\;u}$

666 $u$$u$footnotetext: Underestimated uncertainty due to exceeding our
modeled parameter space

Figure 5: Similar to the bottom right panel of Fig. 4, but also for a suite of
other planets with reported metastable helium observations. The observations
for HD 209458 b are by Alonso-Floriano et al. (2019), for WASP-69 b by
Vissapragada et al. (2022b), for WASP-52 b by Kirk et al. (2022), for HAT-P-11
b by Mansfield et al. (2018), for HAT-P-18 b by Paragas et al. (2021), and for
WASP-107 b by Spake et al. (2018). The constrained mass-loss rates and
temperatures are listed in Table 1.

##  5 Discussion

###  5.1 Comparison with other mass-loss rate estimates

For HD 209458 b, we found a mass-loss rate of
$\dot{M}=10^{9.84^{+0.24}_{-0.27}}$ g s-1 and a temperature of
$T=8200^{+1200}_{-1100}$ K. This is a 2$\sigma$ difference with the mass-loss
rate results of Salz et al. (2016), who found $\dot{M}=10^{10.27}$ g s-1 and a
peak temperature of $T=9100$ K for simulations of an photoevaporative wind.
Murray-Clay et al. (2009) employed a theoretical model of a hydrodynamically
escaping atmosphere driven by photoevaporation and found a peak temperature of
$T\sim 8600$ K and a maximum mass-loss rate of $\dot{M}\leq 10^{10.5}$ g s-1,
dependent on the UV flux of the host star. Both of these works estimated the
mass-loss rate from simulations, while our constraint is based on observations
in the metastable helium line. When comparing to other constraints on the
mass-loss rate based on observations, we find that our estimates are slightly
below the lower limit of $\dot{M}\geq 10^{10}$ g s-1 from Vidal-Madjar et al.
(2003), who observed HD 209458 b in the Ly$\alpha$ line. Koskinen et al.
(2010) later reanalyzed these and other observations in the O I triplet at
1304 Å with a parametrized atmospheric model, treating the thermospheric
temperature as a free parameter. Their best-fitting models indicate a
thermospheric temperature of $T=8000-11000$ K and a mass-loss rate of
$\dot{M}=10^{10-11}$ g s-1, showing reasonable agreement with our results.

For HAT-P-11 b, we found $\dot{M}=10^{9.70^{+0.29}_{-0.15}}$ g s-1 and
$T=4500^{+1300}_{-1100}$ K. This is a 2$\sigma$ difference with
$\dot{M}=10^{10.29}$ g s-1 as reported by Salz et al. (2016). It is also in
disagreement with dos Santos et al. (2022), who used the p-winds code to fit
the high-resolution observations of the helium line by Allart et al. (2018)
with isothermal Parker wind models and found $\dot{M}=10^{10.40\pm 0.12}$ g
s-1 and $T=7200\pm 700$ K. Their temperature is 2$\sigma$ higher than ours,
resulting in a higher mass-loss rate due to the degeneracy between these
parameters. The discrepancy between the temperature inferred by p-winds and
that predicted by Cloudy may be caused by inaccurate assumptions of either
model. Cloudy might underestimate the temperature of HAT-P-11 b due to our
adopted stellar spectrum (we use a template SED instead of the observed
HAT-P-11 spectrum, which is not available) and/or atmospheric composition. On
the other hand, p-winds might overestimate the temperature of HAT-P-11 b by
fitting the line width with a model that does not include some line-broadening
mechanisms which may be important, such as atmospheric circulation. These two
codes use fundamentally different approaches to constrain the atmospheric
temperature; hence, they are highly complementary. Understanding the cause of
this discrepancy may ultimately shed new light on atmospheric properties
and/or dynamics of HAT-P-11 b; however, this task is beyond the scope of this
paper.

For WASP-69 b, we found $\dot{M}=10^{9.91^{+0.19}_{-0.12}}$ g s-1 and
$T=5300^{+1300}_{-300}$ K, but the lower error bars are limited by the
parameter space we ran and is in reality higher than 0.12 dex and 300 K
(visible in the posteriors in Fig. 5). The mass-loss rate we found is
inconsistent with the value of $\dot{M}=10^{11.0}$ g s-1 that Wang & Dai
(2021a) find while fitting the metastable helium line with their 3D
hydrodynamic models. For WASP-107 b, our result of
$\dot{M}=10^{9.96^{+0.23}_{-0.19}}$ g s-1 and $T=3400^{+1400}_{-400}$ K (with
underestimated lower bounds, see Fig. 5) is again not consistent with the
result of Wang & Dai (2021b), who find $\dot{M}=10^{11.3}$ g s-1 with their 3D
hydrodynamic model fits to the helium line. This discrepancy is possibly
caused by the fact that their outflow is always hotter than $T\gtrsim$7000 K
and non-symmetric due to the inclusion of a stellar wind. For WASP-69 b, the
model of Wang & Dai (2021a) has temperatures similar to our constrained value
and the outflow geometry appears rather symmetric, so it is less clear where
the discrepancy stems from. Finally, our self-consistent Parker wind parameter
spaces for HAT-P-11 b, HAT-P-18 b and WASP-69 b are consistent with those
identified in Vissapragada et al. (2022a).

For the observational He 10830 Å data of HD 209458 b, WASP-69 b, HAT-P-11 b
and HAT-P-18 b, line fitting analyses using Parker wind profiles have already
been published (Lampón et al., 2020; Vissapragada et al., 2020; Mansfield et
al., 2018; Paragas et al., 2021). In these works, the assumed constant
temperature $T_{0}$ was used to calculate synthetic spectra. In Appendix A, we
investigate the effect that using this temperature has on the resulting
constraints and we compare our results to the literature. An interesting
finding is that for the low-gravity planets, the blue curves of Fig. 5, which
show the fits to the observations, are much ‘flatter’ than the blue curves of
Fig. 7. This means that using the non-constant temperature structures results
in tighter constraints on the mass-loss rate compared to using the isothermal
profiles.

###  5.2 Comparison to other atmospheric escape modelling efforts

The main advantage of our simulations lies in their applicability to
interpreting observations. Parker wind profiles parametrize the planetary wind
in two main parameters (assuming fixed chemical composition): temperature and
mass-loss rate, allowing us to explore a range of possible outflow profiles
without making assumptions on the physical processes driving the escape (e.g.
photoevaporation versus core-powered mass-loss). This is supported by the
relatively short run time of our code, roughly 15 CPU minutes to converge on
the temperature structure of one Parker wind profile.

Another advantage is that our simulations are NLTE through the use of Cloudy.
The importance of NLTE effects was stressed by Young et al. (2020), who found
that some spectral lines were 40% deeper than in the LTE case. In addition to
this, Cloudy guarantees a detailed treatment of radiative processes. An
example of one such process that is rarely included in simulations that treat
radiative heating/cooling analytically rather than with a dedicated code such
as Cloudy, is cooling by metal species. This was recently shown to be
important in the upper atmosphere of KELT-9 b by Fossati et al. (2021). They
also used Cloudy to solve for the temperature structure but under the
assumption of hydrostatic equilibrium, so our simulations differ by being
applicable to outflowing atmospheres. Figs. 2 and 3 show hydrodynamic effects
to be the dominant heating and cooling agents at high altitudes for both high-
and low-gravity planets, and neglecting them would lead to an overestimation
of the exospheric temperature. Finally, Cloudy supports the use of a full SED
shape instead of representing the incident spectrum with a few flux bins, the
latter being a common approach in atmospheric escape studies in literature
(e.g. Murray-Clay et al., 2009; Shaikhislamov et al., 2020; Caldiroli et al.,
2021a).

A potential drawback of our model is that we extrapolate the 1D sub-stellar
solution to a 3D spherically symmetric structure to calculate the synthetic
spectrum. However, the sub-stellar solution may not be representative of the
part of the atmosphere that is probed with transmission spectroscopy, even if
the outflow itself were spherically symmetric (which is most likely not the
case).

###  5.3 Using different characteristic temperatures

In our Parker wind models, we compared the assumed temperature to the mean
temperature in the helium line-forming region (Eq. 4), such that self-
consistent profiles could be defined and found. This choice was specifically
tailored towards the planets we investigated, as all of them have published
observations in the He 10830 Å line, allowing for a meaningful comparison of
the temperatures. The helium line is particularly suited for constraining the
planetary outflow models, as the line-forming region is also the region where
the bulk of the wind acceleration takes place. This can be seen by comparing
the optical depth and velocity profiles in Figs. 2 and 3. However, our
approach is in not limited to planets that have been observed in the helium
line. For planets with other observed atmospheric escape tracers, we could
define a corresponding characteristic mean temperature weighed by the column
density of the lower energy state.

In some cases, it might even be desirable to compare the assumed temperature
to some characteristic value of the simulated temperature structure that is
not linked to any specific line. Examples are the temperature at some fixed
altitude, or the peak temperature $T_{p}$. The latter ensures an unambiguous
upper limit to the temperature of a planet, since $T_{p}<T_{0}$ implies that
the whole atmosphere is inconsistent with the assumed temperature. We also
tested this criterion for HD 209458 b and again found very similar temperature
constraints to those of Fig. 4. This implies that for this planet, the helium
line forms in the hottest region of the atmosphere, which is indeed what we
see in Fig. 2.

###  5.4 Applicability to lower gravity planets

For the class of planets with lower gravitational potential presented in Sec.
4.1, we lose some of the advantages that Parker wind profiles provide, as the
dependence of the temperature on expansion cooling and, by extension, the
outflow velocity at the base of the atmosphere exposes a stronger dependence
on the mechanism by which the wind is launched. This is not the case for the
class of higher gravity planets, where the resulting temperature structure
does not depend on the velocity assumed at the base of the atmosphere. For the
model of WASP-69 b shown in Fig. 3, the velocity at $1R_{p}$ is $v=540$ m s-1,
while the velocity of the model of HD 209458 b shown in Fig. 2 is $v=10$ m s-1
at $1R_{p}$. This illustrates that the Parker wind model of WASP-69 b does not
transition into a hydrostatic atmosphere at low altitudes. Hydrodynamic
simulations such as those of Salz et al. (2016) and Caldiroli et al. (2021a)
show a more physical picture, where even for low-gravity planets, the velocity
approaches $v\sim 0$ at $\lesssim 1.1R_{p}$. Nevertheless, our method can
provide parameter space constraints to aid in data interpretation because the
simulated temperature structures should still hold for the given Parker wind
profiles. However, because in the low gravitational potential regime Parker
winds do not correctly model the lower altitudes, the contribution of
expansion cooling may be too strong. Therefore, especially for spectral lines
originating at lower altitudes, the temperature we constrain is potentially
too low, in turn also resulting in lower mass-loss rate estimates.

###  5.5 Energy-limited mass-loss efficiencies

An often-used way of estimating planetary mass-loss rates is to assume that
some given fraction of the incoming high-energy flux is converted into
expansion work that drives the atmospheric mass loss. This mechanism has been
dubbed ‘energy-limited’ escape and dictates that the mass-loss rate depends
linearly on the incoming XUV flux (e.g. Owen, 2019):

$\dot{M}=\eta\frac{\pi R_{p}^{3}F_{XUV}}{GM_{p}K},$ (6)

where $\eta$ is the heating efficiency and

$K=1-\frac{3}{2}\frac{R_{p}}{R_{Roche}}+\frac{1}{2}\Bigg{(}\frac{R_{p}}{R_{Roche}}\Bigg{)}^{3}$
(7)

a correction term for the fact that matter need only be lifted to the
planetary Roche lobe to escape (Erkaev et al., 2007). There exist many
variations of Eq. 6, mainly due to the fact that the planet absorbs high-
energy radiation over an area $\pi R_{XUV}^{2}>\pi R_{p}^{2}$ that should then
be averaged over the full planet surface area, which some authors include
explicitly in the energy-limited formula. We choose to instead absorb these
effects into $\eta$. In principle any planet undergoing atmospheric mass loss
will satisfy Eq. 6 for some value of $\eta$, but an ‘energy-limited’ outflow
usually refers to the case where $\eta$ is of order unity.

The energy-limited prescription finds its use mainly in population studies of
escape, when mass-loss rates need to be estimated for a large sample of
planets and performing individual hydrodynamical simulations is too
computationally expensive (e.g. Lammer et al., 2009; Owen & Wu, 2017; Rogers
et al., 2021). A main difficulty for such analyses lies in finding the value
of $\eta$ and how it depends on planetary parameters. Many works have
addressed this problem through different avenues (Owen & Jackson, 2012;
Shematovich et al., 2014; Owen & Alvarez, 2016; Caldiroli et al., 2021b).
Although this manuscript is not a systematic investigation into the value of
$\eta$, we still find it illustrative to plug the mass-loss rates we
constrained for our sample of 6 planets into Eq. 6 and compare the
corresponding heating efficiencies to the literature. This comparison is
especially useful because our mass-loss rates are constrained from
observations, while most literature studies use theoretical predictions to
evaluate $\eta$.

Figure 6: The black scatter points show our constrained mass-loss rates
corrected by $K$ as a function of the ratio of high-energy irradiation to bulk
density. The error bars indicate the results for the simulations at 1/3x and
3x EUV flux. The orange and blue scatter points show the results for the
hydrodynamic simulations of Salz et al. (2016) and Caldiroli et al. (2021b),
respectively. The yellow and purple shaded regions indicate the constraints on
the efficiency parameter of Owen & Jackson (2012) and Shematovich et al.
(2014), respectively. The green line and shaded region mark the best-fit and
1$\sigma$ contour of the efficiency constraint from Vissapragada et al.
(2022b).

Inspired by Fig. 4 of Caldiroli et al. (2021b), in Fig. 6 we show $\dot{M}/K$
versus $F_{XUV}/\rho_{p}$ of our sample together with various estimates from
the literature, where the choice of axes ensures $\eta$ is read off easily.
Including the EUV flux uncertainty, our sample spans efficiencies of roughly
2%-30%. This range generally agrees with the results of Shematovich et al.
(2014), who calculate a heating efficiency of 10%-20% for hydrogen-dominated
upper atmospheres similar to HD 209458 b. Owen & Jackson (2012) explored X-ray
driven outflows of different planetary masses and radii and found efficiencies
in the range 1%-10% for gas giant planets. We marked this region of
$\eta$-space completely in Fig. 6, but note that in principle their
efficiencies depend on the planetary parameters. Furthermore, their definition
of the efficiency is different from Eq. 6 up to a factor of a few. Despite
these considerations, the reported efficiencies seem generally consistent with
our results. However, we do find somewhat lower efficiencies than the
hydrodynamic simulations of Salz et al. (2016) and Caldiroli et al. (2021b),
who (excluding stable atmospheres) find efficiencies in the range 30%-90% for
planets with $F_{XUV}/\rho_{p}<10^{4}$ erg cm s-1 g-1. We also find lower
efficiencies than Vissapragada et al. (2022b), who fit the efficiency
parameter for a sample of five planets whose mass-loss rates they constrained
using Parker wind models. They report an efficiency of $0.41^{+0.16}_{-0.13}$,
but we note that they divided their mass-loss rates by a factor 4 to account
for flux averaging over the planet surface, so that their efficiency
definition implicitly differs from ours. Although the discrepancies with the
latter three works are interesting and potentially worthwhile exploring in
more detail, our sample size is currently not large enough to draw definitive
conclusions about the value of $\eta$.

##  6 Summary

Understanding the atmospheric mass loss of an exoplanet is crucial for
understanding the planet’s evolution. Constraining mass-loss rates is commonly
done by fitting spectroscopic observations with a set of isothermal Parker
wind models, which parametrize the outflow in the mass-loss rate and
temperature. This approach often results in a strong degeneracy between these
parameters, particularly so for spectrally unresolved observations. We aim to
resolve this degeneracy by placing prior constraints on the temperature.

We do this by simulating Parker wind models with the photoionization code
Cloudy. We develop an algorithm around the code to include hydrodynamical
effects on the temperature structure that are normally not supported. The
algorithm works by running Cloudy with a fixed temperature profile and
combining the resulting radiative heating rates with advection and expansion
terms to propose a new temperature structure based on the heating to cooling
ratio. After a few iterations, this yields a converged temperature structure
with the hydrodynamic terms included. Our code is relatively fast to execute,
and by virtue of using Cloudy includes a highly detailed treatment of the
stellar SED and NLTE processes in the upper atmosphere.

We find different thermal structures for planets with relatively high
gravitational potential like HD 209458 b and lower gravity planets like
WASP-69 b. The atmospheres of high-gravity planets ($|\phi|\gtrsim 10^{12.82}$
erg g-1) are characterized by a transition from radiative cooling to expansion
work cooling at higher altitudes. Low-gravity giant planets cool by expansion
work even at low altitudes, making the temperature at the base of the
atmosphere dependent on the assumed outflow velocity.

We run grids of Parker wind profiles for HD 209458 b, WASP-69 b, WASP-52 b,
HAT-P-18 b, HAT-P-11 b and WASP-107 b, which all have metastable helium
observations in the literature. For each Parker wind profile, we calculate the
difference between the mean temperature in the He 10830 Å line and the
isothermal value assumed to create the model. Treating this difference as a
quantifier for the self-consistency of each model, we identify a self-
consistent parameter space for each planet. We then calculate the transmission
spectrum in the He 10830 Å line for every model and fit these to the observed
values. Combining the self-consistent parameter spaces and the best-fit Parker
wind profiles of the observations give much better constraints on both the
temperature and mass-loss rate. In this way, for HD 209458 b, WASP-69 b,
WASP-52 b, HAT-P-11 b, HAT-P-18 b and WASP-107 b, we find mass-loss rates of
$\dot{M}=10^{9.84^{+0.24}_{-0.27}}$ g s-1,
$\dot{M}=10^{9.91^{+0.19}_{-0.12}}$, $\dot{M}=10^{11.06^{+0.26}_{-0.25}}$,
$\dot{M}=10^{9.70^{+0.29}_{-0.15}}$, $\dot{M}=10^{10.33^{+0.29}_{-0.24}}$ and
$\dot{M}=10^{9.96^{+0.23}_{-0.19}}$ g s-1, respectively. In an energy-limited
framework, these mass-loss rates would indicate heating efficiencies of
2%-30%.

###### Acknowledgements.

We thank the anonymous referee for providing valuable comments and suggestions
that improved the manuscript. We are grateful to G. J. Ferland and other
developers for making Cloudy publicly available, as well as providing support
on the help forum777https://cloudyastrophysics.groups.io. We appreciate the
helpful discussions with J. Kirk, C. Dominik, F. Nail and P. Uttley. We thank
SURFsara888www.surfsara.nl for the support in using the Lisa Compute Cluster.
AO gratefully acknowledges support from the Dutch Research Council NWO Veni
grant.

## References

  * Allart et al. (2019) Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58
  * Allart et al. (2018) Allart, R., Bourrier, V., Lovis, C., et al. 2018, Sci, 362, 1384
  * Alonso-Floriano et al. (2019) Alonso-Floriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019, A&A, 629, A110
  * Ankit (2020) Ankit, R. 2020, WebPlotDigizer
  * Caldiroli et al. (2021a) Caldiroli, A., Haardt, F., Gallo, E., et al. 2021a, A&A, 655, A30
  * Caldiroli et al. (2021b) Caldiroli, A., Haardt, F., Gallo, E., et al. 2021b, arXiv:2112.00744 (v1)
  * Cubillos et al. (2020) Cubillos, P. E., Fossati, L., Koskinen, T., et al. 2020, AJ, 159, 111
  * dos Santos et al. (2022) dos Santos, L. A., Vidotto, A. A., Vissapragada, S., et al. 2022, A&A, 659, A62
  * Drake et al. (2020) Drake, J. J., Kashyap, V. L., Wargelin, B. J., & Wolk, S. J. 2020, ApJ, 893, 137
  * Eastman et al. (2013) Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83
  * Erkaev et al. (2007) Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329
  * Ferland et al. (1998) Ferland, G., Korista, K., Verner, D., et al. 1998, PASP, 110, 761
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzman, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Fossati et al. (2020) Fossati, L., Shulyak, D., Sreejith, A. G., et al. 2020, A&A, 643, A131
  * Fossati et al. (2021) Fossati, L., Young, M. E., Shulyak, D., et al. 2021, A&A, 653, A52
  * France et al. (2016) France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Ginzburg et al. (2018) Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759
  * Kirk et al. (2020) Kirk, J., Alam, M. K., López-Morales, M., & Zeng, L. 2020, AJ, 159, 115
  * Kirk et al. (2022) Kirk, J., Dos Santos, L. A., López-Morales, M., et al. 2022, AJ, 164, 24
  * Koskinen et al. (2010) Koskinen, T. T., Yelle, R. V., Lavvas, P., & Lewis, N. K. 2010, ApJ, 723, 116
  * Kruijssen et al. (2020) Kruijssen, J. M. D., Longmore, S. N., & Chevance, M. 2020, ApJ, 905, L18
  * Lamers & Cassinelli (1999) Lamers, H. J. G. L. M. & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge University Press)
  * Lammer et al. (2009) Lammer, H., Odert, P., Leitzinger, M., et al. 2009, A&A, 506, 399
  * Lampón et al. (2020) Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13
  * Lampón et al. (2021) Lampón, M., López-Puertas, M., Sanz-Forcada, J., et al. 2021, A&A, 647, A129
  * Lavie et al. (2017) Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7
  * Lecavelier des Etangs et al. (2010) Lecavelier des Etangs, A., Ehrenreich, D., Vidal-Madjar, A., et al. 2010, A&A, 514, A72
  * Lee et al. (2022) Lee, E. J., Karalis, A., & Thorngren, D. P. 2022, arXiv: 2201.09898
  * Linsky et al. (2010) Linsky, J. L., Yang, H., France, K., et al. 2010, ApJ, 717, 1291
  * Lopez & Rice (2018) Lopez, E. D. & Rice, K. 2018, MNRAS, 479, 5303
  * Mansfield et al. (2018) Mansfield, M., Bean, J. L., Oklopčić, A., et al. 2018, ApJ, 868, L34
  * Matsakos & Königl (2016) Matsakos, T. & Königl, A. 2016, ApJ, 820, L8
  * Mignone et al. (2012) Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23
  * Nortmann et al. (2018) Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388
  * Oklopčić (2019) Oklopčić, A. 2019, ApJ, 881, 133
  * Oklopčić & Hirata (2018) Oklopčić, A. & Hirata, C. M. 2018, ApJ, 855, L11
  * Owen (2019) Owen, J. E. 2019, Annu. Rev. Earth Planet. Sci., 47, 67
  * 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 & Lai (2018) Owen, J. E. & Lai, D. 2018, MNRAS, 479, 5012
  * Owen et al. (2021) Owen, J. E., Murray-Clay, R. A., Schreyer, E., et al. 2021, arXiv: 2111.06094
  * Owen & Wu (2013) Owen, J. E. & Wu, Y. 2013, ApJ, 775, 105
  * Owen & Wu (2017) Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29
  * P. Loyd et al. (2016) P. Loyd, R. O., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102
  * Palle et al. (2020) Palle, E., Nortmann, L., Casasayas-Barris, N., et al. 2020, A&A, 638, A61
  * Paragas et al. (2021) Paragas, K., Vissapragada, S., Knutson, H. A., et al. 2021, ApJ, 909, L10
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664
  * Rogers et al. (2021) Rogers, J. G., Muñoz, C. J., Owen, J. E., & Booth, R. A. 2021, arXiv: 2110.15162
  * Salz et al. (2015) Salz, M., Banerjee, R., Mignone, A., et al. 2015, A&A, 576, A21
  * Salz et al. (2018) Salz, M., Czesla, S., Schneider, P. C., et al. 2018, A&A, 620, A97
  * Salz et al. (2016) Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75
  * Shaikhislamov et al. (2020) Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2020, MNRAS, 491, 3435
  * Shematovich et al. (2014) Shematovich, V. I., Ionov, D. E., & Lammer, H. 2014, A&A, 571, A94
  * Sing et al. (2019) Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91
  * Spake et al. (2021) Spake, J. J., Oklopčić, A., & Hillenbrand, L. A. 2021, AJ, 162, 284
  * Spake et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nat, 557, 68
  * Turner et al. (2020) Turner, J. D., de Mooij, E. J. W., Jayawardhana, R., et al. 2020, ApJ, 888, L13
  * Vidal-Madjar et al. (2003) Vidal-Madjar, A., des Etangs, A. L., Désert, J.-M., et al. 2003, Nat, 422, 143
  * Vissapragada et al. (2022a) Vissapragada, S., Knutson, H. A., dos Santos, L. A., Wang, L., & Dai, F. 2022a, ApJ, 927, 96
  * Vissapragada et al. (2022b) Vissapragada, S., Knutson, H. A., Greklek-McKeon, M., et al. 2022b, arXiv:2204.11865
  * Vissapragada et al. (2020) Vissapragada, S., Knutson, H. A., Jovanovic, N., et al. 2020, AJ, 159, 278
  * Wang & Dai (2021a) Wang, L. & Dai, F. 2021a, ApJ, 914, 98
  * Wang & Dai (2021b) Wang, L. & Dai, F. 2021b, ApJ, 914, 99
  * Yan & Henning (2018) Yan, F. & Henning, T. 2018, Nat Astron, 2, 714
  * Yan et al. (2021) Yan, F., Wyttenbach, A., Casasayas-Barris, N., et al. 2021, A&A, 645, A22
  * Young et al. (2020) Young, M. E., Fossati, L., Koskinen, T. T., et al. 2020, A&A, 641, A47
  * Youngblood et al. (2019) Youngblood, A., France, K., Ksokinen, T., et al. 2019, BAAS, 51, 320
  * Youngblood et al. (2016) Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101
  * Zhang et al. (2022) Zhang, M., Knutson, H. A., Wang, L., Dai, F., & Barragán, O. 2022, AJ, 163, 67

##  Appendix A Using the fixed temperature to fit observations

In Sec. 4.3, we fit a grid of Parker wind models to observations in the He
10830 Å line, resulting in the constraints shown in blue in Fig. 5. We used
the temperature structures and metastable helium densities as given by the
converged Cloudy simulations to calculate the synthetic spectral lines of the
models. Here, we investigate how these constraints change if we fit Parker
wind models with their assumed fixed temperature $T_{0}$. This has been the
standard approach in literature, and for the observational data we use in this
work, such an analysis has been performed by Lampón et al. 2020 for HD 209458
b, Vissapragada et al. 2020 for WASP-69 b, Mansfield et al. 2018 for HAT-P-11
b and Paragas et al. 2021 for HAT-P-18 b.

We used Cloudy to simulate a grid of Parker wind models at a constant
temperature $T_{0}$ for each planet. Fig. 7 shows the constraints we find for
these simulations. For HD 209458 b and WASP-52 b (the planets with the deepest
gravitational potential wells of our sample) the results are very similar to
those of Fig. 5, but for the other planets we find a stronger degeneracy
between the mass-loss rate and temperature. This stems from the fact that the
synthetic helium lines are different for simulations with fixed or converged
temperature structures, which result in different abundances of helium atoms
in the metastable state. We observe that this is particularly the case for the
lower gravity planets, likely because they have a qualitatively different
energy balance than higher gravity planets (see Sec. 4.1) in the region where
the helium line forms (typically $r\lesssim 2R_{p}$).

Figure 7: Similar to Fig. 5, but using the assumed fixed temperature $T_{0}$
instead of Cloudy’s temperature $T(r)$ to calculate the synthetic spectral
lines. As a consequence, the observational fits shown in blue are different.
The thick grey lines indicate the best-fit Parker wind models from similar
analyses in the literature: Lampón et al. 2020 for HD 209458 b, Vissapragada
et al. 2020 for WASP-69 b, Mansfield et al. 2018 for HAT-P-11 b and Paragas et
al. 2021 for HAT-P-18 b (who used a transit depth of $0.46\pm 0.12$ % that is
lower than our used value of $0.70\pm 0.16$ %). The thickness of the grey
lines is arbitrary and does not represent the confidence intervals of those
studies.

Comparing the observational constraints of Fig. 7 to Fig. 5 reveals that
simulating the temperature structure of Parker wind models with Cloudy instead
of using the assumed temperature already results in better constrained mass-
loss rates for lower gravity planets, even when no priors are placed on the
self-consistent parameter space. Further comparing our observational
constraints of WASP-69 b and HAT-P-11 b of Fig. 7 to the corresponding results
in literature (indicated by thick grey lines), we typically find a very
similar shape and offsets on the order of $\lesssim$0.2 dex in mass-loss rate.
We expect these differences to be mostly due to the use of Cloudy. For
HAT-P-18 b, the offset is a bit larger, but also partly stems from the fact
that we used the updated transit depth from Vissapragada et al. (2022b) that
is higher than the value used in Paragas et al. (2021). In the case of HD
209458 b, the difference with the constraints from Lampón et al. (2020) is on
the order or $\sim$0.5 dex, which might be due to fitting the EW instead of
the resolved line shape.