Free Access
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/0004-6361/202142038
Published online 08 March 2022

© ESO 2022

1. Introduction

The evolution of short-period 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., Vidal-Madjar 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 super-Earths (Beaugé & Nesvorný 2013; Owen & Wu 2013; Mazeh et al. 2016; Fulton et al. 2017; Fulton & Petigura 2018; Hardegree-Ullman 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 (Vidal-Madjar 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; Vidal-Madjar 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 transition-region 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 space-based instruments (Spake et al. 2018; Allart et al. 2018). This spectral channel is not photon-starved 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; Alonso-Floriano 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 planet-hosting 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 in-transit absorption. Furthermore, a subtler effect in time-series analyses of transmission spectroscopy is the dilution of a planetary absorption signal when the data are co-added in phase space. Since upper atmospheres are extended, the in-transit absorption is variable with time. This variability dilutes the in-transit absorption because time series of transmission spectra are co-added in phase space to improve the signal-to-noise 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 p-winds, an open-source, fully documented Python implementation of the one-dimensional, isothermal Parker wind1 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 open-source implementation allows for an independent verification of results as well as community contributions to the code. p-winds implements limb darkening and a ray-tracing algorithm that allows the user to change the transit geometry (namely the transit impact parameter and phase in relation to mid-transit).

In this manuscript we describe the overarching implementation of p-winds, 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 HAT-P-11 b and GJ 436 b and their corresponding atmospheric escape rates retrieved by fitting p-winds models to observations. Finally, in Sect. 4 we discuss the conclusions of this work.

2. Methods

The code p-winds 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 p-winds 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 steady-state, spherically symmetric outflow follows the equation of mass conservation:

(1)

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 steady-state momentum equation:

(2)

where G is the gravitational constant, p is the thermal pressure and Mpl 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, vs, 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 T0, then the constant sound speed is calculated as

(3)

According to Lampón et al. (2020), the temperature T0 in this approach corresponds to roughly the maximum of the temperature profile obtained by more comprehensive, self-consistent 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

(4)

where mp is the mass of a proton and nX is the number density of species X. For clarity, we note that nH = nH0 + nH+ and nHe = nHe 11S (singlet) + nHe 23S (triplet) + nHe+. 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 T0, the corresponding average mean molecular weight, , is calculated as in Eq. (A.3) of Lampón et al. (2020):

(5)

It is convenient to normalize the radii, velocities and densities to, respectively, the radius at the sonic point (rs), 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

(6)

(7)

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 X-rays 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 Newton-Raphson 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 one-dimensional 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 p-winds (continuous curves) with a one-dimensional model of the same planet calculated self-consistently using the formulation of Allan & Vidotto (2019 dashed curves, assuming 90% H and 10% He). In order to be comparable, the p-winds 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 self-consistent model when we assume a T0 corresponding to the maximum temperature from the latter. This is the same result that Lampón et al. (2020) obtains in their description.

thumbnail Fig. 1.

One-dimensional structure of the upper atmosphere of the hot Jupiter HD 209458 b computed with p-winds (continuous curves). Velocities are shown in blue and densities in orange. For comparison, we also plot a model for the same planet computed self-consistently using the formulation of Allan & Vidotto (2019) as dashed curves. The circles mark the sonic point.

Naturally, a one-dimensional 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, three-dimensional 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 one-dimensional 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 steady-state distribution of neutral and ionized H in the upper atmosphere. The quantity of interest here is fion, whose radial profile is obtained by calculating the steady-state balance between advection and source-sink terms for H ions. In this case, the source is (photo-) ionization by high-energy 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):

(8)

where nH(r) = xρ(r)/[(x + 4y) mp], 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:

(9)

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

(10)

The velocities v and densities ρ are calculated using the module parker. αrec is the case-B H recombination rate at a given temperature (Osterbrock & Ferland 2006; Tripathi et al. 2015), calculated as

(11)

As seen in Eq. (10), τλ depends on fion, 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 fion profiles iteratively until convergence is achieved (the user can define the convergence criterion).

We solve Eq. (8) using solve_ivp, an explicit Runge-Kutta integrator of hybrid 4th and 5th orders implemented in scipy.integrate. The user inputs an initial guess for fion 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 self-consistent 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 p-winds model overpredicts the ion fraction by a factor of a few when compared to the self-consistent model, likely because of the larger densities (see Fig. 1), which increase the optical depth of the atmosphere to ionizing irradiation.

thumbnail Fig. 2.

Neutral H atom fraction in the upper atmosphere of the hot Jupiter HD 209458 b computed with p-winds for the same setup from Sect. 2.1 (continuous curve). We also show the neutral fraction calculated with a self-consistent escape model for comparison (dashed curve).

2.3. The helium module

The helium module calculates the steady-state distribution of neutral singlet, neutral triplet, and ionized He in the upper atmosphere. The quantities of interest here are f1 = nHe 11S/nHe and f3 = nHe 23S/nHe. The radial profiles df1/dr and ddf3/dr are described by a coupled system of differential equations with source and sink terms:

(12)

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 He2. 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 self-consistent 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 f1 and f3 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 high-energy 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 f1 and f3, and then a first solution is obtained using odeint3, 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, f1 and f3 iteratively until convergence is achieved. The solution can, however, become numerically unstable for large density gradients, which can sometimes happen near the R = 1 Rpl. A practical work-around is to establish a cutoff near 1 Rpl 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 p-winds, 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 × 1010 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 p-winds code.

thumbnail Fig. 3.

Distribution of He in the upper atmosphere of HD 209458 b calculated with p-winds 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 two-dimensional intensity maps containing the host star and a transiting planet at a user-defined phase and impact parameter. The one-dimensional profiles of metastable He volumetric densities are required to calculate the two-dimensional array of column densities mapped to the same geometry as the transit. The output intensity map is normalized in a way that the disk-averaged stellar intensity is 1.0 when the planet is out of transit. Optionally, the user can also input a limb-darkening law.

The most important function in this module is radiative_transfer_2d, which, as the name implies, calculates the in-transit absorption spectrum caused by the opaque disk of the planet and its upper atmosphere. In each cell of the two-dimensional transit mapped by the ij indexes, the resulting attenuated intensity Iij(ν)4 of the stellar light caused by absorption of He in the upper atmosphere is given by

(13)

where Iij,  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. Iij(ν) 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 cell-by-cell in the transit map. Formally, the optical depth is given by

(14)

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

(15)

where f is the oscillator strength of the transition, e is the electron charge, and me is the electron mass. This formula is only valid in the Gaussian-cgs 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 vbulk 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

(16)

where mHe 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 γ = Aij/4π, where Aij 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) database5, and list them in Table 1.

