More Evidence for Variable Helium Absorption from HD 189733b

HD 189733b is currently one of the most extensively studied exoplanets in the literature. Its exceptional observational favorability has made it a popular target for atmospheric characterization studies. It is a large (1.1 R_J ) planet orbiting a somewhat small K dwarf (0.78 R_⊙) on a relatively tight orbit (P = 2.2 days), and is located just 20 pc from Earth (Addison et al. 2019; Rosenthal et al. 2021). Since its discovery in 2005 (Bouchy et al. 2005), it has been observed in wavelengths ranging from the X-ray to the mid-infrared, using a wide variety of observational approaches. HD 189733b is currently the only exoplanet with a possible X-ray transit detection (Poppenhaeger et al. 2013), and was the first planet with a measured infrared phase curve (Knutson et al. 2007). It has an extensively characterized transmission and emission spectrum (e.g., Sing et al. 2011; McCullough et al. 2014; Morello et al. 2016), which provide us with some of the most precise constraints on its atmospheric composition (e.g., Zhang et al. 2020; Harrington et al. 2022). This planet has been the focus of a diverse array of theoretical studies seeking to explain its observed properties, including its atmospheric circulation patterns (e.g., Showman et al. 2009) and the nature of its atmospheric aerosols (e.g., Steinrueck et al. 2021).

HD 189733b is also well suited for mass-loss studies, for many of the same reasons that make it favorable for observations of the lower atmosphere. Atmospheric mass loss is an important phenomenon that could reshape exoplanet demographics, potentially turning mini-Neptunes into smaller and denser super-Earths, as well as creating the Fulton gap (Fulton et al. 2017; Fulton & Petigura 2018) and the Neptunian desert (e.g., Szabo & Kiss 2011; Fulton & Petigura 2018). However, there are many theoretical uncertainties in mass-loss models that make exoplanet atmospheric evolution difficult to model, including the effects of magnetic fields (Adams 2011) and of interaction with the stellar wind. The 10833 Å line has long been proposed as a possible probe of exoplanet atmospheres (Seager & Sasselov 2000; Oklopčić & Hirata 2018), but observations have only recently succeeded in detecting exoplanetary absorption in this line (Nortmann et al. 2018; Spake et al. 2018; Allart et al. 2019). It has been suggested that K-type stars have the most suitable stellar spectrum for populating the metastable helium state (Oklopčić 2019; Poppenhaeger 2022). HD 189733b, aside from being large, nearby, and highly irradiated, orbits a K-type star.

In 2016–2017, Salz et al. (2018) observed three transits of HD 189733b with CARMENES and detected clear He i absorption in all three, with consistent depths. Guilluy et al. (2020) observed five transits with GIANO and obtained inconsistent depths indicating that HD 189733b's outflow may be varying in time. In this paper, we present a new transit observation with Keck/NIRSPEC. Section 2 describes the data acquisition, Section 3 the data analysis, Section 4 compares our results to prior work, Section 5 uses a 1D hydrodynamic model to explore how much the varying high-energy stellar flux should change the signal, and Section 6 summarizes our conclusions.

We observed a transit of HD 189733b with Keck/NIRSPEC on 2020 July 14, from 08:28 UTC to 14:41 UTC, using the NIRSPEC-1 filter and the 0.288 × 12'' slit. These observations cover a wavelength range of 0.947–1.121 μm at a resolution of R = 37,500. Over the course of these 6 hr, we took 252 exposures in a ABBA nod pattern, with each exposure taking slightly over a minute and each ABBA pattern taking roughly 5 minutes. These data cover 1.0 hr of pretransit baseline, the 1.8 hr transit, and 3.4 hr of posttransit baseline. The airmass decreased from 1.2 at the beginning of the night to 1.0 around 0.5 hr after midtransit, before rising again to 1.66 at the end of the observations. The resulting signal-to-noise ration (S/N) per spectral pixel in the continuum surrounding the He i line ranged from 150–280, with a S/N peak 2 hr after transit, and an unfortunate S/N trough coinciding with midtransit.

Two of the raw spectra have multiple traces due to telescope nodding errors: one during egress, and one an hour after the end of egress. We exclude these spectra and their nod companions from the analysis, for a total of four excluded reduced spectra.

The NIRSPEC observation was intended to be simultaneous with a Lyα observation by the Hubble Space Telescope's COS spectrograph (GO 15710). Unfortunately, HST failed to acquire the guide star, so we did not collect any Lyα data.

