The Upper Edge of the Neptune Desert Is Stable Against Photoevaporation

,


INTRODUCTION
The transit and radial velocity techniques have revealed masses and radii for thousands of planets spanning a wide range of orbital periods.These surveys indicate that relatively few sub-Jovian planets reside on close-in orbits (P 5 days).This "Neptune desert" is a robust feature of the planetary census, which is relatively complete in this regime (Szabó & Kiss 2011;Beaugé & Nesvorný 2013;Lundkvist et al. 2016).Mazeh et al. (2016) defined boundaries for the Neptune desert by searching for curves that maximize the contrast between regions of the mass-period plane with and with-out planets.We plot their definition for the Neptune desert in Figure 1 along with the planetary population observed to date.
Given its close-in location, it is natural to wonder whether the desert might have been created by atmospheric mass loss (e.g.Youdin 2011;Beaugé & Nesvorný 2013;Kurokawa & Nakamoto 2014;Ionov et al. 2018;Owen & Lai 2018).Planets with short orbital periods are subjected to intense amounts of high-energy radiation that can drive atmospheric outflows (Owen 2019).The resulting atmospheric mass-loss rates can be large, especially for planets orbiting young stars, which emit proportionally more of their energy in high-energy photons.If close-in gas-giant planets formed in situ (e.g.Batygin et al. 2016) or arrived at their present-day locations relatively early (∼10 Myr) via disk migration (e.g.Ida & Lin 2008), they might have experienced a pe- Planets orbiting stars with 4000 K< T eff < 5400 K are outlined in red, and the seven targets constituting our sample are outlined in blue.The point color denotes the host star's J magnitude, and the point size is proportional to logarithm of the predicted SNR in the WIRC helium bandpass from Section 2.1.
riod of strong photoevaporation (e.g.Murray-Clay et al. 2009).Hot Neptunes typically have lower gravitational potentials than their Jovian counterparts, making them more susceptible to photoevaporation and/or Roche lobe overflow (Kurokawa & Nakamoto 2014;Valsecchi et al. 2015;Koskinen et al. 2022).This means that Jupiters could survive at orbital separations where most Neptunes would be destroyed.
If this population of close-in giant planets instead arrived at their present locations relatively late (∼Gyr) via dynamical interactions with other bodies in their systems (Rasio & Ford 1996), we would expect photoevaporation to play a relatively minor role in their long-term evolution.In this scenario, planets undergoing high-eccentricity migration circularized onto orbits with final semimajor axes of a ≈ 2r peri , where r peri is the pericenter distance when the eccentricity e ∼ 1. Neptune-mass planets with small present-day semimajor axes would have had pericenter passages inside the stellar tidal disruption radius, whereas Jupiter-mass planets would survive at comparable separations (Ford & Rasio 2006;Guillochon et al. 2011;Matsakos & Königl 2016;Owen & Lai 2018).
Our ability to explain the Neptune desert is thus intimately linked to our understanding of the origins of close-in giant planets (Dawson & Johnson 2018).Unfortunately, previous studies have found it difficult to differentiate between atmospheric escape and higheccentricity migration as explanations for this feature.
Depending on their assumptions about the energylimited mass-loss efficiency, stellar X-ray and ultraviolet (XUV) flux, initial star-planet separation, and pressure at the XUV photobase, some published studies suggest that sub-Jovian planets near the upper edge of the desert are stable against evaporation, while other suggest that they experience runaway envelope escape (Kurokawa & Nakamoto 2014;Ionov et al. 2018;Owen & Lai 2018;Rao et al. 2021).High-eccentricity migration models also have their own uncertainties, most notably the assumed tidal quality factor, the initial star-planet separation, and the eccentricity excitation mechanism (Matsakos & Königl 2016;Owen & Lai 2018).Previous observational studies have used the orbital properties and multiplicity of hot Jupiter systems to constrain the fraction of hot Jupiters that might have formed via high-eccentricity migration (e.g.Winn et al. 2010;Knutson et al. 2014;Dawson et al. 2015;Ngo et al. 2016;Rice et al. 2022).Although it appears unlikely that all hot Jupiters formed this way, such planets almost certainly comprise a subset of the hot Jupiter population (see the review by Dawson & Johnson 2018).
Observations of atmospheric escape can provide a complementary basis to differentiate between these scenarios for the origin of the Neptune desert.If atmospheric escape is indeed the mechanism that clears the Neptune desert, planets at its edge would be marginally stable against photoevaporation, potentially exhibiting large outflows today.Recent studies have demonstrated that the metastable helium triplet can be used to detect and characterize planetary outflows for close-in transiting planets (e.g.Spake et al. 2018;Oklopčić & Hirata 2018).Unlike the Lyman-α transit depth -which is controlled not by the planetary mass-loss rate, but by the distance neutral hydrogen atoms can travel before they are photoionized (Owen et al. 2021) -the metastable helium absorption amplitude probes low-velocity thermospheric gas near the wind-launching radius.This means that it can be readily modeled using a one-dimensional isothermal Parker wind to infer mass-loss rates as low as 10 9 − 10 10 g s −1 (e.g.Mansfield et al. 2018;Oklopčić & Hirata 2018;Lampón et al. 2020;Vissapragada et al. 2020b).For comparison, the required average mass-loss rate to clear a hypothetical M J planet from the desert in a few hundred Myr is of order 10 14 g s −1 .The stellar XUV luminosity (on which the planetary mass-loss rate is linearly dependent in the energy limit) decreases by about 2 orders of magnitude over the stellar lifetime (e.g.Johnstone et al. 2021).Therefore, if planets at the upper edge of the desert are only marginally stable against photoevaporation, they would be expected to exhibit relatively strong outflow signatures even given the reduced XUV fluxes of their middle-aged host stars.Although uncertainties in the assumed incident XUV spectrum, outflow temperature, and outflow composition limit our ability to predict the metastable helium population in the Parker wind model, we can still constrain presentday mass-loss rates well enough for planetary evolution studies.This is especially true when we incorporate additional information like line shapes and energetics (dos Santos et al. 2022;Vissapragada et al. 2022).
We recently installed a narrowband filter centered on the metastable helium line (Vissapragada et al. 2020b) in the Wide-field InfraRed Camera (WIRC; Wilson et al. 2003) at the Hale 200-inch Telescope at Palomar Observatory.We commissioned this new observing mode by confirming the outflow of WASP-69b, finding an absorption signal commensurate with spectroscopically resolved observations by Nortmann et al. (2018).As part of this same study we also observed a transit of WASP-52b but did not detect a strong signal.We have since used this same observing mode to detect an atmospheric outflow from HAT-P-18b (Paragas et al. 2021) and found a tentative evidence for an outflow from the young planet V1298 Tau d (Vissapragada et al. 2021).
From 2019-2021, we used our narrowband helium photometer to survey the metastable helium line in a seven-planet sample near the upper edge of the Neptune desert.There are currently ten published helium detections in the literature: WASP-107b (Spake et al. 2018;Allart et al. 2019;Kirk et al. 2020;Spake et al. 2021), HAT-P-11b (Allart et al. 2018;Mansfield et al. 2018), WASP-69b (Nortmann et al. 2018;Vissapragada et al. 2020b), HD 209458b (Alonso-Floriano et al. 2019), HD 189733b (Salz et al. 2018;Guilluy et al. 2020), GJ 3470b (Palle et al. 2020;Ninan et al. 2020), HAT-P-18b (Paragas et al. 2021), HAT-P-32b (Czesla et al. 2022), TOI-560.01 (Zhang et al. 2022), and GJ 1214b(Orell-Miquel et al. 2022).These observations, along with many non-detections and tentative detections in the literature, span wide ranges in planetary mass and radius, stellar spectral type and activity level, system age, observing methodology, and data reduction technique.This heterogeneity can obfuscate underlying trends in the data, because the measured helium absorption signal depends on a variety of factors that all vary across the sample.We therefore focus on a uniformly selected sample of helium observations obtained with Palomar to search for trends in measured mass-loss rates for planets at the upper edge of the Neptune desert.
In this work we present the results of our survey, including a re-analysis of the published transit observations from Vissapragada et al. (2020b) and Paragas et al. (2021).We also present new observations of eight additional transits for five planets: WASP-52b, WASP-80b, WASP-177b, HAT-P-26b, and NGTS-5b.In Section 2, we describe our sample selection methodology as well as our observation and data reduction procedures.In Section 3, we present helium-band light curves for each of our seven targets.In Section 4, we model each of our light curves using a one-dimensional isothermal Parker wind model to constrain the planetary mass-loss rates.We discuss the implications of our results for the population-level picture of atmospheric escape in Section 5, and summarize our conclusions in Section 6.

