The Hubble/STIS Near-ultraviolet Transmission Spectrum of HD 189733b

The benchmark hot Jupiter HD 189733b has been a key target to lay out the foundations of comparative planetology for giant exoplanets. As such, HD 189733b has been extensively studied across the electromagnetic spectrum. Here, we report the observation and analysis of three transit light curves of HD 189733b obtained with {\Hubble}/STIS in the near ultraviolet, the last remaining unexplored spectral window to be probed with present-day instrumentation for this planet. The NUV is a unique window for atmospheric mass-loss studies owing to the strong resonance lines and large photospheric flux. Overall, from a low-resolution analysis ($R=50$) we found that the planet's near-ultraviolet spectrum is well characterized by a relatively flat baseline, consistent with the optical-infrared transmission, plus two regions at $\sim$2350 and $\sim$2600 {\AA} that exhibit a broad and significant excess absorption above the continuum. From an analysis at a higher resolution ($R=4700$), we found that the transit depths at the core of the magnesium resonance lines are consistent with the surrounding continuum. We discarded the presence of \ion{Mg}{ii} absorption in the upper atmosphere at a $\sim$2--4$\sigma$ confidence level, whereas we could place no significant constraint for \ion{Mg}{i} absorption. These broad absorption features coincide with the expected location of \ion{Fe}{ii} bands; however, solar-abundance hydrodynamic models of the upper atmosphere are not able to reproduce the amplitude of these features with iron absorption. Such scenario would require a combination of little to no iron condensation in the lower-atmosphere, super-solar metallicities, and a mechanism to enhance the absorption features (such as zonal wind broadening). The true nature of this feature remains to be confirmed.


