The open-source sunbather code: Modeling escaping planetary atmospheres and their transit spectra

Atmospheric escape is thought to significantly influence the evolution of exoplanets, especially sub-Jupiter planets on short orbital periods. Theoretical models predict that hydrodynamic escape could erode the atmospheres of such gaseous planets, leaving only a rocky core. Deriving atmospheric mass-loss rates from observations is necessary to check these predictions. One of the ways to obtain mass-loss-rate estimates is to fit transit spectra of the 10830 Å helium or UV metal lines with Parker wind models. We aim to provide the community with a tool that enables this type of analysis, and present sunbather, an open-source Python code that can be used to model escaping exoplanet atmospheres and their transit spectra. sunbather incorporates the Parker wind code p-winds and the photoionization code Cloudy, with the ability to calculate any currently known spectral tracer, using an arbitrary atmospheric composition. With sunbather, we investigate how the atmospheric structure of a generic hot-Neptune planet depends on metallicity. We find that the mass-loss rate drops by roughly one order of magnitude as we increase the metallicity from solar to 50 times solar. Line cooling by metal species is already important for a solar composition, and is even more so at higher metallicity. We then demonstrate how sunbather can be used to interpret observations of spectral lines that form in the upper atmosphere. We fit the observed helium spectrum of the mini-Neptune TOI-2134 b and show how, even for helium data, the inferred mass-loss rate can change by a factor of up to three, depending on the assumed metallicity.

1 Introduction

Recent advances in the field of exoplanet research have identified a lack of sub-Jovian planets, as well as a bimodal distribution of the radii of sub-Neptune-sized planets, both at short orbital periods (Szabó & Kiss, 2011; Beaugé & Nesvorný, 2013; Owen & Wu, 2013; Lundkvist et al., 2016; Fulton et al., 2017). Atmospheric escape is thought to play a major role in sculpting these features. Theoretical works aiming to explain the observed demographic characteristics in this way typically prescribe analytical mass-loss-rate estimates based on the bulk parameters of planets (e.g., Owen & Wu, 2017; Ginzburg et al., 2018). However, the relative importance of atmospheric escape compared to other proposed mechanisms such as formation and migration is still a matter of debate (e.g., Owen & Lai, 2018; Lee et al., 2022; Burn et al., 2024). The same is true for the physical details of the atmospheric escape process, which can be due to photoevaporation or core-powered mass loss, or indeed both (e.g., Owen & Schlichting, 2024), with possible contributions from giant impacts (e.g., Chance et al., 2022; Zhong et al., 2024). Therefore, obtaining mass-loss rates that are based on observational measurements for a large sample of individual planets would allow us to benchmark the theoretical predictions and test the atmospheric escape hypothesis (Vissapragada et al., 2022b).

Observational evidence of atmospheric escape for individual planets is gathered through transit spectroscopy of specific spectral lines. In the early days of atmospheric characterization, this consisted of UV observations of the hydrogen Lyman-$\alpha$ line and lines from metals such as oxygen, carbon, and magnesium (e.g., Vidal-Madjar et al., 2003, 2004; Fossati et al., 2010). In more recent years, the 10,830 Å line of metastable helium has additionally been identified as a tracer of atmospheric escape (Spake et al., 2018; Oklopčić & Hirata, 2018; Nortmann et al., 2018; Allart et al., 2018). Consequently, the number of atmospheric escape detections has rapidly grown to a few dozen (dos Santos, 2023), enabling the search for trends at the population level (e.g., Allart et al., 2023; Bennett et al., 2023; Krishnamurthy & Cowan, 2024).

Based on the specific spectral line considered, as well as the quality and complexity of the data, different types of models may be needed for their interpretation. Models ranging in complexity from a relatively simple 1D isothermal Parker wind (e.g., Mansfield et al., 2018) to expensive 3D coupled hydrodynamics-thermochemistry simulations (e.g., Wang & Dai, 2021) have been applied in the literature. With the goal being to infer a mass-loss rate from a spectral line measurement, the isothermal Parker wind prescription has been widely used because it treats the mass-loss rate as a free parameter (e.g., Mansfield et al., 2018; Lampón et al., 2020, 2021a, 2023; Kasper et al., 2020; Palle et al., 2020; Vissapragada et al., 2020, 2021, 2022b; Zhang et al., 2021, 2023a, 2023b; Krishnamurthy et al., 2021, 2023; Paragas et al., 2021; Czesla et al., 2022; dos Santos et al., 2022; Kirk et al., 2022; Spake et al., 2022; Allart et al., 2023; Bennett et al., 2023; Gaidos et al., 2023; Orell-Miquel et al., 2023; Gully-Santiago et al., 2024; Pérez-González et al., 2024). This prompted dos Santos et al. (2022) to release p-winds, an open-source implementation of the isothermal Parker wind model that includes a module to produce metastable helium transit spectra. Although the isothermal Parker wind model has thus far been mainly used for interpreting helium observations, it need not be limited to this spectral line, and the newest version of p-winds includes the ability to make carbon and oxygen spectra (dos Santos et al., 2023b). Indeed, in addition to the upsurge in the number of helium studies in this regard, there have been several recent detections of escaping metals in the UV (e.g., Sing et al., 2019; Cubillos et al., 2020, 2023; García Muñoz et al., 2021; Ben-Jaffel et al., 2022; Sreejith et al., 2023), which provide highly complementary data that would benefit from similar modeling strategies.

Here, we aim to address this need and present sunbather, an open-source code to model the structure and transit spectra of escaping exoplanet atmospheres. The code couples p-winds to the detailed nonlocal thermodynamic equilibrium (NLTE) photoionization code Cloudy (Ferland et al., 1998, 2017; Chatzikos et al., 2023). This enables the calculation of atmospheres of arbitrary composition and of any currently known spectral line tracing escape (Linssen & Oklopčić, 2023). Starting from the isothermal Parker wind as produced by p-winds, sunbather calculates a more realistic nonisothermal temperature profile, which can have important consequences for the predicted transit spectrum. As shown in Linssen et al. (2022), the nonisothermal profile can additionally be used to restrain the model parameter space and yield more stringent mass-loss-rate constraints.

This paper is organized as follows. In Sect. 2, we describe the theory behind the model, and explain how it is implemented in three main Python modules. In Sect. 3, we use sunbather for a theoretical study of the effects of high metallicity on the upper atmospheric structure and transit signatures of a sub-Neptune planet. In Sect. 4, we demonstrate how sunbather can be used to constrain mass-loss rates from spectral-line observations. In Sect. 5, we summarize our findings. A reproduction package for all of the results of this manuscript is available online.

2 The sunbather code

We model an exoplanet atmosphere that is escaping hydrodynamically. In this escape regime, gas particles obtain an average internal energy comparable to that needed to escape from the planet (Gronoff et al., 2020). Hence, thermal expansion of the atmosphere drives a bulk outflow, requiring that the gas is collisional up until the sonic point (Owen, 2019). The thermal energy required can either be supplied by stellar high-energy radiation (called photoevaporation) or the planet’s interior luminosity (called core-powered mass-loss). Hydrodynamic escape is expected to be the dominant mechanism for highly irradiated planets.

Following Murray-Clay et al. (2009), the Euler equations describing a steady-state hydrodynamic outflow in the planetary frame in 1D are as follows. Mass continuity can be written as

$$\frac{\partial}{\partial r}(r^{2}\rho v)=0,$$ (1)

where $r$ is the distance to the planet center, $\rho$ the density and $v$ the velocity. Momentum conservation implies

$$\rho v\frac{\partial v}{\partial r}+\frac{\partial P}{\partial r}+\frac{GM_{p}\rho}{r^{2}}+\frac{3GM_{*}\rho r}{a^{3}}=0,$$ (2)

where $P$ is pressure, $G$ the gravitational constant, $M_{p}$ the planet mass, $M_{*}$ the stellar mass, and $a$ the semi-major axis. The last term on the left-hand side is the stellar tidal gravity term. It can be optionally excluded in sunbather. Energy balance reads

$$-\rho v\frac{\partial}{\partial r}\Bigg{[}\frac{kT}{(\gamma-1)\mu}\Bigg{]}+\frac{kTv}{\mu}\frac{\partial\rho}{\partial r}+\Gamma+\Lambda=0,$$ (3)

where $k$ is the Boltzmann constant, $T$ the temperature, $\gamma=5/3$ the adiabatic index for a perfect gas, $\mu$ the mean particle mass, $\Gamma$ the radiative heating rate, and $\Lambda$ the radiative cooling rate. Solving this system of equations results in a “self-consistent” solution to the problem. One obtains the radial density, velocity, and temperature profiles, as well as the corresponding mass-loss rate of the planet. This exercise is performed by various open-source codes, such as TPCI (Salz et al., 2015), ATES (Caldiroli et al., 2021), and AIOLOS (Schulik & Booth, 2023). A useful asset of these self-consistent codes is that they have predictive power of the atmospheric mass-loss rate. However, the associated predicted transit spectrum may be inconsistent with observational data.

The isothermal Parker wind (Parker, 1958; Lamers & Cassinelli, 1999) is an alternative type of model that parametrizes the mass-loss rate and temperature of the atmosphere. The model thus does not aim to include all physics and solve the atmospheric structure from first principle. Rather, it assumes ignorance of some parameters and aims to infer those from data. Mathematically, the isothermal Parker wind description comes down to ignoring Eq. 3. Instead, the mass and momentum equations (Eqs. 1 and 2) are solved while assuming a constant sound speed

$$v_{s}=\sqrt{\frac{kT}{\mu}}.$$ (4)

One then obtains the radial density and velocity profiles (see e.g., Lamers & Cassinelli, 1999), and a transit spectrum can be calculated. Using this model, an observed transit spectrum can virtually always be matched with some values of the mass-loss rate and temperature, allowing these parameters to be retrieved from observations. This is especially useful for the mass-loss rate, which is the quantity we are often most interested in.

