Exoplanet Aeronomy: A Case Study of WASP-69b's Variable Thermosphere

Aeronomy, the study of Earth's upper atmosphere and its interaction with the local space environment, has long traced changes in the thermospheres of Earth and other solar system planets to solar variability in the X-ray and extreme ultraviolet (collectively,"XUV") bands. Extending comparative aeronomy to the short-period extrasolar planets may illuminate whether stellar XUV irradiation powers atmospheric outflows that change planetary radii on astronomical timescales. In recent years, near-infrared transit spectroscopy of metastable HeI has been a prolific tracer of high-altitude planetary gas. We present a case study of exoplanet aeronomy using metastable HeI transit observations from Palomar/WIRC and follow-up high-energy data from the Neil Gehrels Swift Observatory that were taken within one month of the WASP-69 system, a K-type main sequence star with a well-studied hot Jupiter companion. Supplemented by archival data, we find that WASP-69's X-ray flux in 2023 was less than 50% of what was recorded in 2016 and that the metastable HeI absorption from WASP-69b was lower in 2023 versus past epochs from 2017-2019. Via atmospheric modeling, we show that this time-variable metastable HeI signal is in the expected direction given the observed change in stellar XUV, possibly stemming from WASP-69's magnetic activity cycle. Our results underscore the ability of multi-epoch, multi-wavelength observations to paint a cohesive picture of the interaction between an exoplanet's atmosphere and its host star.

The metastable HeI signal depends sensitively on the high-energy radiation incident on an exoplanet's upper atmosphere.Not only does XUV power photoevaporative outflows (Owen 2019), but the relative amount of XUV to mid-UV also sets the level populations and ionization environments from which the metastable HeI feature is produced (Oklopčić 2019).Therefore, connecting the growing body of metastable HeI data to astrophysical models of exoplanet mass-loss requires XUV measurements of host stars.Extreme ultraviolet (EUV) radiation is undetectable by current facilities, so X-ray fluxes are often used to infer total XUV outputs (Sanz-Forcada et al. 2011;Chadney et al. 2015;Nortmann et al. 2018;King et al. 2018).When X-ray data have been unavailable for systems with planetary metastable HeI measurements, previous studies have used proxy spectra from stars with similar masses and activity indices (Mansfield et al. 2018;Vissapragada et al. 2022a;Bennett et al. 2023).
Nonetheless, different stars of the same mass can vary by an order of magnitude in X-ray flux (Messina et al. 2003;Jackson et al. 2012;Wright et al. 2018;Johnstone et al. 2021); exoplanet mass-loss rates derived with proxy stellar spectra from another star are similarly uncertain (Behr et al. 2023).Compounding the issue, the X-ray output from a given star is often time-variable (Wagner 1988;Woods & Rottman 2005;Chadney et al. 2015) and should alter exoplanetary atmospheric conditions accordingly (Wang & Dai 2021).Thus, singleepoch snapshots of exoplanet thermospheres likely do not reflect the time-averaged values.
Due to this star-to-star variability and time-variability of individual stars' XUV production, a complete understanding of exoplanet aeronomy requires XUV and metastable HeI observations of the same star-planet system across multiple epochs.Individual stellar XUV and planetary thermospheric data are most complementary when obtained simultaneously.In practice, telescope scheduling constraints and transit ephemerides often prohibit this ideal case.Since spot-driven stellar XUV fluxes evolve on timescales a few stellar rotation periods (Woods & Rottman 2002), we define "contemporaneous" multi-wavelength observations that occur within such a window as being particularly helpful for the aeronomy science case.
Recent research has pushed towards this multiwavelength framework, such as a study of HD 189733 by Zhang et al. (2022).That work found 30% changes in metastable HeI absorption by the companion hot Jupiter over years-long timescales, which could be explained by the general level of the host's X-ray variability (Pillitteri et al. 2022).However, the high-energy flux timeseries was obtained before any metastable HeI observations.In other multi-wavelength work, Nortmann et al. (2018), Czesla et al. (2022), and Zhang et al. (2023) collected metastable HeI transmission spectroscopy and X-ray fluxes for various systems with varying degrees of contemporaneity.
Here, we present new metastable HeI results from a transit of WASP-69b viewed by the Wide Field In-fraRed Camera (WIRC; Wilson et al. 2003) at Palomar Observatory and high-energy observations from the Neil Gehrels Swift Observatory (Swift) that were taken within a one-month window.Combining these measurements with archival metastable HeI (Nortmann et al. 2018;Vissapragada et al. 2020Vissapragada et al. , 2022b;;Allart et al. 2023;Tyler et al. 2024;Guilluy et al. 2024) and highenergy (Nortmann et al. 2018;Spinelli et al. 2023b) observations provides a powerful assessment of WASP-69b's outflow variability.The bright (J = 8.0 mag) host star and large scale height of the hot Jupiter companion (0.26 M J , 1.06 R J ; Anderson et al. 2014) make the WASP-69 system an exceptional target for studies seeking to quantify the time-variability of atmospheric outflows.
Metastable HeI variability was first reported for WASP-69b by Nortmann et al. (2018), who uncovered changes in excess absorption within the same month.Then, Allart et al. (2023) analyzed archival SPIRou data from 2019 and found a lower metastable HeI signal compared to either CARMENES dataset.Later, Tyler et al. (2024) documented changes in post-egress absorption by comparing their Keck/NIRSPEC spectra to CARMENES ones from Nortmann et al. (2018).Most recently, Guilluy et al. (2024) presented data from four transits with metastable HeI absorption and simultaneous optical spectroscopy.That study reported and analyzed metastable HeI variability in the context of simultaneous Hα measurements.We also take a multiwavelength approach in our study, but we instead complement our metastable HeI data with observations of stellar XUV flux to connect our research with the geoscientific discipline of aeronomy.
In Section 2, we describe our new metastable HeI transit for WASP-69b.Section 3 presents individual light curve fits for both new and archival WIRC observations, a joint fit of WIRC and TESS data, and a joint fit of WIRC, TESS, and archival spectroscopic data.Then, we analyze the new Swift observations of WASP-69 in Section 4. With multi-wavelength data, we model the stellar XUV spectrum and hot Jupiter's mass-loss rate Table 1.Summary of Palomar/WIRC observations of WASP-69b."Start" and "Finish" correspond to the beginning and end of the science image sequence, respectively."Moon Frac." and "Moon Sep." refer to the illuminated fraction of the Moon and angular on-sky separation between the moon and WASP-69b, respectively.The transit midpoint uncertainty σ mid is calculated from the ephemeris from Kokori et al. (2022), while zst, zmin, and z end represent the airmasses at the beginning, the minimum, and at the end of the science imaging, respectively.The symbols nexp, ntr, texp denote the number of science images, the number of in-transit images, and the exposure time of the science images, respectively.We observed a full transit of WASP-69b with WIRC on the night of 2023 August 25 UT.A summary of the observations for the new light curve and for an archival WIRC light curve from Vissapragada et al. (2020) are provided in Table 1.Conditions were excellent on both observing nights; the humidity was 41% and 9% during the new and archival observations, respectively.Before collecting the new light curve, we obtained dark frames, flat-fielding, and HeI arc lamp reference images.Then, we constructed a background frame with a fourpoint dither immediately prior to taking science data.Our observing setup was identical to that of Vissapragada et al. (2020), permitting a precise comparison between new and archival transit observations.Science images were obtained with the beam-shaping diffuser (fore wheel) to reshape stellar point spread functions (PSFs) into 3. ′′ full-width at half-maximum (FWHM) top-hats (Stefansson et al. 2017) and the ultra-narrowband HeI filter (aft wheel).

