Axisymmetric High Spot Coverage on Exoplanet Host HD 189733 A

Transmission spectroscopy is a powerful tool to study exoplanet atmospheres, which can be affected by the ability of stellar photospheric heterogeneity to mimic or mask exoplanetary spectral signatures. The canonical HD 189733 system provides a textbook example of this spectroscopic discrepancy with features that have been variously interpreted as signatures of scattering by haze in the planetary atmosphere or unocculted spots on the stellar disk. Here, we leverage three archival data sets from the Hubble Space Telescope to directly infer the covering fraction of HD 189733 A and explore the evidence for photospheric heterogeneity in the out-of-transit spectra. We model the stellar spectrum using one to three spectral components in a nested-sampling framework, finding that the two-component model (photosphere and spot) is preferred for all data sets. We find photospheric and spot temperatures of 5295±4143 and 3222±116100 K, respectively, which are consistent across data sets. The spot covering fraction is large and varies between 38% ± 4% and 47% ± 3%. Combined with time-domain insights from Transiting Exoplanet Survey Satellite data revealing HD 189733 A's 1.4% peak-to-peak variability, our findings imply that the majority of the spots must be distributed axisymmetrically, e.g., in a densely filled latitudinal band or at the poles. More work with complementary data sets is necessary to investigate those possible arrangements.


Introduction
Transiting exoplanet discoveries have paved the way for astronomers to study the atmospheres of distant worlds.Light passing through an exoplanet's upper atmosphere allows for constraints on its physical and chemical makeup, owing to wavelength-dependent absorption from different opacity sources (Seager & Sasselov 2000;Brown 2001;Benneke & Seager 2013;de Wit & Seager 2013).However, the confident attribution of spectral features to planetary atmospheres is only possible for a uniform stellar surface.Stellar photospheric heterogeneities, such as bright (faculae) or dark (starspot) features, can mimic the spectral patterns of atmospheric absorption and thus hinder our ability to probe exoplanet characteristics (Pont et al. 2008;Bean et al. 2010;Berta et al. 2011;Sing et al. 2011;Jordán et al. 2013;Rackham et al. 2018Rackham et al. , 2019Rackham et al. , 2023)).
The HD 189733 system is a binary pair of K2V (Gray et al. 2003) and M4V (Bakos et al. 2006) stars of significant interest to the exoplanet community, as the primary K-dwarf star, HD 189733 A, hosts the well-studied, canonical hot Jupiter HD 189733 b (Bouchy et al. 2005;Knutson et al. 2007Knutson et al. , 2008;;Pont et al. 2008;Redfield et al. 2008).Due to its short (∼2.2 days) orbital period (Bonomo et al. 2017) and its large size relative to the host star (R p /R s = 0.16; Boyajian et al. 2015), HD 189733 b is a prime target for atmospheric studies.Yet, due to HD 189733 Aʼs high photospheric surface activity (Aigrain et al. 2012), surface heterogeneities, such as dark spots, can imprint stellar spectral features on the observed transmission spectrum of the exoplanet, a phenomenon known as the "transit light source (TLS) effect" (Rackham et al. 2018(Rackham et al. , 2019)).
Ten years ago, Pont et al. (2013) studied the transmission spectrum with the Hubble Space Telescope (HST) using five transits with the Space Telescope Imaging Spectrograph (STIS), one transit with the Advanced Camera for Surveys (ACS), and two transits with the Wide Field Camera 3 (WFC3) instruments.They found a strongly sloped optical spectrum, which they interpreted as the result of haze.Shortly after, McCullough et al. (2014) revisited the spectrum and offered a differing interpretation of its optical slope.Adopting a spot temperature of 4250 K, as measured from a spot-crossing event (Sing et al. 2011), they note that it can be wholly due to unocculted spots covering 4.3% of the stellar disk.The difference in interpretation hinges on the total spot coverage of HD 189733 A: While Pont et al. (2013) correct for the level of spot coverage that is evident from the rotational variability of the star, McCullough et al. (2014) find that a slightly higher spot coverage, not discernible from rotational variability due to its axisymmetric distribution, can fully explain the optical slope.
Recently, new work has shown that the baseline out-oftransit HST/WFC3 stellar spectra themselves provide a powerful probe of the photospheric heterogeneity of the star.Zhang et al. (2018) used out-of-transit spectral modeling of TRAPPIST-1 to show that its photosphere contained both faculae and starspots, which provide a significant additional noise source in the interpretations of the planets' transmission spectra (de Wit et al. 2016(de Wit et al. , 2018)).Using new HST/WFC3 observations of TRAPPIST-1, Wakeford et al. (2019) compared models of out-of-transit stellar spectra to observations and emphasized the importance of converting the model data to physical units to avoid unphysical inferences from normalized Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
spectra.Similarly, the study of the planet TRAPPIST-1 h (Garcia et al. 2022) illuminated the utility of incorporating constraints from photometry in interpreting relatively narrow wavelength ranges of spectra in order to avoid seemingly good solutions that violate photometric constraints from other spectral wavelengths.
Here, we use HST/WFC3 spectra of HD 189733 A to study the photospheric heterogeneity of this important exoplanet host.Specifically, we fit out-of-transit spectra from three archival HST/WFC3 observations of HD 189733 A using models of uniform and heterogeneous photospheres to understand the evidence for spottedness in the spectra and constrain the temperatures and filling factors of potential spots.We find the model that best describes the stellar spectrum includes spots with consistent temperatures over time but filling factors that vary between visits.
This paper is organized as follows.Section 2 details the archival data sets we used and our data-reduction process.Section 3 presents our modeling approach and parameter retrieval framework.Section 4 introduces the results of the model fitting and comparison.Section 5 discusses our results in the context of previous observations of HD 189733 A, including those from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015).Finally, Section 6 presents our conclusions in relation to the HD 189733 system and the broader stellar contamination discourse, while also outlining future work to be done.

