© H. Parviainen et al. 2022

1 Introduction

KELT-1b is a strongly irradiated transiting brown dwarf orbiting a 6500 K F5-star on a short-period orbit of 1.2 days as discovered by Siverd et al. (2012). It has the third-shortest orbital period of the known transiting brown dwarfs¹ (Carmichael et al. 2021) after TOI-263b (P = 0.56 day, Parviainen et al. 2020; Palle et al. 2021) and NGTS-7Ab (P = 0.7 day, Jackman et al. 2019), but receives a significantly higher stellar insolation since TOI-263b and NGTS-7Ab both orbit low-mass M dwarfs with T_eff ~ 3400 K. The orbit and insolation level of KELT-1b is similar to ultra-hot Jupiters (Parmentier et al. 2018; Kitzmann et al. 2018; Lothringer et al. 2018) such as MASCARA-2b (Lund et al. 2017; Talens et al. 2018), WASP-121b (Delrez et al. 2016), WASP-189b (Anderson et al. 2018), WASP-76b (West et al. 2016), and WASP-33b (Cameron et al. 2010), but its surface gravity is ~22 times larger than that of Jupiter. This makes KELT-1b an interesting case to study how surface gravity affects the atmospheres of strongly irradiated sub-stellar objects.

KELT-1b’s emission spectrum² 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. (2019) observed a KELT-lb 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) and von 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 collision-induced absorption due to H₂–H₂ and H₂–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 GitHub³. The repository also includes additional Jupyter notebooks that work as appendices to the paper.

2 Observations

2.1 CHEOPS

We observed eight eclipses of KELT-1b 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 nadir-locked orbit of the spacecraft, there may be short-term, non-astrophysical 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.

Fig. 1 Detrended photometry from eight KELT-1b eclipses observed with CHEOPS centred around the eclipse centre time (Te). 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.

2.2 TESS

The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) observed KELT-lb 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, 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.

2.3 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 artefact-corrected CBCD (Corrected Basic Calibrated Data) frames we constructed light curves by firstly using the IRAC archive contributed IDL program BOX_CENTROIDER.PRO⁴ 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⁵ tool, and reject data with fluxes or x- and y-centroid positions 5σ away from the respective, normalised mean.

2.4 Large Binocular Telescope H band

Beatty et al. (2017) observed a spectrally resolved eclipse of KELT-lb 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.⁶

2.5 Canada-France-Hawaii Telescope Ks band

Croll et al. (2015) observed an eclipse of KELT- lb 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 total⁷. 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.4).

3 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 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, 2016) to grids of stellar tracks and isochrones pre-computed by the PARSEC⁸ 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 CLES⁹ 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 χ² -based criterion and finally merged the posterior distributions obtaining M_⊙ and t_* = 1.8 ± 0.5 Gyr.

4 Theory

4.1 Phase curve model

We model the light curves using a phase curve model implemented in PYTRANSIT v2.5.21 (Parviainen 2015, 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 (1)

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 emission¹⁰.

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₀)/p, where t is the mid-exposure time, t₀ 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.

4.2 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 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.

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.

4.3 Reflected light

The reflected light¹¹ is approximated using a Lambertian phase function (Russell 1916; Madhusudhan & Burrows 2012). This 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 (2)

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.

4.4 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 (3)

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.

4.5 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 (4)

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 (5)

where β is a factor that depends on the linear limb darkening coefficient, u, and gravity (darkening coefficient, g, as (6)

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, and Sect. 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.

4.6 Doppler beaming

Doppler beaming is calculated following Loeb & Gaudi (2003), Barclay et al. (2012), and (Claret et al. 2020) as (7)

where the Doppler beaming amplitude is (8)

and c is the speed of light in vacuum, G is the gravitation constant, p is the orbital period, and 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 (9)

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.

5 Phase curve analysis

5.1 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 ExternalDataLPF¹² 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 CHEOPSLPF¹³ 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 FinalLPF¹⁴ 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 curves¹⁵ (including optimisation and MCMC sampling), and the derived classes mainly read in the different light curves, set up the noise models, and set up the model parameter priors based on the analysis scenario.

5.2 Priors

5.2.1 Orbit, geometry, and limb darkening

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.

5.2.2 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¹⁶.

5.2.3 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 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)¹⁷. 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 , 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 (10)

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).