Sample Selection
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 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 (SNR).We estimated the SNR 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 resolv- 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  et al. (2019).The log(R HK ) value for NGTS-5b was calculated using an S-value from P. Eigmüller (private communication).
ing 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,b) to the 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;Allart et al. 2019;Guilluy et al. 2020;Kirk et al. 2020;Spake et al. 2021) 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 SNR > 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.

Palomar/WIRC Observations
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 a full-width at half-maximum (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;Paragas et al. 2021;Vissapragada 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 Jun 05.Weather conditions were poor that night, and we chose to slightly defocus the telescope to 1. 2 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 center 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).

Data Reduction
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 timevarying telluric water absorption.This is necessary because comparison stars occupy different parts of the field and thus experience different amounts of OH emission Note-In the column headers, nexp refers to the total number of exposures, texp is the exposure time, ncomp is the number of comparison stars, n dither is the total number of dither frames used in constructing the background, r phot is the radius of the circular aperture used in the photometric extraction, σ is the rms scatter of the residuals to our final light-curve fit, and σ/σ phot is the ratio of the rms scatter to the photon noise limit.
The two σ and σ/σ phot values for the WASP-80b observation indicate the noise before and after the instrument was reset.All dates and times are UT.
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 (0. 75-3.75) 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 datasets except for our observation of WASP-80b on UT 2020 Jul 09, 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.

