A Large and Variable Leading Tail of Helium in a Hot Saturn Undergoing Runaway Inflation

Atmospheric escape shapes the fate of exoplanets, with statistical evidence for transformative mass loss imprinted across the mass–radius–insolation distribution. Here, we present transit spectroscopy of the highly irradiated, low-gravity, inflated hot Saturn HAT-P-67 b. The Habitable Zone Planet Finder spectra show a detection of up to 10% absorption depth of the 10833 Å helium triplet. The 13.8 hr of on-sky integration time over 39 nights sample the entire planet orbit, uncovering excess helium absorption preceding the transit by up to 130 planetary radii in a large leading tail. This configuration can be understood as the escaping material overflowing its small Roche lobe and advecting most of the gas into the stellar—and not planetary—rest frame, consistent with the Doppler velocity structure seen in the helium line profiles. The prominent leading tail serves as direct evidence for dayside mass loss with a strong day-/nightside asymmetry. We see some transit-to-transit variability in the line profile, consistent with the interplay of stellar and planetary winds. We employ one-dimensional Parker wind models to estimate the mass-loss rate, finding values on the order of 2 × 1013 g s−1, with large uncertainties owing to the unknown X-ray and ultraviolet (XUV) flux of the F host star. The large mass loss in HAT-P-67 b represents a valuable example of an inflated hot Saturn, a class of planets recently identified to be rare, as their atmospheres are predicted to evaporate quickly. We contrast two physical mechanisms for runaway evaporation: ohmic dissipation and XUV irradiation, slightly favoring the latter.


INTRODUCTION
Atmospheric escape appears to be a normal phase for exoplanets of a certain size and insolation, statistically * NASA Sagan Fellows imprinted in the dearth of 1.5−2.0R ⊕ planets within 100-day periods (Fulton et al. 2017).The leading candidate physical mechanisms for this rapid transition between mini-Neptunes and super-Earths include photoevaporation (Owen & Wu 2013;Lopez & Fortney 2013;Owen & Wu 2017) and core-powered mass loss (Gupta & Schlichting 2019;Berger et al. 2023).Whatever the cause, some large fraction of planets undergo transformative atmospheric escape, and the signal should be widely discernable.Such signals have been searched for, and increasingly found, in many transiting planet systems with at least 28 detections to date (Dos Santos 2022).
Uncertainty in system ages, evaporation timescales, X-ray/UV radiation, and dominating physical mechanisms degrade our ability to foretell if any given planet will exhibit ongoing signatures of atmospheric escape.Episodic stellar wind gusts and other forms of astrophysical variability could also subdue the appearance of atmospheric escape, even where we expect it most.The over 57 published non-detections of atmospheric escape (Dos Santos 2022; Guilluy et al. 2023) must encode these natural whims in a way that we have not yet disentangled.Nevertheless, we can boost our chances of witnessing active and significant atmospheric escape by targeting sources that seem predisposed to loss.These intrinsic or extrinsic factors may include proximity to the host star, low surface gravity, and high energy incident radiation.
Inflated hot Saturns stand out as an especially extreme category that should exhibit mass loss.This category is defined as having masses comparable to Saturn's (M p ∼ 0.3M Jup ), with equilibrium temperatures high enough to expect radius inflation (T eq > 1000 K).Their low gravitational potentials should let go of their atmospheres more readily than their hot Jupiter counterparts.Lower gravity also implies larger atmospheric scale heights, making them easier to detect in transmission spectroscopy.And their large transit depths and short periods should make them readily detectable in transit searches in large numbers, like hot Jupiters.
But inflated hot Saturns are rare (Thorngren & Fortney 2018).The cause for their underabundance remains an open question, with at least two conceivable explanations.Tidal migration mechanisms-either higheccentricity or disk-based-could hypothetically proceed in a mass-dependent manner, efficiently for Jupitermass planets, but inefficiently for the lower-mass Saturns (Thorngren & Fortney 2018;Dawson & Johnson 2018).In this scenario, sub-Saturn mass planets never make it to the close-in orbital separations that would lead to the conditions needed for inflation in the first place.
Alternatively-and most consequentially for atmospheric escape-another explanation may prevail.Inflated hot Saturns may either form in-situ or effectively migrate to close-in orbital separations (Dawson & Johnson 2018), but once they arrive, the intense irradiation overheats the planet.This heating leads to runaway inflation and, ultimately, complete disintegration.In this scenario, the inflationary half-life becomes so short that the probability of observing members in the class decreases sharply with density, causing the apparent lack of inflated sub-Saturns (Thorngren et al. 2023).Batygin et al. (2011) predicted runaway inflation of hot Saturns as a consequence of the Ohmic dissipation mechanism.Here, lightly thermally ionized atmospheric flows induce drag in a planetary magnetic field, weakly coupling the stellar incident energy into the planetary interior.A key prediction of Ohmic dissipation is that the anomalous heating efficiency ϵ(T eq ) should exhibit a peak around T eq ∼ 1500 − 2000 K (Menou 2012;Rogers & Komacek 2014;Ginzburg & Sari 2016).Thorngren & Fortney (2018) favored Ohmic dissipation as the mechanism responsible for inflating hot Jupiters by showing that the observed sample of inflated planets implies an anomalous heating efficiency peak at equilibrium temperatures of ∼ 1500 K. Therefore Ohmic dissipation stands out as a leading physical mechanism for inflation and mass loss.
Recently, Thorngren et al. (2023) showed that hot Saturns can undergo catastrophic erosion by stellar X-ray and Extreme UV (XUV) photoevaporative mass loss.Planets with densities less than ∼0.1 g cm −3 achieve mass loss rates of up to 10 3 M ⊕ /Gyr, setting up a positive feedback loop: the planet increases in radius, overflowing its Roche lobe, fueling greater mass loss, and increasing in radius even further.This vicious cycle systematically depopulates the mass-radius plane with a cliff defined by the ρ p ∼0.1 g cm −3 dividing line.
A runaway inflation scenario predicts an inevitable and profound mass loss rate for inflated hot Saturns.As the planet's atmosphere overflows its Roche lobe at an ever-increasing pace, instantaneous mass loss rates may exceed Ṁ > 10 13 g/s, over 10× larger than those previously seen in atmospheric escape measurements to date (Dos Santos 2022).Inflated hot Saturns, therefore, make excellent targets for direct measurement of atmospheric escape and for testing the underlying mechanisms of mass loss.
Large mass loss alone does not guarantee detectability.The ability to detect even immense mass loss hinges on its observability in spectral tracers.Lyα, He I 10833 Å, and Hα have emerged as the most amenable to detection (Seager & Sasselov 2000;Vidal-Madjar et al. 2003;Jensen et al. 2012;Yan & Henning 2018;Oklopčić & Hirata 2018;Spake et al. 2018;Dos Santos 2022;Owen et al. 2023), but each of these has its own observational limitations.Lyα can suffer from interstellar medium H I absorption censoring its low-velocity line core, for example.Here we focus on He I 10833 Å, which offers some advantages.In particular, metastable Helium's ability to be observed at high spectral resolution from the ground has been especially valuable for evincing velocity substructure of the escaping gas motion relative to the exoplanet restframe (Alonso-Floriano et al. 2019;Ninan et al. 2020), and convincingly associating the signal to an exoplanetary origin as opposed to stellar contamination (Cauley et al. 2018).He I has resulted in at least 14 systems with detections (Dos Santos 2022)1 .The sample of detections includes both hot Jupiters and lower mass planets but lacks inflated hot Saturns owing to their intrinsic rarity below the 0.1 g cm −3 threshold.An observational picture of mass loss in inflated hot Saturns appears to be lacking for this reason.
The Habitable Zone Planet Finder (HPF) Helium Exospheres program has been conducting a survey of exoplanets to search for atmosphere loss via the He I 10833 Å metastable triplet.The survey's multiple-visit sampling strategy has enabled a search for atmospheric loss at large out-of-transit separations from the planet, recently revealing giant tidal tails of Helium escaping the hot Jupiter HAT-P-32 b (Zhang et al. 2023).Here we present a multi-year observational campaign searching for atmospheric escape in HAT-P-67 b, a transiting inflated hot Saturn orbiting an F5 subgiant at an orbital separation of 0.06 AU and a 4.8-day orbital period (Zhou et al. 2017).Its strong insolation (T eq ∼ 2000 K) combined with HAT-P-67 b's extremely low surface gravity (log g p < 2.3 dex) makes it an exceptional candidate for strong atmospheric mass loss via Roche lobe overflow.Its ∼0.05 g cm −3 density places it below the 0.1 g cm −3 threshold predicted to exhibit runaway inflation (Thorngren et al. 2023).Figures 1 and 2 show how much of an outlier HAT-P-67 b is: large, low mass, and heavily irradiated.
The evolutionary state of HAT-P-67 b offers even more intriguing possibilities.The planet may have undergone re-inflation (Saunders et al. 2022;Grunblatt et al. 2022Grunblatt et al. , 2023)), as the host star evolved through the subgiant phase-with the planet's insolation increasing with the star's rapid luminosity jump in this part of the shortlived HR diagram.The tidal gravity of the nearby massive (M ⋆ ∼ 1.6 M ⊙ ) host star may amplify mass loss rates even further (Erkaev et al. 2007;Thorngren et al. 2023).Its status as a rare inflated hot Saturn makes HAT-P-67 b a promising laboratory, uniquely suited for testing the runaway inflationary predictions of Ohmic dissipation and XUV photoevaporation.
We assemble both archival, previously published, and new observations of the HAT-P-67 system, chronicled in Section 2. In Section 3, we refine the stellar and planet properties based on those observations, including dis-tance, radius, and orbit re-analyses.The Helium excess detection is presented in Section 3.5, with an analysis of the signal's trend with orbital phase.We use this orbital structure to reconstruct the geometry and mass loss of the escaping Helium exosphere with 1D Parker wind models (Section 4).We assess the physical mechanisms (Section 5) giving rise to the escaping atmosphere, weighing the distinctive predictions of XUV irradiation and Ohmic dissipation.Finally, in Section 6, we question the assumptions in our approach, discuss the overall congruence of predictions and observations, and highlight some implications for future exosphere studies.