To analyze these observations, we use the methods outlined in Zhang et al. (2022a). We make a master dark by stacking five darks with an exposure time of 1.5 s and five coadds. We make a master flat by subtracting the master dark from each flat and stacking the 29 flats. We create four difference images out of each A₁ B₁ B₂ A₂ nod group: A₁ − B₁, B₁ − A₁, B₂ − A₂, A₂ − B₂. This subtracts off bias, dark current and skyglow. We then use a custom variant of optimal extraction to extract the 1D spectrum, while masking both bad pixels identified in the dark and flat processing stages, and cosmic rays identified by the optimal extraction algorithm itself. Afterwards, we generate a template spectrum by combining a PHOENIX model (Husser et al. 2013) with a telluric model, accounting for the radial velocity of the star relative to Earth. We use this template spectrum to obtain a wavelength solution for order 70 (containing the helium line) in each individual exposure and then run molecfit to correct for tellurics. This last step is not strictly necessary for these data because–due to a fortuitous combination of the star's radial velocity and Earth's orbital velocity–there are no telluric lines (with the possible exception of microtellurics) that overlap with the helium line.

Having obtained wavelength calibrated and telluric corrected spectra, we interpolate all spectra onto a common wavelength grid with a uniform logarithmic spacing of λ/110,000 and a range of 10,810–10,850 Å. The resolution of 110,000 is chosen so that the spacing of the grid is equal to that of the spectral pixels at 10833 Å. We remove fringing by applying a notch filter twice, using the same parameters as in Zhang et al. (2021). We divide each spectrum by the continuum, take the logarithm of the entire spectral grid, and subtract the mean of every row and column from that row and column. The end result is a N_obs × N_wav grid of numbers representing the relative deviation of a pixel from the mean for that row and column.

After obtaining this residuals grid, for every column (wavelength), we subtract the mean of the out-of-transit part of the residuals image for that column; we then invert the residuals image. The resulting residuals image now shows the excess absorption relative to the out-of-transit baseline. However, there are still continuum variations that contribute structure in this image. We correct for these variations by masking out the strong lines (including the helium line), fitting a third order polynomial to each row (epoch) with respect to wavelength, and subtracting off the polynomial. The result is shown in Figure 1. Finally, we shift all of the in-transit residuals spectra to the planetary rest frame and stack them to obtain the 1D excess absorption spectrum shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare our results to the eight literature values for HD 189733b's helium absorption signal. Salz et al. (2018) observed three transits with the fiber-fed Echelle spectrograph CARMENES (R ∼ 80,400), situated on the 3.5 m Calar Alto telescope. Guilluy et al. (2020) obtained five transits with the slit-fed Echelle spectrograph GIANO-B (R ∼ 50,000) on the 3.6 m Telescopio Nazionale Galileo telescope. The CARMENES observations have a S/N of 160–240 per spectral pixel in the continuum surrounding the helium line (Salz et al. 2018), but their exposure time is closer to five minutes than one minute, and their wavelength dispersion is 70% higher. Binning our observations to their time resolution and binning their observations to our spectral resolution, we find that our S/N is 1.6–2.0× times their single-night S/Ns. The same comparison, performed for the GIANO data, reveals that our S/N is 3.6–4.4× their single-night ratios. The GIANO comparison is not perfect, as Guilluy et al. (2020) report their ratios over the whole order containing the helium line, not the immediate region surrounding the line.

Guilluy et al. (2020) report the depths for each individual transit, which they calculate by averaging from the end of ingress to the beginning of egress for a wavelength range of −20 to 20 km s⁻¹ in planet-centric velocity space. We use their definition to calculate the analogous average transit depths for the CARMENES transits (using values taken from their Figure A.3) and for our NIRSPEC transit. We did not take into account the differing resolutions of the instruments when defining the equivalent pixel range for our wavelength bounds. We quantified our sensitivity to the instrumental resolution by convolving the CARMENES spectrum down to the NIRSPEC resolution of 37,500 and the NIRSPEC pixel sampling before calculating the average transit depth, and found that it did not change the result by more than 0.01% for any of the three transits.