A common practice in the literature is to fit isothermal Parker wind models with different mass-loss rates and temperatures to an observed transit spectrum (e.g., Lampón et al., 2020; Vissapragada et al., 2020; Allart et al., 2023). However, the best-fit models often show a degeneracy between the temperature and mass-loss rate, making it difficult to constrain them. Especially retrievals on data with low spectral resolution suffer from this degeneracy. It is less of a problem in higher resolution spectra, since the temperature can be constrained from the spectral line width.

In low-resolution spectra, an external constraint on the temperature of the atmosphere can help to break the degeneracy and hence constrain the mass-loss rate (Vissapragada et al., 2022a; Linssen et al., 2022). In sunbather, we provide a way to find a sensible temperature of the Parker wind and hence allow constraining the mass-loss rate. After having calculated the density and velocity profiles of the isothermal Parker wind, we plug them into Eq. 3 and solve for a new temperature profile that is nonisothermal. The density, velocity, and temperature profiles are thus not entirely self-consistent, given that the former two are calculated assuming an isothermal profile. That being said, exactly this fact can be leveraged to estimate a sensible Parker wind temperature. If the isothermal and nonisothermal temperature profiles are similar, the atmospheric structure and hence temperature is relatively self-consistent. Limiting the model parameter space to only those Parker wind models (for example by using a prior) means that effectively the mass-loss rate is the only free parameter left to fit to the data. We first demonstrated this method in Linssen et al. (2022). A second reason to calculate a nonisothermal temperature profile is that it should simply be more realistic than an assumed isothermal profile. Apart from affecting the atmospheric structure and line widths, the temperature also affects the chemical state of gas (i.e.,the ionization state and atomic level populations) and hence the transit spectrum.

sunbather includes the 30 lightest elements (up until zinc) and allows for an arbitrary composition. The inclusion of metal species compared to a pure hydrogen/helium composition manifests itself in the outflow properties and transit spectrum in various ways, which we explore in Sect. 3. The composition of the atmosphere is kept constant with radius in the atmosphere. Heavy metal particles are assumed to be dragged along efficiently in the outflow so that there is no mass fractionation. This assumption breaks down when the mass-loss rate is low and only the lighter elements are allowed to escape (Hunten et al., 1987). It is up to the user of the code to consider whether each species is expected to escape and if the chosen composition is sensible. This can for example be estimated using the crossover mass (Hunten et al., 1987). The cases in which heavier particles either escape efficiently or do not escape at all can be modeled well by choosing the appropriate (constant) composition. Only the intermediate case in which heavy particles do escape, but at a lower rate than the lighter elements, cannot be appropriately modeled as it requires abundance patterns that change with radius in the atmosphere.

The transit spectrum can be calculated over any wavelength range between 1 Å and 10 $\mathrm{\mu}$m. It aims to include all atomic and ionic spectral lines present in the NIST Atomic Spectra Database. Not all of these lines can be calculated (see Sect. 2.3 for more details), but the strongest lines are included. To quantify, 29,530 spectral lines are included, equal to 40% of the lines in the NIST database. Sources of continuum opacity, as well as molecular absorption lines cannot be calculated in the current version of sunbather, but molecules are also expected to be dissociated in the outflow. Our model is 1D and we assume spherical symmetry, and so we cannot produce asymmetric spectral line shapes or light curves. Spectral lines will always be symmetric around their rest-frame wavelength.

We now describe the setup of sunbather in more detail. This is divided into three sections that correspond to the modules of sunbather in which they are implemented. Fig. 1 presents a simplified visualization of the computational scheme.

Figure 1: Simplified visualization of the modeling setup of sunbather. The variables are as follows: $R_{p}$ and $M_{p}$ are the planet radius and mass, $a$ is the orbital distance (only circular orbits are possible), $M_{*}$ and $F_{*}(\lambda)$ are the stellar mass and spectral energy distribution, $T_{0}$ and $\dot{M}$ are the (constant) temperature and mass-loss rate of the Parker wind, $Z$ denotes the atmospheric bulk composition, $\mu(r)$, $\rho(r)$, $v(r),$ and $T(r)$ are the atmospheric mean particle mass, density, velocity, and (nonconstant) temperature profiles, $n_{X}(r)$ represents the number density of any “chemical species” (i.e., any energy level of an atom or ion), $b$ and $\phi$ are the planet transit impact parameter and orbital phase, and $R_{*}$, $u_{1}(\lambda)$ and $u_{2}(\lambda)$ are the stellar radius and quadratic limb-darkening parameters.

2.1 Calculating the density and velocity structure: The construct_parker module

In this module, we calculate isothermal Parker wind density and velocity profiles using a combination of p-winds (dos Santos et al., 2022) and Cloudy (Ferland et al., 1998, 2017; Chatzikos et al., 2023). Despite what the name suggests, the wind is not strictly isothermal in this formalism. Instead, the sound speed and thus the ratio of $T(r)/\mu(r)$ is constant. $\mu(r)$ depends on the composition of the atmosphere (which we keep constant with radius), but also the ionization fractions, which do change with radius. $\mu$ therefore typically varies by a factor $\lesssim 2$ throughout the atmosphere. If one calculates an appropriately averaged value $\bar{\mu}$, a corresponding averaged isothermal value $T_{0}$ can also be assigned (Appendix A of Lampón et al., 2020).

We use p-winds to calculate the density and velocity profiles of the atmosphere. This means it solves Eqs. 1 and 2; we refer the reader to dos Santos et al. (2022) for the details of this calculation. To run, p-winds needs the planet mass $M_{p}$, the planet radius $R_{p}$, and the stellar spectrum $F_{*}(\lambda)$ at the planet orbital distance $a$ as input. Additionally, the three free parameters of the model are the isothermal temperature $T_{0}$, the bulk mass-loss rate $\dot{M}$, and the composition of the atmosphere. These parameters cannot be predicted by the model - they have to be assumed. In p-winds, the atmosphere is assumed to be hydrogen/helium only, meaning that choosing the composition is equivalent to choosing a hydrogen abundance fraction $f_{H}$. The code calculates the ionization structures of hydrogen and helium, resulting in a $\mu(r)$ profile and associated $\bar{\mu}$. Together with $T_{0}$, $\bar{\mu}$ determines the speed of sound, which is needed for the calculation of the density and velocity profiles. If the user of sunbather chooses an atmospheric composition that is hydrogen/helium only, p-winds can be used standalone to obtain the density and velocity profiles. However, when the chosen atmospheric composition includes metals, the $\mu(r)$-structure as calculated by p-winds will not be correct. We therefore use Cloudy to calculate $\mu(r)$ when metals are included.

Cloudy is a 1D NLTE photoionization code applicable to many different astrophysical problems. Currently, versions v17.02 and v23.01 are supported by sunbather, and the user may decide which version to install (no significant differences between the two versions were found). Cloudy calculates microphysical processes inside a gas cloud irradiated by a light source, predicting the ionization and energy level populations, as well as the emerging spectrum. It includes extensive databases of energy levels, transition probabilities, and collision rates for different molecular, atomic, and ionic species. This makes it very suitable for our purpose of predicting the physical state of upper exoplanet atmospheres, which due to their low densities are expected to show a complex chemical state governed by NLTE effects. In recent years, Cloudy has been increasingly applied to the study of exoplanet atmospheres (e.g., Salz et al., 2016; Fossati et al., 2020, 2021, 2023; Turner et al., 2020; Young et al., 2020, 2024; Zhang et al., 2021, 2022; Linssen et al., 2022; Linssen & Oklopčić, 2023; Kubyshkina et al., 2024).

The isothermal Parker wind density and velocity profiles are solved in an iterative manner when metals are included (visualized in the left part of Fig. 1). We start by running p-winds to calculate the density and velocity profile with hydrogen and helium only. We then simulate this density profile with Cloudy, while fixing an isothermal temperature profile $T(r)=T_{0}$, and using the full composition that includes metals. Cloudy calculates - amongst other things - the ionization profiles of all elements and the corresponding free electron number density profile. If we ignore molecules, we only need the electron density to calculate the mean particle mass as

$$\mu(r)=\frac{\sum_{s}X_{s}m_{s}}{1+\frac{n_{e}(r)}{\rho(r)}\sum_{s}X_{s}m_{s}},$$ (5)

where $X_{s}$ and $m_{s}$ are the abundance and mass of species $s$, respectively, and $n_{e}$ is the free electron number density. The abundances are expressed as a dimensionless fraction and obey $\sum_{s}X_{s}=1$. We use Eq. A.3 from Lampón et al. (2020) to find $\bar{\mu}$ from $\mu(r)$. This completes the first iteration and we start the second. We run p-winds again and use the new $\bar{\mu}$ value, while bypassing its own calculation of $\bar{\mu}$. This results in a new density and velocity profile. We simulate the density profile with Cloudy to obtain a new $\bar{\mu}$, completing the second iteration. We repeat these steps until $\bar{\mu}$ is converged to $(|\bar{\mu}_{i}-\bar{\mu}_{i-1}|)/\bar{\mu}_{i-1}<0.01$, where $i$ is the iteration number. sunbather also stops if convergence is not reached after seven iterations (the user can change the convergence threshold and maximum iteration number).

The runtime of this module depends on the composition. For a pure hydrogen/helium composition, p-winds can calculate the density and velocity profiles in a second. When the p-winds/Cloudy iterative scheme is applied, the runtime is of order ten minutes. One could also still invoke the iterative scheme when calculating a pure hydrogen/helium composition, as Cloudy may report a different $\mu(r)$-structure from p-winds.

2.2 Calculating the temperature structure: The convergeT_parker module

In this module, we calculate a nonisothermal temperature structure using Cloudy. This means we now introduce Eq. 3 and solve it for $T(r)$ while using $\rho(r)$ and $v(r)$ that were calculated before under the isothermal assumption. The four terms of Eq. 3 represent advection, expansion, radiative heating, and radiative cooling, respectively. The advection and expansion terms depend explicitly on the temperature profile. The radiative heating and cooling terms take a few minutes to calculate using Cloudy and implicitly depend on the temperature in a nontrivial manner. Therefore, we cannot instantly solve Eq. 3 for the temperature profile, but use an iterative algorithm instead. The middle part of Fig. 1 shows a simplified schematic of the iterative algorithm, and Fig. 2 shows an example case.

