The geometric albedo of the hot Jupiter HD 189733b measured with CHEOPS

Context. Measurements of the occultation of an exoplanet at visible wavelengths allow us to determine the reﬂective properties of a planetary atmosphere. The observed occultation depth can be translated into a geometric albedo. This in turn aids in characterising the structure and composition of an atmosphere by providing additional information on the wavelength-dependent reﬂective qualities of the aerosols in the atmosphere. Aims. Our aim is to provide a precise measurement of the geometric albedo of the gas giant HD189733b by measuring the occultation depth in the broad optical bandpass of CHEOPS (350 – 1100 nm). Methods. We analysed 13 observations of the occultation of HD189733b performed by CHEOPS utilising the Python package PyCHEOPS . The resulting occultation depth is then used to infer the geometric albedo accounting for the contribution of thermal emission from the planet. We also aid the analysis by reﬁning the transit parameters combining observations made by the TESS and CHEOPS space telescopes. Results. We report the detection of an 24 . 7 ± 4 . 5 ppm occultation in the CHEOPS observations. This occultation depth corresponds to a geometric albedo of 0 . 076 ± 0 . 016. Our measurement is consistent with models assuming the atmosphere of the planet to be cloud-free at the scattering level and absorption in the CHEOPS band to be dominated by the resonant Na doublet. Taking into account previous optical-light occultation observations obtained with the Hubble Space Telescope, both measurements combined are consistent with a super-stellar Na elemental abundance in the dayside atmosphere of HD189733b. We further constrain the planetary Bond albedo to between 0.013 and 0.42 at 3 σ conﬁdence. Conclusions. We ﬁnd that the reﬂective properties of the HD189733b dayside atmosphere are consistent with a cloud-free atmosphere having a super-stellar metal content. When compared to an analogous CHEOPS measurement for HD209458b, our data hint at a slightly lower geometric albedo for HD189733b (0 . 076 ± 0 . 016) than for HD209458b (0 . 096 ± 0 . 016), or a higher atmospheric Na content in the same modelling framework. While our constraint on the Bond albedo is consistent with previously published values, we note that the higher-end values of ∼ 0.4, as derived previously from infrared phase curves, would also require peculiarly high reﬂectance in the infrared, which again would make it more di ﬃ cult to disentangle reﬂected and emitted light in the total observed ﬂux, and therefore to correctly account for reﬂected light in the interpretation of those phase curves. Lower reported values for the Bond albedos are less a ﬀ ected by this ambiguity.