5.3 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 (11)

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 FinalLPF¹⁸ and ExternalDataLPF¹⁹ 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³). 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.

6 Results and discussion

6.1 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, give a more detailed view of the CHEOPS band dayside flux ratio posterior in Fig. 4, and show the observations and the model 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.

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 day-side flux ratio, while in the TESS band we detect a clear eclipse corresponding to a flux ratio of 6 ± 1%.

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 Dayside flux ratio posterior distribution for the CHEOPS passband. The dashed lines shows the posterior percentile limits.

6.2 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 Tb 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.

6.3 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 night-and dayside, and agree with the previous Spitzer and TESS analyses by Beatty et al. (2020) and von Essen et al. (2021).

6.4 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, 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) (12)

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 (13)

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 (14)

to calculate the flux F_pb inside each passband. Here, ℱ 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 passband-integrated 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 (15)

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 Κ 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 Κ. 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 , 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. 8. The 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 Κ. 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 snowing 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-lb. This is not surprising given its high equilibrium temperature.

The derived geometric albedo of in the TESS pass-band 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.

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 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 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 Tb, 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.

Fig. 8 Posterior distributions for the retrieval calculations using Helios-r2. We show only the four most important parameters. The rest are strongly constrained by priors and do not show significant correlations between other parameters.

Fig. 9 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.

6.5 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₂/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-lb’s albedo spectrum.

We modelled whether the CHEOPS-TESS discrepancy could be explained by chromatic variations in KELT-lb’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, λ₀, 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²⁰.

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 λ₀ = 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.

Fig. 10 Upper panel: values of E, as a function of the heat recirculation parameter ϵ and Bond albedo AB. 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 Tb limit of 1900 K. Lower panel: temperature-pressure profile of a HELIOS model (solid line), computed for the median values of E, and Tint. Besides the temperature-pressure profile, the lower panel also includes a set of condensate stability curves (dashed lines) for a selected set of species.

6.6 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. (2019, 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.

6.7 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)²¹ and DRIFT-PHOENIX spectra (Witte et al. 2011)²² 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 3100 Κ (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.

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).

6.8 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 night-side 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.

6.9 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.

7 Conclusions

We have estimated a self-consistent dayside emission spectrum of KELT-lb 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. (2017, 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 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 S2/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.

Acknowledgements

We thank the anonymous referee for their insight and useful comments. We kindly thank T. G. Beatty and B. Croll for providing their observations to our use, and C. von Essen, M. Mallonn, N. Madhusudhan, and A. Piette for providing their model spectra from von Essen et al. (2021). CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. A.C.C. and T.G.W. acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. M.L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project FOUR ACES. grant agreement no. 724427). It has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). D.E. acknowledges financial support from the Swiss National Science Foundation for project 200021_200726. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalizacão by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987, O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. L.M.S. and D.G. gratefully acknowledge financial support from the CRT foundation under Grant No. 2018.2323 “Gaseous or rocky? Unveiling the nature of small worlds”. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). Y.A. and M.J.H. acknowledge the support of the Swiss National Fund under grant 200020_172746. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiología (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. X.B., S.C., D.G., M.F. and J.L. acknowledge their role as ESA-appointed CHEOPS science team members. A.Br. was supported by the SNSA. A.C.C. acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L.D. is an FRS-FNRS Postdoctoral Researcher. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-190080). M.F. and C.M.P. gratefully acknowledges the support of the Swedish National Space Agency (DNR 65/19, 174/18). M.G. is an FRS-FNRS Senior Research Associate. S.H. gratefully acknowledges CNES funding through the grant 837319. K.G.I. is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. This work was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. P.M. acknowledges support from STFC research grant number ST/M001040/1. G.Sc., G.Pi., I.Pa., L.Bo., V.Na. and R.Ra. acknowledge the funding support from Italian Space Agency (ASI) regulated by “Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C”. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). I.R.I. acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme. Gy.M.Sz acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a PRODEX Institute Agreement between the ELTE Eötvös Loránd University and the European Space Agency (ESA-D/SCI-LE-2021-0025), the Lendület LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. V.V.G. is an FRS-FNRS Research Associate. N.A.W. acknowledges UKSA grant ST/R004838/1. M.F. gratefully acknowledges the support of the Swedish National Space Agency (SNSA; DNR 177/19, 174/18, and 65/19). S.S. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 833925, project STAREX).