We start by running the density profile through Cloudy with an initial temperature profile. The user of sunbather can choose whether this initial temperature profile is

1. the isothermal profile assumed when calculating the density and velocity profiles,

2. the temperature profile of a previously calculated model with a similar $T_{0}$ and $\dot{M}$,

3. unspecified, in which case Cloudy will calculate it assuming radiative equilibrium.

The chosen initial temperature profile does not affect the final nonisothermal profile, but it does affect the number of iterations needed to find the solution. The output of this first Cloudy simulation is $\Gamma(r)$, $\Lambda(r)$, $\mu(r)$ (after applying Eq. 5), as well as $T(r)$ if option iii) was chosen. For this initial temperature profile, Eq. 3 will not yet be satisfied. Therefore, we guess a new temperature structure and run Cloudy again to obtain new $\Gamma(r)$, $\Lambda(r)$, and $\mu(r)$ profiles. This is repeated until Eq. 3 is respected. We use two complementary algorithms to guess the temperature structure of the next iteration, which we describe in the following two subsections.

2.2.1 Algorithm 1: “Relaxing” the temperature profile

In this algorithm, we guess a new temperature profile based on the heating/cooling inequality of the last iteration. We aim to find the temperature profile for which the total heating and cooling rate are equal and hence their ratio equal to 1, as that implies we have solved Eq. 3. After running Cloudy with the temperature profile of the last iteration, we have all the quantities needed to calculate the four different terms of Eq. 3. The advection term can be positive (heating) or negative (cooling), depending on the slope of $T(r)/\mu(r)$, while expansion and $\Lambda(r)$ are always negative and $\Gamma(r)$ is always positive. Fig. 2b shows an example of these rates for the temperature profile of the second iteration. We sum the different heating and cooling rates into a total heating rate $\mathcal{H}(r)$ and total cooling rate $\mathcal{C}(r)$ (which we define as positive here) as a function of radius. Fig. 2c shows the ratio $\mathcal{H}(r)/\mathcal{C}(r)$ as a red line if heating is stronger or $\mathcal{C}(r)/\mathcal{H}(r)$ as a blue line if cooling is stronger. Based on this ratio, we guess a new temperature structure, assigning a higher temperature to radii where $\mathcal{H}>\mathcal{C}$, and a lower temperature to radii where $\mathcal{C}>\mathcal{H}$. Specifically,

$$T_{i+1}(r)=\begin{cases}T_{i}(r)\dfrac{1}{1-\delta_{i}(r)\cdot\mathrm{log}_{10}\left(\frac{\mathcal{H}(r)}{\mathcal{C}(r)}\right)}&\text{if }\mathcal{H}(r)>\mathcal{C}(r)\\ \\ T_{i}(r)\left[1-\delta_{i}(r)\cdot\mathrm{log}_{10}\left(\frac{\mathcal{C}(r)}{\mathcal{H}(r)}\right)\right]&\text{if }\mathcal{H}(r)<\mathcal{C}(r),\end{cases}$$ (6)

where $i$ is the iteration number and $\delta(r)$ is a scaling factor that sets how large the temperature changes are between iterations. In the first iteration, $\delta_{1}(r)=0.3$ at all radii. In some cases, the changes to the temperature structure are too aggressive and $T(r)$ will fluctuate up and down between iterations. We check for radial bins in which such temperature fluctuations happen, and we lower the scaling factor there as $\delta_{i+1}(r)=2/3\cdot\delta_{i}(r)$. This ensures smaller changes in temperature between iterations, preventing fluctuations and allowing the temperature structure to converge. $T_{i+1}/T_{i}$ is furthermore bounded by 0.5 and 2 to prevent extremely large jumps in temperature. In Fig. 2a, the dotted line labeled “iteration 3” shows the newly guessed temperature profile after applying the algorithm described in this section. We do not immediately start the third iteration with this profile, however. First, we apply a second, complementary algorithm described in the next subsection to refine the guessed temperature profile of iteration 3 further.

2.2.2 Algorithm 2: “Constructing” the temperature profile

In this algorithm, we further modify the guessed temperature from the relaxation algorithm described in the previous section, before starting the next iteration. Contrary to the relaxation algorithm, we here do not use the temperature profile of the last iteration. Instead, we use only $\mu(r)$, $\Gamma(r)$ and $\Lambda(r)$ of the last iteration as given by Cloudy. Together with $\rho(r)$ and $v(r)$, we can solve Eq. 3 for $T(r)$, which is the only remaining variable. The reason that this does not immediately yield the correct temperature structure, is that $\Gamma$ and $\Lambda$ implicitly depend on the temperature, meaning that the $\Gamma(r,T)$ and $\Lambda(r,T)$ values from the last iteration are not consistent anymore with the newly constructed $T_{i+1}(r)$. Iterations are therefore still required to find the converged temperature structure.

This algorithm will not work well in parts of the atmosphere where the dominant heating and cooling terms are the radiation rates. This is because we solve Eq. 3 while keeping the radiation rates fixed at those from the last iteration, such that only the expansion and advection terms can change. If the latter are not the dominant terms, the algorithm will make huge adjustments to the temperature structure in order to change the expansion and advection terms enough to achieve heating/cooling balance. However, smaller changes to the temperature profile would actually already lead to energy balance, since $\Gamma(r,T)$ and $\Lambda(r,T)$ would change accordingly as well. Oppositely, this algorithm works well and achieves convergence quickly when radiation is relatively unimportant. In this case, the radiation rates are still fixed at those from the last iteration and will again not be fully consistent with the newly constructed temperature structure. However, this is not important for the total heating/cooling balance, which is dominated by advection and/or expansion. For these reasons, we always employ the relaxation algorithm of Sect. 2.2.1 first, and afterwards only employ the construction algorithm of this section over the region of the atmosphere where radiation is relatively unimportant (provided such a region exists).

We construct the new temperature profile as follows. Given a starting value $T(r_{start})$ specified at one particular radial bin $r_{start}$, we calculate $T$ in the next radial bin by minimizing the left-hand side of Eq. 3 with respect to $T$ (using the minimize_scalar function from the scipy library). We can then calculate $T$ in the following radial bin and proceed in this way until the full profile is found. The temperature profile constructed with this method is very sensitive to $T(r_{start})$, and therefore we only apply the algorithm if there is a starting point where the temperature profile is already almost converged, through the use of the relaxation algorithm described in Sect. 2.2.1. Additionally, we only apply the algorithm if the whole region of the atmosphere that we construct (i.e., $r>r_{start}$) is dominated by expansion or advection, for the reasons discussed above. The exact criterion to find a suitable $r_{start}$ is somewhat complex and we refer the interested reader to the sunbather source code. Roughly speaking, we require $\mathcal{H}(r_{start})/\mathcal{C}(r_{start})$ or $\mathcal{C}(r_{start})/\mathcal{H}(r_{start})$ (whichever is higher) to be smaller than 1.3, but this value is allowed to be higher if the contributions from radiative heating and cooling are sufficiently low.

2.2.3 Reaching convergence

Our complete method to calculate the nonisothermal temperature profile consists of iteratively running the temperature profile through Cloudy, updating the temperature profile using the described algorithms, and radially smoothing the temperature profile over small scales. We smooth the temperature profile mainly because the advection rate shows small discrete jumps that propagate to jumps in the newly guessed temperature profile. These jumps arise from the dependence of the advection rate on the temperature gradient and Cloudy’s limited precision in reporting the temperature. Without smoothing, small discrete jumps in the temperature profile would lead to much larger changes in the advection rate in the next iteration, quickly destabilizing the algorithm. We consider the temperature profile converged when $1/f_{c}<\mathcal{H}(r)/\mathcal{C}(r)<f_{c}$ at all radii. The default value of the convergence threshold $f_{c}$ is 1.1 and can be changed by the user of the code. In some cases, the smoothing of the temperature profile prevents it from reaching a convergence factor of $f_{c}$, even though the profile does not change substantially anymore between iterations. These solutions are thus correct up to the level of the smoothing (typically $\lesssim 50$ K), and we consider them valid as well. The runtime of this module is on the order of ten minutes.

We finally clarify that the advection term of Eq. 3 only accounts for the advective heat flow. In principle, the advective particle flow should also lead to the transport of atoms from the deeper layers of the atmosphere to the ionic upper layers of the atmosphere, thereby influencing the degree of ionization. This would in turn affect the mean particle mass that we calculate in Eq. 5, as well as the radiative heating and cooling rates through different ionization and recombination rates. Properly including the advective velocity component in Cloudy would require editing the source code (Salz et al., 2015) and drastically increase its computation time, and so we ignore these second-order effects on $\mu$, $\Gamma$, and $\Lambda$.

Figure 2: Iterative method to solve for the nonisothermal temperature profile. a) Temperature profiles of the different iterations $i$. Convergence is reached after nine iterations. Upon each iteration, two algorithms are applied to guess the temperature profile. Dotted lines show the profile based on the relax algorithm of Sect. 2.2.1 and are labeled (r). Solid lines then show the profile after the construct algorithm of Sect. 2.2.2 has also been applied and are labeled (c). b) The different heating and cooling rates corresponding to the temperature profile of the second iteration. c) Ratio of the total heating rate to total cooling rate, for the rates shown in panel (b). If the ratio is $<f_{c}=1.1$ everywhere, the temperature profile is considered converged. Based on the $\mathcal{H}/\mathcal{C}$ and $\mathcal{C}/\mathcal{H}$ ratios shown here, the guessed $i$=3 (r) profile will be cooler at $r<5R_{p}$ and slightly hotter at $r>5R_{p}$, as confirmed in panel (a).

2.3 Calculating the transmission spectrum: The RT module