Habitable Zone Planet Finder (HPF)
The Habitable Zone Planet Finder Spectrograph (HPF; Mahadevan et al. 2012Mahadevan et al. , 2014;;Metcalf et al. 2019) on the queue-scheduled 10-meter Hobby-Eberly Telescope (HET; Ramsey et al. 1998) operates in the near-IR from 8100 − 12800 Å spanning the z, Y, and J bands at spectral resolving power R = 55, 000.The HET fixedelevation design (Shetrone et al. 2007) limits the observability of HAT-P-67 to less than 1 hour "tracks" for a fixed range of hour angles before (east track) and after (west track) the star transits the meridian.Whereas conventional steerable telescopes could conduct continuous point-and-stare observations of HAT-P-67 for hours, HET cannot.In practice, this limitation means that in-transit and out-of-transit observational phases were rarely possible on the same night.Instead, we organized the observations into four campaigns to coincide with HAT-P-67 b transits on UT dates 2020 April 28, 2020;May 22, 2020;June 15, 2020;and 2022 April 29.These campaigns have out-of-transit observations at least one night before and one night after, and often two nights on either side of the transit.Two more transit snapshots were obtained in 2022 June-July without the accompanying visits immediately before and after.The in-transit campaigns had up to 14 exposures per HET track, with integration times between 5-8.5 minutes.We also obtained random-in-phase lower-priority "P1-P4" reconnaissance observations (Shetrone et al. 2007).These out-of-transit snapshot observations typically received 4 or fewer individual exposures.We observed HAT-P-67 with HPF for a total of 41 visits on 39 unique nights, with two of those nights observing both the east and west HET tracks.The total on-source integration time exceeds 13.8 hours, with N = 152 individual spectra possessing a typical signal-to-noise ratio of 80 per pixel.Table 1 summarizes the log of observations.
The observation schedule was strategically restricted to the spring season when the Barycentric Earth Radial Velocity (BERV) would Doppler shift telluric absorption lines sufficiently far away from the core of the Helium 10833 Å feature (Spake et al. 2022).The small BERV still means the redward Helium line wing has significant telluric contamination between 10834−10836 Å, but the core and blue line wings appear relatively pristine.Sky emission lines land in this region but are more easily mitigated by HPF's simultaneous sky reference fiber.
Only 19 out of the 39 nights possess A0V telluric calibration standard stars.These standard stars were used to spot-check our telluric masking.Figure 3 indicates the epochs of select HPF visits overlaid on the TESS time-series lightcurve.

TESS Light Curves
HAT-P-67 was observed with the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) in Sectors 24, 26, 51, 52, 53 with 2-minute cadence, and in Sector 25 with 30-minute (FFI) cadence.The Sector 25 FFI data fell just barely off the TESS detectors, in collateral pixels where no starlight ever hit the detector.
We also assembled a comparison sample of about 1000 lightcurves to interpret the prevalence of lightcurve modulation and stellar activity among broadly F subgiant-like stars.We selected sources based on Gaia DR3 T eff estimates and similar log g, and availability of at least one sector of TESS 2-minute cadence data.We visually spot-checked these lightcurves to understand artifacts and windowing effects.

Gaia DR3
The stellar system consists of a binary with an M dwarf companion HAT-P-67B (Gaia DR3 1358614983131339904 ) separated on-sky by 9. ′′ 0 (Mugrauer 2019), well-separated from the planet-host star and not a source of contamination for the HPF observations.HAT-P-67A (Gaia DR3 1358614983131339392 ) has a parallax of 2.69 ± 0.01 mas in Gaia DR3, placing it at about 372±1.4pc.Minor corrections to the parallax uncertainty (El-Badry et al. 2021) and bias (Lindegren et al. 2021) appear negligible for the G = 9.98 mag source.The DR3 parallax places the system about 8.7% farther than previously estimated by Zhou et al. (2017), which adopted a Gaia DR1-informed parallax of 2.92 ± 0.23 mas, including a systematic −0.325 mas bias term (Stassun & Torres 2016).The wide companion HAT-P-67B has nearly identical parallax (2.58±0.05mas) and proper motions, confirming its interpretation as co-moving, with a projected separation of about 3400 AU.The IAU naming convention would demand HAT-P-67Ab to refer to the planet, which orbits the primary.Hereafter, we simply drop the A designation for notational simplicity since the wide companion will not factor into our analysis.

ASAS-SN, DASCH, and ZTF
We retrieved ground-based photometry with the All-Sky Automated Survey for Supernovae (ASAS-SN) using the Sky-Portal (Shappee et al. 2014;Kochanek et al. 2017).The precision of ASAS-SN was too low to perceive stellar variability, so we can place a relatively uninformative limit of < 5% stellar variability on years-long timescales.
Similarly, HAT-P-67 appears in the Harvard Plate Archive, with 5809 measurements digitized through the Digital Access to a Sky Century at Harvard (DASCH) program, spanning over 120 years of coarse photometric monitoring.In principle, these datasets could inform long-term variability trends such as stellar cycles.In practice, the 0.15 magnitude jitter appears too coarse to perceive any genuine astrophysical variability, with no conspicuous trend seen.We can therefore place a relatively mild constraint that the star appears stable at the ∼ 30% level over periods of tens to hundreds of years.
HAT-P-67A was saturated in ZTF (Bellm et al. 2019), but the M-dwarf companion HAT-P-67B had up to thousands of visits across several years.The data quality appeared too poor to perceive any genuine astrophysical variability, with the indication of some lunar background signals in the periodogram.Zhou et al. (2017) previously derived stellar radius estimates of 2.1-2.7 R ⊙ through Spectral Energy Distribution (SED) and isochrone fitting as part of a joint orbit fit.We systematically increase those stellar radius estimates by 8.7% to match the greater Gaia DR3 distance ( §2.3).For a fixed R p /R ⋆ from the measured transit depth, the larger R ⋆ implies a proportionally larger planet radius.This update systematically decreases the estimate for the already-low density of HAT-P-67 by 28%, to a mere <0.035 g cm −3 , albeit with significant uncertainties from the weak mass constraint.The luminosity increases by 18%, to about 10.3 L ⊙ .