Light-curve Modeling
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. 2021a,b).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 0 , the scaled semi-major axis a/R , the impact parameter b, and the quadratic limb darkening coefficients (u 1 , u 2 ).
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 our data to constrain the limb darkening themselves.
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.MacLeod & Oklopčić 2022; Wang & Dai 2021a,b).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 non-zero orbital eccentricities (Hartman et al. 2011a,b;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 BJD times as a decorrelation vector, which functions as a linear baseline.Each weight was assigned a uniform prior of U(−2, 2).We list the rest of the priors for the transit light-curve modeling in Table 3.
After defining the model, we sampled the posterior distributions for each model parameter using the No U-Turn Sampler (NUTS Hoffman & Gelman 2011) implemented in pymc3 (Salvatier et al. 2016).In each lightcurve fit, we ran four chains for 1,500 tuning steps each before taking 2,500 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 R 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 repeat-ing the fit with every combination of these three parameters and calculated the Bayesian Information Criterion (BIC; Schwarz 1978) for each fit: where k is the number of fit parameters, n is the number of data points, and 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 (MAP) 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 dataset, 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 log(σ extra ), which quantifies the discrepancy between the photon noise and the true variance in the data (i.e.σ 2 = σ 2 photon + σ 2 extra ).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 lightcurve 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 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 order ∼ 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 nm-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.

WASP-69b
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 Aug 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. We find the excess depth at mid-transit to be 0.512 +0.049  −0.048 % 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.

WASP-52b
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 +300 −200 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;Alam et al. 2018;Chen et al. 2020).
We first searched for helium in the atmosphere of this planet on UT 2019 Sep 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 Aug 04.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.14  −0.13 % on the first night, with a 95th-percentile upper limit of 0.44%, and 1.02 +0.39  −0.41 % on the second night with a 95th-percentile upper limit of 1.68%.When fitting the datasets 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 3.

HAT-P-18b
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. 2011a).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  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.8nm-1141.6 nm) is narrower and closer to the helium filter bandpass than the TESS bandpass (600 nm-1000 nm), and the depth in this bin is more precise than the TESS depth (∼ 100 ppm precision for WFC3 vs. ∼ 500 ppm for TESS ).Additionally, in the previous work we reported (R p /R ) 2 WIRC − (R p /R ) 2 comparison 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 dataset independently, we ob-tained (R p /R ) 2 values of 2.14 +0.23  −0.24 % and 2.56 +0.17 −0.16 % 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 p /R ) 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 datasets jointly and show the results in Figure 4.In the joint fit, we find (R p /R ) 2 = 2.46 ± 0.15%, which is about 1σ larger than the result from Paragas et al. (2021)

WASP-80b
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(Triaud et al. , 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;Parviainen et al. 2018;Kirk et al. 2018).Multiple groups running self-consistent wind-launching simulations have reported WASP-80b to be an excellent candidate for massloss observations (Salz et al. 2015(Salz et al. , 2016a;;Caldiroli et al. 2022).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 Jul 09.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 δ mid = 0.02 +0.20  −0.20 with a 95th-percentile upper limit of 0.35% excess absorption in our bandpass.

WASP-177b
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 Oct 2. TESS also observed this planet (TOI-4521.01)at 2 minute cadence in Sector 42 from UT 2021 Aug -2021 Sep.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 semi-major 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 6.The fitted TESS lightcurve is given in Appendix D. Using the TESS transit depth as our baseline, we constrain the excess depth to be 0.53 +0.23  −0.28 %, a non-detection with a 95th-percentile upper limit of 0.90%.

HAT-P-26b
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. 2011b).The Neptune-mass planet has a well-constrained heavy element abundance thanks to the presence of strong water bands in the red-optical and near-infrared (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)  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 Feb 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 p /R = 0.034 +0.028 −0.023 .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.18  −0.17 % 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 7.We obtained a mid-transit excess depth of 0.31 +0.10 −0.10 , evidence for helium absorption at 3.1σ confidence.

NGTS-5b
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 dataset, and the issue stemmed purely from our approach to modeling the systematics.We therefore kept this dataset 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.99  −0.95 % and 0.94 +0.55 −0.54 % for the individual fits to the first and second nights, respectively.We proceeded to fit the datasets jointly, and the results are shown in Figure 8.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.48  −0.46 %, a tentative detection at 2.2σ confidence.