Introduction
The reflective properties of an exoplanetary atmosphere are quantified by the geometric albedo (Russell 1916), which is defined as the albedo of the planet at full phase.The global energy budget of a planet is dependent on how much light from the host star is able to enter the atmosphere without being reflected at its top.The reflectivity of the atmosphere can also be wavelength-dependent as different atmospheric constituents absorb and reflect light at different wavelengths differently.This results in the geometric albedo being dependent on the spectral range of the observed light (Sudarsky et al. 2000;Heng & Demory 2013;Parmentier et al. 2016).The measurement of the albedo therefore provides an additional observational constraint when modelling the atmospheric structure and composition of an exoplanet.
In practice, the geometric albedo can be determined by observing the secondary eclipse (occultation) of the planet at optical wavelengths.While observations of the occultation at infrared wavelengths are limited to the thermal contribution of the planetary emission, measurements at optical wavelengths also allow for the characterisation of the reflective contribution.The total occulted flux is a combination of thermal emission from the planet and reflected stellar light.By accounting for the thermal contribution to the total occulted flux (e.g.Cowan & Agol 2011;Heng & Demory 2013;Wong et al. 2020Wong et al. , 2021)), the geometric albedo can then be inferred from the remaining reflective contribution.
Such measurements have been done with a variety of space telescopes.Both Angerhausen et al. (2015) and Esteves et al. (2015) report geometric albedos for 20 and 14 planets, respectively.They determined the albedos with phase curve observations performed by the Kepler Space telescope.Both studies find that geometric albedos of gas giants in the Kepler bandpass (420-910 nm; Koch et al. 2010) are usually < 0.1.
These findings are also supported by Heng & Demory (2013), who obtained A g < 0.15 for 9 out of 11 planets looked at in their study.The two exceptions are HAT-P-7b and Kepler-7b, for which they report A g = 0.225 ± 0.004 and A g = 0.352 ± 0.023, respectively.The case of the high albedo of Kepler-7b, however, has been widely discussed, and Heng et al. (2021) revised the value to A g = 0.25 +0.01  −0.02 .The high albedos for these planets are attributed to the presence of clouds or condensates in the atmosphere.All of these observations also prove the capability of space-grade photometry to detect low amplitude occultation signals with, for example, a measured occultation depth of 10.9 ± 2.2 ppm for TrES-2b and 16.5 ± 4.5 ppm for Kepler-8b (Kipping & Spiegel 2011;Angerhausen et al. 2015).Wong et al. (2020Wong et al. ( , 2021) ) report detections of secondary eclipses at optical wavelengths for 15 planets from data acquired with the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014).Due to design and purpose, the photometric precision of TESS is smaller than that of Kepler.The smallest occultation depth detected with at least 3σ confidence in these studies is that of WASP-100b at 94 ± 17 ppm (Wong et al. 2021).However, the reported uncertainties are often an order of magnitude higher for other systems.They also infer geometric albedos for the planets with a detected occultation, while being careful to account for the thermal contribution to the total occulted flux.For planets with dayside temperatures below 1500 K, they find that with increasing temperature the geometric albedo decreases because fewer condensates are created in the atmosphere.However, for planets with dayside temperatures between 1500 and 3000 K, they find a weak positive correlation of the dayside brightness temperature and the geometric Fig. 1.Wavelength-dependent response functions of CHEOPS (dashed red), Kepler (dotted blue), and TESS (solid green).The response functions are from the Spanish Virtual Observatory Filter Profile Service (Rodrigo et al. 2012;Rodrigo & Solano 2020).
albedo.This means that the reflective qualities of planets with very high temperatures improve again despite the temperature increase.They attribute this to high-temperature condensates, which cause higher atmospheric reflectivity.Alternatively they also suggest that opacity sources may contribute additional unaccounted emission at visible wavelengths (Cowan & Agol 2011), which lead to biased high albedos.
The Characterising Exoplanet Satellite (CHEOPS; Benz et al. 2021) has been used to detect occultations of several planets.Lendl et al. (2020) reported the first detection of an occultation (WASP-189b) with CHEOPS.Brandeker et al. (2022) used CHEOPS to observe occultations of the hot Jupiter HD 209458b at optical wavelengths.They prove the capability of this space telescope to also detect shallow secondary eclipses by reporting a detection of an occultation in the CHEOPS bandpass with a depth of 20.4 +3.2 −3.3 ppm.Apart from the higher photometric precision, CHEOPS also provides a wider bandpass (350-1100 nm) than TESS (580-1120 nm).The response functions for CHEOPS, TESS, and Kepler are shown in Fig. 1.The CHEOPS and Kepler bandpasses are very similar.The TESS bandpass does not cover the bluer parts of the visible spectrum, but is much more sensitive at infrared wavelengths and therefore to thermal emission from planets, while Kepler and CHEOPS are more sensitive to reflection.Bouchy et al. (2005) reported the discovery of the hot Jupiter HD 189733b, which was among the first exoplanets for which both radial velocity and photometric transit measurements became available.The combination of these two complementary methods allowed for a first characterisation of the planet using its density.It has since become a benchmark for studying hot Jupiters and their atmospheres.Deming et al. (2006) were the first to report the detection of an occultation of HD 189733b in data obtained with the Spitzer Space Telescope.Knutson et al. (2007Knutson et al. ( , 2009) ) used phase curves observed by Spitzer in the 8 and 24 µm channels to study the day-night contrast of the planet.They reported a maximum brightness temperature of 1212 ± 11 K at a wavelength of 8 µm.They also find a small day-night contrast and a bright-spot offset in both wavelength channels, indicating that the energy absorbed at the day-side of the planet is efficiently redistributed across the atmosphere.Agol et al. (2010) used a total of seven occultation observations by Spitzer to study phase variations of the infrared flux emitted by the planet.Most A24, page 2 of 16 Krenn,A. F.,et al.: A&A proofs, notably they reported no significant eclipse depth variations over the 2-yr baseline of their observations.The set of Spitzer phase curve and occultation observations was complemented by additional observations in the 3.6 and 4.5 µm channels (Knutson et al. 2012).They also reported no significant eclipse depth variations and present enhanced absorption in the 4.5 µm channel as evidence for the existence of vertical mixing resulting in excess of CO.
Thermal phase curves can be used to derive estimates of the Bond albedo of a planet by measuring day-and night-side temperatures and computing the equilibrium temperature of the planet (Cowan & Agol 2011;Schwartz & Cowan 2015).The Bond albedo is a measure of the total reflectance of a planet across all wavelengths and in all directions.It is specifically important when determining the energy budget of a planet.This has been done several times for the case of HD 189733b using the Spitzer thermal phase curve observations.Schwartz & Cowan (2015) and Schwartz et al. (2017) report consistent measurements of A B = 0.37 +0.04  −0.05 and A B = 0.41 ± 0.07 respectively.Notably, Zhang et al. (2018) show that applying the same approach used in Schwartz et al. (2017) to only Spitzer 3.6 and 4.5 µm phase curves (discarding the 8 and 24 µm phase curves), yields a lower value (but still consistent within 1σ) of A B ≈ 0.3 ± 0.1.All of these works use similar approaches (error weighted means) to convert broad-band brightness temperatures to day-and nightside effective temperatures.Keating et al. (2019), on the other hand, retrieve a significantly lower estimate of A B = 0.16 +0.11   −0.10 from the same datasets using a Gaussian process regression developed by Pass et al. (2019) to create a wavelength-dependent brightness map, which in turn is used to estimate effective temperatures.Evans et al. (2013) studied the occultation of HD 189733b at optical wavelengths using observations by the Hubble Space Telescope (HST).They were able to detect the occultation signal in the 290-450 nm wavelength channel at 126 +37 −36 ppm, but were unable to detect it in the redder 450--570 nm wavelength channel where the signal is expected to be significantly smaller.They provide an upper limit of the geometric albedo of A g < 0.12 across 450-570 nm at 1σ confidence, while determining a geometric albedo across 290-450 nm of A g = 0.40 ± 0.12.They interpreted the decrease in the albedo towards longer wavelengths as an indication of the presence of optically thick reflective clouds on the dayside.In particular, they suspect that sodium absorption suppresses the scattering of light from the red part of the visible spectrum.
Using the 390-435 nm and 435-480 nm channels of the HST data published in Evans et al. (2013), Wiktorowicz et al. (2015) computed a geometric albedo of HD 189733b in the B-band (390-480 nm) of A g = 0.226 ± 0.091.They supported their computation with polarimetric observations, which allowed for the determination of a 3σ confidence upper limit of A g < 0.37 in the B-band (in this case 391-482 nm).Combining these photometric and polarimetric observations they reject a B-band geometric albedo of A g = 0.61 ± 0.12, which was previously reported by Berdyugina et al. (2011) using different polarimetric observations.Their result is consistent with models including the presence of Rayleigh scattering in the atmosphere.
In this work, we present observations of the occultation of HD 189733b at optical wavelengths performed by CHEOPS with the aim of utilising the exquisite precision of the space telescope for bright stars to determine the geometric albedo of the planet.Additionally, we also present transit observations of the target performed by TESS and CHEOPS that aid in constraining the planetary parameters in the course of the analysis of the occultation observations.Section 2 contains a description of the acquired data, and in Sect. 3 the analysis of the data is described in detail.In Sect.4, using the retrieved occultation depth, we provide an estimate of the geometric albedo of the planet and discuss it in the context of other geometric albedo measurements of similar planets.Finally, we also derive lower and upper boundaries for the Bond albedo of HD 189733b and compare them with previous measurements of the Bond albedo from thermal phase curves.