Hubble Space Telescope Observations
In this study, we utilize three archival data sets from HST/ WFC3 observations made in spatial scanning mode from 2013 June 5, 2013 June 24 (HST program 12881, PI: McCullough), and 2022 June 21 (HST program 16180, PI: Kataria; see Table 1).The spatial scanning mode was pioneered to help with spectral observations of bright stars, like HD 189733 A (McCullough & MacKenty 2012).All data sets were collected using the G141 grism, covering a spectral range from 1.1 to 1.7 μm on the 512 × 512 pixel subarray.We take the publicly available Intermediate MultiAccum images (ima.fits) and separate the images according to the scanning direction,5 resulting in two sets of data per observation date, similar to the initial process of previous studies (Crouzet et al. 2014;McCullough et al. 2014).Nominally, we do not expect differences in spectral features between scan directions, though we split the data sets in this way to investigate any potential systematic differences.

Data Reduction and Processing
We reduced each of the three, five-orbit data sets using the HST data-reduction Python package PACMAN (Zieba & Kreidberg 2022). 6We carried out Stages 00 to 20 of the pipeline to arrive at the extracted wavelength-dependent spectral data.For the observations taken in the forward scan direction on 2013 June 5, the spectra were off the subarray on the first two of the six total samples, so we remove those samples and adjust the total exposure time in the later processing stage to reflect the amount of samples used.For Stage 10, we ensured that the star location on the direct image is found by updating the di_rmin, di_rmax, di_cmin, and di_cmax parameters in the obs_par.pcffile, as done in the example. 7For Stage 20, we set the rmin, rmax, and window parameters in the obs_par.pcffile to a wide value in order to extract the full 2D spectra from the scans, as done in the example. 8In order to preview the scan images and extract the full 2D 1st-order spectra, we utilized the visualization software package SAOImageDS99 (Joye & Mandel 2003).After extracting each scan's spectra, we interpolated all of the spectra to the same wavelength grid, binned the uniformly sampled spectra into 50 μm bins, and finally took the median over all observations per visit/scan direction.The final output of our data-reduction process is a set of six spectra, reflecting the median out-of-transit observation of HD 189733 A for each pair of visits and scan directions.Taking the median of the spectral flux over the five orbit collections per observation date allows us to negate the effects of the transit/eclipse.Additionally, each of the data sets also has a wavelengthdependent spectral error flux value that is Poisson distributed, representing the photon noise error.
As shown previously (Wakeford et al. 2019;Garcia et al. 2022), it is important to convert the data and model to the same physical units, turning the spectral observations (given in electron counts, e − ) into units of spectral flux density F λ at the stellar surface (erg s −1 cm −2 Å −1 ) in order to facilitate comparisons to model spectra.To do so, first we divide the binned median observed electron counts by the WFC3/G141 1st-order sensitivity curve.10Next, we divide the resultant data by the integration time for each observation date, 5.97 s for 2013 June 5 and 2013 June 24 and 23.77 s for 2022 June 21.As noted previously (McCullough et al. 2014), to adjust the exposure time for the removed samples in the forward scan observations on 2013 June 5, we take the exposure time to be

= ´=
. Lastly, we again divide by the size of the wavelength bins, 50 μm.
At this point, our data are in the same units as the model spectra (erg s −1 cm −2 Å −1 ), though in terms of spectral flux density at the detector rather than at the stellar surface.As detailed in Section 3, we fit for the stellar radius and distance in our retrievals using priors from literature estimates in order to perform the final conversion from flux received at the detector to F λ at the stellar surface.Finally, the WFC3/G141 sensitivity quickly drops to zero on the edges of its bandpass, so to avoid large errors in the subsequent modeling, we trim the processed spectra outside the 1.125-1.650μm bands, as done by Garcia et al. (2022).
Figure 1 shows the median spectra divided by the exposure time for each scan direction per visit.From this, we note that the overall flux magnitude recorded during the 2013 June 5 visit is ≈11% (Reverse) and ≈30% (Forward) lower than the median spectral fluxes from the other two visits due to varying amounts of the spectral scan falling off the detector in each scan direction.In addition, there is ≈20% less flux in the forward scan data relative to the reverse scan direction.To understand the curvature of the observed spectra relative to each observation scan data set, we subtract the overall mean flux between the six spectra from each individual median spectrum in Figure 1.In doing so, we see that the spectra are curved with respect to the median, with the higherflux spectra and lower-flux spectra showing different curvatures.The offsets between spectra from forward and reverse scans of the same visits, along with the curvature dependence on overall flux, indicate that the offsets relate to observational systematics and that a model which utilizes a multiplicative scale factor (a scale ) for each spectrum would be appropriate for the parameter retrieval process, which we explore further in Section 3.