We compare the resulting transit depths for each observation in Figure 3. This figure shows that the consistency of the three CARMENES transits over a year appears to have been a coincidence, and is not representative of the broader sample. The GIANO transits show widely varying depths, while our NIRSPEC depth is lower than all but one of the previous transit observations. Combining the depths for every instrument with a weighted average, we find an average transit depth of 0.617% ± 0.017% for CARMENES, 0.508% ± 0.015% for GIANO, and 0.420% ± 0.013% for NIRSPEC. The error bars are the standard deviation of the mean, not the sample. The CARMENES average differs from the NIRSPEC result by more than 9σ, and from the GIANO result by 4.4σ.

4.1. Effect of Stellar Variability on Observed He Transit Depths

We consider two ways in which stellar variability can affect the observed transit depths. First, if the stellar disk is inhomogenous and the planet does not transit a perfectly representative chord, the transit depth will be biased. Second, the stellar He line may change during the observations, making the transit appear deeper or shallower than it really is.

Guilluy et al. (2020) used the five GIANO transit observations, which also included simultaneous Hα measurements, to extensively explore the first type of variability. Hα is an excellent tracer of stellar activity, making it a useful proxy for the spot coverage fraction during each visit. In active regions on the star, Hα is seen in emission while the stellar He i line becomes deeper. Guilluy et al. (2020) detected changes in the Hα line during the two transit observations with the highest and lowest He i depths, suggesting that the planet was transiting an inhomogenous star. The remaining three transits had consistent He i depths (∼0.5%) and no detectable variations in the stellar Hα line. The deep He i transit could be explained if the planet transited an especially quiet region of the star (thus preferentially leaving the active regions unocculted). However, it is difficult to explain why the other transit that took place during a period of enhanced activity is shallower than expected. Guilluy et al. (2020) argue that this observation is unlikely to be caused by a starspot occultation, as they do not observe a concurrent weakening in the Si 10830 Å line, which is sensitive to changes in temperature. They tentatively suggest that the planet might have occulted a filament instead. Our new transit depth is also shallower than expected, and we also do not see any evidence for a concurrent change in the Si 10830 Å line during the transit. However, unlike Guilluy et al. (2020), we do not see a midtransit spike in the He band-integrated light curve, nor any other abnormalities that would indicate the occultation of a filament or starspot.

We can quantify the effects of stellar active region contamination on the helium transmission spectrum using the techniques outlined in Cauley et al. (2018). In order to reproduce the strength and shape of the stellar He i 10833 Å triplet we found that an active region filling factor of ≈65% and an optical depth at line center of τ ≈ 0.7 are required. We note that this is the facular or plage filling factor and not a spot filling factor. The large filling factor is in line with previous estimates for helium absorption on active K-dwarfs (Andretta et al. 2017) and is similar to the value of 75% estimated by Guilluy et al. (2020). After estimating the filling factor, we simulated the transit of HD 189733 b across the active stellar surface and tuned the atmospheric parameters to match the observed transmission spectrum. We then tested various filling factors and geometries (e.g., uniform distribution, active latitudes) to see how the transmission spectrum changes while the atmospheric parameters are held constant. We find that even drastic changes in the filling factor (≈30%–40%) are not enough to reproduce the difference between the atmospheric absorption observed by Salz et al. (2018) and our NIRSPEC transmission spectrum.

The strength of the stellar helium line does not change much from epoch to epoch: Guilluy et al. (2020) found that the line core flux varied by only ≈4% across a 1.5 yr baseline. The NIRSPEC spectrum displays a similar line depth (≈0.7) compared to Guilluy et al. (2020) and Salz et al. (2018), suggesting that the surface activity level, at least as measured by absorbing metastable helium regions, is fairly constant. This implies that the filling factor cannot change dramatically from epoch to epoch, otherwise we would have observed large variations in the stellar line strength (of order the percent change in filling factor). The relative stability of the star's filling factor and the large size of the changes in the filling factor required to account for the epoch to epoch variability in the observed helium absorption signal make it unlikely that this variability is due to inhomogeneities on the stellar surface.

On the other hand, variability in the stellar He line is a plausible contributing factor to the differing transit depths. Figure 4 shows the light curve of the helium line, centered slightly blueward of the main peak. The planetary absorption is evident, but so is substantial variability after transit, including a sudden drop 2 hr after midtransit and a slow 0.6% rise in the 2 hr following the drop. This variability can also be seen from a close examination of Figure 1, but the band-integrated light curve makes it clearer. If similar stellar variability occurred during the transit or in the immediate out-of-transit baseline, the planetary excess absorption could appear ∼0.2% larger or smaller than it actually is.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Potential Evidence for Variability in the Planetary Outflow

