The Upper Edge of the Neptune Desert Is Stable Against Photoevaporation

To construct our survey sample, we searched the NASA Exoplanet Archive (Akeson et al. 2013; NASA Exoplanet Archive 2022) for transiting giant (R > 4R_⊕) planets with measured masses orbiting K-type host stars (T_⋆ = 4000–5400 K). Planets orbiting K-type stars are predicted to have the largest fractional metastable helium populations, and therefore to exhibit the strongest helium absorption signatures during transit (Oklopčić 2019; Wang & Dai 2021a). This is because K stars have relatively high extreme-UV (EUV) fluxes, which populate the metastable state via ground-state ionization and subsequent recombination, and relatively low mid-UV fluxes, which ionize and depopulate the metastable state. A large fraction of the EUV flux also comes from coronal iron lines, and Poppenhaeger (2022) highlighted that this consideration also favors K stars for helium studies. The stellar spectral type gives only an approximate expectation for the high-energy spectrum of the star, however, and the discovery of a strong helium signature for HAT-P-32b (a gas-giant planet orbiting an F star) is a reminder that other stars with favorable spectral energy distributions (SEDs) can also exhibit strong signals regardless of their spectral type (Czesla et al. 2022). Nevertheless, the bulk of planets with detected outflows in the literature orbit K stars.

We used this initial sample to create a rank-ordered list based on predicted signal-to-noise ratio (S/N). We estimated the S/N by first calculating lower atmospheric planetary scale heights using masses, radii, and equilibrium temperatures from the NASA Exoplanet Archive with a mean molecular weight of 2.3 amu. We then assumed that a typical helium absorption extends 85.5 lower atmospheric scale heights at the core of the helium feature, corresponding to the magnitude of the measured excess absorption signal for WASP-69b (Nortmann et al. 2018). We selected this planet as a benchmark simply because it was among the first published measurements of a planetary helium line at high resolving power, but a different choice for the signal amplitude would not have affected the relative rankings. Assuming a line shape similar to that observed by Nortmann et al. (2018), we convolved the expected absorption features with our measured filter transmission profile to get a Palomar/WIRC signal estimate. We calculated the expected photometric noise for each observation by scaling on-sky noise statistics from previous observations (Vissapragada et al. 2020a, 2020b) to the Two Micron All Sky Survey (2MASS) J magnitudes for each host star. From the resulting rank-ordered list, we removed two targets with extensive previous observations (WASP-107b and HD 189733b; Salz et al. 2018; Spake et al. 2018, 2021; Allart et al. 2019; Guilluy et al. 2020; Kirk et al. 2020) and one target that lacked nearby comparison stars, which are required for our photometric observation technique, and had an exceptionally long transit duration (KELT-11b).

We were able to obtain observing time for eight of the highest-priority targets with predicted S/N > 3, one of which (HAT-P-12b) we did not ultimately observe due to telescope closures and poor weather conditions. In Table 1, we present the stellar and planetary parameters for the remaining seven targets constituting our sample. All of the new targets in our sample except for NGTS-5b (a relatively recent discovery) were independently identified by Kirk et al. (2020) as good candidates for metastable helium observations. These seven planets lie near the edge of the Neptune desert, visualized in Figure 1. HAT-P-26b is close to the lower edge of the desert, and the other six planets trace the upper edge.

Table 1. Stellar and Planetary Parameters for Sample Selected in This Work

Star	J mag	T⋆ (K)	R⋆ (R☉)	M⋆ ( M☉)	$\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })$	Mp (MJ)	Rp (RJ)	ρp (g cm−3)	P (days)	a (au)	Reference
WASP-69	8.0	${4715}_{-50}^{+50}$	${0.813}_{-0.028}^{+0.028}$	${0.826}_{-0.029}^{+0.029}$	−4.54	${0.260}_{-0.017}^{+0.017}$	${1.057}_{-0.047}^{+0.047}$	${0.291}_{-0.041}^{+0.041}$	3.87	0.045	A14
WASP-52	10.6	${5000}_{-100}^{+100}$	${0.79}_{-0.02}^{+0.02}$	${0.87}_{-0.03}^{+0.03}$	−4.4	${0.459}_{-0.021}^{+0.022}$	${1.27}_{-0.03}^{+0.03}$	${0.29}_{-0.03}^{+0.03}$	1.75	0.027	H13
HAT-P-18	10.8	${4803}_{-80}^{+80}$	${0.749}_{-0.037}^{+0.037}$	${0.770}_{-0.031}^{+0.031}$	−4.73	${0.197}_{-0.013}^{+0.013}$	${0.995}_{-0.052}^{+0.052}$	${0.25}_{-0.04}^{+0.04}$	5.51	0.056	H11a, V22
WASP-80	9.2	${4143}_{-94}^{+92}$	${0.586}_{-0.018}^{+0.017}$	${0.577}_{-0.054}^{+0.051}$	−4.04	${0.538}_{-0.036}^{+0.035}$	${0.999}_{-0.031}^{+0.030}$	${0.717}_{-0.032}^{+0.039}$	3.07	0.034	T15, F21
WASP-177	10.7	${5017}_{-70}^{+70}$	${0.885}_{-0.046}^{+0.046}$	${0.876}_{-0.038}^{+0.038}$	⋯	${0.508}_{-0.038}^{+0.038}$	${1.58}_{-0.36}^{+0.66}$	${0.172}_{-0.110}^{+0.203}$	3.07	0.040	T19
HAT-P-26	10.1	${5079}_{-88}^{+88}$	${0.788}_{-0.043}^{+0.098}$	${0.816}_{-0.033}^{+0.033}$	−4.99	${0.059}_{-0.007}^{+0.007}$	${0.565}_{-0.032}^{+0.072}$	${0.40}_{-0.10}^{+0.10}$	4.23	0.048	H11b
NGTS-5	12.1	${4987}_{-41}^{+41}$	${0.739}_{-0.012}^{+0.014}$	${0.661}_{-0.061}^{+0.068}$	−4.63	${0.229}_{-0.037}^{+0.037}$	${1.136}_{-0.023}^{+0.023}$	${0.193}_{-0.034}^{+0.032}$	3.36	0.038	E19

Note. J magnitudes are taken from 2MASS (Cutri et al. 2003). Orbital periods and semimajor axes are truncated with errors omitted for clarity. In the references column, A14 is Anderson et al. (2014), H13 is Hébrard et al. (2013), H11a is Hartman et al. (2011), V22 is Vissapragada et al. (2022), T15 is Triaud et al. (2015), F21 is Fossati et al. (2022), T19 is Turner et al. (2019), H11b is Hartman et al. (2011), and E19 is Eigmüller et al. (2019). The $\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })$ value for NGTS-5b was calculated using an S-value from P. Eigmüller (2022, private communication).

Download table as:  ASCIITypeset image

From 2019–2021, we observed transits for the seven targets in our sample over 12 nights. We summarize these observations in Table 2. All transits were observed in a custom metastable helium filter centered at 1083.3 nm with an FWHM of 0.635 nm, previously introduced in Vissapragada et al. (2020b). Each observation followed a similar sequence to that described in our previous papers (Vissapragada et al. 2020b, 2021; Paragas et al. 2021). We began each night by observing a helium arc lamp through the helium filter, allowing us to position the target star where the effective transmission function of the filter was centered on the helium feature. The exact positioning within this region was selected to maximize the number of suitable comparison stars; the number of comparison stars for each night is noted in Table 2. On most nights, after acquiring and positioning the target star, we used a custom near-infrared (NIR) beam-shaping diffuser to mold the PSF into a top-hat 3'' in diameter, allowing for precise control of time-correlated systematics (Stefansson et al. 2017; Vissapragada et al. 2020a). The only exception was our observation of HAT-P-18b on UT 2020 June 5. Weather conditions were poor that night, and we chose to slightly defocus the telescope to 12 instead of using the diffuser in order to minimize the sky background flux in our photometric aperture. As described in Vissapragada et al. (2020b), the central wavelength of the filter shifts across the field of view, and the resulting sky background is highly structured with bright rings corresponding to OH emission features. To correct for the structured background, we obtained a number of dithered background frames on each night (noted in Table 2).