MASS-LOSS MODELING
We can convert the excess depths from Section 3 into constraints on planetary mass-loss rates by modeling the outflows with one-dimensional isothermal Parker winds (Oklopčić & Hirata 2018;Lampón et al. 2020).In this section, we use the p-winds code (dos Santos et al. 2022) to model our observations.Given a mass-loss rate assuming irradiation conditions at the substellar point Ṁsubstellar , a thermosphere temperature T 0 , and a hydrogen fraction for the outflow f H , we first set the density and velocity profiles for the wind.Then, we use a high-energy stellar spectrum to calculate the ionization structure and level populations throughout the outflow.
The level populations are used to compute the expected wavelength-dependent absorption signal at high resolving power, which we then convolve with our filter bandpass to compare to the measurements in Section 3. The free parameters in this model are the mass-loss rate, thermosphere temperature, the hydrogen fraction, and the high-energy spectrum.We compute the model over a 50 × 50 grid of mass-loss rates Ṁsubstellar = 10 10 g s −1 − 10 13 g s −1 and thermosphere temperatures T 0 = 5000 K−15000 K.The grid was uniformly spaced in log( Ṁsubstellar ) and T 0 .We set the hydrogen fraction for all outflows to 0.9, corresponding to a 90-10 hydrogen-helium outflow.Finally, for the high-energy spectra we select proxy stars from the v2.2 panchromatic SEDs from the MUSCLES survey (at 1 Å binning; France et al. 2016;Loyd et al. 2016;Youngblood et al. 2016).We identify three K stars in MUSCLES that appear to be good proxies (i.e. they are the most similar in both spectral type and activity index log(R HK )) for our targets: HD 85512, HD 40307, and Eridani.We model the late K stars (T 4800 K) in our sample using HD 85512's spectrum, and use HD 40307's spectrum for the early K stars.We use Eridani's spectrum to model the most active stars in our sample, includ- We list the spectrum used for each star in Table 5.We then scale each MUSCLES spectrum by the stellar radius and planetary orbital distance to obtain the highenergy spectrum incident at the top of the planet's atmosphere.We integrate each spectrum up to the hydrogen photoionization threshold 91.2 nm to obtain an estimate of the top-of-atmosphere XUV flux F XUV for each planet, which are included in Table 5.
After running an initial grid of models, we noticed that many of our winds had sonic points located outside the Roche radius.As discussed in Murray-Clay et Note-FXUV corresponds to the XUV flux at the planet.We quote the 95% confidence intervals for Ṁ , ε, and Menv/ Ṁ as the marginalized Ṁ distributions are highly non-Gaussian (see Figure 10).Upper limits correspond to 95% on the corresponding distributions.
al. ( 2009), this happens when the tidal gravity term in the momentum equation for the wind is not taken into account.We show in Appendix E that the tidal gravity term appreciably alters the velocity and density profiles for the planets in our sample, and we update the p-winds code (in version 1.3) to account for the effects of the stellar gravity.We then re-calculate our model grids for each planet and use these grids to determine the range of mass-loss rates and thermosphere temperatures that match our observations.Next, we assess which combinations of Ṁsubstellar and T 0 are energetically self-consistent using the methodology outlined in Vissapragada et al. (2022).If photoionization is the only heat source, certain combinations do not satisfy energy balance in the isothermal Parker wind, and we omit them from our final results.
Finally, we divide the substellar-point mass-loss rates by a factor of 4 to obtain estimates for the overall massloss rate Ṁ .This correction factor accounts for the fact that the planet is only irradiated over π steradians rather than the 4π steradians assumed by our 1D models, so the true mass-loss rates will be about a factor of 4 smaller (Salz et al. 2016a;Vissapragada et al. 2022).Other prescriptions to correct for 2D effects exist (e.g.Odert et al. 2020), but they all agree well for outflows from low-gravity planets (like the ones in our sample) that are dominated by adiabatic cooling (Caldiroli et al. 2021).Performing the radiative transfer using the 1D model rather than a 2D or 3D model also introduces error, but comparisons between 1D and 3D models suggest that this is a minor error term (relative to the uncertainty from the Ṁ − T 0 degeneracy) except in the presence of strong stellar winds (MacLeod & Oklopčić 2022; Wang & Dai 2021a,b).
We show the resulting mass-loss constraints for each planet in Figure 9.We used a rejection sampling scheme with this grid of results to generate a list of 100,000 samples for each planet corresponding to the underlying joint distribution of mass-loss rates and thermosphere temperatures.We began by uniformly sampling temperatures from T 0 ∼ U(5000, 150000) and log( Ṁsubstellar ) ∼ U(10, 13).For each sample, we found the closest point on the grid in Figure 9 with corresponding xσ discrepancy between the model and the data.We accepted the sample only if y > erf(x/ √ 2), where y is a random draw from U(0, 1).
With the final list of 100,000 samples in hand for each planet, we marginalize over the thermosphere temperature to obtain distributions for the mass-loss rates, which are summarized in Table 5.We evaluate the significance of these mass-loss rates for the expected atmospheric lifetime of each planet by calculating the envelope loss timescale M env / Ṁ .For HAT-P-26b we adopt an envelope mass fraction of 31.7% (Lopez & Fortney 2014).For the other higher-mass gas giants in the sample, we take M env ∼ M .We give the resulting envelope loss timescales in Table 5.