It is plausible to think that the observed variability might instead be due to changes in the planet's atmospheric outflow, which might be caused by variations in the star's high-energy radiation or stellar wind environments. In this case, we might also expect to see visit-to-visit variations in the shape of the wavelength-dependent planetary absorption signal. Figure 5 compares the excess absorption spectrum obtained by the three instruments, which tells the same story as the average depths: NIRSPEC sees a similar signal as GIANO-B (all transits combined), which in turn sees a weaker signal than CARMENES.

Zoom In Zoom Out Reset image size

One interesting difference between the GIANO and NIRSPEC excess absorption spectra is that despite having a very similar primary peak, the height of the secondary peak at 10832 Å is lower for GIANO by ∼0.1%. Adopting the formal errors on the two spectra (0.07% for GIANO and 0.05% for NIRSPEC) and assuming statistical independence, we calculate that the difference is significant at the ∼3.5σ level. However, the significance of this feature is likely overestimated, as noise in high-resolution spectra is often correlated across adjacent wavelengths. If the difference is real, we can use it to constrain the optical thickness of the outflow. A completely optically thin outflow would have a primary-to-secondary peak ratio of 8:1 (derived from the g_i f_ik of the three lines), whereas a completely optically thick outflow would have a peak ratio of 1:1, as the gas would absorb all light in both peaks. The higher peak ratio in the GIANO data could indicate that the absorption is coming from a more optically thick region than in our NIRSPEC observations. However, the fact that we see significant absorption extending between the two peaks makes it difficult to determine how much of the absorption at 10832 Å is due to the line there and how much is due to absorption from strongly blueshifted gas absorbing at a rest wavelength of 10833.3 Å. Given the right kinematic structure–namely a long tail of gas being accelerated away from the star by the stellar wind–it is conceivable that even the CARMENES observations, for which the apparent peak ratio is 2.8 (Salz et al. 2018), might be consistent with an optically thin outflow.

We can use atmosphere models to quantify the effect of stellar variability on the helium absorption signal. Here we focus on the effect of variations in the stellar XUV flux, which can be captured using relatively simple 1D hydrodynamic models. The code we use is The PLUTO-CLOUDY Interface (TPCI; Salz et al. 2015a), a combination of two sophisticated and widely used codes: the hydrodynamic solver PLUTO (Mignone et al. 2007), and the plasma simulation and spectral synthesis code CLOUDY (Ferland et al. 2013). TPCI has been extensively used to study photoevaporation (e.g., Salz et al. 2015b, 2016; Kasper et al. 2020; Zhang et al. 2022a).