We dark-corrected and flat-fielded each image, and then corrected detector cosmetics including bad pixels and residual striping (for a more detailed description, see Vissapragada et al. 2020a). Because the effective bandpass of the filter changes across the field of view, we also corrected for telluric emission lines and time-varying telluric water absorption. This is necessary because comparison stars occupy different parts of the field and thus experience different amounts of OH emission and time-varying water absorption (Vissapragada et al. 2020b). To correct the sky background and telluric emission lines from each science frame, we first constructed a background frame from the aforementioned dither sequence. Then, we sigma-clipped the science frame to remove sources, and finally we median-scaled the dither frame to match the science frame in 10 px radial steps from the "zero-point," which is where light encounters the filter at normal incidence (Vissapragada et al. 2020b). We demonstrated in Paragas et al. (2021) that the scaling factors from this procedure can also be used to decorrelate the time-varying telluric water absorption, which can otherwise contaminate the transit signal. The OH line at 1083.4 nm overlaps with nearby telluric water features at the resolution of our filter, and its flux evolves differently over the night relative to the uncontaminated OH lines. By taking the ratio of scaling factors in the contaminated and uncontaminated OH lines, we can construct a proxy for the time-varying telluric water absorption. This proxy is included as a covariate in our light-curve fits.

We performed aperture photometry on each source using the photutils package, following the procedure detailed in Vissapragada et al. (2020b). We tested apertures from 3 to 15 pixels (075–375) in radius for each transit observation. We allowed the locations of the apertures to shift (separately for each star) from image to image to account for variations in telescope pointing. We found that the pointing variations were smaller than 1'' across all data sets except for our observation of WASP-80b on UT 2020 July 9, where the instrument crashed during the observation and had to be reset. We included the centroid offsets as covariates in our light-curve fits along with the airmass curves, which we found to be important for a subset of the observations. To select the optimal aperture for use in light-curve modeling, we first normalized the target star's light curve by an average of the comparison star light curves, and subsequently removed 4σ outliers using a moving-median filter. We then selected the aperture size that minimized the per-point rms scatter in the averaged and moving-median-corrected target star photometry. The optimized aperture sizes are noted in Table 2.

After calculating the optimized photometry for each transit observation along with the corresponding decorrelation vectors (airmass, centroid offset, and absorption proxy), we proceeded to fit the transit light curves using exoplanet (Foreman-Mackey et al. 2021). Our procedure is similar to that described in Paragas et al. (2021) and Vissapragada et al. (2021). Briefly, each target light curve is modeled with a limb-darkened starry transit light curve (Luger et al. 2019), which is multiplied by a systematics model. The transit light curve is parameterized by the radius ratio R_p/R_⋆, the orbital period P, the epoch T₀, the scaled semimajor axis a/R_⋆, the impact parameter b, and the quadratic limb-darkening coefficients (u₁, u₂).

We experiment with leaving the limb-darkening coefficients free (sampling them using the approach from Kipping 2013) and fixing them to values computed from ldtk (Husser et al. 2013; Parviainen & Aigrain 2015). Both approaches have advantages and drawbacks. Using the computed values can be advantageous for faint stars, where the precision on the transit shape is too low for solely our data to constrain the limb darkening. On the other hand, the PHOENIX models (Husser et al. 2013) underlying ldtk may not accurately describe the stellar brightness profile in the helium line because the line is formed in the chromosphere. For some stars in our sample, we observe mismatches between the ldtk calculations and the retrieved limb-darkening parameters; this could either be due to strong stellar limb darkening or an optically thin component of the outflow altering the transit shape and biasing the inference for the observed limb-darkening parameters (e.g., Wang & Dai 2021a, 2021b; MacLeod & Oklopčić 2022). Therefore, we repeat the fits for each planet leaving the limb-darkening coefficients free and holding them fixed to the calculated values, and adopt the fixed solution when it agrees with the free retrieval or when the free retrieval is otherwise uninformative on the limb-darkening coefficients. We show how the posteriors on the limb-darkening coefficients compare to the calculated values in Appendix A.

We fix the eccentricities to zero in our fits; none of the planets in our sample are known to have eccentric orbits. Although we note that there is tentative evidence (1σ–2σ) that HAT-P-18b and HAT-P-26b may have small but nonzero orbital eccentricities (Hartman et al. 2011, 2011; Wallack et al. 2019), the eccentricities are too small to have an impact on the light-curve modeling. We modeled the systematics as a linear combination of comparison star light curves and decorrelation vectors, with the weights left as free parameters in each fit. We also included the mean-subtracted Barycentric Julian Date (BJD) times as a decorrelation vector, which functions as a linear baseline. Each weight was assigned a uniform prior of ${ \mathcal U }(-2,2)$. We list the rest of the priors for the transit light-curve modeling in Table 3.

Table 3. Priors for the Transit Light-curve Fits

Planet	P (days)	T0 (BJD—2,450,000)	b	a/R⋆	WIRC Rp/R⋆	References
WASP-69b	${ \mathcal N }(3.8681390,0.0000006)$	${ \mathcal N }(7176.17789,0.00017)$	${ \mathcal N }(0.686,0.023)$	${ \mathcal N }(12.00,0.46)$	${ \mathcal U }(0.,0.25)$	A14, K22
WASP-52b	${ \mathcal N }(1.74978119,0.00000010)$	${ \mathcal N }(6770.05972,0.00004)$	${ \mathcal N }(0.60,0.02)$	${ \mathcal N }(7.38,0.11)$	${ \mathcal U }(0.,0.25)$	H13, K22
HAT-P-18b	${ \mathcal N }(5.5080300,0.0000008)$	${ \mathcal N }(7276.25646,0.00010)$	${ \mathcal N }(0.352,0.057)$	${ \mathcal N }(16.39,0.24)$	${ \mathcal U }(0.,0.25)$	K17, K22
WASP-80b	${ \mathcal N }(3.06785271,0.00000019)$	${ \mathcal N }(6671.49615,0.00004)$	${ \mathcal N }(0.215,0.022)$	${ \mathcal N }(12.63,0.13)$	${ \mathcal U }(0.,0.25)$	T15, K22
WASP-177b	${ \mathcal N }(3.071722,0.000001)$	${ \mathcal N }(7994.37140,0.00028)$	${ \mathcal N }(0.980,0.092)$	${ \mathcal N }(9.610,0.530)$	${ \mathcal U }(0.,0.4)$	T19
HAT-P-26b	${ \mathcal N }(4.2345002,0.0000007)$	${ \mathcal N }(6892.59046,0.00010)$	${ \mathcal N }(0.303,0.122)$	${ \mathcal N }(13.06,0.83)$	${ \mathcal U }(0.,0.25)$	H11, K22
NGTS-5b	${ \mathcal N }(3.3569866,0.0000026)$	${ \mathcal N }(7740.35262,0.00026)$	${ \mathcal N }(0.653,0.048)$	${ \mathcal N }(11.111,0.32)$	${ \mathcal U }(0.,0.4)$	E19

Note. ${ \mathcal U }(a,b)$ indicates a uniform prior between a and b, and ${ \mathcal N }(a,b)$ indicates a normal prior with mean a and standard deviation b. In the References column, A14 is Anderson et al. (2014), K22 is Kokori et al. (2022), H13 is Hébrard et al. (2013), K17 is Kirk et al. (2017), T15 is Triaud et al. (2015), T19 is Turner et al. (2019), H11 is Hartman et al. (2011), and E19 is Eigmüller et al. (2019).

Download table as:  ASCIITypeset image

After defining the model, we sampled the posterior distributions for each model parameter using the No U-Turn Sampler (Hoffman & Gelman 2011) implemented in pymc3 (Salvatier et al. 2016). In each light-curve fit, we ran four chains for 1500 tuning steps each before taking 2500 draws. This resulted in a total of 10,000 draws, which we used to derive the final results. When performing joint fits across multiple nights of WIRC data, we doubled the number of tuning steps and draws. We verified the convergence of each fit by visually inspecting the trace plots to ensure that the chains were well mixed, and also by checking that the Gelman-Rubin statistic (Gelman & Rubin 1992) was $\hat{R}\ll 1.01$ for every sampled parameter.

In addition to the comparison star photometry, each night of data has three additional covariates as described in Section 2.3: the airmass curve, the water absorption proxy, and the centroid offsets. We determined which of these covariates to include in the final model by repeating the fit with every combination of these three parameters and calculated the Bayesian Information Criterion (BIC; Schwarz 1978) for each fit:

$$\begin{eqnarray}&&\mathrm{BIC}=k\mathrm{ln}n-2\mathrm{ln}\hat{{ \mathcal L }},\end{eqnarray} \tag{ 1 }$$

