CHEOPS finds KELT-1b darker than expected in visible light: Discrepancy between the CHEOPS and TESS eclipse depths

Recent TESS-based studies have suggested that the dayside of KELT-1b, a strongly-irradiated brown dwarf, is significantly brighter in visible light than what would be expected based on Spitzer observations in infrared. We observe eight eclipses of KELT-1b with CHEOPS (CHaracterising ExOPlanet Satellite) to measure its dayside brightness temperature in the bluest passband observed so far, and model the CHEOPS photometry jointly with the existing optical and NIR photometry from TESS, LBT, CFHT, and Spitzer. Our modelling leads to a self-consistent dayside spectrum for KELT-1b covering the CHEOPS, TESS, H , Ks, and Spitzer IRAC 3.6 and 4.5 ${\mu}$m bands, where our TESS, H , Ks, and Spitzer band estimates largely agree with the previous studies, but we discover a strong discrepancy between the CHEOPS and TESS bands. The CHEOPS observations yield a higher photometric precision than the TESS observations, but do not show a significant eclipse signal, while a deep eclipse is detected in the TESS band. The derived TESS geometric albedo of $0.36^{+0.12}_{-0.13}$ is difficult to reconcile with a CHEOPS geometric albedo that is consistent with zero because the two passbands have considerable overlap. Variability in cloud cover caused by the transport of transient nightside clouds to the dayside could provide an explanation for reconciling the TESS and CHEOPS geometric albedos, but this hypothesis needs to be tested by future observations.

KELT-1b's emission spectrum 2 has been studied using ground-and space-based eclipse and phase curve observations covering wavelengths from the red visible (TESS) to 4.5 µm (Spitzer).Soon after KELT-1b's discovery, Beatty et al. (2014) observed one secondary eclipse using Spitzer (in IRAC 3.6 and 4.5 µm passbands) and four eclipses in z ′ band using the 0.6 m RCOS telescope at Moore Observatory.Later, Croll et al. (2015) observed an eclipse in the Ks band using the WIRCam on the Canada-France-Hawaii Telescope, Beatty et al. (2017) observed an eclipse in the H band using the LUCI1 spectrograph on the Large Binocular Telescope, and Beatty et al. (2019Beatty et al. ( ) observed a A&A 668, A93 (2022) ) KELT-1b phase curve with Spitzer.Finally, the Transiting Exoplanet Survey Satellite (TESS) observed KELT-1b continuously for 25 days in 2019, and the photometry spanning 15 orbits was studied by Beatty et al. (2020) andvon Essen et al. (2021).
The previous eclipse and phase curve studies have shown that KELT-1b's eclipse spectrum in NIR and mid-infrared agrees with a constant dayside brightness temperature of ∼2900 K (Beatty et al. 2019).However, a recent study of the TESS observations by Beatty et al. (2020) has presented that the brown dwarf's dayside brightness temperature in the TESS band is significantly higher than in the Spitzer bands.They propose that this higher-than-expected brightness in visible light could be due to a high albedo caused by silicate clouds forming on the nightside of the brown dwarf.KELT-1b's dayside is too hot for clouds to form, but the nightside temperature is expected to be cool enough for the formation of silicate clouds.These nightside clouds could then be blown over to the dayside by winds, where they could survive all the way to the local noon, boosting the otherwise low dayside albedo significantly.
However, the cloud hypothesis is not the only viable explanation for the unexpected dayside brightness in the TESS band.von Essen et al. (2021) made an independent analysis of the TESS phase curve and proposed that KELT-1b's dayside emission spectrum can be explained with thermal emission alone.von Essen et al. (2021) suggest that the differences in the TESS and Spitzer brightness temperatures can be explained by collisioninduced absorption due to H 2 -H 2 and H 2 -He decreasing the brightness temperature in the Spitzer bands.So, rather than the brightness of the TESS passband being boosted by reflection from clouds, the Spitzer band brightness would be suppressed by molecular absorption.
CHaracterising ExOPlanet Satellite (CHEOPS; Benz et al. 2021) is an ESA S-class mission aiming to characterise transiting exoplanets around bright (V < 12) stars based on ultra-high quality photometry (Broeg et al. 2013).The satellite is equipped with a 30 cm telescope that allows it to reach a photometric precision of 20 ppm in 6 h of integration time for a V = 9 star (as a design requirement); this is sufficient to detect Earth-sized planets around G5 dwarf stars (Benz et al. 2021).CHEOPS was launched on 18 December 2019, the science observations started in April 2020, and the observations have already been used to improve the characterisation of new and known exoplanet systems (Bonfanti et al. 2021;Delrez et al. 2021;Szabó et al. 2021;Barros et al. 2022;Lacedelli et al. 2022;Wilson et al. 2022), improve the characterisation of eclipsing binaries with M-dwarf companions (Swayne et al. 2021), search for transit timing variations (TTVs, Borsato et al. 2021), search for planets transiting white dwarfs (Morris et al. 2021b), and study secondary eclipses and full phase curves of super-Earths and hot Jupiters (Lendl et al. 2020;Morris et al. 2021a;Hooton et al. 2022;Deline et al. 2022).By now, the satellite has shown it can surpass its design requirements by reaching 10-20 ppm photometric precision for 1 h of integration (Lendl et al. 2020).
In this paper, we present the phase curve analysis of KELT-1b based on new CHEOPS secondary eclipse photometry and existing photometry from TESS, Spitzer, LUCI1, and WIRCam.While the CHEOPS and TESS passbands have a significant overlap, the CHEOPS band emphasises bluer wavelengths than the TESS band.This difference allows us to test the two hypotheses about the TESS-Spitzer dayside brightness temperature discrepancy.If the high-albedo hypothesis is correct, both the CHEOPS and TESS bands are expected to be dominated by reflected light, and the CHEOPS band can be expected to also show a high albedo (albeit this cannot be taken for granted since there are scenarios where KELT-1b could have a low albedo in the CHEOPS band and high in the TESS band, as will be discussed later).If the IR absorption hypothesis is correct, CHEOPS passband brightness temperature could be expected to follow closely the model presented in Fig. 9 in von Essen et al. (2021).
All of the data and code for the analysis and figures are publicly available from GitHub3 .The repository also includes additional Jupyter notebooks that work as appendices to the paper.

CHEOPS
We observed eight eclipses of  with the CHEOPS spacecraft between 2020-10-02 and 2020-12-01 as a part of the Guaranteed Time Observers (GTO) programme (Fig. 1).Each eclipse observation consisted a CHEOPS visit of ≈7 h of near-continuous observations with an exposure time of 60 s.The length of the visits combined with the high precision of the KELT-1b ephemeris (thanks to recent TESS observations) ensured that we obtained at least 2 h of pre-and post-eclipse baseline for each visit.As CHEOPS is in a low-Earth orbit, sections of the observations are unobtainable as they may occur during Earth-target occultations, passages through the South Atlantic Anomaly (SAA), or at times when the stray light level passes an acceptable threshold.These effects manifest as gaps in the light curves and result in lower efficiency observations.For the data presented here, we obtain average efficiency of 56%.
The data were automatically processed using the latest version of the CHEOPS Data Reduction Pipeline (DRP v13;Hoyer et al. 2020) that performs image calibration including bias, gain, non-linearity, dark current, and flat fielding corrections, and amends of environmental and instrumental effects, for example, cosmic-ray hits, smearing trails, and background variations.Subsequently, aperture photometry is conducted on the corrected images using a set of standard apertures; R = 22.5 ′′ (RINF), 25.0 ′′ (DEFAULT), and 30.0 ′′ (RSUP), with an additional aperture that selects the radius based on contamination level and instrumental noise (OPTIMAL), which for these observations yields 20.5 ′′ .Furthermore, the DRP calculates contamination estimates of the visits that comes from background sources, as detailed in Sect.6.1 of Hoyer et al. (2020), which we remove from the data.
As the field of view of CHEOPS rotates due to the nadirlocked orbit of the spacecraft, there may be short-term, nonastrophysical flux trends on the order of a CHEOPS orbital timescale due to nearby contaminants, background variations, or changes in the instrumental environment.In previous studies, these trends have been corrected via linear decorrelation with instrumental basis vectors (Bonfanti et al. 2021;Delrez et al. 2021) or Gaussian process regression (Lendl et al. 2020).However, a recent study has identified that these roll angle trends can be corrected by a novel PSF detrending method (Wilson et al. 2022).To summarise, this tool determines PSF shape changes of a visit by computing vectors associated with these variations by conducting principal component analyses (PCA) of the auto-correlation functions of the CHEOPS subarray images.The components that contribute most to the PSF shape changes are subsequently selected by a leave-one-out-cross-validation and are used as covariates in the light curve model.For this study, we  perform this method on all visits using the DEFAULT aperture photometry fluxes as these yield the lowest RMS (root mean square) noises.