Table 1.

Spectral line properties of the metastable He triplet in the near-infrared.

2.4.1. Line broadening by the planetary outflow

In reality, φν depends on the three-dimensional position in relation to the planet because each position has a different line-of-sight 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 line-of-sight velocity, vLOS, as a function of distance, r, from the planet center is calculated using the formulation of Seidel et al. (2020):

(17)

where vver 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

(18)

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 vw in the line of sight, and add it in quadrature to the square-velocity term of Eq. (16), yielding

(19)

The wind-broadening velocity term vw is calculated as the average of vLOS weighted by the metastable He number density:

(20)

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):

(21)

The turbulence velocity term is added quadratically in Eq. (19).

p-winds allows the user to decide on which method to use to calculate the Voigt profile: the formal calculation (Eqs. (13) and (14)) or the density-weighted 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 density-weighted average-velocity 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 average-velocity broadening methods for HD 209458 b and HAT-P-11 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 average-velocity 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 HAT-P-11 b, the average-velocity 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.

thumbnail 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 density-weighted average-velocity broadening (red). In the case of a more extended atmosphere of HAT-P-11 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 phase-space to build enough signal to noise and produce a detectable signal. However, upper atmospheres are so extended that the in-transit absorption signature is variable with phase, and phase-averaging 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 phase-averaging for both HD 209458 b and HAT-P-11 b. In the case of the hot-Jupiter, the mid-transit spectrum (phase = 0.0) is comparable to the spectrum phase-averaged between second and third contacts (T2–T3), but it differs more significantly when the phase-averaging is taken between first and fourth contacts (T1–T4). This is due to a more compact upper atmosphere than HAT-P-11 b, whose extent produces a larger difference between the phase-averaged and the mid-transit spectra. Since the planet-to-star 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.