where k is the number of fit parameters, n is the number of data points, and $\hat{{ \mathcal L }}$ is the maximum likelihood estimate, or MLE. We report the differences in BIC (ΔBIC) between each fit and a fit with no additional covariates in Appendix B, and select the model with the minimal BIC. Because pymc3 takes a fully Bayesian approach, we actually obtain an estimate of the maximum a posteriori solution, and the likelihood at this point can differ from the MLE when the priors are informative. To ensure that our model selection is not strongly impacted by this difference, we also consider the two Bayesian model selection methodologies provided by the arviz package: the Pareto-smoothed importance sampling leave-one-out statistic (PSIS-LOO; Vehtari et al. 2017) and the Watanabe-Akaike Information Criterion, sometimes referred to as the Widely Applicable Information Criterion (WAIC; Watanabe 2010). Both are methods for evaluating the out-of-sample predictive accuracy of a model, and more detailed comparisons between these statistics and the BIC can be found in, e.g., Gelman et al. (2014). For each data set, we verified that the detrending model with the lowest BIC also had the highest PSIS-LOO and WAIC values (within the uncertainties reported by arviz), indicating the highest predictive accuracy.

After selecting optimized detrending models for each night of data, we obtained final fits for each planet, in some cases jointly fitting multiple nights of data. In the final version of the fits, we also included a jitter term on each night $\mathrm{log}({\sigma }_{\mathrm{extra}})$, which quantifies the discrepancy between the photon noise and the true variance in the data (i.e., σ² = ${\sigma }_{\mathrm{photon}}^{2}$ + ${\sigma }_{\mathrm{extra}}^{2}$). Final σ values for each light curve are provided in Table 2, as well as the ratio σ/σ_phot between the achieved rms and the photon noise in the unbinned residuals. Allan deviation plots for each light curve are presented in Appendix C.

We summarize the resulting posteriors for our light-curve fits in Table 4. In this table, we also report the excess absorption at mid-transit, δ_mid. To derive this parameter, we first construct a comparison light curve at each step in the fit using the radius ratio measured in a nearby bandpass by another instrument along with our transit shape and limb-darkening parameters. We then subtract the minimum value of the comparison light curve from the minimum value of our helium light curve at each step in the fit to obtain a distribution of δ_mid values, which is summarized in Table 4. To construct the comparison light curves for WASP-69b, WASP-52b, HAT-P-18b, WASP-80b, and HAT-P-26b, we used the radius ratios from the Hubble Space Telescope (HST)/WFC3 G141 between 1110.8 nm and 1141.6 nm obtained by Tsiaras et al. (2018), i.e., the first point of the corresponding spectra in their Figure 4. We chose this wavelength bin because it was closest to 1083.3 nm. We also note that the water absorption features from the lower atmosphere are of the order of ∼100 ppm for these planets, far smaller than our precision on the helium transit depth, making the choice of baseline a negligible source of error in our analysis. WASP-177b and NGTS-5b are relatively recent discoveries and have not been observed by HST; for the former planet, we fit jointly with existing TESS (bandpass between 600 and 1000 nm) observations of this system, and for the latter planet, we use the reported radius ratio from NGTS (an average across optical and NIR wavelengths; Eigmüller et al. 2019) since TESS has not yet observed the system.

Table 4. Posteriors for the Light-curve Fits

Planet	LDC	P (days)	T0 (BJD-2,450,000)	b	a/R⋆	WIRC Rp/R⋆	u1	u2	δmid (%)
WASP-69b	free	${3.86813899}_{-0.00000054}^{+0.00000053}$	${7176.17789}_{-0.00016}^{+0.00017}$	${0.679}_{-0.019}^{+0.019}$	${11.64}_{-0.28}^{+0.28}$	${0.1462}_{-0.0025}^{+0.0026}$	${0.47}_{-0.30}^{+0.33}$	${0.20}_{-0.42}^{+0.40}$	${0.506}_{-0.066}^{+0.068}$
	fixed	${3.86813940}_{-0.00000052}^{+0.00000051}$	${7176.17797}_{-0.00016}^{+0.00016}$	${0.691}_{-0.017}^{+0.017}$	${11.79}_{-0.28}^{+0.27}$	${0.1467}_{-0.0017}^{+0.0017}$	0.38	0.12	${0.512}_{-0.048}^{+0.049}$
WASP-52b	free	${1.749781181}_{-0.000000099}^{+0.000000099}$	${6770.05972}_{-0.000040}^{+0.000040}$	${0.596}_{-0.018}^{+0.018}$	${7.40}_{-0.10}^{+0.10}$	${0.1728}_{-0.0039}^{+0.0038}$	${0.26}_{-0.19}^{+0.28}$	${0.12}_{-0.25}^{+0.31}$	${0.28}_{-0.14}^{+0.14}$
	fixed	${1.749781183}_{-0.000000099}^{+0.000000099}$	${6770.059719}_{-0.000040}^{+0.000039}$	${0.594}_{-0.016}^{+0.017}$	${7.373}_{-0.096}^{+0.096}$	${0.1729}_{-0.0036}^{+0.0035}$	0.35	0.12	${0.29}_{-0.13}^{+0.13}$
HAT-P-18b	free	${5.50802974}_{-0.00000077}^{+0.00000078}$	${7276.256448}_{-0.000098}^{+0.000099}$	${0.365}_{-0.049}^{+0.047}$	${16.44}_{-0.23}^{+0.23}$	${0.1567}_{-0.0047}^{+0.0046}$	${0.67}_{-0.32}^{+0.29}$	${0.08}_{-0.37}^{+0.41}$	${0.70}_{-0.16}^{+0.16}$
	fixed	${5.50802973}_{-0.00000077}^{+0.00000077}$	${7276.256446}_{-0.000097}^{+0.000099}$	${0.412}_{-0.041}^{+0.040}$	${16.54}_{-0.23}^{+0.23}$	${0.1624}_{-0.0036}^{+0.0035}$	0.38	0.12	${0.84}_{-0.13}^{+0.13}$
WASP-80b	free	${3.06785259}_{-0.00000018}^{+0.00000019}$	${6671.49614}_{-0.00004}^{+0.00004}$	${0.216}_{-0.022}^{+0.022}$	${12.66}_{-0.12}^{+0.12}$	${0.1708}_{-0.0055}^{+0.0055}$	${0.40}_{-0.27}^{+0.36}$	${0.21}_{-0.40}^{+0.38}$	$-{0.01}_{-0.22}^{+0.23}$
	fixed	${3.06785259}_{-0.00000018}^{+0.00000019}$	${56671.496140}_{-0.000040}^{+0.000040}$	${0.219}_{-0.022}^{+0.021}$	${12.68}_{-0.11}^{+0.11}$	${0.1716}_{-0.0051}^{+0.0051}$	0.31	0.15	${0.02}_{-0.20}^{+0.20}$
WASP-177b	free	${3.07172220}_{-0.00000068}^{+0.00000068}$	${7994.37141}_{-0.00024}^{+0.00025}$	${0.980}_{-0.057}^{+0.071}$	${9.51}_{-0.30}^{+0.37}$	${0.230}_{-0.043}^{+0.056}$	${0.62}_{-0.45}^{+0.59}$	$-{0.03}_{-0.42}^{+0.43}$	${0.70}_{-0.35}^{+0.35}$
	fixed	${3.07172221}_{-0.00000069}^{+0.00000068}$	${7994.37142}_{-0.00024}^{+0.00025}$	${0.975}_{-0.054}^{+0.069}$	${9.52}_{-0.31}^{+0.36}$	${0.213}_{-0.034}^{+0.050}$	0.36	0.12	${0.53}_{-0.28}^{+0.23}$
HAT-P-26b	free	${4.23450038}_{-0.00000067}^{+0.00000068}$	${6892.590466}_{-0.000097}^{+0.000100}$	${0.288}_{-0.111}^{+0.099}$	${13.19}_{-0.56}^{+0.51}$	${0.0878}_{-0.0058}^{+0.0055}$	${0.26}_{-0.19}^{+0.32}$	${0.15}_{-0.27}^{+0.36}$	${0.30}_{-0.11}^{+0.11}$
	fixed	${4.23450039}_{-0.00000066}^{+0.00000067}$	${6892.590466}_{-0.000099}^{+0.000099}$	${0.290}_{-0.110}^{+0.099}$	${13.22}_{-0.52}^{+0.48}$	${0.0882}_{-0.0054}^{+0.0051}$	0.35	0.12	${0.31}_{-0.10}^{+0.10}$