Model
In the following, we describe the components of our modeling effort, including the model stellar spectra we used (Section 3.1), the details of our retrieval framework (Section 3.2), our priors (Section 3.3), and the details of our model inference process (Section 3.4).

Model Stellar Spectra
To generate model stellar spectral data for this study, we employ the PHOENIX stellar spectral model grid (Husser et al.  2013). 11For our purposes, the PHOENIX model grid provides stellar surface spectra for temperatures of T ä {3000, 6000} K, surface gravities of g log 0.0, 6.0 10 ( ) { } Î , and metallicities of [Fe/H] ä { −4.0, 1.0} in steps of 100 K, 0.5, and 1.0 dex, respectively.To ensure that our spectra parameters are not significantly stellar model dependent, we perform additional retrievals using the MPS-ATLAS (Witzke et al. 2021;Kostogryz et al. 2023)

Parameter Retrieval Model
For each of the three, five-orbit data sets, we retrieve stellar parameters using models that account for one, two, and three spectral components (referred to as "1comp," "2comp," and "3comp" models hereafter).Comparing these models allows us to determine if there is statistically significant evidence for more than one photospheric contribution to the HST spectra of HD 189733 A, e.g., from the quiescent photosphere and spots.First, the model converts the HST observed spectral flux to the flux output at the surface of the star, based on the sample-drawn stellar radius and distance parameters.This is calculated as where F obs , Surf l is the spectral flux at the stellar surface, F obs , (0.818 ± 0.012 R e ; Boyajian et al. 2015), and a scale is a multiplicative scale factor.Simultaneously, we calculate the model spectral energy distribution (SED) in the form of Equation ( 1), but without using the scale factor a scale , as the absolute flux of the SED is well-constrained for HD 189733 A in the Bessell B and V and Two Micron All Sky Survey J, H, and Ks bands (Stassun et al. 2018). 13We include this SED constraint in order to better anchor the overall magnitude of the model spectra and ensure that the parameters we retrieve do not lead to spectra that disagree with SED constraints.
For the flux conversion for a given model, the scale factor applied to the HST spectra, a scale , is included as a reflection of the discussion of the median spectra in Section 2. We find that the multiplicative discrepancies between the data and nominal model could arise from a systematic error in the sensitivity curve flux calibration or from errors in the reduction process due to idiosyncrasies of the data sets, such as the spectra scanning off the subarray in the 2013 June 5 forward scan direction data set (McCullough et al. 2014).We note that Wakeford et al. (2019) also added a model offset term in their analysis, but interpreted it in terms of an astrophysical origin (a systematic offset to the stellar radius).
For this retrieval process, we aim to maximize the loglikelihood function, which relates the observed data F obs,λ (at the stellar surface) to the model flux F model,λ (with the subscript λ denoting the wavelength dependence).We calculate Equation (2) for both the HST spectral data and the SED separately and then add them.In this work, we take the natural logarithm of the likelihood function, as given by in which the total noise σ λ is the quadrature sum of the photon noise σ photon,λ and an additional "jitter" noise term σ jitter,λ , i.e., .
We calculate σ jitter,λ as f F , 4 jitter, var model, in which f var is a scaling parameter that defines this additional noise in terms of a fraction of the model flux.This noise term procedure follows the methods of Rackham & de Wit (2023), adopted in turn from Foreman-Mackey et al. (2013, 2019). 14he free parameters in all models are the primary photospheric temperature of the star (T 1 ), the stellar radius (R star ), the Earth to HD 189733 A distance (D star ), the multiplicative scale factor (a scale ), and the additional scaled model uncertainty (given as f ln var ).In the 2comp model, additionally, there is the secondary photospheric temperature (T 2 ) and its covering fraction ( f 2 ).In the 3comp model, there are also the secondary and tertiary photospheric temperatures (T 2 , T 3 ) with their respective covering fractions ( f 2 , f 3 ).Following Rackham et al. (2023), we define the covering fractions in terms of the ratios Reworking these equations, we arrive at the formulations for the covering fractions in terms of the ratio parameters: All together, these formulations result in five, seven, and nine free parameters in the 1comp, 2comp, and 3comp models, respectively.