Description of acquired data
2.1.Transit data HD 189733 was observed by TESS in Sector 41 of cycle 4 of the extended mission at 2 min cadence.In our analysis we used the Pre-search Data Conditioning Simple Aperture Photometry (PDC-SAP) provided by the Science Processing Operations Center (SPOC) pipeline (Smith et al. 2012;Stumpe et al. 2014;Jenkins et al. 2016).CHEOPS also performed two transit observations of the target in August 2021 (see Table 1).The CHEOPS observations are available as two different data products (Benz et al. 2021): sub-arrays, which contain a circular region around the target with a radius of 100 pixels and are a product of combining three individual 10.45 s exposures resulting in an effective cadence of 31.35 s, and imagettes, which contain circular subsections of a 25-pixel radius around the target and are available at a cadence equal to the exposure time of 10.45 s.Aperture photometry is available only for the sub-arrays via the official CHEOPS Data Reduction Pipeline (DRP; Hoyer et al. 2020).However, PSF photometry can be performed on the imagettes using PIPE 1 (see also Morris et al. 2021;Szabó et al. 2021;Brandeker et al. 2022).
We processed all of the CHEOPS transit observations with the DRP using an aperture radius of 30 pixels (RSUP aperture radius option).This option was chosen as it ensures that the whole PSF is included in the aperture radius, and minimises the number of nearby contaminating stars.Notably, we also performed both the transit and occultation analysis using data processed with PIPE, and achieved consistent results with both reduction alternatives.
In general, CHEOPS observations are affected by instrumental noise such as stray light from the Earth and the Moon (Moon glint), smearing effects, or spacecraft jitter.The flux measurements show a particularly strong correlation with the spacecraft roll angle (see also Lendl et al. 2020;Bonfanti et al. 2021).The spacecraft is designed to rotate around itself exactly once every orbit.Therefore, the roll angle parameter is directly linked to the orbital position of the spacecraft.Instrumental noise must be accounted for during the data analysis in order to identify the transit and occultation signals of the planet (see Sect. 3).Prior to performing the transit analysis, we removed all of those points that were flagged by the DRP; this includes those points that are contaminated, for example by cosmic rays.Additionally, we performed a sigma clipping and removed all points with the median absolute deviation (MAD) higher than 5 to discard outliers.

Occultation data
CHEOPS observed HD 189733 ten times in a time period from 30 June 2021 until 16 August 2021, and three times from 10 July 2022 until 2 August 2022 (see Table 1).Each visit contains a single occultation event at the middle of the observation with out-of-occultation data being acquired both before and after the occultation.Each individual visit comprises either six or seven CHEOPS orbits.A single CHEOPS orbit covers roughly 100 min.However, the target cannot be observed throughout the entire orbit, due to Earth occultations and South Atlantic Anomaly (SAA) crossings.This leads to gaps in the observations, with a width that varies from visit to visit, spanning values between 20 and 40 min, depending on target coordinates and observation time.
For the occultation observations, the dataset of each individual visit was also processed with the standard CHEOPS DRP version 13.1.0(Hoyer et al. 2020).Again an aperture radius of 30 pixels (RSUP aperture radius option) was used during processing.To remove outliers, we applied a MAD clipping with a clipping factor of σ = 2.5 on all datasets.Since the observed roll angle trends in the data are partially caused by changes in the amount of background light, we opted for a clipping of data points with high background values.All points for which the background estimate >700 ADU were removed.This procedure especially removes data points shortly before and after an Earth occultation, when stray light from Earth's atmosphere hits the telescope, as well as data points with increased flux levels due to a crossing of the SAA.Finally, we also removed data points with a smear estimate larger than 9 × 10 −5 times the median flux value of the whole dataset as well as all data points with at least one centroid coordinate being shifted more than 1 pixel from the median centroid position.Both the threshold for background clipping and the threshold for smear estimate clipping were determined by inspecting plots showing the correlation of the corresponding parameters with the observed fluxes.The thresholds were chosen in such a way that they represent a tradeoff between ensuring that data points affected very strongly by instrumental systematics are clipped and minimising the number of clipped points.