thumbnail Fig. 5.

Metastable He transmission spectrum of HD 209458 b (left panel) and HAT-P-11 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 (Rp/Rs)2 was removed, as in actual ground-based observations. Phase 0.0 and 0.5 represent, respectively, the times of mid-transit and of first (or fourth) contact. The dashed red spectrum is the average of all phases between the first and fourth contact. The dot-dashed 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 one-dimensional descriptions of the metastable He transmission spectrum did not take into account the transit geometry, phase-averaging, and limb darkening. The p-winds code allows the user to set the transit impact parameter, phase in relation to the first and fourth contacts, and set a limb-darkening law. To this end, we utilize the auxiliary open-source code flatstar6 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 p-winds by performing retrievals of the atmospheric escape rate, temperature, line-of-sight bulk velocity, and the H fraction (nH/natoms) of the warm Neptunes HAT-P-11 b and GJ 436 b. Both planets were observed in transmission spectroscopy by the CARMENES spectrograph (Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs), but only the first had a strong in-transit signal and the second had only a non-detection. 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 online7.

3.1. Fitting the He signature of HAT-P-11 b

The metastable He signature of HAT-P-11 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 platform8. 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 in-air using the following formula (Morton 2000, IAU standard):

(22)

In general, we fit three free parameters: the atmospheric escape rate , the upper atmosphere temperature T, and the bulk line-of-sight velocity vbulk 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 log-likelihood is given by

