© D. C. Linssen et al. 2022

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-α line (e.g. Vidal-Madjar et al. 2003; Lecavelier des Etangs et al. 2010; Lavie et al. 2017), the Balmer-α line (e.g. Yan & Henning 2018; Yan et al. 2021) the helium line at 10 830 Å (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 10 830 Å 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.¹ 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α observations, but these analyses may be complicated by the fact that Lyα 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 Cloudy² (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 10 830 Å, 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 Sect. 2 we review the basic assumptions underlying isothermal Parker wind models and how we can test them to assess the model self-consistency. In Sect. 3, we present an algorithm to model the temperature structure of upper exoplanet atmospheres with Cloudy. In Sect. 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 Sect. 5, we compare our results and methods to similar studies and discuss the advantages and potential drawbacks of our modeling approach. We summarize in Sect. 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{\text{d}v}{\text{d}r}+\frac{1}{\rho }\frac{\text{d}P}{\text{d}r}+\frac{G{M}_{\text{p}}}{r^2}=0,$$(2)

$$-\rho v\frac{\text{d}}{\text{d}r}\left[\frac{kT}{\left(\gamma -1\right)\mu }\right]+\frac{kTv}{\mu }\frac{\text{d}\rho }{\text{d}r}+\text{Γ}+\text{Λ}=0\text{.}$$(3)

In these equations, $\dot{M}$ is the mass-loss rate, r the distance from the center of the planet, ρ the density, v the outflow velocity, P the gas pressure, G the graviational constant, M_p the planet mass, k the Boltzmann constant, T the temperature, γ the adiabatic index (5/3 for a perfect gas), µ the mean molecular weight and Γ and Λ 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, p(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 p, 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-α emission are the main contributors to Γ and Λ, respectively. Recently, Caldiroli et al. (2021) 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/µ 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/µ as r → ∞ 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 Oklopcic & Hirata (2018). Even though the constant T(r)/µ(r) does not immediately imply an isothermal structure, a fixed temperature T₀ can be assigned if a fixed value for µ(r) is calculated. This $\overline{\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₀ 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₀ 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₀. 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 ρ and v of the Parker wind profile. The comparison between T(r) and T₀ then acts as a ‘sanity check’ on the Parker wind assumption that Eq. (3) is well approximated by the parametrized T0. By quantifying the discrepancy between T(r) and T₀, 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/µ to the T(r)/µ(r) structure we get from solving Eq. (3). In literature however, Parker wind fits to observations are usually presented in terms of T₀ rather than T/µ, and thus benefit more from constraints on T₀. Furthermore, for a few test models we found that the relative differences were similar for T and T/µ.

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 Sect. 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 magnetohydrodynamics 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α and Hβ 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 ≲ 10⁻⁴ 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 (Γ and Λ 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.3log₁₀(H/C)), where H and 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.3log₁₀(H/C)). We then pass the new temperature structure on to Cloudy to get Γ and A 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 H/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 $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$1.1$}\right.<H/C<1.1$, which for a typical temperature structure means ≲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 two and ten 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 H/C ratio as described above, we construct the new temperature structure by solving Eq. (3) for T, using the µ, Γ and Λ of the previous iteration. Since Γ and Λ depend on the temperature, this approach only works well when Γ and Λ 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 ~2 iterations to converge. In radiation-dominated regimes, it is exactly the dependence of Γ and Λ on T that determines the temperature, and using Γ and Λ of the last iteration to construct the new temperature structure with Eq. (3) does not work, which is why we use the H/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 Sect. 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 WebPlotDigitizer³ 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. 1⁴.

Our results agree quite well, but there is a constant ~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.

Fig. 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.

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.35 R_J and a mass of M_p = 0.71 M_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₀ < 12 000 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. Figure 2 shows the result for the Parker wind profile with T₀ = 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.

This behavior changes for planets with lower gravitational potential (|ϕ| = 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.26 M_J, a radius of R_p = 1.057 R_J and used the version 2.2 SED of the K6 star HD 85512 from the MUSCLES survey (France et al. 2016; Youngblood et al. 2016; 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₀ < 11 000 K. Figure 3 shows a typical temperature structure, in this case for T₀ = 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)).

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 |ϕ| > 10^12.82 erg g⁻¹ 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. (2021).

Fig. 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 T0 = 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. Bottom: relevant heating and cooling rates shown in red and blue, respectively. We observe that at lower altitudes (≲2Rp), radiative heating and cooling dominate, while at higher altitudes (≳3Rp), advection heating and expansion cooling are the dominant thermal processes.

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 10 830 Å 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 10 830 Å. 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 10 830 Å (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_{\text{He}}=\frac{{\displaystyle \sum T\left(r\right)·N_{\text{He*}}\left(r\right)}}{{\displaystyle \sum N_{\text{He*}}\left(r\right)}}.$$(4)

The difference between the temperature of the He 10 830 Å line-forming region T_He and the temperature assumed to create the Parker wind profile T₀ 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 ≈ T₀. 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{{\displaystyle \sum {\left(T\left(r\right)-T_{\text{He}}\right)}^2·N_{\text{He*}}\left(r\right)}}{{\displaystyle \sum N_{\text{He*}}\left(r\right)}}.$$(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₀.

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

4.3 Mass-loss rate constraints

For each of our simulated Parker wind profiles for HD 209458 b, we calculated T_He − T₀ and σ_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 ≈ T₀, at an almost constant temperature of T ~ 8000 K for the range of mass-loss rates explored. Our code always points us towards this temperature: models with T₀ > 8000 K result in helium line formation temperatures lower than T₀, 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₀ < 8000 K reveals that Cloudy prefers hotter models.

We next fit the Parker wind models to the He 10 830 Å 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 Sect. 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 ± 0.06% in a 0.915 Å bandpass could be translated to an equivalent width of EW = 4.03 ± 0.55 mÅ.

For every simulated Parker wind profile, we calculate the transmission spectrum following the procedure described in Sect. 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 EXOFAST⁵ (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₀) to a p-value assuming it is normally distributed with standard deviation σ_T. For example, a temperature difference of T_He − T₀ = 2σ_T would be considered a 2σ discrepancy and thus translate to a prior p-value of ~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 3× and 1/3× flux at λ < 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σ contours of these two posteriors are shown in gray in the bottom right panel of Fig. 4. In total then, we find a temperature of $T={8200}_{-1100}^{+1200}$ K and a mass-loss rate of $\dot{M}={10}^{{9.84}_{-0.27}^{+0.24}}$ g s−1, where we quote the fiducial simulations for our central value and the 1σ 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 ± 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 ± 0.25% in a 2.44 Å bandpass, indicating an EW of 40.3 ± 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 ± 15 ppm and a 3560 ± 34 ppm transit depth in the 10 809–10 858 Å channel. Adding the uncertainties in quadrature, we obtain EW = 9.26 ± 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 ± 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 ± 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; 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 є 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.

Fig. 4 Top left: helium line-forming temperatures THe (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 gray dots). The gray 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 THe (top left panel) and the assumed isothermal value T0. The dashed and dotted lines indicate the 1σ and 2σ standard deviations of the self-consistent region based on Eq. (5), respectively. Top right: compatibility between the He 10 830 Å 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: joint constraint on the mass-loss rate and temperature. The orange line and dark and light areas are the self-consistent models, 1σ and 2σ 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σ and 2σ contours from the top right panel, respectively, which are our likelihood values. The 1σ contour of the posterior is shown in black. The 1σ posteriors of the 3× and 1/3× EUV flux simulations are shown in gray.

Fig. 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.27}^{+0.24}}$ g s−1 and a temperature of $T={8200}_{-1100}^{+1200}$ K. This is a 2σ 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 ~ 8600 K and a maximum mass-loss rate of $\dot{M}\le {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}\ge {10}^{10}$ g s−1 from Vidal-Madjar et al. (2003), who observed HD 209458 b in the Lyα 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.15}^{+0.29}}$ g s−1 and $T={4500}_{-1100}^{+1300}$ K. This is a 2σ 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 ± 700 K. Their temperature is 2σ 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.12}^{+0.19}}$ g s−1 and $T={5300}_{-300}^{+1300}$ 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.19}^{+0.23}}$ g s−1 and $T={3400}_{-400}^{+1400}$ 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 ≳ 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 10 830 Å 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₀ 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. A.1. 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. Figures 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. 2021).

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 10 830 Å 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 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₀ 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 Sect. 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 1 R_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. (2021) show a more physical picture, where even for low-gravity planets, the velocity approaches v ~ 0 at <1.1 R_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_{\text{p}}^3{F}_{\text{XUV}}}{G{M}_{\text{p}}K},$$(6)