TESS
The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) observed KELT-1b during Sector 17 for 25 days (2019-10-08-2019-11-02) with a 2-min cadence.The observations cover 15 orbital periods with a high-enough precision to measure all the significant phase curve components in the TESS passband.
We use the Presearch Data Conditioning (PDC-SAP) light curve (Stumpe et al. 2014(Stumpe et al. , 2012;;Smith et al. 2012) produced by the SPOC pipeline (Jenkins et al. 2016).As in Beatty et al. (2020), we remove all the exposures with a non-zero quality flag, but we do not apply any additional detrending or outlier rejection.Instead, we use a CELERITE-based (Foreman-Mackey et al. 2017) time-dependent Gaussian process (GP) to calculate the likelihood as explained in more detail later.Furthermore, we ignore the residual light contamination of ≈0.82% reported by von Essen et al. ( 2021), since it will not have any practical effect on the results.

Spitzer
KELT-1 was observed with the Spitzer space telescope using the Infrared Array Camera (IRAC; Fazio et al. 2004) in both S1 (3.6 µm) and S2 (4.5 µm) passbands across multiple programs with the observations covering a secondary eclipse (Beatty et al. 2017) and a full phase curve (Beatty et al. 2019).
Using these data, we extract light curves for all visits using our own custom pipeline.Using the flux-calibrated and artefactcorrected CBCD (Corrected Basic Calibrated Data) frames we constructed light curves by firstly using the IRAC archive contributed IDL program BOX_CENTROIDER.PRO 4 to locate the target in the images.In addition to utilising these values as decorrelation basis vectors during the joint fitting, these x-and y-centroid positions were used as inputs into the APER.PRO code to perform aperture photometry using an aperture radius of 3 pixels and a sky annulus 12-20 pixels.Lastly, we perform the aperture corrections, correct intrapixel gain using the IRACPC_PMAP_CORR.PRO 5 tool, and reject data with fluxes or x-and y-centroid positions 5σ away from the respective, normalised mean.

Stellar characterisation
We revise the stellar radius, mass, and age and show the updated stellar properties in Table 1.The stellar age and mass are compatible with the previous estimates by Siverd et al. (2012) within 1σ, but our new radius is slightly larger than the original estimate of 1.46 +0.039 −0.030 R ⊙ .Using the stellar parameters derived via our spectral analysis as priors, we use a Markov chain Monte Carlo (MCMC) modified infrared flux method (IRFM;Blackwell & Shallis 1977;Schanche et al. 2020) to compute the stellar radius of KELT-1.We construct spectral energy distributions (SEDs) using the ATLAS Catalogue's stellar atmospheric models (Castelli & Kurucz 2003) and subsequently compute synthetic fluxes by integrating the SEDs over passbands of interest with attenuation of the SED included as a free parameter to account for extinction.To calculate the stellar apparent bolometric flux, hence the stellar angular diameter and effective temperature, we compare the synthetic fluxes to the observed broadband photometry retrieved from the most recent data releases for the following passbands; Gaia G, G BP , and G RP , 2MASS J, H, and K, and WISE W1 and W2 (Gaia Collaboration 2021; Skrutskie et al. 2006;Wright et al. 2010).We convert the angular diameter into stellar radius using the offset-corrected Gaia EDR3 parallax (Lindegren et al. 2021) and determine a R ⋆ = 1.530 ± 0.009 R ⊙ .
We then use T eff , [Fe/H], and R ⋆ to compute the stellar mass M ⋆ and age t ⋆ through two different sets of theoretical evolutionary models.In detail, we retrieved a first pair of mass and age estimates by applying the isochrone placement technique (Bonfanti et al. 2015(Bonfanti et al. , 2016) ) to grids of stellar tracks and isochrones pre-computed by the PARSEC8 v1.2S code (Marigo et al. 2017).A second pair of mass and age estimates, instead, was derived by directly inputting T eff , [Fe/H], and R ⋆ into the CLES9 code (Scuflaire et al. 2008), which computes the best-fit evolutionary track following the Levenberg-Marquadt minimisation scheme, as presented in Salmon et al. (2021).As described in Bonfanti et al. (2021), once the two pairs of outcomes are available, we checked their mutual consistency through a χ 2based criterion and finally merged the posterior distributions obtaining M ⋆ = 1.370 +0.047 −0.070 M ⊙ and t ⋆ = 1.8 ± 0.5 Gyr.

Phase curve model
We model the light curves using a phase curve model implemented in PYTRANSIT v2.5.21 (Parviainen 2015(Parviainen , 2020;;Parviainen & Korth 2020) that includes a primary transit, secondary eclipse, thermal emission, reflection, Doppler beaming, and ellipsoidal variation, as illustrated in Fig. 2. The observed flux is presented as a sum of stellar and planetary components divided by a baseline flux model as where f b represents the baseline flux, T is the transit model, f EV is the ellipsoidal variation component, f DB is the Doppler beaming component, E is the eclipse model, k is the planet-star radius ratio, f R is the reflected flux, and f E is the thermal emission10 .The model parametrisation is shown in Table 2 and the parameter priors are discussed in more detail in Sect.5.2.We assume a circular orbit since KELT-1b's orbital eccentricity agrees with zero (Siverd et al. 2012), which allows us to model the phase curve simply as a function of the orbital phase, ϕ = 2π(t − t 0 )/p, where t is the mid-exposure time, t 0 is the mid-transit time, and p is the orbital period.
The light curves are modelled jointly, and each light curve has its own baseline and noise model.We use a simple constant baseline model (that is, a normalisation factor), and model the systematics as a Gaussian process, as detailed later.All the baseline and noise model parameters are free (or loosely constrained) in the posterior optimisation and sampling.

Transits and eclipses
The transits are modelled using the quadratic limb darkening model by Mandel & Agol (2002) as implemented in PYTRANSIT (Parviainen 2015), and the limb darkening is parameterised using the triangular parameterisation by Kipping (2013).The stands for the ratio of the radiation emitted towards the observer by the brown dwarf and the star per projected unit area.That is, the flux ratio is not the ratio of the disk-integrated fluxes, but this can be obtained by multiplying the flux ratio by the companion-star area ratio, k 2 .secondary eclipses are also modelled using PYTRANSIT, using an eclipse model that returns the fraction of the visible planetary disk area to the stellar disk area.We do not include light travel time into the eclipse modelling since it is at maximum ≈25 s.

Reflected light
The reflected light 11 is approximated using a Lambertian phase function (Russell 1916;Madhusudhan & Burrows 2012).This 11 We fix the reflected light component to zero for the main analyses due to its degeneracy with the thermal emission, but describe the component here for consistency.
Notes.The global parameters are independent of passband or light curve, the passband-dependent parameters are repeated for each passband, and the light-curve-dependent parameters are repeated for each separate light curve.N(µ, σ) stands for a normal prior with a mean µ and standard deviation σ, and U(a, b) stands for a uniform distribution from a to b. (a) The zero epoch and orbital period priors are based on the values reported in Siverd et al. (2012) with uncertainties multiplied by three. (b) The limb darkening coefficients have normal priors calculated using LDTK. (c) Ellipsoidal variation is constrained based by EV amplitude ratios as explained in the text and listed in Table 3. (d) The Doppler beaming priors are given in Table 3. (e) The baseline levels have wide normal priors that are based on the light curve variability. (f ) The average log white noise parameters have normal priors centred around a numerical white noise estimate. (g) The Gaussian process log input scales have wide normal priors while the log output scales have normal priors based on light curve variability.
choice is based on simplicity since our data is not precise enough to justify a more realistic (and flexible) phase function.The planet-to-star flux ratio from reflected light is where a g is the geometric albedo, a s the scaled semi-major axis, a s = a/R ⋆ , and α is the phase angle, α = |ϕ − π|.The reflected light reaches its maximum near the eclipse and goes to zero near the transit.

Thermal emission
The thermal emission from the companion is modelled as a sine wave between the minimum flux ratio, f n , and maximum flux ratio, f d , as where o is a phase offset that sets the location of the hot spot, moving the minimum and maximum emission phases marked in Fig. 2.

Ellipsoidal variation
When a star is orbited by a massive companion on a short-period orbit, the companion's gravitational pull distorts the star's shape from a sphere into an ellipse.This causes an ellipsoidal variation (EV) signal (Morris 1985;Faigler & Mazeh 2011;Barclay et al. 2012;Lillo-Box et al. 2014) that can be approximated to a first order as where o is an angular offset, the amplitude, A EV , depends on the companion-host mass ratio, orbital inclination, i, and scaled semi-major axis, a s , as where β is a factor that depends on the linear limb darkening coefficient, u, and gravity darkening coefficient, g, as The EV signal is usually inconsequential when studying planetary phase curves, but it can dominate over the other phase curve components when a star is orbited by a massive companion on a short-period orbit, such as in the case of KELT-1b.Including the EV component into an eclipse and phase curve model is especially important because the EV signal has its minimum close to the mid-eclipse time (Bell et al. 2019;Beatty et al. 2020, andSect. 6.3 in von Essen et al. 2021).Unaccounted for EV signal can artificially enhance the eclipse depth estimated from secondary eclipse observations, as happened with some of the previous KELT-1b observations.

Doppler beaming
Doppler beaming is calculated following Loeb & Gaudi (2003), Barclay et al. (2012), and (Claret et al. 2020) as where the Doppler beaming amplitude is and c is the speed of light in vacuum, G is the gravitation constant, p is the orbital period, and B is the photon-weighted passband-integrated beaming factor (Bloemen et al. 2010).The passband-integrated beaming factor can be calculated from a stellar spectrum for any given passband as where T (λ) is the passband transmission, F λ is the stellar flux at wavelength λ, and B = 5 + d log F λ /d log λ is the beaming factor.Unlike ellipsoidal variations, Doppler beaming has only a minimal effect on the observed eclipse depth.This is because the beaming signal behaves as a linear slope at the vicinity of the eclipse.In addition, beaming is expected to have a much lower amplitude for KELT-1b than EV.

Overview
The goal of our study is to estimate the day-and nightside flux ratios between KELT-1b and its host star in the CHEOPS, TESS, H, Ks, and Spitzer 3.6 and 4.5 µm passbands.We do this by modelling the light curves described in Sect. 2 jointly using the phase curve model defined in Sect.4, which will also yield improved orbital and geometric parameter estimates as a by-product.The analysis follows the standard procedures of Bayesian inference (Gelman et al. 2013;Parviainen 2018), where we aim to estimate the posterior probability distributions (posteriors) for the model parameters given a model, observations, and prior probability distributions (priors) for the model parameters.
We carry out three analyses that are each further divided into separate scenarios (summarised in Appendix B) with different priors on the phase curve model parameters.The analyses are: 1. External data analysis models the TESS, LBT, CFHT, and Spitzer observations jointly, and is implemented in ExternalDataLPF12 Python class.The external data analysis was carried out to create priors for the CHEOPS analysis, and also allows us to test how including the CHEOPS observations affects the final parameter estimates.
2. CHEOPS analysis models only the CHEOPS-observed eclipses and is implemented in the CHEOPSLPF13 class.The main motivation for the CHEOPS analysis is to create detrended CHEOPS light curves to be used in the final analysis.As mentioned in Sect.2.1, the CHEOPS light curves contain strong instrument-related systematics that are captured as basis vectors (covariates) using the approach described in Wilson et al. (2022).Modelling the CHEOPS observations jointly using a linear baseline model leads to a model with 180 free parameters, most of which are baseline model coefficients.This is a lot considering that the external data analysis has 79 free parameters, most of which are physically interesting.After some tests, we decided to carry out the CHEOPS analysis with a full baseline model and simplify the final analysis by using the CHEOPS analysis posterior baseline model to remove the systematics from the CHEOPS light curves included in the final analysis.
3. Final analysis models all the data jointly and is implemented in the FinalLPF14 class.
The posterior estimates from the final analysis are adopted as the final results, but the differences between the analyses and scenarios are discussed in Appendix B.
All the analyses are implemented as Python classes that inherit pytransit.lpf.PhaseCurveLPF, where the base class implements the functionality to model phase curves jointly for a set of heterogeneous light curves15 (including optimisation and For the final and external dataset analyses, we set uninformative or only slightly informative priors on the geometric and orbital parameters, such as the orbital period, impact parameter, and companion-star radius ratio, as listed in Table 2.We also assume that the radius ratio is constant from passband to passband because of a small atmospheric scale height caused by KELT-1b's high surface gravity (Beatty et al. 2020).Further, we assume a circular orbit and force eccentricity to zero, since both the RV and eclipse observations agree with a circular orbit (Siverd et al. 2012).
The transit model adopts a quadratic limb darkening law with the Kipping (2013) parametrisation.We set uninformative priors on the TESS band limb darkening coefficients since the TESS light curve does not feature any strong systematics and the transit signal-to-noise ratio is high.However, we set LDTK-calculated normal priors on the limb darkening coefficients for the two Spitzer passbands because the strong systematics in the Spitzer light curves will not allow us to constrain limb darkening well.

Doppler beaming
Since we have good estimates for the orbital parameters and the brown dwarf to star mass ratio (from Siverd et al. 2012 RV observations), we can set strict priors on the Doppler beaming amplitudes for each passband.We calculate the beaming factors and beaming amplitude priors for each passband with PYTRANSIT (namely using the doppler_beaming_factor and doppler_beaming_amplitude functions), using the BT-Settl (Allard 2013) model spectra to derive the beaming factors rather than a black body approximation.The beaming amplitude estimates are listed in Table 3, and they are computed in the Doppler beaming notebook 16 .

Ellipsoidal variation
As with Doppler beaming, we can calculate theoretical EV amplitudes for all the passbands by combining our prior knowledge about KELT-1b and its host star with stellar spectrum models.We use the companion-star mass ratio estimate from model supersampling and others not), and each possibly with their own noise and systematics models. 16https://github.com/hpparvi/cheops_kelt_1/tree/master/A2_doppler_beaming.ipynb Siverd et al. (2012) and the semi-major axis estimate from Beatty et al. (2020).We estimate the linear limb darkening coefficients for each passband using LDTK (Parviainen & Aigrain 2015), and estimate the gravity darkening coefficients for the CHEOPS passband from the tables in Claret et al. (2020), for the TESS passband from the tables in Claret (2017), and for the H, Ks, and Spitzer passbands from the tables in Claret & Bloemen (2011)  17 .As shown in Table 3, the EV signal amplitude is 10-20 times larger than the DB signal amplitude, being roughly equal to the amplitude of the brown dwarf's phase curve signal in the TESS passband (Beatty et al. 2020;von Essen et al. 2021).
The EV amplitude scales with the semi-major axis as a −3 s , which makes the theoretical EV amplitudes sensitive on our prior a s choice.This causes a potential problem because EV directly affects our eclipse depth measurements.A biased prior a s estimate could also significantly bias the eclipse depths measured for passbands with only near-eclipse photometry (CHEOPS, H, Ks) if we were to base the EV amplitude priors on the theoretical amplitudes.Thus, it would be better to try to find a way to constrain the EV amplitudes in a way that is independent of prior mass ratio and semi-major axis estimates.
Fortunately, the TESS light curve covers the whole orbital phase with a precision that allows us to estimate the EV amplitude in the TESS passband.The EV amplitude ratio between a passband i and the TESS passband t is simply where the β are the factors in Eq. ( 6) that depend on stellar gravity and limb darkening (which both are passband-dependent), but have no dependence on KELT-1b's mass or orbit.Thus, we can constrain the EV signal by constraining the EV amplitudes relative to TESS passband rather than constraining the absolute EV amplitudes themselves.This approach is less likely to introduce biases to EV amplitudes, and so also less likely to bias the secondary eclipse signal.In the final analysis, we set uninformative priors on the absolute EV amplitudes on all passbands, and set informative priors on the EV amplitude ratios between the TESS passband any every other passband where the β factors are calculated using the LDTK and gravity darkening tabulations described earlier and σ is the uncertainty in the β ratio estimate.Further, instead of using the uncertainties listed in Table 3, we set σ to 0.01 to reduce our sensitivity to any biases in the β factors (the chosen value is rather arbitrary but in general 3-10 times larger than the theoretical uncertainties). 17The EV amplitudes are calculated in the EV notebook (https://github.com/hpparvi/cheops_kelt_1/tree/master/A1_ellipsoidal_variation.ipynb) available from GitHub.

Noise and systematics model
In addition to the planetary signal and the EV and Doppler boosting signals from the star, the photometry contains systematics from astrophysical and instrumental sources, as well as from changes in the observing conditions.These systematics need to be taken into account during the phase curve modelling.
We model the variability in the photometry not explained by the phase curve model as a Gaussian process (Rasmussen & Williams 2006;Gibson et al. 2012;Roberts et al. 2013).We use time as the only GP input for the TESS data, which allows us to use CELERITE (Foreman-Mackey et al. 2017) to evaluate the GP log likelihood for the whole unbinned TESS data set.The other data sets use more complex covariance functions with multiple input variables (e.g.x-and y-centroids, airmass and seeing), and for those we evaluate the GP log likelihood using GEORGE (Ambikasaran et al. 2015).
Each GP adds a set of hyperparameters into the joint model.We use simple covariance kernels (Matern-3/2 and squared exponential) which add only a log input scale parameter per input variable and a log output scale parameter per GP (detailed later separately for each data set).These GP hyperparameters are constrained only loosely during the optimisation and MCMC sampling.We standardise the GP input variables (that is, we divide each input variable time series with its standard deviation and remove its mean) and constrain the log input scale parameters using normal priors centred at one with a standard deviation of one.The priors for the log output scale parameters are also normal and centred at log variance of the light curve with a standard deviation of 1.5.
We also fit the average white noise for each light curve.The white noise is parametrised using its log variance, and constrained with loose normal prior centred around a white noise estimate computed from the standard deviation of the flux point-to-point differences as where f is a vector containing the fluxes, sd stands for standard deviation, and diff stands for discrete difference.
The noise models are defined in the _init_lnlikelihood method of the FinalLPF18 and ExternalDataLPF19 analysis classes, and the noise model parameter priors are set in the classes' _post_initialisation method.The posterior MCMC sampling marginalises over the whole GP hyperparameter space allowed by the priors and the data to ensure we are not overfitting the data.
TESS.We model the variability in the TESS data as a Gaussian process with time as the only input variable using the CELERITE package.The GP uses basic Matern-3/2 covariance kernel that yields three additional parameters to the joint model, the log input and output scales and the log white noise variance, and these parameters are constrained with relatively loose normal priors as described earlier.
LBT H.The systematics in the H band light curve are modelled as a GP with eight input variables (including, for example, the airmass and the target and comparison star locations).The covariance kernel is a product of eight squared exponential kernels where each kernel provides an independent input scale parameter into the joint model.This brings the total number of parameters from the H band data set noise model to ten.
CFHT Ks.The systematics in the Ks band light curve are modelled as a GP with two input variables, the x-and y-centroids.The covariance kernel is a product of two squared exponential kernels where each kernel provides an independent input scale parameter into the joint model, leading in total to four additional parameters.
Spitzer.The systematics in the Spitzer data are also modelled using a GP with the x and y centroids as the input variables.The Beatty et al. (2019) Spitzer phase curve observations were split into three 12 h long stares and the telescope was repointed between the stares.This repointing led to strong discontinuities in the raw photometry, so we decided to treat each 12 h stare as a separate light curve.This gives us four light curves per Spitzer passband (one covering the Beatty et al. 2014 eclipse and three covering the phase curve).However, all the eight light curves use the same GP kernel and hyperparameters (we assume the centroid-related systematics in both Spitzer passbands have similar input and output scales), so the Spitzer GP model adds only four parameters to the joint model.
The Spitzer light curves consist of 10 000 (eclipse) and 64 000 (phase curve) flux measurements.This is too much for a standard brute-force GP that requires the inversion of the covariance matrix and scales as O(n 3 ).Splitting the phase curve observations helps, but we still needed to bin the Spitzer observations to a 2-min cadence to make a GP noise model computationally feasible.
CHEOPS.The CHEOPS data is modelled differently in the final analysis and in the initial CHEOPS-only analysis.In the CHEOPS analysis, we use a linear model to model the systematics with basis vectors determined by the pipeline of Wilson et al. (2022) used as covariates.The number of basis vectors used per visit varies from 13 to 29, and the total number of basis vectors (and thus free parameters from the baseline model) is 156.
The noise model for the CHEOPS data assumes i.i.d.(independent and identically distributed) normally distributed noise (that is, white noise), and introduces only the logarithm of the average standard deviation of the noise distribution for each visit as an additional free parameter into the full model (that is, we assume the noise does not vary significantly from observation to observation within a single light curve, and marginalise over the average observation uncertainty).
In the final analysis, we detrend the CHEOPS data using the posterior median baseline model from the CHEOPS analysis and give each visit only a free normalisation term.We do this to simplify the final model, since using the full CHEOPS baseline model would add 156 free parameters into the final model.There is a danger that this approach could lead to underestimated uncertainties in some of our quantities of interest, but a comparison of the CHEOPS and final analysis results (see Appendix B) shows that any effects from not marginalising over the baseline model coefficients in the final analysis are insignificant.

Results from the phase curve analysis
We show the dayside and nightside planet-star flux ratio posteriors from the final joint modelling for all passbands in Fig. 3  Notes.The nightside estimates give an upper limit corresponding to the 99th posterior percentiles.The dayside estimates correspond to posterior medians with uncertainties based on the 16th and 84th posterior percentiles (except for CHEOPS, for which we only give the upper limit).The reported dayside brightness temperature values are lower and upper limits corresponding to 2.5th and 97.5th posterior percentiles.
posteriors for all the passbands in Fig. 5.The final flux ratio estimates are reported in Table 4, and the rest of the parameter estimates in Table 5.The results from the external dataset and CHEOPS-only analysis scenarios are presented in Appendix B. Except for the CHEOPS band, the dayside flux ratio posteriors are approximately normal with a well-defined non-zero mode, and for these posteriors we report their median values with uncertainties based on the 16th and 84th central percentiles in Table 4.The CHEOPS band flux ratio posterior has its mode at zero, and we report its posterior median and 99th percentile upper limits.
The nightside flux ratio posteriors all have their modes at zero, and we report only the 99th posterior percentile upper limits.The nightside flux ratio upper limit for the TESS band is 4.6%, while for the Spitzer bands it is ≈12%.
The dayside flux ratios for the TESS and Spitzer passbands (Fig. 3) agree fairly well with the previous studies by Beatty et al. (2019, 2020), and von Essen et al. (2021).The flux ratios for the H and Ks bands are lower than what was reported by Beatty et al. (2017) and Croll et al. (2015), but this is expected since neither of the original studies included an ellipsoidal variation signal in their models.
The dayside flux ratio for the 3.6 µm Spitzer band is higher than the estimate by Beatty et al. (2019).The Spitzer light curves use relatively flexible GP systematics models and all the Spitzer light curves are well-fit visually (Fig. 5).Even then, we cannot rule out poorly modelled systematics as the source of the anomalously high 3.6 µm flux ratio since the pre-eclipse baseline is relatively short for the full-phase Spitzer observations.0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

CHEOPS dayside flux ratio [%]
Posterior probability P > 99% P > 90% P > 50% Curiously, the CHEOPS and TESS bands show a strong discrepancy in their dayside flux ratios.Panel a of Fig. 6 shows the CHEOPS observations with the median posterior model and panel b shows the same for TESS.The combined CHEOPS observations have a noise estimate of ≈50 ppm over a 20-min bin, half of the TESS noise of ≈97 ppm, but we do not detect a significant eclipse signal in the CHEOPS band.Instead, we are only able to set an upper limit of 1.4% for the CHEOPS dayside flux ratio, while in the TESS band we detect a clear eclipse corresponding to a flux ratio of 6 ± 1%.Notes.The estimates correspond to the posterior median with the uncertainties based on the 16th and 84th posterior percentiles. (a) The equilibrium temperature is calculated using the stellar T eff estimate, scaled semi-major axis distribution, heat redistribution factor distributed uniformly between 0.25 and 0.5, and albedo distributed uniformly between 0 and 0.4.

Dayside brightness temperatures
We estimate the day-and nightside brightness temperatures, T b , separately for each passband using the flux ratio estimates from the final joint modelling and BT-Settl stellar spectrum models by Allard (2013).We calculate the T b estimates as averages over three approaches: (a) physical-physical, where both the brown dwarf and its host star spectra are modelled with BT-Settl models; (b) blackbody-physical, where KELT-1b is modelled as a black body and the star uses a BT-Settl model; and (c) blackbody-blackbody, where both KELT-1b and the star are modelled as black bodies.
Figure 7 shows the approach-averaged brightness temperatures together with the best-fitting constant dayside temperature (the T b estimates are also included in Table 4, and the perapproach estimates are shown in Appendix C).The dayside temperatures mostly agree with each other within the uncertainties.In agreement with Beatty et al. (2020), the TESS passband has a somewhat higher T b than the Spitzer bands, but the disagreement is less significant for us (for example, the TESS and Spitzer 3.6 µm band 68% central posterior intervals overlap).Our T b estimate for the TESS passband agrees well with the Beatty et al. (2020) result, but our Spitzer T b estimates are higher than what was reported in Beatty et al. (2019).
As expected from the differences in the dayside flux ratios, the dayside T b estimate for CHEOPS is significantly lower than for TESS.

Nightside brightness temperatures
We can measure the nightside brightness temperatures only for the passbands where we have observations over the full orbital phase (TESS and Spitzer), and the estimates are rather poorly constrained.Both Spitzer passbands support cool night sides, with 99th T b posterior percentiles ≈ 1850 K, while in the TESS passband we can only set a 3000 K upper limit.These results support a significant temperature difference between the nightand dayside, and agree with the previous Spitzer and TESS analyses by Beatty et al. (2020) and von Essen et al. (2021).

Atmospheric modelling
For the atmospheric modelling, we follow the same basic approach as done in Lendl et al. (2020).In particular, we employ the self-consistent atmosphere model HELIOS (Malik et al. 2017(Malik et al. , 2019) ) to generate a grid of possible atmospheres for KELT-1b.As discussed in Lendl et al. (2020), in the context of 1D models, the flux of the host star impinging upon the top of the atmosphere (TOA) of KELT-1b is (Cowan & Agol 2011) where F * is the stellar photometric flux, R * the radius of the host star KELT-1, a the orbital distance of KELT-1b, A B its Bond albedo, and ϵ the heat redistribution efficiency (0 ≤ ϵ ≤ 1).In terms of the atmosphere model HELIOS, the term is fully degenerate with respect to A B and ϵ, such that we only vary E * rather than A B and ϵ independently.Since KELT-1b is a brown dwarf, its emission might also have a considerable contribution from an internal energy flux, described by an internal temperature T int .In the case of an isolated object without exterior insolation, T int would be the object's effective temperature.For the construction of the model grid, we  therefore use T int as an additional free parameter.Given the lack of any spectrally resolved data for KELT-1b and the fact that the host star's [Fe/H] value is consistent with 0 (see Table 1), we use solar element abundances for the model grid, for simplicity.The general modelling parameters for KELT-1b (surface gravity, stellar radius, orbital distance) are taken from Tables 1 and 5.
Using HELIOS, we generate a model grid of self-consistent forward atmosphere models for KELT-1b in terms of the free parameters T int and E * .In total we compute 165 different atmospheric model structures for internal temperatures between 0 K and 3000 K and E * values from 0.667 to 0.06.The former value is the highest possible one for E * , with a Bond albedo and heat redistribution ϵ of 0. The latter value represents a case with a high Bond albedo and/or very efficient heat redistribution.
In addition to the atmospheric structure, we also obtain the brown dwarf's emission spectrum for each model in a postprocess procedure.For this postprocess, we exclude any scattering contribution, such that the obtained emission spectra describe only the atmosphere's thermal emission.These spectra are then integrated over the passbands shown in Table 4 F pb = F F λ dλ (14) to calculate the flux F pb inside each passband.Here, F is the corresponding filter transmission function and F λ the calculated emission spectrum of KELT-1b.
Since these F pb are a smooth function of T int and E * , we parameterise them as a two-dimensional, fifth-order polynomial in each passband.This allows us to evaluate the passbandintegrated fluxes F pb (T int , E * ) as a continuous function of T int and E * directly, rather than having to perform an interpolation in two dimensions.
These passband-integrated fluxes are then added to a forward model in the Bayesian retrieval framework Helios-r2 (Kitzmann et al. 2020) to constrain the geometric albedos in the TESS and CHEOPS passbands (A g,T and A g,C ).This forward model calculates the secondary eclipse depths dF se as a function of the free parameters A g , R p /a, R p /R * , T int , and E * in each photometric passband.The stellar photospheric spectrum F * is taken from the grid of PHOENIX stellar atmosphere models (Husser et al. 2013), based on the stellar parameters given in Table 1.
Within the retrieval forward model we assume that only the CHEOPS and TESS passbands have contributions by a geometric albedo, whereas the other ones are dominated by thermal emission rather than reflected light.This is based on the fact that the stellar spectrum peaks at lower wavelengths, where the CHEOPS and TESS passbands are located.The TESS band receives 64% of the stellar flux received by the CHEOPS band, and the H, Ks, 3.6 µm, and 4.5 µm bands receive 7%, 3%, 1%, and 0.5% of the flux received by the CHEOPS band, respectively.The four photometric measurements (H, Ks, Spitzer) in the infrared are dominated by the brown dwarf's own thermal emission rather than by scattering of its host star's incident flux.
The posterior distributions of our retrieval calculations are shown in Fig. 8 while the posterior photometric points are depicted in Fig. 9.The posterior distributions show a strong correlation between the parameters E * and T int .When the retrieval is performed only on the TESS and Spitzer data points, bimodal solutions for these two quantities are obtained (see Appendix D).This bimodality is broken primarily by the H and Ks data points, while the CHEOPS data point further narrows the posterior distribution on E * .The retrieved value of T int = 2670 +190 −150 K implies that the brown dwarf has a very strong contribution from an internal heat flux to the total energy budget of the atmosphere and corroborates the findings of Beatty et al. (2020).As already noted by Beatty et al. (2019), an isolated brown dwarf with the same age and mass would have an internal temperature of about 850 K.The relatively low value of E * = 0.22 ± 0.08, on the other hand, suggests that the heat recirculation ϵ must be very high, or that the atmosphere has a large Bond albedo A B .
For the TESS passband, we obtain a geometric albedo of about 0.36 +0.12  −0.13 , while in the CHEOPS passband the results are consistent with a geometric albedo of 0.
The resulting photometry posteriors in comparison to the measured data are shown in the lower panel of Fig.  Posterior model spectrum for the retrieval calculations using Helios-r2.The light grey line shows the high-resolution version of the posterior thermal emission model median, and the darker grey line shows the same spectrum in medium resolution.The black points show the retrieval eclipse depths with their corresponding 1σ intervals, and the orange shaded areas show the eclipse depth posteriors estimated from the photometry as given in Table 4.The black labelled lines show the passbands, where the Spitzer passbands are simplified for visualisation.
results imply that atmosphere forward models can explain most of the measured data points while especially the Spitzer 3.6 µm band shows larger deviations.This mismatch in the Spitzer 3.6 µm band might point to non-equilibrium chemistry effects.
The predicted CHEOPS photometry secondary eclipse depth also deviates from the measured one.In fact, even without including an additional geometric albedo contribution, the measured eclipse depth cannot be described by the expected thermal emission of the planet.The upper panel of Fig. 10 depicts the heat recirculation efficiencies ϵ and Bond albedos A B that are consistent with the retrieved value of E * = 0.22 ± 0.08 and the upper nightside T b limit of 1900 K.If the heat redistribution is very efficient (ϵ ≈ 1), then the Bond albedo has to be close to 0. On the other hand, if only very little heat gets transported to the nightside (ϵ ≈ 0), then the Bond albedo has to be in the range of 0.7, implying that potentially clouds scatter almost 70% of the incident stellar radiation back to space.In principle, given constraints on T int , the dayside temperature and the nightside temperature, the degeneracy between A B and ϵ could be broken.
Additionally, we calculate a separate HELIOS model for the above stated retrieved median values of E * and T int .The resulting temperature-pressure profile is shown in the bottom panel of Fig. 10.Due to the very high T int , the temperature increases very strongly in the lower atmosphere, while being close to isothermal around 0.1 bar, and showing a temperature inversion in the upper atmosphere, originating from the absorption of stellar light.The plot also depicts a set of stability curves for selected high-temperature condensates.These stability curves in combination with the atmosphere's temperature profile imply that most likely condensates do not form on the dayside of KELT-1b.This is not surprising given its high equilibrium temperature.
The derived geometric albedo of 0.36 +0.12 −0.13 in the TESS passband is challenging to explain since even the most refractory mineral (corundum) does not condense at the derived dayside temperatures.However, it agrees with the scenario where silicate clouds form on the nightside of the brown dwarf and are then transported to the dayside by winds, as well as with geometric albedos for silicate cloud-covered planets of ∼0.4 predicted by Sudarsky et al. (2000).It is nevertheless difficult to reconcile this high value of the TESS geometric albedo with a CHEOPS geometric albedo of essentially zero, because the two passbands have considerable overlap.

Chromatic albedo variability as a possible source of discrepancy between CHEOPS and TESS
The albedo spectrum of a planet (or, as in our case, a brown dwarf) can feature clear distinct regions of low and high reflectivity in wavelength space (Sudarsky et al. 2000;Burrows et al. 2008).In theory, a low blue-optical albedo combined with a high NIR albedo could be explained with silicate clouds and a strong optical absorber above the silicate cloud deck (such as gaseous TiO/VO or S 2 /HS, Schwartz & Cowan 2015), and there is a possibility that the discrepancy between the CHEOPS and TESS dayside flux ratios could be explained by a high-contrast feature in KELT-1b's albedo spectrum.
We modelled whether the CHEOPS-TESS discrepancy could be explained by chromatic variations in KELT-1b's albedo using a simple two-level albedo spectrum model parametrised by the geometric albedo in blue, A g,0 , geometric albedo in red, A g,1 , and step location, λ 0 , dictating where the blue albedo changes to red.The modelling is detailed in Appendix A and carried out in the chromatic albedo variation test notebook 20 .
As a conclusion, the discrepancy is challenging to explain with a simple two-level albedo spectrum model due to the large overlap between the CHEOPS and TESS bands.The largest TESS to CHEOPS reflected light contrast, C = (F T − F C )/(F T + F C ), of 0.41 is obtained for a step model with A g = 0.01 for λ < 860 nm and A g = 0.4 for λ ≥ 860 nm, as illustrated in Fig. A.1.A geometric albedo of 0.01 is not unheard of for a highly irradiated body, but is similar to what has been estimated for TrES-2b (Kipping & Spiegel 2011), while the geometric albedo of 0.4 corresponds to the silicate cloud models by Sudarsky et al. (2000), and is in line with the previous KELT-1b study by Beatty et al. (2020).However, the maximum contrast of 0.4 is at the lower tail of the measured C posterior distribution that has its mode at 0.9 (Fig. A.2).Even a blue A g of 0.001 with λ 0 = 940 nm would lead to C of only 0.6, so it is unlikely that a physically plausible albedo spectrum would be able to reproduce the observed CHEOPS-TESS discrepancy.Pressure (bar) Fig. 10.Upper panel: values of E * as a function of the heat recirculation parameter ϵ and Bond albedo A B .The solid white line represents the retrieved median value of E * , the dashed lines are its corresponding 1σ interval, and the shaded region corresponds to the parameter space excluded by the upper nightside T b limit of 1900 K. Lower panel: temperature-pressure profile of a HELIOS model (solid line), computed for the median values of E * and T int .Besides the temperature-pressure profile, the lower panel also includes a set of condensate stability curves (dashed lines) for a selected set of species.

Temporal variability as a possible source of discrepancy between CHEOPS and TESS
Since the upper limits on the nightside temperature are low enough that clouds could conceivably form, it is plausible that nightside clouds could be transported to the dayside (Beatty et al. 2020), where they could survive transiently.This raises the possibility of dayside cloud cover being transient.Variability in cloud cover could be one possible explanation for the dayside flux ratio discrepancy between the CHEOPS and TESS passbands.von Essen et al. ( 2021) reported possible temporal variation in the eclipse depths in the TESS data but concluded that the variations were most likely associated with stellar variability.We likewise carried out per-eclipse analyses for the TESS and CHEOPS data, but also need to conclude that the per-eclipse precision is not high enough to study whether the eclipse depths show systematic variations that would securely not be caused by stellar variability or systematics.
The spatial distribution of transient dayside clouds transported from the brown dwarf's nightside can be expected to be strongly nonuniform, which would lead to an asymmetric reflection phase curve as with Kepler-7b (Demory et al. 2013).While our analyses simplify the KELT-1b's flux contribution by using only the emission component, the emission component has a loosely constrained phase offset that allows it to catch any existing phase curve asymmetries.
The emission offset (see Table 5) agrees with zero for the TESS and Spitzer 4.5 µm bands, but is significantly non-zero for the Spitzer 3.6 µm band.Of these three, only the TESS band phase curve is expected to have a significant contribution from reflection, and the lack of measurable phase curve asymmetry speaks against the cloud transport hypothesis.Then again, were the dayside cloud coverage to be variable, the variability could blur the asymmetry when averaged over the TESS observations of 25 days.
The TESS and Spitzer phase offsets disagree with the ones estimated by Beatty et al. (2019Beatty et al. ( , 2020)).They observe slight eastward shifts for the TESS and Spitzer 4.5 µm bands and a strong eastward shift for the Spitzer 3.6 µm band, while our shifts are westwards.
The phase curve asymmetry in Spitzer 3.6 µm band is curious.The Spitzer phase curve observations were carried out eight days apart of each other, so the differences could be explained by temporal variability in emission and reflection.This could also explain the higher-than expected dayside eclipse depth for the 3.6 µm band (Fig. 9).Our eclipse depth estimate differs from the result from Beatty et al. (2019; whose estimate agrees with the theoretical models) and we extracted the Spitzer light curves using our own pipeline, so we cannot completely rule out untreated systematics.However, the timescale of the variability (smooth over KELT-1b phase) combined with our use of a flexible GP to model the Spitzer systematics lead us to believe the differences in the two Spitzer phase curves arise from true variability in the KELT-1b's phase curve.

Comparison between emission spectrum models
We compare our retrieved emission spectrum with the spectra by von Essen et al. ( 2021), the BT model spectra (Allard 2013) 21and DRIFT-PHOENIX spectra (Witte et al. 2011) 22 for isolated brown dwarfs in Fig. 11.The "Helios-r2" results correspond to the modelling described earlier in Sect.6.4 and the "VE21a" and "VE21b" results correspond to the modelling detailed in von Essen et al. (2021).The rest of the results correspond to best-fitting BT-SETTL, BT-NextGen, BT-COND, BT-Dusty, and DRIFT-PHOENIX models where the host star was modelled using a BT-SETTL spectrum with T Eff = 6500 K, log g = 4.0 and z = 0.0.
The best fitting BT and DRIFT-PHOENIX models correspond to a dayside temperature of 3000 K (DRIFT-PHOENIX) and Fig. 11.Comparison between different emission spectrum models.The shaded areas show the 68%, 95%, and 99.7% central posterior intervals for the eclipse depths as estimated from the joint modelling, and the black points show the passband-integrated emission spectrum models.Helios-r2 corresponds to the model presented here, VE21a and VE21b to the models in von Essen et al. ( 2021), BT-CN, BT-DS, BT-NG, and BT-ST to the best-fitting BT-Cond, BT-Dusty, BT-NextGen, and BT-Settl models by Allard (2013), respectively, and DRIFT-PH to the best-fitting Drift-PHOENIX model (Witte et al. 2011).3100 K (BT models) and fit even the CHEOPS passband surprisingly well.However, none of the isolated brown dwarf models can simultaneously explain both the CHEOPS and TESS eclipse depths.

Sensitivity on parametrisation and priors
The choice of parametrisation and priors generally affects the parameter posteriors, and this effect can be significant when studying weak signals.We parameterise the day-and nightside emission using log flux ratios and set uniform priors on these, which equals to setting log-uniform (reciprocal) priors on the day-and nightside emission flux ratios.This differs from previous studies using uniform priors on the flux ratios.
We repeated our final joint analysis parametrising it with flux ratios rather than log flux ratios to study how sensitive our posteriors are on the choice of parametrisation.The most notable difference was in nightside flux ratios, where uniform priors led to distributions with non-zero modes and the maximum flux ratios (99th posterior percentile) were somewhat higher than when using log-uniform priors.The second notable difference was in the CHEOPS band dayside flux ratio posterior, which again had a non-zero mode with slightly larger maximum values.However, the difference in the maximum flux ratio was not large enough to change the main outcome of the analysis.The dayside flux ratio posteriors for other passbands were not significantly affected.
We adopt the log flux ratio parametrisation since this should be more suitable for a "scale" parameter with an unknown magnitude.This ensures that we do not give too much weight on apparent posterior modes if in reality we can only estimate the upper limit securely.

The role of ellipsoidal variation in the CHEOPS band
Our CHEOPS band eclipse depth estimate is very sensitive on the amplitude of the ellipsoidal variation signal, as illustrated in Fig. B.3.However, the two approaches leading to physically expected EV amplitudes yield similar CHEOPS band eclipse depths, and only nonphysically low EV amplitudes can lead to eclipse depths that would agree with the TESS measurement.

Conclusions
We have estimated a self-consistent dayside emission spectrum of KELT-1b covering the CHEOPS, TESS, H, Ks, and Spitzer IRAC 3.6 and 4.5 µm passbands by modelling new CHEOPS secondary eclipse photometry of KELT-1b jointly with the existing ground-and space-based light curves.The study adds so far the bluest point to KELT-1b's emission spectrum, improves the two ground-based NIR measurements by using a phase curve model that includes a stellar ellipsoidal variation signal, and improves the Spitzer measurements by modelling the Beatty et al. (2017Beatty et al. ( , 2019) ) observations jointly with all the other light curves.
The emission spectrum largely agrees with the previous studies.The H and Ks dayside flux ratios are lower than the prior estimates but this was expected since the previous modelling ignored the stellar ellipsoidal variation signal.The TESS and Spitzer 4.5 µm bands also agree with the results by Beatty et al. (2019, 2020), and von Essen et al. (2021).The 3.6 µm Spitzer band dayside flux ratio estimate deviates both from previous results by Beatty et al. (2019) and theoretical expectations.This difference can either be due to insufficiently modelled systematics or non-equilibrium chemistry effects in the A93, page 15 of 23 A&A 668, A93 (2022) brown dwarf atmosphere, but we cannot securely determine its cause.
Another, significantly more difficult-to-explain discrepancy occurs between the CHEOPS and TESS passbands: the TESS observations show a deep secondary eclipse with a depth of 360 ± 60 ppm, but the CHEOPS observations strictly rule out an eclipse with a depth larger than 80 ppm (Fig. 6).The CHEOPS observations yield nearly twice higher photometric precision than the TESS observations, so an eclipse similar to the TESS band one would be easily detected.
Atmospheric models can reproduce KELT-1b's emission spectrum fairly well in most of the passbands considered, but generally manage to manage to fit only either the CHEOPS or TESS band well.Atmospheric modelling with HELIOS leads to a solution where KELT-1b's geometric albedo is ≈0.3 in the TESS passband but consistent with 0 in the CHEOPS band.A contrast like this in the geometric albedo is difficult to explain due to the overlap between the two passbands.The models by von Essen et al. ( 2021) explain the TESS band without reflection but cannot explain the CHEOPS band eclipse depth (Fig. 11).Finally, the BT and DRIFT-PHOENIX models (Allard 2013;Witte et al. 2011) for isolated brown dwarfs also reproduce the observed KELT-1b emission spectrum fairly well, but have trouble reconciling the CHEOPS and TESS eclipse depths.
It is challenging to explain the discrepancy between the CHEOPS and TESS band eclipse depths.Temporal variability in cloud cover caused by the transport of transient nightside clouds to the dayside (as suggested by Beatty et al. 2020) could provide a potential explanation.High contrast features in KELT-1b's albedo spectrum could also play a role in explaining a part of the discrepancy, and could be caused by a layer of strong optical absorber such as gaseous TiO/VO or S 2 /HS residing above a silicate cloud layer (Schwartz & Cowan 2015), but are unable to explain the discrepancy fully.
All in all, instead of shedding light on which of the previously proposed theories might work best to explain KELT-1b's atmosphere, our additional blue-optical eclipse depth measurement has introduced a new open question that is challenging to explain with our current knowledge.Further observations combined with comprehensive modelling are required to study whether KELT-1b's dayside brightness spectrum is explained by transient silicate clouds (Beatty et al. 2020), thermal emission (von Essen et al. 2021), or something else, and how the discrepancy between the CHEOPS and TESS band eclipse depths can be explained.
Chromatic albedo variation test notebook.

Appendix B: Comparison between EV scenarios
In addition to the final joint analysis, we carried out a set of analyses for the external dataset (ED) consisting of the TESS, H, Ks, and Spitzer photometry, and another set of analyses for the CHEOPS photometry alone.The ED analysis was carried out to provide priors on the orbital and geometric parameters for the CHEOPS analysis and to study how different approaches to constrain the EV amplitudes affect our parameter estimates.
The CHEOPS analysis was carried out to also study how the EV amplitude constraints affect the eclipse depth estimate, and to provide detrended CHEOPS light curves for the final joint analysis.
Both the ED and CHEOPS-only analyses consider three scenarios that differ in the priors set on the EV amplitude: a) Theoretical EV: the EV amplitudes are given normal priors based on theoretical estimates in Table 3.This is the most constraining scenario but can lead to biased eclipse depth estimates if the semi-major axis estimate used to calculate the amplitudes is biased.b) Constrained EV ratios: the EV amplitudes are constrained relative to the EV amplitude in the TESS passband and the TESS passband EV amplitude is given an uninformative prior.This scenario works as a stepping stone between the strongly constrained and completely unconstrained EV cases and is not sensitive to our prior semi-major axis estimate.c) Unconstrained EV : all the EV amplitudes are given uninformative priors to see whether the EV amplitudes from the TESS and Spitzer passbands (where we have coverage over the full orbital phase) agree with the theoretical expectations.
The CHEOPS-only analysis considers also an additional scenario with EV amplitude forced to zero: d) No EV: the EV amplitude is forced to zero.The motivation for this scenario is to see how ignoring EV would affect our CHEOPS eclipse depth estimate.
Except for the data included, the models are identical to the full joint model.The baseline and noise models for the ED analyses are the same as in the final analysis.However, the baseline in the CHEOPS-only analysis is modelled differently than in the final analysis.In the CHEOPS analysis, we use the basis vectors determined by the pipeline by Wilson et al. as covariates and each visit is given its own set of baseline coefficients.The number of basis vectors per visit varies from 13 to 29, and the total number of basis vectors (and thus free parameters from the baseline model) is 156.This approach allows us to detrend the CHEOPS light curves together with the phase curve model to be used in the final joint analysis (this is safe and does not affect the final parameter estimates).The CHEOPS analysis uses the posteriors from the external dataset analysis as priors for the zero epoch, orbital period, stellar density, impact parameter, and planet-star area ratio.
We show the EV amplitude posteriors for all the analyses in