Appendix A Chromatic albedo variability

Appendix A.1 Albedo spectrum model

We use a simple two-level albedo spectrum model to study whether the discrepancy between the CHEOPS and TESS band dayside flux ratios could be caused by chromatic variability in albedo. The model is parametrised by two geometric albedo levels, A_g,0 (blue) and A_g,i (red), and a step location λ₀, so that (A.1)

as illustrated in the middle panel of Fig. A.1.

Our main quantity of interest is the contrast between the reflected light in the CHEOPS and TESS bands, (A.2)

S is the stellar spectrum (a BT-Settl spectrum by Allard 2013), A_g is the albedo spectrum model, and T is the transmission function for either TESS (F_T) or CHEOPS (F_c).

We assume A_g,1 that matches the geometric albedo for a “roaster” planet with silicate clouds from Sudarsky et al. (2000), compute C for a set of A_g,0 values between 0.01 (similar to TrES-2b, Kipping & Spiegel 2011) and 0.05, and λ₀ values between 400 and 1000 nmn and show the results in the lower panel of Fig. A.1. The maximum C is achieved for λ₀ ~ 800 nm and varies from 0.17 to 0.41. We obtain C = 0.41 for the lowest included blue A_g (as expected), A_g,0 = 0.01, and C decreases quickly with increasing A_g,0.

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 Ag,0 and λ0. The step locations corresponding to the maximum flux ratios are marked as dashed vertical lines.

Appendix A.2 Comparison between modelled and observed reflection ratio

Before we can compare our theoretical C values to the C distribution from the light curve analysis, we need to estimate and remove the contribution from thermal emission in the latter. This contribution should be small for the CHEOPS and TESS bands but is not necessarily negligible.

We calculate a set of “reflection only” flux ratio samples from the original dayside flux ratio samples using a Monte Carlo approach. We draw a set of KELT-1b temperatures from a uniform distribution from 2800 to 3200 Κ (corresponding to the KELT-1b Spitzer 4.5 µm T_b ranges given in Table 4), we calculate the emission contribution in the CHEOPS and TESS bands using BT-Settl model spectra by Allard (2013) for each temperature sample, and finally we remove these contributions from a random sample of CHEOPS and TESS dayside flux ratios. The full process is carried out in the chromatic albedo variation test notebook.²³ Figure A.2 shows the original dayside flux ratios with emission and reflection, and the “reflection-only” flux ratios. The correction is rather small for the CHEOPS band but significant for the TESS band. The bottom panel of Fig. A.2 shows the final observed R distribution.

Our maximum C of 0.41 is at the lower tail of the observed C distribution. The maximum C corresponds to a low but still physically realistic red A_g of 0.01, which is similar to what has been measured for TrES-2b (Kipping & Spiegel 2011). However, while both A_g,0 and A_g,1 are physically plausible separately, it is unclear if this is the case for an albedo spectrum combining the two.

Fig. A.2 Top and middle: original CHEOPS and TESS dayside flux ratio distribution estimates together with versions of the distributions with the thermal emission removed as described in Appendix A.2. Bottom: the final TESS to CHEOPS contrast (C) distribution and the maximum C of ≈ 0.4 from the albedo spectrum modelling is marked with a dashed vertical line.

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 Fig. B.1, the dayside flux ratio posteriors in Fig. B.2, and the CHEOPS dayside posteriors in Fig. B.3. All in all, the dayside flux ratio posteriors all agree with each other for the final joint modelling and EV prior scenarios a and b. The flux ratio posteriors also agree for scenario d for passbands with photometry covering a full orbital phase, but the free- and no-EV scenarios (c and d) lead to overestimated flux ratios for CHEOPS photometry.

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. B.2 As in Fig. B.1, but for the dayside flux ratio.

Fig. B.3 Dayside flux ratios the CHEOPS-only analysis scenarios. The scenario labels are the same as in Fig. B.1 with an additional no-EV scenario d.

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 Fig. C.1. All three approaches agree with each other within uncertainties.