where η is the heating efficiency and

$$K=1-\frac{3}{2}\frac{R_{\text{p}}}{R_{\text{Roche}}}+\frac{1}{2}{\left(\frac{R_{\text{p}}}{R_{\text{Roche}}}\right)}^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_{\text{XUV}}^2>\pi R_{\text{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 η. In principle any planet undergoing atmospheric mass loss will satisfy Eq. (6) for some value of η, but an ‘energy-limited’ outflow usually refers to the case where η 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 η 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. 2022). Although this manuscript is not a systematic investigation into the value of η, we still find it illustrative to plug the mass-loss rates we constrained for our sample of six 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 η.

Inspired by Fig. 5 of Caldiroli et al. (2022), in Fig. 6 we show $K\dot{M}$ versus F_XUV/ρ_p of our sample together with various estimates from the literature, where the choice of axes ensures η 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 η-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. (2022), who (excluding stable atmospheres) find efficiencies in the range 30–90% for planets with F_XUV/ρ_p < 10⁴ erg cm s−1 g⁻¹. 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.13}^{+0.16}$, but we note that they divided their mass-loss rates by a factor of four 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 η.

Fig. 6 Mass-loss rate corrected by K as a function of the ratio of high-energy irradiation to bulk density. Our constraints are indicated by black scatter points, where the error bars indicate the results for the simulations at 1/3× and 3× EUV flux. The orange and blue scatter points show the results for the hydrodynamic simulations of Salz et al. (2016) and Caldiroli et al. (2022), 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σ contour of the efficiency constraint from Vissapragada et al. (2022b).

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 (|ϕ| ≳ 10^12.82 erg g⁻¹) 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 10 830 Å 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 10 830 Å 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.27}^{+0.24}}$ g s−1, $\dot{M}={10}^{{9.91}_{-0.12}^{+0.19}},\dot{M}={10}^{{11.06}_{-0.25}^{+0.26}},\dot{M}={10}^{{9.70}_{-0.15}^{+0.29}},\dot{M}={10}^{{10.33}_{-0.24}^{+0.29}}$ and $\dot{M}={10}^{{9.96}_{-0.19}^{+0.23}}$ 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 forum (https://cloudyastrophysics.groups.io). We appreciate the helpful discussions with J. Kirk, C. Dominik, F. Nail and P. Uttley. We thank SURFsara (www.surfsara.nl) for the support in using the Lisa Compute Cluster. AO gratefully acknowledges support from the Dutch Research Council NWO Veni grant.

Appendix A Using the fixed temperature to fit observations

In Sect. 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₀. 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₀ for each planet. Fig. A.1 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 Sect. 4.1) in the region where the helium line forms (typically r ≲ 2R_p).