(23)

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 (Foreman-Mackey 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 HAT-P-11 b, we ran in total four different models: (1) no limb darkening, H number fraction fixed to 0.90; (2) a quadratic limb-darkening law with coefficients c1 = 0.63 and c2 = 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 average-velocity broadening approximation. For model (4), instead of a full MCMC, we perform only a maximum-likelihood (Eq. (23)) using the Nelder-Mead algorithm implemented in scipy.optimize.minimize. The reason is because we simply want to assess the accuracy of the average-velocity broadening approximation for HAT-P-11 b in comparison to the formal, more computationally costly radiative transfer. HAT-P-11 does not have a full high-energy spectrum measurement, so we use the spectrum of a similar star from the MUSCLES Treasury Survey9 (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 HAT-P-11.

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 × 1010 g s−1, 800 K, and −2.0 km s−1, respectively, for , T, and vbulk. We remove a total of 2t burn-in 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.

thumbnail Fig. 6.

Posterior distributions of the mass loss rate, upper atmospheric temperature, and line-of-sight bulk velocity of HAT-P-11 b using p-winds models (no limb darkening included) as a retrieval tool against a CARMENES transmission spectrum (Model 1).

thumbnail Fig. 7.

Transmission spectrum of HAT-P-11 b measured with CARMENES (black) and a sample of 100 p-winds models (Model 1) fit to the data (red).

Table 2 contains the retrieved atmospheric escape parameters for HAT-P-11 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 HAT-P-11 b, the limb darkening of the star does not significantly affect the retrieved escape parameters when fitting a ground-based 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 HAT-P-11 b in Fig. 9 based on the best-fit model to the CARMENES data. Finally, we show that the retrieved escape parameters of Models 1-3 (average-velocity 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.

thumbnail Fig. 8.

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

thumbnail Fig. 9.

Distribution of He in the upper atmosphere of HAT-P-11 b based on the best-fit solution obtained by fitting p-winds models (no limb darkening included) to the CARMENES transmission spectrum (Model 1).

Table 2.

Upper atmosphere properties retrieved for HAT-P-11 b from the CARMENES transmission spectrum.

In order to compare our results with those obtained by the three-dimensional 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 × 1010 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 metastable-He escape rate of 1.2 × 105 g s−1, which is compatible with the upper-limit rate of ∼3 × 105 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 self-consistent, assuming that the heating comes solely from the available high-energy 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 high-altitude wind flowing from the day- to the night-side of the planet. This net blueshift is not predicted by the one-dimensional 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 non-detection in GJ 436 b

GJ 436 b is a high-profile case of atmospheric escape because it possesses the deepest transmission spectrum feature detected to date: a repeatable 50% in-transit 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 non-detection (Nortmann et al. 2018).

In this section we attempt to measure an upper limit of atmospheric escape rate to the non-detection 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 pipeline-reduced spectral time series from Nortmann et al. (2018) is available in the CARMENES data archive10.

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 HAT-P-11 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 [107, 1012] 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 high-energy 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 × 109 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 × 1010 g s−1 and 2400 K. This result is consistent with the escape rate of ∼2.5 × 108 g s−1 inferred by Bourrier et al. (2016), and with the mass loss rate of (6 − 10)×109 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 H-rich outflow for GJ 436 b due to selective escape, leading to a He-rich lower atmosphere (Hu et al. 2015). We cannot, however, draw strong conclusions on this matter because GJ 436 b had a non-detection 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 open-source Python code p-winds 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 p-winds 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 (nH/natoms) as a free parameter in the retrieval.

We used p-winds to fit the escape rate, outflow temperature, and the H fraction of the warm Neptune HAT-P-11 b based on CARMENES transmission spectroscopy previously reported in Allart et al. (2018). For a model without limb darkening and with nH/natoms fixed at 0.90, we find that the escape rate of HAT-P-11 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 HAT-P-11 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 non-detection with the CARMENES spectrograph reported by Nortmann et al. (2018). We find an upper limit of 3.4 × 1010 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 nH/natoms, with the posterior distribution peaking near 0.99. These results are fully compatible with the escape rate of ∼2.5 × 108 g s−1 inferred by Bourrier et al. (2016) based on Lyman-α transmission spectroscopy. For both HAT-P-11 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 one-dimensional, isothermal Parker wind model are: (1) It does not capture the three-dimensional nature of very extended atmospheres, particularly when they have both a thermospheric and an exospheric contributions (see the case of WASP-107 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 self-consistent models of escape (e.g., Salz et al. 2016; Allan & Vidotto 2019) and (3) it does not self-consistently consider the sources of heating and cooling that control the atmospheric escape process. The usefulness of simple models such as p-winds lies in an efficient exploration of the parameter space that defines atmospheric escape (scalability) and ease of use (open-source, fully documented code) when more sophisticated models are not yet warranted.

As for the next steps, we aim to improve p-winds 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 Vidal-Madjar et al. (2013) and Sing et al. (2019). We shall also add day-to-nightside winds to the atmospheric modeling, similar to Seidel et al. (2020). Another avenue to explore p-winds in the future consists in coupling it with more complex tridimensional hydrodynamic escape models.


1

In this manuscript, we use the terms “wind” and “outflow” interchangeably, but the first should not be confused with horizontal winds in the lower atmosphere.

2

We note that, in Table 2 of Lampón et al. (2020), the units for the recombination and collisional processes are cm3 s−1, and not cm−3 s−1 as the authors list in their manuscript.

3

When calculating the steady-state distribution of He, we opted to use odeint instead of solve_ivp in this case because the first is more stable and 2.6 times faster than the second. solve_ivp is faster than odeint when solving the distribution of H.

4

The radiative transfer routine uses input in wavelength space, but the actual calculations are performed in frequency space for code clarity and brevity. The grid size is defined by the user.

6

The code is freely available at https://flatstar.readthedocs.io.

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ópez-Puertas, 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 p-winds code makes use of the open source software NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Pillow (https://python-pillow.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.open-mpi.org), Jupyter (Kluyver et al. 2016), MPI for Python (mpi4py; Dalcin et al. 2011), emcee (Foreman-Mackey et al. 2013), and schwimmbad (Price-Whelan & Foreman-Mackey 2017). Finally, the authors also extend a special thanks to the platforms GitHub, Conda-Forge, Read the Docs, and Travis.ci for the valuable support of open-source initiatives.

References

  1. Allan, A., & Vidotto, A. A. 2019, MNRAS, 490, 3760 [Google Scholar]
  2. Allart, R., Bourrier, V., Lovis, C., et al. 2018, Science, 362, 1384 [Google Scholar]
  3. Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Alonso-Floriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019, A&A, 629, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  6. Bean, J. L., Raymond, S. N., & Owen, J. E. 2021, J. Geophys. Res. (Planets), 126 [Google Scholar]
  7. Beaugé, C., & Nesvorný, D. 2013, ApJ, 763, 12 [Google Scholar]
  8. Bourrier, V., & Lecavelier des Etangs, A. 2013, A&A, 557, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. 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]
  10. Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017, A&A, 599, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., et al. 2018, A&A, 620, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bourrier, V., dos Santos, L. A., Sanz-Forcada, J., et al. 2021, A&A, 650, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Carolan, S., Vidotto, A. A., Plavchan, P., Villarreal D’Angelo, C., & Hazra, G. 2020, MNRAS, 498, L53 [NASA ADS] [CrossRef] [Google Scholar]
  14. Chaffin, M. S., Chaufray, J. Y., Deighan, J., et al. 2015, Geophys. Res. Lett., 42, 9001 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chamberlain, J. W. 1963, Planet. Space Sci., 11, 901 [NASA ADS] [CrossRef] [Google Scholar]
  16. Claret, A. 2000, A&A, 363, 1081 [NASA ADS] [Google Scholar]
  17. Claret, A., & Hauschildt, P. H. 2003, A&A, 412, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017, MNRAS, 466, 1868 [Google Scholar]
  19. 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]
  20. Diaz-Cordoves, J., & Gimenez, A. 1992, A&A, 259, 227 [Google Scholar]
  21. Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2019, A&A, 629, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2020a, A&A, 634, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2020b, A&A, 640, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Dos Santos, L. A., Bourrier, V., Ehrenreich, D., et al. 2021, A&A, 649, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459 [Google Scholar]
  26. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  27. Fossati, L., Haswell, C. A., Froning, C. S., et al. 2010, ApJ, 714, L222 [Google Scholar]
  28. France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89 [NASA ADS] [CrossRef] [Google Scholar]
  29. Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [Google Scholar]
  30. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
  31. Gaidos, E., Hirano, T., Wilson, D. J., et al. 2020, MNRAS, 498, L119 [NASA ADS] [CrossRef] [Google Scholar]
  32. García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36 [Google Scholar]
  33. Hardegree-Ullman, K. K., Zink, J. K., Christiansen, J. L., et al. 2020, ApJS, 247, 28 [NASA ADS] [CrossRef] [Google Scholar]
  34. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  35. Hu, R., Seager, S., & Yung, Y. L. 2015, ApJ, 807, 8 [Google Scholar]
  36. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  37. Indriolo, N., Hobbs, L. M., Hinkle, K. H., & McCall, B. J. 2009, ApJ, 703, 2131 [NASA ADS] [CrossRef] [Google Scholar]
  38. Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86 [Google Scholar]
  39. Kameda, S., Ikezawa, S., Sato, M., et al. 2017, Geophys. Res. Lett., 44, 11706 [Google Scholar]
  40. Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, AJ, 160, 258 [Google Scholar]
  41. King, G. W., Wheatley, P. J., Bourrier, V., & Ehrenreich, D. 2019, MNRAS, 484, L49 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kirk, J., Alam, M. K., López-Morales, M., & Zeng, L. 2020, AJ, 159, 115 [Google Scholar]
  43. Klinglesmith, D. A., & Sobieski, S. 1970, AJ, 75, 175 [Google Scholar]
  44. Kluyver, T., Ragan-Kelley, 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]
  45. Knutson, H. A., Charbonneau, D., Noyes, R. W., Brown, T. M., & Gilliland, R. L. 2007, ApJ, 655, 564 [NASA ADS] [CrossRef] [Google Scholar]
  46. Kopal, Z. 1950, Harvard College Obs. Circular, 454, 1 [NASA ADS] [Google Scholar]
  47. Koskinen, T. T., Yelle, R. V., Lavvas, P., & Lewis, N. K. 2010, ApJ, 723, 116 [CrossRef] [Google Scholar]
  48. Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
  49. Kulow, J. R., France, K., Linsky, J., & Loyd, R. O. P. 2014, ApJ, 786, 132 [Google Scholar]
  50. Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds [Google Scholar]
  51. Lampón, M., López-Puertas, M., Lara, L. M., et al. 2020, A&A, 636, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lampón, M., López-Puertas, M., Sanz-Forcada, J., et al. 2021, A&A, 647, A129 [EDP Sciences] [Google Scholar]
  53. Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Lecavelier Des Etangs, A. 2007, A&A, 461, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Lecavelier des Etangs, A., Ehrenreich, D., Vidal-Madjar, A., et al. 2010, A&A, 514, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Loyd, R. O. P., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102 [NASA ADS] [CrossRef] [Google Scholar]
  57. MacLeod, M., & Oklopčić, A. 2021, AAS J., submitted [arXiv:2107.07534] [Google Scholar]
  58. Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Morton, D. C. 2000, ApJS, 130, 403 [Google Scholar]
  60. Nortmann, L., Pallé, E., Salz, M., et al. 2018, Science, 362, 1388 [Google Scholar]
  61. Oklopčić, A. 2019, ApJ, 881, 133 [Google Scholar]
  62. Oklopčić, A., & Hirata, C. M. 2018, ApJ, 855, L11 [Google Scholar]
  63. Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei [Google Scholar]
  64. Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
  65. Paragas, K., Vissapragada, S., Knutson, H. A., et al. 2021, ApJ, 909, L10 [NASA ADS] [CrossRef] [Google Scholar]
  66. Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
  67. Price-Whelan, A. M., & Foreman-Mackey, D. 2017, J. Open Source Softw., 2, 388 [NASA ADS] [CrossRef] [Google Scholar]
  68. Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916 [Google Scholar]
  70. Seidel, J. V., Ehrenreich, D., Pino, L., et al. 2020, A&A, 633, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Sing, D. K. 2010, A&A, 510, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. 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]
  73. Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91 [Google Scholar]
  74. Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68 [Google Scholar]
  75. Spake, J. J., Oklopčić, A., & Hillenbrand, L. A. 2021, AJ, 162, 284 [NASA ADS] [CrossRef] [Google Scholar]
  76. Tripathi, A., Kratter, K. M., Murray-Clay, R. A., & Krumholz, M. R. 2015, ApJ, 808, 173 [NASA ADS] [CrossRef] [Google Scholar]
  77. Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J.-M., et al. 2003, Nature, 422, 143 [Google Scholar]
  78. Vidal-Madjar, A., Huitson, C. M., Bourrier, V., et al. 2013, A&A, 560, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Villarreal D’Angelo, C., Vidotto, A. A., Esquivel, A., Hazra, G., & Youngblood, A. 2021, MNRAS, 501, 4383 [CrossRef] [Google Scholar]
  80. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
  81. Vissapragada, S., Knutson, H. A., Jovanovic, N., et al. 2020, AJ, 159, 278 [Google Scholar]
  82. Vissapragada, S., Knutson, H. A., Dos Santos, L. A., Wang, L., & Dai, F. 2022, ApJ, submitted [arXiv:2201.09889] [Google Scholar]
  83. Waalkes, W. C., Berta-Thompson, Z., Bourrier, V., et al. 2019, AJ, 158, 50 [NASA ADS] [CrossRef] [Google Scholar]
  84. Wang, L., & Dai, F. 2021a, ApJ, 914, 98 [NASA ADS] [CrossRef] [Google Scholar]
  85. Wang, L., & Dai, F. 2021b, ApJ, 914, 99 [NASA ADS] [CrossRef] [Google Scholar]
  86. Wyttenbach, A., Ehrenreich, D., Lovis, C., Udry, S., & Pepe, F. 2015, A&A, 577, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Wyttenbach, A., Mollière, P., Ehrenreich, D., et al. 2020, A&A, 638, A87 [EDP Sciences] [Google Scholar]
  88. 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 p-winds. 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 (Nx, Ny) and stellar radius Rs in number of pixels, and optionally setting a limb-darkening law. The limb-darkening laws currently implemented in the code are: linear, quadratic (Kopal 1950), square-root (Diaz-Cordoves & Gimenez 1992), logarithmic (Klinglesmith & Sobieski 1970), exponential (Claret & Hauschildt 2003), the three-parameter law of Sing et al. (2009), and the four-parameter law of Claret (2000). Finally, the user can also set a custom limb-darkening law. The star is always centered to the grid.

In addition, flatstar can add a planetary transit with user-defined planet-to-star ratio, Rp/Rs, 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 (yp) in pixel space is calculated as

(A.1)

The x coordinate of the planetary center (xp) 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:

(A.2)

The distance β from the planet center to the y axis at first contact is given by

(A.3)

Thus, the distance x0 of the planet center from the border of the simulation at first contact is

(A.4)

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 xp of the planet from the border of the simulation is thus

(A.5)

This formulation assumes that the arc that the planet follows during the transit can be approximated to a chord.

The grid can be super-sampled 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 p-winds). 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 HAT-P-11 b and GJ 436 b data