NGTS-5b	free	${3.3569892}_{-0.0000024}^{+0.0000024}$	${7740.35267}_{-0.00026}^{+0.00026}$	${0.675}_{-0.046}^{+0.044}$	${11.16}_{-0.31}^{+0.31}$	${0.188}_{-0.014}^{+0.014}$	${0.81}_{-0.55}^{+0.57}$	$-{0.10}_{-0.46}^{+0.52}$	${1.07}_{-0.50}^{+0.52}$
	fixed	${3.3569890}_{-0.0000024}^{+0.0000024}$	${7740.35267}_{-0.00026}^{+0.00026}$	${0.691}_{-0.044}^{+0.041}$	${11.20}_{-0.31}^{+0.31}$	${0.187}_{-0.013}^{+0.013}$	0.36	0.12	${1.02}_{-0.46}^{+0.48}$

Note. The "LDC" column denotes the fitting strategy of leaving the limb-darkening coefficients free or fixing them to a calculation from ldtk (Husser et al. 2013; Parviainen & Aigrain 2015); the adopted strategy for each planet is noted in bold. The mid-transit excess depth δ_mid is a derived parameter. Posteriors on the detrending weights are omitted for brevity.

Download table as:  ASCIITypeset image

WASP-69b is a Saturn-mass (0.26M_J), Jupiter-sized (1.06R_J) planet in a 3.87 day orbit around its K5 host star (Anderson et al. 2014). Its low gravitational potential makes it an excellent target for atmospheric observations. To date, observations of this planet's lower atmosphere have revealed water absorption, sodium absorption, and a Rayleigh scattering slope indicating the presence of hazes (Tsiaras et al. 2018; Murgas et al. 2020; Estrela et al. 2021; Khalafinejad et al. 2021). Additionally, Nortmann et al. (2018) detected helium escaping from the upper atmosphere of WASP-69b with two nights of CARMENES observations. We confirmed the strong absorption reported by these authors in our observation of this planet on UT 2019 August 16, which was initially published in Vissapragada et al. (2020b).

Our transit modeling methodology has changed since that initial work, and we therefore re-fit the data to be consistent with the other planets in the sample. We found that the calculated limb-darkening coefficients agreed well (within 2σ) with the joint distribution of retrieved quadratic limb-darkening coefficients (Appendix A), and the choice does not make much of a difference for the retrieved parameters. Therefore, we adopted the fixed limb-darkening coefficient fit for this planet. The best-fit light curve, residual, and Allan deviation plot are shown in Figure 2(a). We find the excess depth at mid-transit to be ${0.512}_{-0.048}^{+0.049}$% in our bandpass, which is consistent with our previous work albeit with slightly larger uncertainty. In our previous work, the detrending model used the best-fit linear combination of comparison star vectors, effectively marginalizing over the detrending vector weights that we fit for in this work, so the previously reported uncertainties were slightly underestimated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

WASP-52b is a 0.46M_J, 1.27R_J planet in a short 1.75 day orbit around its K2 host star (Hébrard et al. 2013). Notably, the system appears to be relatively young: the discovery paper quotes a gyrochronological age of ${400}_{-200}^{+300}$ Myr (Hébrard et al. 2013), though this may be an underestimate (Mancini et al. 2017). Its precarious position—just past the edge of the Neptune desert—was also noted by Owen & Lai (2018). In the high-eccentricity migration framework, which is their preferred scenario for planets at the upper edge of the desert, this planet's position past the edge of the desert suggests that it should have been tidally disrupted. To resolve this discrepancy, Owen & Lai (2018) suggest that WASP-52b's radius may have been smaller during migration (with a correspondingly smaller tidal disruption radius), and that it underwent radius inflation only recently. Other published observations of WASP-52b indicate that this planet has a relatively flat transmission spectrum (Kirk et al. 2016; Chen et al. 2017; Louden et al. 2017; Mancini et al. 2017; Alam et al. 2018; May et al. 2018), although sodium, potassium, and Hα absorption have all been detected at high spectral resolution (Chen et al. 2017, 2020; Alam et al. 2018).

We first searched for helium in the atmosphere of this planet on UT 2019 September 17, and published the results of these observations in Vissapragada et al. (2020b). In this work, we re-fit that light curve jointly with a new observation taken on UT 2020 August 4. Observing conditions were worse on the second night, which rendered two faint comparison stars from the first night unobservable, but better positioning of the target star allowed us to add one extra comparison star that was not in the field of view on the first night.

For this planet, we adopted a fixed limb-darkening coefficient fit, as we found that the calculated limb-darkening coefficients agreed well with the joint distribution of retrieved quadratic limb-darkening coefficients (Appendix A), and the choice did not make a difference for the retrieved parameters. In separate retrievals, we constrained the excess absorption to be ${0.22}_{-0.13}^{+0.14}$% on the first night, with a 95th percentile upper limit of 0.44%, and ${1.02}_{-0.41}^{+0.39}$% on the second night with a 95th percentile upper limit of 1.68%. When fitting the data sets jointly, we found an excess absorption of ${0.29}_{-0.13}^{+0.13} \%$, corresponding to a tentative detection at 2.2σ confidence. The joint fit results are shown in Figure 2(b).

HAT-P-18b is a Jupiter-sized (1.00R_J), low-mass (0.20M_Jup) planet in a 5.5 day orbit around its K2 host star (Hartman et al. 2011). Despite its faintness (J = 10.8), the low gravitational potential of this planet makes it amenable to atmospheric characterization, and a Rayleigh scattering slope in HAT-P-18b's transmission spectrum was detected by Kirk et al. (2017). We previously detected a metastable helium signature for this planet in Paragas et al. (2021). The fitting method we employ in this survey is somewhat different from that in Paragas et al. (2021), so we re-fit the observations here. Rather than fitting jointly with TESS observations of the system as we did previously, we fit the WIRC data alone and compare to the precise transit depth from HST/WFC3 G141 (Tsiaras et al. 2018). We prefer this approach because the spectrophotometric bin we use from WFC3 (1110.8–1141.6 nm) is narrower and closer to the helium filter bandpass than the TESS bandpass (600–1000 nm), and the depth in this bin is more precise than the TESS depth (∼100 ppm precision for WFC3 versus ∼500 ppm for TESS). Additionally, in the previous work we reported $({R}_{{\rm{p}}}/{R}_{\star }{)}_{\mathrm{WIRC}}^{2}-({R}_{{\rm{p}}}/{R}_{\star }{)}_{\mathrm{comparison}}^{2}$ as the excess absorption, but this is not necessarily equal to the difference in transit depth between the WIRC light curve and comparison light curve, especially when there is strong stellar limb darkening in the helium line (as we infer for HAT-P-18b). To account for this, the comparison light curve should be constructed taking into account limb darkening in the helium bandpass, which we do when obtaining the excess absorption δ_mid in this work (see Section 2.4).

For HAT-P-18b, we adopted the free limb-darkening coefficient fit. The fixed limb-darkening coefficients disagreed with the retrieved coefficients by more than 2σ in the joint posterior (Appendix A), and the choice was consequential for the final fitted parameters, so we choose to adopt the free limb-darkening coefficient solution. When fitting each data set independently, we obtained ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$ values of ${2.14}_{-0.24}^{+0.23}$% and ${2.56}_{-0.16}^{+0.17}$% for the first and second nights, respectively. These are consistent to 1σ with the values obtained by Paragas et al. (2021), who fit the WIRC data together with data from TESS to obtain ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$ values of 2.11% ± 0.25% for the first night and 2.35% ± 0.14% for the second night. We then fit the WIRC data sets jointly and show the results in Figure 2(c). In the joint fit, we find ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}=2.46\pm 0.15$%, which is about 1σ larger than the result from Paragas et al. (2021) of 2.29% ± 0.12%. Our result for the excess depth is 0.70% ± 0.16%, again slightly larger than our previously reported signal of 0.46% ± 0.12%.

WASP-80b is a Jupiter-sized (1.00R_J), 0.54M_J planet in a 3.07 day orbit around its K7/M0 host star (Triaud et al. 2013, 2015). At low resolving power, the transmission spectrum of this planet is quite flat (Fukui et al. 2014; Mancini et al. 2014; Turner et al. 2017; Kirk et al. 2018; Parviainen et al. 2018). Multiple groups running self-consistent wind-launching simulations have reported WASP-80b to be an excellent candidate for mass-loss observations (Caldiroli et al. 2022; Salz et al. 2015, 2016b). However, Fossati et al. (2022) recently obtained four nights of high-resolution data on WASP-80b and did not detect any metastable helium signature, setting an upper limit of 0.7% absorption in the line core.