Fig. C.1 Dayside and nightside brightness temperatures estimated from the flux ratios using the p) physical-physical, x) blackbody-physical, and b) blackbody-blackbody approaches. The shading marks the 68% and 99.7% central posterior intervals.

Appendix D Additional retrievals to explore the dependence of outcomes on data

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 by Beatty et al. (2019). The left panel of D.2 shows a a retrieval with the full prior ranges. This yields a bimodal solution for T_int and E_*: a high value of E_* and a low value of T_int, or a low value of E_* and a high value of T_int. The two different solutions are isolated in the other two panels of Fig. D.2.

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.

Fig. D.1 Eclipse depths in ppm as a function of the internal temperature Tint 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 Table 4, dashed lines refer to their 1σ error bars. The red lines for the two Spitzer passbands refer to the eclipse depths reported by Beatty et al. (2019).

Fig. D.2 Top: Posterior distributions for a retrieval with Helios-r2 using only the TESS secondary eclipse measurement and the two Spitzer points from Beatty et al. (2019). The left panel shows the results for the full range of prior distributions discussed in Sect. 6.4, which yields a bimodal solution with respect to Tint and the parameter E*. The two other panels isolate the two solution modes. Bottom: Posterior distributions for a retrieval with Helios-r2 using all secondary eclipse measurement from Table 4 but with the two Spitzer points from Beatty et al. (2019).

References

All Tables

All Figures

	Fig. 1 Detrended photometry from eight KELT-1b eclipses observed with CHEOPS centred around the eclipse centre time (Te). 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.
In the text

	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.
In the text

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.

In the text

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

	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.
In the text

	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.
In the text

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 Tb, 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.

In the text

	Fig. 8 Posterior distributions for the retrieval calculations using Helios-r2. We show only the four most important parameters. The rest are strongly constrained by priors and do not show significant correlations between other parameters.
In the text

Fig. 9 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.

In the text

Fig. 10 Upper panel: values of E, as a function of the heat recirculation parameter ϵ and Bond albedo AB. 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 Tb limit of 1900 K. Lower panel: temperature-pressure profile of a HELIOS model (solid line), computed for the median values of E, and Tint. Besides the temperature-pressure profile, the lower panel also includes a set of condensate stability curves (dashed lines) for a selected set of species.

In the text

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).

In the text

	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 Ag,0 and λ0. The step locations corresponding to the maximum flux ratios are marked as dashed vertical lines.
In the text

	Fig. A.2 Top and middle: original CHEOPS and TESS dayside flux ratio distribution estimates together with versions of the distributions with the thermal emission removed as described in Appendix A.2. Bottom: the final TESS to CHEOPS contrast (C) distribution and the maximum C of ≈ 0.4 from the albedo spectrum modelling is marked with a dashed vertical line.
In the text

	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.
In the text

	Fig. B.2 As in Fig. B.1, but for the dayside flux ratio.
In the text

	Fig. B.3 Dayside flux ratios the CHEOPS-only analysis scenarios. The scenario labels are the same as in Fig. B.1 with an additional no-EV scenario d.
In the text

	Fig. C.1 Dayside and nightside brightness temperatures estimated from the flux ratios using the p) physical-physical, x) blackbody-physical, and b) blackbody-blackbody approaches. The shading marks the 68% and 99.7% central posterior intervals.
In the text

Fig. D.1 Eclipse depths in ppm as a function of the internal temperature Tint 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 Table 4, dashed lines refer to their 1σ error bars. The red lines for the two Spitzer passbands refer to the eclipse depths reported by Beatty et al. (2019).

In the text

Fig. D.2 Top: Posterior distributions for a retrieval with Helios-r2 using only the TESS secondary eclipse measurement and the two Spitzer points from Beatty et al. (2019). The left panel shows the results for the full range of prior distributions discussed in Sect. 6.4, which yields a bimodal solution with respect to Tint and the parameter E*. The two other panels isolate the two solution modes. Bottom: Posterior distributions for a retrieval with Helios-r2 using all secondary eclipse measurement from Table 4 but with the two Spitzer points from Beatty et al. (2019).

In the text