Night
We positioned the centroid of WASP-69's PSF within 7 pixels (1.′′ 75) of where the star was placed for the archival WIRC observations.This consistency is particularly important because HeI filter observations have position-dependent wavelength sensitivity (Vissapragada et al. 2020).Imaging was done with an exposure time of 60s, the same integration cadence as Night 1. Guiding was stable throughout the observations.We captured the entire transit as well as approximately 40 and 120 minutes of pre-ingress and post-egress baseline, respectively.For comparison, the archival WIRC observations contain nearly equal pre-ingress and post-ingress baselines (about 120 minutes for each).

Light Curve Data Reduction
We reduced data from the archival (Night 1) and new (Night 2) observations with the pipeline implemented by Vissapragada et al. (2022b).Raw science images were dark-corrected, flat-fielded, and cleaned of detector cosmetic imperfections.Then, we corrected for telluric emission lines by exploiting the fact that the ultra-narrowband filter has position-dependent absorption.Following the procedure from Vissapragada et al. (2020), we removed sources in the dither sequence images via sigma clipping and median-scaled the combined dither image for each science frame in 10 pixel radial steps from the point on the filter where the incident light rays are normal to the filter surface.In the light curve modeling, we used the resultant scaling factors from this procedure to detrend on the time-varying telluric water absorption (Paragas et al. 2021).
Next, we extracted WASP-69's light curve with aperture photometry.First, we selected comparison stars on which to perform differential photometry.Four and two comparison stars in the field were detected beyond 50σ significance on Nights 1 and 2, respectively.Previously published Night 1 reductions by Vissapragada et al. (2020) and Vissapragada et al. (2022b) used four comparison stars that were detected beyond 50σ significance.Even though the Night 2 data had the same exposure time, only two of those comparison stars met that threshold on Night 2; we attribute this result to minor differences in seeing between the nights (Figure 1) that led to a slightly larger stellar PSF size on Night 2. For all analyses in this paper, we elected to keep only the two comparison stars that were detected beyond 50σ significance on both nights.We found, however, that selecting two versus four comparison stars on Night 1 did not meaningfully change the derived transit model in Section 3.
We determined the optimal aperture size by testing apertures with radii ranging between 3 − 20 pixels (0. ′′ 75 − 5. ′′ 00) in radius.For each candidate size, we constructed WASP-69's light curve and detrended by the mean of the comparison star fluxes.Local background subtraction around the stars was done with an annulus of 25 − 50 pixels (6.′′ 25 − 12. ′′ 50).We then chose the aperture that minimized the root-mean-square (rms) scatter after clipping 5σ outliers with a moving median filter.We found optimal aperture sizes of 11 and 12 pixels for Nights 1 and 2, respectively.To ensure that the stellar PSF was stable and consistent during the imaging, we examined the Gaussian FWHM of WASP-69's radial seeing profile for each science exposure.Specifically, we constructed seeing profiles for each image from 100 uniformly-spaced radii between 2-40 pixels and calculated the FWHM with photutils.Figure 1 plots the FWHM of WASP-69's PSF versus WASP-69b's orbital phase, assuming the ephemeris published by Kokori et al. (2022).
Both PSFs were stable across their respective light curves, likely due to the excellent and consistent weather conditions within each night.Although the Night 2 PSF was slightly larger on average than the Night 1 PSF, the FWHM of the Night 2 PSF still remained within 1 pixel of the median value for most exposures.The smaller average FWHM of the Night 1 PSF matches the fact that we found a smaller best-fit aperture (11 pixels) for this timeseries than for the Night 2 PSF (12 pixels).

MULTI-EPOCH TRANSIT FITTING
3.1.Palomar/WIRC: Individual Nights With these reduced light curves, we first fit each WIRC dataset independently.Model setups for each night were identical, yielding a direct comparison between the best-fit transit depths.
Our procedure was similar to the one described by Vissapragada et al. (2022b), which uses exoplanet to fit a starrygenerated light curve to the data.Our procedure simultaneously fit a light curve model with a model of the systematics.Free parameters in the light curve model were the planet-to-star radius ratio R p /R ⋆ , the exoplanet's orbital period P , the reference transit midpoint t 0 , the scaled semi-major axis a/R ⋆ , the impact parameter b, a jitter j to account for the difference between Poisson noise and the realized scatter, and weights of the comparison stars and covariates for linear systematics detrending.We fixed the orbital eccentricity to zero, in agreement with published radial velocity constraints on WASP-69b that were consistent with a circular orbit (e < 0.11; Bonomo et al. 2017).
We fixed limb-darkening coefficients to values that were calculated with ldtk for WASP-69 by Vissapragada et al. (2022b): [u 1 , u 2 ] = [0.38,0.12].Using these constant values, we determined an optimal set of covariates on which to detrend in the final fits.We considered airmass, the centroid offsets, and the water absorption proxy introduced by Paragas et al. (2021).For each combination of these covariates, we optimized the light curve model parameters to find the best-fit solution and removed outlier points more than 4σ from the best-fit solution.We iterated this optimization and sigma-clipping until no 4σ outliers remained.Then, we sampled the posterior with the No U-Turn Sampler (NUTS) in pymc3 (Hoffman & Gelman 2011) using four chains with 1500 tuning steps and 3000 draws.We set a target acceptance fraction of 0.99, a high value that is necessary to properly sample topologically complex posteriors.
All priors were informed by the ephemeris from Kokori et al. (2022) and planetary data from Anderson et al. (2014).Except for the following changes, priors for astrophysical parameters were the same as the ones used for the joint fit in Section 3.2 and provided in Table 2. Instead of the wide priors for b, a/R * , P , and t 0 that are listed in Table 2, we used tighter ones based on literature values selected by Vissapragada et al. (2022b): b = N (0.686, 0.023), a/R * = N (12.00,0.46), P = N (3.8681390, 0.0000006), and t 0 (BJD−2450000) = N (7176.17789,0.00017).Reference star weights and covariate weights were all assigned uniform U(−2, 2) priors.We used these prior values to check our results versus the analysis of Vissapragada et al. (2022b); we later widened these priors when doing joint fits as another robustness check.
After finding the best-fit light curve for a given set of covariates, we calculated the Bayesian Information Criterion (BIC) for that model.Once we obtained BICs for all detrending possibilities, we selected the covariate set for our final fit that provided the best BIC.To fairly compare sets of covariates, we kept a constant jitter term in the fits for covariate selection and ensured that the sigma-clipping procedure in the light curve fit yielded timeseries of the same length.We validated our covariate selection by ensuring that neither the Watanabe-Akaike Information Criterion nor the Pareto-Smoothed Importance Sampling Leave-One-Out (PSIS-LOO) statistic preferred different covariates within their arviz-reported (Kumar et al. 2019) 1σ uncertainties.
With our preferred covariates chosen -water proxy and airmass for Night 1 and airmass only for Night 2 -we did a final fit with free values of the limbdarkening coefficients and jitter.Although fixing the limb-darkening coefficients could result in smaller uncertainties on the planet-to-star radii ratio, we decided to remove this external model dependence from our light curve fits.We optimized the solution and sampled the posterior with the same sampling hyperparameters as the aforementioned covariate-finding methodology.Visual inspection of the posterior corner plot indicated that our chains had converged, which we statistically validated by finding a Gelman-Rubin diagnostic R < 1.01 for each sampled parameter.
We found best-fit (R p /R ⋆ ) = 0.1470 +0.0027 −0.0025 for Night 1 and (R p /R ⋆ ) = 0.1360 +0.0032 −0.0032 for Night 2. These depths differ by 2.6σ.We confirmed that our results for all derived parameters of the WIRC Night 1 light curve were consistent with the result reported by Vissapragada et al. (2022b).In Appendix A, we analyze how each covariate affected the measured transit depth, assess the degree of correlated noise, and check the fit robustness to changes in aperture sizes and out-of-transit baseline.Comparing the night-to-night results is suggestive of astrophysical variability in WASP-69b's metastable HeI signal, but this procedure cannot rule out the possibility that different best-fit transit shapes masqueraded as two different transit depths.Each timeseries was fit independently, albeit with identical procedures, but a more rigorous test would be to fit the data jointly.