We observed a full transit of WASP-80b on UT 2020 July 9. During the first half of observations, the detector started reading large numbers of negative counts relatively infrequently until the issue became more persistent near ingress. The instrument was reset near the middle of the transit, which fixed the problem. Prior to analysis, frames in the first half with >0.1% of pixels reading negative counts were discarded. Additionally, the pre-reset data exhibited stronger systematics than the post-reset data, so we treated each of the two halves as separate observations (i.e., allowing each half of the light curve to have different coefficients for the systematics model) and fit them jointly. We optimized the aperture size using the second half of data only, finding an optimized radius of 8 px.

We adopted the fixed limb-darkening coefficient fit for WASP-80b, as the calculated limb-darkening coefficients agreed well with the joint distribution of retrieved quadratic limb-darkening coefficients (Appendix A), and the choice did not make a difference for the retrieved parameters. We did not detect any excess absorption in our data, finding ${\delta }_{\mathrm{mid}}={0.02}_{-0.20}^{+0.20}$ with a 95th percentile upper limit of 0.35% excess absorption in our bandpass. The fit results are shown in Figure 2(d). Our nondetection agrees with the more sensitive observations from Fossati et al. (2022).

WASP-177b is a highly inflated (1.6R_J), 0.5M_J planet in a 3.07 day orbit around its K2 host star (Turner et al. 2019). This planet has an exceptionally low density (similar to WASP-107b), making it a good target for transmission spectroscopy. However, the planetary transit is grazing, which makes it more difficult to extract precise constraints on the wavelength-dependent planet–star radius ratio.

We observed a full transit of this planet on UT 2020 October 2. TESS also observed this planet (TOI-4521.01) at 2 minutes cadence in Sector 42 from UT 2021 August–2021 September. Given the grazing nature of the transit and the relatively small number of comparison light curves in the literature, we decided to fit our Palomar/WIRC light curve jointly with the TESS light curve. This allowed us to better capture covariances between the radius ratio, impact parameter, scaled semimajor axis, and limb-darkening parameters while also giving us a reference transit depth to which we could compare the WIRC measurement. The TESS portion of the fit was carried out similarly to our fit for HAT-P-18b in Paragas et al. (2021); the only additional parameters we included were separate TESS limb-darkening parameters, the radius ratio in the TESS bandpass, and an error re-scaling term.

For a grazing geometry it is difficult to differentiate excess absorption associated with the planet from strong limb darkening, so we first tried a fit with free limb-darkening coefficients. However, we found that the data were uninformative on the quadratic limb-darkening coefficients (see Appendix A), so we defaulted to using fixed limb-darkening coefficients instead. The results of the joint fit are shown in Figure 2(e). The fitted TESS light curve is given in Figure 2(f). Using the TESS transit depth as our baseline, we constrain the excess depth to be ${0.53}_{-0.28}^{+0.23}$%, a nondetection with a 95th percentile upper limit of 0.90%.

HAT-P-26b resides near the lower edge of the Neptune desert, and is the smallest planet in our sample. With a radius of a 0.57R_J and a mass of 0.059M_J, this planet orbits its K1 host star every 4.23 days (Hartman et al. 2011). The Neptune-mass planet has a well-constrained heavy element abundance thanks to the presence of strong water bands in the red-optical and NIR (Stevenson et al. 2016; Wakeford et al. 2017) as well as the potential presence of metal hydride features (MacDonald & Madhusudhan 2019) and Spitzer secondary eclipse depth ratios (Wallack et al. 2019). At optical wavelengths, HAT-P-26b's alkali features appear to be obscured by a cloud layer (Wakeford et al. 2017; Panwar et al. 2021).

This planet is an excellent target for observations of atmospheric escape. It is larger and less massive than HAT-P-11b (a benchmark planet for metastable helium studies; Allart et al. 2018; Mansfield et al. 2018), and as a result has an inferred H/He envelope fraction approximately twice that of HAT-P-11b (Lopez & Fortney 2014). The slightly earlier spectral type of HAT-P-26 (HAT-P-11 is a K4 dwarf) also puts it at the "sweet spot" for metastable helium observations as identified by Oklopčić (2019). We observed this high-priority target three times: on UT 2021 February 18, 2021 May 1, and 2021 May 27. The weather on the first night was relatively poor, and the target was not optimally positioned. Between the first and second nights of observation, we shifted the placement of the target on the detector in order to include an additional comparison star, which greatly increased the final precision.

We do not detect the transit in the first night of data, for which we obtained an overall ${R}_{{\rm{p}}}/{R}_{\star }={0.034}_{-0.023}^{+0.028}$. As expected, this first night is much noisier than the latter two (Table 2): the rms scatter in the unbinned residuals was 1.79% (2.7× the photon noise) in the first observation, whereas it was 0.82% (1.6× the photon noise) and 0.61% (1.4× the photon noise) on the second and third nights, respectively. In order to avoid biasing the final joint fit, we therefore only included the latter two nights of photometry.

After optimizing the set of detrending vectors selected for each of the latter two nights, we obtained mid-transit excess depths of ${0.43}_{-0.17}^{+0.18} \%$ and ${0.23}_{-0.13}^{+0.13} \%$. We adopted the fixed limb-darkening coefficient fit for HAT-P-26b, as the calculated limb-darkening coefficients agreed well with the joint distribution of retrieved quadratic limb-darkening coefficients (Appendix A), and the choice did not make a significant difference for the retrieved parameters. The results of our joint fit to the latter two nights are shown in Figure 2(g). We obtained a mid-transit excess depth of ${0.31}_{-0.10}^{+0.10}$, evidence for helium absorption at 3.1σ confidence.

NGTS-5b is a Jupiter-sized (1.1R_J), sub-Saturn-mass (0.23M_J) planet orbiting a K2 dwarf every 3.36 days (Eigmüller et al. 2019). The planet is similar to WASP-69b in mass, radius, orbital period, and stellar host type, so despite the rather faint host star (J = 12.1), we identified it as a high-priority target. We observed two full transits of the planet: one on UT 2021 Apr 30 and one on UT 2021 May 27.

This target had only one comparison star available in the field of view, and this along with the faintness of the target star created challenges for our standard systematics model approach. When testing different covariate combinations for the first night of data, the detrending models failed. With relatively poor weather conditions on that first night, the noisy comparison star light curve was severely underweighted in the systematics model when other covariates were included. When included in the joint fit, this underweighted model severely increased the magnitude of the correlated noise. We stress that unlike the first night for HAT-P-26b, which we discarded because there was no transit detected and the noise was nearly 3× the photon noise limit, the transit is clearly detected in this data set, and the issue stemmed purely from our approach to modeling the systematics. We therefore kept this data set in the final model, but to avoid the sharp increase in correlated noise, we did not use any additional covariates in the systematics model for the first night, and as such, we do not report the model selection statistics for this night in Appendix B. Conditions on the second night were more favorable, so the fit was more typical.

We obtained mid-transit excess absorption values of ${1.02}_{-0.95}^{+0.99}$% and ${0.94}_{-0.54}^{+0.55}$% for the individual fits to the first and second nights, respectively. We proceeded to fit the data sets jointly, and the results are shown in Figure 2(h). We adopted the fixed limb-darkening coefficient fit for NGTS-5b as the calculated limb-darkening coefficients agreed well with the joint distribution of retrieved quadratic limb-darkening coefficients (Appendix A) and the choice did not make a significant difference for the retrieved parameters. We obtained a final mid-transit excess absorption of ${1.02}_{-0.46}^{+0.48}$%, a tentative detection at 2.2σ confidence.

Self-consistent, 1D mass-loss models like the Pluto-CLOUDY Interface (Salz et al. 2015, 2016b) and the ATmospheric EScape code (ATES; Caldiroli et al. 2021, 2022) are crucial for interpretation of mass-loss-related phenomena in the exoplanet population. But to date, there have been relatively few attempts to benchmark their predictions against measurements of present-day mass-loss rates for transiting planets. Here, we compare the inferred mass-loss rates from our seven-planet sample to predictions from the ATES code as a function of the incident XUV flux. Because Caldiroli et al. (2022) provide analytical approximations for the efficiency as a function of incident XUV flux and gravitational potential, comparisons between the model and the data are simple and do not require expensive computations for each planet.