We obtain the nominal stellar spectrum in two ways. In the first method, we use the XUV spectrum derived by Lampón et al. (2021b; their Figure 2, data provided by M. Lamp'on) for the wavelength range of 5–1450 Å. We chose a cutoff of 1450 Å because their spectrum redward of that wavelength is not based on observations. In the second method, the X-ray spectrum is taken from the thermal plasma model of Poppenhaeger et al. (2013), the Lyα flux at 1 au is taken to be 11.8 erg s⁻¹ cm⁻² (France et al. 2013), and the EUV spectrum is derived from the Lyα flux using the scaling relations of Linsky et al. (2014). For both methods, fluxes at longer wavelengths are taken from a PHOENIX model (Husser et al. 2013; T_eff = 5000, log(g) = 4.5, solar metallicity). Of particular importance is the spectrum at 1230–2588 Å, which we call midultraviolet, because this radiation ionizes metastable helium, but cannot contribute to producing metastable helium because it cannot ionize ground state hydrogen. Fortunately, we have observational constraints for the flux at these wavelengths from XMM-Newton's Optical Monitor, a photometer which operates simultaneously with the X-ray instruments. We reviewed all archival observations taken by the OM in the UVW2 (212 nm, width 50 nm) and UVM2 (231 nm, width 48 nm) filters, and found that they do not differ by more than a few percent. (For a summary of all XMM-Newton observations of HD 189733, see Pillitteri et al. 2022.) To roughly match the observations, we boosted the model spectrum between 1230–2588 Å by 30% and used it to predict the count rates that OM should see. The prediction was only 6% too low for UVW2 and 3% too high for UVM2, indicating our model MUV flux is accurate. The final spectra constructed using the two methods is shown in Figure 6, while Table 1 shows the band-integrated fluxes.

Zoom In Zoom Out Reset image size

Comparing the two methods, the X-ray flux is in excellent agreement: Lampón et al. (2021b) reports an X-ray luminosity that corresponds to a flux of 7.8 erg s⁻¹ cm⁻² at 1 au, while we obtain 7.6 erg s⁻¹ cm⁻². However, the EUV (100–912 Å) flux is dramatically discrepant. Lampón et al. (2021b) imply a flux of 57 erg s⁻¹ cm⁻² at 1 au, while we obtain 6.5 erg s⁻¹ cm⁻², a difference of 9×. Other authors have obtained 11 (Sanz-Forcada et al. 2011), 22 (Bourrier et al. 2020), and 36 erg s⁻¹ cm⁻² (Poppenhaeger et al. 2013). This large uncertainty in the EUV flux is ultimately because the stellar EUV is only measurable from space and no space telescopes currently have EUV capabilities (France et al. 2022). The uncertainty has been noted by many previous publications. For example, Oklopčić (2019) compared the results of two different methods of EUV reconstruction for the same star that differed by an order of magnitude, which produced helium excess absorption depths that differed by a factor of 3. France et al. (2022) compared four different methods of reconstructing the EUV spectrum for Proxima Centauri and found that they were discrepant by 3–100×, depending on the wavelength.

Since the EUV flux is uncertain to an order of magnitude, we analyze the behavior of the escaping atmosphere over a wide range of EUV fluxes. For the first stellar spectrum we constructed, which has low-EUV flux compared to literature values, we set up four 1D hydrodynamic simulations: at 1×, 1.3×, 2×, and 3× the nominal stellar XUV flux as seen by the planetary dayside (which is half the flux seen by the substellar point). For the second spectrum, which has high EUV flux compared to other literature values, we set up five 1D hydrodynamic simulations: at 0.5×, 1×, 1.3×, 2×, and 3× the nominal stellar XUV flux. Following Zhang et al. (2022a, 2022b), we ran the simulation for 100 time units with advection off, where the time unit is calculated as the planetary radius divided by 10 km s⁻¹ (roughly the sound speed). Unlike in this previous study, we did not run the models for another 100 time units with advection on, due to numerical problems. Also unlike in the previous study, we irradiate the atmosphere with the average flux experienced by the planetary dayside, and not the flux experienced by the substellar point. The former is half of the latter, and is more representative of the outflow as a whole.

Figure 7 shows the results of the TPCI simulations as a function of EUV flux. The simulation results are similar to those of other 1D models. In agreement with Lampón et al. (2021a), we find that the the outflow is hot (>∼10,000 K) and in the recombination regime, which is characterized by a sharp transition from neutral to ionized hydrogen. For the 1/2 nominal XUV model using the first stellar spectrum, which is the model with median EUV flux (see Figure 7), the outflow is 95% neutral at 1.12 R_p and drops to 5% neutral by 1.32 R_p . In the recombination regime, most of the incident XUV is radiated away via recombination, leading to low mass-loss efficiency. In fact, the mass-loss rate we obtain for the aforementioned model, 5 × 10⁹ g s⁻¹, is substantially lower than even the 10¹¹ g s⁻¹ found by Lampón et al. (2021a, 2021b), and similar to the 10¹⁰ g s⁻¹ found by Caldiroli et al. (2022) for similar XUV flux. This corresponds to a mass-loss efficiency of 1.3%, similar to the 2%–3% found by Caldiroli et al. (2022).

Zoom In Zoom Out Reset image size

The absorption depth is higher than observed even in the simulation with the lowest EUV. Increasing the EUV flux 34% above this lowest point increased the helium absorption by 60%, but the depth increases roughly with the square root of EUV flux thereafter, so that further 34% increases in flux only produce 15% increases in depth.

HD 189733 has been extensively observed at XUV wavelengths, providing us with good empirical constraints on the magnitude of the star's variability. Pillitteri et al. (2022) analyzed 25 XMM-Newton observations spanning a total of 8 yr, the last of which was taken before the first helium observation. They found that the 25%–75% range of quiescent flux varies by 23%, and that 44% of the total observing time was occupied by flares, which typically lasted several kiloseconds and increased the XUV flux by a few tens of percent (see their Figure B.2). Using the data behind their Figure B.3 (I. P. Pillitteri 2022, private communication), we calculate that the typical variability in the X-ray flux, as defined by the gap between the 84th and 16th percentiles, is 60%. EUV flux increases with X-ray flux as ${F}_{\mathrm{EUV}}\propto {F}_{{\rm{X}}}^{0.86}$ (Sanz-Forcada et al. 2011), so a 60% increase in X-rays corresponds to a 50% increase in EUV, and a 23% increase in X-rays corresponds to a 19% increase in EUV. These numbers suggest that the observational variability plotted in Figure 3 can plausibly be due to variability in stellar EUV output if the EUV flux is low (∼6 erg s⁻¹ cm⁻² at 1 au), but is less likely if the EUV flux is high (>∼10).

It is worthwhile to briefly consider other sources of variability that we could not explore with our 1D model. The effect of differing stellar wind conditions on the outflow has been extensively studied with 3D simulations (e.g., McCann et al. 2019; Vidotto & Cleary 2020), but generally without modeling the helium line. MacLeod & Oklopčić (2022) does model the helium line, and find that given the same planetary mass-loss rate, the equivalent width of helium absorption changes by 10% when the stellar wind is increased 10× from weak to moderate, but changes far more drastically when the wind is strong enough to "break through" the sonic surface and confine the outflow. Zhang et al. (2022a) study the effect of different stellar wind conditions on the helium absorption from TOI 560b, a young mini-Neptune, by alternately halving the stellar wind density and velocity. They find changes in peak absorption of up to 25%, but smaller changes in equivalent width of up to 12%. These studies show that changes in stellar wind conditions can plausibly explain a portion of the variability we see in HD 189733b.

Another source of variability is shear instability, which exists even with constant irradiation and wind from the star. For the inflated sub-Saturn WASP-107b, Wang & Dai (2021a) uses a 3D hydrodynamic model to predict ∼10% fluctuations in helium absorption over hour-long timescales. These fluctuations, however, change the planetary mass-loss rate by less than 1%. Shear instability, like the impact of the stellar wind, is not possible to model in 1D.

One final source of variability that has been modeled in the literature is flares. Flares increase the star's XUV output, and a sufficiently long-lasting flare would increase helium absorption, as predicted by our models. However, depending on the spectrum of the flare, the immediate effect of a flare might be to ionize metastable helium and reduce helium absorption. After a dynamical timescale, when the surge of photoevaporative mass loss generated by the flare has reached higher altitudes, helium absorption should increases above baseline (Wang & Dai 2021b). For WASP-69b, the planet simulated by Wang & Dai (2021b), this takes about an hour; helium absorption then stays elevated for 2 days after the end of the flare. Thus, even though flares can both decrease and increase helium absorption, the latter is far more likely to be observed than the former. To explore the effect of flares that raise the XUV flux by tens of percent on a high-gravity planet like HD 189733b–rather than a 10x flare on a low gravity planet like WASP-69b, as Wang & Dai (2021b) did–new 3D simulations are desirable.

Metastable helium observations provide us with a sensitive probe of planetary outflow characteristics. In this paper, we present a new helium transit observation of HD 189733b, a hot Jupiter that has been observed over a wide range of wavelengths in the fifteen years since its discovery. This planet will be also observed by 5 James Webb Space Telescope programs in Cycle 1 alone; these programs will doubtlessly further expand our knowledge of its atmospheric properties.

When combined with the set of eight published helium transit observations of HD 189733b, our new observations add to the growing body of evidence suggesting that this planet's outflow properties may vary in time. HD 189733 is a relatively active K dwarf and it has been extensively observed at high energies by XMM-Newton and Chandra. Our simulations suggest that stellar XUV variability, by itself, can plausibly explain the observed variations in the planet's helium transit depth. However, this does not mean that other sources of variability we do not model, such as stellar wind variations, shear instability, and short flares, cannot also explain the variability.

Understanding the high-energy environment and hydrostatic atmosphere of the planet is crucial to understanding mass loss, and the quantity of atmospheric observations available for HD 189733b make it one of the most promising planets for understanding mass-loss physics in detail. In order to better understand the amplitude, timescale, and likely cause of this planet's observed variability, it is important to continue monitoring of the planet in the helium line, preferably with simultaneous Hα and/or X-ray measurements. This will in turn inform how much we should rely on single-epoch observations of helium absorption to quantify the mass-loss rates and population-level properties of other close-in planets.

The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. A.O. gratefully acknowledges support from the Dutch Research Council NWO Veni grant.

Software: numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), TPCI (Salz et al. 2015a).