Benchmarking Mass-Loss Models
Self-consistent, one-dimensional mass-loss models like the Pluto-CLOUDY Interface (TPCI; Salz et al. 2015Salz et al. , 2016a) ) and the ATmospheric EScape code (ATES; Caldiroli et al. 2021Caldiroli et al. , 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 sevenplanet 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-toplanet 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): 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 : Planets in the energy-limited regime should therefore exhibit a linear relation between the K-corrected massloss rate and the ratio F XUV /ρ p (Caldiroli et al. 2022).
As the incident XUV flux increases, however, selfconsistent models predict that the outflow begins to lose substantial energy to radiative cooling, leading to a sub-linear dependence on F XUV /ρ p beyond a threshold of F XUV /ρ p ∼ 10 4 erg cm g −1 s −1 (Murray- Clay et al. 2009;Salz et al. 2016a,b;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 −11 erg) exceeds the particle binding energy, which for a hydrogen atom (m H ∼ 10 −24 g) occurs roughly at a gravitational potential of φ p ∼ ∆E/m H ∼ 10 13 erg/g (for a more thorough explanation, see Caldiroli et al. 2022).More sophisti-cated numerical experiments place this threshold closer to φ p 10 13.2 erg/g (Salz et al. 2016a;Caldiroli et al. 2022).
In Figure 10 we plot our marginalized and K-corrected Ṁ 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, comprised of 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.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 Ṁ 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 ε = 0.41 +0.16  −0.13 .We overplot this average ε model in Figure 10.
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 massloss rate and thermosphere temperature, are too large to differentiate between this model and a constant-ε model.Resolving the line profiles of these planets with highresolution 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;Allart et al. 2019;Kirk et al. 2020;Spake et al. 2021) is well-matched by the three-dimensional hydrodynamical models of Wang & Dai (2021b), who report a mass-loss rate of approximately 2 × 10 11 g s −1 .Using the scaled MUSCLES spectrum of the similarly active star Eridani (log(R HK ) = −4.44;Noyes et al. 1984) to calculate the XUV flux for WASP-107b (log(R HK ) = −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 MUS-CLES 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.7 −0.5 × 10 10 g s −1 (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.052 −0.037 , 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. (2021a).For their 90/10 H/He composition model, they found a mass-loss rate at T = 12000 K of about 4 × 10 9 g s −1 , with significant spread due to the Ṁ − 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. 2016a;Caldiroli et al. 2021Caldiroli et al. , 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 (log (R HK ) = −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 predicted by 1D hydrodynamical models, which give mass-loss rates of about 10 9.5−10.0g s −1 for this planet (Salz et al. 2016a;Caldiroli et al. 2021Caldiroli et al. , 2022)).Lampón et al. (2021a) resolve 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.Figure 10.K-corrected mass-loss rate Ṁ as a function of FXUV/ρp for all the planets in our sample, following Caldiroli et al. (2022).The shaded violins indicate the marginalized distribution of K-corrected mass-loss rates with the median and 95% confidence interval for the distribution given by the black points and error bars, respectively.The triangle denotes the 95% upper limit for WASP-80.The orange curve indicates the average inferred mass-loss efficiency ε, with the 1σ uncertainty given by the shaded orange region.The blue curve indicates predictions from the ATES code (Caldiroli et al. 2021(Caldiroli et al. , 2022)).

Low Mass-Loss
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 log(R HK ) 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, strongerthan-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 ten 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.064 −0.048 .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 ten 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 X = 4.85 +0.12 −0.23 ×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 2 /s.Two other papers independently derive this star's X-ray flux and find values that are approximately a factor of three smaller (Monsch et al. 2019) and a factor of two larger (Foster et al. 2022, who obtain a flux close to our adopted value in Table 5) than that obtained by Fossati et al. (2022).Despite these methodological differences, factorof-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 X = 3.1 +1.2 −1.0 × 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 2 /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 ≈3.6 km/s with a photometric rotation period of 16.4 days, Hébrard et al. 2013) and WASP-80 (v sin i ≈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 log R HK values that also suggest ages as young as a few hundred Myr (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 three-dimensional hydrodynamical model (Shaikhislamov et al. 2021;Khodachenko 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 13 g s −1 ∼ 8 Ṁ , where the solar mass-loss rate is Ṁ ≈ 2 × 10 −14 M yr −1 (e.g.Wood et al. 2005), but at early times the wind can be even stronger.Wood et al. (2005) used astrospheric Lyman-α absorption to infer mass-loss rates of 30 Ṁ and 100 Ṁ 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 massloss rate for WASP-80.We obtain a mass-loss rate of ∼ 4 Ṁ , well within the range modeled by Fossati et al. (2022).Using the same relation, we obtain a massloss rate of ∼ 35 Ṁ 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. 2020Lampón et al. , 2021a) 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 Lyman-α and metastable helium lines to Parker wind models (Lampón et al. 2020(Lampón et al. , 2021a,b) ,b) or 3D hydrodynamical models (Shaikhislamov et al. 2021).However, Lyman-α 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 Lyman-α 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 log(R HK ) 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 in a few systems (SPI; Cauley et al. 2019).Other observations suggest that these planets may have more modest field strengths (1-10 G), including ultraviolet observations of Lyman-α 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 ∼ 1/ √ 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 G 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 windlaunching 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 non-detection 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 non-detection 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.

The Upper Neptune Desert Is Stable Against Mass Loss
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 observationallyconstrained 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 env / Ṁ ∼ 4 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 Ṁ − T 0 degeneracy: if the outflow is hot (T 0 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 Ṁ − 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 0 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, highresolution 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 energylimited 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: which is the typical expression for the energy-limited mass-loss rate (Equation 2) with the scaled luminosity L XUV /(4πa 2 ) 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 & Kipping 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 ∝ ρ 2 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: where we have introduced the initial planetary mass M 0 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 XUV = t 0 L XUV dt, we obtain the equation: Finally, writing the current mass of the planet as a fraction f < 1 of the initial mass (such that M 0 = M p /f ), we obtain: 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 46 erg using the fast rotator track from Johnstone et al. (2021) as a guide (see their Figure 18).In Figure 11, we draw the f = 0.5 boundary on the M p − a plane for a 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).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.WASP-52b and NGTS-5b; and did not detect (< 2σ) absorption from WASP-80b and WASP-177b.
2. When interpreting these signals with a onedimensional Parker wind model, we found that the six planets in our sample near the upper edge of the desert (WASP-69b, WASP-52b, HAT-P-18b, WASP-80b, WASP-177b, and NGTS-5b) have atmospheric lifetimes M env / Ṁ much larger than 10 Gyr.Our result for HAT-P-26b (near the lower edge of the desert) is inconclusive due to the Ṁ − T 0 degeneracy.
3. We compared our empirically measured mass-loss rates to predictions from the one-dimensional, selfconsistent hydrodynamics code ATES (Caldiroli et al. 2021(Caldiroli et al. , 2022)).We found that five planets were in good agreement with the ATES predictions: HAT-P-18b, WASP-69b, HAT-P-26b, NGTS-5b, and WASP-177b.The mass-loss rates for these planets are all consistent with a mean outflow efficiency of ε = 0.41 +0.16 −0.13 .4. We found that WASP-52b and WASP-80b have much lower inferred mass-loss rates than predicted by these models.The measured helium absorption signals for these two planets may be affected by stellar wind confinement, helium depletion, the inverse first ionization potential effect, or magnetic field confinement, although stellar wind confinement appears unlikely for WASP-80b.
5. Our empirically measured mass-loss efficiencies are too small for photoevaporation to sculpt the population of giant planets at the upper edge of the Neptune desert, in agreement with previous work by Owen & Lai (2018) and Ionov et al. (2018).
We conclude that another mechanism besides photoevaporation must be responsible for carving the upper edge of the Neptune desert.Candidate explanations include high-eccentricity migration (Matsakos & Königl 2016;Owen & Lai 2018) and in situ formation near the magnetospheric truncation radius of the natal protoplanetary disk (Bailey & Batygin 2018).Although our observations cannot distinguish between these two scenarios, future studies of the functional dependence of the desert's upper boundary on stellar properties may yield more insight into the problem.This will require a large, statistically uniform sample of hot Jupiters, which will soon be provided by TESS (Yee et al. 2021(Yee et al. , 2022)).Additional observations aimed at constraining the presence of long-period planetary companions, for instance with Gaia astrometry, could also help to differentiate between these hypotheses.
Our survey also revealed that 1D mass loss models overpredict the metastable helium signatures from WASP-52b and WASP-80b, in agreement with previous work by Fossati et al. (2022).Although it cannot explain WASP-80b, our estimates suggest that a strong stellar wind might confine the outflow of WASP-52b, similar to what has been proposed for WASP-107b (Spake et al. 2021;Wang & Dai 2021b) and TOI-560.01 (Zhang et al. 2022).If the helium line for WASP-52b can be detected with more sensitive observations at higher resolving power, we would expect to see a characteristic blueshifted line profile and/or an extended tail of postegress absorption in this scenario.If the outflow is instead confined by magnetic fields, the line profile would not be strongly Doppler shifted.We could also test the magnetic confinement explanation for both planets by searching for evidence of star-planet interactions.Cauley et al. (2019) recently demonstrated that close-in planets with strong magnetic fields can affect the stellar chromospheric emission.These magnetic star-planet interactions (SPI) are detected in the Ca II K line, where variations in stellar emission occur on the planetary orbital timescale.It is also possible that the stellar XUV templates we have used for these stars are inadequate as the coronal iron abundances may differ between WASP-80/WASP-52 and the template star Eridani (Poppenhaeger 2022); this motivates future X-ray monitoring of these sources.Lastly, observations of the Hα absorption signal from both planets could also provide constraints on the outflow composition (Czesla et al. 2022).Hα absorption has already been detected and used to constrain composition for WASP-52b (Chen et al. 2020;Yan et al. 2022); similar constraints should also be obtained for WASP-80b.
During the review process for this paper, Kirk et al. (2022) announced a strong detection of helium absorption in WASP-52b and tentative evidence for helium absorption in WASP-177b, both with Keck/NIRSPEC.For WASP-177b, their result is consistent with our upper limit, but is more constraining on the mass-loss rate: they obtain a 3σ upper limit of 7.9 × 10 10 g s −1 on the mass-loss rate.The Kirk et al. (2022) result matches the posterior distribution shown in Figure 10 quite well after accounting for the K = 0.52 correction factor.For WASP-52b, Kirk et al. (2022) observe a slightly stronger signal than ours: 0.66±0.14%(in a 0.635 nm bin) compared to our result of 0.29±0.13%,and they obtain a mass-loss rate of 1.2 ± 0.5 × 10 11 g s −1 .This is still an order of magnitude smaller than the prediction from ATES in Figure 10 (K Ṁ = 4.7 × 10 11 g s −1 ) after accounting for the K = 0.43 correction factor.Both Kirk et al. (2022) and the subsequent modeling effort by Yan et al. (2022) find similar inferred outflow rates that are not very sensitive to the assumed composition.Additionally Kirk et al. (2022) find that the outflow is not strongly confined by stellar winds as they do not observe a Doppler-shifted line profile or a tail.Of the factors we considered for the low inferred mass-loss rate for WASP-52b compared to model predictions, this leaves the inverse first ionization potential effect and/or confinement by magnetic fields as the most likely explanations.
Intriguingly, the TESS survey is also beginning to discover planets within the Neptune desert, including TOI-849b (Armstrong et al. 2020) and LTT 9779b (Jenkins et al. 2020).Dai et al. (2021) show that these desertdwellers tend to orbit metal-rich stars and lack planetary companions, suggesting that planets within the desert are more similar to gas giant planets than the rocky, ultra-short period planets below the desert.However, our result confirm that these desert-dwellers are unlikely to be the photoevaporated cores of more massive gas giant planets.In the high-eccentricity migration scenario, these unique planets could be the results of partial tidal disruption as they circularized onto their current orbits (Faber et al. 2005;Guillochon et al. 2011).Alternatively, Pezzotti et al. (2021) suggest that if TOI-849b had a large initial mass and radius, it may instead have lost its envelope to Roche-lobe overflow (RLO).This process can lead to mass-loss rates far exceeding those from photoevaporation, and is predicted to be consequential for planets on exceptionally short-period orbits (Jackson et al. 2017).While it likely plays a role for the closest-in Neptune desert planets, the RLO massloss rate drops off precipitously with orbital distance, so it cannot explain the entirety of the desert (Koskinen et al. 2022).
Although these high density desert-dwelling Neptunes likely have very low present-day mass-loss rates, the picture for lower density Neptunes at the lower edge of the desert is somewhat murkier.Previous modeling studies (Kurokawa & Nakamoto 2014;Lundkvist et al. 2016;Ionov et al. 2018;Owen & Lai 2018) concluded that mass loss can be significant for planets near the lower edge of the desert.Our result for HAT-P-26b appears to be consistent with this prediction.TESS has identified a large sample of new planets near the lower edge of the desert that are favorable targets for mass-loss measurements, making this a promising area for future investigation.If planetary evolution is dominated by photoevaporation in this mass regime, we also expect the lower boundary of the desert to shift to smaller insolations for less luminous late-type stars (e.g.McDonald et al. 2019;Kanodia et al. 2021).This appears to be confirmed observationally (Petigura et al. 2022), constituting additional evidence for the importance of mass loss in sculpting the lower part of the desert.With the continued success of TESS and the ability to probe mass-loss rates with metastable helium, we now have the means to unveil the divergent evolutionary pathways of planets on either side of the Neptune desert.

ACKNOWLEDGMENTS
We thank the Palomar Observatory telescope operators, support astronomers, and directorate for their support of this work, especially Tom Barlow, Andy Boden, Carolyn Heffner, Paul Nied, Joel Pearman, Kajse Peffer, and Kevin Rykoski.We also thank Aaron Bello-Arufe, Yayaati Chachan, Philipp Eigmüller, Akash Gupta, Julie Inglis, Shubham Kanodia, Katja Poppenhaeger, Hannah Wakeford, and Nicole Wallack for insightful conversations.We acknowledge the referees for thorough reports that improved the quality of this paper.SV is supported by an NSF Graduate Research Fellowship.HAK acknowledges support from NSF CA-REER grant 1555095.AO gratefully acknowledges support from the Dutch Research Council NWO Veni grant.In Figure 12, we present the posteriors on the limb darkening coefficients for all stars analyzed in this work alongside calculations for the coefficients with ldtk (Husser et al. 2013;Parviainen & Aigrain 2015).

B. DETRENDING VECTOR SELECTION
In Table 6, we summarize the ∆BIC values obtained when comparing models with different combinations of detrending vectors.-26.87 -25.87 -22.53 -21.66 -21.31 -18.08 Note-All ∆BIC values are relative to that for no additional detrending vectors, with negative values indicating favored models.Vector 1 corresponds to the centroid offsets, vector 2 corresponds to the water absorption proxy, and vector 3 corresponds to the airmass.Bolded values indicate the adopted detrending model for a given planet; if there are no bolded values then the model with no additional detrending vectors was selected.
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.

C. ALLAN DEVIATION PLOTS
In Figure 13, we present the Allan deviation plots for all light curves analyzed in this work.2.

D. TESS LIGHT-CURVE FIT
In Figure 14, we present the light-curve fit for the WASP-177b TESS data.We found a radius ratio in the TESS bandpass of 0.187 +0.054  −0.036 .The rest of the fit parameters are given in Table 4 as this was a joint fit with the WIRC data.

E. THE PARKER WIND MODEL WITH TIDAL GRAVITY
We demonstrate that a one-dimensional transonic Parker wind model should achieve the sound speed: within the planetary Roche lobe (R s < R Roche ) when the tidal gravity term is considered (Murray-Clay et al. 2009).Here, M p is the planet mass, G is the gravitational constant, k B is the Boltzmann constant, T 0 is the isothermal outflow temperature, μ is the average molecular weight as defined by Lampón et al. (2020), and m p is the proton mass.The sonic point without considering the stellar gravity is defined as: As in the original Parker wind solution, we rewrite the pressure gradient term with a density gradient using the derivative of the ideal gas law.We then rewrite the density gradient term as a velocity gradient term using the derivative of the continuity equation (a step-by-step walkthrough of this method is provided in Chapter 3 of Lamers & Cassinelli 1999).The result is: The sonic point is the singular point of Equation (E4) where the numerator and denominator of the right-hand side both go to zero.The latter condition gives v = c s , identically to the original Parker wind solution.Requiring the numerator to go to zero gives the cubic: 3GM r 3 a 3 + 2c 2 s r − GM p = 0.
(E5) This is the major difference from the Parker wind considered by Oklopčić & Hirata (2018).Neglecting the tidal gravity term, we obtain the original sonic point solution in Equation ( E2), but the general solution behaves differently: where M 1 is further defined as: and v K = GM /a is the Keplerian velocity.Equation (E6) agrees with Equation (E2) in the limit of large sound speed, but the crucial point here is to consider the case where the sound speed is very small, corresponding to cold outflows.Whereas the original sonic point solution would have these outflows fall outside the Roche radius, in this limit M 1 = M p /2, and substitution into Equation (E6) yields R s = R Roche .A similar result is obtained in Appendix C of Tang & Young (2020).Therefore, when the sonic point in our original models (which neglected the tidal gravity term) began to venture outside the Roche lobe, it meant that our model was not providing an adequate representation of the density and velocity profiles, which are dominated by the stellar gravity.In general, Equation (E7) shows that stellar gravity begins to dominate the outflow behavior when: 32 81 Expanding out the Keplerian velocity term on the lefthand side, this condition scales with c 6 s a 3 /(M M 2 p ); therefore the tidal gravity term is important in planets with cool outflows and planets in close proximity to more massive stars.For our sample, the relevant quantities are approximately v K ∼ 100 km/s, c s ∼ 10 km/s, and M p /M ∼ 3 × 10 −4 , and Equation (E8) is nearly satisfied, indicating that the tidal gravity term in the momentum equation appreciably alters our results.The new momentum equation (Equation E4) is readily integrated from the sonic point (R s , c s ) to an arbitrary point (r, v(r)) to obtain the velocity profile: and using the continuity equation we can obtain the density profile: The updated Parker wind structure has been included into the latest release of the p-winds code (dos Santos et al. 2022).

Figure 1 .
Figure 1.Transiting exoplanets with fractional mass uncertainties of less than 30%, with the dashed lines indicating the Neptune desert boundaries from Mazeh et al. (2016).Planets orbiting stars with 4000 K< T eff < 5400 K are outlined in red, and the seven targets constituting our sample are outlined in blue.The point color denotes the host star's J magnitude, and the point size is proportional to logarithm of the predicted SNR in the WIRC helium bandpass from Section 2.1.

Figure 2 .
Figure 2. Narrowband light curve (top) and residuals (bottom) for WASP-69b.Grey points are detrended data, which are binned to 10 minute cadence in black.The red curve indicates our best-fit model, with 1σ uncertainty indicated by the red shading.The blue curve indicates the nominal light curve model (outside the helium bandpass).

Figure 3 .
Figure 3. Same as Figure 2, but for WASP-52b.Points with circular markers indicate the first night of data collection and points with triangular markers indicate the second night.

Figure 4 .
Figure 4. Same as Figure 2, but for HAT-P-18b.Points with circular markers indicate the first night of data collection and points with triangular markers indicate the second night.
The fit results are shown in Figure5.Our non-detection agrees with the more sensitive observations fromFossati et al. (2022).

Figure 5 .
Figure 5. Same as Figure 2, but for WASP-80b.Points with circular markers indicate the first half of the observations (prior to the instrument reset) and points with triangular markers indicate the second half.

Figure 6 .
Figure 6.Same as Figure 2, but for WASP-177b.The blue shading indicates the non-negligible uncertainty on the baseline transit model for this planet.

Figure 7 .
Figure 7. Same as Figure 2, but for HAT-P-26b.Points with circular markers indicate the second night of data collection and triangular markers indicate the third night.For both the second and third nights, the dashed red lines indicate 1.5× the photon noise.

Figure 8 .
Figure 8. Same as Figure 2, but for NGTS-5b.Points with circular markers indicate the first night of data collection, and points with triangular markers indicate the second night.

Figure 9 .
Figure 9. Mass-loss models for the seven planets in our sample.The shading gives the σ agreement between the observations and a 1D transonic Parker wind model characterized by the mass-loss rate Ṁ and thermosphere temperature T0.Models that violate energy balance (Vissapragada et al. 2022) are rejected (blue shaded regions).

F
Rates for WASP-80b and WASPXUV /ρ p (erg cm g −1 s

Figure 11 .
Figure11.Mp − a plane with marginal stability curves.Planets above the red curve cannot have lost more than 50% of their initial mass M0 to photoevaporation, even assuming an energy-limited outflow for the entire planetary lifetime.The red shaded region corresponds to the uncertainty on our empirical estimate for the mass-loss efficiency ε with the uncertainty in the empirical radius-temperature relation fromOwen & Lai (2018) included as well.The orange curve shows the same for 10% mass loss.Points indicate transiting exoplanets with fractional mass uncertainties of less than 30%; those orbiting stars with 4000 K< T eff < 5400 K are colored blue, and the seven targets constituting our sample are labeled and outlined with red stars.The gray shaded region indicates Mp < 0.2MJ wherein our assumed radii are incorrect.

Figure 13 .
Figure 13.Allan deviation plots for all light curves analyzed in this work.The rms of the binned residuals for each light-curve fit are shown with the black curves, the red noise indicates the expectation from photon noise statistics, and the dashed red line is the red line scaled up to match the first point of the black curve.The photon noise values and scaling factors are given in Table2.

Figure 14 .
Figure 14.Same as Figure 2, but for the WASP-177b TESS data.
derive the sonic point with tidal gravity considerations, we first write the full momentum equation for the wind (e.g.Equation (2) from Murray-Clay et al. 2009):

Table 1 .
Stellar and planetary parameters for sample selected in this work .

Table 2 .
Summary of Palomar/WIRC observations analyzed in this work .

Table 3 .
Priors for the transit light-curve fits.

Table 4 .
(Husser et al. 2013;Parviainen & Aigrain 2015)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.

Table 6 .
Summary of ∆BIC between detrending models in fits with free limb-darkening coefficients .