In order to place all seven planets on the same plot, we need to convert our measured mass-loss rates and XUV fluxes into scaled units that correct for planet-to-planet variations in gravitational potential and Roche radius. Scaled units can be obtained by starting from the energy-limited approximation from the planetary mass-loss rate (e.g., Caldiroli et al. 2022):

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{\varepsilon \pi {R}_{{\rm{p}}}^{3}{F}_{\mathrm{XUV}}}{{{KGM}}_{{\rm{p}}}},\end{eqnarray} \tag{ 2 }$$

where ε is the mass-loss efficiency, F_XUV is the incident high-energy flux, and K is the Erkaev et al. (2007) Roche lobe correction factor. Then, we can rewrite the expression in terms of the planetary density ρ_p:

$$\begin{eqnarray}&&K\dot{M}=\displaystyle \frac{3\varepsilon }{4G}\displaystyle \frac{{F}_{\mathrm{XUV}}}{{\rho }_{{\rm{p}}}}.\end{eqnarray} \tag{ 3 }$$

Planets in the energy-limited regime should therefore exhibit a linear relation between the K-corrected mass-loss rate and the ratio F_XUV/ρ_p (Caldiroli et al. 2022).

As the incident XUV flux increases, however, self-consistent models predict that the outflow begins to lose substantial energy to radiative cooling, leading to a sublinear dependence on F_XUV/ρ_p beyond a threshold of F_XUV/ρ_p ∼ 10⁴ erg cm g⁻¹ s⁻¹ (Murray-Clay et al. 2009; Salz et al. 2016a, 2016b; Caldiroli et al. 2022). Higher-gravity planets are also predicted to have low overall outflow efficiencies, as the atmospheres are more tightly bound and tend to re-emit more energy locally (Owen & Alvarez 2016). To an order of magnitude, this happens when the energy liberated by a single photoionization event (ΔE ∼ 10 eV ∼10⁻¹¹ erg) exceeds the particle binding energy, which, for a hydrogen atom (m_H ∼ 10⁻²⁴ g), occurs roughly at a gravitational potential of ϕ_p ∼ ΔE/m_H ∼ 10¹³ erg g⁻¹ (for a more thorough explanation, see Caldiroli et al. 2022). More sophisticated numerical experiments place this threshold closer to ϕ_p ≳ 10^13.2 erg g⁻¹ (Salz et al. 2016b; Caldiroli et al. 2022).

In Figure 4 we plot our marginalized and K-corrected $\dot{M}$ distributions as a function of F_XUV/ρ_p, along with the corresponding ATES model prediction for each planet. Immediately, two distinct populations of planets emerge in this parameter space. The first, comprising HAT-P-18b, WASP-69b, HAT-P-26b, NGTS-5b, and WASP-177b, appear to agree quite well with the ATES predictions. The second group of planets, which includes WASP-80b and WASP-52b, appear to have mass-loss rates that are more than an order of magnitude lower than the model predictions. We discuss possible explanations for the surprisingly low inferred mass-loss rates of the second group in Section 5.2.

Zoom In Zoom Out Reset image size

We use a bootstrap averaging method to estimate the mean mass-loss efficiency of planets in the first group. We omit WASP-52b and WASP-80b from this exercise because their helium signals appear to be affected by factors that are not included in our model, such as magnetic fields or stellar winds, as discussed in the next section. For the other planets, we first converted our marginalized $\dot{M}$ distribution for each planet into a distribution on ε using Equation (3), and these distributions are summarized in Table 5. We then took a random draw from each planet's ε distribution and averaged the draws to get an estimate for the mean efficiency. We repeated this procedure 100,000 times to obtain a final estimate of $\varepsilon ={0.41}_{-0.13}^{+0.16}$. We overplot this average ε model in Figure 4.

The measured mass-loss constraints for all five planets in the first group overlap with the average-ε model at 2σ confidence, and the model is a reasonable approximation to the ATES (Caldiroli et al. 2022) predictions for our sample as well. ATES also predicts that the mass-loss efficiency should level off with increasing F_XUV/ρ_p as compared to a fixed-ε curve. Unfortunately, the uncertainties in the inferred mass-loss rates for our sample, which are dominated by the degeneracy between mass-loss rate and thermospheric temperature, are too large to differentiate between this model and a constant-ε model. Resolving the line profiles of these planets with high-resolution transmission spectroscopy can help break the degeneracy (e.g., dos Santos et al. 2022).

We can compare our measured efficiency to other giant planets (R_p > 4R_⊕) orbiting K stars with reported helium detections in the literature. The reported absorption for WASP-107b (Spake et al. 2018, 2021; Allart et al. 2019; Kirk et al. 2020) is well-matched by the 3D hydrodynamical models of Wang & Dai (2021b), who report a mass-loss rate of approximately 2 × 10¹¹ g s⁻¹. Using the scaled MUSCLES spectrum of the similarly active star Eridani ($\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })=-4.44;$ Noyes et al. 1984) to calculate the XUV flux for WASP-107b ($\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })=-4.43;$ calculated using the S-value from Spake et al. 2018), we obtained an outflow efficiency of 0.34, in good agreement with our measured efficiency. The signal observed for HAT-P-11b (Allart et al. 2018; Mansfield et al. 2018) has been modeled using p-winds using the scaled MUSCLES spectrum of HD 40307 as a proxy for the stellar XUV spectrum, and dos Santos et al. (2022) report a 1D mass-loss rate of ${2.3}_{-0.5}^{+0.7}\times {10}^{10}$ g s⁻¹ (with stellar limb darkening taken into account). Dividing this result by 4 to get the overall mass-loss rate (for consistency with the method outlined in Section 4), we obtain an outflow efficiency of ${0.170}_{-0.037}^{+0.052}$, somewhat smaller than our measured efficiency but consistent within 2σ.

Finally, the signal observed for HD 189733b (Salz et al. 2018; Guilluy et al. 2020) was modeled using a Parker wind methodology by Lampón et al. (2021b). For their 90/10 H/He composition model, they found a mass-loss rate at T = 12,000 K of about 4 × 10⁹ g s⁻¹, with significant spread due to the $\dot{M}-{T}_{0}$ degeneracy. Again dividing by 4 to account for the multidimensional outflow, this corresponds to an efficiency of 0.0023, much smaller than our measured efficiency. A small efficiency is expected: the relatively large mass of HD 189733b puts it in a regime where radiative losses dominate the cooling budget rather than the adiabatic expansion of the atmosphere (Salz et al. 2016b; Caldiroli et al. 2021, 2022), so the outflow is more similar to WASP-52b and WASP-80b than the other planets in our sample. Similarly to WASP-52b and WASP-80b, the planet has a smaller semimajor axis (a = 0.031; Bouchy et al. 2005) and more active host star ($\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })=-4.501;$ Knutson et al. 2010) than the other planets in our sample. The inferred mass-loss rate for HD 189733b at a 90/10 composition is still at least a factor of a few smaller than those predicted by 1D hydrodynamical models, which give mass-loss rates of about 10^9.5−10.0 g s⁻¹ for this planet (Salz et al. 2016b; Caldiroli et al. 2021, 2022). Lampón et al. (2021b) resolved this discrepancy by invoking a hydrogen-rich composition for the outflow, but this is not the only factor that can achieve a reduction of the helium signal. In the next section, we consider a range of explanations for the smaller-than-expected helium signals on WASP-52b and WASP-80b that could apply to HD 189733b as well.

We now consider potential explanations for WASP-80b and WASP-52b's low apparent mass-loss rates, which are discrepant with both the ATES model predictions and the rest of the planets in our sample. We note that these two planets have the smallest semimajor axes and the largest $\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })$ values (Table 1) in our sample; this may or may not have anything to do with their low observed mass-loss rates. In this section, we consider five possible explanations: overestimated XUV fluxes, the XUV spectral shapes of their host stars, stronger-than-predicted stellar winds, outflow composition, and magnetic fields.

If the XUV flux of the host star was overestimated, it would cause us to overpredict the magnitude of a planet's mass-loss rate. To test this effect, we decreased the incident XUV flux by a factor of 10 in the p-winds models for WASP-52b and WASP-80b and recomputed the mass-loss efficiencies. Averaging the two efficiency parameters using the same bootstrapping method from Section 5.1, we found a mean efficiency of ${0.130}_{-0.048}^{+0.064}$. This is still smaller than our inferred efficiency for the other five planets, but the distributions overlap at the 2σ level. We conclude that the XUV fluxes would have to be overestimated by at least a factor of 10 for the mass-loss efficiencies to come into reasonable agreement with the rest of the sample.