Priors
The full set of priors used in the 1comp, 2comp, and 3comp models is given in Table 2.We place normal priors on D star using the Gaia Data Release 3 (DR3) measurement (19.7758 ± 0.0063 pc; Gaia Collaboration et al. 2023) and on R star using the interferometric limb-darkened angular diameter measured by Boyajian et al. (2015) and the Gaia DR3 distance, which yields R star = 0.818 ± 0.012 R ☉ .We choose to place wide, uniform priors on the photospheric component temperature parameters, T 1 , T 2 , and T 3 , in order to allow for a broad range of spectral components.In the 1comp model, the covering fraction of the spectral temperature component T 1 is implicitly defined as f 1 = 1.The covering fraction ( f 2 ) corresponding to the second spectral temperature component is set to a uniform prior ranging between (0.0, 0.5), spanning the full range of possible covering fractions while also not allowing for degeneracies between the two spectral temperature components.Therefore, in the 2comp model, the total model generated spectral flux is calculated as with F T , 1 l and F T , 2 l being the PHOENIX spectral models for temperatures T 1 and T 2 .In the 3comp model, we place uniform Prior distributions for the 1comp, 2comp, and 3comp parameter retrieval modeling. refers to a uniform distribution where the first and second inputs are the parameter space range. refers to a normal distribution where the first input is the mean and the second input is the standard deviation.The prior distribution on the temperature applies to all temperature components in each of the three models.
priors on r 1 and r 2 over ranges that ensure f 1 > f 2 > f 3 .In the 3comp model, the total model generated spectral flux is calculated as which can be related from r 1 and r 2 using Equations (7), (8), and (9).
The last two free parameters in the model relate to any systematic offset in the observed flux (a scale ) or underestimate in the uncertainties f ln var .Again, we place a wide, uniform prior on a scale to account for systematic offsets that are suggested by the data (see Equation (1)).Likewise, we set a wide, uniform prior on f ln var between [−50, 0], noting the median value of log F 10 photon, obs, ( ) is −3, i.e., the typical photon noise is 0.1% of the flux.The need for including the f ln var parameter stems from the current state-of-the-art stellar models not being able to adequately fit stellar spectra at the exquisite precision afforded by HST, as demonstrated with retrieval modeling in previous works (Wakeford et al. 2019;Garcia et al. 2022;Rackham & de Wit 2023).

Model Inference and Comparison Metric
We derive posterior probability distributions and the Bayesian evidence associated with the model parameters with the nestedsampling Monte Carlo algorithm MLFriends (Buchner 2016(Buchner , 2019) ) using the UltraNest package (Buchner 2021). 15We employ the SliceSampler functionality to explore the parameter space and set the maximum number of improvement loops to three to arrive at a balance with computational runtime, as discussed in Rackham & de Wit (2023).We take the number of steps to be 2× the number of parameters in the given model, but also verified that increasing the steps to 10× the parameter space provided no significant change in the posterior distribution.At each sampling stage, the speclib (see footnote 12) library generates a model spectra based on the sampled photospheric stellar temperature component drawn, utilizing the design efficiency of the SpectralGrid object.This allows for a computationally faster sampling procedure as the stellar model grid is saved in memory derived from the fixed stellar surface gravity and metallicity inputs, while linearly interpolating between different temperature values.
In order to compare the suitability of the 1comp, 2comp, and 3comp models for the data, we take the Bayesian evidence (ln  ) of each model and calculate the Bayes factor ( ln  D ) in favor of the higher-evidence model.We then convert the Bayes factors to frequentist p-values following Equations (10) and (11) of Benneke & Seager (2013).Table 2 of Benneke & Seager (2013) translates the calculation of the ln  D between two models and the p-value or "significance" of the higherevidence model over the lower-evidence model.For example, ln  D = 5.0 translates to a "p-value" or 0.0003 (3.6σ), indicating a strong preference of the higher-evidence model over the lower.

Results
First, we explore the results of running the 2013 and 2022 HST observations of HD 189733 A through the parameter retrieval sampling.To keep the analysis consistent, we only compare the results from each observing date using the data observed utilizing the same scan direction.

Model Complexity
We first compare the results of the 1comp, 2comp, and 3comp models for each visit and scan direction.We note that in all six cases (three visits × two scan directions) the model retrievals resulted in a ln  D that indicated a strong preference of the 2comp model over the 1comp model.
Second, the ln  D between the 3comp and 2comp models showed a weak preference favoring the 2comp formulation.Based on these findings, we take the 2comp model to be the best spectral description of the HST data (Figure 2), which can be seen in Figure 3, where the 2comp model spectra show a better fit to the HST data compared to the 1comp model spectra.