Stellar properties
To aid our transit and occultation analyses, we provided a prior on the stellar density of the host star HD 189733.The stellar density was calculated from stellar radius and mass estimates computed specifically for this work.To derive the radius we used a modified Markov chain Monte Carlo (MCMC) infrared flux method (IRFM; Blackwell & Shallis 1977;Schanche et al. 2020) with spectral priors taken from values available in the SWEET-Cat catalogue (Sousa et al. 2021; for more details see Wilson et al. 2022).To retrieve a mass estimate we assumed the stellar T eff , [Fe/H], and radius as the basic set of input parameters to derive the isochronal mass and age values from two different stellar evolutionary models.To retrieve the first pair of mass and age values we used the interpolation capability of the isochrone placement algorithm (Bonfanti et al. 2015(Bonfanti et al. , 2016) ) to fit the input parameters within pre-computed grids of PARSEC2 v1.2S (Marigo et al. 2017) isochrones and tracks.In particular, we added v sin i = 3.5 ± 1.0 km/s (Bouchy et al. 2005) to the basic set of input parameters because the synergy between isochrones and gyrochronology improves the routine convergence, as detailed in Bonfanti et al. (2016).A second pair of mass and age values was computed by the Code Liégeois d'Évolution Stellaire (CLES; Scuflaire et al. 2008), which generates the best-fit stellar track according to the basic set of input parameters following the Levenberg-Marquadt minimisation scheme (Salmon et al. 2021).As discussed in Bonfanti et al. (2021), we finally checked the mutual consistency of the two respective pairs of outcomes through a χ 2 -based criterion, and we merged the mass and age distributions to obtain the results used in this study.All adopted stellar parameters including the stellar density are presented in Table 2.
To aid our interpretation of the measured geometric albedo we also determined the stellar sodium abundance from the publicly available optical spectrum obtained with the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph using the derived and adopted stellar parameters.We adopted the classical curve-of-growth analysis method assuming local thermodynamic equilibrium.We used the ARES v2 code3  Notes.The Gaia ID of HD 189733b is 1827242816201846144.(Sousa et al. 2015) to measure the equivalent widths of the spectral lines.Then we used a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973) to convert the EWs into abundances closely following the methods described in Adibekyan et al. (2012Adibekyan et al. ( , 2015)), among others.
Given the low temperature of the star, it was not possible to determine the abundances of C, N, and O directly from the HARPS spectrum.The abundances of these elements have been estimated using the available abundance datasets of solar  2) into absolute abundances we adopted the solar reference values from Asplund et al. (2021).This yielded the following elemental abundances of carbon, oxygen, nitrogen, and sodium: C/H = (2.4 ± 0.7) × 10 −4 , (O/H = 6.6 ± 2.7) × 10 −5 , C/H = (5.1 ± 1.3) × 10 −4 , Na/H = (1.6 ± 0.2) × 10 −6 .These values are used as input to our computation of the planetary geometric albedo (see Sect. 4.1).