For WASP-80, Fossati et al. (2022) combined archival ROSAT and XMM-Newton observations with a new Gaia distance to calculate a precise X-ray flux measurement of ${L}_{{\rm{X}}}={4.85}_{-0.23}^{+0.12}\times {10}^{27}$ erg s⁻¹. They used the King et al. (2018) relations to scale this value into the UV, obtaining a final F_XUV ≈ 6300 erg cm⁻² s⁻¹. Two other papers independently derive this star's X-ray flux and find values that are approximately a factor of 3 smaller (Monsch et al. 2019) and a factor of 2 larger (Foster et al. 2022 obtained a flux close to our adopted value in Table 5) than that obtained by Fossati et al. (2022). Despite these methodological differences, factor-of-a-few variations in the XUV flux are not enough to reconcile the tight upper limit on this planet's helium absorption signal with model predictions. For WASP-52, which was observed by Chandra (ACIS), Monsch et al. (2019) found ${L}_{{\rm{X}}}={3.1}_{-1.0}^{+1.2}\,\times {10}^{28}$ erg s⁻¹. Using the X-ray to EUV scaling relation for Chandra (ACIS) from King et al. (2018), we obtain an XUV irradiance of F_XUV ≈ 38,000 erg cm⁻² s⁻¹. This is quite close to our assumed value (see Table 5). We therefore conclude that, just as with WASP-80, differences in assumed XUV flux are unlikely to explain this planet's low helium absorption signal.

While variations in the integrated XUV flux are unlikely to explain the helium observations for WASP-80b and WASP-52b, the assumed spectral shape may play a role. In order to create a significant population of metastable helium, there must be enough EUV photons to ionize ground-state helium in the outflow, which then recombines into the metastable state (Oklopčić & Hirata 2018; Oklopčić 2019). As demonstrated by Poppenhaeger (2022), many of the stellar EUV photons that ultimately ionize helium come from coronal iron lines. Iron has a relatively low first ionization potential (FIP), and in very active, low-mass stars, the abundance of species with low FIPs often appear to be diminished (termed the inverse first ionization potential effect, or iFIP; Brinkman et al. 2001; Güdel et al. 2001; Wood et al. 2018). WASP-52 and WASP-80 are by far the most active stars in our sample, so they may be affected by this phenomenon. Our assumed proxy star for these two planets, Eridani, does not exhibit the iFIP effect and has strong coronal iron lines (Poppenhaeger 2022); this might cause us to overestimate their metastable helium populations in our models.

It is also worth noting that the models we are using in this study assume spherical symmetry, but not all outflows are spherically symmetric. If the outflow is confined into a comet-like shape by a strong stellar wind, it could reduce the magnitude of the helium absorption (e.g., McCann et al. 2019; Carolan et al. 2020; MacLeod & Oklopčić 2022). This phenomenon has been directly detected for the hot Jupiter WASP-107b (Khodachenko et al. 2021; Spake et al. 2021; Wang & Dai 2021b), and may also explain the relatively weak helium signal measured for the mini-Neptune TOI-560.01 (Zhang et al. 2022). Importantly, both WASP-52 and WASP-80 appear to be relatively young, suggesting that they may have correspondingly enhanced stellar winds. WASP-52 ($v\sin i\,\approx$ 3.6 km s⁻¹ with a photometric rotation period of 16.4 days; Hébrard et al. 2013) and WASP-80 ($v\sin i\,\approx$ 3.55 km s⁻¹; Triaud et al. 2013) both rotate rapidly, with estimated gyrochronological ages of less than 1 Gyr (Barnes 2007; Hébrard et al. 2013; Triaud et al. 2013). They are also the two most active stars in our sample, with B − V colors and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ values that also suggest ages as young as a few hundred megayears (Mamajek & Hillenbrand 2008). However, some caution is warranted as stellar angular momenta are known to be affected by giant planets on close-in orbits (e.g., Lanza 2010; Poppenhaeger & Wolk 2014; Mancini et al. 2017).

Fossati et al. (2022) modeled the impact of stellar wind on WASP-80b using a 3D hydrodynamical model (Khodachenko et al. 2021; Shaikhislamov et al. 2021), and showed that it does not suppress the helium signal enough to match their nondetection. They modeled stellar mass-loss rates up to 10¹³ g s${}^{-1}\sim 8\,{\dot{M}}_{\odot }$, where the solar mass-loss rate is $\,{\dot{M}}_{\odot }\approx 2\times {10}^{-14}\,{M}_{\odot }$ yr⁻¹ (e.g., Wood et al. 2005), but at early times, the wind can be even stronger. Wood et al. (2005) used astrospheric Lyα absorption to infer mass-loss rates of $30\,{\dot{M}}_{\odot }$ and $100\,{\dot{M}}_{\odot }$ for Eridani and 70 Oph, both relatively young K dwarf stars. We can use the aforementioned X-ray fluxes along with the Wood et al. (2005) relation to estimate the stellar mass-loss rate for WASP-80. We obtain a mass-loss rate of $\sim 4\,{\dot{M}}_{\odot }$, well within the range modeled by Fossati et al. (2022). Using the same relation, we obtain a mass-loss rate of $\sim 35\,{\dot{M}}_{\odot }$ for WASP-52. We conclude that a strong stellar wind might suppress the helium signal for WASP-52b, but it is unlikely to explain WASP-80b.

The composition of the outflow could also affect the magnitude of the observed signal. It has been suggested (e.g., Lampón et al. 2020, 2021b) that outflows from giant planets may be depleted in helium, which would suppress both the helium signal and the corresponding inferred mass-loss rate. Fossati et al. (2022) modeled the composition of WASP-80b's outflow and showed that helium depletion (with n_He/n_H ∼ 0.01) might plausibly explain these data. Previous studies of HD 209458b, HD 189733b, and GJ 3470b have sought to quantify the fractional abundances of hydrogen and helium in planetary outflows by comparing the magnitude of absorption measured in both the Lyα and metastable helium lines to Parker wind models (Lampón et al. 2020, 2021a, 2021b) or 3D hydrodynamical models (Shaikhislamov et al. 2021). However, Lyα is a poor tracer for the density distribution in the planetary thermosphere, typically probing above the exobase at high velocities (e.g., Owen et al. 2021). This suggests that anchoring the Parker wind model composition near the wind-launching radius to models of Lyα absorption may be inadvisable. A more reasonable route would be to jointly model Hα and metastable helium absorption signals, as suggested by Czesla et al. (2022). Such joint models would be invaluable in resolving this question.

Strong magnetic fields can also confine planetary outflows. For a strongly ionized wind launched in a dipolar planetary magnetic field, equatorial field lines are closed (ionized material traveling along field lines would loop back to the planet), leading to an equatorial "dead zone," while polar field lines remain open (Adams 2011; Trammell et al. 2011; Owen & Adams 2014). For close-in planets orbiting exceptionally active stars, the outflows are largely ionized, so material should follow field lines if the planets are magnetized. WASP-80b and WASP-52b have the smallest semimajor axes, and their host stars have the largest $\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })$ values in our whole sample, making this an attractive explanation.

The ratio of magnetic pressure to ram pressure is a good diagnostic for the influence of magnetic fields (Owen & Adams 2014), but this ratio depends on the magnetic field strength, which is quite uncertain. There are several lines of evidence indicating that the magnetic fields of hot Jupiters may be relatively strong (10–100 G), including the inflated radii of these planets (Batygin & Stevenson 2010; Yadav & Thorngren 2017) and evidence for magnetic star–planet interactions (SPIs) in a few systems (Cauley et al. 2019). Other observations suggest that these planets may have more modest field strengths (1–10 G), including ultraviolet observations of Lyα and ionized carbon lines in HAT-P-11b (Ben-Jaffel et al. 2022), and a small inferred hotspot offset in the phase curve of WASP-76b, which might be caused by Lorentz drag (May et al. 2021; Beltz et al. 2022). Elsasser number scalings also predict more moderate magnetic field strengths (Stevenson 2003). If the Elsasser number is of order unity, the field strength scales with the square root of the rotation rate; for tidally locked hot Jupiters, this implies $B\sim 1/\sqrt{P}$. Assuming the fluid densities and conductivities are similar to that of Jupiter, we have B ∼ B_J(P/10 hr)^−1/2, of order a few Gauss for the planets in consideration here. Even at the lower end of this range, magnetic pressure is expected to dominate over ram pressure near the wind-launching region (Owen & Adams 2014), suggesting that magnetic fields may indeed play a role in confining the outflows of WASP-80b and WASP-52b.