In this module, we use the Cloudy output to calculate the transmission spectrum. We do not use the spectrum predicted by Cloudy. Cloudy is a hydrostatic code and thus does not include spectral line broadening due to the outflow velocity. Furthermore, instead of one simulation along the substellar direction, multiple simulations at different impact parameters from the planet would be needed to produce the correct transit geometry (see e.g., Young et al., 2020). Therefore, we instead postprocess the chemical state of the atmosphere that Cloudy calculates with a ray-tracing calculation as follows.

The user provides a wavelength region and a set of atomic and ionic species for which to calculate the transit spectrum. For each species, sunbather includes a downloaded table of all spectral line coefficients of the NIST database. This table is filtered to only those lines that fall within the requested wavelength range. Each absorption line corresponds to a transition from a lower energy level of an atom/ion to a higher energy level. To calculate a spectral line, the number density of the lower energy level as a function of radius in the atmosphere is needed. For a given line in the NIST database, sunbather therefore first tries to identify if the lower energy level can be matched to any of the energy levels used by Cloudy. This is not always possible, because we use the default Cloudy energy level selection, which for most species includes the 15 lowest levels. Additionally, Cloudy will sometimes “fold” multiple energy levels together into one, or include energy levels that simply cannot be matched to the levels reported by NIST. In those cases, the spectral line cannot be calculated. In practice, however, the strongest lines in a transit spectrum originate from the lowest few energy levels of an atom or ion, which are typically included in the default Cloudy selection. Therefore, we choose to not increase the number of energy levels included in Cloudy beyond the default selection, as this would increase the computation time. For the energy levels that can be matched to NIST, Cloudy provides the number density of that energy level as a function of radius in the atmosphere.

We then project the 1D number density of the energy level $n_{X}(r)$, temperature $T(r),$ and velocity $v(r)$ profile onto a 2D semi-circle with the $x$-coordinate along the direction of the stellar light rays and the $y$-coordinate indicating the impact parameter of the light rays with respect to the planet surface. The velocity is directed radially outward, but only the $v_{x}$-component affects the Doppler shift of the observed line. By default, sunbather uses 100 rays, but this number can be changed if higher or lower precision is needed. The code then loops over each individual spectral line and calculates the cumulative optical depth of each of the 100 rays as a function of frequency $\nu$ as

$$\tau_{y}(\nu)=\int_{x}n_{X}(x,y)\sigma_{0}\Phi(\Delta\nu)\mathrm{d}x,$$ (7)

where $\sigma_{0}=2.654\times 10^{-2}\times f_{12}$ is the absorption cross-section of the given spectral line. $f_{12}$ is the oscillator strength of the line given by the NIST database. $\Phi$ is the Voigt line profile, which is a convolution of a Lorentzian line profile due to natural broadening and a Gaussian line profile due to thermal broadening. The Lorentzian profile has a (constant) half-width at half-maximum of $\gamma=A_{21}/4\pi$, where $A_{21}$ is the Einstein coefficient of the line given by the NIST database. The Gaussian profile has a (temperature-dependent) half-width at half-maximum of

$$\alpha=\sqrt{\frac{2\mathrm{ln}2\;kT(x,y)}{m_{X}}}\frac{\nu_{0}}{c},$$ (8)

where $m_{X}$ is the mass of the absorbing species and $v_{0}$ is the rest-frame frequency of the line (i.e., the line center). The line profile is dominated by thermal broadening around the line center, but natural broadening can dominate the absorption in the far wings of the line. If the Voigt profile itself is centered around zero, it is evaluated at

$$\Delta\nu=(\nu-\nu_{0})-\frac{v_{x}(x,y)}{c}\nu_{0}$$ (9)

in order to take into account the deviation from the line center at the current frequency, as well as the Doppler shift due to the outflow velocity. We then sum the $\tau_{y}(\nu)$ contributions of all spectral lines together.

The last step is to convert these optical depths of the 100 rays into a transit spectrum. First, sunbather adds 50 more rays between $y=0$ and $y=R_{p}$ with $\tau_{y}(\nu)=\infty$ to take into account the opaque planet core. Each ray $i$ then corresponds to a projected annulus on the stellar disk (the $y-z$-plane) of area $S_{i}=\pi(y_{i+1}-y_{i})^{2}$. However, the background stellar flux varies across such an annulus, based on the stellar limb darkening and the projected planet position (determined by its impact parameter $b$ and orbital phase $\phi$). Each annulus is thus divided into 500 cells distributed spherically along the annulus. In each of these cells, the stellar intensity $I(\nu)$ is evaluated according to the (frequency-dependent) quadratic limb darkening law. Taking the average of the 500 cells results in an average intensity $\bar{I}_{i}(\nu)$ per ray. We then calculate the average intensity of the complete stellar disk $\bar{I}_{*}(\nu)$ by sampling it in 1000 annuli. Finally, the total transit spectrum can be calculated, defined as the ratio of the flux in-transit $F_{in}(\nu)$ and the flux out-of-transit $F_{out}(\nu)$:

$$\frac{F_{in}(\nu)}{F_{out}(\nu)}=1-\sum_{i}\frac{\bar{I}_{i}(\nu)\cdot S_{i}}{\bar{I}_{*}(\nu)\cdot\pi R_{*}^{2}}\Big{[}1-\mathrm{exp}(-\tau_{y_{i}}(\nu))\Big{]},$$ (10)

where $R_{*}$ is the stellar radius.

Having discussed the main modules of sunbather, we provide a clarification of how the planet radius $R_{p}$ is defined. The isothermal Parker wind keeps the mass-loss rate constant at all radii, from a singular point at $r\to 0$ until $r\to\infty$. In reality, the planetary outflows are being launched at some finite altitude between the lower atmosphere, which is in hydrostatic equilibrium, and the upper atmosphere, which is well described by a Parker wind profile. Our models do not capture this transition, and so the Parker wind profile is truncated at a lower boundary of $1R_{p}$. In calculating the transmission spectrum, we assume that the planet is fully opaque at all wavelengths for $r<R_{p}$.

3 Escaping atmospheres at supersolar metallicity

We employ sunbather to investigate how atmospheric escape depends on the composition. The dependence on the helium abundance in hydrogen/helium outflows has been studied before, and the degeneracy with the extreme UV flux of the host star is well known (Lampón et al., 2021b; Khodachenko et al., 2021; Czesla et al., 2022; Fossati et al., 2022; Xing et al., 2023; Biassoni et al., 2024; Allan et al., 2024). Even though atmospheric escape models that incorporate metals have been employed in the literature (e.g., Salz et al., 2016; Shaikhislamov et al., 2020; Huang et al., 2023), a systematic exploration of the influence of metals has so far only been performed by Zhang et al. (2022) and Kubyshkina et al. (2024).

Studying how metals impact atmospheric mass-loss is especially important in the context of small gaseous planets. In the solar system, a clear anticorrelation between planet mass and atmospheric metallicity has been observed (Guillot et al., 2022). However, establishing a similar relationship for exoplanets proves more challenging. While there is evidence supporting an inverse correlation between planet mass and bulk metallicity (Thorngren et al., 2016), it is not yet clear how this relates to the atmospheric metallicity (e.g., Welbanks et al., 2019; Bean et al., 2023).

In this section, we investigate how the atmospheric bulk properties such as density, velocity and temperature respond to an increasing metal content. We also study how the changing atmospheric structure and abundances in turn affect the transit spectrum. A front-to-back reproduction script for the results of this section is available on Zenodo${}^{\ref{fn:zenodo}}$. The simulations presented here neglected the stellar tidal gravity term of Eq. 2, but we checked that the results are qualitatively the same as simulations including it. We simulate a generic hot-Neptune planet with $M_{p}=0.1$ $M_{J}$, $R_{p}=0.5$ $R_{J}$, orbiting a Sun analog at $a=0.05$ AU. The solar spectrum that we use is a combination of different datasets. The 0.5-189.5 nm wavelength region consists of data from the SEE instrument aboard the TIMED spacecraft, which has been measuring the full solar disk irradiance since 2002 (Woods et al., 2005). We use the Level 3 data product, which constitutes one-day average spectra with flares removed. The data we use is from January 1st 2016. This date is in-between the maximum and minimum of the 11-year solar cycle and was thus assumed to represent a somewhat average spectrum. The 190-310 nm and 310-2412 nm wavelength regions consist of data from the SOLSTICE and SIM instruments, respectively, aboard the SORCE spacecraft that also continuously monitored the solar spectrum (Rottman, 2005). These components are again daily averages at January 1st 2016. We finally extrapolate the red end of the spectrum with a Rayleigh-Jeans power-law until 10,000 nm. The spectrum is shown in Fig. A.1 of Linssen & Oklopčić 2023.

We simulate the planet using six different atmospheric compositions: a pure hydrogen/helium mixture in the ratio 10:1 by number, a complete solar composition including metals, and three, ten, 30, and 50 times solar metallicity. A scaled solar metallicity $Z$ in this case means we scale the metal/hydrogen abundance ratio relative to the solar ratio as

$$Z=\Bigg{(}\frac{X_{m}}{X_{H}}\Bigg{)}\;\Bigg{/}\;\Bigg{(}\frac{X_{m,\odot}}{X_{H,\odot}}\Bigg{)},$$ (11)

where $m$ runs over every metal element. With this definition, $Z$ can take any value, but the increase in metal abundance becomes nonlinear at high $Z$. To illustrate, while the solar carbon abundance is $X_{C}=2.225\times 10^{-4}$, at 50 times solar metallicity, the abundance increases to $X_{C}=1.063\times 10^{-2}$, which is only a factor $\sim$48 increase. The carbon abundance converges towards $X_{C}=0.233$ as $Z\to\infty$.

3.1 Assessing the validity of a constant composition

As explained in Sect. 2, sunbather does not treat mass fractionation but instead has a constant composition with altitude. To clarify, the ionization fractions and energy level populations of each element do change as a function of radius, but the overall abundance of the element does not. We can assess the validity of the assumption of constant composition in the high-metallicity outflows by using the concept of the crossover mass, as introduced by Hunten et al. (1987). In an atmosphere consisting of two constituents, the crossover mass is defined as