Observation Date and Scan Direction Variation
We are further interested in whether the best-fit model parameters differ between scan directions and visits.Figure 4 illustrates the posterior samples of the photospheric parameters (T 1 , T 2 , and f 2 ) for the preferred 2comp model in all six retrievals (see Figure B1 in Appendix B for the posteriors of all parameters).The posterior distribution results from Figures 4  and B1 show that the physically relevant parameters are fully consistent within each visit, while a scale varies between scan directions.This underscores both that the flux-level differences apparent in Figure 1 are due to systematics and that our scaling approach accounts for this while not impacting our physical inferences.
Figure 2. Bayesian evidence for different model fits over the three observation dates (rows) and two scan directions (columns).Each panel has the marginal Bayesian evidence for the specified model ("1comp," "2comp," and "3comp) relative to the reference (2comp), which is highlighted in green.Triangles show Bayesian evidences that represent models with >5σ preference against, demonstrating a poor fit.The 2comp model is preferred over all observation dates and scan directions; it provides a significant increase in the marginal Bayesian evidence over the 1comp model, while the additional complexity of the 3comp model does not provide a similar improvement.
From Figure 4, it is also evident that all data sets give evidence for primary and secondary spectral components that agree within 1σ of 5295 41 43  and 3222 K 116 100  , respectively, with a generally consistent (within 2σ) covering fraction of 42% ± 4%.The posterior distributions on the median primary and secondary temperatures are consistent at 1σ across all six observations, suggesting two underlying photospheric components contribute to the spectra from all visits.Importantly, between the six data sets, the spread between the median covering fractions varies from 38% ± 4% to 47% ± 3% (seen The flux level varies between models because the conversion from raw counts depends on the free parameter a scale .Notably, the simplest 1comp retrievals settle on a distinct value of a scale in order to optimize the fit to the spectral slope with a single spectral component.In addition to the mean (gray line) and standard deviation (gray shaded region) of total model flux, contributions from individual spectral components are shown as colored lines (means) and shaded regions (standard deviations).Below each panel with spectral fluxes, another panel shows the residuals between the data and models in units of the standard deviation of the set of posterior models.Over all observation dates, the 2comp model fits the data better than the 1comp model, while the 3comp model does not significantly improve the fit.Note that the higher precision of the 2022 June 21 observations leads to better constraints on the contributions from the individual components (right middle panel).
in Table 3), and this variation on f 2 is discussed in Section 5.The preferred mean parameter retrieval results are summarized in Table 4, where we take the average over all of the median parameter values from the 2comp models.
Lastly, as previously described, we include the f ln var parameter to account for additional scaled model uncertainty, relaying the fact that the current preferred stellar spectral models cannot reproduce observational data from HST to high precision.As seen in the results of our preferred 2comp model (Table 4), we retrieve a f ln 4.06 var = -, indicating the need to artificially inflate the errors on HST data by a factor of ≈50 to match the level of uncertainty that is inherent to the stellar models.

Discussion
We find that the HST/WFC3 spectra of HD 189733 A are best described by a model that includes two spectral components, suggesting that the disk-integrated stellar spectrum comprises the emission from a quiescent photosphere and a substantial coverage of spots.Here, we first contextualize this result by comparison to studies of spot-crossing events on HD 189733 A as well as spectral studies of a comparably active K dwarf.We then discuss whether our results in combination with the light curve of HD 189733 A from TESS can enable constraints on the distribution of spots on HD 189733 A.

Comparison to Results from Spot Crossings
Our best-fit model calls for spots with temperatures of 3222 K 116 100  (Table 4).This result is in tension with those from studies of spot-crossing events on HD 189733 A.  Berdyugina (2005), suggests that a star like HD 189733 A with a photospheric temperature of 5050 ± 50 K (Bouchy et al. 2005) should have ∼3700 K spots.Thus, while estimates from spot-crossing events may be ∼500 K higher than expected, our results are likewise ∼500 K lower than expectations.
This discrepancy may owe to the generally narrow wavelength coverage (1.1-1.7 μm) and low resolving power (R ∼ 130) of the HST/WFC3/G141 spectra we study.In the K-dwarf regime, spectra covering this wavelength range and resolution are essentially sloped spectra with few strong features (Figure 3), which may limit the ability of our retrievals to key into the right temperatures.However, it is unclear why this would lead to a specific bias toward lower spot temperatures.
Considering this broader context, we suggest that our results should be interpreted as a lower limit on the spot temperature for HD 189733 A. By extension, our estimate of the spot covering fraction (0.42 ± 0.04) then also represents a lower limit: Warmer spots have less of a flux contrast with respect to the photosphere than cooler spots, so they would need to cover more of the star for a given quiescent photospheric temperature and hemisphere-averaged flux.Future work with JWST spectra with higher resolution and broader wavelength coverage may shed light on this issue.Additionally, joint analyses of the existing out-of-transit HST spectra and the in-transit spotcrossing events, while outside of the scope of the current study, could help to reconcile the out-of-transit and in-transit inferences.4.
With this caveat noted, we consider the implications of our results at face value in the following sections.

Comparison to σ Draconis
Chromospheric activity is intrinsically linked to photospheric starspots, since both phenomena depend on the strength and distribution of the stellar magnetic field (see review by Hall 2008).In short, concentrated magnetic fields both suppress convection, producing dark starspots in the photosphere, and lead to enhanced heating in the chromosphere, giving rise to bright plage regions (Osterbrock 1961;Skumanich et al. 1975;Spruit 1977).Emission in the Ca II H and K lines traces this chromospheric activity, with a commonly adopted measure being the chromospheric R log 10 HK ¢ index (Linsky et al. 1979), a measurement of the magnitude of emissions from the H and K lines normalized to the total stellar flux across all wavelengths (bolometric flux).
In this context, the star σ Draconis provides an interesting analog for HD 189733 A, as a well-studied, active, early K dwarf (K0V; Keenan & McNeil 1989) =  K with a spot covering fraction of f 2 = 0.22 ± 0.06.This result is similar to our parameter retrieval with the 2comp model of T 1 ≈ 5300 K, T 2 ≈ 3200 K, and f 2,avg = 0.42 ± 0.04 (Table 4).R log 10 HK ¢ scales with spot covering fraction for at least G and K dwarfs (Morris et al. 2019b), so HD 189733 Aʼs higher chromospheric activity suggests a higher spot coverage.Likewise, the same survey also found that the G/K-type stars that had relatively high inferred spot coverages ( f 2  0.3) also had short stellar rotation periods (Morris et al. 2019b).HD 189733 Aʼs rotation period (11.4 days, as suggested by a Lomb-Scargle periodicity analysis of its TESS light curve; Lomb 1976;Scargle 1982) is short compared to the rotation period of σ Draconis (>29 days, based on its projected rotational velocity; Absil et al. 2013), further suggesting that our derived covering fraction has a physical basis.

Constraints on Spot Sizes and Distribution Using TESS
Our analysis of the HST/WFC3 spectra provides an understanding of the contributions to the full-disk-integrated stellar spectrum of HD 189733 A at specific epochs.This analysis, however, cannot tell us about the typical spot size or the spatial distribution of spots on the star.Beyond providing insights into the stellar dynamo, these properties are interesting because they govern the impact of the stellar heterogeneity on the transmission spectrum of HD 189733 b.
To understand what spot sizes and distributions are consistent with our inferences from the HST spectra, we turn to the light curve of HD 189733 A from TESS, which monitored the star over a wider time range.TESS observed HD 189733 A at both 120 and 20 s cadence during Sectors 41 (2021 July 23-2021 August 20) and 54 (2022 July 9-2022 August 5).Both light curves show evolving rotational modulation with a maximum peak-to-peak amplitude of 1.4% (Figure 5).
We adapted the approach of Rackham et al. (2018Rackham et al. ( , 2019) ) to constrain what spot sizes and total coverages are consistent with the TESS rotational variability amplitude.The main modification for our adapted approach is that, in lieu of a rectangular grid, we used fleck (Morris 2020a) to model the stellar surface.fleck is a Python software package for simulating light curves of spotted, rotating stars.It models spots as circular dark regions with a set contrast and accounts for limb darkening and foreshortening toward the limb.More details of the fleck algorithm are provided by Morris (2020b).
Briefly, our approach involves iteratively adding spots to a fleck model photosphere until its peak-to-peak rotational variability reaches the reference level, repeating the process many times to place a statistical constraint on the spot coverage consistent with the reference variability for the given spot parameters.We set the spot contrast to 0.0990, determined by integrating PHOENIX spectra with our best-fit temperatures (Table 4) over the TESS bandpass. 16We tested spot radii from 0.02 R s to 0.20 R s in 0.01 R s increments, and repeated the simulation 10,000 times for each spot radius.
When a simulation reaches the reference variability level, we use shapely (Gillies et al. 2022) to measure full-disk and transit-chord spot coverages.shapely is a Python package for manipulation and analysis of geometric features, underscoring some of the capabilities of fleck.For each simulation, we use shapely to measure at the model zero phase the full-disk spot coverage, the fraction of the projected stellar disk covered by spots, and the spot coverage within the transit chord, the region of the projected stellar disk covered by HD 189733 b over the course of a transit.
Figure 6 shows the results of this analysis.Two effects are important to note.First, both whole-disk and transit-chord covering fractions increase with decreasing spot radius.This is because smaller, more numerous spots contribute less to rotational variability than larger, less numerous spots, as they are less likely to be concentrated in rotational phase.The trend suggests that the spot filling factor approaches a maximum of ∼7% as the spot radius approaches 0. Second, the variance of F chord is larger than that of F disk for larger spot radii but approaches the smaller variance of F disk as spot radii decrease.This suggests that, for randomly distributed spots, larger spots tend to lead to more of a mismatch between the full-disk and transit-chord spot filling factors, which results in TLS signals.
By contrast, as spots become smaller and more numerous, mismatches and thus TLS signals are less pronounced.The primary conclusion from this analysis is that the spottedness we infer from the HST/WFC3 stellar spectra of HD 189733 A cannot be due to spots randomly distributed over the whole stellar surface.Given our inferred spot and photosphere temperatures, we cannot distribute randomly (or uniformly) spots with a filling factor of ∼7% (or larger) without leading to variability that would be higher than what has been observed in the TESS light curve.Thus, the majority of the ∼30% spot coverage suggested by our HST/WFC3 analysis implies that HD 189733 Aʼs spot distribution is strongly axisymmetric.We note that if spots on HD 189733 A are actually warmer than our estimates, the total spot coverage would be higher (given the discussion in Section 5.1), further strengthening the case for an axisymmetric spot distribution.
An axisymmetric distribution could originate from a densely filled latitudinal band of spots or one or more large polar spots (Figure 7).Leveraging the Kepler light curve of its misaligned transiting planet, Morris et al. (2017) showed that the active K4 dwarf HAT-P-11 has spots preferentially at latitudes of ∼16°± 1°, much like the Sun near its activity maximum.Given its similar spectral type (K2V), we should expect a similar dynamo and thus latitudinal concentration of spots for HD 189733 A. However, the chromospheric activity index of HD 189733 A has been observed to vary from     et al. 2016), while that of HAT-P-11 averages −4.35 and has varied from roughly −4.5 to −4.3 over 10 yr (Morris et al. 2017, their Figure 6).Thus, while it is unclear whether R log 10 HK ¢ has been measured over the full activity cycle of either star, the generally lower R log 10 HK ¢ index of HD 189733 A relative to that of HAT-P-11 does not clearly support the expectation of higher spot coverage and more densely filled latitudinal bands on HD 189733 A.
The other possible axisymmetric distribution we consider is a polar spot or spots, the likes of which are widely seen on stars studied through Doppler imaging (see review by Strassmeier 2009, and references therein).However, with a rotation period of 11.4 days, suggested by a Lomb-Scargle periodicity analysis (Lomb 1976;Scargle 1982) of its TESS light curve, HD 189733 A rotates more slowly than a typical star amenable to Doppler imaging studies, and it is unclear whether its rotation rate could produce a long-lived polar spot.Thus, beyond the axisymmetric constraint from TESS, the true distribution of the spot coverage suggested by the HST/WFC3 spectra analysis presented here remains a puzzle.