Appendix C: Comparison between brightness temperature calculation approaches
We calculate the brightness temperatures discussed in Sect.6.2 using three approaches: physical-physical, where both the star and the brown dwarf are modelled using the BT-Settl spectra by Allard (2013); blackbody-physical, where KELT-1b is modelled as a black body and its host star using a BT-Settl spectrum; and blackbody-blackbody, where both the brown dwarf and its host star are modelled as black bodies.Our reported brightness temperature values correspond to averages over the three approaches, but we also tested whether the different approaches agree with each other or not.
We present the day-and nightside brightness temperatures separately for each approach in  When the H and Ks data points are added to the retrieval, the bimodality is broken (not shown).This can be understood by looking at the ranges of T int and E * values that are consistent with the measurements in the passbands depicted in D.1.The results presented in this figure suggest that the H and Ks secondary eclipse measurements are only consistent with high values of T int Finally, in the second row of Fig. D.2 we perform a retrieval using all available measurements but with the Spitzer secondary eclipse depths from Beatty et al. (2019).The addition of the CHEOPS data point now further narrows the posterior on E * down to a median value of 0.2.this behaviour can also be easily understood by comparison with the theoretically calculated eclipse depths in Fig. D.1.For CHEOPS, the results indicate that the measured eclipse depth is only compatible with a maximum E * of about 0.2.
We also note that the results using the Spitzer data from Beatty et al. (2019) essentially provide the same posterior distribution as using the reduction of the Spitzer of this work (see Fig. 8). .In particular, the very high geometric albedo in the TESS passband is consistently obtained with either dataset.Despite the difference in the reported Spitzer 3.6 µm eclipse depth, the retrieval outcome is, thus, unaffected because the solution with T int ≈ 2700 K has a blackbody peak at shorter wavelengths, within the H and Ks passbands.A93, page 23 of 23