thumbnail Fig. B.1.

Posterior distribution of parameters of HAT-P-11 b fit to the CARMENES transmission spectrum using a model with a quadratic limb-darkening law and coefficients c1 = 0.63 and c2 = 0.09.

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

thumbnail Fig. B.3.

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

All Tables

Table 1.

Spectral line properties of the metastable He triplet in the near-infrared.

Table 2.

Upper atmosphere properties retrieved for HAT-P-11 b from the CARMENES transmission spectrum.

All Figures

thumbnail Fig. 1.

One-dimensional structure of the upper atmosphere of the hot Jupiter HD 209458 b computed with p-winds (continuous curves). Velocities are shown in blue and densities in orange. For comparison, we also plot a model for the same planet computed self-consistently using the formulation of Allan & Vidotto (2019) as dashed curves. The circles mark the sonic point.

In the text
thumbnail Fig. 2.

Neutral H atom fraction in the upper atmosphere of the hot Jupiter HD 209458 b computed with p-winds for the same setup from Sect. 2.1 (continuous curve). We also show the neutral fraction calculated with a self-consistent escape model for comparison (dashed curve).

In the text
thumbnail Fig. 3.

Distribution of He in the upper atmosphere of HD 209458 b calculated with p-winds assuming the same input parameters as Oklopčić & Hirata (2018).