Conclusion
The HD 189733 A system provides a prime target for studying the atmosphere of a hot giant exoplanet.Constraining the spot coverage of HD 189733 A is crucial for enabling atmospheric studies, and previous works have presented differing interpretations of the planetary transmission spectrum that hinge on the total spot coverage of HD 189733 A. Here, we leveraged the out-of-transit HST/WFC3/G141 spectra of HD 189733 A from three epochs to gain insights into the stellar photospheric heterogeneity and its variability.The main conclusions from our study are as follows.
1. Testing models with one to three spectral components, we consistently find that the Bayesian evidence favors a model with two spectral components.We interpret this as evidence for emission from the quiescent photosphere and a persistent spot coverage on HD 189733 A. Across all observations, we estimate a mean temperature of 5295 41

43
 K for the dominant component (i.e., the quiescent photosphere) and 3222 116 100  K for the minor component (i.e., the spotted area), with a mean covering fraction of the minor component of f 2 = 0.42 ± 0.04.Our inferred spot temperature is lower than estimates from spot-crossing events on HD 189733 A and expectations from Doppler imaging studies of K stars, and we suggest that it should be interpreted as a lower limit on the true spot temperature.2. We find consistent values for T 1 and T 2 across all visits and forward/reverse scan subsets of data.However, while consistent between forward/reverse scan subsets of data, our inferences of f 2 vary at the 2σ level between visits, from 38% ± 4% to 47% ± 3%.We interpret this as a suggestion of variability in the spottedness of HD 189733 A in addition to a baseline level of persistent spottedness.3. Our current inferences are largely dominated by uncertainties stemming from stellar models.The values of f ln var in our fits, ranging from 4.200 0.078 0.071 - to 3.785 0.080 0.078 - , show that the data uncertainties need to be inflated to ∼1%-2% (representing a factor of 28 increase) in order to reach reduced chi-squared values of ∼1 in the model fit.4. Considering models with randomly distributed small spots, we find that the maximal rotational variability of HD 189733 A in its TESS light curve suggests a spot filling factor of no larger than ∼7%, given the spot and photosphere temperatures inferred from our HST/WFC3 spectral analysis.HD 189733 Aʼs spot distribution must thus be primarily axisymmetrical, e.g., densely packed latitudinal bands of spots and/or a large polar spot.
