Issue 
A&A
Volume 659, March 2022



Article Number  A62  
Number of page(s)  12  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202142038  
Published online  08 March 2022 
pwinds: An opensource Python code to model planetary outflows and upper atmospheres^{⋆}
^{1}
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
email: ldsantos@stsci.edu
^{2}
Observatoire astronomique de l’Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland
^{3}
Leiden Observatory, Leiden University, Postbus 9513, 2300 RA Leiden, The Netherlands
^{4}
School of Physics, Trinity College Dublin, the University of Dublin, College Green, Dublin2, Ireland
^{5}
Division of Geological and Planetary Sciences, California Institute of Technology, 1200 East California Blvd, Pasadena, CA 91125, USA
^{6}
Center for Astrophysics  Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
^{7}
Earth and Planets Laboratory, The Carnegie Institution for Science, 5241 Broad Branch Road, Washington, DC 20015, USA
^{8}
Department of Physics, and Institute for Research on Exoplanets, Université de Montréal, Montréal H3T 1J4, Canada
^{9}
European Southern Observatory, Alonso de Córdova 3107, Santiago, Casilla 19001, Chile
Received:
17
August
2021
Accepted:
30
November
2021
Atmospheric escape is considered to be one of the main channels for evolution in subJovian planets, particularly in their early lives. While there are several hypotheses proposed to explain escape in exoplanets, testing them with atmospheric observations remains a challenge. In this context, highresolution transmission spectroscopy of transiting exoplanets for the metastable helium triplet (He 2^{3}S) at 1083 nm has emerged as a reliable technique for observing and measuring escape. To aid in the prediction and interpretation of metastable He transmission spectroscopy observations, we developed the code pwinds. This is an opensource, fully documented, scalable Python implementation of the onedimensional, purely H+He Parker wind model for upper atmospheres coupled with ionization balance, raytracing, and radiative transfer routines. We demonstrate an atmospheric retrieval by fitting pwinds models to the observed metastable He transmission spectrum of the warm Neptune HATP11 b and take the variation in the intransit absorption caused by transit geometry into account. For this planet, our best fit yields a total atmospheric escape rate of approximately 2.5 × 10^{10} g s^{−1} and an outflow temperature of 7200 K. The range of retrieved mass loss rates increases significantly when we let the H atom fraction be a free parameter, but its posterior distribution remains unconstrained by He observations alone. The stellar host limb darkening does not have a significant impact on the retrieved escape rate or outflow temperature for HATP11 b. Based on the nondetection of escaping He for GJ 436 b, we are able to rule out total escape rates higher than 3.4 × 10^{10} g s^{−1} at 99.7% (3σ) confidence.
Key words: methods: numerical / planets and satellites: atmospheres
The source code can be freely obtained in https://github.com/ladsantos/pwinds. Documentation, installation instructions, and tutorials are available in https://pwinds.readthedocs.io/. Contributions to the project are welcome.
© ESO 2022
1. Introduction
The evolution of shortperiod exoplanets is thought to be dictated by atmospheric escape. This conclusion is supported by two different approaches: (i) the detection of planetary outflows and large escape rates in hot exoplanets (e.g., VidalMadjar et al. 2003; Ehrenreich et al. 2015; Bourrier et al. 2018) and (ii) the observation of demographic features possibly carved by atmospheric escape in the population of Neptunes and superEarths (Beaugé & Nesvorný 2013; Owen & Wu 2013; Mazeh et al. 2016; Fulton et al. 2017; Fulton & Petigura 2018; HardegreeUllman et al. 2020). These discoveries have led the community to attempt to combine the theoretical descriptions of escape based on demographic features to predict observable atmospheric signatures in transiting exoplanets (e.g., Lecavelier Des Etangs 2007; Salz et al. 2016; King et al. 2019; Carolan et al. 2020). This experiment has been challenging, mainly because of limitations in our instruments and our theories (e.g., Cubillos et al. 2017; Kasper et al. 2020; Gaidos et al. 2020; Bean et al. 2021).
There are four known spectroscopic windows for observing atmospheric escape: the Lymanα line at 121.57 nm (VidalMadjar et al. 2003), metallic chromospheric lines and continuum in the ultraviolet (Fossati et al. 2010; Sing et al. 2019), the metastable helium triplet at 1083 nm (Seager & Sasselov 2000; Oklopčić & Hirata 2018), and the Balmer series of H lines in the blue optical (Jensen et al. 2012; Wyttenbach et al. 2020). Each one of them has its own set of challenges. While UV observations have classically been used to this end with a variable degree of success (e.g., Lecavelier des Etangs et al. 2010; Fossati et al. 2010; VidalMadjar et al. 2013; Waalkes et al. 2019; Dos Santos et al. 2020a, 2021; Bourrier et al. 2021; García Muñoz et al. 2021), they are particularly complicated because only the Hubble Space Telescope (HST) can access this wavelength range at high spectral resolution; in addition, cool stars usually do not have UV continuum, limiting transmission spectroscopy only to chromospheric or transitionregion emission lines whose count rates are very low (Bourrier et al. 2017; Dos Santos et al. 2019).
One of these techniques, He transmission spectroscopy, has been shown to be reliable and attainable using ground and spacebased instruments (Spake et al. 2018; Allart et al. 2018). This spectral channel is not photonstarved and is devoid of interstellar medium absorption (Indriolo et al. 2009), the main limitations of Lymanα spectroscopy. The disadvantage is that the formation of metastable He in the upper atmospheres of exoplanets depends on a specific level of irradiation arriving at the planet (e.g., Nortmann et al. 2018; Oklopčić 2019; Dos Santos et al. 2020b). Nevertheless, He spectroscopy has the potential to become the main technique of atmospheric escape observations (Allart et al. 2019; AlonsoFloriano et al. 2019; Kirk et al. 2020; Vissapragada et al. 2020; Lampón et al. 2021; Paragas et al. 2021).
Upper atmospheres extend to several planetary radii and can dwarf the size of planethosting stars depending on the properties of the system (e.g., Chamberlain 1963; Chaffin et al. 2015; Lavie et al. 2017; Kameda et al. 2017; Bourrier et al. 2018). For this reason, when observing the upper atmospheres of exoplanets, the transit geometry can have important effects on the interpretation of the data. For example, if a transiting planet has a nonzero impact parameter, a large portion of its exosphere may not transit and thus not contribute to the observed intransit absorption. Furthermore, a subtler effect in timeseries analyses of transmission spectroscopy is the dilution of a planetary absorption signal when the data are coadded in phase space. Since upper atmospheres are extended, the intransit absorption is variable with time. This variability dilutes the intransit absorption because time series of transmission spectra are coadded in phase space to improve the signaltonoise ratio of the combined transmission spectrum (e.g., Wyttenbach et al. 2015).
There are currently no publicly available tools to predict and interpret metastable He transmission spectroscopy. Considering that there is a broad community interest in these observations, we developed pwinds, an opensource, fully documented Python implementation of the onedimensional, isothermal Parker wind^{1} description (Parker 1958), to model exoplanet atmospheres. This code is timely because many data sets used to study metastable He spectroscopy have recently become public. Furthermore, an opensource implementation allows for an independent verification of results as well as community contributions to the code. pwinds implements limb darkening and a raytracing algorithm that allows the user to change the transit geometry (namely the transit impact parameter and phase in relation to midtransit).
In this manuscript we describe the overarching implementation of pwinds, discuss the design decisions, and illustrate the usage of the code. In Sect. 2 we describe the several modules implemented in the code to forward model the metastable He signature in a transiting exoplanet. In Sect. 3 we present case studies of the warm Neptunes HATP11 b and GJ 436 b and their corresponding atmospheric escape rates retrieved by fitting pwinds models to observations. Finally, in Sect. 4 we discuss the conclusions of this work.
2. Methods
The code pwinds is largely based on the formulations of Oklopčić & Hirata (2018) and Lampón et al. (2020). In its current version, the code has four core modules (and two support modules) to model the upper atmospheres and ionization balance of H and He around planetary bodies, which we describe below. In principle, these modules can be used independently of one another depending on the objective of the user. The code to reproduce the examples shown in this section can be obtained via the pwinds documentation.
2.1. The parker module
The parker module calculates the structure of the upper atmosphere following the theoretical description of the solar wind by Parker (1958). In this model, a steadystate, spherically symmetric outflow follows the equation of mass conservation:
where ṁ is the mass loss rate, r is the radius, ρ is the gas density, and v is the outflow velocity. This model also follows the steadystate momentum equation:
where G is the gravitational constant, p is the thermal pressure and M_{pl} is the planetary mass. The isothermal Parker solar wind model assumes that the outflow is completely ionized, yielding a constant mean molecular weight, μ, and consequently a constant sound speed, v_{s}, as a function of radial distance. This allows for a significant simplification of the problem. However, as pointed out by Lampón et al. (2020), the upper atmosphere of a hot planet differs from the solar wind in that μ(r) is not necessarily constant with radial distance. But if we assume that the ratio T(r)/μ(r) is constant over r, then the assumption of a constant sound speed profile still holds. If we let be the value of mean molecular weight corresponding to a given temperature T_{0}, then the constant sound speed is calculated as
According to Lampón et al. (2020), the temperature T_{0} in this approach corresponds to roughly the maximum of the temperature profile obtained by more comprehensive, selfconsistent models (see Sect. 3.1 in their manuscript). As we see in the following paragraphs, we arrive at a similar conclusion in our calculations as well.
We calculate μ(r) as
where m_{p} is the mass of a proton and n_{X} is the number density of species X. For clarity, we note that n_{H} = n_{H0} + n_{H+} and n_{He} = n_{He 11S} (singlet) + n_{He 23S} (triplet) + n_{He+}. We assume that the electrons coming from He ionization do not significantly contribute to changes in μ. According to Oklopčić & Hirata (2018), who also make this same assumption, including electrons from He ionization increases their number density by up to ∼10%.
For a given temperature T_{0}, the corresponding average mean molecular weight, , is calculated as in Eq. (A.3) of Lampón et al. (2020):
It is convenient to normalize the radii, velocities and densities to, respectively, the radius at the sonic point (r_{s}), the constant speed of sound, and the density at the sonic point (ρ_{s}). Based on the formulation of Lamers & Cassinelli (1999), the resulting equations describing the radial profiles of velocity and density are
where , , and are the normalized radial distance, velocity, and density, respectively.
Calculating the structure of the upper atmosphere requires as input: the planetary parameters and the stellar spectrum from Xrays to ultraviolet (XUV) impinging at the top of the atmosphere, as well as values for the atmospheric temperature and escape rate; the latter two are free parameters in the Parker wind model. Equation (6) is transcendental and requires a numerical approach to determine its solutions. To this end, we utilize a NewtonRaphson method implemented in scipy.optimize, which requires an initial guess for the optimization. Equation (6) has many solutions, but we are only interested in the solution that represents an escaping atmosphere (i.e., a transonic solution). In order to guarantee we converge to the correct solution, we enforce that the initial guess is below (above) the speed of sound when calculating the velocities below (above) the sonic point. The end product of the parker module is the structure of the upper atmosphere from Eqs. (6) and (7).
Since the structure of the upper atmosphere (densities and velocities) and the H ionization fraction are interdependent of each other, the code performs a loop that iteratively calculates all of these onedimensional profiles until convergence is achieved (see Sect. 2.2). As an example, we calculated the structure the hot Jupiter HD 209458 b and the result is shown in Fig. 1. We compare the structure computed with pwinds (continuous curves) with a onedimensional model of the same planet calculated selfconsistently using the formulation of Allan & Vidotto (2019 dashed curves, assuming 90% H and 10% He). In order to be comparable, the pwinds model was computed using the same mass loss rate, composition, and a monochromatic ionizing flux (energy > 13.6 eV) of 450 erg s^{−1} cm^{−2} with an average photon energy of 20 eV. The isothermal Parker wind predicts a similar structure as the selfconsistent model when we assume a T_{0} corresponding to the maximum temperature from the latter. This is the same result that Lampón et al. (2020) obtains in their description.
Fig. 1. Onedimensional structure of the upper atmosphere of the hot Jupiter HD 209458 b computed with pwinds (continuous curves). Velocities are shown in blue and densities in orange. For comparison, we also plot a model for the same planet computed selfconsistently using the formulation of Allan & Vidotto (2019) as dashed curves. The circles mark the sonic point. 
Naturally, a onedimensional model does not capture outflow asymmetries that are sometimes observed in Lymanα transit spectroscopy (e.g., Lecavelier des Etangs et al. 2010; Lavie et al. 2017; Bourrier et al. 2018). More complex, threedimensional models are necessary to completely describe these features (e.g., Bourrier et al. 2016; Villarreal D’Angelo et al. 2021; Wang & Dai 2021a,b; Allart et al. 2019; MacLeod & Oklopčić 2021). Simple onedimensional models are nevertheless capable of retrieving atmospheric escape parameters with the assumption that the mass loss process takes place spherically and homogeneously throughout the surface of the planet (e.g., Lampón et al. 2020, 2021). Models that are faster to calculate are also useful when there is a need to explore a large parameter space, which is what we discuss in Sect. 3 and in an upcoming manuscript (Vissapragada et al. 2022).
2.2. The hydrogen module
The hydrogen module calculates the steadystate distribution of neutral and ionized H in the upper atmosphere. The quantity of interest here is f_{ion}, whose radial profile is obtained by calculating the steadystate balance between advection and sourcesink terms for H ions. In this case, the source is (photo) ionization by highenergy photons, and the sink is recombination into neutral atoms. This radial distribution can be calculated with the following differential equation (see Sect. 3.2 in Oklopčić & Hirata 2018):
where n_{H}(r) = x ρ(r)/[(x + 4y) m_{p}], with x being the H atoms number fraction in the outflow, y = 1 − x is the He atoms number fraction, and Φ is the photoionization rate:
where λ_{0} is the wavelength corresponding to the ionization energy of H (911.65 Å) and f_{λ} is the flux density (in units of energy ⋅ time^{−1}⋅ area^{−1}⋅ wavelength^{−1}) arriving at the top of the atmosphere. σ_{λ} is the photoionization cross section, which we calculate in the support module microphysics, following Eq. (10) in Oklopčić & Hirata (2018), which is based on Osterbrock & Ferland (2006). The optical depth of neutral H is given by
The velocities v and densities ρ are calculated using the module parker. α_{rec} is the caseB H recombination rate at a given temperature (Osterbrock & Ferland 2006; Tripathi et al. 2015), calculated as
As seen in Eq. (10), τ_{λ} depends on f_{ion}, which is what we want to calculate in the first place. However, the optical depth depends more strongly on the densities of H than the ion fraction. So instead of solving a system of coupled nonlinear differential equations, a first solution can be achieved by assuming that the whole atmosphere is neutral at first. Later, we relax this assumption by recalculating the τ_{λ} and f_{ion} profiles iteratively until convergence is achieved (the user can define the convergence criterion).
We solve Eq. (8) using solve_ivp, an explicit RungeKutta integrator of hybrid 4th and 5th orders implemented in scipy.integrate. The user inputs an initial guess for f_{ion} at the innermost layer of the upper atmosphere. The code also takes as input the stellar host spectrum arriving at the planet, or the monochromatic flux between 0 and 911.65 Å, and the planetary parameters. The solution for the H distribution in 500 points including the relaxation takes approximately 400 ms on a CPU with frequency 3.1 GHz and four computing threads. Continuing the example for HD 209458 b from Sect. 2.1, we calculated the ion and neutral fraction of H in the upper atmosphere, and the resulting distribution is shown in Fig. 2 (continuous curve). We compare this result with the ion fraction calculated with the selfconsistent escape model from Sect. 2.1 (dashed curve); in order to be comparable, both models are calculated assuming an impinging XUV monochromatic flux of 450 erg s^{−1} cm^{−2}. The pwinds model overpredicts the ion fraction by a factor of a few when compared to the selfconsistent model, likely because of the larger densities (see Fig. 1), which increase the optical depth of the atmosphere to ionizing irradiation.
Fig. 2. Neutral H atom fraction in the upper atmosphere of the hot Jupiter HD 209458 b computed with pwinds for the same setup from Sect. 2.1 (continuous curve). We also show the neutral fraction calculated with a selfconsistent escape model for comparison (dashed curve). 
2.3. The helium module
The helium module calculates the steadystate distribution of neutral singlet, neutral triplet, and ionized He in the upper atmosphere. The quantities of interest here are f_{1} = n_{He 11S}/n_{He} and f_{3} = n_{He 23S}/n_{He}. The radial profiles df_{1}/dr and ddf_{3}/dr are described by a coupled system of differential equations with source and sink terms:
We refer the reader to Oklopčić & Hirata (2018) and Table 2 of Lampón et al. (2020) for detailed equations of all the source and sink terms for He^{2}. In our code we do include the He charge exchange terms pointed out by Lampón et al. (2020).
We assume that the He ionization and the excited helium triplet do not significantly change the structure of the upper atmosphere. This allows us to decouple the helium module from the parker and hydrogen modules. This is advantageous because the user can enter as input a H structure that was calculated by more complex and selfconsistent models than isothermal Parker winds ones. It is important, however, that these models do include He in their calculation of the structure in order to produce consistent results for the metastable He distribution.
The procedure to solve the distribution of He (Eq. (12)) is similar to that for H. The user inputs an initial guess for f_{1} and f_{3} at the innermost atmospheric layer, the stellar host spectrum from 0 to 2600 Å (or monochromatic fluxes in the bands 0–1200 Å and 1200–2600 Å), the structure of the upper atmosphere (profiles of density and velocity), and the planetary parameters. It is important to emphasize that neutral H also contributes to the optical depth between wavelengths 0–911 Å, attenuating the amount of highenergy flux that ionizes and populates the He levels. The code takes this contribution into account, as in Oklopčić & Hirata (2018).
Initially, the code assumes that the entire upper atmosphere has constant f_{1} and f_{3}, and then a first solution is obtained using odeint^{3}, a Python wrapper for the LSODA solver from the Fortran library odepack implemented in scipy.integrate. This solution is then relaxed by updating the optical depths, f_{1} and f_{3} iteratively until convergence is achieved. The solution can, however, become numerically unstable for large density gradients, which can sometimes happen near the R = 1 R_{pl}. A practical workaround is to establish a cutoff near 1 R_{pl} that removes this large density gradient and ignore this layer of the atmosphere in the modeling. However, the user should be aware that this solution could affect the interpretation of more compressed thermospheres, such as that of HD 189733 b (Lampón et al. 2021); we have not yet attempted to model this planet with pwinds, and leave this for future work.
We show the distribution of He in the upper atmosphere of HD 209458 b calculated as described above in Fig. 3. For comparison purposes, this time we assumed a model with the same input parameters as the one described in Oklopčić & Hirata (2018 namely an escape rate of 8 × 10^{10} g s^{−1}, a temperature of 9000 K, a H fraction of 0.9, and a solar irradiating spectrum). Our results match the models of Oklopčić & Hirata, as seen in Fig. 3 of their publication. The solution for the He distribution in 500 points including the relaxation takes approximately 2.5 s on a CPU with frequency 3.1 GHz and four computing threads. This is the main computational bottleneck of the pwinds code.
Fig. 3. Distribution of He in the upper atmosphere of HD 209458 b calculated with pwinds assuming the same input parameters as Oklopčić & Hirata (2018). 
2.4. The transit module
The transit module has two independent functions that can be used to calculate the spectral signatures of the upper atmosphere in transmission. The first function, draw_transit, calculates twodimensional intensity maps containing the host star and a transiting planet at a userdefined phase and impact parameter. The onedimensional profiles of metastable He volumetric densities are required to calculate the twodimensional array of column densities mapped to the same geometry as the transit. The output intensity map is normalized in a way that the diskaveraged stellar intensity is 1.0 when the planet is out of transit. Optionally, the user can also input a limbdarkening law.
The most important function in this module is radiative_transfer_2d, which, as the name implies, calculates the intransit absorption spectrum caused by the opaque disk of the planet and its upper atmosphere. In each cell of the twodimensional transit mapped by the ij indexes, the resulting attenuated intensity I_{ij}(ν)^{4} of the stellar light caused by absorption of He in the upper atmosphere is given by
where I_{ij, 0}(ν) is the intensity emerging from the host star before filtering through the atmosphere and τ_{ij, He} is the optical depth due to metastable He. I_{ij}(ν) is set to zero in the cells corresponding to the opaque disk of the planet. From here onward, we drop the ij indexes for the sake of brevity, but the reader should implicitly assume that the radiative transfer is carried out cellbycell in the transit map. Formally, the optical depth is given by
where φ_{ν} is the Voigt profile and σ_{He} is the cross section of metastable helium lines near 1.083 μm. Following Oklopčić & Hirata (2018) (see also, e.g., Allan & Vidotto 2019), the He cross section is calculated as
where f is the oscillator strength of the transition, e is the electron charge, and m_{e} is the electron mass. This formula is only valid in the Gaussiancgs unit system, where e is given in units of esu or statC (see, for example, Koskinen et al. 2010 for a formula that can be used in other unit systems).
The Voigt profile φ_{ν} is calculated using the voigt_profile implementation of scipy.special, which takes three parameters: the bulk velocity shift v_{bulk} of the profile in relation to the rest wavelength, the standard deviation α of the Gaussian (in our case Doppler) term, and the Lorentzian half width at half maximum (HWHM). Similar to Lampón et al. (2020), the Gaussian width α is calculated as
where m_{He} is the mass of a He atom, T is the temperature of the gas, ν_{0} is the central frequency of the transition. The Lorentzian HWHM is γ = A_{ij}/4π, where A_{ij} is the Einstein coefficient of the transition. We took the properties of the metastable He transitions near 1.083 μm from the National Institute of Standards and Technology (NIST) database^{5}, and list them in Table 1.
Spectral line properties of the metastable He triplet in the nearinfrared.
2.4.1. Line broadening by the planetary outflow
In reality, φ_{ν} depends on the threedimensional position in relation to the planet because each position has a different lineofsight velocity, which broadens the absorption line. Thus, the formal calculation of the Voigt profile is performed for each pencil of light between the star and the observer. In a given position z along the pencil, the lineofsight velocity, v_{LOS}, as a function of distance, r, from the planet center is calculated using the formulation of Seidel et al. (2020):
where v_{ver} is the outflow velocity obtained from the Parker wind model.
This calculation has to be performed for three spectral lines, and it adds an extra dimension for wavelength. For these reasons, the formal calculation of φ_{ν} taking into account all the four dimensions is computationally costly. In order to accelerate the radiative transfer, instead of calculating the Parker wind broadening in full dimensionality, we can optionally assume that it contributes to the Gaussian broadening term of the Voigt profile uniformly through the line of sight. With the dependence on the z axis dropped, we can remove φ_{ν} from the integrand in Eq. (14), yielding the approximation
where η_{He} is the column density of He. In order to validate this approximation, we need to assume that the Gaussian wind broadening has a constant velocity v_{w} in the line of sight, and add it in quadrature to the squarevelocity term of Eq. (16), yielding
The windbroadening velocity term v_{w} is calculated as the average of v_{LOS} weighted by the metastable He number density:
In this approximation, the user can, optionally, include an additional source of broadening: microturbulence. We implement the same formulation of Lampón et al. (2020):
The turbulence velocity term is added quadratically in Eq. (19).
pwinds allows the user to decide on which method to use to calculate the Voigt profile: the formal calculation (Eqs. (13) and (14)) or the densityweighted average broadening parameter (Eqs. (18)–(20)). The turbulent broadening (Eq. (21)) can optionally be included at the discretion of the user. This is done with the optional parameters wind_broadening_method and turbulence_broadening when calling the radiative_transfer_2d function; the default is the densityweighted averagevelocity implementation, which is a good compromise of speed and accuracy.
We assess the validity of the assumption we made above by calculating the formal and the averagevelocity broadening methods for HD 209458 b and HATP11 b. The first planet has a more compact upper atmosphere and lower outflow velocities than the second by a factor of ∼2. In the case of HD 209458 b, the averagevelocity method produces an accurate approximation to the formal calculation (see the left panel of Fig. 4). In the case of the more extended atmosphere of HATP11 b, the averagevelocity method is accurate when turbulence broadening is also included (right panel of Fig. 4). In both cases, the average method is one order of magnitude faster in computation time than the formal method.
Fig. 4. Comparison of the spectral line broadening in the He triplet lines using two different methods: formal definition of the optical depth (black) and the densityweighted averagevelocity broadening (red). In the case of a more extended atmosphere of HATP11 b, a better match between the formal and average methods is obtained when we include turbulence broadening (orange). These are forward models, and we are not yet attempting to fit them to observed signatures. 
We emphasize that, at this point, we are only producing a forward model and not making an attempt to fit it to actual existing observations of this planet) (e.g., Lampón et al. 2020). The results we obtain from the example of HD 209458 b throughout this section, from the Parker wind structure to the predicted metastable He transmission spectrum, are consistent with those obtained by Oklopčić & Hirata (2018).
2.4.2. Dilution of the transit signature
Usually, the absorption of light by upper atmospheres in transiting exoplanets is in the order of several percent or less in a narrow bandpass. Thus, transmission spectroscopy observations sometimes rely on averaging time series in phasespace to build enough signal to noise and produce a detectable signal. However, upper atmospheres are so extended that the intransit absorption signature is variable with phase, and phaseaveraging them dilutes the observed signature. In addition, inhomogeneities in the stellar surface, such as limb darkening, may become important, particularly when fitting (spectro) photometric light curves.
We illustrate this effect in Fig. 5, where we simulated the phaseaveraging for both HD 209458 b and HATP11 b. In the case of the hotJupiter, the midtransit spectrum (phase = 0.0) is comparable to the spectrum phaseaveraged between second and third contacts (T2–T3), but it differs more significantly when the phaseaveraging is taken between first and fourth contacts (T1–T4). This is due to a more compact upper atmosphere than HATP11 b, whose extent produces a larger difference between the phaseaveraged and the midtransit spectra. Since the planettostar ratio is smaller than HD 209458 b, phase averaging between T1–T4 or T2–T3 does not make as much difference as it does for the hot Jupiter.
Fig. 5. Metastable He transmission spectrum of HD 209458 b (left panel) and HATP11 b (right panel) for uniformly sampled phases, transit impact parameter b = 0.499 and b = 0.132, respectively, and limb darkening based on the results of Knutson et al. (2007) and Sing (2010). The baseline (R_{p}/R_{s})^{2} was removed, as in actual groundbased observations. Phase 0.0 and 0.5 represent, respectively, the times of midtransit and of first (or fourth) contact. The dashed red spectrum is the average of all phases between the first and fourth contact. The dotdashed black curve is the average between second and third contact. These are forward models, and we are not yet attempting to fit them to observed signatures. 
Previous onedimensional descriptions of the metastable He transmission spectrum did not take into account the transit geometry, phaseaveraging, and limb darkening. The pwinds code allows the user to set the transit impact parameter, phase in relation to the first and fourth contacts, and set a limbdarkening law. To this end, we utilize the auxiliary opensource code flatstar^{6} to simulate transit grids (see a brief description in Appendix A). In the current implementation, this transit grid only allows for circular orbits and it neglects the curvature of the transit chord.
3. Atmospheric escape retrievals
We further benchmarked the code pwinds by performing retrievals of the atmospheric escape rate, temperature, lineofsight bulk velocity, and the H fraction (n_{H}/n_{atoms}) of the warm Neptunes HATP11 b and GJ 436 b. Both planets were observed in transmission spectroscopy by the CARMENES spectrograph (Calar Alto highResolution search for M dwarfs with Exoearths with Nearinfrared and optical Échelle Spectrographs), but only the first had a strong intransit signal and the second had only a nondetection. For the latter, we attempted to fit upper or lower limits of the atmospheric escape rate and outflow temperature. The Python algorithms to reproduce our retrievals are freely available online^{7}.
3.1. Fitting the He signature of HATP11 b
The metastable He signature of HATP11 b was measured with the CARMENES spectrograph installed in the 3.5 m telescope at the Calar Alto Observatory (Allart et al. 2018). The transmission spectrum is openly available in the DACE platform^{8}. The central wavelengths of the metastable He transitions retrieved from the NIST database are listed as measured in air, but the wavelengths of the CARMENES spectrum are in vacuum. We converted the wavelengths of the latter to inair using the following formula (Morton 2000, IAU standard):
In general, we fit three free parameters: the atmospheric escape rate ṁ, the upper atmosphere temperature T, and the bulk lineofsight velocity v_{bulk} of the upper atmosphere. It is also possible to run fits with additional parameters (such as the H fraction). The fit is performed by maximizing the likelihood 𝒫 of a given transmission spectrum model ℱ_{model} to represent the observed transmission spectrum ℱ. Such a loglikelihood is given by
where n stands for a given bin of the spectrum, σ is the uncertainty of the measurement, and p is the vector containing the free parameters. To determine the uncertainties of the fit, we use the Markov chain Monte Carlo (MCMC) ensemble sampler emcee (ForemanMackey et al. 2013). The uncertainties we report here represent the confidence interval that encompasses the 16th to the 84th percentile of the posterior distribution of the free parameters.
For HATP11 b, we ran in total four different models: (1) no limb darkening, H number fraction fixed to 0.90; (2) a quadratic limbdarkening law with coefficients c_{1} = 0.63 and c_{2} = 0.09 (Sing 2010), H number fraction fixed to 0.90; (3) no limb darkening, H fraction as a free parameter with a uniform prior of [0.80, 0.99]; and (4) same as (1), but using the formal implementation of the radiative transfer instead of the averagevelocity broadening approximation. For model (4), instead of a full MCMC, we perform only a maximumlikelihood (Eq. (23)) using the NelderMead algorithm implemented in scipy.optimize.minimize. The reason is because we simply want to assess the accuracy of the averagevelocity broadening approximation for HATP11 b in comparison to the formal, more computationally costly radiative transfer. HATP11 does not have a full highenergy spectrum measurement, so we use the spectrum of a similar star from the MUSCLES Treasury Survey^{9} (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016) as a proxy. We chose the star HD 40307, which has similar effective temperature, mass, radius, and surface gravity as HATP11.
We ran the MCMC for 7000 steps and 10 walkers in 10 cores of a computer cluster with an average frequency of 3.0 GHz per core. The autocorrelation time t of the MCMC was, on average, 45 steps when we started from a first guess of 1 × 10^{10} g s^{−1}, 800 K, and −2.0 km s^{−1}, respectively, for ṁ, T, and v_{bulk}. We remove a total of 2t burnin steps from the beginning of the MCMC and take a sample thinned by t/2, resulting in a flat chain of approximately 3600 samples. The computation of the MCMC chains took approximately 6.5 h of computing time. Different planets will likely yield different computing times because the numerical bottleneck (calculating the He distribution) is highly dependent on the input parameters. We show the posterior distributions of the fit parameters for Model 1 in Fig. 6 (see also Appendix B), and a sample of corresponding transmission spectrum models fit to the data in Fig. 7.
Fig. 6. Posterior distributions of the mass loss rate, upper atmospheric temperature, and lineofsight bulk velocity of HATP11 b using pwinds models (no limb darkening included) as a retrieval tool against a CARMENES transmission spectrum (Model 1). 
Fig. 7. Transmission spectrum of HATP11 b measured with CARMENES (black) and a sample of 100 pwinds models (Model 1) fit to the data (red). 
Table 2 contains the retrieved atmospheric escape parameters for HATP11 b based on the CARMENES transmission spectrum. All models we tested yield results consistent with one another within their uncertainties, based on the marginalized posterior distribution of the retrieved parameters. A comparison between the results for Models 1 and 4 reveals that, in the case of HATP11 b, the limb darkening of the star does not significantly affect the retrieved escape parameters when fitting a groundbased transmission spectrum. In the case of Model 3, where we allowed the H fraction to vary between 0.80 and 0.99, the retrieval slightly favors fractions > 0.96, but the 3σ upper limit of > 0.80 is not constraining. We show the resulting posterior distributions of the fit to Model 3 in Fig. 8. For fractions above 0.92, the retrieved escape rate tends to increase by a factor of several percent. The retrieved upper atmosphere temperature T is highly anticorrelated with the H fraction. This degeneracy increases the uncertainties of T by a factor of at least two. We show the resulting distribution of He in the upper atmosphere of HATP11 b in Fig. 9 based on the bestfit model to the CARMENES data. Finally, we show that the retrieved escape parameters of Models 13 (averagevelocity approximation for the wind broadening) are fully consistent with that of Model 4 (formal radiative transfer calculation). Hence, we demonstrate that this approximation, which saves an order of magnitude in computation time, does not significantly affect the retrieved escape parameters.
Fig. 9. Distribution of He in the upper atmosphere of HATP11 b based on the bestfit solution obtained by fitting pwinds models (no limb darkening included) to the CARMENES transmission spectrum (Model 1). 
Upper atmosphere properties retrieved for HATP11 b from the CARMENES transmission spectrum.
In order to compare our results with those obtained by the threedimensional model EVE (Bourrier & Lecavelier des Etangs 2013) used by Allart et al. (2018), we need to calculate the escape of metastable He only instead of total mass loss. For a total escape rate of 2.5 × 10^{10} g s^{−1} and an upper atmosphere temperature of 7200 K obtained from the retrieval described above, we calculate an average metastable helium fraction of 4.8 × 10^{−6} and a T/μ fraction of 8000 K amu^{−1}. This result translates into a metastableHe escape rate of 1.2 × 10^{5} g s^{−1}, which is compatible with the upperlimit rate of ∼3 × 10^{5} g s^{−1} determined by Allart et al. (2018). Our retrieved T/μ fraction when assuming a H fraction of 0.9 is discrepant with the results of Allart et al. (2018), who find T/μ = 24 000 K amu^{−1}. Some of the solutions of our retrieval with the H fraction as a free parameter do allow for high values of T/μ up to 12 000, but they are nevertheless incompatible with Allart et al. (2018); the authors, however, do propose that a high T/μ may correspond to a low mean atomic weight, which can be obtained with a large fraction of ionized gas and free electrons. In an upcoming manuscript, Vissapragada et al. (2022) will discuss how solutions with high temperatures can be ruled out because they are not energetically selfconsistent, assuming that the heating comes solely from the available highenergy irradiation budget. A possible explanation for the disagreement can be due to modeling difference, since Allart et al. (2018) use a hydrostatic model while we use a hydrodynamic one; since a hydrostatic thermosphere is less extended, it needs a higher temperature to increase the density of He enough to be detectable at high altitudes.
The bulk velocity of −1.9 ± 0.8 km s^{−1} is consistent with the net blueshift of 3 km s^{−1} reported by Allart et al. (2018), which was previously interpreted as a highaltitude wind flowing from the day to the nightside of the planet. This net blueshift is not predicted by the onedimensional Parker wind model, which is the reason for fitting it as a free parameter in our models. More complex, tridimensional models that take into account other physical processes may be necessary to determine the exact mechanism that causes this bulk velocity shift in the metastable He absorption signature.
3.2. Escape rate upper limit for a nondetection in GJ 436 b
GJ 436 b is a highprofile case of atmospheric escape because it possesses the deepest transmission spectrum feature detected to date: a repeatable 50% intransit absorption in Lymanα (Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017; Dos Santos et al. 2019), which is explained by a large volume of exospheric neutral H fed by escape (Bourrier et al. 2016; Villarreal D’Angelo et al. 2021). In fact, Oklopčić & Hirata (2018) predict a metastable He signature as deep as 9% in the core of the strongest line of the triplet. However, when GJ 436 b was observed by CARMENES, the results yielded only a nondetection (Nortmann et al. 2018).
In this section we attempt to measure an upper limit of atmospheric escape rate to the nondetection of He in GJ 436 b and compare it to the result derived from Lymanα transmission spectroscopy and modeling. The metastable He transmission spectrum of GJ 436 b is rather unfortunately not publicly available, but the pipelinereduced spectral time series from Nortmann et al. (2018) is available in the CARMENES data archive^{10}.
We ran an MCMC of 10 000 steps and 10 walkers with three free parameters: total escape rate, upper atmospheric temperature, and the H fraction. We increased the number of steps compared to the HATP11 b retrieval in order to better explore the parameter space, since we are expecting to obtain only upper/lower limits. We did not include limb darkening. Based on previous theoretical predictions for GJ 436 b (e.g., Salz et al. 2016), we set uniform priors of [10^{7}, 10^{12}] g s^{−1} for the mass loss and [1000, 10 000] K for the temperature, and [0.40, 0.99] for the H fraction. We used the highenergy spectrum of GJ 436 measured in the MUSCLES Treasury Survey as a source of irradiation.
The resulting posterior distributions of the free parameters for GJ 436 b yield, at 99.7% (3σ) confidence, an upper limit of 4.5 × 10^{9} g s^{−1} for the escape rate and a lower limit of 2600 K for the upper atmospheric temperature (given the uniform priors above). In broad terms, mass loss rates above this value or temperatures below the lower limit would yield a detectable metastable He signature. With the H fraction as a free parameter, these 3σ limits become 3.4 × 10^{10} g s^{−1} and 2400 K. This result is consistent with the escape rate of ∼2.5 × 10^{8} g s^{−1} inferred by Bourrier et al. (2016), and with the mass loss rate of (6 − 10)×10^{9} g s^{−1} inferred by Villarreal D’Angelo et al. (2021), both based on the same Lymanα transmission spectroscopy data set.
Given the flat prior of [0.40, 0.99] for the H fraction, the resulting posterior distribution of this parameter is not constraining; however, it seems to favor higher values and peaks near 0.99 (see Fig. B.3). Interestingly, this result could be seen as an agreement with the prediction of a Hrich outflow for GJ 436 b due to selective escape, leading to a Herich lower atmosphere (Hu et al. 2015). We cannot, however, draw strong conclusions on this matter because GJ 436 b had a nondetection of He. More detailed descriptions that fit both H and He simultaneously in the upper atmosphere of this planet are likely going to yield more definitive answers. For example, Lampón et al. (2021) used H densities derived from Lymanα observations to inform metastable He models, and determine that the warm Neptune GJ 3470 b has .
4. Conclusions
We demonstrate in this manuscript the usage of the opensource Python code pwinds to forward model the distribution of He atoms in the upper atmospheres of exoplanets as well as their corresponding metastable He transmission spectra of exoplanets. The code also enables the retrieval of atmospheric escape rates and temperatures based on observations at high resolution when coupled to an optimization algorithm, such as a maximum likelihood estimation or an MCMC sampler. A typical retrieval takes several hours to compute, depending on the setup.
As an implementation of the method originally described by Oklopčić & Hirata (2018), the forward models produced by pwinds are fully compatible with that study. We also implement changes proposed by Lampón et al. (2020), such as the inclusion of charge exchange of He and H particles. Our implementation includes further improvements, such as the addition of transit geometry and the limb darkening of the host star, as well as allowing the H fraction (n_{H}/n_{atoms}) as a free parameter in the retrieval.
We used pwinds to fit the escape rate, outflow temperature, and the H fraction of the warm Neptune HATP11 b based on CARMENES transmission spectroscopy previously reported in Allart et al. (2018). For a model without limb darkening and with n_{H}/n_{atoms} fixed at 0.90, we find that the escape rate of HATP11 b is g s^{−1} and the planetary outflow temperature is 7200 ± 700 K. This temperature is in disagreement with the value of T/μ calculated by Allart et al. (2018), and it is likely caused by a key difference between our models – theirs contains a hydrostatic thermosphere, while ours is hydrodynamic. Including limb darkening does not have a significant impact on the retrieved parameters of HATP11 b. Allowing the H fraction as a free parameter has a stronger impact because it yields an anticorrelation with the retrieved outflow temperature. It also increases the uncertainty of the retrieved atmospheric escape rate. We find that the H fraction is unconstrained, but with a preference for higher values. These results are in agreement with those of Lampón et al. (2020, 2021), although those authors can constrain the H fraction by analyzing He transmission spectra in conjunction with H escape using Lymanα observations.
Finally, we also attempted to fit limits for the escape rate, outflow temperature, and H fraction of GJ 436 b based on a nondetection with the CARMENES spectrograph reported by Nortmann et al. (2018). We find an upper limit of 3.4 × 10^{10} g s^{−1} for the first and a lower limit of 2400 K for the second at 99.7% confidence. Our upper and lower limit determinations show a preference for high values of n_{H}/n_{atoms}, with the posterior distribution peaking near 0.99. These results are fully compatible with the escape rate of ∼2.5 × 10^{8} g s^{−1} inferred by Bourrier et al. (2016) based on Lymanα transmission spectroscopy. For both HATP11 b and GJ 436 b, we find a slight preference for high values of H in the atomic fraction, which is in line with the results of Lampón et al. (2020, 2021) for other hot gas giants.
The main limitations of a onedimensional, isothermal Parker wind model are: (1) It does not capture the threedimensional nature of very extended atmospheres, particularly when they have both a thermospheric and an exospheric contributions (see the case of WASP107 b in Allart et al. 2019; Spake et al. 2021); (2) it does not take into account the variable profile of temperature with radial distance from the planet, which is seen in selfconsistent models of escape (e.g., Salz et al. 2016; Allan & Vidotto 2019) and (3) it does not selfconsistently consider the sources of heating and cooling that control the atmospheric escape process. The usefulness of simple models such as pwinds lies in an efficient exploration of the parameter space that defines atmospheric escape (scalability) and ease of use (opensource, fully documented code) when more sophisticated models are not yet warranted.
As for the next steps, we aim to improve pwinds by including the escape of heavier atomic species, such as C, N, O, Mg, Si, and Fe. This will allow us to use the code to predict and interpret observations of metals escaping hot gas giants, such as the signatures reported by VidalMadjar et al. (2013) and Sing et al. (2019). We shall also add daytonightside winds to the atmospheric modeling, similar to Seidel et al. (2020). Another avenue to explore pwinds in the future consists in coupling it with more complex tridimensional hydrodynamic escape models.
We note that, in Table 2 of Lampón et al. (2020), the units for the recombination and collisional processes are cm^{3} s^{−1}, and not cm^{−3} s^{−1} as the authors list in their manuscript.
The code is freely available at https://flatstar.readthedocs.io.
Available at https://archive.stsci.edu/prepds/muscles/.
Acknowledgments
L.A.D.S. acknowledges the helpful input of A. Wyttenbach, M. Stalport, A. Oklopčić, J. Stürmer, and M. Zechmeister to the development of this project. The authors also thank the referee, Manuel LópezPuertas, for the helpful and detailed review. S.V. is supported by an NSF Graduate Research Fellowship and the Paul & Daisy Soros Fellowship for New Americans. R.A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation, and his contribution was supported in part through a grant from Fonds de recherche du Québec – Nature et technologies. This research was enabled by the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (projects: FOUR ACES grant agreement No 724427; SPICE DUNE grant agreement No 947634; ASTROFLOW grant agreement No 817540), and it has been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). The pwinds code makes use of the open source software NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pillow (https://pythonpillow.org), and Astropy (Astropy Collaboration 2018). The results of this manuscript were also made possible by the open source software Matplotlib (Hunter 2007), OpenMPI (https://www.openmpi.org), Jupyter (Kluyver et al. 2016), MPI for Python (mpi4py; Dalcin et al. 2011), emcee (ForemanMackey et al. 2013), and schwimmbad (PriceWhelan & ForemanMackey 2017). Finally, the authors also extend a special thanks to the platforms GitHub, CondaForge, Read the Docs, and Travis.ci for the valuable support of opensource initiatives.
References
 Allan, A., & Vidotto, A. A. 2019, MNRAS, 490, 3760 [Google Scholar]
 Allart, R., Bourrier, V., Lovis, C., et al. 2018, Science, 362, 1384 [Google Scholar]
 Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AlonsoFloriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019, A&A, 629, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Bean, J. L., Raymond, S. N., & Owen, J. E. 2021, J. Geophys. Res. (Planets), 126 [Google Scholar]
 Beaugé, C., & Nesvorný, D. 2013, ApJ, 763, 12 [Google Scholar]
 Bourrier, V., & Lecavelier des Etangs, A. 2013, A&A, 557, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017, A&A, 599, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., et al. 2018, A&A, 620, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., dos Santos, L. A., SanzForcada, J., et al. 2021, A&A, 650, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carolan, S., Vidotto, A. A., Plavchan, P., Villarreal D’Angelo, C., & Hazra, G. 2020, MNRAS, 498, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Chaffin, M. S., Chaufray, J. Y., Deighan, J., et al. 2015, Geophys. Res. Lett., 42, 9001 [NASA ADS] [CrossRef] [Google Scholar]
 Chamberlain, J. W. 1963, Planet. Space Sci., 11, 901 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A. 2000, A&A, 363, 1081 [NASA ADS] [Google Scholar]
 Claret, A., & Hauschildt, P. H. 2003, A&A, 412, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017, MNRAS, 466, 1868 [Google Scholar]
 Dalcin, L. D., Paz, R. R., Kler, P. A., & Cosimo, A. 2011, in Advances in Water Resources, 34, 1124, New Computational Methods and Software Tools [NASA ADS] [CrossRef] [Google Scholar]
 DiazCordoves, J., & Gimenez, A. 1992, A&A, 259, 227 [Google Scholar]
 Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2019, A&A, 629, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2020a, A&A, 634, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2020b, A&A, 640, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dos Santos, L. A., Bourrier, V., Ehrenreich, D., et al. 2021, A&A, 649, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Fossati, L., Haswell, C. A., Froning, C. S., et al. 2010, ApJ, 714, L222 [Google Scholar]
 France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
 Gaidos, E., Hirano, T., Wilson, D. J., et al. 2020, MNRAS, 498, L119 [NASA ADS] [CrossRef] [Google Scholar]
 García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36 [Google Scholar]
 HardegreeUllman, K. K., Zink, J. K., Christiansen, J. L., et al. 2020, ApJS, 247, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Hu, R., Seager, S., & Yung, Y. L. 2015, ApJ, 807, 8 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Indriolo, N., Hobbs, L. M., Hinkle, K. H., & McCall, B. J. 2009, ApJ, 703, 2131 [NASA ADS] [CrossRef] [Google Scholar]
 Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86 [Google Scholar]
 Kameda, S., Ikezawa, S., Sato, M., et al. 2017, Geophys. Res. Lett., 44, 11706 [Google Scholar]
 Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, AJ, 160, 258 [Google Scholar]
 King, G. W., Wheatley, P. J., Bourrier, V., & Ehrenreich, D. 2019, MNRAS, 484, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Kirk, J., Alam, M. K., LópezMorales, M., & Zeng, L. 2020, AJ, 159, 115 [Google Scholar]
 Klinglesmith, D. A., & Sobieski, S. 1970, AJ, 75, 175 [Google Scholar]
 Kluyver, T., RaganKelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Scmidt (Netherlands: IOS Press), 87 [Google Scholar]
 Knutson, H. A., Charbonneau, D., Noyes, R. W., Brown, T. M., & Gilliland, R. L. 2007, ApJ, 655, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Kopal, Z. 1950, Harvard College Obs. Circular, 454, 1 [NASA ADS] [Google Scholar]
 Koskinen, T. T., Yelle, R. V., Lavvas, P., & Lewis, N. K. 2010, ApJ, 723, 116 [CrossRef] [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
 Kulow, J. R., France, K., Linsky, J., & Loyd, R. O. P. 2014, ApJ, 786, 132 [Google Scholar]
 Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds [Google Scholar]
 Lampón, M., LópezPuertas, M., Lara, L. M., et al. 2020, A&A, 636, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lampón, M., LópezPuertas, M., SanzForcada, J., et al. 2021, A&A, 647, A129 [EDP Sciences] [Google Scholar]
 Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lecavelier Des Etangs, A. 2007, A&A, 461, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lecavelier des Etangs, A., Ehrenreich, D., VidalMadjar, A., et al. 2010, A&A, 514, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Loyd, R. O. P., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, M., & Oklopčić, A. 2021, AAS J., submitted [arXiv:2107.07534] [Google Scholar]
 Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morton, D. C. 2000, ApJS, 130, 403 [Google Scholar]
 Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388 [Google Scholar]
 Oklopčić, A. 2019, ApJ, 881, 133 [Google Scholar]
 Oklopčić, A., & Hirata, C. M. 2018, ApJ, 855, L11 [Google Scholar]
 Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei [Google Scholar]
 Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
 Paragas, K., Vissapragada, S., Knutson, H. A., et al. 2021, ApJ, 909, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
 PriceWhelan, A. M., & ForemanMackey, D. 2017, J. Open Source Softw., 2, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916 [Google Scholar]
 Seidel, J. V., Ehrenreich, D., Pino, L., et al. 2020, A&A, 633, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sing, D. K. 2010, A&A, 510, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sing, D. K., Désert, J. M., Lecavelier Des Etangs, A., et al. 2009, A&A, 505, 891 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91 [Google Scholar]
 Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68 [Google Scholar]
 Spake, J. J., Oklopčić, A., & Hillenbrand, L. A. 2021, AJ, 162, 284 [NASA ADS] [CrossRef] [Google Scholar]
 Tripathi, A., Kratter, K. M., MurrayClay, R. A., & Krumholz, M. R. 2015, ApJ, 808, 173 [NASA ADS] [CrossRef] [Google Scholar]
 VidalMadjar, A., Lecavelier des Etangs, A., Désert, J.M., et al. 2003, Nature, 422, 143 [Google Scholar]
 VidalMadjar, A., Huitson, C. M., Bourrier, V., et al. 2013, A&A, 560, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Villarreal D’Angelo, C., Vidotto, A. A., Esquivel, A., Hazra, G., & Youngblood, A. 2021, MNRAS, 501, 4383 [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
 Vissapragada, S., Knutson, H. A., Jovanovic, N., et al. 2020, AJ, 159, 278 [Google Scholar]
 Vissapragada, S., Knutson, H. A., Dos Santos, L. A., Wang, L., & Dai, F. 2022, ApJ, submitted [arXiv:2201.09889] [Google Scholar]
 Waalkes, W. C., BertaThompson, Z., Bourrier, V., et al. 2019, AJ, 158, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., & Dai, F. 2021a, ApJ, 914, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., & Dai, F. 2021b, ApJ, 914, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Wyttenbach, A., Ehrenreich, D., Lovis, C., Udry, S., & Pepe, F. 2015, A&A, 577, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyttenbach, A., Mollière, P., Ehrenreich, D., et al. 2020, A&A, 638, A87 [EDP Sciences] [Google Scholar]
 Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101 [Google Scholar]
Appendix A: The flatstar code
The code implemented in flatstar was originally written as a part of the transit module of pwinds. However, we decided to transform it into a separate package because this implementation can be useful for other astrophysical applications not necessarily related to transmission spectroscopy.
The typical usage of flatstar involves setting a grid shape (N_{x}, N_{y}) and stellar radius R_{s} in number of pixels, and optionally setting a limbdarkening law. The limbdarkening laws currently implemented in the code are: linear, quadratic (Kopal 1950), squareroot (DiazCordoves & Gimenez 1992), logarithmic (Klinglesmith & Sobieski 1970), exponential (Claret & Hauschildt 2003), the threeparameter law of Sing et al. (2009), and the fourparameter law of Claret (2000). Finally, the user can also set a custom limbdarkening law. The star is always centered to the grid.
In addition, flatstar can add a planetary transit with userdefined planettostar ratio, R_{p}/R_{s}, transit impact parameter, b, and phase, ϕ. The first and fourth contact of the transit are defined as the phases −0.5 and +0.5, respectively, independent of b. The y coordinate of the planetary center (y_{p}) in pixel space is calculated as
The x coordinate of the planetary center (x_{p}) is not as straightforward to calculate, since we define it based on ϕ and b. Let θ be the angle between the y axis and the vector that connects the center of the star and the planet at first contact:
The distance β from the planet center to the y axis at first contact is given by
Thus, the distance x_{0} of the planet center from the border of the simulation at first contact is
As the planet moves from the first contact to fourth contact, it covers a distance of 2β. For arbitrary phases between ϕ = −0.5 (first contact) and ϕ = +0.5 (fourth contact), the distance x_{p} of the planet from the border of the simulation is thus
This formulation assumes that the arc that the planet follows during the transit can be approximated to a chord.
The grid can be supersampled in order to avoid "hard" pixel edges when the grid size is coarse. This is useful to save computation time in cases where we need to mass produce grids while conserving the precision of intensities (which is the case of atmospheric retrievals with pwinds). By default, the resampling algorithm is the "box" method, which takes the value of each pixel with fixed weights to compute the average flux of the resampled pixel. We do not recommend using flatstar to fit wide band photometric light curves, since the computation time is much longer than other codes that implement analytical equations to calculate light curves, such as batman (Kreidberg 2015).
Appendix B: Other posterior distributions for the fits to HATP11 b and GJ 436 b data
Fig. B.1. Posterior distribution of parameters of HATP11 b fit to the CARMENES transmission spectrum using a model with a quadratic limbdarkening law and coefficients c_{1} = 0.63 and c_{2} = 0.09. 
Fig. B.2. Same as Fig. B.1, but for GJ 436 b with the H fraction fixed to 0.90 and no limb darkening. 
All Tables
Upper atmosphere properties retrieved for HATP11 b from the CARMENES transmission spectrum.
All Figures
Fig. 1. Onedimensional structure of the upper atmosphere of the hot Jupiter HD 209458 b computed with pwinds (continuous curves). Velocities are shown in blue and densities in orange. For comparison, we also plot a model for the same planet computed selfconsistently using the formulation of Allan & Vidotto (2019) as dashed curves. The circles mark the sonic point. 

In the text 
Fig. 2. Neutral H atom fraction in the upper atmosphere of the hot Jupiter HD 209458 b computed with pwinds for the same setup from Sect. 2.1 (continuous curve). We also show the neutral fraction calculated with a selfconsistent escape model for comparison (dashed curve). 

In the text 
Fig. 3. Distribution of He in the upper atmosphere of HD 209458 b calculated with pwinds assuming the same input parameters as Oklopčić & Hirata (2018). 

In the text 
Fig. 4. Comparison of the spectral line broadening in the He triplet lines using two different methods: formal definition of the optical depth (black) and the densityweighted averagevelocity broadening (red). In the case of a more extended atmosphere of HATP11 b, a better match between the formal and average methods is obtained when we include turbulence broadening (orange). These are forward models, and we are not yet attempting to fit them to observed signatures. 

In the text 
Fig. 5. Metastable He transmission spectrum of HD 209458 b (left panel) and HATP11 b (right panel) for uniformly sampled phases, transit impact parameter b = 0.499 and b = 0.132, respectively, and limb darkening based on the results of Knutson et al. (2007) and Sing (2010). The baseline (R_{p}/R_{s})^{2} was removed, as in actual groundbased observations. Phase 0.0 and 0.5 represent, respectively, the times of midtransit and of first (or fourth) contact. The dashed red spectrum is the average of all phases between the first and fourth contact. The dotdashed black curve is the average between second and third contact. These are forward models, and we are not yet attempting to fit them to observed signatures. 

In the text 
Fig. 6. Posterior distributions of the mass loss rate, upper atmospheric temperature, and lineofsight bulk velocity of HATP11 b using pwinds models (no limb darkening included) as a retrieval tool against a CARMENES transmission spectrum (Model 1). 

In the text 
Fig. 7. Transmission spectrum of HATP11 b measured with CARMENES (black) and a sample of 100 pwinds models (Model 1) fit to the data (red). 

In the text 
Fig. 8. Same as Fig. 6, but including the H fraction as a free parameter to be fit (Model 3). 

In the text 
Fig. 9. Distribution of He in the upper atmosphere of HATP11 b based on the bestfit solution obtained by fitting pwinds models (no limb darkening included) to the CARMENES transmission spectrum (Model 1). 

In the text 
Fig. B.1. Posterior distribution of parameters of HATP11 b fit to the CARMENES transmission spectrum using a model with a quadratic limbdarkening law and coefficients c_{1} = 0.63 and c_{2} = 0.09. 

In the text 
Fig. B.2. Same as Fig. B.1, but for GJ 436 b with the H fraction fixed to 0.90 and no limb darkening. 

In the text 
Fig. B.3. Same as Fig. B.2, but with the H fraction as a free parameter. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.