TESS Light Curve
Previously, HAT-P-67 b transits had only been detected with HATNet (Bakos et al. 2004) and followed up with KeplerCam on the FLWO 1.2 m telescope (Zhou et al. 2017).These ground-based photometers were not intended to measure weak, long-term stellar variability signals.We, therefore, examined the TESS lightcurves for out-of-transit photometric variability.The revised precision and continuous coverage of TESS can also refine the orbital solution.

Revised exoplanet orbital parameters
We assembled a composite lightcurve by stitching TESS Sectors 24, 26, 51, 52, 53, which were reduced with the default SPOC pipeline (Caldwell et al. 2020), and lightly post-processed with lightkurve (Barentsen et al. 2019).These TESS data exhibit an RMS scatter of better than 1 part-per-thousand (ppt) at native 2minute sampling.We fit a Keplerian orbit model to the TESS lightcurve using the exoplanet framework (Foreman-Mackey et al. 2021).We obtained revised orbital properties shown in Table 2.These properties are broadly consistent with the previously reported values from Zhou et al. (2017) and the updated ephemeris of Ivshina & Winn (2022).Figure 3 shows an overview of all the TESS Sectors with a Gaussian Process trendline in green highlighting the out-of-transit modulations, and the exoplanet transit model in orange.Figure 4 shows the best-fit orbit overlaid in purple on the detrended TESS lightcurve, which is binned in the orange dots.Table 2 lists only one of the previous or-bit determinations-also assuming a circular orbitcompared with the orbital solution reported here.

Stellar rotation rate from periodogram analysis
The TESS Sector 26 lightcurve exhibits a weak ∼3 ppt peak-to-valley out-of-transit modulation.TESS Sectors 51-53 do not show as conspicuous a modulation signal but still show some ostensibly stellar variability.The ebb and flow of the modulation can be seen in the minimally processed TESS lightcurve in Figure 3.We fit the TESS lightcurve modulation with a quasiperiodic Gaussian Process (GP) model using celerite (Foreman-Mackey et al. 2017;Foreman-Mackey 2018).The GP model fit to the entire composite lightcurve yields a 4.7day period; when fit to individual sectors alone, the periods hover around 5.9 days.The signal is weak enough, especially in Sectors 51-53, that TESS instrumental systematics may contribute some of the out-of-transit modulation that we see in Figure 3.
We can independently constrain the expected stellar rotation period based on the stellar radius estimate, v sin i ⋆ , and the observation of the planet's orbital inclination i p ∼ 90 • from orbit fitting and Doppler tomography (Zhou et al. 2017).We assume spin-orbit alignment, i ⋆ ∼ i p .We adopt a high limit of v sin i = 35.8± 1.1 km/s and a lower v sin i = 30.9± 2 km/s value if macroturbulence is accounted for.We obtain a range of 3.2 < P rot < 4.8 days.This range is typical for F stars that have not yet evolved too far into the subgiant branch (Avallone et al. 2022).The geometrical constraint comports with the 4.7-day GP-based modulation period derived from the stitched TESS lightcurve but is lower than the 5.9-day per-Sector fits, slightly preferring the lower 4.7-day modulation as the stellar rotation period.The 5.9-day value would require a stellar radius over 3.2R ⊙ , which appears implausible.