Future work should aim to derive further constraints on the photospheric heterogeneities of HD 189733 A in other data sets, including high-resolution spectroscopy (e.g., Morris et al. 2019a).Additionally, simultaneous modeling of the heterogeneous stellar photosphere and planetary spectrum in transit studies will provide the most consistent constraints of stellar heterogeneity and its impact on transmission spectra for active stars like HD 189733 A.
and SPHINX (Iyer et al. 2023) spectral models, which we detail in Appendix A. In our model retrievals, we fix [Fe/H] = −0.03and g log 4.53 10 ( ) = (Bouchyet al. 2005), as the low-resolution WFC3 spectra are not sensitive to these parameters.We use the speclib Python library(Rackham 2023) to linearly interpolate for spectra within the model grid.12

Figure 1 .
Figure1.Left: median spectra divided by the exposure time for each scan direction observed flux within each observation date data set.Right: median spectra divided by the exposure time for each scan direction observed flux subtracted by the overall (six spectra) mean spectra within each observation date data set.Error bars for each spectra are smaller than the width of the lines.The left figure highlights the offset between the six scan observation spectra.Similarly, the right figure identifies the difference in curvature between the observed spectra, suggesting that a multiplicative scale factor is employed in the model spectra formulation (Equation (1)).

Figure 3 .
Figure 3. 1comp (top row), 2comp (middle row), and 3comp (bottom row) posterior data and models for the 2013 June 5 (left column), 2013 June 24 (middle column), and 2022 June 21 (right column) HST reverse scan observations of HD 189733 A. Spectra are given in units of flux at the stellar surface.The flux level varies between models because the conversion from raw counts depends on the free parameter a scale .Notably, the simplest 1comp retrievals settle on a distinct value of a scale in order to optimize the fit to the spectral slope with a single spectral component.In addition to the mean (gray line) and standard deviation (gray shaded region) of total model flux, contributions from individual spectral components are shown as colored lines (means) and shaded regions (standard deviations).Below each panel with spectral fluxes, another panel shows the residuals between the data and models in units of the standard deviation of the set of posterior models.Over all observation dates, the 2comp model fits the data better than the 1comp model, while the 3comp model does not significantly improve the fit.Note that the higher precision of the 2022 June 21 observations leads to better constraints on the contributions from the individual components (right middle panel).
Pont et al. (2008) find a spot-crossing event in an HST/ACS transit of HD 189733 b is well-modeled by a T ∼ 4000 K spot.In a similar analysis, Sing et al. (2011) find a spot temperature of T = 4250 ± 250 K from an HST/STIS spot occultation.Turning to a wider set of stars, Equation (1) of Rackham et al. (2019), derived from the Doppler imaging studies summarized by Median of the retrieved model parameters for all six data sets (reverse + forward scan data from the three observation dates) using the 2comp model, presented in Section 4.