$$m_{c}=m_{1}+\frac{kTF_{1}^{0}}{g_{0}X_{1}b},$$ (12)

where $m_{1}$ is the mass of the first constituent, $k$ is the Boltzmann constant, $T$ is the temperature, $F_{1}^{0}$ is the number outflow flux of the first constituent, $g_{0}=GM_{p}/R_{p}^{2}$ is the planet surface gravitational acceleration, $X_{1}$ is the abundance of the first constituent, and $b$ is the diffusion parameter, given by (Cussler, 2009):

$$b=n\cdot\frac{AT^{3/2}}{P\sigma_{12}^{2}\Omega}\sqrt{\frac{1}{M_{1}}+\frac{1}{M_{2}}},$$ (13)

where $n$ is the number density, $P$ is the pressure, $\sigma_{12}=(\sigma_{1}+\sigma_{2})/2$ is the average collision diameter of the two constituents, $\Omega$ is a dimensionless constant of order unity that we take to be 1 here (Hirschfelder et al., 1967), $M_{1}$ and $M_{2}$ are the molar masses of the two constituents, and

$$A=\frac{3}{8}k^{3/2}\sqrt{\frac{N_{A}}{2\pi}},$$ (14)

where $N_{A}$ is the Avogadro constant. If we take atomic hydrogen as constituent 1, assuming it drives the outflow and drags along the second constituent (a metal species), the crossover mass can be understood as the particle mass above which the metal species is no longer effectively dragged along. In other words, if $m_{2}<<m_{c}$, the metal species is dragged along efficiently, resulting in a constant composition with altitude. However, as $m_{2}$ approaches $m_{c}$, the metal escape flux approaches 0.

We can thus check our model assumption of efficient metal drag by estimating the crossover mass at $r=R_{p}$. We do this for the model that requires the most extreme metal drag, namely the 50 times solar metallicity case with Parker wind parameters $T_{0}=5000$ K and $\dot{M}=10^{9.9}$ g s^-1 (see Sect. 3.2 for the choice of parameters). With such a composition, $X_{1}=0.87$ and the hydrogen number flux is $F_{1}^{0}=1.3\times 10^{13}$ s^-1. In principle, our atmosphere is a mixture of many different constituents, but Eq. 12 is derived under the assumption of only two constituents. We take the second constituent to be atomic iron, with $M_{2}=55.845$ g mol^-1. The collision diameters are not readily available, and so we assume $\sigma_{12}=2.58$ Å (similar to the value of various molecular gases, see Hirschfelder et al. 1967). Using the temperature, pressure and number density of the Parker wind profile at $R_{p}$, we find $b=1.44\times 10^{20}$ cm^-1 s^-1 and $m_{c}=43$ amu. This means that the assumption of efficient metal drag and a constant composition with altitude is questionable for the heaviest metal species in this particular simulation (zinc is the most massive particle included, with a a mass of 65 amu). The simulations at lower metallicities have higher crossover masses ($m_{c}=127$ amu at 30 times solar metallicity) and metal drag can thus be assumed to be efficient in them. Even for the 50 times solar metallicity case however, our calculation used hydrogen and iron in atomic form as the two constituents. For ions, the collision diameter is much stronger (Hunten et al., 1987), and consequently the crossover mass is much larger. The Cloudy simulation shows that both hydrogen and iron (and other metal species for that matter) are indeed mostly ionized already at $r=R_{p}$. If we re-calculate the crossover mass for an ionized outflow using the more generalized description given in Eq. 33 of Xing et al. (2023), we find $m_{c}\approx 5\times 10^{5}$ amu. Hence, we consider the assumption of efficient metal drag to be reasonable for the 50 times solar metallicity simulations as well.

3.2 Effect of metallicity on the outflow structure

3.2.1 The model parameters

When comparing models with different metallicities, we have to choose appropriate values for the model free parameters $T_{0}$ and $\dot{M}$. We first explore how the atmospheric structure changes when using the same mass-loss rate of $\dot{M}=10^{11}$ g s^-1 in each model. To find a suitable $T_{0}$, we perform a model self-consistency check. This means we run multiple models with different values for $T_{0}$. For each, we calculate the average of the nonisothermal temperature profile between $r=R_{p}$ and $r=2R_{p}$ and compare it to $T_{0}$. We select the model where the difference is the smallest, because the nonisothermal temperature profile is then most consistent with the density and velocity profiles which were derived under the isothermal approximation. We refer the interested reader to Linssen et al. (2022) and Linssen & Oklopčić (2023) for a more in-depth discussion on assessing the self-consistency of a Parker wind model. We find self-consistent values for $T_{0}$ of 5100, 5000, 5000, 4900, 4700, 4100 K for the six models ordered by ascending metallicity.

Figure 3: Atmospheric profiles for a generic hot-Neptune planet at different metallicities. The left column shows a set of models with a fixed mass-loss rate. The right column shows another set of models with a fixed density at the planet radius, leading to different mass-loss rates at different metallicities. For every model, the value of $T_{0}$ is chosen such that it is consistent with the mean of the nonisothermal temperature profile of panel (d) or (i) between 1$R_{p}$ and 2$R_{p}$. Those values for $T_{0}$ can be found in the text in Sect. 3.2.1 and Sect. 3.2.4. In panels (e) and (j), solid lines show the radiative heating rates and dashed lines show the radiative cooling rates.

3.2.2 Effect of metallicity on the hydrodynamics

The left column of Fig. 3 shows the bulk atmospheric properties of the six models. Panels a, b, and c show the density, velocity and mean particle mass profiles, respectively, which are calculated under the assumption of an isothermal Parker wind. Unsurprisingly, the mean particle mass increases with the metallicity of the atmosphere. This difference in mean particle mass then propagates into different density and velocity profiles through the sound speed (Eq. 4). As the sound speed becomes lower, the outflow velocity decreases and depends more steeply on radius. Since the mass-loss rate $\dot{M}=4\pi r^{2}\rho(r)v(r)$ is kept constant, the density shows the reverse behaviour and increases while also depending more steeply on radius. These changes only become substantial at metallicities higher than three times solar. The conclusion that the mean particle mass is the main driver of structural changes is in agreement with the findings of Kubyshkina et al. (2024), who compare a hydrogen, a hydrogen/helium, and a solar composition outflow.

3.2.3 Effect of metallicity on the temperature

Panel d shows the nonisothermal temperature profiles that we obtain by solving Eq. 3. We generally observe a radial flattening of the temperature profile with increasing metallicity. We can understand this by considering the “thermostat temperature”, which is the temperature that the gas would have in radiative equilibrium (in the absence of the expansion and advection terms). The thermostat temperature has a rather weak dependence on density (and hence atmospheric radius) and is $\sim$10,000 K for the pure hydrogen/helium atmospheric density profile of panel d (agreeing with the findings of Murray-Clay et al., 2009). Due to metal cooling, the thermostat temperature steadily decreases towards $\sim$5000 K at 50 times solar metallicity. For the pure hydrogen/helium and lower metallicity atmospheres, the thermostat temperature is not actually reached because of expansion cooling, and we instead find a temperature profile that decreases with radius. On the other hand, for the high metallicity atmospheres, the expansion cooling rates are much lower because of the lower outflow velocities. The high metallicity outflows therefore do reach their thermostat temperature, and we find a temperature profile that is $\sim$5000 K with little dependence on radius. This explanation is supported by the radiative heating and cooling rates shown in panel e. For the lower metallicity models, the radiative heating rate is roughly an order of magnitude higher than the radiative cooling rate. The radiative heating is instead balanced by expansion cooling. For the 50 times solar metallicity model, we see that radiative heating and cooling balance each other throughout the radial domain. The models at ten and 30 times solar metallicity are (almost) in radiative equilibrium only at lower radii, and this is indeed where their temperature profiles are relatively flat. We can thus comment on the thermal properties of the gas by looking at the difference between the radiative heating and cooling rates. At high metallicity, thermal balance between radiative processes indicates that chemistry plays a dominant part in governing the outflow. At lower metallicity, however, the onset of expansion cooling highlights the relevance of dynamics in affecting the atmosphere, especially at higher radii.

We also note that it is somewhat coincidental that the low metallicity models have a temperature of $\sim$5000 K at $r\lesssim 2R_{p}$, which is similar to the high metallicity thermostat temperature. This is because expansion cooling dominates throughout the atmosphere of the hot-Neptune planet that we simulate here. As identified in for example Salz et al. (2016); Caldiroli et al. (2021); Linssen et al. (2022), a higher gravity planet like a hot Jupiter is in radiative equilibrium at lower radii, and would therefore reach a temperature close to the thermostat value of $\sim$10,000 K.

3.2.4 Effect of metallicity on the mass-loss rate

It is interesting to examine how the mass-loss rate of the planet changes as a function of atmospheric metallicity. Unfortunately, we cannot easily answer this question, since the mass-loss rate is assumed, rather than predicted, in our model. However, under the assumption that the density at the planet radius stays the same irrespective of the outflow metallicity, we can adjust the mass-loss rate of the models with different compositions in such a way that we respect that density. This assumption can be debated, as modeling efforts have shown that metallicity affects the conditions at the photosphere (e.g., Mollière et al., 2015; Drummond et al., 2018). However, we are only aiming to do a simple exploration of the problem here. With this exercise, we are effectively imposing a boundary condition on $\rho(R_{p})$ that our isothermal Parker wind model normally does not have. Other atmospheric escape models that treat the hydrodynamics and thermal balance self-consistently also have such a boundary condition. In those models, the chosen boundary conditions can affect the resulting mass-loss rate ranging from a few tens of percent up to a factor of a few, depending on the planet considered (e.g., Salz et al., 2016; Kubyshkina et al., 2024). We should thus take into account that changes to the mass-loss rate at that level could be attributed to the boundary condition rather than metallicity.