In the text
thumbnail 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 density-weighted average-velocity broadening (red). In the case of a more extended atmosphere of HAT-P-11 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
thumbnail Fig. 5.

Metastable He transmission spectrum of HD 209458 b (left panel) and HAT-P-11 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 (Rp/Rs)2 was removed, as in actual ground-based observations. Phase 0.0 and 0.5 represent, respectively, the times of mid-transit and of first (or fourth) contact. The dashed red spectrum is the average of all phases between the first and fourth contact. The dot-dashed 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
thumbnail Fig. 6.

Posterior distributions of the mass loss rate, upper atmospheric temperature, and line-of-sight bulk velocity of HAT-P-11 b using p-winds models (no limb darkening included) as a retrieval tool against a CARMENES transmission spectrum (Model 1).

In the text
thumbnail Fig. 7.

Transmission spectrum of HAT-P-11 b measured with CARMENES (black) and a sample of 100 p-winds models (Model 1) fit to the data (red).

In the text
thumbnail Fig. 8.

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

In the text
thumbnail Fig. 9.

Distribution of He in the upper atmosphere of HAT-P-11 b based on the best-fit solution obtained by fitting p-winds models (no limb darkening included) to the CARMENES transmission spectrum (Model 1).

In the text
thumbnail Fig. B.1.

Posterior distribution of parameters of HAT-P-11 b fit to the CARMENES transmission spectrum using a model with a quadratic limb-darkening law and coefficients c1 = 0.63 and c2 = 0.09.

In the text
thumbnail 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
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.