Introduction
The hot Jupiter HD 189733 b is one of the earliest transiting exoplanets to be detected (Bouchy et al. 2005), being one of the closest transiting planets to Earth. HD 189733 b's large mass of 1.13 M Jup , large radius of 1.13 R Jup , and short-period 2.2day orbit (Stassun et al. 2017) provide some of the most favorable conditions for atmospheric characterization. Consequently, HD 189733 b has become one of the most extensively observed exoplanets, now constituting a benchmark target for atmospheric characterization via follow-up observations and theoretical modeling.
From space, the Hubble Space Telescope (HST) has been instrumental in probing the optical-to-infrared transmission spectra of this planet. Collective observations with the Advanced Camera for Surveys (ACS), the Space Telescope Imaging Spectrograph (STIS), and the Wide Field Camera 3 (WFC3) have obtained nearly continuous spectra over the 0.3-1.7 µm range Sing et al. 2011;Huitson et al. 2012;McCullough et al. 2014). Complementary infrared (IR) transmission photometry from the Spitzer Space Telescope between 3.6 and 8.0 µm (Désert et al. 2009;Agol et al. 2010;Knutson et al. 2012), and posterior re-analyses of the cumulative data (Pont et al. 2013; have revealed a steep Rayleigh scat-tering slope dominating the optical and near-infrared spectrum. These results are consistent with the presence of high-altitude atmospheric hazes, ruling out a clear atmosphere. Ground-based facilities have further allowed the detection of several spectral features; in particular, high-resolution observations have lead to the detection of narrow absorption features above the haze layer from the H Balmer lines (Jensen et al. 2012;Cauley et al. 2015) and the Na D lines (Redfield et al. 2008;Wyttenbach et al. 2015). These observations constrained the temperature of the thermosphere, revealed strong wind velocities consistent with an eastward equatorial jet (Louden & Wheatley 2015), and detected a transit early ingress, possibly caused by a hydrogen bow-shock ahead of the planet (Cauley et al. 2015). We note however that these detections have been contested in subsequent studies (see Barnes et al. 2016;Guilluy et al. 2020), arguing that the observed signal could be stellar in nature, rather than planetary. Similarly, stellar effects are also suspected to affect the steep UV slope observed by HST for this planet, although the magnitude of such contribution is not fully understood (see, e.g., McCullough et al. 2014;Rackham et al. 2019).
Infrared observations of the day-side emission during secondary eclipse has revealed a lack of thermal inversion and a nearly solar H 2 O abundance in the lower planetary atmosphere Todorov et al. 2014). Furthermore, infrared photometry along the entire orbital phase has measured the brightness temperature variation, constraining the energy budget and the global circulation regime (Knutson et al. 2007(Knutson et al. , 2009(Knutson et al. , 2012. These observations unveiled an efficient energy redistribution of the incident stellar irradiation and an eastward offset of the peak emission, relative to the sub-stellar point. As with the conclusions from optical observations, these measurements are consistent with circulation models exhibiting an eastward advection by an equatorial super-rotating jet with wind speeds of the order of ∼km s −1 at pressures around the mbar level (e.g., Showman et al. 2009Showman et al. , 2013.

Upper-atmosphere Characterization of HD 189733 b
Given the close planet-to-star distance, it is expected that HD 189733 b exhibits atmospheric escape. In fact, the host star, HD 189733, is an active K-dwarf star (radius of 0.75R , mass of 0.79M , effective temperature of 5052 K, Stassun et al. 2017), showing clear evidence of activity and spots during transit events (e.g., Sing et al. 2011;Pont et al. 2013), which makes it a particularly interesting study case for star-planet interactions (e.g., Poppenhaeger et al. 2013;Pillitteri et al. 2015;Cauley et al. 2018).
Lyα transit observation in the far-ultraviolet with HST/STIS detected an extended hydrogen envelope that expands hydrodynamically and escapes to space, powered by absorption of the high-energy stellar flux (X-ray and extreme ultraviolet; Lecavelier des Etangs et al. 2012;Bourrier et al. 2013). Strong temporal variations and the high velocity of the escaping gas suggest the presence of significant star-planet interactions (consistent with the active nature of the host star), which requires an acceleration mechanism on top of stellar radiation pressure, such as the stellar wind. Although later, Guo & Ben-Jaffel (2016) found that the Lyα variability can be simply explained by invoking a thermal population of H i in the planet's upper atmosphere, without a need for complex dynamics involving stellar winds or radiation pressure. In this case, changes to the star's EUV output would explain the different absorption signals. Far-ultraviolet observations with the HST Cosmic Origins Spectrograph have also led to the detection of O i and C ii in the planetary upper atmosphere (Ben-Jaffel & Ballester 2013), indicating that heavy atoms are possibly entrained and dragged upwards in the hydrogen hydrodynamic outflow of the atmosphere.
The Hα absorption feature detection by Jensen et al. (2012) indicates that the upper atmosphere of HD 189733 b is neither in radiative equilibrium nor thermodynamic equilibrium, due to the large ultraviolet flux received by its active host star. In addition, Salz et al. (2018), Zhang et al. (2022), and Guilluy et al. (2020) detected a strong absorption signal from the He i triplet at 10830 Å by means of ground-based high-resolution transmission spectroscopy with CARMENES and GIARPS, respectively. The observed He absorption favors a compact helium atmosphere with a substantial column density and no evidence of He escape. This interpretation of the He observations would seem to be consistent with the aforementioned explanation of the Lyman alpha transit observations by Guo & Ben-Jaffel (2016). The many observations of HD 189733 b across the electromagnetic spectrum have enabled an in-depth characterization of the composition and physical state of both the lower and upper atmosphere of the planet. However, until recently, the nearultraviolet (NUV) window, probing both the exosphere and thermosphere, escaped the spectral coverage probed by these observations. King et al. (2021) analyzed 20 NUV photometric transits of HD 189733 b obtained from the XMM-Newton optical monitor. They detected a transit-depth signal consistent with that in the optical, ruling out extended broad-band absorption towards or beyond the Roche lobe. However, to detect the NUV absorption signature from individual metallic lines, higher spectral-resolution observations are required, like those from HST/STIS, with a resolving power of R ≈ 30, 000 (see Sing et al. 2019;. For the slighly more strongly irradiated planet HD 209458 b,  detected Fe ii absorption in a 100 Å-wide range around 2370 Å, lying beyond the planetary Roche lobe, while finding no evidence for absorption by Mg ii, neutral magnesium (Mg i), nor neutral iron (Fe i). These results suggest that hydrodynamic escape is strong enough to carry heavy atoms beyond the planetary Roche lobe; however, condensate formation sequesters primarily Mg (rather than Fe) in the lower atmosphere of the planet, consistent with microphysical cloudformation models (Gao et al. 2020). For the ultra-hot Jupiter WASP-121b, Sing et al. (2019) reported the detection of strong narrow NUV absorption features, consistent with the singly ionized magnesium (Mg ii) doublet and several lines of ionized iron (Fe ii), emanating from altitudes that exceed the Roche lobe radius of the planet. These results indicate that Mg ii and Fe ii species are not trapped into condensates at depth, and can hydrodynamically escape the atmosphere.
To complete the spectral survey of this benchmark planet and aid the quantitative assessment of the metallicity, ionization state, and outflow rate of the upper atmosphere of HD 189733 b, here we analyzed 15 HST orbits obtained from three STIS nearultraviolet transit light curves. In Section 2 we describe the observations in more detail. In Section 3 we describe our dataanalysis procedure. In Section 4 we present the main results from the NUV data analysis. In Section 5 we place the observations in context by contrasting them with theoretical models of the planet. In Section 6 we speculate on the origin of the features seen in the NUV observations. Finally, in Section 7 we summarize our conclusions.
Each transit observation (one HST visit) consists of five HST orbits. For each visit, the fourth orbit occurs during transit, whereas the first three and last HST orbits occur out of transit. The orbits are broken into exposures (frames) of ∼332 seconds each. The first two visits contain nine frames per orbit (except for the first orbit, which contains one fewer frame to accommodate acquisition observations). The third visit contains one fewer frame per orbit, possibly due to a longer pointing adjustment time after HST returned to operations from a gyroscope failure. This spectrum has been constructed from the mean time-series spectra obtained during the first visit, after discarding bad pixels. The spectrum is color coded by echelle order (see color bar on the right). The inset panels zoom-in around the location of the Mg ii and Mg i resonance lines.
Each frame consists of an echelle spectrum comprising 23 partially overlapping orders, covering the 2300-3100 Å spectral range. Each order contains 1024 wavelength samples, with a resolving power of R = λ/∆λ = 30 000, where the resolution element spans approximately two pixels (or equivalently ∼0.09 Å, or ∼10 km s −1 ). Figure 1 shows a sample spectrum from visit 1. As expected, the flux decreases toward shorter wavelengths, with a clear drop in signal-to-noise ratio (S/N) at wavelengths shorter than ∼ 2550 Å. The stellar spectrum shows many emission features, most of them correspond to Fe ii lines rising from low energy levels, but also the Mg i and Mg ii resonance lines (the strongest line-emission feature in the data, as shown in the inset panels of Fig. 1). In contrast, stellar Fe i features appear in absorption.
The input for our analysis consisted of the x1d CALSTISreduced fits files 1 (version 3.4.2). The reduction comprises dark and flat-field correction, 1D spectral extraction of each echelle order, background subtraction, and wavelength and flux calibration. The results are the 1D spectra (erg s −1 cm −2 Å −1 ) and their uncertainties.
In addition, following Sing et al. (2019), we also considered the Engineering Data Processing System "jitter" files, which contain 28 different engineering measurements of HST's Pointing Control System for each exposure during an observation.

Data Analysis
The core of our analysis is based on the methodology of . We separated our light-curve analysis into two main steps. First we constructed "white" lightcurves by integrating the flux per transit visit and per echelle order to boost the astrophysical and instrumental signals, from which we characterized the instrumental systematics. Then, we divided out the systematics model from the raw data, and combined the systematics-1 https://hst-docs.stsci.edu/stisdhb/chapter-3-stis-calibration corrected data from all visits into spectral light curves binned at different resolving powers, from which we extracted the transmission spectra of the astrophysical signal.

Analysis of Individual Visits
We characterized the instrumental systematics to detrend them from the astrophysical signal and remove them from further analyses. To increase the S/N of the data, we integrated the flux over each echelle spectral order. The implicit assumption is that the instrumental systematics behave similarly over the integrated range (Kreidberg et al. 2014), in this case, over each echelle order.
The analysis begins by masking bad and overly noisy data values. Following the standard practice for exoplanet HST/STIS time-series analyses, we discarded the first orbit from each visit and the first frame from each orbit, since they show significantly stronger instrumental systematics (see, e.g., Sing et al. 2011). Then, we used the overlap between echelle orders to discard data points where the fluxes differ from each other by more than 3σ. Finally, we manually inspected the edges of the echelle spectra to exclude data points with abnormally low fluxes (typically, a handful of data points at the ends of the echelle orders). We then produced raw "white" light curves by summing the flux over each spectral order and propagating the uncertainties accordingly.
We fit the raw light curves with parametric transit and systematics models at each echelle order (denoted by their wavelength λ) as a function of time (t), HST orbital phase (φ), and jitter vector ( j i ): (1) where T λ (t) is a Mandel & Agol (2002) transit model and S λ (t, φ, j i ) is a model of the instrumental systematics. To obtain statistically robust parameter estimations, we computed best-fitting parameters and uncertainties from a Levenberg-Marquardt optimization and a Markov-chain Monte Carlo (MCMC) sampling, respectively. For this matter, we employed the open-source mc3 statistical package 2 , which implements the Snooker Differential-evolution MCMC algorithm of ter Braak & Vrugt (2008).
The transit model parameterizes the astrophysical signal via the orbital parameters, the planet-to-star radius ratios R p /R s (λ), the out-of-transit stellar fluxes F s (λ), and limb-darkening coefficients C i (λ). For the orbital parameters we adopted previously measured values for the orbital period P = 2.21857567 ± 1.5 × 10 −8 day, inclination i = 85 • .71, and semi-major-axis to stellar radius ratio a/R s = 8.84 ± 0.27 (Stassun et al. 2017), which we kept fixed during the fitting and MCMC process. We fit the transit mid-time epoch, applying a Gaussian prior according to the measured optical value (T 0 = 53955.025551 ± 9 × 10 −6 MJD, Bonomo et al. 2017), propagating uncertainties according to the epoch of the HST observations.
For the limb-darkening parameters, we computed theoretical coefficients based on the stellar properties, since the low S/N ratios and gaps in the HST observations do not allow for a good estimation. For this calculation we used the open-source routines of Espinoza & Jordán (2015) to compute limb-darkening coefficients over each wavelength range based on the PHOENIX stellar model (Husser et al. 2013) that best matches the properties of HD 189733: effective temperature T eff = 5052 K, surface gravity log g = 4.49, and metallicity [M/H] = −0.02 (Stassun et al. 2017). We note that, in absence of any better known guess, the adopted limb-darkening coefficents correspond to the photosphere of the star. If a significant fraction of the stellar flux comes from its chromosphere and transition region, as suggested by the stellar line emission, then the used limb-darkening coefficients may not be as representative as expected.
The HST/STIS time-series observations exhibit a number of well documented instrumental systematics: (1) periodic variations in phase with the HST orbital period that are suspected to originate from the telescope's thermal cycle (Brown et al. 2001), and (2) visit-long variations that might be attributed to stellar activity (e.g., Wakeford et al. 2016;Sing et al. 2019). Following the example of previous studies, we model these systematics considering a family of polynomials in time (t) and HST orbital phase. Furthermore, Sing et al. (2019) showed that additional data from HST's Pointing Control System (so called "jitter" data) can help decorrelate instrumental systematics from the astrophysical signal. For each science exposure, the jitter data describe the HST's fine guidance sensor (FGS) coordinates, FGS jitter, right ascension, declination, position angle, latitude, longitude, angle between HST zenith and target, Earth-limb angle, terminator angle, and magnetic field. HST stores the jitter measurements into time-tag ancillary engineering files, recorded at a high temporal resolution. For each science exposure we computed a value for each of the jitter parameters by calculating the median of the jitter values during the exposure. Therefore, our systematics model consists of a polynomial expression up to a quadratic degree in time, up to a quartic degree in HST phase, and a up to a linear degree in one of the jitter vectors: where a k , b k , and c k are the polynomial coefficients to be fit. Note that we considered an individual set of fitting parameters for each echelle order (we dropped the λ-dependence of these coefficients 2 https://mc3.readthedocs.io/ for clarity). Regarding the jitter vectors j i , we considered each of the 21 vectors listed in Table 1 of Sing et al. (2019). t 0 and φ 0 are reference values for the time and phase, which are fixed at the transit mid-time t 0 = T 0 and at the HST mid-phase φ 0 = 0.2, respectively. j i denotes the mean value of the jitter vector j i along the visit.

Systematics Model Selection
Given our limited understanding of the systematics affecting HST observations, we determined the optimal functional form for the systematics model S λ by comparing fits considering all possible combinations of polynomial orders in t and φ up to the degrees shown in Eq.
(2), and each of the 21 jitter vectors, and fits with no jitter vectors. To this end we employed a Bayesian model selection approach (see, e.g., Trotta 2007), where one compares the posterior probability of two competing modelsgiven the data-by computing the ratio of their Bayesian evidences (also known as Bayes factor). The main problem of this approach is that the evidence of a model is not readily calculable. Numerically, calculating the evidence entails evaluating the model likelihood throughout its parameter space (e.g., via a nested-sampling algorithm). To avoid the computational cost of this task, we adopt a simpler model selection approach where we approximate the evidence ratio via the Bayesian Information Criterion (BIC, Liddle 2007): where χ 2 min is the goodness-of-fit parameter corresponding to the maximum likelihood for the given model, k is the number of free parameters, and N is the number of data points. The Bayes ratio is then simply calculated from the difference in BIC values between the competing models B = exp{−(BIC 1 − BIC 2 )/2} (Raftery 1995), which only requires to find the maximum likelihood for the given models. A Bayes factor B > 1 favors model 1 over model 2, or alternatively, the favored model minimizes BIC.
Certainly, these approximations are based on various assumptions, which may not be fulfilled to different extents, such as independent and identically distributed data or (nearly) Gaussian posterior distributions or that the forward model is correct, for example. Therefore, we also considered a second model-selection criterion based on information theory, that is the Akaike Information Criterion corrected for small sample sizes (AICc, Liddle 2007): Just like before, the favored model minimizes AICc.

Model Selection Applied to the HST Observations
For the HD 189733 b data we treated each HST visit independently of each other, where for each visit we fit simultaneously the astrophysical and systematics parameters. Also, for each visit we fit all 23 echelle orders simultaneously, jointly fitting the transit mid-time epoch and adopting the same configuration of the systematics model for all orders. Table 1 reports the BIC and AICc of the global fit to all orders for the favored systematics model (with and without jitter decorrelation). Figure 2 shows an example of the light-curve fits for a selected echelle order. As seen from the raw light curves (top row of Fig. 2  ). The top row shows the raw light curves (one column for each of the three visits). The markers with error bars denote the fluxes and their uncertainties integrated over the echelle order for each frame. The markers are color-coded according to the HST orbit, where the gray colors denote the frames discarded from the light-curve fitting. The bottom rows show the systematics-corrected light curves corresponding to the best-BIC (colored circles) and best-AICc statistics (empty black diamonds), including the jitter decorrelation. The solid lines denote the astrophysical model of the best-BIC (colored lines) and best-AICc models (black lines). Notes. The first two terms in the systematics model column show the polynomial degree for time and orbital phase, the third term shows the name of the jitter vector (which is always modeled as a linear polynomial). The absence of a third term indicates that no jitter vector was included.
variations. Overall, the second visit shows stronger systematics than the first visit, and the third visit shows stronger systematics than the second, which is reflected in the best-fitting BIC, AICc, and reduced χ 2 values (χ 2 red ). In particular, the second and third visit systematics vary on the order of 10% of the raw flux. If we only included the more "traditional" polynomials in t and φ in the fit, we found clear systematic residual trends, which vary in strength for each echelle order. The relatively new practice of jitter decorrelation introduced by Sing vidual visits). The circle markers with error bars denote the transmission planet-to-star radius ratios from each echelle order and visit (see legend), corresponding to the best-AICc fits. The colored dashed lines denote the mean value for each visit. The square marker denotes the XMM-Newton/UVM2 photometric measurement by King et al. (2021). improves the fit for the second and third visits, minimizing both information criteria statistics and residual systematics (bottom row of Fig. 2). In these fits, the AICc imparts a weaker penalty per fitting parameter than the BIC, thus favoring more complex models, which notably improves the χ 2 red values. However, although the BIC and AICc statistics favor different models, the two criteria usually produced consistent transit depths. All following analyses presented here will focus on the best AICc analysis including jitter decorrelation. Figure 3 shows the resulting transmission spectra according to the best-AICc models for each visit and at each echelle order. One clear feature seen in the transmission spectra is an offset between the third visit and the first two. This offset is consistent with the signature of unocculted stellar spots (see, e.g., Rack-ham et al. 2018;Zellem et al. 2017), which is a likely cause given the well known active nature of HD 189733. For each visit, we found that the transit radius ratios vary closely proportional to the square root of the bolometric fluxes of the observations (R p /R s ∝ F s −1/2 ), which is a direct consequence of observations arising from a heterogeneous stellar photosphere including spots (Rackham et al. 2018). The weighted mean of the R p /R s values for each individual visit are 0.1585, 0.1595, and 0.1701; whereas the relative bolometric fluxes between visits are F v1 s /F v2 s = 1.02 and F v1 s /F v3 s = 1.16. Later, the combined analysis in Section 4.1 will show an independent confirmation that we are observing stellar spot contamination in the data.
Once corrected for stellar activity, the analysis of each of the three visits produces a consistent transmission spectrum across the NUV. We obtained a mostly flat transmission spectrum between 2300 Å and 3100 Å, with a planet-to-star radius ratio of ∼ 0.16, which becomes significantly noisier at wavelengths shorter than 2550 Å, due to the much lower stellar flux. These results are also consistent with the previous NUV observation by King et al. (2021), which measured the transit depth of HD 189733 b in the UVM2 photometric band of the XMM-Newton Space Observatory (effective wavelength 2310 Å, width 480 Å).
As a reference, the estimated altitude of the Roche lobe radius as probed during transit is R L1 = 0.44 R s (i.e., perpendicular to the star-planet axis). This value is approximated by 2/3 the distance to the L1 Lagrange point (see, e.g., Vidal-Madjar et al. 2008). All transit depths at each echelle order are well under the Roche lobe radius.

Analysis of Combined Visits
To search for the narrow metal line-absorption features in the transmission spectrum, we need to analyze the observations at a high resolving power, finding a compromise between sufficient S/N per wavelength channel and diluting the narrow absorption features. Therefore, we treated the data at multiple spectral resolving powers ranging from R = 10 − 70 (coarse sampling analysis), and then at a R = 4700 (fine sampling analysis). We proceeded in the same fashion as in , that is, calibrate the wavelength solution to ensure that the spectra are well aligned, use the modeling results from the previous analysis to divide out the instrumental systematics, combine the data from all three visits into a single spectrum, and fit the transit depths at each bin. All future analyses use the best-AICc dataset since it attained χ 2 red values closer to one.

Wavelength Calibration
We used a cross-correlation procedure to calibrate the wavelength solution of each echelle order in each frame and visit. First, we created a master stellar template (enhancing the S/N) by co-adding the spectra from all frames within an echelle order and visit. We then Doppler-shifted the individual-frame spectra and cross-correlated them with their template, searching for the relative velocity with respect to the template given by the shift that maximized the correlation function. Once we obtained the wavelength correction of the individual frames within an echelle order, we repeated the procedure between the template spectra of each visit to find the relative wavelength correction between all spectra. Finally, we compared our wavelength solution with a theoretical model spectrum of HD 189733 (Shulyak et al. 2004) to obtain an absolute wavelength solution at a rest reference frame. Once calibrated, we Doppler-shifted all frames within a given echelle order to a common wavelength solution using a cubic spline interpolation. To minimize the impact of the interpolation, we selected a reference point that minimized the shifts of the set of frames. For a given visit, we found offsets of less than half a pixel between the frames, and in general the wavelength solution for each echelle order follows a similar trend along the visit, with no outliers.

Divide-white Spectral Extraction
To extract a high-resolution transmission spectrum of HD 189733 b, we followed the "divide-white" spectral analysis by Kreidberg et al. (2014), as described in .
In this approach, we first constructed a non-parametric model of the instrumental systematics by dividing the white-light curves (for each echelle order and visit) by their best-fit astrophysical transit model (see Sect. 3.1). Assuming that the instrumental systematics vary weakly with wavelength, each raw spectral data point is divided by the systematics model corresponding to their orbital phase, echelle order, and visit. Data uncertainties are properly accounted for, i.e., using the white-light model's posterior distribution to estimate the uncertainties systematics model and propagating the errors.
In a second step, we split the data into constant resolvingpower wavelength bins, where for each spectral channel we coadded all flux contributions from all echelle orders, and propagated the errors accordingly.
To analyze the spectral light curves we fit a Mandel & Agol (2002) transit model to the systematics-corrected spectral lightcurves in an MCMC run. As previously, we fixed the orbital a/R s , cos(i), and P parameters, and applied a Gaussian prior for the mid-transit epoch based on the white-light fitting results. We also computed theoretical limb-darkening curves at each spectral channel, which we kept fixed during the fit. At each wavelength bin we fit a unique transit depth common to all three visits, but we fit a separate out-of-transit flux value for each visit. Additionally, we incorporated a stellar-spot correction model that we detail in the following section.

Stellar Spot Correction
To account for the spectral transit-depth variations due to unocculted stellar spots, we modeled the stellar flux as a linear combination of a nominal photospheric region and a spotted region. We assumed that a fraction f spot of the observed photosphere is covered by spots with a characteristic temperature T spot , whereas the rest of the photosphere has the nominal effective temperature of the star T eff = 5052 K (with T eff ≥ T spot ). We used the PHOENIX library of high-resolution synthetic stellar spectra (Husser et al. 2013) to compute the flux of the nominal stellar photosphere F s (T eff ) and the spotted region F spot (T spot ). This leads to the following correction term for the observed transit depth d obs (λ): where d(λ) is the true transit depth and (λ) is the contamination spectrum due to stellar hetereogeneities (as defined in Rackham et al. 2018). To evaluate the correction factor at any given T spot , we interpolated in temperature between all available PHOENIX models (ranging from 2300 K to 5100 K) at a fixed metallicity ([M/H]=0.0) and surface gravity (log(g) = 4.5). The stellar-spot correction model thus contains two free parameters for each visit ( f spot and T spot ), which we fit simultaneously along with the other astrophysical parameters during the MCMC. We considered a uniform prior for the f spot values between 0.0 and 0.5 (it would be unlikely to find unocculted spot covering fractions greater than 50%). For the spot temperature we considered a logarithmic-uniform prior in ∆T spot = T eff −T spot because the spectral variation of the spot correction is more pronounced at small ∆T spot values . When ∆T spot 1500 K the correction becomes nearly achromatic. Regardless of the parameterization, the main conclusions of our analysis do not change when we consider a uniform prior in T spot .

Coarse-binning Spectra
To boost the characterization of the astrophysical signal, we started with the analysis of the combined data binning over a relatively coarse wavelength sampling. Figure 4 shows the NUV transmission spectrum of HD 189733 b sampled at a resolving power of R = 50. This resolving power resulted in the best compromise between having a good S/N per bin and having sufficient spectral resolution to identify variability with wavelength. Appendix A shows the corresponding systematics-and stellar-spot corrected lightcurves for this analysis. At this resolving power, we obtained a relatively flat transmission spectrum with the exception of two regions where the transit depth increases (by more than ∼1σ relative to the surrounding wavelengths) located at λ ∼ 2350 Å and λ ∼ 2600 Å.
For completeness, we also analyzed the data at other similar resolving powers of R = ∼10, 30, 70 finding consistent results. Likewise, analyses with the bins shifted by half a bin-width returned consistent results as well.
The stellar spot contamination was well characterized, as shown in Figure 5. The first two visits show a negligible spot contamination ( < 1.05), with the models finding either a small spot covering fraction or a spot temperature close to the nominal stellar effective temperature (∆T spot 10). In contrast, the third visit shows a clearly significant spot contamination ranging from (λ) ≈ 1.20 to 1.15 from short to long wavelengths. The contamination shows a weak variation with wavelength, which is not larger than the 1σ credible intervals of the constraints. Accordingly, we constrained the spot covering fraction to f spot > 0.11 and the spot temperature to 100 K ∆T spot 500 K, values which are strongly correlated (see posterior pair-wise distribution in Fig. 5). Analyzing the combined data without a stellarspot model produced slightly larger transit depths (Fig. 4); however, the shape of the spectrum remains consistent with that of the spot-corrected analysis, which is consistent with the nearly wavelength-independent stellar spot correction.
All in all, the ratio between contamination spectra of the visits (e.g., v3 / v1 = 0.15 − 0.2) matches precisely the ratio between bolometric fluxes (e.g., F v1 s /F v3 s = 0.16), which reinforces the hypothesis that the variability between visits is caused by a varying stellar spot coverage at different epochs.

Fine-binning Spectra: Magnesium Resonance Lines
We attempted to analyze the data binning at much higher resolving power such that we can capture absorption of individual metallic lines, however, the high-resolution lightcurves were largely dominated by noise, and not much could be concluded from the resulting spectra. The only region where we were able to place significant constraints at high resolution is at the core of the magnesium resonance lines, which have fluxes nearly ten times larger than at other wavelengths, and thus a much better S/N. For this, we selected a bin resolving power of R = 4700 (65 km s −1 wide bins), which encapsulates the core of the Mg i and Mg ii resonance lines. Figure 6 shows the high-resolution light curves centered at the Mg i and Mg ii resonance (top), and the stellar spectra (middle) and transit radius-ratio spectra (bottom) around the magnesium lines. For the Mg ii bins, the light curves show a clear transit signature (albeit with some residual noise for the Mg ii h For context, in the bottom panels of Fig. 6 we also show the theoretical upper-atmospheric model for this planet (see Sect. 5.1 for details). While HD 189733 b's escape rate is not high enough to generate magnesium absorption signatures close to the the Roche-lobe altitude (R L1 /R s = 0.44), the model predicts a noticeable excess absorption above the continuum. Thus, we used the binned model estimate to place non-detection limits for the Mg ii measurements. We found that the observations sit 2-4σ below the prediction for a solar-abundance atmospheric model containing Mg in the upper atmosphere (Table 2). Due to the higher noise at the Mg i line, we cannot reject nor confirm the presence of Mg i absorption from our observations.

The NUV-to-Infrared Transmission Spectrum of HD 189733 b
For context, the Hei infrared metastable triplet and Hi Lyman-α transit observations of HD 189733 b can be interpreted in terms of the planet hosting a compact upper atmosphere; extended, but not as much as that of other exoplanets observed in the NUV, <1σ † number of standard deviations below the expected R p /R s from the solar-abundance theoretical model. e.g., HD 209458 b. Metals like magnesium or iron may condense in the lower atmosphere and even if they make it to the upper atmosphere, a lower escape rate would not allow their density profiles to stretch out to high altitudes. At longer wavelengths, the optical-IR observations show a strong slope towards the blue, which can be interpreted as absorption from high-altitude hazes.
To test these hypotheses and to get a better understanding of the interaction between the upper and lower atmosphere of HD 189733 b, we compare the observations to theoretical models. Since it is not possible to model all required upper-and lower-atmosphere properties at once, we model individually the upper and lower regions, and then study how they relate to each other.

Upper-atmosphere Modeling
To interpret the NUV transmission spectrum of HD 189733 b, we employed the physical model of the upper atmospheres of closein planets (Koskinen et al. 2013a(Koskinen et al. ,b, 2022. Since the parts of the NUV transmission spectrum that probe the upper atmosphere are dominated by signatures of ionized metals Fe ii and Mg ii (Sing et al. 2019;, we added the heavy elements C, O, N, Mg, Si, Fe, S, Ca, Na, and K with solar abundances to the model, in addition to H and He, with related chemistry and physics adapted from Huang et al. (2017). We included Mg, Si, Fe, Ca with first and second ionization states, due to the relatively low ionization potential of the first ionization state, and the rest of the heavy elements with their first ionization state. The electron density is equal to the sum of the ion densities under the assumption of quasineutrality. The lower boundary of the escape model is at 0.1 µbar where the temperature is sufficient to dissociate molecules that are not included in the model.
We used the results from the photochemical model of Lavvas & Koskinen (2017) to obtain the mean molecular weight profile in the lower and middle atmosphere to set the lower boundary altitude and temperature of the escape model. The temperature profile in the lower and middle atmosphere was not calculated self-consistently by the photochemical model and was instead adapted from Moses et al. (2011). Here, we used the MK profile from Lavvas & Koskinen (2017). The application of the photochemical model to HD189733b did not include Mg and Fe that are detectable in the NUV. In order to crudely model the full profiles of the absorption lines, we added Mg and Fe to the lower and middle atmosphere, assuming solar abundances and thermal ionization according to the Saha equation (Menou 2012). This approximation does not affect our conclusions because our focus in this section is on the possible upper atmosphere signatures.
Given the temperature profile assumed for HD 189733 b here and assuming equilibrium condensation, Mg, Fe, Si, and Ca are expected to condense to form mineral clouds and rain out from the atmosphere (e.g., Wakeford & Sing 2015). We included Mg and Fe in the escape model simply to explore if the modeled ab- Fig. 6. Transit light curves (top) and spectrum (bottom) of HD 189733 b of the combined HST/STIS transits at the location of the magnesium resonance lines, and binned at a resolving power of R = 4700 (i.e., roughly the width of the magnesium lines). The top panels show the systematicscorrected light curves for the spectral channels centered at the core of the Mg ii and Mg i lines (see labels). The solid red curve and orange area denote the best-fitting transit depths and span of their 1σ uncertainties. The middle panels show the system flux around the magnesium lines, with the vertical dashed lines marking the edges of the spectral bins. The bottom panels show the transit planet-to-star radius ratio for each of the bins around the magnesium resonance lines (red markers with 1σ error bars). The blue curve shows the upper-atmosphere model (the diamond markers show the model binned at R = 4700). The peach-colored area shows the span of the low-resolution transit radius ratio (R = 50). sorption lines of Mg II and Fe II would be detectable in the NUV data. The rest of the heavy elements were retained for consistency with the no-condensation scenario. Although efficient mixing in the lower atmosphere could interfere with mineral cloud formation and allow for higher abundances of the relevant heavy elements in the upper atmosphere (Spiegel et al. 2009;Koskinen et al. 2013b), we do not hold a view on whether this is possible on HD 189733 b or not at this point. Figures 4 and 7 show the transmission spectrum from the upper atmosphere model in blue. Here, we calibrated the model spectrum so that the measured base R p /R s (in the near-IR) corresponds to 2.3 mbar in the atmosphere (Lavvas & Koskinen 2017). So as not to bias the fit to any particular interpretation at this point, we did not include extinction by the high-altitude haze modeled by Lavvas & Koskinen (2017) in this forward model. Instead, the transit continuum is due to Rayleigh scattering by H 2 , H, and He. It is immediately clear that the model continuum falls below the extrapolation based on the optical-IR transmission spectrum of the planet, due to the lack of the high-altitude haze in the forward model. Since transit depth is not a cumulative quantity, however, the model still correctly represents absorption by Fe II and Mg II in the upper atmosphere in the absence of condensation, even if high-altitude hazes are present.
We ran the upper atmosphere model assuming globally averaged radiative forcing with the XUV spectrum constructed from Eridani as an input, due to the well known similarity to HD 189733 (Lavvas & Koskinen 2017). The reference simula-tion, on which the spectrum in Figure 4 is based on, assumes an isolated planet. We also ran a second simulation that includes Roche lobe overflow (see, Koskinen et al. 2022) but this had negligible effect on the modeled spectrum, in line with the expectation that Roche lobe overflow does not significantly enhance the escape rate or upper atmosphere densities on HD 189733 b.
We note that NUV observations have the potential to offer valuable new insights to the properties of possible highaltitude hazes and/or clouds. Unfortunately, S/N is an issue for HD 189733 b that precludes clear conclusions on the origin of the continuum in the transit observations. The visible spectrum shows a strong slope indicative of small particle hazes Lavvas & Koskinen 2017) that appears to be supported to some degree by observations by the Stratospheric Observatory for Infrared Astronomy (SOFIA) (Angerhausen et al. 2015); however, the latter found a slope offsetted to lower depths by ∆R p /R s =∼ 0.001 (see Appendix A). The slope can, however, be enhanced by stellar activity (star spots) . Also, a general circulation model (GCM) that simulates the production and transport of hazes predicts a shallower slope unless mixing by sub-grid scale eddies exceeds mixing by global circulation (Steinrueck et al. 2021). The NUV observations are statistically consistent with the model that includes no haze, and generally the observed transmission falls below the extrapolation based on the visible slope (Figure 7), with the exception of the 2350 Å region. Due to low S/N in the NUV, however, the observations are also consistent with this extrapolation. Thus, these NUV observations neither confirm or exclude the presence of high-altitude hazes, especially if the visible slope is enhanced by star spots.

Lower-atmosphere Modeling
To constrain the lower-atmosphere properties of HD 189733 b we used the open-source Pyrat Bay framework 3 (Cubillos & Blecic 2021). We modeled the atmosphere from 100 to 10 −9 bar, adopting a parametric temperature profile from Madhusudhan & Seager (2009), hydrostatic equilibrium, and a H/He-dominated composition considering constant-with-altitude volume mixing ratios.
The atmospheric retrieval employed the differentialevolution Markov-chain Monte Carlo sampler (ter Braak & Vrugt 2008), implemented via the open-source code mc3 (Cubillos et al. 2017), guided by the transit observations from HST's STIS G430 and G750, WFC3 G102 and G141 spectrographs; and from Spitzer's IRAC 3.6, 4.5, 5.8, 8.0 µm and MIPS 24 µm photometric bands (Pont et al. 2013;. Figure 7 shows the fit to the optical-IR transit observations. Appendix B presents the retrieval posterior-distribution results for all model parameters. The optical-IR observations point to a prominent haze slope, steeper than a Rayleigh slope (α haze ≈ −9.3±2.3). This is a similar conclusion to that found by previous retrieval studies on this planet (e.g., Pinhas et al. 2019;Barstow 2020). From a physical perspective, the steep haze slope can be created by a combination of particle growth and efficient eddy (turbulent) mixing (Ohno & Kawashima 2020) or by stellar inhomogeneities (Rackham et al. 2018), effects that were captured by the parametric haze model (rather than being explicitly accounted for). The optical absorption is also stronger than that produced by hydrogen (as in the upper-atmosphere model), leading to a higher continuum level at optical wavelengths, as well when extrapolated to NUV wavelengths. The NUV observations are consistent with the optical-IR continuum model, with the exception of the few data points that extend above the continuum (mainly the 2350 Å region). Lastly, we can see that the constraint on the temperature profile becomes less precise above the µbar level since optical/IR data is mostly probing the lower layers of the atmosphere. At the same time, at this high altitudes the retrieved temperature profile starts to differ from the higher temperatures predicted by self-consistent models that account for high-energy flux that heats the thermosphere of the planet (Lav-3 https://pyratbay.readthedocs.io/ vas & Koskinen 2017), i.e., the models used as the base of our upper-atmosphere model.

Fe ii Absorption
It is interesting that the excess NUV absorption matches precisely the locations of the two strongest Fe ii bands in the region, while there does not seem to be any absorption corresponding to the Fe i band shortly blueward of the magnesium doublet. Unfortunately, the low S/N achieved when binning the data at higher resolution does not allow us to identify absorption from individual lines.
In fact, a detection of Fe ii would be very controversial in this planet for multiple reasons. In first place, iron is not expected to survive the cold trap at the lower layers. Even if we neglected the iron condensation, as is the case for our upper-atmosphere model, we cannot match the large amplitude of the observed features when we bin the model at the resolution of the observation. Another argument against Fe ii is that given the relatively cool expected temperature of the planet and the late spectral type of the star, one would also expect to observe Fe i if Fe ii were detected.
If the NUV absorption indeed corresponded to Fe ii, it would imply that our hydrodynamic model severely underestimates the absorption signature, and we would be in need of a mechanism to enhance the lines such that they become detectable in broader bins. This would require a significant broadening of the iron lines and super-solar abundances, even when assuming no condensation in the lower atmosphere. However, global circulation models (GCM) predict winds of the order of ±8 km s −1 at the terminator (Steinrueck et al. 2021). These models include regions of significant size with flows both away from and towards the Earth due to the complexity of the circulation and rotation effects. These kinds of flows could effectively broaden metal lines in spatially unresolved transit observations.

Absorption from an Unknown Absorber or other Source
The broad-band absorption suggests the presence of an extra absorber, but the fact that the absorption matched the location of Fe ii bands might be deceptive. It is certainly possible that the absorption comes from another unidentified species, possibly from either a haze or from the electronic band system of a molecular species, given the width of the absorption and the transition energies involved.
Following Lothringer et al. (2022), we explored whether SiO, say vaporized from strongly mixed cloud particles, could appear in the upper atmosphere. However, the atmospheric retrievals do not favor significant SiO absorption when it was included in the model. From a further heuristic exploration, we could not find any configuration of temperatures and abundances producing SiO bands at the observed strength nor location. We finally tested for other candidates available in the Ex-oMol database and that extended into the NUV, which led us to discard both HS and OH.
One might consider that the absorption originates from an external source such as accreted material from trojan or circumplanetary satellites (Kislyakova et al. 2016;Oza et al. 2019). However, this is a highly speculative scenario, since to date there is no model that demonstrates that material from an external source can produce such large absorption features, nor that it could generate exclusively an Fe ii signature but not that from other particles like magnesium.
It seems also unlikely that the absorption features come from a yet-undetermined aerosol. To our knowledge there is no selfconsistent GCM that couples haze production and circulation to a degree that it can state with confidence that particles could form and exist at the required high altitudes (∼0.1 nbar). It is hard then to explain how aerosols would form and remain aloft at such altitudes.
We also tested whether the stellar spot correction could generate the observed features (Rackham et al. 2018(Rackham et al. , 2019, but on close inspection of the spectral dependency of the spot contamination at different temperatures, we found no indication that it could mimic extra absorption at the location of the Fe ii bands. Finally, we considered the impact of stellar center-to-limb flux variation. As the stellar Fe ii lines show in emission (Fig.  1), the observations should be probing the (hotter) stellar chromosphere and transition region rather than the (cooler) photosphere. Therefore the lines may be showing limb brightening instead of limb darkening, leading to a quite different transit signature (see, e.g., Schlawin et al. 2010). Unfortunately, the limb-darkening/brightening law for these lines is unknown, and it cannot be reasonably estimated with current instrumentation. In any case, while the unknown limb-darkening law of the emission lines introduces a degree of uncertainty in the NUV transit depths, unaccounted limb brightening would lead to underestimates transit depths, and thus does not explain the excess absorption seen at 2350 Å and 2600 Å.

Exoplanet NUV Observations in Context
To date three exoplanets have been spectroscopically observed during transit in the NUV with HST/STIS: WASP-121b (Sing et al. 2019), HD 209458 b (Cubillos et al. 2020, and HD 189733 b (this work). These three planets span a wide range of irradiation regimes, receiving a stellar flux of 364× (HD 189733 b), 769× (HD 209458 b), and 5624× (WASP-121b) that of the one received by the Earth. The analysis of this aggregated data provides a first glance at the physical processes that shape the properties of the planets' upper atmospheres.
The upper panel of Fig. 8 shows the transmission spectra of these three planets at a resolution of R = 1300 in units of r/R p , with R p being the optical radius. We can see that WASP-121b exhibits enhanced transit depths coinciding with Fe ii and Mg ii lines, whereas the other two planets show subtler features. When comparing the transmission spectra at a resolving power of R = 130 in units of r/R RL , with R RL being the Roche-lobe terminator radius, (lower panel of Fig. 8) we see that the pseudocontinuum extends to the Roche-lobe boundary, which indicates that the middle atmosphere/lower thermosphere of WASP-121b undergoes Roche-lobe overflow. In clear contrast, this is not the case for HD 189733 b nor HD 209458 b. WASP-121b is therefore unique among these planets, and must have a much higher mass loss rate than the others, as expected. This broad behavior is qualitatively consistent with expectations based on model predictions.
Article number, page 11 of 17 A&A proofs: manuscript no. nuv_HD189733b The lower panel of Fig. 8 illustrates the dramatic difference between WASP-121b and the two other planets. The NUV spectra of HD 209458 b and HD 189733 b look nearly featureless when compared to that of WASP-121b. Certainly, WASP-121b is not an ordinary planet-it must have substantially larger massloss rates than the other two planets. Neither HD 209458 b nor HD 189733 b show the Mg ii absorption feature that survives in the spectrum of WASP-121b. As we already know, the spectrum of HD 189733 b shows a strong feature close to the blue Fe ii absorption band. This feature does not extend to the Roche lobe, but the peak is ∼150 surface scale heights above the optical radius. Given that heating in the thermosphere increases the scale height by a factor of ∼20 above the ∼1 µbar level, the peak probably probes the 1 nbar level.
Finally, we note that the feature in the spectrum of HD 189733 b is not identical to the Fe ii feature in the spectrum of WASP-121b though. In particular, the point at ∼2400 Å coincides with a strong Fe ii band and should show maximum absorption towards the longer-wavelength edge. This is true for WASP-121b, but in the spectrum of HD 189733 b absorption drops instead. This poses a conundrum, however, given the uncertainties in the data analysis, noise could be responsible for this apparent inconsistency.

Conclusions
HD 189733 b is the coolest hot Jupiter that has been spectrosopically observed in the NUV so far. Transmission observations indicate that this is a very hazy planet, but the composition of the hazes are not really know, some of them may present additional absorption in the NUV.
By analyzing the combined NUV HST/STIS transmission spectra (3 transits) at a relatively coarse resolution (R = 50), we obtained better S/Ns and spectral resolution than previous NUV observation of HD 189733 b in the NUV by King et al. (2021). The transit signature in the light curves is clear and with no significant residual systematics. We were able to constrain a continuum level that sits in between the extrapolation of the optical haze slope and the continuum of a pure-H 2 /He upper atmosphere model. Thus, given the uncertainties of our continuum data points, the NUV data is statistically consistent with either a hazy or clear atmosphere scenario. In addition, the NUV transmission spectrum presents strong and broad absorption features that coincide with the location of two strong Fe ii bands, while there is no apparent excess absorption correlated with Fe i bands. When analyzing the data at a higher resolution (R = 4700), we found that most of the lightcurves are dominated by noise, which does not allow us to confirm nor discard iron as the source of the absorption. However, thanks to the much higher stellar flux at the location of the magnesium resonance lines, we were able to rule out the presense of Mg ii excess absorption above the continuum. When compared to a solar-metallicity model with magnesium in the upper atmosphere, the measured transmission R p /R s values at the core of the Mg ii lines are 2-4σ below the model predictions. The S/Ns of the data at the Mg i line is not high enough to reject or confirm excess absorption due to magnesium.
HD 189733 b has a surface gravitational acceleration twice that of HD 209458 b and unlike WASP-121b, Roche-lobe overflow effects are not expected to be that significant. Its atmosphere is also cooler, suggesting that the formation of clouds is expected. This is well in agreement with optical observations, where the data are consistent with a high-altitude haze, which could have a soot-like precursor. HD 189733 b upperatmosphere models assuming that magnesium did not condense in the lower atmosphere show considerable excess absorption at the core of the Mg ii resonance lines. Thus, the absense of excess absorption suggests that magnesium is not escaping on HD 189733 b and that Mg-Si clouds likely form in the lower atmosphere.
At this point, the nature of the broad NUV absorption features is not known. An iron origin raises significant challenges for existing upper-atmosphere models. In first place, given the relatively low temperatures of HD 189733 b, condensation is expected to sequester heavy metals into the lower layers of the atmosphere, though it is possible that condensates form primarily from other metals instead, such as magnesium (Gao et al. 2020;Woitke et al. 2020). However, even if iron is not strongly depleted, one would still need a mechanism to enhance its absorption beyond that predicted by our upper-atmosphere models. The two most immediate possibilities would be to consider higher metallicities and zonal-wind velocity broadening of the absorption lines. GCM estimations of HD 189733 b predict wind velocities at ∼8 km s −1 at the terminator, the question for future studies is whether these mechanisms are enough to boost the signal. Heavy metals play a key role in the modelling and interpretation of NUV transit observations and hence in our understanding of exoplanetary atmospheres as a whole. The detection of metals at Roche lobe distances would be surprising instead of expected for HD 189733 b. Here we presented the detection of a broad absorption feature in the NUV that matches the location of Fe ii bands, but cannot fully confirm the nature of the absorber, both future observational and theoretical work are needed to understand more precisely the origin of this absorption. Fig A.1 shows the NUV transmission spectrum binned at a coarse resolving power of R = 10 and shows a SOFIA-like offset of the transit depth for comparison. Figure A.2 shows the systematics-and stellar-spot corrected lightcurves of the combined HST/STIS HD 189733 b observations, when binned over a coarse wavelength sampling with a resolving power of R = 50. We used the posterior-distribution median and the central 68% percentiles statistics to estimate the parameter values and their credible intervals.
Tables A.1, A.2, A.3 present the estimated values of the NUV transmission planet-to-star radius ratio when analyzed at resolving powers of 50, 33, and 10, respectively. combined HST/STIS observations. The red markers with error bars denote the systematics-and stellar-spot corrected transmission spectrum, their 1σ uncertainties, and the span of the spectral bins for a resolving power of R = 10. The blue spectrum shows a theoretical model of the planet's upper atmosphere. The orange curve shows the fit to the optical-IR transmission spectrum extrapolated into the NUV (the shaded areas denote the 1σ and 2σ uncertainties). The black dashed curve shows the optical-IR fit shifted downwards according to a SOFIA-like offset of ∆R p /R s = 0.001. Table B.1 shows the parameterization, priors, and retrieved posterior values for the analysis of the HD 189733 b optical-IR transmission spectrum. Figure B.1 shows the posterior distributions of the model parameters and temperature profile. Transit light curves of the combined HST/STIS observations of HD 189733 b binned at a resolving power of R = 50. The colored markers with error bars denote the systematics-and stellar-spot corrected measurements, color coded for each visit (see legend). The solid red curve and orange area denote the best-fitting transit depths and span of their 1σ uncertainties. The lightcurve of each bin has been shifted in the vertical axis for visualization. The labels next to each ligtcurve denotes the mean wavelength of the spectral bin. Note that the scale of the vertical axis changes from left to right as the S/N of the data improves with increasing wavelength.