We keep the Parker wind parameters of the pure hydrogen/helium composition the same as before, that is, a self-consistent temperature of $T_{0}=5100$ K and a mass-loss rate of $\dot{M}=10^{11}$ g s^-1. For the simulations with different metallicity, we then find the mass-loss rate for which $\rho(R_{p})$ is the same as that of the hydrogen/helium profile, while using a self-consistent $T_{0}$. This yields the following values: $T_{0}=5100$ K and $\dot{M}=10^{10.95}$ g s^-1 for the model with a solar composition, $T_{0}=5100$ K and $\dot{M}=10^{10.95}$ g s^-1 for the model with three times solar metallicity, $T_{0}=5100$ K and $\dot{M}=10^{10.8}$ g s^-1 for the model with ten times solar metallicity, $T_{0}=5000$ K and $\dot{M}=10^{10.3}$ g s^-1 for the model with 30 times solar metallicity, and $T_{0}=5000$ K and $\dot{M}=10^{9.9}$ g s^-1 for the model with 50 times solar metallicity. Hence, for a fixed density $\rho(R_{p})$, the mass-loss rate is roughly an order of magnitude lower at high metallicity than at solar metallicity.

The right column of Fig. 3 presents the atmospheric structure of these models. Since the mass-loss rates of the lower metallicity models are close to $10^{11}$ g s^-1, their atmospheric structures are almost the same as in the left column of Fig. 3. The differences are larger for the higher metallicity models. The temperature profiles shown in panel i are now all remarkably similar. Like in Sect. 3.2.3, we can explain this with an argument based on the expansion cooling rate. The velocity and hence expansion cooling rate of the higher metallicity models is relatively low at radii $r\lesssim 2R_{p}$. As panel j confirms, the result is that these models are in radiative equilibrium at their thermostat temperatures of 5000 K in this region of the atmosphere. Contrary to the previous case at fixed $\dot{M}=10^{11}$ g s^-1, their velocities are now close to the velocities of the solar composition model at radii $r\gtrsim 2R_{p}$. Due to expansion cooling, their temperature profiles are therefore not flat, but decrease with radius.

3.2.5 Radiative heating and cooling by metals

We take a closer look at the different heating and cooling processes that make up the total radiation rates. For the simulations with a fixed density at $R_{p}$ (right column of Fig. 3), Fig. 4 shows the contribution of each process, as given by Cloudy. We present these as a fraction of the total radiative heating or cooling rate. At solar composition, heating is dominated by hydrogen ionization, but helium ionization contributes significantly, especially at smaller radii. Metals only slightly affect the heating through line heating, which makes up $\lesssim 5$%. Cooling is dominated by hydrogen recombination and free-free emission at radii $r\gtrsim 2R_{p}$. At smaller radii, line cooling from iron, magnesium, and calcium ions becomes important and comparable to the hydrogen recombination cooling rate. Interestingly, cooling through hydrogen Lyman-$\alpha$ emission seems to be unimportant, contrary to the findings of other works both employing Cloudy (e.g., Salz et al., 2016; Zhang et al., 2022; Kubyshkina et al., 2024) and not employing Cloudy (e.g., Murray-Clay et al., 2009; Allan et al., 2024). It is not immediately clear to us why this is the case. As we go to ten times solar metallicity, heating is still mostly dominated by hydrogen and helium ionization, but metal line heating and oxygen ionization is becoming more important. Cooling below $\lesssim 2R_{p}$ is now completely dominated by iron, magnesium, and calcium line cooling. At 50 times solar metallicity, heating occurs through wealth of different agents. At smaller radii, hydrogen, helium, and carbon ionization still dominate, but above $2R_{p}$, ionization of various metal atoms and ions all contribute a few % to 20 % of the heating. The cooling rate qualitatively shows roughly the same behaviour as at ten times solar metallicity. Iron, magnesium, and calcium line cooling dominate at smaller radii, while at larger radii cooling occurs through the lines of various metal ions, the strongest being Ne²⁺. These results show that metals are important in the radiative budget of the planet already at solar abundance through their line cooling. At enhanced metallicities of approximately ten times solar and above, their ionization and line heating also become important.

Since Fig. 4 shows the fractional contribution of each heating and cooling agent, we need to compare with Fig. 3j to get a sense for the absolute change in the heating and cooling rates of the different agents. This reveals that at radii $\gtrsim 1.5R_{p}$, the lower hydrogen and helium ionization heating fractions at high metallicity indeed also correspond to a lower radiative heating rate. In other words, in those regions, the decrease in hydrogen/helium heating is not replaced by similar amounts of metal ionization heating. At radii $\lesssim 1.5R_{p}$ however, the radiative heating rate does increase with metallicity due to O⁺ ionization. For the cooling rates, we see similar behaviour, but now at a dividing radius of $\sim 3R_{p}$. At radii $\lesssim 3R_{p}$, the loss of hydrogen recombination cooling is not counterbalanced by metal cooling. However, at radii $\lesssim 3R_{p}$, the onset of iron, magnesium, and calcium line cooling does lead to an overall much higher cooling rate.

Figure 4: Radiative heating and cooling agents as given by Cloudy for models with a fixed density at the planet radius (right column of Fig. 3). The rates are expressed as a fraction of the total radiative heating or cooling rate. The grey dotted line labeled “other” is the sum of all other agents that are not plotted as they do not exceed a fraction of 0.01. For some processes such as heating by “metal lines”, we unfortunately cannot obtain more details from the Cloudy output about the specific elements responsible. Please note the log-scale when interpreting the relative importance of the different agents.

3.3 Effect of metallicity on the transit spectrum

To investigate the effect of the outflow metallicity on the transit spectrum, we calculated the transit depths of various atmospheric escape tracers. We first calculated the transit depths of all spectral lines identified and listed in Table C.1 of Linssen & Oklopčić (2023). From these, we selected the 18 strongest lines and supplemented this with the metastable helium line and calcium infrared triplet, to obtain a sample of 20 spectral lines. Fig. 5 presents the transit depths of these spectral lines as a function of outflow metallicity, for the simulations with a fixed photospheric density (right column of Fig. 3). The spectra were made at mid-transit and without stellar center-to-limb variations.

The hydrogen Lyman-$\alpha$ line appears to have a transit depth independent of metallicity. This is an artifact of our atmospheric radial grid, which extends to 8$R_{p}$ so that the maximum transit depth is $(8R_{p}/R_{*})^{2}=16.9\%$. Even though there is certainly less absorption in the Lyman-$\alpha$ line at high metallicity, the line is still optically thick and hence its transit depth does not decrease. A model that extends to larger radii (and takes into account interactions with the stellar environment) would be needed to properly explore the dependence of the Lyman-$\alpha$ line on the outflow metallicity. Since the metastable helium line forms at smaller radii, it is seen to become weaker at high metallicity. This is not due to the decrease in helium abundance, which is only $\sim 4\%$, but due to structural changes in the outflow. The dominant factor is that the total density (hence also helium density) decreases because the mass-loss rate gets lower. At the same time, Fig. 3 shows that the temperature decreases at radii $r\lesssim 2R_{p}$ and increases at radii $r\gtrsim 2R_{p}$ at higher metallicity, leading to an increase and decrease in the fraction of helium atoms in the metastable state, respectively (see Fig. 1 of Biassoni et al., 2024). The combined result turns out to be a weaker helium line.

For all the other lines that originate from metals, we notice an interesting response to an increasing metallicity. Spectral lines originating from atomic or singly ionized metals become stronger until ten to 30 times solar metallicity and then become weaker again at 50 times solar metallicity. In contrast, spectral lines originating from doubly or higher ionized metals keep growing stronger until 50 times solar metallicity. This behavior is the result of a decreasing number density of atoms and singly charged ions on the one hand and an increasing number density of doubly and higher charged ions on the other hand. We can understand why the ionization fractions change in this way as follows. If we temporarily consider a total density profile that is kept constant, increasing the metal content will lead to a lower free electron density. This is due to the fact that a metal particle generates fewer free electrons per unit mass than hydrogen or helium - assuming the metal particle is not completely ionized, which is indeed not the case in the conditions of an upper exoplanet atmosphere. The lower electron density leads to a lower recombination rate of the ions, as the radiative recombination rate depends linearly on the electron density (Rybicki & Lightman, 1991) and the three-body recombination rate depends on the electron density squared. The photoionization rate stays roughly the same, and as a consequence the degree of ionization increases. At 50 times solar metallicity, this effect is so strong that we get a lower number density of metal atoms and singly charged ions, even though the overall metal abundance is higher than in the 30 times solar metallicity case. Vice versa, the number density of doubly and higher charged ions grows faster than the increase in overall metal abundance. Finally, we can consider the case where the total density profile is not kept constant. As the top-right panel of Fig. 3 shows, the total density in fact decreases with metallicity, which adds to the weakening of the lines.

Figure 5: Depths of the 18 strongest lines in the transit spectra, as well as the metastable helium and calcium infrared triplet lines. The results are shown for the models with a fixed density at the planet radius (right column of Fig. 3). Our computational radial domain extends to 8$R_{p}$, corresponding to a maximum absorption of 16.9%, which is reached by various spectral lines. Lines originating from metal atoms (solid) and singly charged ions (dashed) tend to become stronger up to ten times solar metallicity, and then become weaker. Lines originating from doubly or higher charged metal ions (dotted and dash-dotted) keep growing stronger as the metallicity increases.

4 Constraining mass-loss rates from spectral line observations