TESS + WIRC Joint Fitting
Towards our goal of determining whether WASP-69b's thermosphere is time-variable, we implemented a simultaneous fit of the two WIRC light curves along with the light curve from the Transiting Exoplanet Survey Satellite (TESS ; Ricker et al. 2015).This space-based Table 2. Summary of input priors for fitting data from individual nights (Section 3.1) and joint fitting (Section 3.2, 3.3) along with the resultant posteriors from joint fitting WIRC Night 1, WIRC Night 2, and TESS Sector 55.Priors on the ephemeris and planetary parameters were informed by Kokori et al. (2022) and Anderson et al. (2014), respectively.On the planet-to-star radii ratios (Rp/R * ), the subscripts T, (W,1), and (W,2) correspond to the TESS, WIRC Night 1, and WIRC Night 2 datasets, respectively.Quadratic limbdarkening coefficients are given by u1 and u2, with subscripts T and W for TESS and WIRC, respectively.Priors for the limb-darkening coefficients were implemented using the methodology by Kipping (2013) via the exoplanet package.The jW,1, jW,1, and ϵT are the jitter terms for WIRC Night 1, jitter term for WIRC Night 2, and the TESS error scaling parameter. .By fitting the WIRC data with TESS broadband photometry, we can self-consistently constrain the transit timing and shape.Thus, this procedure interrogates differences in the WIRC transit depths while controlling for a cause that would be unrelated to the planetary thermosphere.Importantly, our goal in using the TESS photometry is not to provide a reference by which to compute excess metastable HeI absorption; other bandpasses that are closer to the 10830 Å feature are more appropriate for that task (see Section 5).

Parameter
We downloaded the TESS 20s cadence Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; STScI/MAST 2021) times, fluxes, and flux uncertainties from the mission's Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) via lightkurve, which queries the Mikulski Archive for Space Telescopes (MAST) server.While seven transits were expected to occur within the sector based on the ephemeris published by Kokori et al. (2022), two of those events occurred during gaps in data coverage.
The PDCSAP photometry exhibits secular trends which are likely attributed to stellar variability and/or instrumental effects.To detrend photometry while preserving transits, we fit a spline with break point spacing of 5 days to the PDCSAP points with Keplerspline v2 code1 (Vanderburg & Johnson 2014;Shallue & Vanderburg 2018).Dividing the original data by the spline trend flattened the light curve.We illustrate these preprocessing steps in Figure 2. WASP-69b's transits are deeper than any stellar variability or instrumental systematics, so we did not sigma-clip outliers until the jointfitting itself.To reduce the computational burden of the joint fit, we filtered to keep only points that were within 0.15 d of any observed transit midpoints.
In our joint model, we required that all three light curves share the same orbital period P , transit midpoint reference time t 0 , scaled semi-major axis a/R * , impact parameter b, and stellar radius R * .We allowed different quadratic limb-darkening coefficients between the TESS (u 1,T , u 2,T ) and WIRC (u 1,W , u 2,W ) bandpasses.To assess the time-variability of the metastable HeI signal, we independently fit for the planet-to-star radii ratio in the TESS bandpass (R p /R ⋆ ) T , the first WIRC night (R p /R ⋆ ) W,1 , and the second WIRC night (R p /R ⋆ ) W,2 .Each WIRC time series could have different values for jitter, weights for comparison stars, and weights for shared covariates.We also included an errorscaling parameter for the TESS photometry to account for deviations from pure photon noise; this term is similar in function to the WIRC jitter terms.Finally, we kept the same optimal sets of covariates for each night as were determined in Section 3.1.
When comparing ultra-narrowband photometry from multiple epochs, transit depth variability can emerge from the changing topocentric radial velocity of the target star.If the metastable HeI triplet falls on parts of the bandpass with different transmission properties in two different observations, then the recorded transit depths could differ even if the planetary spectrum is unchanged.This issue was negligible between our two WIRC timeseries, as the difference in WASP-69's topocentric radial velocity was less than 5 km s −1 at the respective transit midpoints.The 10830 Å wavelength was shifted in the observatory frame by 0.2 Å, much smaller than the WIRC filter's FWHM of 6.35 Å.Nevertheless, we accounted for this shift when we analyzed variability among a larger set of observations in Section 3.3.
With this model setup, we first removed 10σ outliers from the reduced WIRC photometry.Then, we maximized the log-likelihood of the joint model with pymc3 and executed sequential rounds of sigma-clipping to eliminate points from both the WIRC and TESS timeseries whose residuals were at least 5σ discrepant from the maximum a posteriori (MAP) solution.When no 5σ outliers remained, we re-optimized the joint model.Finally, we sampled the posteriors with NUTS to obtain the posterior probability distributions.We ran 4 chains with 1000 tuning steps and 5000 draws.As in Section 3.1, we set a target acceptance fraction of 0.99 and verified that R < 1.01 for each parameter in the posteriors.
Derived parameters of the joint fit are presented in Table 2 along with best-fit light curves and residual plots in Figure 3.The Allan deviation plots in Figure 4 show that the residuals do not contain significant correlated noise.Another indicator of fit robustness is that the posterior distributions of parameters where we used wider priors versus the ones from Section 3.1 -namely b, a/R * , P , and t 0 -converged in these fits to being consistent with those literature-derived tighter priors.Although the WIRC transit depths would be 3.21σ discrepant if their posteriors were uncorrelated, inspecting the samples drawn from the distribution (Figure 5) reveals that the (R p /R ⋆ ) W values are not independent.To account for this positive correlation, we computed the covariance matrix between the vectors of posterior samples for the (R p /R * ) values and calculated the variance while including the off-diagonal terms.With this procedure, we found that the actual difference is 3.96σ.This result reinforces our conclusion from analyzing each night separately: WASP-69b's metastable HeI absorption signal differed between the two WIRC-observed transits.
The correlation between the WIRC planet-to-star radius ratios is explained by our requirement that the transit shape remain constant between the two transits.Changes in a/R * , b, or the limb-darkening coefficients that affect one transit depth affect the other one in the same direction.The naive assumption of independent error bars underestimates the true statistical significance of transit depth variability between the WIRC nights, and the joint model reveals more significant variability in the metastable HeI absorption.