Comparing the observational constraints of Fig. A.1 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. A.1 to the corresponding results in literature (indicated by thick gray lines), we typically find a very similar shape and offsets on the order of ≲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 ~0.5 dex, which might be due to fitting the EW instead of the resolved line shape.

Fig. A.1 Similar to Fig. 5, but using the assumed fixed temperature T0 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 gray 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 ± 0.12% that is lower than our used value of 0.70 ± 0.16%). The thickness of the gray lines is arbitrary and does not represent the confidence intervals of those studies.

References

All Tables

All Figures

Fig. 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.

In the text

Fig. 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 T0 = 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. Bottom: relevant heating and cooling rates shown in red and blue, respectively. We observe that at lower altitudes (≲2Rp), radiative heating and cooling dominate, while at higher altitudes (≳3Rp), advection heating and expansion cooling are the dominant thermal processes.

In the text

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

Fig. 4 Top left: helium line-forming temperatures THe (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 gray dots). The gray 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 THe (top left panel) and the assumed isothermal value T0. The dashed and dotted lines indicate the 1σ and 2σ standard deviations of the self-consistent region based on Eq. (5), respectively. Top right: compatibility between the He 10 830 Å 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: joint constraint on the mass-loss rate and temperature. The orange line and dark and light areas are the self-consistent models, 1σ and 2σ 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σ and 2σ contours from the top right panel, respectively, which are our likelihood values. The 1σ contour of the posterior is shown in black. The 1σ posteriors of the 3× and 1/3× EUV flux simulations are shown in gray.

In the text

Fig. 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.

In the text

Fig. 6 Mass-loss rate corrected by K as a function of the ratio of high-energy irradiation to bulk density. Our constraints are indicated by black scatter points, where the error bars indicate the results for the simulations at 1/3× and 3× EUV flux. The orange and blue scatter points show the results for the hydrodynamic simulations of Salz et al. (2016) and Caldiroli et al. (2022), 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σ contour of the efficiency constraint from Vissapragada et al. (2022b).

In the text

Fig. A.1 Similar to Fig. 5, but using the assumed fixed temperature T0 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 gray 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 ± 0.12% that is lower than our used value of 0.70 ± 0.16%). The thickness of the gray lines is arbitrary and does not represent the confidence intervals of those studies.

In the text