In this section, we demonstrate one of the main use cases of sunbather: fitting an observed transit spectrum with a grid of models in order to constrain the planet’s mass-loss rate. Measuring the mass-loss rate is especially interesting as it is a fundamental parameter in the atmospheric evolution. Obtaining accurate mass-loss rates for individual planets will help to reveal the role of atmospheric escape in shaping the sub-Jovian desert and radius valley (e.g., Owen & Lai, 2018; Owen & Wu, 2013). Constraining the mass-loss rate with sunbather is enabled by the fact that it is a free parameter in the Parker wind description, and can thus be adjusted until the model fits the observations. With this in mind, many works have indeed used an isothermal Parker wind model to fit observations of the metastable helium line (see references in the introduction). In most cases, the fit is performed only for $T_{0}$ and $\dot{M}$ while considering a hydrogen/helium outflow with some assumed helium abundance. However, some of these works have additionally fit for the helium abundance (e.g., Lampón et al., 2020; dos Santos et al., 2022), but it is typically degenerate with other model parameters. In Sect. 4.1, we start by demonstrating the method on a mock spectrum, which serves as a simple example. We then move on to an observed spectrum of a real planet in Sect. 4.2, which serves as a more complicated example. In Sect. 4.3, we explore the influence of the assumed metal abundance on the fit, and show that even though we keep the hydrogen/helium ratio the same, the presence of metals in the outflow still affects the mass-loss rate that is inferred from helium observations. A front-to-back reproduction script for the results of this section is available on Zenodo${}^{\ref{fn:zenodo}}$, and it can easily be modified to interpret other observations.

4.1 Fitting a mock spectrum

We start by generating a fake observed metastable helium spectrum. We take the generic hot-Neptune planet that we introduced in Sect. 3 and choose a solar composition atmosphere, $T_{0}=5100$ K and $\dot{M}=10^{10.95}$ g s^-1. We first use sunbather to calculate the metastable helium spectrum between 10,830 and 10,836 Å, shown as the purple line in Fig. 6a. We then turn this synthetic spectrum into a mock observed spectrum by introducing noise. We sample the synthetic spectrum at a spectral resolution of 80,000. At each wavelength point, we add random values drawn from a normal distribution with standard deviation 0.0025 to the model $F_{in}/F_{out}$ value. We then also assign an observed errorbar of 0.0025 to the drawn data points. This results in the mock spectrum shown in black in Fig. 6a.

Figure 6: Demonstration of a retrieval of the mass-loss rate and temperature from the helium spectrum of the generic hot-Neptune introduced in Sect. 3. a) True model metastable helium spectrum and the mock observed spectrum. The best retrieved model from the resolved line fit is also shown. b) Results of the fit to the spectrally resolved line, without using a prior. c) Results of the fit to the EW of the line. As visible from the blue contours, the fit offers no tight constraints on its own. We thus place a normal prior on the parameter space based on the model temperature self-consistency, shown in orange. The posterior distribution shown in black then yields a precise and accurate constraint.

We now pretend that this is an observed spectrum of the planet and aim to fit for the mass-loss rate and temperature. We fix all other parameters (see Fig. 1 for the parameters that are needed) at the values that we used to create the mock spectrum, thus effectively assuming we know them precisely. For $R_{p}$, $a$, and $R_{*}$, this is a reasonable assumption as for real planets they are usually known well enough to not dominate the uncertainties of the analysis. The value of $M_{p}$ typically has larger uncertainties. To get a sense for when the uncertainty in $M_{p}$ may become significant, we ran an identical model as that used to generate the mock spectrum, but increased the planet mass by 30%. The resulting helium spectrum was within $1\sigma$ of the original spectrum. When increasing the planet mass by 60%, the difference became more significant at $\sim 3\sigma$. The final two “parameters” that we fix at the real values are the stellar spectral energy distribution (SED) and the planet atmospheric composition, which are typically not known accurately for real planets, but may also significantly influence the results of the analysis. One may thus wish to also explore different SEDs and compositions, in addition to $T_{0}$ and $\dot{M}$ (we explore this in more detail in Sect. 4.2 for the atmospheric composition).

We run a grid of models, with temperatures ranging from 4000 to 7000 K in steps of 100 K, and mass-loss rates ranging from 10^10.5 to 10^11.5 g s^-1 in steps of 0.05 dex (totaling 651 models). Each model calculation consists of first calling the construct_parker.py routine to calculate the isothermal Parker wind density and velocity profiles (Sect. 2.1), then calling the convergeT_parker.py routine to calculate the nonisothermal temperature profile (Sect. 2.2), and finally using the RT.py module to calculate the metastable helium spectrum (Sect. 2.3).

With all model spectra calculated, we can turn to the statistics. We first resample each model spectrum onto the same wavelength grid as the observed spectrum. Assuming normally distributed errors, we can fit each model to the data using a $\chi^{2}$-statistic. In a Bayesian framework, these $\chi^{2}$-values are then related to the likelihood $L$ as $L=\mathrm{exp}(-\chi^{2}/2)$, which we show in Fig. 6b. One could choose to also use a (nonflat) prior, but we abstain from doing so here since the fit itself is already quite constraining. The purple cross in Fig. 6b marks the true parameters that were used to generate the mock spectrum. It falls within the 2$\sigma$ contour of the likelihood, confirming that our retrieval has been successful.

To demonstrate a case in which we have spectrally unresolved data, we now repeat the analysis but fit only the equivalent width (EW) of the line. Observations with the James Webb Space Telescope (dos Santos et al., 2023a), the Hubble Space Telescope (Spake et al., 2018), and the Palomar/WIRC ultranarrowband filter (Vissapragada et al., 2020) would produce such low-resolution or unresolved metastable helium spectra. We integrate our true model spectrum to obtain its EW and then again randomly sample a value around it, obtaining a mock observed EW=$13.5\pm 1.8$ mÅ (this value is thus not related to the mock resolved spectrum we generated before). We can then simply compare the EW of each model to this observed value, which directly gives the likelihood (shown in blue in Fig. 6c). As visible from the figure, the fit itself is not very constraining, and we instead find a degeneracy between $T_{0}$ and $\dot{M}$.

To mitigate this problem (which is common when fiting Parker wind models to spectrally unresolved data), we follow the approach outlined in Linssen et al. (2022). This comes down to now using a nonflat prior to break the degeneracy. We choose a prior that is based on the self-consistency of the atmospheric structure. Specifically, for each Parker wind model we compare the nonisothermal temperature profile $T(r)$ to the isothermal value $T_{0}$. A large difference between the two indicates that the model is not very self-consistent. To compare the $T(r)$ profile to a single value $T_{0}$, we need to somehow collapse $T(r)$ into a single value as well. Since we are fitting metastable helium observations, we care most about the self-consistency of the model in the region where this line forms. Therefore, we calculate the mean and standard deviation of $T(r)$ weighted by the metastable helium density, resulting in $T_{He}$ and $\sigma_{T}$, respectively (Eqs. 4 and 5 of Linssen et al. 2022). Our prior value for the model is then given by a normal distribution centered on $T_{He}$ with standard deviation $\sigma_{T}$ evaluated at $T_{0}$. The orange line in Fig. 6c show models for which $T_{He}=T_{0}$, while the dark and light orange contours indicate models where $T_{0}$ falls within 1$\sigma_{T}$ and 2$\sigma_{T}$ of $T_{He}$.

The posterior distribution now follows by multiplying the prior and likelihood and normalizing it. In Fig. 6c, the solid and dashed black lines indicate the 1$\sigma$ and 2$\sigma$ contours, respectively. We see that the retrieval has been successful, which is in large part due to using the prior. Without it, the fit would have been quite unconstraining for both the temperature and mass-loss rate. We finally note that in the case of this mock spectrum, the temperature self-consistency prior that we used in the equivalent width fit is also consistent with the temperature that we constrained directly from the resolved line fit. As we show in the following section, this is not always the case for real spectra.

4.2 Fitting the helium spectrum of TOI-2134 b

We now turn to a real observed spectrum, which will prove to be more complicated to interpret. We use this example to illustrate how we can use sunbather to derive insights into outflow physics when our current modeling struggles to self-consistently match the data.

We model the helium line of the mini-Neptune TOI-2134 b as observed with Keck/NIRSPEC by Zhang et al. (2023a). The original spectrum showed a $\sim$5 km s^-1 redshift from the rest-frame helium triplet wavelength. Such offsets can have valuable physical implications (e.g., Nail et al., 2024), but cannot be reproduced by sunbather, which by virtue of its 1D nature yields lines centered on the rest-frame wavelength. So as to not obtain biased fit results, we first shift the spectrum back to the rest-frame wavelength (shown in Fig. 7c). We then adopt the system parameters from Zhang et al. (2023a): $R_{p}=0.24$ $R_{J}$, $M_{p}=0.0287$ $M_{J}$, $R_{*}=0.709$ $R_{\odot}$, $a=0.078$ AU, $b=0.2$, and use their stellar SED that is reconstructed from XMM-Newton observations.

We first run models assuming a pure hydrogen/helium composition in the solar ratio of 10:1 by number. The helium abundance is thus $X_{He}=0.091$. We run a grid with $2000<T_{0}<8000$ K and $10^{8.5}<\dot{M}<10^{10}$ g s^-1. The simulations excluded the stellar tidal gravity term of Eq. 2, but we verified that the differences with models including it are minimal for this planet. We calculate the helium spectrum for each model, convolve it to the instrument resolution of $R=32,000$, and fit this to the observed spectrum using a $\chi^{2}$-statistic. The associated likelihood contours are shown in blue in Fig. 7a. We find $\dot{M}=10^{9.28^{+0.06}_{-0.04}}$ g s^-1 and $T_{0}=5000\pm 400$ K and the best-fit model is shown as the solid blue line in Fig. 7c. Zhang et al. (2023a) also modeled their own observations using p-winds and found $\dot{M}=10^{9.74\pm 0.08}$ g s^-1 and $T_{0}=4640\pm 230$ K. The difference with our constraints likely stems from the fact that sunbather calculates a nonisothermal temperature structure. The best-fit models with $T_{0}\approx 5000$ K have nonisothermal profiles that are quite a bit colder than 5000 K. At those lower temperatures, the metastable helium fraction becomes higher, and a lower mass-loss rate is therefore needed to match the observed line (e.g., Biassoni et al., 2024).