Assessing Variability With Spectroscopic Data
In addition to the ultra-narrowband photometry from WIRC, WASP-69b's metastable HeI signal has been probed spectroscopically by numerous works.Two transits were observed using CARMENES (Nortmann et al. 2018), one transit was observed using SPIRou (Allart et al. 2023), one transit was observed using Keck/NIRSPEC (Tyler et al. 2024), and four transits were observed using GIANO-B (Guilluy et al. 2024).These spectroscopic timeseries2 provide further data through which variability can be probed.We therefore implemented another joint model, this time with these spectroscopic metastable HeI data, the WIRC data, and TESS photometry.Since Guilluy et al. (2024) noted that the transit from 2021 October 28 UT was contaminated by telluric OH emission, we excluded that dataset from our analysis; thus, we only examined three GIANO-B transits.
Unlike spectroscopic data with resolution R ≳ 40k, the WIRC HeI filter has an effective R ∼ 1700 and can only detect the underlying triplet convolved with the filter's transmission profile (Vissapragada et al. 2020).Thus, comparing spectroscopic data and WIRC light curves requires translating the velocity-resolved observations into simulated WIRC-like photometry (Vissapragada et al. 2020).Alternatively, we could have attempted to translate the WIRC observations into a velocity-resolved profile, but this procedure would have required undue assumptions on the underlying spectral feature.
WASP-69's topocentric radial velocity among the spectroscopic datasets differs by more than 30 km s −1 , which would Doppler shift the wavelengths on the narrowband filter by more than 1 Å.To eliminate this non-astrophysical cause of transit depth variability, we Doppler shifted all spectroscopic data as though the observations were being conducted from WIRC on Night 1. Shifting to Night 2 instead does not tangibly affect the final fit since the WIRC radial velocities were close.
Starting from each of the reduced spectroscopic timeseries, we calculated the excess absorption f ex for each spectrum in the WIRC bandpass using Eq. 3 from Vissapragada et al. ( 2020): where T W,λ is the WIRC ultranarrowband filter's transmission and f ex,λ is the excess absorption, both at wavelength λ.Both integrals in Equation 1 were computed over the same wavelength range.Then, we constructed WIRC-like light curves by combining the WIRC-like excess depth timeseries and a ref-erence transit depth of (R p /R * ) ref = 0.1286 ± 0.0002 in the 11108-11416 Å bandpass from the Hubble Space Telescope's (HST 's) Wide Field Camera 3 presented by Tsiaras et al. (2018).These wavelengths are closer to the 10830 Å triplet than the TESS bandpass, making HST a more appropriate reference.This benchmark transit depth was also used for WASP-69b by Vissapragada et al. (2020).Although planetary absorption from water is present in this HST bandpass, this spectral feature is small in magnitude versus our transit depth precision and does not affect our analysis (Vissapragada et al. 2022b).
We estimated uncertainties on the WIRC-like light curve points through a Monte Carlo procedure, where we generated 1000 trial excess absorption timeseries for each spectral timeseries.In each trial, we perturbed each pixel-by-pixel normalized flux by a Gaussian centered on zero and with standard deviation of the uncertainty for that pixel.Pixel-level uncertainties were only provided with the CARMENES and SPIRou data in each spectrum, so we estimated the uncertainty in each spectrum from the Keck/NIRSPEC and GIANO-B timeseries by calculating the rms scatter of the pixels from 10820-10830 Å and 10836-10839.5Å, respectively.
Next, we added these WIRC-like data into another joint fit model based on the procedure from Section 3.2.We required that all light curves -TESS, WIRC, and WIRC-normalized spectroscopic data -share global parameters of the orbit: period P , impact parameter b, and transit midpoint reference time t 0 .All metastable HeI light curves were required to use the same limbdarkening coefficients, but the planet-to-star radius ratio for each metastable HeI transit was allowed to vary independently.Each excess depth timeseries from the spectroscopic observations was given a free jitter parameter but was otherwise not detrended.For computational efficiency, we binned the TESS data from the original 20 s cadence into 200 s intervals for only this fit; this coarser broadband timeseries reduced the posterior sampling runtime.
To fit the light curve model, we followed the same iterative optimization, sigma-clipping, and posterior sampling procedures as described in Section 3.2 and treated all metastable HeI data as WIRC-like light curves.During this procedure, we noticed that the post-egress tails of escaping gas reported from Keck/NIRSPEC by Tyler et al. (2024) and from GIANO-B by Guilluy et al. (2024) may have been detectable by WIRC.The light curve model that we implemented fits a symmetric transit shape, so including post-midpoint data could have spuriously depressed the best-fit transit depth.To treat all spectroscopic data fairly and identically, we therefore Finally, we sampled the posterior using NUTS with the same procedure as Section 3.2.Table 3 summarizes the results from this joint model, and Figure 6 (bottom right panel) shows a timeline of the WIRC-normalized transit depths.Light curves and residuals for all metastable HeI timeseries are located in Appendix B. We found that each of the actual WIRC transit depths were consistent within the 1σ uncertainties of the transit depth results for those timeseries from the joint fit with only TESS data in Section 3.2.A corner plot with detailed information on the fits is provided in Appendix B. In the posterior, the two most discrepant metastable HeI transit depths are the two WIRC transit depths (4.0σ).The first and second CARMENES transit depths are 2.8σ and 1.3σ discrepant with the transit depth from WIRC Night 2, respectively.Taken together, this comprehensive light curve model with archival spectroscopic data on the planet demonstrates that its metastable HeI signal is variable over timescales of several years.

X-Ray Telescope Data
To understand the space weather environment in which the August 2023 metastable HeI signal was produced, we collected high-energy data on WASP-69 with the Neil Gehrels Swift Observatory (target ID: 34059).This satellite, previously known as the Swift Gamma-Ray Burst Explorer, is a NASA spacecraft that was launched into low-Earth orbit in 2004 with the primary purpose of studying gamma-ray bursts and their afterglows (Gehrels et al. 2004).Although exoplanet astronomy was only emerging during the mission development, Swift's ability to measure high-energy radiation has since contributed to the understanding of host stars (Bourrier et al. 2020;Becker et al. 2020;Spinelli et al. 2023a) and the escape of short-period exoplanet atmospheres (Salz et al. 2019;Corrales et al. 2021).
We collected Swift data on WASP-69 between 2023 September 13 UT and 2023 September 19 UT as a Target of Opportunity (ToO; as outlined in Table 4) with 11.9ks of total integration time.This timeframe is contemporaneous with the WIRC Night 2 observations of metastable HeI by our definition from Section 1 since WASP-69b's rotation period is approximately 23 d (Bonomo et al. 2017).All observations consisted of simultaneous observations with the X-Ray Telescope (XRT; Burrows et al. 2005) and monitoring with the Ultra-Violet Optical Telescope (UVOT; Roming et al. 2005) in the UVM2 filter.Swift had previously visited WASP-69 once in 2015 for approximately 2ks (observation ID: 00034059001).That short integration cannot provide meaningful constraints in the context of our scientific objective, so we did not analyze those archival data.
Upon creating a stacked XRT image of all five Swift visits from 2023, we did not see any readily apparent sources by-eye.There was a small cluster of counts within a 20" radius near the location of WASP-69, indicating a possible marginal detection.To check our intuition, we ran a source detection algorithm from Ximage in HEASoft (version 6.28).We found that WASP-69 was detected at 2σ confidence but not at the 3σ level, reinforcing the notion of a marginal detection.Given the marginal detection, the most robust constraint on WASP-69's X-ray flux requires an upper limit analysis.Using the astropy.statsmodule, we calculated a Poissonian confidence interval based on the Bayesian approach of Kraft et al. (1991).This yielded a 95% confidence interval on the number of source counts to be [0-11.5],with the upper end of this interval corresponding to a count rate of 0.9 ks −1 .To convert this result to a X-ray flux from WASP-69, we scaled the twotemperature APEC model fit to the 2016 XMM-Newton observation by Spinelli et al. (2023b) so that the model Swift count rate indicated by Xspec matched our statistical upper limit for the measured count rate.While X-ray spectral shapes can change with varying flux, especially during flares, we did not detect any flares (see Section 4.2).We found the 95% upper limit on the 0.3-2.4keV unabsorbed flux to be 2.6×10 −14 erg cm −2 s −1 .For comparison, the best-fit model published by Spinelli et al. (2023b) from XMM-Newton observations in 2016 predicts a flux of 5.19 +0.23  −0.17 × 10 −14 erg cm −2 s −1 in the same 0.3-2.4keV energy range.The archival XMM-Newton flux is two times larger than the Swift XRT upper limit from 2023, and we show these data on the timeline in Figure 6.

Ultraviolet Optical Telescope Results
We examined the Swift/UVOT light curve of WASP-69 for flares.All data were obtained with the UVM2 filter, which covers the 200-250 nm wavelength range.This bandpass is particularly important because the ionization energy of metastable HeI 2 3 S 1 is 4.767 eV, mak-ing wavelengths close to and shortward of 260 nm efficient at removing metastable HeI through photoionization.Since 2023 August, a malfunction in one of the Swift inertial reference units has degraded its attitude control, introducing a source of instrumental jitter.A visual inspection of the UVOT images demonstrated that the pointing instability did not exceed a few arcseconds.We followed standard UVOT data reduction procedures to stack images and ran uvotsource to do aperture photometry.To account for the increased apparent size of the UVOT image point sources caused by the degraded attitude control, we increased the photometric aperture region to have a 15 ′′ radius with a 20 − 35 ′′ background annulus.
We found that WASP-69 has an apparent magnitude of 18.16 ± 0.02 in the UVM2 band, equivalent to an apparent flux of 1.20(±0.02)× 10 −15 erg s −1 cm −2 A −1 .We used uvotmaghist to extract an NUV light curve for each UVOT observation, using the same photometric aperture described above.No significant NUV flaring activity was detected.The consolidated light curve hints at a slight downward trend in the NUV flux over the course of our observing campaign, but the change in flux is around 5%, similar in magnitude to the error estimates for individual flux values in the uvotmaghist light cruve.Thus, the observed light curve is relatively consistent with a constant NUV flux.With contemporaneous data from 2023 on WASP-69's X-ray output and WASP-69b's metastable HeI transmission, we can model WASP-69b's upper atmosphere without fully relying on proxy methods to estimate the stellar irradiation.Then, we can assess whether the magnitude of the observed metastable HeI variability is consistent with the magnitude of the stellar XUV variability.We used p-winds v1.4.5 (dos Santos et al. 2022) for this exercise, which determines the HeI level populations with an input incident stellar spectrum.Metastable HeI is populated when the ground-state is ionized by EUV photons (λ < 504 Å) and subsequently recombines into the metastable state.Metastable HeI is depleted both via collisions with electrons and via ionization by mid-UV (λ < 2600 Å) photons.For giant planets orbiting K-type stars like WASP-69b, electron collisions are believed to dominate (Oklopčić 2019).

XUV Spectrum Estimation
We estimated WASP-69's XUV spectral energy distribution in four physically-motivated ranges: • F x : photons from 5-41 Å observed by Swift.
• F euv,2 : photons from 504-911 Å that cannot ionize the HeI singlet but do ionize ground-state H and the metastable HeI triplet.
• F uv : photons from 911-2593 Å that cannot ionize the HeI singlet or ground-state H but can ionize the metastable HeI triplet.
To estimate the total flux in the latter three bands, we determined and applied power-law scaling relations based on solar XUV measurements from TIMED/SEE mission (Woods & Rottman 2005), a methodology previously applied by Chadney et al. (2015) and King et al. (2018).Notably, we used scaling relations based on the Sun to infer XUV bandpass fluxes for WASP-69, a Ktype star.Chadney et al. (2015) and King et al. (2018) previously considered the applicability of extrapolating these power-laws to other stellar masses, finding that the Sun's F euv /F x versus F x correlation predicted values that were consistent with observations of active Mdwarfs.This evidence supports our assumption of using these power-laws to infer WASP-69's EUV flux.
For this calculation, we downloaded the Version 12 file of the mission's merged Level 3 data products that provides the daily-averaged solar flux observed by the satellite from low-Earth orbit in 10 Å bins over the range 0-1950 Å. Notably, reported fluxes have been slightly modified since the Version 11 release used by Chadney et al. (2015) and King et al. (2018) due to an adjustment in the detector's estimated degradation.From these data, we filtered to only the most recent complete solar cycle: 2008 December 01 through 2019 December 01.
The possible wavelength coverage of the TIMED/SEE data reaches 1950 Å, but values above 1900 Å are not provided for most days in our selected range.Because the data do not cover the longest wavelengths of our previously defined F uv , we broke-up this bandpass into F uv,1 ∈ [911, 1900] Å and F uv,2 ∈ [1900, 2593] Å. Next, we ignored days where any bin in F x , F euv,1 , F euv,2 , or F uv,1 did not have a reported flux.This strict quality cut eliminated about 1% of days in the original range.
For each spectrum in this set of days, we summed the flux in F x , F euv,1 , F euv,2 , and F uv,1 .After scaling to what would be detected at the Sun's surface (Chadney et al. 2015), we determined power-law relationships: where F bp is the flux in a given UV bandpass and the best-fit α and γ are constants determined for each fit.
Results for each bandpass are shown in Table 5.As an example, we show the best-fit power-law between F x and F euv,2 in the inset of Figure 7. Importantly, the correlation between bandpassess shows that the Sun's EUV changes in concert with X-ray output across the activity cycle; we should expect X-ray variability to reliably trace EUV variability for exoplanet hosts.
To estimate the spectrum of WASP-69 in the epoch of the Swift data and Night 2 metastable HeI data, we started from the HD 85512 spectrum published by the MUSCLES survey (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016).This star is an active K-type whose observed spectrum was first used to proxy WASP-69's XUV output by Vissapragada et al. (2020).
We scaled the entire HD 85512 spectrum to what would be observed at the surface of the star, then scaled fluxes in each of the aforementioned ranges to be consistent with the Swift F x upper bound from the X-ray observations and the power-law relationships in Table 5.First, we multiplied all fluxes with wavelengths in the F x range such that the total flux in that band summed-up to the upper confidence bound from Swift.For F euv,1 , F euv,2 , and F uv,1 , we multiplied fluxes within each range by values such that the total fluxes in those bandpasses summed to the respective power-law estimates with the parameters in Table 5.We left F uv,2 unchanged from the original HD 85512; the longest wavelength fluxes in TIMED/SEE barely change as F x varies, so we would not expect that changes in F x affect mid-UV fluxes.
This procedure derived an input spectrum for our atmospheric modeling routine, which we display in Figure 7.Although spectral shapes are inconsistent between F x and the UV ranges, the shapes within each range are consistent with an observed star of similar mass and activity level.More importantly, this input spectrum for atmospheric modeling has fluxes in each relevant bandpass that are consistent with contemporaneous X-ray observations and scaling-relations derived from high-fidelity solar data.

Atmospheric Variability Modeling
Since contemporaneous X-ray data only exists for the recent epoch of metastable HeI absorption data, we constructed an atmospheric model to estimate whether the XUV change observed between 2016 and 2023 by XMM-Newton and Swift could explain the metastable HeI variability.We focused on comparing the GIANO Night 1 and WIRC Night 2 datasets because the magnitude of those two best-fit WIRC-like transit depths differed the most among all the metastable HeI absorption data on WASP-69b's thermosphere; if the observed stellar XUV variability could explain that observed atmospheric variability, then all metastable HeI variability could be attributed to changes across WASP-69's activity cycle.While the CARMENES Night 1 transit was collected closest in time to the XMM-Newton visit, the XMM-Newton observation was too time-separated from any metastable HeI data to meet our definition of contemporaneous.Thus, we used the observation with the largest variability as a basis for comparison.
We used p-winds, which implements the isothermal Parker wind outflow formulation by Oklopčić & Hirata (2018) and Lampón et al. (2020).This 1D model describes an outflow for an exoplanet at the substellar point given the planetary parameters, a mass-loss rate from the substellar point Ṁsub , the thermospheric temperature T 0 , the H/He fraction, and an incident irradiating spectrum.To convert to a planet-averaged mass-loss rate Ṁ , we divide the value inputted to p-winds by 4 (Vissapragada et al. 2022a).
For a reference planet-to-star radii ratio when calculating excess depths, we used the aforementioned transit depth from HST (Tsiaras et al. 2018) that was used in Section 3.3.We assumed an atmosphere composed of 90% H and 10% He by number, a common modeling setup for giant planets (i.e.Oklopčić & Hirata 2018;Mansfield et al. 2018), and planetary parameters from Anderson et al. (2014) except for the radius R p = 1.014R J from Tsiaras et al. (2018).In addition, we adjusted WASP-69b's gravity field in the momentum equation by the tidal effects elucidated by Murray-Clay et al. (2009) and implemented in p-winds by Vissapragada et al. (2022a).Given those inputs, we found the outflow velocity and density profiles with p-winds.
We used the radiative transfer functionality in p-winds to obtain the H and He ionization fractions and level populations, assuming an incident spectrum scaled to a/R ⋆ = 12.21 (from Table 2) from the one shown on Figure 7.After calculating the wavelengthresolved in-transit metastable HeI triplet, we Doppler shifted the model spectrum by the planet's radial velocity versus Palomar Observatory on Night 2. Finally, we convolved the simulated metastable HeI signal with the ultranarrowband filter's bandpass to determine the theoretical WIRC transit depth that would be observed for the given exoplanetary outflow parameterization.This convolution procedure is similar to the one used to create WIRC-like light curves from spectroscopic data in Section 3.3, but we are only finding the WIRC-like transit depth for one spectrum in this case.
With this model setup, we selected a fiducial value of T 0 = 10000 K and found the best-fit Ṁ via the optimize.minimizefunctionality in scipy with the Powell solver that matched the planet-to-star radius ratio from WIRC Night 2 from the joint fit in Section 3.3.For an energy-limited outflow, a given (T 0 , Ṁ ) combination implies an efficiency ϵ via the formulation for the mass-loss rate where M p is the planet mass, G is the gravitational constant, and K corrects the gravity field for the Roche lobe (Erkaev et al. 2007).
To illustrate the effect of changing XUV flux on the observed metastable HeI signal, we found the efficiency ϵ at the fiducial T 0 point with F XUV inferred from the Swift upper bound F x .Then, we held ϵ and T 0 constant while adjusting F x and recalculating the expected WIRC transit depths with the same p-winds procedure.Theoretical modeling has suggested that the constant ϵ approximation is valid for planets like WASP-69b (Caldiroli et al. 2021(Caldiroli et al. , 2022;;Wang & Dai 2021).Notably, the terrestrial thermosphere increases in temperature with increasing solar XUV irradiation (Mlynczak et al. 2018), but self-consistently treating that effect is beyond the scope of p-winds.Specifically, we varied the independent variable F x by an order-of-magnitude from the Swift upper bound, then estimated the EUV flux in each band from Table 5 for the input F x value according to Equation 2 and the procedure in Section 5.1.Next, we obtained adjusted Ṁ values from Equation 3 given the modified high-energy fluxes.With these new input parameters for our outflow model, we found the expected WIRC transit depth with the aforementioned p-winds procedure.Finally, we checked to see if these newly-constructed outflows were energetically selfconsistent since excluded outflows could constrain T 0 .
Figure 8 shows the results from this approach: the theoretical WIRC planet-to-star radius ratio as a function of X-ray flux.We ran multiple models with various assumptions on T 0 to probe the effect of changing the thermospheric temperature.Nearly all trial outflows were found to be energetically feasible, and the only inadmissible ones corresponded to low F x and high T 0 .A ver-sion of Figure 8 that includes all WIRC-like metastable HeI transit depths instead of only the extreme values is located in Appendix C.
We found that the direction of the change in metastable HeI variability -decreasing signal for decreasing stellar X-ray flux -is consistent with the hypothesis that the observed changes in the planetary thermosphere resulted from years-long changes in stellar XUV output.While Figure 8 implies that the Xray flux during the GIANO Night 1 observations in 2019 must have been at least 2-5x larger than what XMM-Newton observed in 2016, a critical assumption in our procedure was anchoring our p-winds models at the X-ray flux from Swift 95% confidence bound.Had we assumed a lower value for the X-ray flux on Night 2, then the intersection of the XMM-Newton and GI-ANO Night 1 uncertainty ranges could have been closer to overlapping with the p-winds results.Importantly, the demand on WASP-69's high-energy variability implied by Figure 8 is consistent with broad expectations for years-long variability in stellar X-ray luminosities (Woods & Rottman 2005).Thus, changes in WASP-69b's metastable HeI signal are plausibly explained by changes in its host star's XUV output across the activity cycle.Future epochs with contemporaneous high-energy and metastable HeI data can add more points to Figure 8, probing the correlation's strength and better testing the outflow model.
Notably, Guilluy et al. (2024) analyzed the velocityresolved GIANO data from their years-long baseline (2019-2022) and found a 1.9σ significant drop in the signal from GIANO Night 1 to GIANO Night 3.That change in metastable HeI absorption is in the expected direction given a decrease in stellar XUV luminosity.While our analysis in Section 3.3 found that the WIRClike transit depths from GIANO would have appeared consistent, the nature of our convolution procedure decreases the resolution and information content of spectroscopic timeseries for the purpose of treating data consistently.Guilluy et al. (2024) also presented simultaneous Hα measurements for the 2019 and 2020 transit observations that differed.From that trend, the study suggested that the change in the metastable HeI signal between those two transits may have stemmed from WASP-69b occulting stellar regions with different amounts of activity.While that hypothesis can plausibly explain the difference between the 2019 and 2020 results, no Hα data were taken with the 2022 transit.Thus, no similar hypothesis could be proposed to explain the additional drop in metastable HeI signal in the 2022 GIANO data.Our results, which found a tentative decrease in the metastable HeI signal from that final GIANO epoch to WIRC Night 2 in August 2023, suggest an alternative hypothesis: that the overall drop in the metastable HeI signal stems from years-long changes in stellar XUV flux across an activity cycle.In this scenario, WASP-69b's thermosphere would have changed in response to the changing highenergy radiation balance.Only by adding more epochs with contemporaneous multi-wavelength observationsmetastable HeI as well as XUV, Hα, and other stellar activity indicators -to the dataset can these hypotheses be definitively tested.

DISCUSSION & CONCLUSIONS
Along with archival data, we have analyzed new metastable HeI absorption data from a transit of WASP-69b observed by WIRC and contemporaneous stellar Xray data from Swift.The recent metastable HeI transit depth differed by approximately 4σ versus data from the same instrument four years earlier.High-resolution spectroscopy of the metastable HeI triplet in WASP-69b also shows variable absorption versus other spectroscopic data and versus the WIRC data.To interrogate the cause of this change in the planetary thermosphere, we gathered X-ray and UV data with Swift within a few weeks of the metastable HeI observation.These data revealed that WASP-69's X-ray flux had dropped by at least a factor of two versus a snapshot in 2016, a change that is consistent with the expected magnitude of variability across stellar activity cycles and in a direction consistent with the change in stellar XUV flux required to explain the changes in metastable HeI absorption.
This study highlights the importance of two emerging frontiers in exoplanet research: multi-epoch and multiwavelength campaigns.Looking ahead, the opportunity exists to extend the timeseries of stellar activity and planetary thermospheric state for the WASP-69 system when TESS is scheduled to observe this target in Sector 81.Contemporaneous high-energy data from the host star and metastable HeI absorption data from the hot Jupiter in this upcoming epoch along with simultaneous broadband photometry would provide an unprecedented means to probe atmospheric outflow variability.Additional opportunity to probe outflow variability will be unlocked by forthcoming purpose-built instrumentation (Jentink et al. 2023), and pilot studies like this one can set expectations for these future campaigns.
More broadly, our results contribute to the growing literature of metastable HeI data from exoplanet thermospheres and place WASP-69's atmosphere within the broader geoscientific discipline of aeronomy.Solar system aeronomy has flourished for decades (Liu et al. 2011), as studies of Earth (Mlynczak et al. 2015), Mars (Lee et al. 2017), and the gas giants (Jackman & Arridge 2011), have traced their variable thermospheric and magnetospheric conditions to changing levels of XUV irradiation over the solar magnetic activity cycle.
In this study, we have established the feasibility of probing the same cause-and-effect process in extrasolar environments.Figure 6 summarizes the data that we have analyzed from the WASP-69 system and plots data from terrestrial aeronomic experiments.Solar activity is proxied as the 60-day rolling mean sunspot count3 .Thermospheric conditions are given by the Thermospheric Climate Index (TCI; Mlynczak et al. 2018) 4 , an estimate of the total power radiated by the thermopshpere and a reliable indicator of this layer's temperature.
It is apparent that the terrestrial thermosphere heatsup during solar maximum when the Sun's XUV output reaches its zenith and that the reverse occurs during solar minimum.While the WASP-69 data's cadence is sparser and the uncertainties are larger, hints of these same trends can be discerned in this star-planet system.The thermospheres of WASP-69b and Earth are in different physical regimes (Owen 2019), but this contrast is instrumental to examining the role of stellar radiation and gravitational potential wells in shaping atmospheric dynamics across various environments.
A deeper understanding of exoplanet aeronomy can reveal analogous processes driven by planets' responses' to their stars, possibly answering longstanding questions on exoplanet magnetospheres (Oklopčić et al. 2020) and the distribution of exoplanet radii (Howard et al. 2010).An important input to theories of photoevaporation is the stellar XUV irradiation (e.g.Owen 2019), so observing changes in planetary thermospheres as stellar XUV flux changes across the activity cycle provides a natural experiment by which to calibrate these models.Thus, studying planetary aeronomy on human timescales can help to illuminate mass-loss processes that may occur on orders-of-magnitude longer timescales.
Besides providing a means to test demographic hypotheses, exoplanets can expand the parameter space over which comparative aeronomy operates and provide critical tests of generalized atmospheric models (Bauer & Lammer 2004).Answering these questions requires combining insight from atmospheric physics and astrophysics with additional data on WASP-69 and other planetary systems.In conclusion, the nascent field of exoplanet aeronomy can illuminate these distant worlds in the context of well-studied processes from geoscience and provide hypothesis tests for widely-invoked mechanisms that may explain distributions of fundamental planetary properties.Here, we examine light curve modeling results for a range of aperture sizes and combinations of possibly relevant covariates.Figure 11 plots the best-fit planet-to-star radii ratios at various aperture sizes, allowing free limb-darkening coefficients for each model.To assist with the interpretation, Figure 12 shows the BIC for each fit relative to the BIC at each aperture size where no covariates other than time are considered (∆BIC).To facilitate the BIC comparison, fixed limb-darkening coefficients were used in Figure 12 (see Section 3.1).Together, these plots demonstrate changes in the best-fit transit depth and model preference and over a broad range of methodological parameter space.Night 1 results are insensitive to aperture size and the covariates on which the light curve is detrended; all (R p /R ⋆ ) values on Figure 11 are consistent within their 1σ bounds.The globally preferred light curve from the methodology of Section 3.1 was detrended on the water proxy and airmass with an 11 pixel aperture.In contrast, the derived Night 2 transit depths do change for different covariates.This effect is especially apparent for small aperture sizes, but all those light curves had large rms scatter in the residuals and were not globally preferred.Even though covariate choice is less important for the transit depth determination when using large apertures, those models also result in large rms scatter in the residuals.For 1-2 pixel changes in aperture size, keeping a constant covariate set does not change (R p /R * ) beyond the 1σ uncertainties.At the optimized Night 2 aperture size of 12 pixels, either detrending B. FIGURES FROM JOINT FITTING ALL METASTABLE HE I DATA Here, we display the best-fit light curve models and flux residuals (Figure 13) as well as a corner plot (Figure 14) of posterior samples for all planet-to-star radii ratios from the joint fitting procedure in Section 3.3.In Figure 14, we annotate each panel with the sigma difference between the two (R p /R * ) values from the NUTS runs and denote the line where the two values would be equal.All WIRC-like metastable HeI light curves except the SPIRou data (with large uncertainties) are at least 2σ discrepant from the TESS broadband light curve.
It is apparent that the second WIRC transit, taken nearly four years after any other metastable HeI data, is shallower than all other metastable HeI light curves.The best-fit SPIRou (R p /R * ) is closest to the result for WIRC Night 2, but the uncertainty on SPIRou is large enough that this result is effectively uninformative.Aside from the two WIRC transit depths, the posterior distributions for most of the (R p /R * ) pairs are not strongly correlated.This result indicates that changes in the best-fit ephemeris and shape affect WIRC transit depths more than other timeseries.8 that includes all the WIRC-like metastable HeI transits instead of only the extreme values.For each WIRC-like metastable HeI lightcurve without contemporaneous XUV data, we display a horizontal line at the best-fit planet-to-star radii ratio from the joint fit in Section 3.3.Uncertainty bounds on the WIRC-like metastable HeI transit depths are not displayed on the plot but are presented in Table 3.

Figure 1 .
Figure 1.Gaussian FWHM of WASP-69's PSF plotted versus orbital phase during the WIRC science imaging sequence for Night 1 (blue line) and Night 2 (red line).Normalized radial seeing profiles were constructed and fitted with a Gaussian using photutils (Bradley et al. 2023).Although the Night 2 PSF varied more than the Night 1 PSF, both timeseries were stable over the observations.

Figure 2 .
Figure 2. Photometry of WASP-69 from TESS Sector 55.The lightest grey points are the PDCSAP flux that were downloaded from MAST, and the magenta line shows the best-fit spline from Keplerspline.Detrended flux, obtained by dividing the PDCSAP flux by the spline, is shown with two colors.Dark grey dots are detrended points far from the observed transits that are removed by the mask.Black dots are detrended points that are included by the mask.Upward-pointing arrows indicate expected transit midpoints from orbital parameters published by Kokori et al. (2022).

Figure 3 .
Figure 3. Results from joint-fitting WIRC Night 1 (blue), WIRC Night 2 (red), and TESS (grey) data.Panel (a) shows the phase-folded light curves for all three timeseries.Individual points and their uncertainties are shown with small markers, and Palomar/WIRC data binned to ten-minute intervals are shown with larger markers.Best-fit light curve models are shown as solid lines.Panels (b), (c), and (d) show the residuals for the WIRC Night 1, WIRC Night 2, and TESS timeseries, respectively.

Figure 4 .
Figure 4. Allan deviation plots for each detrended WIRC light curve in the joint fitting with both time series and TESS data.The top and bottom panels show results from the detrended Night 1 and Night 2 light curves, respectively.Colored lines with errorbars denote the rms error of the binned residuals for each light curve.The solid black lines denote the expected noise from purely Poisson statistics, and the dotted line shows the photon noise scaled-up to the rms error of the binned residuals.

Figure 5 .
Figure 5. Two-dimensional histogram of the posterior distributions for (Rp/R⋆)W,2 versus (Rp/R⋆)W,1.The dotted line represents equal planet-to-star ratios for the two metastable HeI transits obtained with WIRC.

Figure 6 .
Figure 6.Summary of comparative aeronomy between the terrestrial thermosphere and WASP-69b.Left-Hand Side: Timeseries of solar activity indicator (top panel) and terrestrial thermospheric conditions (bottom panel) from 1960-2020.Numbers on the top panel label the solar cycles by their numbers.Right-Hand Side: Analogous plots for the WASP-69 system.The top panel shows stellar activity proxied by X-ray flux observed at Earth, and the bottom panel shows WASP-69b's thermospheric conditions as traced in metastable HeI from Section 3.3.Multi-epoch, multi-wavelength campaigns are necessary to extend the well-developed field of solar system aeronomy to the exoplanets.

Figure 7 .
Figure 7.The main plot shows the estimated WASP-69 spectrum (black line), which was established from Swift observations and scaling HD 85512's spectrum (grey line) to solar X-ray/UV flux ratios.Each bandpass-of-interest is shown with background shading.The inset shows the TIMED/SEE data from Solar Cycle 24 from the Fx and Feuv,2 bands along with the best-fit power-law.

Figure 8 .
Figure 8. WIRC planet-to-star radius ratio versus WASP-69's X-ray output, normalized by the Swift 95% upper confidence bound.The top x-axis shows the X-ray fluxes from WASP-69 in terms of what would be measured at Earth.The black point represents the constraint from contemporaneous multi-wavelength data, where the WIRC Night 2 transit depth is taken from Section 3.3.The grey horizontal shaded bar and the blue vertical shaded bar show the WIRClike metastable HeI transit depth from GIANO Night 1 in 2019 (best-fit value from Section 3.3) and the XMM-Newton X-ray flux from 2016, respectively.Lines represent theoretical p-winds model results for various T0 values, where the mass-loss rates are calculated with the same efficiency ϵ that is required to match the contemporaneous data from 2023.WASP-69's X-ray flux from 2016, although not contemporaneous with any metastable HeI transmission, changes in a direction and magnitude consistent with the interpretation that the drop in planetary absorption data over time stems from long-term changes in the stellar XUV across the magnetic activity cycle.

Figure 10 .
Figure 10.Analogous plots to Figure 9 for the Night 2 WIRC photometry.

Figure 11 .
Figure 11.Derived planet-to-star radius ratios (Rp/R⋆) versus aperture size from fitting WIRC photometry with various sets of detrending covariates.Error bars indicate 1σ uncertainties.Vertical shading highlights the aperture size that minimized the rms scatter as described in Section 3.1.In addition to the covariates listed in the legend, all light curves are detrended by time.
C. ADDITIONAL FIGURES FROM P-WINDS MODELINGHere, we show another version of Figure8with all the metastable HeI data that was used in Section 3.3 for comparison to the results of the p-winds modeling from Section 5.2.Only the WIRC Night 2 observations in August 2023 were accompanied by contemporaneous X-ray measurements of WASP-69's high-energy luminosity; this epoch is the only one that can be represented by a point on this plot.

Figure 13 .
Figure 13.Panel (a) shows the best-fit, phase-folded light curve models from the procedure in Section 3.3.Panels (b), (c), (d), (e), (f), (g), (h), (i), and (j) show the residuals of the metastable HeI light curve models versus the best-fit light curve from panel (a).These residual plots correspond to WIRC Night 1, WIRC Night 2, CARMENES Night 1, CARMENES Night 2, SPIRou, Keck, GIANO Night 1, GIANO Night 2, and GIANO Night 3, respectively; the legend in panel (a) also corresponds to the color-coding in panels (b)-(j).This figure is the analog of Figure 3 from Section 3.2 for the joint fit from Section 3.3.

Figure 14 .
Figure 14.Corner plot of the planet-to-star radii ratios for all light curves from the joint fit in Section 3.3 that included broadband TESS observations with metastable HeI data from both ultra-narrowband photometry and spectroscopy.Each of the 2D histograms is annotated with the difference (in σ) between the two parameters as derived from the posterior samples themselves.Dotted black lines on each 2D histogram indicates where the two radii are equal.

Figure 15 .
Figure 15.Analogous plot to Figure8that includes all the WIRC-like metastable HeI transits instead of only the extreme values.For each WIRC-like metastable HeI lightcurve without contemporaneous XUV data, we display a horizontal line at the best-fit planet-to-star radii ratio from the joint fit in Section 3.3.Uncertainty bounds on the WIRC-like metastable HeI transit depths are not displayed on the plot but are presented in Table3.

Table 4 .
Summary of Swift ToO observations of WASP-69.All UT dates and times correspond to September 2023.

Table 5 .
Power-law parameters and scaling factors, per the variables in Equation2, derived from solar data and used to estimate flux in XUV bandpasses for WASP-69.Bandpass α [erg cm −2 s −1 ] γ[unitless]