Figure 4 .
Figure 4. 2comp posterior parameters inferred from the model sampling derived from the forward and reverse scan data sets of the 2013 June 5 (orange), 2013 June 24 (gray), and 2022 June 21 (blue) HST observations of HD 189733 A. Both the forward and reverse scan data sets retrieve similar parameter distributions, discussed in Section 4. These corner plots compromise the model parameters used to describe the best-fit model data, summarized in Table4.
. It has been observed to be a relatively active star, with a R log 4.68 10 HK ¢ = -(Morris et al. 2019b).By contrast, HD 189733 A has been shown to have R log 4.54 10 HK ¢ = -(Melo et al. 2006), indicating that it might be more active than σ Draconis.Morris et al. (2019b), studying σ Draconis, perform a 2comp modeling process similar to our work, focused on fitting observational data to segments of the stellar spectra around the absorption band of TiO.They arrive at estimates of T

Figure 5 .
Figure 5.The TESS 120 s light curve of HD 189733 A. TESS data from Sectors 41 and 54 (black) are shown alongside five examples of rotational light curves produced with fleck (gray) for spot radii of 0.05 R s .The blue bar gives the maximum peak-to-peak amplitude of the observed rotational modulation.We stress that the aim of this analysis is not to fit the evolving light curve entirely but to constrain probabilistically the spot properties consistent with the maximum TESS rotational amplitude.Transits of HD 189733 b have been masked.

Figure 6 .
Figure 6.Maximum spot coverage consistent with HD 189733 Aʼs peak-topeak variability of 1.4% in TESS data as a function of spot radius.Median (points) and 68% confidence intervals (error bars) are shown for the spot coverage of the whole disk and transit chord.Points are offset along the x-axis for clarity.

Figure 7 .
Figure 7. Examples of hypothetical spot distributions with ∼40% spot coverage and ∼1.4% rotational variability in the TESS bandpass.The left panel provides an example of densely filled latitudinal spot bands in which overdensities of spots lead to rotational variability and apparent spot-crossing events, while the right shows large, nonvariable polar spots with a few additional spots that produce the variability and occasional spot-crossing event.The region occulted by HD 189733 b is outlined by dotted white lines in both panels.Made with spotter (Garcia & Rackham 2023).

Figure A1 .
Figure A1.All posterior parameters from the 2comp model sampling derived from the forward and reverse scan data sets of the 2013 June 5 (orange), 2013 June 24 (gray), and 2022 June 21 (blue) HST observations of HD 189733 A, using the (left) PHOENIX and (right) MPS-ATLAS/SPHINX (Witzke et al. 2021; Kostogryz et al. 2023; Iyer et al. 2023) spectral models.Utilizing different spectral grids gives confidence that the parameter retrieval is independent of stellar model.

Table 2
Free Parameters and their Priors for the Retrieval Models Mean of the retrieved model parameters of the six data sets (reverse + forward scan data from the three observation dates) using the 2comp model, presented in Section 4.

Table A1 2comp
Model Parameters over All Dates and Scans Median of the retrieved model parameters for all data sets + forward scan data from the three observation dates) using the 2comp model, comparing the usage of the PHOENIX spectral model (presented in main text) with the MPS-ATLAS (Witzke et al. 2021; Kostogryz et al. 2023) and SPHINX (Iyer et al. spectral models used to verify that the results are independent of spectral model. Note.