Figure 7: Two fits to the metastable helium data of the mini-Neptune TOI-2134 b from Zhang et al. (2023a) assuming different atmospheric compositions. In panel (a), we assume a pure hydrogen/helium atmosphere in the solar ratio. In panel (b), we assume a 100 times solar metallicity atmosphere. The colors are the same as in Fig. 6, but here we do not use the model temperature self-consistency as a prior as we did in Fig. 6c. c) Observations from Zhang et al. (2023b) shifted back to the rest-frame helium triplet wavelength (black). Solid lines show the best-fit models from the blue contours in panels (a) and (b). Dashed lines show models with self-consistent $T_{0}$-values from the orange contours in panels (a) and (b). The self-consistent models produce lines that are too narrow, indicating that an additional line broadening mechanism may be at play.

We also calculate the model self-consistency in the same way as in Sect. 4.1. We find a self-consistent temperature of $T_{0}\sim 2600$ K, independent of the mass-loss rate (orange contours in Fig. 7a). This value is significantly lower and inconsistent with the temperature obtained from the spectral line fit. Because the difference is so large and the line fit itself already yields tight constraints, we do not use the model self-consistency as a prior here. However, it is interesting to discuss the implications and potential causes of the $T_{0}$-discrepancy. The dashed blue line in Fig. 7c shows the spectrum of a (self-consistent) model with $T_{0}=2600$ K and $\dot{M}=10^{8.8}$ g s^-1 (the value of the mass-loss rate is manually chosen so that the peak absorption is matched). We see that the $T_{0}$-discrepancy in practice means that the observed spectral line is broader than expected from our self-consistent models.

The $T_{0}$-discrepancy can stem from a variety of factors. It may be partly due to the way in which we assess the self-consistent $T_{0}$ in the first place. Using a different criterion than the metastable helium density could lead to a different temperature constraint. We can explore this line of thinking a little further by noting that we would obtain the highest possible self-consistent temperature if we were to compare the peak of the nonisothermal temperature profile to $T_{0}$. With this criterion, we obtain a hard upper limit of $T_{0}\lesssim 4000$ K, which is still somewhat inconsistent with the temperature preferred by the data. The next possibility to consider is that our models underestimate the temperature of the outflow. This may be due to uncertainties in the stellar SED, or additional heating from, for example, Joule heating (Cohen et al., 2024). If applicable, these additional heat sources might increase the self-consistent temperature of the outflow, leading to larger outflow velocities that broaden the line proportionately. Alternatively, it is also possible that the thermodynamics are modeled appropriately and the temperature is indeed around 2600 K, but there are additional kinematic components not considered by sunbather that broaden the velocity distribution of the line. Such extra velocity components could arise due to the stellar tidal gravity, the stellar radiation pressure (e.g., Bourrier & Lecavelier des Etangs, 2013), the stellar wind ram pressure (e.g., MacLeod & Oklopčić, 2022) and its interaction with a planetary magnetic field (e.g., Carolan et al., 2021), atmospheric turbulence and circulation, and the tidally locked rotation of the planet (e.g., Huang et al., 2023).

Ultimately, observing more spectral lines than only the metastable helium line would help to discriminate between the various hypotheses. Other spectral lines can probe different elements, complementary regions of the atmosphere where the temperature and velocity is different, and have varying natural and thermal line widths, all of which may help us understand the apparent inconsistency (Linssen & Oklopčić, 2023).

4.3 Assuming a 100 times solar metallicity composition

We repeat the previous analysis of TOI-2134 b while assuming a 100 times solar metallicity atmosphere. The helium abundance is $X_{He}=0.083$ in this case, which is roughly 10% lower than in Sect. 4.2. The results are presented in Fig. 7b. We find $\dot{M}=10^{9.68\pm 0.04}$ g s^-1 and $T_{0}=9700\pm 700$ K. We again also calculate the model self-consistent temperature, which is $T_{0}\sim 4000$ K.

We can use the findings of Sect. 3 to explain the changes relative to the pure hydrogen/helium fit. As visible from Fig. 7, the fit at 100 times solar metallicity prefers a much higher $T_{0}$. A higher $T_{0}$ leads to higher outflow velocities and it thus compensates for the otherwise lower velocities at high metallicity (see Sect. 3.2.2). The velocity must remain roughly the same because it is an important contributor to the spectral line broadening that we are fitting. The fit at 100 times solar metallicity also prefers a higher $\dot{M}$. This is due to the well-known $T_{0}-\dot{M}$-degeneracy mechanism: at higher $T_{0}$, the overall density and hence metastable helium density decreases, and a higher mass-loss rate is needed to compensate. Focusing on the self-consistent temperature, we see it has increased from $T_{0}\sim 2600$ K for pure hydrogen/helium, to $T_{0}\sim 4000$ K at 100 times solar metallicity. This is because the temperature profile becomes flatter at higher metallicity and approaches the thermostat temperature (see Sect. 3.2.3) of $\sim$4000 K. We thus find a self-consistent $T_{0}$ of that value. We plot the helium spectrum of a model with $\dot{M}=10^{9.1}$ g s^-1 in Fig. 7c (this mass-loss rate is manually chosen to match the peak absorption of the data).

The results of this and the previous section show that the inferred mass-loss rate differs by up to half a dex depending on the assumed outflow metallicity. The helium abundance itself is hardly different, but changes in the outflow structure caused by the high metallicity still affect the metastable helium line. As shown by Lampón et al. (2021b), changing the assumed helium abundance can furthermore change the inferred mass-loss rate by orders of magnitude. Concluding, the retrieved mass-loss rate is degenerate with the helium abundance, and to a lesser extent, the metallicity. Given that we usually do not know the upper atmospheric composition of a given exoplanet, we again pose that we can achieve better mass-loss-rate constraints by combining observations of multiple spectral lines that trace atmospheric escape.

5 Summary

We present sunbather, an open-source${}^{\ref{fn:github}}$ code that can be used to model escaping exoplanet atmospheres and their transit spectra. The code calculates the atmospheric density and velocity profiles assuming an isothermal Parker wind using the p-winds code. A refined nonisothermal temperature profile is then calculated using the photoionization code Cloudy. A radiative transfer module finally allows the calculation of the transit spectrum over an arbitrary wavelength range.

The assets of sunbather are:

• It has a highly detailed NLTE treatment of the ionization and excitation state and temperature of the atmosphere through the use of Cloudy. The calculation of a nonisothermal temperature profile should produce transit spectra that are more accurate than those provided by an isothermal Parker wind model.

• It allows an arbitrary composition that includes the 30 lightest elements. Metals are completely accounted for in both the hydrodynamics and thermodynamics. The transit spectrum includes metal lines, and enables joint modeling of multiple spectral lines.

• It is highly suited to the interpretation of spectral line observations through a parametrized description of the atmosphere. The mass-loss rate is a free parameter of the model, allowing one to retrieve it from observations.

• It has a reasonably fast runtime of a few minutes to an hour, depending on the composition and density of the atmosphere. This runtime supports parameter space exploration for individual planets and/or studies at the population level (but may require a computer cluster to do so).

The limitations of sunbather are:

• It is a 1D code, while atmospheric escape is an inherently 3D process. For example, day-to-night variations, magnetic fields, and interaction with a stellar wind can significantly alter the outflow and observational signatures (e.g., Nail et al., 2024; Schreyer et al., 2024).

• It keeps the composition of the atmosphere constant with radius. Every element escapes at the same rate, and it thus implicitly assumes efficient drag of heavier elements. This assumption may be invalid, particularly at low mass-loss rates.

• It does not treat the hydrodynamics and thermodynamics self-consistently, but instead parametrizes the outflow with a chosen mass-loss rate. As pointed out when describing the assets of the code above, it is this feature that makes sunbather well suited to retrievals. The downside is that it is less suited to forward modeling and theoretical studies designed to predict the mass-loss rate.

One may prefer to use sunbather over p-winds when fitting spectrally unresolved data and/or metal lines that are not included in p-winds. Conversely, one may prefer to use p-winds when fitting lines included in it over a large parameter space, due to its significantly faster runtime.

We use sunbather to explore the effect of atmospheric metallicity on the escape process. We simulate a generic hot-Neptune planet orbiting the Sun, with different atmospheric compositions ranging from hydrogen/helium-only to 50 times solar metallicity. If we keep the mass-loss rate fixed, a higher metallicity leads to a denser but slower wind, while the temperature profile starts to become isothermal. However, if we instead keep the density at the planet radius fixed, we can gain insight into the effect of metallicity on the mass-loss rate. We then find that a higher metallicity leads to a more tenuous and slower wind, with little effect on the temperature profile. The mass-loss rate decreases by an order of magnitude as we move from solar-like to 50 times solar metallicity. We furthermore find that metals significantly contribute to the radiative cooling rate already at solar abundances. They become important in the radiative heating rate at approximately ten times solar metallicity and above. Regarding the planet’s transit spectrum, all metal lines become stronger up to ten times solar metallicity. At even higher metallicity, lines of atoms and singly charged ions generally become weaker owing to the lower mass-loss rate, while lines of doubly and higher charged ions generally continue to grow stronger. We do not investigate whether all of the described effects of high metallicity also pertain to other types of star–planet systems.

After validating sunbather’s ability to retrieve a known model from a mock spectrum, we used the code to interpret the metastable helium spectrum of mini-Neptune TOI-2134 b as observed by Zhang et al. (2023a). If we assume a pure hydrogen/helium composition in the ratio of 10:1 by number, we find a mass-loss rate of $\dot{M}=10^{9.28^{+0.06}_{-0.04}}$ g s^-1 and a temperature of $T_{0}=5000\pm 400$ K. However, a model self-consistency check indicates a much cooler temperature of $T_{0}\sim 2600$ K, which is hard to reconcile with the data, unless other types of line broadening not accounted for here (e.g., atmospheric turbulence and circulation) play a significant role. If we instead assume a 100 times solar metallicity composition, we find a mass-loss rate of $\dot{M}=10^{9.68\pm 0.04}$ g s^-1 and a temperature of $T_{0}=9700\pm 700$ K. This temperature is significantly different from the self-consistent temperature of $T_{0}\sim 4000$ K, which again potentially indicates additional physics. In order to constrain the composition, instead of assuming it, additional observations of other spectral lines would be valuable. Such observations would also help to elucidate the interesting temperature differences that we find.