In summary, inaccuracies in the integrated XUV fluxes are unlikely to explain the helium nondetection for WASP-80b and the small inferred mass-loss rate for WASP-52b, but inaccuracies in the spectral shape could decrease the helium signals for these planets. The stellar wind might also be a contributing factor for WASP-52b, but it is unlikely to explain the nondetection for WASP-80b. Composition and magnetic fields could both play a role in explaining the observations. We propose some observational tests in Section 6 that could help to distinguish between these competing hypotheses.

At the beginning of this work, we hypothesized that if the Neptune desert is indeed cleared by mass loss, then planets at its boundaries should be marginally stable against photoevaporation. We can now test this hypotheses using our measurements and some simple calculations. First, we inspect our observationally constrained mass-loss rates to evaluate our sample's present-day stability to atmospheric erosion. We expect that close-in giant planets will be destroyed while their host stars are on the main sequence (Hamer & Schlaufman 2019), so a reasonable timescale for atmospheric stability is ${M}_{\mathrm{env}}/\dot{M}\sim 4\,\mathrm{Gyr}$. In Table 5, all of the planets in our sample near the upper edge of the desert have predicted atmospheric lifetimes that are much greater than this timescale.

For HAT-P-26b near the lower edge of the desert, our result is inconclusive due to the $\dot{M}-{T}_{0}$ degeneracy: if the outflow is hot (T₀ ≳ 10,000 K), the planet's predicted atmospheric lifetime may be as short as 3.3 Gyr. This is substantially shorter than the envelope mass-loss timescale for HAT-P-11b: using the mass-loss rate from (dos Santos et al. 2022) and the envelope mass fraction of 15.1% from (Lopez & Fortney 2014), the mass-loss timescale for HAT-P-11b is ∼30 Gyr. Importantly, the $\dot{M}-T$ degeneracy for HAT-P-11b is broken by the precise line shape measurement (Allart et al. 2018; dos Santos et al. 2022) and also by energetic arguments (Vissapragada et al. 2022). Both methods suggest the outflow temperature for HAT-P-11b is T₀ ≲ 8000 K, and if this is the case for HAT-P-26b as well, it would have a longer inferred envelope loss timescale. In the future, high-resolution spectroscopy of HAT-P-26b could resolve this degeneracy by measuring the helium line shape and corresponding outflow temperature.

Although the present-day mass-loss rates of the planets in our sample appear to be low, they experienced stronger XUV irradiation (and thus suffered higher mass-loss rates) in the past. We can reconstruct their past mass-loss histories by integrating the energy-limited mass-loss rate back in time (see, e.g., Lecavelier Des Etangs 2007; Jackson et al. 2012; Mordasini 2020). We first start with the expression:

$$\begin{eqnarray}&&\dot{M}=-\displaystyle \frac{\varepsilon {R}_{{\rm{p}}}^{3}{L}_{\mathrm{XUV}}}{4{a}^{2}{KGM}},\end{eqnarray} \tag{ 4 }$$

which is the typical expression for the energy-limited mass-loss rate (Equation (2)) with the scaled luminosity L_XUV/(4π a²) substituted for the flux F_XUV. For this estimate, we assume that the planetary radius is approximately invariant to small changes in the mass. We can do this because the mass–radius relation is nearly flat for the giant planets at the upper edge of the Neptune desert (e.g., Chen & Aipping 2017; Owen & Lai 2018; Thorngren et al. 2019). This is because giant planets at the upper edge of the desert are well approximated by P ∝ ρ² polytropes for which the radius does not depend on mass (Stevenson 1982). Additionally, the variations in K with mass are also small and can be ignored. The integral then reads:

$$\begin{eqnarray}&&{\int }_{{M}_{0}}^{{M}_{{\rm{p}}}}{MdM}={\int }_{0}^{t}-\displaystyle \frac{\varepsilon {R}_{{\rm{p}}}^{3}{L}_{\mathrm{XUV}}}{4{a}^{2}{KG}}{dt},\end{eqnarray} \tag{ 5 }$$

where we have introduced the initial planetary mass M₀ and the planet's current age t. To allow for maximum possible mass loss, we assume that the planet has maintained the same mass-loss efficiency for its whole life, even though the efficiency at early times is likely much lower in the recombination limit (Murray-Clay et al. 2009; Owen & Alvarez 2016). Additionally, we assume the planet has remained at its current semimajor axis for its whole life, i.e., that it did not migrate from a more distant location where it experienced lower instellation. Then, defining the integrated stellar XUV luminosity ${E}_{\mathrm{XUV}}={\int }_{0}^{t}{L}_{\mathrm{XUV}}{dt}$, we obtain the equation:

$$\begin{eqnarray}&&{M}_{{\rm{p}}}^{2}-{M}_{0}^{2}=-\displaystyle \frac{\varepsilon {R}_{{\rm{p}}}^{3}{E}_{\mathrm{XUV}}}{2{a}^{2}{KG}}.\end{eqnarray} \tag{ 6 }$$

Finally, writing the current mass of the planet as a fraction f < 1 of the initial mass (such that M₀ = M_p/f), we obtain:

$$\begin{eqnarray}&&{M}_{{\rm{p}}}=\displaystyle \frac{1}{a}\sqrt{\left(-\displaystyle \frac{1}{1-\tfrac{1}{{f}^{2}}}\right)\displaystyle \frac{\varepsilon {R}_{{\rm{p}}}^{3}{E}_{\mathrm{XUV}}}{2{KG}}}.\end{eqnarray} \tag{ 7 }$$

This is a line in the M_p − a plane above which planets have lost less than a fraction f of their initial mass, and below which they have lost more. We hereafter take f = 0.5—a planet losing half of its initial mass—as a metric for marginal stability. To calculate this boundary, we require a prescription for the planetary radii and the integrated stellar XUV flux in addition to our empirically constrained distribution for ε. For the planetary radii, we use the empirical radius–temperature relation from Equation (1) of Owen & Lai (2018) for a typical K dwarf temperature T_⋆ = 4700 K and radius R⋆ = 0.8R_☉. The spread in the relation is incorporated into our final uncertainties. We note that this relation is only valid for M_p ≳ 0.2M_J; below this mass, the planetary radius will change as the planet undergoes photoevaporation. For the integrated XUV flux, a reasonable upper limit for a K dwarf is E_XUV ∼ 10⁴⁶ erg using the fast rotator track from Johnstone et al. (2021) as a guide (see their Figure 18).

In Figure 5, we draw the f = 0.5 boundary on the M_p − a plane for an R = R_J and K = 0.5 planet using our empirical efficiency from Section 5.1. For comparison, we also draw the f = 0.9 boundary (corresponding to 10% mass loss over the planetary lifetime). We conclude that planets at the upper edge of the Neptune desert have lost less than 10% of their initial masses to photoevaporation. This empirically benchmarked result is in good agreement with previous theoretical calculations by Ionov et al. (2018) and Owen & Lai (2018). HAT-P-26b, near the lower edge of the Neptune desert, is a potential exception—it could have lost ∼50% of its initial mass if its outflow was energy-limited for the entire planetary lifetime. Although planets like HAT-P-26b along the lower edge of the desert could still be substantially affected by mass loss, we caution that our simplifying assumption that the radius is insensitive to small changes in envelope mass will not hold true for these smaller planets, and a more detailed framework taking into account the core-envelope structure is needed for understanding this part of the population (e.g., Owen & Lai 2018; Mordasini 2020).

Zoom In Zoom Out Reset image size

Throughout this section, we made a number of assumptions to try to estimate the maximum possible mass-loss endured by planets in and near the upper edge of the Neptune desert. In reality, not all planetary host stars will have XUV luminosities as large as those of the most rapidly rotating stellar models from Johnstone et al. (2021). We also assumed an energy-limited outflow for the entire planetary lifetime; this is certainly an overestimate at early times, and the efficiency will be much lower during the period of maximum XUV irradiation (Murray-Clay et al. 2009; Owen & Alvarez 2016). Despite this, we still found that the boundary for marginal stability is located well below the actual edge of the desert. We conclude that photoevaporation cannot explain the upper boundary of the Neptune desert.