Fig. 1 .
Fig.1.Detrended photometry from eight KELT-1b eclipses observed with CHEOPS centred around the eclipse centre time (T e ).Each eclipse is observed as a single continuous CHEOPS visit of ≈7 h (labelled as V1-V8 in the figure), but the satellite's orbital configuration leads to periodic gaps in the photometry, as described in Sect.2.1.
Beatty et al. (2017) observed a spectrally resolved eclipse of KELT-1b in the H band using the LUCI1 multiobject spectrograph on the 2 × 8.4 m Large Binocular Telescope (LBT).The observations were carried out on the night of 2013-10-26 and lasted ≈5 h.The exposure time was 60 s, leading to 245 exposures in total.6   2.5.Canada-France-Hawaii Telescope Ks bandCroll et al. (2015) observed an eclipse of KELT-1b in the Ks band (∼2.15 µm) using the Wide-field Infrared Camera (WIRCam) on the 3.6 m Canada-France-Hawaii Telescope (CFHT).The observations were carried out in 2012-10-10 and covered 6.8 h centred around the expected eclipse centre with an exposure time of 8 s, leading to 1440 exposures in total7 .The reduction of the photometry is detailed extensively in Croll et al. (2015, see especially the discussion about the sensitivity to different data reduction approaches in Sect. 5

Fig. 2 .
Fig. 2. Components of a phase curve model (transit, eclipse, emitted light, reflected light, Doppler beaming, and ellipsoidal variation) and the full model.The figure is centred around the secondary eclipse, the horizontal dotted line shows the baseline level (flux = 1), the vertical dashed line shows the eclipse centre, and the solid grey line shows the transit centre.

Fig. 3 .
Fig. 3.The dayside (orange) and nightside (blue) planet-star flux ratio posterior distributions for all passbands.The central 68% posterior interval is marked with a dark shade and the previous dayside flux ratio estimates are shown as diamonds with errorbars showing their reported uncertainties.For TESS, the Beatty et al. (2020) estimate is shown on the left and the von Essen et al. (2021) estimate on the right.

Fig. 4 .
Fig. 4. Dayside flux ratio posterior distribution for the CHEOPS passband.The dashed lines shows the posterior percentile limits.

Fig. 5 .
Fig. 5. Observed light curves from CHEOPS, TESS, LBT, CFHT, and Spitzer (points) and the posterior median phase curve model (black line).The systematics model has been removed from the observations, the TESS and CHEOPS observations have been folded over the phase, and the data has been binned to 10 min for visualisation purposes.

Fig. 6 .
Fig. 6.CHEOPS eclipse observations binned to a 20-min cadence with the phase curve model (panel a) and TESS observations with a similar binning (panel b).The median posterior phase curve model is shown as a black line and its 68% central posterior interval is shown as a shaded blue region.

Fig. 7 .Fig. 8 .
Fig. 7. Dayside and nightside brightness temperatures estimated from the flux ratios using the BT-Settl stellar model spectra by Allard (2013).The shading marks the 68% and 99.7% central posterior intervals, from the darkest to the lightest shade, and the dashed horizontal line shows the best-fitting constant brightness temperature.The previous dayside T b estimates for the TESS and Spitzer passbands by Beatty et al. (2019, 2020), and von Essen et al. (2021) are shown as points with errorbars.

40 A
Fig. A.1.Top: a model spectrum of KELT-1 (light blue line), the CHEOPS and TESS passband transmission functions (black lines).Middle: the step-function albedo spectrum model.Bottom: TESS to CHEOPS reflected light contrast as a function of A g,0 and λ 0 .The step locations corresponding to the maximum flux ratios are marked as dashed vertical lines.
Fig. B.1.Ellipsoidal variation amplitude posteriors for f) final joint analysis with relative prior on the EV amplitude, a) theoretical prior on EV amplitude, b) relative prior on EV amplitude, and c) uninformative prior on EV amplitude.Case c with an uninformative prior on EV is shown only for the TESS and Spitzer passbands with full phase coverage.
Fig. C.1.All three approaches agree with each other within uncertainties.
Figure D.1 shows the degeneracy between T int and E * when the model accounts for the TESS secondary eclipse.Within this family of solutions, the solution of von Essen et al. (2021), which explains the dayside emission spectrum without the need for reflected light, is included.In Fig.D.2, we perform a retrieval with just the TESS, Spitzer 3.6 µm and 4.5 µm data points reported byBeatty et al. (2019).The left panel of D.2 shows a a retrieval with the full Fig. D.1.Eclipse depths in ppm as a function of the internal temperature T int and the parameter E * based on the HELIOS model calculations described in Sect.6.4.The six panels show the eclipse depths in the various passbands.White, solid lines depict the measured occultation depths from Table4, dashed lines refer to their 1σ error bars.The red lines for the two Spitzer passbands refer to the eclipse depths reported byBeatty et al. (2019).

Table 2 .
Phase curve model parameters and priors used in the final joint analysis.

Table 3 .
Theoretical Doppler beaming and ellipsoidal variation amplitudes, and ellipsoidal variation amplitude ratios relative to the TESS passband.

Table 4 .
Final day-and nightside planet-star flux ratios and eclipse depths with their uncertainties and upper limits.

Table 5 .
Final estimates for the geometric and orbital parameters of KELT-1b and its host star.