Transit analysis
Before jointly analysing the transit observations of TESS and CHEOPS, we analysed these datasets individually.This was mainly to select and constrain the systematic and astrophysical noise models.In this analysis we did not assume any prior knowledge of planetary parameters except the period and the transit time, for which we used Gaussian priors based on values from Baluev et al. (2015).We chose this option because the Baluev et al. (2015) analysis uses a longer photometric baseline and radial velocity measurements to constrain these parameters, and therefore using these priors provides a more precise value than an analysis using only photometric observations by TESS and CHEOPS.We also adopted a Gaussian prior on stellar density (see Sect. 3.1 and Table 2).For the rest of the planetary parameters, we used wide uninformative priors (see Table 3).To parametrise the limb-darkening effect, we used the quadratic limb-darkening law with the efficient coefficient parametrisation proposed by Kipping (2013).We used the code of Espinoza & Jordán (2015) to compute Gaussian priors on the limb darkening coefficients (see Table 3).
The individual TESS light curve was modelled using juliet 4 (Espinoza et al. 2019), which uses a transit model from batman 5 (Kreidberg 2015).In addition to the transit model, we added a jitter term, the mean out-of-transit flux, and a Gaussian process (GP) model to account for systematic and/or temporal astrophysical trends.The GP model was built using the Exponential-Matérn kernel in juliet, based on celerite 6 models (Foreman-Mackey et al. 2017).In practice, we first analysed only out-of-transit data with instrumental and GP models in order to better constrain the nuisance parameters and to obtain a better and fast convergence.We then used posteriors from this analysis as priors when we modelled the whole TESS dataset along with transit parameters.
CHEOPS data are known to correlate with many instrumental parameters, mainly the roll angle, as the CHEOPS field of view rotates around the target during each of its orbits (see Sect. 2).To correct for these trends we decorrelated the dataset against several instrumental parameters with the PyCHEOPS7 Python package, which was specifically developed to detrend and fit CHEOPS data (Maxted et al. 2022).We chose a set of detrending basis vectors to be included in the analysis by individually adding them one by one and retaining an additional vector only when supported by a higher Bayes factor computed from the Bayesian information criterion (BIC).This way we ended up using up to third-order harmonics of the roll angle and a second-order polynomial in the PSF centroid position.Both of the CHEOPS visits showed the well-known ramp effect caused by thermal effects in the instrument (Maxted et al. 2022).This effect manifests itself in the reduced data as an increased flux at the beginning of the observation, which quickly decreases to the same normalised flux level as the rest of the observation.We detrended this effect using the telescope-tube temperature as a correlated variable by adding a linear detrending model for the thermfront2 parameter provided by the DRP to the detrending basis vectors.We then used juliet to actually fit the data with linear decorrelation against these parameters, and also a GP model to fit for the temporal trends.We followed the same twostep procedure to model the individual CHEOPS visits as we did with the TESS data.
In the end, we modelled one sector of TESS data and two visits of CHEOPS jointly.In addition to the joint transit model, the model that we used contained a jitter term, an out-of-transit offset and a GP model for systematic and astrophysical trends for each dataset.The GP models were again built using the Exponential-Matérn kernel in juliet.We also included linear models to take care of instrumental correlations in two CHEOPS visits.The priors on these noise parameters were Gaussian and were based on our earlier analysis of individual datasets.We used nested sampling methods (Skilling 2004(Skilling , 2006;;Higson et al. 2019), as implemented in juliet via dynesty8 (Speagle 2020) to sample from posteriors.The dataset and the median model, along with other randomly drawn models from posteriors, for CHEOPS and TESS are illustrated in Fig. 2. The residuals show no significant long-term trend indicating the quality of our fit.The median posterior parameters, along with the 1σ confidence interval, from our analysis are listed in Table 3.We opted for the juliet Python framework to fit all of the transit data as it allows for the simultaneous fitting of light curves from different instruments and the use of GPs with a variety of different kernels.Within the transit analysis PyCHEOPS was only used to determine the best detrending basis vectors for the CHEOPS data as it is especially designed to do so.

PyCHEOPS analysis
We used the PyCHEOPS python package (Maxted et al. 2022) to perform a combined analysis of all occultations.The package was chosen as it is optimised for fitting CHEOPS data.There was no need for the functionality of a multi-instrument analysis as the occultation dataset consists of only CHEOPS observations and we did not need the option of choosing between a variety of GP kernels as we did not use GPs in the occultation analysis (see Sect. 3.3.2).We defined Gaussian priors according to the results of the transit analysis (see Sect. 3.2) on timing, period, depth, width, and impact parameter.For the stellar density we assumed the same prior as for the transit analysis (see Sect. 3.1 and Table 2).
The data of each individual visit were corrected for the ramp effect (see Sect. 3.2) with the remove_ramp function built into PyCHEOPS (Maxted et al. 2022).We removed long-term time trends, with characteristic timescales much larger than one CHEOPS orbit, with a second-order polynomial.For the polynomial fit we excluded the in-occultation data in order to not remove the planetary signal.Periodic flux changes correlated to the roll angle parameter (and consequently to the orbital position of the spacecraft) were treated with a two-step process.Firstorder linear models of the sines and cosines of the roll angle θ and the angle 2θ were used to remove large-scale trends.Additionally, a 30-segment spline was fit on top of these models.The spline fits small-scale flux changes as a function of the roll angle to remove glint effects.We again exclude the in-occultation data from the spline fit in order to not remove the occultation signal.The glint-model is scaled by fitting an individual scale factor for A24, page 6 of 16 each visit during the analysis.Analogously to the transit observations, we added additional linear models to detrend for smear noise, which is caused by electrons being smeared over several pixels during the read-out process of the CCD, x-and y-centroid offsets, which quantify the spacecraft jitter, contamination from background stars, and changes in background flux only if their addition was supported by a higher Bayes factor for the specific model.
Since we still observed flux changes correlated to time at characteristic timescales of about one CHEOPS orbit, we decided to apply additional time detrending.To do so we first used a subset of each individual visit containing only out-ofoccultation data and ran the MCMC method implemented in PyCHEOPS to fit for the parameters of the corresponding detrending model.The flux values of these subsets were then corrected for the instrumental systematics using the median values of the posterior distributions of these MCMC fits.Then we fitted a spline to the corrected fluxes.The spline fitted to the out-ofoccultation data of each individual visit was set up to have four segments; it used a fourth-order polynomial in each segment to account for the observed temporal flux-variations.Subsequently, the spline was added to the second-order polynomial that was applied to the dataset prior to fitting.Values of the spline during the occultation were interpolated from values before and after the occultation to ensure that the planetary signal was not removed.
We decided to remove the sixth visit, performed on 20 July 2021, from our multi-visit analysis because as a result of Earth occultations this visit does not cover either ingress or egress.Both the second-order polynomial to remove long-term time trends and the spline fit to account for time trends at characteristic timescales of one CHEOPS orbit are fitted only on out-of-occultation data in order to not remove the planetary signal.Both functions are interpolated during the occultation event.Because this dataset is missing both ingress and egress, there are gaps before and after the in-occultation data.This results in a longer interpolation window than for other visits, making the fitted time trends less reliable.
We then fitted each individual visit for the planetary parameters including the occultation depth, a white noise term, and the parameters of the detrending models using again the MCMC method implemented in PyCHEOPS.The fitted values of the detrending parameters of the individual visits were used as initial guesses of a multi-visit MCMC fit of all 12 remaining occultation observations.We fitted for a single occultation depth, the remaining planetary parameters, a white-noise term, and individual detrending parameters for the detrending models and the glint scale for each visit.The parameters of the roll angle trend models were not explicitly calculated, but implicitly marginalised over in the MCMC following Luger et al. (2017).The multi-visit fit used a total of 53 free parameters.All fits also accounted for a light travel time delay of 31 s (Agol et al. 2010;Knutson et al. 2012).
We measured an occultation depth in the CHEOPS bandpass (330-1100 nm) of 24.7 ± 4.5 ppm.The close-up detrended, phase-folded, and fitted light curve showing the occultation is plotted in Fig. 3.The entire phase-folded light curve and the individual light curves can be found in the Appendix.To provide further evidence in favour of our result, we repeated the whole procedure assuming a model without an occultation (i.e. an occultation depth equal to 0).When comparing the two models we find that Baysian evidence supports the model with an occultation over the model without an occultation.

Residual noise
After the analysis of the acquired data, we observed some remaining residual noise: mainly short-term variations in stellar flux.These flux variations were particularly challenging to account for as their amplitudes and timescales were similar to that of the occultation signal, which has a duration of 108.6 min.We attempted to fit GPs (see e.g.Gibson et al. 2012;Foreman-Mackey et al. 2017), including the use of different kernels, on both the time and roll angle parameters as an alternative to the models previously used.We found that using a GP was not a suitable detrending approach for these specific observations because of the similarity of the amplitudes and timescales of the residual noise and the ingress-egress of the occultation.The GP model has no well-constrained GP hyperparameters, and we found that it tends to overfit the occultation events because it is not able to adequately distinguish possible stellar variations and the occultation signal.We decided not to remove flux variations with typical timescales smaller than one CHEOPS orbit and no significant correlation to any of the other instrumental parameters.
We further investigated if the residual noise could be due to the stellar granulation signals, as in Delrez et al. (2021) and Sulis et al. (2023).Unfortunately, our CHEOPS observations contain too many gaps and are too short (see Sect. 2.2) to properly identify and characterise this stellar noise source when studying periodograms 9 .We then only compare the amplitude (RMS) of our residuals with the expected amplitude of this stellar noise.The RMS of the entire time series is about 95 ppm, while if we bin the data over 4 min (since the power spectrum on <4 min timescales is dominated by photon noise), we get a typical RMS of around 32 to 45 ppm.These amplitudes are similar to the amplitudes predicted by 3D hydrodynamic (HD) simulations of stellar granulation, which are 39.3 ppm (Rodríguez Díaz et al. 2022).However, without accessing the frequency constant of this residual noise, we cannot fairly conclude on a stellar origin.A white Gaussian noise of 95 ppm RMS could also lead to an RMS of around 40 ppm when binned over 4 min.The nature 9 According to Sulis et al. (2023) we would expect an increase in the power spectrum from high to low frequency with a cutoff frequency f g < 600 µHz.
A24, page 7 of 16 A&A 672, A24 (2023) (instrumental or stellar origin) of this residual noise is therefore not well understood with the data at hand.

Geometric albedo
The total occulted flux is composed of two major contributions: the thermal emission of the planet, which is determined by its equilibrium temperature, and the stellar light that is reflected by the planetary atmosphere.Therefore, to determine the reflective contribution the thermal emission has to be estimated.Following Brandeker et al. (2022), we used the dayside temperature of 1192 ± 9 K, estimated with occultation measurements in the Spitzer 4.5 µm channel (Knutson et al. 2012), and extrapolate it to the CHEOPS bandpass.In the course of the extrapolation we assumed an irradiated atmosphere model provided by Mollière et al. (2015) with the parameters T eff = 1250 K, [Fe/H] = 0.0, [C/O] = 0.70, planetary log g = 3.0, and stellar spectral type K5.To model the stellar emission spectrum, we used a synthetic spectrum drawn from the PHOENIX 10 library (Husser et al. 2013).The library contains synthetic spectra for a grid of stellar parameters.We chose the spectrum computed for a star with the parameters T eff = 5000 K, [M/H] = 0.0 and stellar log g = 4.5 as it most closely resembles HD 189733 in the PHOENIX sample.The thermal contribution to the occultation depth is determined to be 1.42 ± 0.03 ppm.
The geometric albedo A g is then computed as (Brandeker et al. 2022) with a being the semi-major axis of the planet, R p the planetary radius, and L the measured occultation depth corrected for the thermal contribution.In particular, a/R p was determined from the result of the transit fit (Table 3).The geometric albedo of HD 189733b is calculated as The measured value is consistent with the upper limit of A g < 0.12 measured in the 450-570 nm Hubble Space Telescope, Space Telescope Imaging Spectrograph (HST STIS) band reported by Evans et al. (2013).The achieved uncertainty of the measurement is identical to the uncertainty of the geometric albedo in the CHEOPS bandpass of HD 209458b for which Brandeker et al. (2022) report A g = 0.096 ± 0.016.

Atmospheric Na content
The measured low geometric albedo of HD 189733b in the optical wavelength range makes this planet part of the hot Jupiter population (T eq = 1000-1800 K) for which low geometric albedos were measured in the Kepler and CHEOPS optical wavelength range (Fig. 4).This would make HD 189733b like HD 209458b a typical dark and cloud-free 'Class IV roaster' (Sudarsky et al. 2000;Seager & Sasselov 2000), a planet for which the low albedo in the optical is consistent with a model that assumes that atmospheres are cloud-free at the pressure level where Rayleigh scattering occurs on the dayside 11 and also 10 https://phoenix.astro.physik.uni-goettingen.de 11Rayleigh scattering in hot Jupiters occurs typically at ≈10 mbar (Vidal-Madjar et al. 2011).Angerhausen et al. 2015, blue) and CHEOPS (Brandeker et al. 2022 and this work, orange) vs effective temperature.The green line at A g = 0.12 approximates the 'dark roaster' regime with very low optical geometric albedos (Sudarsky et al. 2000).Kepler-7b is a prominent outlier, even when taking into account a recent update on its albedo (alt, Heng et al. 2021).Its unusually bright albedo may be explained by reflection of clouds at the planetary limbs (Adams et al. 2022).
assumes for the absorption in the CHEOPS (and in the Kepler) band to be dominated by the resonant Na doublet.
Following the same approach as in Brandeker et al. ( 2022) taken from Heng et al. (2021), we computed the theoretical geometric albedo of HD 189733b as a function of atmospheric Na content using the above model.We assumed a dayside temperature of 1200 K consistent with Spitzer secondary eclipse observations by Knutson et al. (2007), Agol et al. (2010), andTodorov et al. (2014).Furthermore, a cloud-free hydrogen dominated atmosphere was assumed, where we considered H 2 O and Na as dominant absorbers in the CHEOPS bandpass as outlined in Brandeker et al. (2022).Using elemental abundances for the host star as listed in Table 2, we derived for H 2 O a volume mixing ratio of 5.4 × 10 −4 in chemical equilibrium using Heng & Tsai (2016).
Comparing the theoretical A g calculations for different X Na contents integrated over the CHEOPS, HST STIS, and TESS bandpass (Fig. 5) shows that our CHEOPS albedo measurement is consistent with sodium elemental abundances between 1 − 3× its stellar value.When combining the measured CHEOPS and HST STIS albedos (Evans et al. 2013), a super-stellar sodium elemental abundance (X ≈ 3) is required for the observations to be consistent with the calculated model values.The corresponding albedo in the TESS bandpass would be very low, and thus probably not detectable.

Bond albedo
Occultation observations at optical wavelengths can also be used to infer lower and upper boundaries of the Bond albedo (Rowe et al. 2006).They can therefore complement Bond albedo measurements obtained from thermal phase curves.Schwartz & Cowan (2015) show that by assuming either zero or perfect reflection for all wavelengths outside of the observed bandpass, lower and upper limits for the Bond albedo can be derived from the spherical albedo (A S ) in the observed bandpass.The spherical albedo also covers all directions, but in contrast to the Bond A24, page 8 of 16 Krenn, A. F., et al.: A&A proofs, manuscript  albedo, it is defined at specific wavelengths.Following Schwartz & Cowan (2015), the Bond albedo can be computed as where A opt S is the spherical albedo at optical wavelengths, f opt is the ratio of stellar flux emitted at optical wavelengths to the total emitted flux, and A oob S is the spherical albedo for all wavelengths outside of the observed bandpass.To derive lower and upper limits on the Bond albedo, A oob S is assumed to be 0 and 1, respectively.In the case of the CHEOPS bandpass, this results in where f CH is the ratio of stellar flux emitted in the CHEOPS bandpass to the total emitted flux, and A CH S is the spherical albedo of the planet in the CHEOPS bandpass.Using the same PHOENIX spectrum as in Sect.4.1 to model the stellar emission, we find f CH = 0.637 for HD 189733.The geometric and spherical albedo are linked to each other by with q depending on the exact reflective qualities of the atmosphere (Pollack et al. 1986;Burrows & Orton 2010).Heng et al. (2021) derive q = 0.761 for Rayleigh single scattering and isotropic multiple scattering under the assumption of the single scattering albedo to be unity.As the low geometric albedo measured in the CHEOPS bandpass implies Rayleigh scattering to be dominant at the height in the atmosphere where optical light is scattered (see Sect. 4.1), we apply q = 0.761 to transfer our measured geometric albedo in the CHEOPS bandpass to a spherical albedo: A CH S = 0.058 ± 0.012.
Accounting for 3σ confidence intervals of our geometric albedo measurement from CHEOPS, this results in All previous measurements of the Bond albedo of HD 189733b using thermal phase curves are consistent with the lower and upper boundaries derived from the CHEOPS geometric albedo.Using the same approach, Schwartz & Cowan (2015) report a slightly higher value for the lower limit of A min B = 0.043 from secondary eclipse observations in the Hubble 290-450 nm channel.However, this value does not take into account confidence intervals, and assumes q = 5/4 and a black-body spectrum for the stellar emission.In comparison, when assuming our median value for A CH g and A oob S = 0, we retrieve A min B = 0.037.Although the measured geometric albedo derived from shortwavelength Hubble data was high (Evans et al. 2013), those data add little information when constraining the Bond albedo.The bandpasses mostly overlap with only the wavelength range of 290-350 nm being covered by Hubble alone.The integrated stellar emission flux at these wavelengths amounts to only about 1.4% of the total stellar emission and is therefore negligible.We note that the PHOENIX model predicts less emission at very short optical and long UV wavelengths compared to a black-body spectrum.
Schwartz & Cowan (2015) also suggest assuming A oob S = 0.5 for hot Jupiters, which would result in A 0.5 B = 0.22 ± 0.08.On the other hand, it is also possible to take Bond albedos inferred from thermal phase curves and calculate the corresponding A oob S values.Since the stellar fluxes at wavelengths shorter than the CHEOPS band are mostly negligible for HD 189733, A oob S can also be used as a proxy of the spherical albedo at infrared wavelengths.Assuming A B = 0.37 ± 0.05 from Schwartz & Cowan (2015) results in a peculiarly high out-of-band albedo of A oob S = 0.92 ± 0.16.Using A B = 0.16 ± 0.11 from Keating et al. (2019) results in a more moderate out-of-band albedo of A oob S = 0.34 ± 0.32.Since the measured light from a planet is composed of both light emitted by the planet and reflected stellar light, observations of thermal phase curves have to account for both reflected and emitted light.Figure 6 shows the estimated contribution of reflected light to the total planetary flux for different grey geometric albedos as a function of wavelength for the HD 189733 system for a given planetary dayside temperature.Here we used the same PHOENIX and irradiated atmosphere models as in Sect.4.1 to model the stellar emission and planetary atmosphere respectively.Analogously to determining the thermal emission of the planet in the CHEOPS bandpass, we again extrapolated thermal emission at a given wavelength from the dayside temperature of 1192 ± 9 K, estimated with occultation measurements in the Spitzer 4.5 µm channel (Knutson et al. 2012).Contrary to the usual assumption that the vast majority of infrared flux is contributed by thermal emission from the planet, in the case of a high geometric albedo the reflective component at longer wavelengths remains significant.If this contribution to the occultation depth is neglected or underestimated when interpreting infrared observations, this leads to overestimated thermal flux, and thus to overestimated temperatures.Earlier models that retrieve higher Bond albedo estimates at values of roughly 0.4, as reported by Schwartz & Cowan (2015) and Schwartz et al. (2017), imply peculiarly high reflectivity at infrared wavelengths.In this case, reflected light needs to be properly accounted for in the phase A24, page 9 of 16 curve analysis, which itself becomes more complex, as the absolute amount of reflected light is larger, and therefore small inaccuracies in its treatment lead to more significant deviations from the true value.On the other hand, models resulting in more moderate estimates, like the measurement A B = 0.16 +0.11  −0.10 by Keating et al. (2019), seem to suffer much less from this problem.

Conclusion
We reported the measurement of the geometric albedo of the hot Jupiter HD 189733b in the CHEOPS bandpass (A g = 0.076 ± 0.016) by analysing 13 observations of its occultation.We started by refining the transit parameters using both TESS and CHEOPS observations (see Table 3).Subsequently, we fitted for the occultation depth (L = 24.7 ± 4.5 ppm) and estimated the contribution of thermal emission from the planet to the total occulted flux.Finally, we inferred the geometric albedo in the CHEOPS bandpass.The measured value is consistent with previous measurements of the geometric albedo of this target and other gas giants, which are often found to be 0.1 (Angerhausen et al. 2015;Esteves et al. 2015;Heng & Demory 2013).The achieved precision of the measurement further proves the capability of CHEOPS to detect low-amplitude occultation signals, just as was done for HD 209458b (Brandeker et al. 2022).
We interpreted the measurement with atmospheric models assuming a cloud-free atmosphere at the pressure level where Rayleigh scattering occurs on the dayside, as well as varying Na enhancement of the planet compared to its host star.We compared these models to the CHEOPS detection as well as the HST STIS upper limit (Evans et al. 2013).Both measurements are consistent with a super-stellar Na elemental abundance (X ≈ 3).In this case the geometric albedo in the TESS band would be very low, making a secondary eclipse of HD 189733b in TESS data unlikely to be detectable.Although consistent within 1σ confidence intervals, CHEOPS measurements hint at HD 189733b having a slightly lower geometric albedo in the CHEOPS band than HD 209548b (Brandeker et al. 2022).This is matched, in the framework of the applied models, by a higher required Na abundance as the planetary brightness at optical wavelengths is keenly sensitive to absorption in the broad Na doublet.Finally, following Schwartz & Cowan (2015), we present lower and upper limits at 3σ confidence on the Bond albedo of the planet at 0.013 and 0.42, respectively.When compared with Bond albedos derived from thermal phase curve observations, we note that earlier higher estimates of the Bond albedo would also require peculiarly high reflectance at infrared wavelengths, while the more recent lower estimates lead to more consistent values.
neighbourhood stars: Delgado Mena et al. (2021) for C, Suárez-Andrés et al. (2016) for N, and Bertran de Lis et al. (2015) for O.We determined the mean value and standard deviation of the abundances of these elements for stars with metallicities similar to that of HD 189733 ([Fe/H] = -0.07± 0.02 dex).Expanding the metallicity range by factors of 5 ([Fe/H] = -0.07± 0.10 dex) and 10 ([Fe/H] = -0.07± 0.20 dex) has a minor impact on the mean abundance and its standard deviation.Finally, to transform the relative atmospheric abundances (see Table

Fig. 2 .
Fig. 2. Top panel: phase-folded and detrended transit light curve from two CHEOPS visits.The light and dark blue points are the original data and the 3-min binned data points, respectively.The black line is the median fitted model and the orange curves are the random models computed from the posterior samples.The residuals of the fit are shown below the light curve.Bottom panel: Same, but using TESS data observed during Sector 41.

Fig. 3 .
Fig. 3. Close-up of the phase-folded, detrended, and fitted occultations of HD 189733b observed by CHEOPS (top) and the residuals of the fitted model (bottom).The light blue points represent the individual data points, the dark blue points the 5 min bins, and the red line the model.

Fig. 4 .
Fig.4.Geometric albedo of hot Jupiters measured by instruments sensitive in the mid-optical, i.e.Kepler (as reported byAngerhausen et al.  2015, blue)  and CHEOPS(Brandeker et al. 2022 and this work, orange)   vs effective temperature.The green line at A g = 0.12 approximates the 'dark roaster' regime with very low optical geometric albedos(Sudarsky et al. 2000).Kepler-7b is a prominent outlier, even when taking into account a recent update on its albedo (alt,Heng et al. 2021).Its unusually bright albedo may be explained by reflection of clouds at the planetary limbs(Adams et al. 2022).

Fig. 6 .
Fig. 6.Contribution of reflected light to the total planetary flux for different grey geometric albedos as a function of wavelength.For the stellar emission spectrum a PHOENIX model (Husser et al. 2013) with the parameters T eff = 5000 K, [M/H] = 0.0, and log g = 4.5 was used.To model the planetary atmosphere an irradiated atmosphere model provided by Mollière et al. (2015) with the parameters T eff = 1250 K, [Fe/H] = 0.0, [C/O] = 0.70, planetary log g = 3.0, and stellar spectral type K5 was adopted.The thermal emission from the planet was extrapolated from the dayside temperature of 1192 ± 9 K, estimated with occultation measurements in the Spitzer 4.5 µm channel (Knutson et al. 2012).

Table 1 .
Overview of the log details of CHEOPS observations.
1 http://github.com/alphapsa/PIPENotes.Time notation follows the ISO-8601 convention.The file keys can be used to retrieve data from the CHEOPS archive.

Table 2 .
Adopted stellar parameters of HD 189733.

Table 3 .
Retrieved planetary and stellar parameters.
-Notes.The Gaussian priors with mean µ and standard deviation σ are displayed as N(µ, σ).U(a, b) shows the uniform prior between a and b.
Evans et al. (2013)red and theoretical geometric albedos of HD 189733 b.The blue dotted line denotes the upper limit for the geometric albedo A g = 0.12 as measured byEvans et al. (2013)in the HST STIS bandpass.The red shaded region denotes our CHEOPS albedo measurement (±1σ).The blue dashed, red, and green curves respectively denote for HST STIS, CHEOPS, and TESS the theoretically calculated albedos A g integrated over the respective bandpasses for different X Na abundances.HST STIS and CHEOPS measurements are consistent with each other if X Na ≥ 6 × 10 −6 .The corresponding TESS albedo would be very low, and thus undetectable according to the model used here.