Revised system evolutionary state
The slightly increased distance and radius imply an 18% more luminous host star than previously estimated.Figure 5 shows evolutionary tracks using the solar-metallicity MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015) ) Isochrones and Stellar Tracks (MIST; Dotter 2016; Choi et al. 2016) for HAT-P-67, which has [Fe/H]= −0.08 ± 0.05 (Zhou et al. 2017).The tracks show that a 1.8 Gyr, 1.64 M ⊙ HAT-P-67 could be at the end of its main sequence lifetime (teal track), yielding a gradual rise in incident radiation over the lifetime of HAT-P-67 b.Alternatively, a 1.54 M ⊙ evolutionary track would send HAT-P-67 along the subgiant branch at an age over 2.2 Gyr, with a rapid 50% increase in luminosity, potentially leading to "re-inflation" of the planet (Thorngren  Either evolutionary track appears consistent with the observed SED, stellar surface gravity, and available constraints from spectral fitting.We, therefore, adopt a MIST-based 2.0±0.2Gyr age for the system, in-between the previous Geneva and Dartmouthbased scenarios (Zhou et al. 2017).Additional metallicity and rotation effects could slightly alter the evolutionary tracks and therefore increase uncertainty in inferred ages and masses.We conducted orbit fitting via precision radial velocity (PRV) measurements following the procedures described in Tran et al. (2021) and updated for joint lightcurve and RV fits (Tran et al. 2022), but without a joint Gaussian Process modeling procedure (Tran et al. 2023).Here we included the new TESS lightcurves and the existing Keck HIRES points reported in Table 3 of Zhou et al. (2017).The HPF points exhibited an RV jitter of ∼ 164 m s −1 , a few times larger than the typical Keck HIRES measurements, due to the lower information content in the near-IR than the visible for this relatively rapidly rotating F star.So even though the HPF points were more numerous, their marginal value for RV orbit determination was subdued.The period and T 0 were fixed from the TESS transits, and we assume a circular orbit with K > 10 m s −1 .The semiamplitude prior acts as a mass constraint, excluding planets with masses so low that the observed radius would exceed the Hill radius (Roche lobe overflow).Figure 6 shows the fitted RV semiamplitude K = 33 +21 −15 m s −1 , a relatively weak constraint but consistent with broadly Saturn-mass planets.Table 2 lists the revised properties, with minor differences from the existing values.The semiamplitude prior accounts for the apparent improvement in the mass constraint.Hypothetically, a truly Roche lobe overflow planet could be consistent with the available data, so our mass constraint could be considered an upper limit.

HPF analysis II: He I 10833 Å
The preparation for analyzing the Helium feature followed a similar but slightly different procedure than that used for the RV.Namely, we used Goldilocks2 for 2D échellogram reduction.This tool outputs 1D extracted spectra for the target and two reference fibers: blank sky and a laser frequency comb (LFC, Metcalf et al. 2019).The observations were acquired with the LFC turned off, so this unilluminated spectrum was discarded.
The sky fiber and target fiber have slightly different throughputs and illumination properties, with the target fiber receiving 93% of the flux of the sky reference fiber on average.This ratio depends on wavelength and season at the few percent level.We quantified the wavelength dependence by acquiring calibration observations of blank sky in both the target and sky fiber.We applied this wavelength-based scale factor to each target spectrum's associated reference sky fiber to achieve skyline subtraction residuals typically less than the photon noise (Gully-Santiago et al. 2022).A few lines exhibit residuals that may arise from genuine differences in the local atmospheric conditions between the target and sky fiber.
The HET's fixed-altitude design means that the airmass remains relatively constant across all pointings, lowering the telluric line variability compared to fully steerable telescopes, which may sample a wider range of airmasses.Variability in atmospheric conditions still makes it difficult to mitigate telluric lines to within the photon noise limit, so we masked spectral regions predicted to have significant telluric lines with a template generated by telfit (Gullikson et al. 2014).We then shifted the spectral coordinates to their common barycentric-corrected reference frame (Wright & Eastman 2014), as implemented in astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)).The spectral continuum was flattened and normalized to two pre-selected continuum indices highlighted as vertical blue bands in Figure 7 shows a zoom-in on the He 10833 Å region of interest, with all 152 individual exposures overlaid.We see variability in the Helium line of up to 10%, much greater than the < 1% pixel-to-pixel variation.The feature width spans about 3 Å.
Figure 8 shows the spectra for four campaigns, with before and after spectra showing conspicuous excess absorption during transit and 1 day before transit but negligible excess absorption after transit.The individual transits show significant morphological variation, with substructure in the bulk line-of-sight velocity distribution.
The 13.8 hours of exposure, combined with HAT-P-67 b's short (P ∼ 4.8 day) period, means that a large fraction of the orbit has been collected, with some phases (especially near transit) heavily sampled, and some other out-of-transit phases sampled more sparsely.The largest gap spans merely 0.18 in phase, as seen in the Equivalent Width time series in Figure 9.For visualization purposes, we constructed a 2-D intensity phase scan, binning in phase and wavelength, and spanning the entire orbit of HAT-P-67 b. Figure 10 overlays the planet's approximately ±150 km/s orbital Doppler velocity, with the stellar rest frame velocity demarcated by the vertical line.The vanishing < 36 m/s reflex motion of the star is imperceptible at this scale.The horizontal dashed lines indicate the moments of transit ingress (-0.03), mid-center (0.0), and egress (+0.03), demarcated as TRANSIT in the figure.
We define 4 additional distinct groups of phases with adjectival qualifiers in anticipation of the need to discuss bulk trends: PRE, POST, CONTROL, and BASELINE.The BASELINE phases define the out-of-transit baseline.The "control" spectra designate phases just before the baseline but are not used to define the baseline, to inspect minor correlation structure in the spectra without dividing it out.Table 1 lists the adjectival qualifiers for each HPF visit.
The in-transit phases exhibit over 10% excess absorption peaks.A large absorption signal can be seen preceding transit ingress, with some significant absorption evident before −0.2 in phase.The egress drop-off is extremely sharp, with almost no excess absorption directly after transit, as seen in the Equivalent Width (EW) time series (Figure 9).We construct a fractional residual absorption spectrum by subtracting off and re-normalizing to the non-varying baseline.We define the non-varying baseline spectrum as the average over phases 0.4 − 0.5, which exhibited stable spectra with the least absorption.Figure 10 shows that significant absorption can be seen as a few percent at −0.37 in phase, rising to 10% just before and during transit.The abrupt dropoff in helium excess at planet egress is conspicuous.
The underlying structure of the metastable He I triplet consists of three quantum components, two of which (J = 1 and 2) are typically blended owing to Doppler broadening from the finite temperature of the gas T 0 .HAT-P-67 b exhibits blending of all three components (J = 0, 1, and 2) into a single broad Gaussian-like feature seen in Figure 7.The mere observation of this broadening implies that HAT-P-67 b probes a larger velocity dispersion than typical measurements and may be associated with either a higher gas temperature, complex planetary wind flows, or some mix of both.Whatever the cause, we treat the feature as a single Gaussian line for the purpose of estimating bulk characteristics of the Helium excess absorption feature.We restricted the fitting to the 10828−10838 Å region of the residual spectra.The model constructed in this way has a total of 4 parameters: amplitude A, location λ c , Gaussian width σ, and constant offset.We repeated a similar process for the nearby Si line at 10830 Å as a control sample.
Excess is confidently detected from phases −0.3 to +0.03.The full-width at half maximum (FWHM) of the feature is about 1.8−2.6Å, equivalent to a line-of-sight velocity distribution in the range of 50-75 km/s.Individual fits to the He feature show a typical line center uncertainty ±0.1 Å or better.The out-of-transit phases exhibit a line center position of 10833.4Å, consistent with the stellar He I rest frame zero velocity.The excess absorption resides slightly blueward of the stellar zero restframe velocity reference.As mentioned previously, the telluric masking near 10836 Å censors our view of the existence (or not) of any redshifted lobe in this region.The excess detections exhibit a bulk blueshift of up to 30 km/s.
The excess absorption spectrum exhibits an increasing line-of-sight velocity distribution as the planet transits, going from 60 km/s at ingress, up to 100 km/s at egress.
Finally, we detect a much weaker but still significant post-egress blue-shifted absorption.This lobe appears to decrease linearly in wavelength, corresponding to a characteristic blueshift increasing from near zero at transit midpoint to 80 km/s at 12 hours after transit.It exhibits a slightly lower line-of-sight velocity distribution of 50 km/s, about half of the peak during transit.
We examined other spectral lines using the same methodology as above, finding no detectable variability in the Ca II infrared triplet, Hydrogen lines, or other neutral metal lines.We show some of the diagnostic plots for these non-detections in Appendix Section A.

Keck HIRES
We retrieved 22 epochs of archival Keck HIRES spectra via the Keck Observatory Archive (KOA).Of these, 19 were acquired through the iodine cell as reported by Zhou et al. (2017), providing RV orbit constraints.The 14 ′′ slit height of the C2 decker appears to cause order overlap in Ca II H and K lines, which contaminated the spectral extraction for 14 of the 22 spectra.The 3. ′′ 5 B2 decker allowed the faithful extraction of this region for 6 spectra, which exhibit no perceptible Ca II H and K variability.The Hα region exhibited no perceptible variability, but the timing of these spectra coincidentally missed the orbital phases (−0.3 to +0.03) where we see the greatest He I 10833 Å absorption excess, leaving open the possibility that detectable Hα absorption could be present at more favorable orbital phases.New spectra or a more careful extraction of the archival C2 decker spectra would be needed.

Velocity substructure
Figure 11 shows the centroid positions of the Helium excess feature.The bulk velocity drifts of the planet's outflow can be seen with a few conspicuous trends.The pre-transit phases (−0.3 to +0.03) show a slight tendency towards blueshift relative to stellar restframe, with an accelerating blueshift from phases −0.3 to −0.15, albeit with some visit-to-visit scatter.The post-transit observations (+0.03 to +0.15) exhibit weak-albeit-significant equivalent widths, as denoted by their smaller marker size.The centroids dramatically accelerate to larger bulk blueshifts, as low as 10831 Å.
At the moment of transit midcenter, the centroids exhibit a redshifted absorption relative to the planet restframe.Measurements from April 2022 and May 2020 transits partially overlap in phase coverage while yielding slightly different line centroid locations.Both time series sequences appear consistent with near-zero bulk restframe velocity but with significant variability in the line profile substructure.The May 2020 sequence shows a slight trend towards alignment with the planet restframe.The April 2020 and June 2020 partial transits sampled phases close to planet egress, with a clear blueshift relative to both the star and planet.
Overall this pattern of velocity centroids appears consistent with the majority of the gas rapidly obeying the stellar restframe but with a large spread in the line-ofsight velocity distribution, as we would expect from a range of launch angles.We revisit the interpretations and caveats for this substructure in the next sections.

Signal Inconsistent with Stellar Activity Interpretation
The planet's orbital period is close to the stellar rotational period, so the prospect of stellar activity contamination arises.Overall, the HPF spectra disfavor a stellar activity interpretation for a few reasons.First, a hypothetical stellar origin of Helium 10833 Å variability should be accompanied by variability in other tracers.Instead, the Ca II Infrared Triplet (IRT) shows stable line profiles in our HPF spectra.The Appendix Section A discusses these and other non-detections.Second, the velocity substructure and phasing appears inconsistent with stellar activity: the gradual rise and then abrupt truncation of the Helium excess at the moment of planet egress (Figure 9) appears inconsistent with a more smoothly varying stellar variability.Finally, the weak post-egress excess appears blueshifted to wavelengths as short as 10830 Å, which would fall entirely outside of the original Helium spectral line's rotationally broadened profile, ruling out a heritage from the star's surface.
The recent discovery of giant tidal tails in HAT-P-32 b (Zhang et al. 2023) provides an additional anchor.HAT-P-32 b is in many ways the most analogous known system to HAT-P-67 b, having comparably low mass and a comparably large radius for similar insolation around an F star.In HAT-P-32 b the planetary interpretation is unambiguous owing to the different stellar rotation and orbital periods.In Section 5.4 we show that radius inflation mechanisms expect similar and large mass loss rates for both of these systems.

Size and geometry of outflow
Helium absorption occurs predominantly before the planet is in-transit, indicating a leading tail escaping the planet.The leading tail is detected 0.1 AU away from the planet, or 130 planetary radii, well outside of the Roche lobe at < 2.6 planetary radii.We are detecting Roche lobe overflow, escaping material that is entirely unbound from the planet.The bulk of the detectable unbound gas may be understood as following its own Keplerian trajectory at a slightly shorter period orbit, pulling away from the planet parallel to the orbital path.Such a scenario may arise from preferential emission on the planet dayside.
The large size of the gas stream overfills the transit chord, which means both the advancing blue side and receding red side of the rotating star are nearly continuously obscured in this spin-orbit-aligned system (Zhou et al. 2017).If the planet transits near the stellar equator (Zhou et al. 2017), the transit chord covers about 10% of the projected stellar surface area.Hypothetically a nearly 100% opaque optically thick stream of gas covering only the transit chord could reproduce the observed 10% absorption excess.Alternatively-and more likely-an optically thin gas stream has some additional vertical extent perpendicular to the plane of orbit, which would act to overfill the entire projected stellar disk.Assuming the gas has a spatially uniform optical depth of about 0.1 would reproduce the 10% absorption excess.
The leading tail geometry seen in HAT-P-67 b contrasts with the morphology of the HAT-P-32 b system, which exhibits giant symmetric tidal tails with comparable absorption depth preceding and trailing the planet (Zhang et al. 2023).These systems represent two of the most extended features associated with exoplanets, and we compare and contrast them in Section 6.

1D Parker wind models and mass loss rates
We explore one-dimensional (1D) models, which offer ease of interpretation and rapid calculation.We employ the open-source p-winds code (Dos Santos et al. 2022), which implements a transonic 1D Parker wind model with radiative transfer of the Helium 10833 Å triplet (Oklopčić & Hirata 2018;Lampón et al. 2020).
The predicted Helium ionization depends sensitively on the XUV flux (Oklopčić 2019).Ideally, we would have a measurement of the XUV spectrum of HAT-P-67, but the limited facilities and distance of the source preclude these challenging observations.Instead, we constructed panchromatic X-ray to visible SEDs by scaling and stitching together synthetic spectra to provide a range of high and low L X /L bol .We chose a 6400 K PHOENIX (Husser et al. 2013) solar metallicity photospheric spectral template with log g = 3.5-the closest PHOENIX grid point to published and updated values available-scaled to the solid angle seen by HAT-P-67 b.The F7IV-V τ Boo serves as the closest analog with available synthetic X-ray coronal spectra (Sanz-Forcada et al. 2011).But HAT-P-67 resides closer to the Kraft break than does τ Boo, and the move to hotter F stars may be associated with weakening of the stellar wind and other atmospheric changes leading to the prospect of much lower XUV luminosity for HAT-P-67 (Avallone et al. 2022).So there remains significant uncertainty about the applicability of the τ Boo XUV spectrum to HAT-P-67.We scaled the synthetic SED of τ Boo retrieved from X-exoplanets (Sanz-Forcada et al. 2011) so that the integral of ground-state ionizing photons (λ < 504 Å) exhibited L X /L bol ∈ [10 −6 , 10 −4 ]. Figure 12 shows two conceivable SEDs with these high and low radiation hardnesses.
The SED shows that the XUV corona model does not extend redward to 2600 Å photons capable of ionizing out of the 2 3 S Helium metastable state.This apparent deficit should be negligible since the F spectral type obtains most of its NUV photons (λ ∼ 2600 Å) from the Wein side of the photospheric spectrum, so the hardness of the radiation stems mostly from the assumptions of the corona flux level.The wavelength-dependent cross section for absorption (not shown) is largest just blueward of the ionization thresholds, so the differences in the spectrum between 504 and 1000 Å have relatively little impact on the overall ionization out of the metastable state.
Equipped with these two SEDs of differing radiation hardness, we explored different mass loss rates and exosphere temperatures.Figure 13 shows one example model spectrum with Ṁ = 2 × 10 13 g/s, and T 0 = 14 000 K, with an XUV spectrum possessing L x /L bol = 10 −5 .This mass loss rate would imply a characteristic lifetime of <1000 Myr.The model spectrum exhibits approximately the same equivalent width as our HPF observations, making it a hypothetical scenario among a family of partially degenerate solutions.
At least a few shortcomings limit the applicability of this 1D Parker wind model.First, the model-dependent line profile (i.e.width and depth) is degenerate with XUV flux, Ṁ , and T 0 , and so a range of these parameters can be fine-tuned to obtain a large range of mass loss rates consistent with the data and our limited understanding of the XUV flux.These known degeneracies have been pointed out previously (Vissapragada et al. 2022;Oklopčić 2019), but the problem for HAT-P-67 b appears somewhat more acute since XUV data are scarce for F-type stars near the Kraft break.Second, the 1D model breaks down when attempting to explain the inherently 3D leading tail geometry.

Direct Evidence for Preferential Dayside Mass Loss
Several physical phenomena could conceivably control the geometry and extent of the escaping material.The stellar potential controls the overall geometry through tides and the Coriolis force, resulting in lobe morphologies that lead and trail the planet (McCann et al. 2019).Here we explore the three predominant dynamical effects: orbital shear, stellar wind confinement, and day-/night-side mass loss asymmetries.
Figure 14 shows an illustration of Keplerian shear adapted to the system properties of HAT-P-67 b.In this shear-dominated scenario, the planetary wind launches primarily from the dayside, with relatively little or no wind launched from the nightside.The planet wind ini- tially launches radially outward from the exobase, with the strongest wind located near the sub-stellar point, the line connecting the planet to the star along the vertical axis in the figure.Inefficient grazing incidence heating near the terminator subdues the mass loss in the x−direction, meaning that this wind exhibits not a hemispherical shape, but more concentration along the star-planet line.The gas increasingly experiences the star's Keplerian potential past the Roche lobe, accelerating in the direction of orbital motion, +x.The accelerating column eventually overtakes the planet completely, with the prospect of extending to hundreds of planetary radii.The mass loss has to be large enough to shield the material to make it observable in the metastable He I 10833 Å triplet.
The observation of such a prominent leading tail leads us to the inescapable conclusion that HAT-P-67 b is predominantly losing mass on the planet's dayside.An isotropic mass loss would manifest a comparably conspicuous trailing tail, which we do not observe.
Importantly, orbital shear, stellar wind confinement, preferential dayside mass loss, and radiative transfer must all conspire to create the conditions of high enough column density to populate the He I metastable triplet to detectable levels while imbuing the velocity substructure that we see.The detailed 3D modeling of the inter- play of these effects is beyond the scope of the current work, but is feasible with adaptations to existing 3D simulations (MacLeod & Oklopčić 2022).
Here we interpret the bulk velocity substructure in Figure 11 under these conceptual gas dynamics mechanisms.We attribute the mild pre-transit arc to stellar wind acceleration.The stellar wind is initially too weak to plow the dense escaping gas until a separation of 60-130 planetary radii (+0.15 orbital phase), the inflection point in an enormous bow shock.Gas at these separations has had enough time to diffuse, both out of the orbital plane and to slightly larger radii shells, in turn lowering the column density along the line-of-sight.The lower column density provides both fewer Helium atoms to participate in absorption, and less overall NUV shielding, leading to greater fractional ionization and further depopulation of the He I metastable state.
The weak post-egress tail can be understood as follows.A weak nightside mass loss means that there is both less overall column density and less NUV photoionization shielding, subduing the overall He I metastable state's signal strength.The blueshift arises from the Keplerian shear that lags the planet, and the ever-blueshifting centroid represents the stellar wind's greater ability to carry away the lower total inertia of less nightside material.
In order to examine the in-transit velocity structure, we have to consider the interplay of two subtle geometrical effects.First, the planet's Roche lobe is small enough (2.7 R p ) that the majority of the projected stellar disk should be filled with escaping material unbound from the planet, even at the time of mid-transit.In other words, only a small fraction of the He I excess signal would be expected to trace the planet's motion, as opposed to the hypothetical in-transit signal for WASP-107 b that definitively traces out the planet's orbital path (MacLeod & Oklopčić 2022).
Second, the star's non-negligible rotational velocity sets up a configuration analogous to Doppler Tomography in the Rossiter-McLaughlin (RM) effect, but distinct in a subtle way.Whereas traditional Doppler Tomography treats the scanning reticle as an opaque planetary disk, here the reticle resembles a transmissive filter with a wavelength center and width that varies in space and time.The interplay of spatial and spectral illumination and absorption can yield minor second-order effects, and overall we anticipate those effects to be secondary compared to the mere existence of absorbing material spanning the entire stellar disk.

Variable planetary wind
Figure 8 and 11 show differences in the Helium line profiles over months-long and years-long timescales.The line profiles are qualitatively similar, but show slightly different tendencies towards redshifting.For example, the centroid of the Helium line during the April 2022 transit appears slightly blueshifted relative to the stellar rest frame.In comparison, the May 2020 transit probes the same planetary phases, yet resides slightly redshifted at transit midcenter.We interpret this line profile variability as indicative of genuine planetary wind variability, which can arise from the interplay of stellar and planetary winds (Murray-Clay et al. 2009) and weatherlike feedbacks in the planetary upper atmosphere.

XUV Irradiation-driven mass loss
XUV photoevaporation heats the upper layers in the atmosphere, driving a hydrodynamic wind (Murray-Clay et al. 2009).Photoevaporation stands out as of-fering a natural cause of the observed leading tail attributable to dayside/nightside differences in mass loss: the greatest supersonic motions arise from the sub-solar point on the planetary dayside.The effect may be boosted if dayside/nightside energy transport proceeds inefficiently (Murray-Clay et al. 2009).Such a scenario is illustrated in Figure 14.
The combined effects of photoevaporation, anomalous heating, and tidal gravity are predicted to have especially drastic outcomes for inflated hot Saturns (Thorngren et al. 2023) such as HAT-P-67 b, where the low densities lead to large mass loss rates: where K t is the tidal gravity term, ρ XUV is the planet density determined using the radius at which XUV radiation is deposited, and η ∼ 0.4 comes from Caldiroli et al. (2022).Assuming log L XUV /L bol = −4.2,we would expect HAT-P-67 b to exhibit Ṁ ∼ 10 2 M ⊕ /Gyr (2 × 10 13 g/s).At a current M p ∼ 95 M ⊕ and an accelerating mass loss rate, its lifetime would measure in the <100 Myr range.The Thorngren et al. ( 2023) simulations focused on 0.75-1.25 M ⊙ host stars; the higher mass ∼1.6 M ⊙ HAT-P-67 would deliver greater tidal gravity for a given semi-major axis, and possibly lower log L XUV /L bol compared to more active G stars.Nevertheless, we can project the time series trends in their Figures 3 and 4 to recreate a qualitative history and fate of HAT-P-67 b under the assumptions of XUV photoevaporation.Broadly, the planet would have started with an initial mass up to 50% greater than at present, with R ∼ 1.4 R Jup , for an initial density of ∼0.2 g/cm −3 .It would lose mass at a rate of tens of Earth masses per Gyr for its 1 Gyr lifetime, expanding modestly until it reaches the critical ∼0.1 g/cm −3 threshold, at which point the mass loss rate increases to its current value.It will last only tens of Myr in its current state before losing almost all of its envelope and settling as a 5 − 15 M ⊕ core with a final radius of 0.2−0.3R Jup .

Ohmic dissipation-driven mass loss
The Ohmic dissipation mechanism made two key observable predictions.Anomalous heating efficiency, ϵ, should initially increase as a function of equilibrium temperature, then degrade as irradiation exceeds equilibrium temperatures of about 2000 K. Second, hot Saturns (≲ 0.5M Jup ) should undergo runaway evaporation, whereas hot Jupiters (≳ 1M Jup ) should reach stable equilibrium-albeit inflated-radii after Gyr timescales.Both of these outcomes predated data that could validate them, and these outcomes differ from the behavior of other heating mechanisms, such as photoevaporation or tides.
Ohmic dissipation may therefore be responsible for both dramatic atmospheric escape and significant radius inflation.Batygin et al. (2011) showed that Ohmic heating acts to inflate hot Jupiters, with planets ≲ 0.5M Jup overflowing their Roche lobes and leading to evaporation on Gyr timescales.The Ohmic heating scenario requires only a modest planetary magnetic field (≳1 G) and an equilibrium temperature great enough to thermally ionize some modest fraction of neutral metals, such as the low ionization species of Na I and K I.The "sweet spot" for this phenomenon appears to prefer equilibrium temperatures in the range of 1500 < T eq < 2000 K (Batygin et al. 2011;Menou 2012;Ginzburg & Sari 2016;Thorngren & Fortney 2018;Knierim et al. 2022), where the conductivity is strong enough to cause an effective drag without being so strong that magnetic braking slows the planetary wind.HAT-P-67 b's ∼2000 K sits to the higher end, but still in a region of high Ohmic dissipation heating efficiency.A proposed order-of-magnitude scaling law predicts an inflation timescale, Eq. 20 in Batygin et al. (2011): yielding an incredibly short <5 Myr timescale for HAT-P-67 b, assuming a typical ϵ ∼ 0.01.The order of magnitude of this inflation timescale is so fleetingly short that-according to this scenario-HAT-P-67 b must be in the runaway stage of inflation, rapidly losing mass and growing in surface area to fuel a positive feedback loop.Under this interpretation, HAT-P-67 b would represent an example of a new category of planet system that is doomed to evaporate entirely due to the Ohmic dissipation mechanism, as predicted by Batygin et al. (2011).At Roche lobe overflow, a 5 Myr inflation timescale may imply an instantaneous Ṁinfl > 10 3 M ⊕ /Gyr.A few caveats complicate the unambiguous causal interpretation of Ohmic dissipation.First, the scaling law arguments that produced Equation 2 were only proposed as coarse estimates, with numerical simulations needed to quantify inflation timescales for individual systems.Accordingly, a 10× higher inflation timescale of 50 Myr-allowed by the coarse scaling law-would yield Ṁinfl ∼ 10 14 g/s, still a very large mass loss rate, and only a factor of 5 away from the baseline 1D model.Factors of a few uncertainties in the Ohmic dissipation efficiency ϵ may also contribute.
Second, numerical simulations (Wu & Lithwick 2013) and analytic theory (Ginzburg & Sari 2016) indicate that Ohmic dissipation can stall the contraction of hot Jupiters but cannot easily "re-inflate" them after having undergone a traditional cooling curve.Heat transfer from the atmosphere to a cooled core appears to proceed too slowly, on the order of tens of Gyr (Ginzburg & Sari 2016).This path dependence of Ohmic dissipation would restrict the allowed evolutionary histories of HAT-P-67 b, to have arrived at its current location within a few to tens of Myr.This short time window prefers a physical mechanism such as disk migration, which would be faster than secular eccentric migration with an outer companion.In-situ formation of a hot Saturn at these close-in separations may be implausible (Dawson & Johnson 2018).
Other caveats like unknowns in planetary magnetic fields, zonal band geometries, and dayside/nightside temperature differences make it impossible to uniquely prescribe Ohmic dissipation, and instead, several additional factors may also be at play (Sarkis et al. 2021).

Reinflation
Evolved stars increase in luminosity, delivering greater insolation to planets at a fixed separation.The heightened equilibrium temperature can cause mature planets to inflate, the phenomenon known as reinflation.Such re-inflated hot planets around red giant stars have been recently found by TESS (Grunblatt et al. 2022(Grunblatt et al. , 2023)), though not all planets around evolved stars appear to re-inflate (Saunders et al. 2022).Zhou et al. (2017) estimated that HAT-P-67 b received about twice the incident flux as a Zero Age Main Sequence (ZAMS) HAT-P-67, based on comparison to the Geneva isochrones.In Section 3.3 MIST evolutionary tracks showed two equally plausible states: the tail-end of the main sequence or a recently evolved sub-giant.Figure 15 employs these tracks to quantify the prospects for re-inflation of HAT-P-67 b.
We further assume that the orbital location has not changed over the system lifetime, that the planet's response to anomalous heating stimulus is instantaneous, and that the anomalous heating efficiency εG peaks at T eq = 1750 K as derived by Thorngren & Fortney (2018).
We conclude that reinflation likely does not have a significant effect on the evolution of HAT-P-67 b: inflation timescales remain relatively unchanged over the planet's lifetime, despite a 50% spike in incident radiation in the subgiant scenario.The reason is subtle, as we describe next.
The second panel from the top illustrates a countervailing effect: the increase in flux actually triggers a decrease in anomalous heating efficiency, ϵ.The greater stellar energy couples into the planet less effectively than the weaker radiation did, roughly balancing out.Together the effects nearly cancel when computing the inflationary lifetimes, yielding nearly identical curves in the bottom panel: a secularly evolving main sequence history gives almost the same inflation timescale as a rapidly increasing subgiant.It is important to emphasize that the anomalous heating efficiency found by Thorngren & Fortney (2018) is agnostic to what the actual heating mechanism is: either XUV Irradiation or Ohmic dissipation both have to obey the trend of εG (T eq ).
The third panel recreates the numerically computed curve for 0.5 M Jup , T eq =1800 K from (Batygin et al. 2011), which should be considered as representative since it was not necessarily tailored to the evolutionary history of HAT-P-67 b.Nevertheless, the similarity of the curve is remarkable since it arrives at approximately the correct planet radius at the right age for about the right initial mass.

Weighing the causes for a lack of sub-Saturns
When applied to an ensemble of planet systems, both the XUV irradiation and Ohmic dissipation mechanisms expect a void in the planet mass-radius plane.However, the theories make different quantitative predictions for the placement of the dividing lines between stable and unstable populations.We illustrate these differences in Figure 16, which shades the mass-radius plane with predictions for the inflation timescale under HAT-P-67blike conditions.The XUV irradiation timescale better matches the distribution of planets, with the 0.1g cm −3 density contour setting a conspicuous dividing line in density.The Ohmic dissipation shading predicts short inflation timescales extending into a region of Jupitermass planets with radii commonly observed.The better match of observed exoplanet demographics to the XUV irradiation shading disfavors Ohmic dissipation.
Under either scenario, HAT-P-67b resides in an extremely short-lifetime region.Its nearest analog, HAT-P-32b, offers an interesting test case since it has recently been shown to exhibit a large helium excess (Zhang et al. 2023).It resides just slightly denser than the 0.1 g cm −3 dividing line, with XUV predicting over 10× longer inflationary timescale for HAT-P-32b than for HAT-P-67b, while Ohmic dissipation expects merely a factor of 3 difference.HAT-P-32b exhibits a more symmetric leading and trailing tail, whereas HAT-P-67b stands out as exhibiting evidence for preferential mass loss on the highly irradiated planet dayside.These differences may make this pair an especially valuable laboratory for developing theories of atmospheric escape.The second panel from the top shows the corresponding estimate for anomalous heating efficiency from Thorngren & Fortney (2018).The Ohmic dissipation model makes predictions for runaway inflation, depending on the incident stellar flux.The heuristic scaling relation for the inflation timescale predicts vanishingly short inflationary lifetimes in this extreme regime but illustrates the countervailing effects of the top two panels.
Finally, we consider the prospect that the planetary mass of HAT-P67b is much lower than the 1-σ estimate, such that the white-light planetary radius is the Hill radius (we have previously assumed the Hill radius is 2.7 planetary radii).This extreme scenario could manifest enormous mass loss rates, with the atmosphere's reservoir of mass fleeing the gravitational potential well directly, without the cushion of an exosphere.This terminal Roche Lobe overflow should have observational consequences.In particular, heavy elements would easily leak into the planetary wind, yielding possibly many observable metal lines in the UV.

Other planets likely to be evaporating
The key figure-of-merit can be distilled to τ infl , which can be thought of as an atmospheric escape spectroscopy metric (Kempton et al. 2018) for inflated planets in the shaded regime of Figure 16.Table 3 presents the rankordered list of planets by this metric, indicating that they are mostly smaller and more massive than HAT-P-67 b.Only KELT-11 b shows a shorter nominal inflationary timescale than HAT-P-67 b, owing to its low mass and equilibrium temperature residing closer to the peak of ϵ G (T eq ).We predict that these few other inflated HAT-P-67 b-analogs should show evidence for significant atmospheric escape, comparable to what we see for HAT-P-67 b.These sources make excellent targets for He I 10833 observations, and/or other atmospheric escape diagnostics.As emphasized, there are significant uncertainties in the inflation timescale, but at least its order-of-magnitude value gives us a quantitative and justified way to prioritize target selection.

Cause for stellar modulation
As discussed in Section 3.2.2, the TESS lightcurve exhibits up to 0.36% peak-to-valley modulation amplitude, with a characteristic timescale of P = 4.7 − 5.9 days.Some sectors exhibit lower amplitudes.The conventional interpretation would consign these cyclical modulations to the familiar stellar activity: surface featureseither starspots, faculae, or plage-entering and exiting the projected stellar disk on the stellar rotation period comparable to the 4.8-day orbital period of planet b.The cause for such orbital and rotational synchronization would then be either coincidence or secular starplanet tidal interactions.For the latter, disk migration could have naturally ceased near the co-rotation radius, leaving the planet to reside naturally near the stellar P rot .
Alternatively, Star Planet Magnetic Interaction (SPMI) could be at play (Strugarek 2018).In this scenario, the magnetic field of the planet permeates the space between the star and the planet.Magnetic perturbations propagate via planetary magnetic fields with the Alfvén speed.The planet can interact with the star if the Alfvén speed exceeds the stellar wind speed controlling the bulk motion of the intervening medium.The non-detection of variability in the Ca II H and K lines and Hα lines appears to disfavor this SPMI interpretation since SPMI could be expected to cause variations in these diagnostics.It is hypothetically possible-albeit unlikely-that mass loss could be great enough to produce variability in the wide band TESS lightcurve.The HPF spectra reveal up to 10% signal depth over a few Angstroms.The TESS bandpass barely includes 10833 Å, at a location of diminishing throughput.The Helium signal alone would manifest as a mere ∼ 4 ppm flux loss when integrated over a TESS-throughput-weighted F-star spectrum-negligible compared to the observed 0.36% peak-to-valley modulation.An ensemble of additional lines in the red-optical cannot realistically resemble the TESS modulation.The inventory of such conceivable atomic lines detectable from the planet in the TESS bandpass numbers merely a few, with Ca II infrared triplet and Hα chief among them (Linssen & Oklopčić 2023).A putative Hα line would have to be about 30% deep over 10 Å wide to manifest perceptibly in the TESS data.Such an implausibly deep and wide line likely would have been observed as perturbations in the Keck HIRES spectra, even with its limited phasing.Hypothetically, dust dredged up in the mass loss process could cause a large enough broadband continuum flux loss to manifest in TESS.We may expect to see some variable reddening in that case.

Implications for exosphere detection in non-transiting planets
We measured a large extent of Helium escape along the arc of the orbit, but we necessarily can place only coarse constraints on the vertical extent-in the direction perpendicular to the orbital plane, ±z.We estimate the detectable vertical extent must be at least a few stellar radii, such that we still could have detected Helium escape in a hypothetical HAT-P-67b-like system even if it were non-transiting, in a "near miss" configuration.We propose a new category of semi-transiting planet, dubbed "exospheric grazers", in which the planet does not produce a detectable white-light transit depth, but does produce measurable line-based exospheric absorption.This category appears to have been neglected due to the assumption that exospheres were only easily detectable at several planetary radii.While large tails have been seen previously in Lyα, the scarcity of UV resources prohibited searches of this kind.The discovery of large Helium tails in HAT-P-32 b (Zhang et al. 2023) and now HAT-P-67 b suggest that the identification of these exospheric grazers may be achievable with current instrumentation and possibly existing archival data from near-IR RV-planet searches.Non-transiting, strongly irradiated planets with well-constrained orbits may make ideal targets for the detection of this phenomenon.

CONCLUSIONS
We have presented a multi-year spectroscopic survey of HAT-P-67 b, a low-density, heavily irradiated Saturnmass (or lower) planet.We identified a large leading tail, with a much weaker trailing tail, which we found to be direct evidence of preferential dayside mass loss.HAT-P-67 b stands out as an outlier in the mass-radius plane, which we examine through the lens of different mechanisms for anomalous heating.Both XUV irradiation and Ohmic dissipation predict runaway inflation for such inflated hot Saturns, and we quantitatively weigh these two scenarios, finding XUV irradiation as more straightforwardly predictive of the overall demographic population of exoplanets observed to date.
We report in-transit line profile variability, which we attribute to the delicate interplay of planetary and stellar winds.We identify several avenues for future work, including additional monitoring of the line profile variability to probe the stellar-and-planetary wind interaction.The large signal should be perceptible in other spectral tracers, such as metal lines in the UV.We offer a list of other planets that may be likely to exhibit mass loss under the Ohmic dissipation and XUV irradiation scenarios.Based on observations obtained with the Hobby-Eberly Telescope (HET), which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximillians-Universitaet Muenchen, and Georg-August Universitaet Goettingen.The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly.(Ninan et al. 2018)

Figure 1 .
Figure 1.Overview of exoplanet demographics.The pixel bins reflect the observed density of over 5000 planets accessed from the NASA Exoplanet Archive.Detections of atmospheric escape are common among large planets with strong insolation.HAT-P-67 bis among the most inflated planets known.

Figure 2 .
Figure 2. Mass-radius trends for inflated hot sub-Saturns and hot Jupiters, with layout following Figure 2 of Thorngren & Fortney (2018) and updated with NASA Exoplanet Archive confirmed planets.The trend lines show the massradius relationship for equilibrium temperatures of 500, 1000, 1250, 1500, and 2000 K, assuming the mean composition and mean heating efficiency from the original figure.HAT-P-67 b (⋆ symbol) stands alone in a region defined by the lack of inflated sub-Saturns and explained by short lifetimes from runaway inflation and mass loss.

Figure 3 .
Figure 3. Overview of all available TESS Sectors showing 24 full or partial transits with 34 visits with HPF (vertical gray bars), 7 of which coincide with transits.The rest of the visits sample out-of-transit phases.The 8 HPF visits between days 2050 and 2690 are not shown.Note the large time breaks between some consecutive panels, which correspond roughly to the duration of a TESS sector.

Figure 4 .
Figure 4. Best fit orbit overlaid on the composite TESS lightcurve.Revised planet properties are consistent with Zhou et al. (2017), with a systematic shift from a revised Gaia DR3 distance.

aFigure 5 .
Figure 5. MIST evolutionary model tracks for HAT-P-67.The T eff and luminosity of HAT-P-67 (gray circle) are consistent with either the late main sequence or recently evolved subgiant.The labeled numbers indicate age in Gyr.

Figure 6 .
Figure6.RV orbit fit including both HPF and Keck HIRES data points.The RV information content in HPF is much less than in Keck HIRES for this F5 spectral type, and so the joint fit constrains the mass to roughly a Saturn mass.
Figure 7.The sky-subtraction, telluric masking, barycentric correction, flattening, and all of our other standard pre-processing steps were carried out in the open-source Python interface muler (Gully-Santiago et al. 2022).

Figure 7 .
Figure 7. Overlay of all 152 individual HPF exposures of HAT-P-67, spanning 2020-2022.Variability is seen in the He I 10833 Å triplet near the vertical orange shaded band but not in the adjacent Si line.These snapshot spectra were barycentric corrected and continuum flattened with a linear fit to the regions in the blue vertical bands.Sharp telluric absorption lines have been masked in regions near 10835 and 10837.5 Å.

Figure 8 .
Figure 8. Four observing campaigns centered on an intransit epoch with before-and-after visits typically separated by 1 night.The after-transit spectra tend to show negligible absorption.Line-of-sight velocity substructure can be seen in the before and during transit epochs.

Figure 9 .
Figure 9. Equivalent Width lightcurve of Helium absorption in HAT-P-67b.The system exhibits a characteristic absorption leading up to transit, followed by a sharp decline after transit passage.The EWs were computed in the orange shaded band (10832.3−10834.2Å) in Figure 7.

Figure 10 .
Figure 10.Absorption depth phase scan over the entire HAT-P-67b orbit from 41 visits with HPF.The planet orbital rest frame is shown as the black sine wave.Gaps in data at 10834-10836 Å arise from telluric masking.Unmasked telluric lines are perceptible at 10825 Å.

Figure 11 .
Figure 11.Line centroid positions of the He I 10833 Å feature.A slight blueshift relative to the stellar restframe (vertical line) can be seen, with a weak excess blueshift trend after planet egress.The line-of-sight velocity distribution of 2-3 Å means that the gas distribution exhibits both advancing and receding lobes at all phases other than post-egress.

Figure 12 .
Figure 12.Synthetic XUV radiation scenarios for HAT-P-67 b.The low and high XUV scenarios correspond to LX /L bol of 10 −6 and 10 −4 , where LX is defined as He 1 1 S Ionizing photons that trigger the recombination cascade needed for observing the He 2 3 S ionization metastable state.He 2 3 S Ionizing photons depopulate the state and suppress the observability of He I 10833 Å.The XUV flux of HAT-P-67 is uncertain.

Figure 13 .
Figure 13.Simulated 1D Parker wind model of Helium absorption in HAT-P-67 b.The spectrum was generated with a mass loss rate of 2×10 13 g/s, exosphere gas temperature of 14 000 K, and with an XUV luminosity Lx/L bol = 10 −5 intermediate between those shown in Figure12.The p-winds model yields the He I 10833 triplet lines (colored thin lines), with a cumulative feature in black resembling the data, shown as the May 2020 in-transit mean.Additional 3D velocity dispersion of the gas can explain differences between the line shapes of the model and data.Overall the 1D models may be too simplistic to represent the inherently 3D structure of the outflow.

Figure 14 .
Figure 14.Schematic of Keplerian orbital shear.Locations outside the Roche lobe experience orbital shear from the Keplerian potential.A wind launched primarily from the dayside will tend to form a leading tail.

Figure 15 .
Figure 15.Conceivable evolutionary scenarios for HAT-P-67 b.The planet's equilibrium temperature increases as the stellar luminosity gradually increases on the main sequence.The second panel from the top shows the corresponding estimate for anomalous heating efficiency fromThorngren & Fortney (2018).The Ohmic dissipation model makes predictions for runaway inflation, depending on the incident stellar flux.The heuristic scaling relation for the inflation timescale predicts vanishingly short inflationary lifetimes in this extreme regime but illustrates the countervailing effects of the top two panels.

Figure 16 .
Figure 16.Runaway inflation timescale in the exoplanet mass-radius diagram.The shading of the left panel shows τ infl from Ohmic dissipation; the right panel shows τ infl ≡ M/ Ṁ from photoionization-driven mass loss.Individual planets with confident mass and radius detections are shown as small black dots.Helium non-detections are shown in red squares, atmospheric escape detections from any origin in green circles (Dos Santos 2022).Under either scenario, HAT-P-67 b resides in a sparsely populated region with a short inflationary timescale, making it unstable to runaway evaporation with large mass loss expected.Photoionization better predicts the depopulation of sources in the upper left low-density region of the diagram, defined by the 0.1 g cm −3 iso-density dashed white line.

Figure 17 .
Figure 17.Non-detection of variability in the Calcium IR Triplet line at 8862 Å.No other HPF lines appear to show any variability.
based on observations obtained with the Habitable-zone Planet Finder Spectrograph on the HET.The HPF team acknowledges support from NSF grants AST-1006676, AST-1126413, AST-1310885, AST-1517592, AST-1310875, ATI 2009889, ATI-2009982, AST-2108512, AST-2108801, and the NASA Astrobiology Institute (NNA09DA76A) in the pursuit of precision radial velocities in the NIR.The HPF team also acknowledges support from the Heising-Simons Foundation via grant 2017-0494.
Exoplanets and Habitable Worlds is supported by the Pennsylvania State University and the Eberly College of Science. 17 18GS acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51519.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI).Funding for the TESS mission is provided by the NASA Explorer Program.STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration.
made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.made use of NASA's Astrophysics Data System.

Table 2 .
Revised orbital and planetary parameters