# The stable climate of KELT-9b,

K. Jones, B. M. Morris, B. -O. Demory, K. Heng, M. J. Hooton, N. Billot, D.

Ehrenreich, S. Hoyer, A. E. Simon, M. Lendl, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

K. Jones, B. M. Morris, B. -O. Demory, K. Heng, M. J. Hooton, et al. The stable climate of KELT-9b, *Astronomy and Astrophysics*, 2022, 666, 10.1051/0004-6361/202243823. insu-03839239

HAL Id: insu-03839239

https://insu.hal.science/insu-03839239

Submitted on 4 Nov 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
The stable climate of KELT-9b*,★★

K. Jones¹, B. M. Morris¹, B.-O. Demory¹, K. Heng¹,², M. J. Hooton³,⁴, N. Billot⁵, D. Ehrenreich⁵,⁶, S. Hoyer⁷, A. E. Simon⁴, M. Lendl⁵, O. D. S. Demangeon⁸,⁹, S. G. Sousa⁸, A. Bonfanti¹⁰, T. G. Wilson¹¹, S. Salmon⁵, Sz. Csizmadia¹², H. Parviainen¹³, G. Bruno¹⁴, Y. Alibert⁴, R. Alonso¹³,¹⁵, G. Anglada¹⁶,¹⁷, T. Bárczy¹⁸, D. Barrado¹⁹, S. C. C. Barros⁸,⁹, W. Baumjohann¹⁰, M. Beck⁵, T. Beck⁴, W. Benz⁴,¹, X. Bonfils²⁰, A. Brandeker²¹, C. Broeg⁴,¹, J. Cabrera¹², S. Charnoz²², A. Collier Cameron¹¹, M. B. Davies²³, M. Deleuil⁷, A. Deline⁵, L. Delrez²⁴,²⁵, A. Erikson¹², A. Fortier⁴,¹, L. Fossati¹⁰, M. Fridlund²⁶,²⁷, D. Gandolfi²⁸, M. Gillon²⁴, M. Güdel²⁹, K. G. Isaak³⁰, L. L. Kiss³¹,³², J. Laskar³³, A. Lecavelier des Etangs³⁴, C. Lovis⁵, D. Magrin³⁵, P. F. L. Maxted³⁶, V. Nascimbeni³⁵, G. Olofsson²¹, R. Ottensamer³⁷, I. Pagano¹⁴, E. Pallé¹³, G. Peter³⁸, G. Piotto³⁵,³⁹, D. Pollacco², D. Queloz⁴⁰,³, R. Ragazzoni³⁵,³⁹, N. Rando⁴¹, F. Ratti⁴¹, H. Rauer¹²,⁴²,⁴³, C. Reimers³⁷, I. Ribas¹⁶,¹⁷, N. C. Santos⁸,⁹, G. Scandariato¹⁴, D. Ségransan⁵, A. M. S. Smith¹², M. Steller¹⁰, Gy. M. Szabó⁴⁴,⁴⁵, N. Thomas⁴, S. Udry⁵, V. Van Grootel²⁵, I. Walter³⁸, N. A. Walton⁴⁶, and W. Wang Jungo¹

(Affiliations can be found after the references)

Received 20 April 2022 / Accepted 8 August 2022

ABSTRACT

Even among the most irradiated gas giants, so-called ultra-hot Jupiters, KELT-9b stands out as the hottest planet thus far discovered with a dayside temperature of over 4500 K. At these extreme irradiation levels, we expect an increase in heat redistribution efficiency and a low Bond albedo owed to an extended atmosphere with molecular hydrogen dissociation occurring on the planetary dayside. We present new photometric observations of the KELT-9 system throughout 4 full orbits and 9 separate occultations obtained by the 30 cm space telescope CHEOPS. The CHEOPS bandpass, located at optical wavelengths, captures the peak of the thermal emission spectrum of KELT-9b. In this work we simultaneously analyse CHEOPS phase curves along with public phase curves from TESS and *Spitzer* to infer joint constraints on the phase curve variation, gravity-darkened transits, and occultation depth in three bandpasses, as well as derive 2D temperature maps of the atmosphere at three different depths. We find a day-night heat redistribution efficiency of ~0.3 which confirms expectations of enhanced energy transfer to the planetary nightside due to dissociation and recombination of molecular hydrogen. We also calculate a Bond albedo consistent with zero. We find no evidence of variability of the brightness temperature of the planet, excluding variability greater than 1% (1σ).

**Key words.** techniques: photometric – instrumentation: photometers – planets and satellites: atmospheres – planets and satellites: gaseous planets – eclipses – occultations

**1. Introduction**

Understanding the climate of an exoplanet involves quantifying
its global thermal and chemical structure. To date there have
been a plethora of methods used to aid in this characterisation,
including observing a phase curve. This is a lightcurve of the
system taken as the planet orbits the star. The planet can add
to the total flux of the system with either thermal emission or
reflected light from the star. Similar to phases of the Moon, dif-
ferent fractions of the dayside of the exoplanet will be visible to
us at different times in its orbit. There are three main observ-
ables within a phase curve: phase curve amplitude and shape,
the occultation depth (when the planet goes behind the star)
and the planet's hotspot phase offset (from the substellar point).
From this information, one can derive 2D temperature maps
of the planet's photosphere. For planets with a rotation period
significantly less than its orbital period, it could be assumed
that the planet would absorb the light from the star and then
re-radiate it as thermal emission relatively uniformly across the

entire surface. However, for tidally locked planets, where a sin-
gle hemisphere is always facing the star, it can be seen that the
heat from the star is not so evenly distributed. The level of this
'heat redistribution' can provide insights into the dynamics and
chemistry of the planet's atmosphere.

Ultra-hot Jupiters (UHJs) are a newly emerging class of short-period exoplanets with temperatures exceeding ~2500 K. Crucially, several chemical properties distinguish them from regular hot Jupiters. Firstly, their dayside atmospheres are hot enough that hydrogen is predicted to exist in its atomic form (Bell & Cowan 2018). The photon energies involved in H⁻ opacities become the dominant source of the spectral continuum over Rayleigh scattering (Arcangeli et al. 2018), resulting in a low optical geometric albedo. As these planets are tidally locked with their host star, the hydrogen will recombine to molecular hydrogen at cooler longitudes, assisting with the heat transport around the planet. The prediction of higher heat redistribution in UHJs as opposed to moderately hot Jupiters follows this (Bell & Cowan 2018). Secondly, the high dayside temperatures result in the thermal dissociation of most molecules, leaving only water and carbon monoxide as the main opacity sources (see also, e.g. Lothringer et al. 2018; Kitzmann et al. 2018; Parmentier et al. 2018). Following this, metals are predicted to exist in their

<footnote>* The photometric time series data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/666/A118</footnote>

<footnote>** The CHEOPS program ID is CH_PR100036.</footnote>
atomic form rather than being bound in molecules, a prediction that has been observationally verified for KELT-9b (Hoeijmakers et al. 2018, 2019; Pino et al. 2020; Bello-Arufe et al. 2022). Due to the high levels of irradiation, UHJs are expected to have extended atmospheres which also lends itself to hydrogen atmospheric escape, a phenomenon also detected for KELT-9b (Yan & Henning 2018; Wyttenbach et al. 2020). Chemically, UHJs are objects in between gas-giant exoplanets and stars.

Not only is it important to characterise the thermal and reflective properties of a planet, but also any transient departures from the mean global state (weather). Observed variability may be caused by dynamical processes in the atmosphere of the exoplanet (see review by Heng & Showman 2015); for example, baroclinic instabilities in the Earth's atmosphere occur on timescales of 3–7 days (Peixóto & Oort 1984). To date, quantifying the climate of an exoplanet and any associated transient phenomena has only been attempted for hot Jupiters. This is due to their short orbital periods (which facilitate repeated observations) and large atmospheric pressure scale heights (due to their hot, hydrogen-dominated atmospheres). Armstrong et al. (2016) claimed variability from HAT-P-7b in the form of a shift in the peak offset of its *Kepler* phase curve from either side of the substellar point. However this was reassessed in Lally & Vanderburg (2022) and they concluded that the apparent variations were artifacts most likely caused by stellar supergranulation. Furthermore, Agol et al. (2010) analysed 7 transits and 7 secondary eclipses of HD 189733b, measured by the *Spitzer* Space Telescope, and set an upper limit of 2.7% on the variability of flux from its dayside. In addition, Owens et al. (2021) scrutinised 27 days of TESS data and found no detectable variability from WASP-12b (another UHJ). Other works focussed on the aspects of a climate that can be probed with phase curves include Kreidberg et al. (2018), who used *Hubble-WFC3* and *Spitzer* phase curves of the UHJ WASP-103b to infer inefficient dayside-to-nightside heat redistribution, an absence of water and a carbon-to-oxygen ratio less than 0.9. Arcangeli et al. (2019) used *Hubble-WFC3* phase curves of the UHJ WASP-18b to find that atmospheric drag is needed to explain the observed dayside-nightside flux contrast and they speculated that this drag may be magnetic in nature.

KELT-9b is the hottest known exoplanet and the most extreme member of the UHJ class with an equilibrium temperature of $4050 \pm 180$ K (Gaudi et al. 2017). It has a near-polar, 1.48-day orbit around a rapidly rotating A0/B9 star ($T_{\text{eff}} = 10\,170$ K) (Gaudi et al. 2017). Its dayside brightness temperature has been measured to be 4600 K in both the z-band (Gaudi et al. 2017) and the TESS bandpass (Wong et al. 2020), making the dayside of KELT-9b as hot as the photosphere of a K4 star. This dayside brightness temperature implies a peak blackbody emission wavelength of 0.63 µm, which sits in the centre of the CHEOPS bandpass. In Hooton et al. (2018), they found a near-ultraviolet geometric albedo upper limit of < 0.14 (3$\sigma$) and in Sudarsky et al. (2000) at wavelengths applicable to CHEOPS, they predict a geometric albedo around 0.05, which corresponds to a percentage of reflected light in the CHEOPS bandpass of $\le 8\%$ of the total dayside flux (thermal emission plus reflected light). These properties imply that the CHEOPS space telescope is well-positioned to monitor the thermal flux of KELT-9b.

The Characterising Exoplanet Satellite (CHEOPS) mission is a 30 cm space telescope in a low-Earth orbit since December 2019 (Benz et al. 2021). A major component of the guaranteed-time observations planned in the nominal mission involve a thorough atmospheric classification of a wide range of transiting exoplanets, using precise full-phase curve and

occultation observations. Published detections include a hot dayside atmosphere for the UHJ WASP-189b (Lendl et al. 2020; Deline et al. 2022) variations in the phase curve of the super-Earth 55 Cnc e (Morris et al. 2022), and a hint of dayside reflection for another UHJ, MASCARA-1b (Hooton et al. 2022). In this work, we embark on a comprehensive observational campaign to quantify the climate and variability of KELT-9b. Using the CHEOPS space telescope, we observed 9 secondary eclipses and 4 full phase curves. Additionally, we re-analyse the 4.5 µm *Spitzer* and TESS phase curves of KELT-9b. We interpret the CHEOPS, TESS, and *Spitzer* phase curve variation jointly using *kelp*, a recently published framework that describes a 2D temperature map using parabolic cylinder functions (Morris et al. 2022).

In Sect. 2, we detail the technical aspects of the CHEOPS, TESS, and *Spitzer* observations and include a description of the CHEOPS space telescope and its data reduction pipeline. Section 3 explains the different sections of the data analysis and includes a detailed description of the models used to model the phase curves in each bandpass and the CHEOPS occultations. In Sect. 4, we report the results of our phase curve fit and occultation eclipse depths. Within this section, in Sects. 4.1 and 4.2 we show the bandpass-dependent 2D thermal maps derived from the best-fit phase curve parameters and discuss the resulting day and nightside integrated temperatures. The results of the CHEOPS occultation analysis is detailed in Sect. 4.3 and in Sect. 4.4 we report the Bond albedo and heat redistribution efficiency. Finally in Sect. 4.5, we compare our results with previous results using the same datasets. A discussion of the impact and importance of this work for the study of the climate on UHJs and future multi-wavelength JWST phase curves can be found in Sect. 5.

## 2. Observations

### 2.1. CHEOPS observations

CHEOPS is an on-axis Ritchey-Chrétien telescope (Benz et al. 2021). Its nadir-locked orientation keeps the bottom of the spacecraft pointed towards Earth throughout each orbit, causing the field of view of CHEOPS to rotate during science observations. This frequently results in systematic noise in phase with the spacecraft roll angle, as neighbouring stars contribute varying levels of contamination into the CHEOPS aperture over one spacecraft orbit (e.g.: Lendl et al. 2020). Detailed systematic investigations have revealed that CHEOPS observations occasionally contain a ramp feature whereby the flux of the first few orbits can show a significant correlation with the temperature fluctuations of the telescope assembly and is most likely caused by the change in temperature of the telescope as it adjusts to a new pointing position (Morris et al. 2021).

CHEOPS observed nine occultations of KELT-9 between 25 July 2020 and 24 July 2021 and four phase curves; observed on 31 August 2020, 10 September 2020, 31 July 2021 and 22 August 2021, each observation lasting around 2.4 days. The individual occultation observations lasted between 5.9 and 6.9 h, distributed over 4–5 CHEOPS orbits. Further details of the individual observations can be found in Table 1.

CHEOPS data is automatically processed by the CHEOPS Data Reduction Pipeline (DRP, Hoyer et al. 2020). The CHEOPS DRP performs the basic calibration of the science images (i.e. bias, dark, flat-field corrections) and also performs background correction, cosmic-rays hits removal, correction of bright stars' smear trail contamination and provide estimations of the flux contamination of background stars (see details in
<table><caption>**Table 1.** CHEOPS observation logs, corresponding to the occultation-only observations in the first 9 rows and the phase curve observations in the last 4 rows.</caption><thead><tr><th>Date start<br>(UT)</th><th>Date stop<br>(UT)</th><th>File key</th><th>Duration<br>(dd:hh:mm)</th><th>Exposure<br>Time (s)</th><th>Exposures<br>per stack</th><th>Efficiency<br>(%)</th></tr></thead><tbody><tr><td>2020-07-25 07:15</td><td>2020-07-25 14:05</td><td>CH_PR100036_TG001201_V0200</td><td>00:06:50</td><td>36.7</td><td>3 (×12.2 s)</td><td>61</td></tr><tr><td>2020-08-04 16:02</td><td>2020-08-04 21:52</td><td>CH_PR100036_TG001202_V0200</td><td>00:05:50</td><td>36.7</td><td>3 (×12.2 s)</td><td>66</td></tr><tr><td>2020-08-17 23:54</td><td>2020-08-18 05:49</td><td>CH_PR100036_TG001203_V0200</td><td>00:05:55</td><td>36.7</td><td>3 (×12.2 s)</td><td>73</td></tr><tr><td>2020-08-20 23:28</td><td>2020-08-21 05:57</td><td>CH_PR100036_TG001204_V0200</td><td>00:06:28</td><td>36.7</td><td>3 (×12.2 s)</td><td>65</td></tr><tr><td>2020-08-23 23:20</td><td>2020-08-24 05:21</td><td>CH_PR100036_TG001205_V0200</td><td>00:06:00</td><td>36.7</td><td>3 (×12.2 s)</td><td>72</td></tr><tr><td>2020-08-28 10:32</td><td>2020-08-28 16:30</td><td>CH_PR100036_TG001206_V0200</td><td>00:05:58</td><td>36.7</td><td>3 (×12.2 s)</td><td>71</td></tr><tr><td>2020-09-03 08:19</td><td>2020-09-03 14:17</td><td>CH_PR100036_TG001207_V0200</td><td>00:05:58</td><td>36.7</td><td>3 (×12.2 s)</td><td>72</td></tr><tr><td>2020-09-13 16:15</td><td>2020-09-13 22:44</td><td>CH_PR100036_TG001208_V0200</td><td>00:06:28</td><td>36.7</td><td>3 (×12.2 s)</td><td>50</td></tr><tr><td>2021-07-24 15:45</td><td>2021-07-24 21:37</td><td>CH_PR100036_TG001209_V0200</td><td>00:05:52</td><td>36.7</td><td>3 (×12.2 s)</td><td>69</td></tr><tr><td>2020-08-31 13:29</td><td>2020-09-03 00:53</td><td>CH_PR100036_TG001001_V0200</td><td>02:11:24</td><td>36.7</td><td>3 (×12.2 s)</td><td>61</td></tr><tr><td>2020-09-10 03:31</td><td>2020-09-12 14:11</td><td>CH_PR100036_TG001002_V0200</td><td>02:10:39</td><td>36.7</td><td>3 (×12.2 s)</td><td>59</td></tr><tr><td>2021-07-31 04:25</td><td>2021-08-02 15:06</td><td>CH_PR100036_TG000901_V0200</td><td>02:10:40</td><td>36.7</td><td>3 (×12.2 s)</td><td>63</td></tr><tr><td>2021-08-22 11:04</td><td>2021-08-24 22:15</td><td>CH_PR100036_TG000902_V0200</td><td>02:11:10</td><td>36.7</td><td>3 (×12.2 s)</td><td>62</td></tr></tbody></table>

**Notes.** The File Key is useful for uniquely identifying the visits used in this work.

Hoyer et al. (2020). Finally, the DRP extracts the photometric signal of the target in 3 standard apertures called RINF, DEFAULT and RSUP (at radius of $R=21.5''$, $25''$ and $30''$), in addition to the OPTIMAL aperture with a radius set to minimise the effect of the contamination by close-by background stars. In the case of KELT-9 observations, this aperture was set at $R=40''$. In our analysis we use the light curves obtained with the DEFAULT aperture.

## 2.2. TESS observations

The TESS satellite (Ricker et al. 2014) observed more than 20 phase curves of KELT-9 during Sectors 14 and 15 of the telescope’s operation (months of July and August 2019). These observations were first published in Wong et al. (2020) using their own reduction techniques. For our analysis we used the PDCSAP flux measurements at 2 min cadence (pre-reduced using the TESS SPOC pipeline to remove long-term trends and systematics; Jenkins et al. 2016). We downloaded the data and stitched the light curves into a single array using lightkurve (Lightkurve Collaboration et al. 2018).

## 2.3. Spitzer observations

In this work, we analyse the Spitzer archival data of KELT-9 b that have already been published (Mansfield et al. 2020). We downloaded KELT-9 b archival IRAC data from the Spitzer Heritage Archive<fnref n="1" />. The data consist of one full phase curve at 4.5 µm split in two Astronomical Observation Requests (AORs) obtained under program ID 14059 (PI J. Bean). The reduction and analysis of these datasets are similar to Demory et al. (2016a). We model the correlated noise associated with IRAC intra-pixel sensitivity (Ingalls et al. 2016) using a modified implementation of the BLISS (BiLinearly Interpolated Sub-pixel Sensitivity) mapping algorithm (Stevenson et al. 2012).

In addition to the BLISS mapping, our baseline model includes a linear function of the Point Response Function (PRF) Full Width at Half-Maximum (FWHM) along the *x* and *y* axes, which significantly reduces the level of correlated noise as shown in previous studies (see, e.g.: Lanotte et al. 2014; Demory et al. 2016b,a; Gillon et al. 2017; Mendonça et al. 2018). Our baseline model does not include time-dependent parameters. Our implementation of this baseline model is included in a Markov chain Monte Carlo (MCMC) framework already presented in the literature (Gillon et al. 2012). We run two chains of 200 000 steps each to determine the phase curve properties at 4.5 µm based on

<table><caption>**Table 2.** Stellar parameters used to derive the stellar radius and mass used in this paper.</caption><thead><tr><th>Stellar parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Metallicity</td><td>dex</td><td>0.14 ± 0.30</td></tr><tr><td>Surface gravity</td><td>dex</td><td>4.1 ± 0.3</td></tr><tr><td>Radius</td><td>R<sub>☉</sub></td><td>2.379 ± 0.038</td></tr><tr><td>Mass</td><td>M<sub>☉</sub></td><td>2.45<sup>+0.19</sup><sub>-0.17</sub></td></tr><tr><td>Age</td><td>Gyr</td><td>0.3 ± 0.1</td></tr></tbody></table>

**Notes.** This table also includes these calculated parameter values.

the entire dataset described above. From our BLISS mapping + FWHM baseline model, we obtain a median RMS of 723 ppm per 23 s integration time for that dataset.

## 2.4. Stellar parameters

The stellar parameters used in this analysis are shown in Table 2. As a key stellar prior in our analysis, we determine the stellar radius of KELT-9 using a Markov-chain Monte Carlo (MCMC) modified infrared flux method (IRFM; Blackwell & Shallis 1977; Schanche et al. 2020). To achieve this, we compute synthetic fluxes from constructed spectral energy distributions (SEDs). These were built from stellar atmospheric models and stellar parameters derived via the spectral analysis in Borsa et al. (2019) that were integrated over bandpasses of interest with the SED attenuated to determine the extinction within the model fitting. We compared these fluxes to observed broadband photometry retrieved from the most recent data releases for the following bandpasses; Gaia G, G<sub>BP</sub>, and G<sub>RP</sub>, 2MASS J, H, and K, and WISE W1 and W2 (Gaia Collaboration 2021; Skrutskie et al. 2006; Wright et al. 2010) to calculate the apparent bolometric flux, and hence the stellar angular diameter and effective temperature. In our analysis, we use stellar atmospheric models from the ATLAS Catalogues (Castelli & Kurucz 2003). By converting the angular diameter to the stellar radius using the offset-corrected Gaia EDR3 parallax (Lindegren et al. 2021), we obtain $R_{\star} = 2.379 \pm 0.038 R_{\odot}$.

The stellar radius $R_{\star}$ together with the effective temperature $T_{\text{eff}}$ and the metallicity [Fe/H] constitute the input set to determine the isochronal age $t_{\star}$ and mass $M_{\star}$. To make the derivation process more robust we employed two different stellar evolutionary models, namely PARSEC v1.2S<fnref n="2" /> (Marigo

<footnote n="2">² PAdova and TRieste Stellar Evolutionary Code: http://stev.oapd.inaf.it/cgi-bin/cmd</footnote>

<footnote n="1">¹ http://sha.ipac.caltech.edu</footnote>
<table><caption>**Table 3.** Global priors and best-fit values for the model and detrending parameters as described in Sect. 3, along with the derived parameters.</caption><thead><tr><th>Parameter</th><th></th><th>Unit</th><th>Prior</th><th>Best-fit value</th></tr></thead><tbody><tr><td colspan="5"><strong>Global</strong></td></tr><tr><td colspan="5"><em>Fitted parameters</em></td></tr><tr><td>Period</td><td>P</td><td>days</td><td>N(1.4811235, 0.0000011) <sup>(a)</sup></td><td>1.48111949 ± 0.00000034</td></tr><tr><td>Transit duration</td><td>T<sub>14</sub></td><td>days</td><td>U(0.13, 0.2) <sup>(a)</sup></td><td>0.16552<sup>+0.00016</sup><sub>-0.00015</sub></td></tr><tr><td>Impact parameter</td><td>b</td><td>deg</td><td>N(0.168, 0.017) <sup>(b)</sup></td><td>0.195<sup>+0.016</sup><sub>-0.015</sub></td></tr><tr><td>Sky-projected spin orbit angle</td><td>λ</td><td>deg</td><td>N(-85.78, 0.46) <sup>(b)</sup></td><td>-85.67<sup>+0.46</sup><sub>-0.45</sub></td></tr><tr><td>Gravity-darkening coefficient</td><td>β</td><td></td><td>U(0.01, ∞), N(0.237, 0.01) <sup>(d)</sup></td><td>0.2270<sup>+0.0077</sup><sub>-0.0080</sub></td></tr><tr><td>GP periodic timescale<sup>(*)</sup> (SHO)</td><td>ρ<sub>SHO</sub></td><td>days</td><td>U(0.2, 0.4)</td><td>0.3386<sup>+0.0037</sup><sub>-0.0036</sub></td></tr><tr><td>Stellar inclination</td><td>i<sup>(*)</sup></td><td>deg</td><td>U(0, 180)</td><td>47.1 ± 1.1</td></tr><tr><td colspan="5"><em>Derived parameters</em></td></tr><tr><td>Semi major axis</td><td>a</td><td>R<sub>*</sub></td><td></td><td>3.0914<sup>+0.0090</sup><sub>-0.0100</sub></td></tr><tr><td>Orbital inclination</td><td>i</td><td>deg</td><td></td><td>86.38<sup>+0.29</sup><sub>-0.30</sub></td></tr><tr><td>Stellar density</td><td>ρ</td><td>g cm<sup>-3</sup></td><td>N(0.256, 0.33) <sup>(c)</sup></td><td>0.2548<sup>+0.0022</sup><sub>-0.0025</sub></td></tr><tr><td>GP damping timescale<sup>(*)</sup> (SHO)</td><td>τ<sub>SHO</sub></td><td>days</td><td></td><td>0.6772<sup>+0.0074</sup><sub>-0.0072</sub></td></tr><tr><td>Stellar rotational period</td><td>P<sub>rot</sub></td><td>hrs</td><td>U(4.61, ∞) <sup>(e)</sup></td><td>18.96 ± 0.34</td></tr><tr><td>True spin orbit angle</td><td>Ψ</td><td>deg</td><td></td><td>84.36<sup>+0.37</sup><sub>-0.38</sub></td></tr><tr><td colspan="5"><em>Fixed parameters</em></td></tr><tr><td>Stellar radius</td><td>R<sub>*</sub></td><td>R<sub>☉</sub></td><td>Source</td><td>Value</td></tr><tr><td>Effective/polar stellar temperature</td><td>T<sub>eff</sub></td><td>K</td><td>This work (Sect. 2.4)</td><td>2.379</td></tr><tr><td>Projected rotational velocity</td><td>ν sin i<sup>(*)</sup></td><td>km s<sup>-1</sup></td><td>Gaudi et al. (2017)</td><td>10170</td></tr><tr><td>Eccentricity</td><td>e</td><td></td><td>Borsa et al. (2019)</td><td>111.8</td></tr><tr><td>Argument of periastron</td><td>ω</td><td>deg</td><td>assumed</td><td>0</td></tr><tr><td>Fluid number</td><td>α</td><td></td><td>assumed</td><td>90</td></tr><tr><td>Drag frequency</td><td>ω<sub>drag</sub></td><td></td><td>Morris et al. (2022)</td><td>0.6</td></tr><tr><td>Highest present spherical mode</td><td>l<sub>max</sub></td><td></td><td>Morris et al. (2022)</td><td>4.5</td></tr><tr><td>GP Amplitude<sup>(*)</sup> (SHO)</td><td>σ<sub>SHO</sub></td><td>ppm</td><td>Morris et al. (2022)</td><td>1</td></tr><tr><td>GP Amplitude<sup>(*)</sup> (Matérn)</td><td>σ<sub>Matern</sub></td><td>ppm</td><td></td><td>100</td></tr><tr><td>GP timescale<sup>(*)</sup> (Matérn)</td><td>ρ<sub>Matern</sub></td><td>days</td><td></td><td>200</td></tr><tr><td></td><td></td><td></td><td></td><td>12</td></tr></tbody></table>

**Notes.** The reported errors are the 16th and 84th percentile interval for every parameter. (*) is for parameters that are only shared between CHEOPS and TESS. Priors based on previous results from (a) Wong et al. (2020), (b) Borsa et al. (2019), (c) this work (see Sect. 2.4), (d) Claret (2016), (e) restricted by break-up velocity.

et al. 2017) and CLES (Code Liègeois d'Évolution Stellaire, Scuflaire et al. 2008). In detail, we used the capability of the isochrone placement technique (Bonfanti et al. 2015, 2016) to fit the input parameters within pre-computed PARSEC grids of tracks and isochrones so to retrieve a first pair of age and mass estimates. A second pair of age and mass values, instead, was directly computed by the CLES code which generates the best evolutionary track that is compatible with the input parameters following a Levenberg-Marquadt minimisation scheme (see Salmon et al. 2021, for further details). After checking the mutual consistency of the two respective pairs through a χ²-based criterion, we finally merged our outcomes as described in Bonfanti et al. (2021) and we obtained $t_{\star} = 0.3 \pm 0.1$ Gyr and $M_{\star} = 2.45^{+0.19}_{-0.17} M_{\odot}$.

# 3. Analysis

## 3.1. Phase curves

Over its first and second year, CHEOPS has observed 4 full phase curves of KELT-9; TESS has observed over 20 and Spitzer has observed 1. We present here a joint analysis of all these observations.

As already discussed in Sect. 1, KELT-9b is the hottest known exoplanet with a dayside temperature of ~4600 K and a thermal emission peak in the CHEOPS bandpass as well as efficient optical absorbers. We therefore assume that all phase curves are completely dominated by thermal emission, an assumption justified in Schwartz et al. (2017) and Morris et al. (2022).

We fit each instrument's observations with a self-consistent transit, eclipse and thermal phase curve model. Table 3 shows the parameters common to all phase curves and Table 4 shows

parameters that are instrument-specific. In the following subsections, we detail the different components of the model used in the fitting procedure, along with the detrending and other models specific to each bandpass.

To fit our global model, we used an affine-invariant ensemble sampling to estimate the parameter posterior distributions with 900 walkers, 5000 burn-in steps and then 80 000 steps (Foreman-Mackey et al. 2013). We confirmed that the chains converged by analysing the autocorrelation time. The integration length is 50 times the number of samples.

### 3.1.1. Gravity-darkened transit model

Gravity can darken stellar photospheres due to rapid rotation, shaping the star into an oblate spheroid. von Zeipel (1924) showed that this oblateness causes a temperature gradient across the surface of the star and consequently reduces the emitted flux near the equator. Depending on the geometry of the star-planet system, this causes significant deviations from a simple symmetric transit light curve (Barnes 2009). This effect is in addition to limb darkening (Claret 2017). We must accurately model the transits in our observations in order to calculate the planetary radius and model the phase curve amplitude accurately, which in turn gives us information about the temperature map of the planet.

We model the transit of KELT-9b with `pytransit` (version 2.5.17, Parviainen 2015), an open-source Python package which implements a gravity-darkened transit model based on Barnes (2009). For describing this model we will use the same notation as in Hooton et al. (2022). The model is characterised by

- $R_*$, the stellar radius;

- $R_p$, the planetary radius (in units of stellar radii, $R_*$);
<table><caption>**Table 4.** Bandpass-specific priors and best-fit values for the model and detrending parameters as described in Sect. 3, along with the derived parameters.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Prior</th><th>Best-fit value</th></tr></thead><tbody><tr><td colspan="4"><strong>CHEOPS</strong></td></tr><tr><td colspan="4"><em>Fitted parameters</em></td></tr><tr><td>Zero transit epoch</td><td>t<sub>0</sub></td><td>BJD time - 2459095.2</td><td>U(-0.01, 0.01)</td><td>-0.003751<sup>+0.000060</sup><sub>-0.000059</sub></td></tr><tr><td>Transit depth</td><td>depth</td><td></td><td>U(0.0049, 1)</td><td>0.006212<sup>+0.000019</sup><sub>-0.000020</sub></td></tr><tr><td>Transit scaling factor</td><td>f<sub>0</sub></td><td></td><td>U(0.9, 1.1)</td><td>0.99998<sup>+0.00012</sup><sub>-0.00011</sub></td></tr><tr><td>1st quadratic limb darkening component</td><td>q<sub>1</sub></td><td></td><td>N(0.34052, 0.01)<sup>(a)</sup></td><td>0.3299<sup>+0.0078</sup><sub>-0.0074</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>q<sub>2</sub></td><td></td><td>N(0.22030, 0.01)<sup>(a)</sup></td><td>0.2210<sup>+0.0086</sup><sub>-0.0087</sub></td></tr><tr><td>Hotspot offset</td><td>Δφ</td><td>deg</td><td>U(-180, 29)</td><td>-14.1 ± 2.4</td></tr><tr><td>C<sub>ml</sub> power coefficient</td><td>C<sub>1,1</sub></td><td></td><td>U(0, 1)</td><td>0.205<sup>+0.039</sup><sub>-0.028</sub></td></tr><tr><td>Mean background temperature</td><td>T&#773;</td><td>K</td><td>U(2890, 5780)<sup>(d)</sup></td><td>4060<sup>+120</sup><sub>-150</sub></td></tr><tr><td colspan="4"><em>Derived parameters</em></td></tr><tr><td>Planetary radius</td><td>R<sub>p</sub></td><td>R<sub>*</sub></td><td></td><td>0.07882<sup>+0.00012</sup><sub>-0.00013</sub></td></tr><tr><td>1st quadratic limb darkening component</td><td>u<sub>1</sub></td><td></td><td></td><td>0.2541<sup>+0.0094</sup><sub>-0.0097</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>u<sub>2</sub></td><td></td><td></td><td>0.320 ± 0.012</td></tr><tr><td colspan="4"><strong>TESS</strong></td></tr><tr><td colspan="4"><em>Fitted parameters</em></td></tr><tr><td>Zero transit epoch</td><td>t<sub>0</sub></td><td>BJD time - 2458693.8</td><td>U(-0.01, 0.01)</td><td>0.013185<sup>+0.000044</sup><sub>-0.000045</sub></td></tr><tr><td>Transit depth</td><td>depth</td><td></td><td>U(0.0049, 1)</td><td>0.006263 ± 0.000016</td></tr><tr><td>Transit scaling factor</td><td>f<sub>0</sub></td><td></td><td>U(0.9, 1.1)</td><td>1.00017 ± 0.00012</td></tr><tr><td>1st quadratic limb darkening component</td><td>q<sub>1</sub></td><td></td><td>N(0.1690, 0.01)<sup>(b)</sup></td><td>0.1541<sup>+0.0062</sup><sub>-0.0064</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>q<sub>2</sub></td><td></td><td>N(0.2082, 0.01)<sup>(b)</sup></td><td>0.2175<sup>+0.0094</sup><sub>-0.0097</sub></td></tr><tr><td>Hotspot offset</td><td>Δφ</td><td>deg</td><td>U(-180, 29)</td><td>-12.32 ± 0.97</td></tr><tr><td>C<sub>ml</sub> power coefficient</td><td>C<sub>1,1</sub></td><td></td><td>U(0, 1)</td><td>0.1884<sup>+0.0088</sup><sub>-0.0083</sub></td></tr><tr><td>Mean background temperature</td><td>T&#773;</td><td>K</td><td>U(2890, 5780)<sup>(d)</sup></td><td>3955<sup>+40</sup><sub>-42</sub></td></tr><tr><td colspan="4"><em>Derived parameters</em></td></tr><tr><td>Planetary radius</td><td>R<sub>p</sub></td><td>R<sub>*</sub></td><td></td><td>0.079142 ± 0.000099</td></tr><tr><td>1st quadratic limb darkening component</td><td>u<sub>1</sub></td><td></td><td></td><td>0.1705<sup>+0.0076</sup><sub>-0.0075</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>u<sub>2</sub></td><td></td><td></td><td>0.2216 ± 0.0096</td></tr><tr><td colspan="4"><strong>Spitzer</strong></td></tr><tr><td colspan="4"><em>Fitted parameters</em></td></tr><tr><td>Zero transit epoch</td><td>t<sub>0</sub></td><td>BJD time - 2458415.36</td><td>U(-0.01, 0.01)</td><td>0.00206 ± 0.00015</td></tr><tr><td>Transit depth</td><td>depth</td><td></td><td>U(0.0049, 1)</td><td>0.006152 ± 0.000068</td></tr><tr><td>Transit scaling factor</td><td>f<sub>0</sub></td><td></td><td>U(0.9, 1.1)</td><td>0.99993 ± 0.00020</td></tr><tr><td>1st quadratic limb darkening component</td><td>q<sub>1</sub></td><td></td><td>N(0.0097, 0.01)<sup>(c)</sup></td><td>0.0118<sup>+0.0065</sup><sub>-0.0054</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>q<sub>2</sub></td><td></td><td>N(0.1463, 0.01)<sup>(c)</sup></td><td>0.1464<sup>+0.0099</sup><sub>-0.0101</sub></td></tr><tr><td>Hotspot offset</td><td>Δφ</td><td>deg</td><td>U(-180, 29)</td><td>-18.2<sup>+1.7</sup><sub>-1.6</sub></td></tr><tr><td>C<sub>ml</sub> power coefficient</td><td>C<sub>1,1</sub></td><td></td><td>U(0, 1)</td><td>0.2453 ± 0.0092</td></tr><tr><td>Mean background temperature</td><td>T&#773;</td><td>K</td><td>U(2890, 5780)<sup>(d)</sup></td><td>3934<sup>+54</sup><sub>-58</sub></td></tr><tr><td colspan="4"><em>Derived parameters</em></td></tr><tr><td>Planetary radius</td><td>R<sub>p</sub></td><td>R<sub>*</sub></td><td></td><td>0.07843 ± 0.00043</td></tr><tr><td>1st quadratic limb darkening component</td><td>u<sub>1</sub></td><td></td><td></td><td>0.0318<sup>+0.0082</sup><sub>-0.0084</sub></td></tr><tr><td>2nd quadratic limb darkening component</td><td>u<sub>2</sub></td><td></td><td></td><td>0.077<sup>+0.019</sup><sub>-0.020</sub></td></tr></tbody></table>

**Notes.** The reported errors are the 16th and 84th percentile interval for every parameter. Priors based on previous results from (a)Claret (2021), (b)Claret (2017), (c)Claret & Bloemen (2011), (d)Morris et al. (2022).

- *P*, planetary orbital period;

- *ρ*, the stellar density;

- *a*, the semi-major axis in units of $R_*$;

- *i*, the planetary orbital inclination;

- *e*, orbital eccentricity;

- *ω*, argument of periastron;

- *P<sub>rot</sub*, the rotation period of the star;

- *T<sub>pole</sub*, the stellar temperature at its pole;

- *λ*, the sky-projected spin orbit angle (the sky-projected angle between the planetary orbital plane and the stellar equatorial plane);

- *i*\*, stellar inclination, defined as the angle between the observer's line of sight and the stellar rotation axis. This is related to the stellar obliquity *φ*, specified in Ahlers et al. (2020), by *i*\* = 90° – *φ*. We note this is not the same *Δφ* that we use to describe the thermal phase curve model;

- *β*, the gravity-darkening coefficient (defines how strong the temperature gradient is), defined in von Zeipel (1924);

- *u*₁ and *u*₂, the quadratic limb-darkening coefficients, which we reparametrise as *q*₁ and *q*₂ (Kipping 2013);

- *t*₀, the time of mid-transit;

- *f*₀, transit scaling factor (scales the out-of-transit baseline);

- *filter*, the telescope's bandpass transmission efficiency;

- stellar spectrum.

We fixed $T_{\text{pole}} = 10170$ K, $e=0$ and $\omega=90^\circ$. We also fixed $R_* = 2.379 R_\odot$ (see Sect. 2.4). For CHEOPS and TESS we use a PHOENIX model stellar spectrum (Husser et al. 2013), and for Spitzer we use a blackbody spectrum.

Due to the large degeneracy between the limb darkening and gravity darkening parameters, we had to use previous studies to inform our priors, instead of just letting them all float freely – something we tested and found no clear unique solution. Therefore we used results from previous Doppler tomography studies
(Gaudí et al. 2017; Borsa et al. 2019) to inform priors on b, λ and v sin i* (where v is the rotational velocity of the star), and theoretical limb-darkening tables in Claret & Bloemen (2011); Claret (2017, 2021) to inform priors on q₁ and q₂. We also used Claret (2016) to inform priors on β. See Tables 3 and 4 for the priors used for all the gravity-darkened model parameters.

Using the information obtained from the transit fit we can calculate the true 3D spin-orbit angle (Ψ), given by the equation:

$$ \cos \Psi = \sin i^* \cos i + \cos i^* \sin i \cos \lambda. \quad (1) $$

Gravity-darkened transits have also been analysed in other CHEOPS work (see, e.g.: Lendl et al. 2020; Hooton et al. 2022; Deline et al. 2022).

### 3.1.2. Thermal phase curve model

To model the shape of the thermal phase curve we use kelp (Morris et al. 2022), a Python package that models the surface temperature of the planet as a 2D thermal map constructed by modified spherical harmonics (parabolic cylinder functions) (Heng & Workman 2014), given by Eq. (1) of Morris et al. (2022):

$$ T(\theta, \phi) = \bar{T} \left( 1 + \sum_{m,l}^{\ell_{\max}} h_{ml}(\theta, \phi) \right), \quad (2) $$

where θ and φ are latitude and longitude. $\bar{T}$ is a background temperature upon which 2D perturbations exist and are quantified by the “alphabet” or basis functions $h_{ml}$ as defined by Eq. (258) of Heng & Workman (2014). These functions are also dependent on the dimensionless parameters α and ω<sub>drag</sub>. The dimensionless fluid number α is constructed from the Rossby, Reynolds and Prandtl numbers; it quantifies the competition between stellar heating (forcing) and sources of friction (drag). The normalised drag frequency ω<sub>drag</sub> quantifies the strength of friction. In practice, we find that it is sufficient to truncate the series in the preceding equation at $l_{max} = 1$ (Morris et al. 2022). The entire temperature map may be shifted back and forth in longitude by Δφ, which effectively fits for the hotspot offset if ω<sub>drag</sub> ≥ 3 (Morris et al. 2022). Full details of the theoretical formalism may be found in Sect. 2 of Morris et al. (2022), which we will not repeat here.

Overall, our thermal phase curve model is parametrised by the following variables:

- Δφ, hotspot offset;

- α, dimensionless fluid number;

- ω<sub>drag</sub>, dimensionless drag frequency;

- C<sub>ml</sub>, the power of individual harmonic modes;

- $\ell_{max}$, describes the highest spherical harmonic mode in the model;

- planetary parameters including a, R<sub>p</sub>, T<sub>pole</sub>;

- telescope bandpass transmission efficiency;

- $\bar{T}$, scaling term of the mean temperature field.

kelp integrates over the passband-weighted thermal map visible at a specific time in its orbit and converts the map into a flux measurement using the following equation from Cowan & Agol (2011), which can be compared to the observed phase curve:

$$ \frac{F_p}{F_*} = \frac{1}{\pi I_*} \left( \frac{R_p}{R_*} \right)^2 \int_0^\pi \int_{-\xi-\pi/2}^{-\xi+\pi/2} I_p(\theta, \phi) \cos(\phi + \xi) \sin^2(\theta) \, d\phi \, d\theta \quad (3) $$

where I, the intensity, is defined by

$$ I = \int \lambda F_{\lambda} B_{\lambda}(T(\theta, \phi)) d\lambda \quad (4) $$

where $F_\lambda$ is the instrument-specific filter response function and $B_\lambda(T(\theta, \phi))$ is the Planck function of each temperature element T(θ, φ).

The initialisation of this model at every time-step required is a bottleneck in the joint MCMC fit for the phase curves. To increase the speed we evaluated the model at 200 orbital phases, spaced equally between −π and π. We then linearly interpolated the model flux values back to the original time resolution. We investigated the effect of increasing the number of samples on the fitted parameters and found that the median parameter values converged to a constant value at around 150, and so we choose 200 samples to increase the speed of the code whilst not decreasing the accuracy of the phase curve fit.

In the joint fit, we followed Morris et al. (2022) and set α = 0.6, ω<sub>drag</sub> = 4.5 and $\ell_{max} = 1$, since latitudinal variations in temperature and flux are not constrained by thermal phase curves and it was shown in that study that it is sufficient to describe the temperature map only with the first mode. We also set all C<sub>ml</sub> = 0 apart from C<sub>11</sub> which is a free parameter. See Table 3 for the priors of the phase curve model.

### 3.1.3. Secondary eclipse model

The batman package models the flux decrement as the planet is occulted by the star (Kreidberg 2015). We multiply this normalised model (where the planetary flux is null during the eclipse) by the thermal phase curve model to produce the full phase curve model.

### 3.1.4. Stellar pulsation model

After examining the residuals of an initial CHEOPS phase curve fit, without any stellar pulsation model, it was clear that there was a periodic signal present in the phase curves at a period of around 7.5 hours (see Fig. 1). This signal was identical to the one observed in Wong et al. (2020) and Mansfield et al. (2020) and has been attributed to stellar pulsations. In Wyttenbach et al. (2020), they analysed these pulsations and concluded that they are compatible with p-mode oscillations present in a δ Scuti-type star.

To correct for this signal we used a Gaussian process with a simple harmonic oscillator (SHO) kernel implemented by celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). We fixed the amplitude (100 ppm) and the damping timescale, τ, to 2× the Gaussian process period (GP period), which is a fitted parameter. This damping timescale provides coherence of pulsations over several cycles while allowing for evolution of the pulsation signal on longer timescales. Because of this we calculate the Q value of this SHO kernel to be 6.

We found evidence of the same stellar pulsations in the TESS phase curves as well. Although the PDCSAP flux light curves have already been corrected for long-term trends, there were still significant trends present in the KELT-9 light curve which we believe to be of instrument systematic origin. Therefore we implemented the same kernel described above for these TESS phase curves, along with an additional Matérn 3/2 kernel with an amplitude of 200 ppm and a timescale of 12 days to remove the long-term trends in this data.

The pulsations are also detectable in the Spitzer data, albeit at less than 10% of the total amplitude of the phase curve (a lower percentage than in the other two bandpasses). After further
<figure><img src="image_3.png" /></figure>

**Fig. 1.** Periodograms reveal the presence of stellar pulsations in the residuals of the phase curve fitting. In red are the four CHEOPS phase curves, in green are the TESS phase curves and in blue is the *Spitzer* phase curve. It is clear the stellar pulsation period (near 0.3 days) is present in all phase curves and matches that signal seen also in Wong et al. (2020) and Mansfield et al. (2020) (dashed black line). The short period peaks are only present in CHEOPS but not TESS or *Spitzer* so we assume they are CHEOPS systematics and disregard the signal.

analysis showed modelling the pulsations within this dataset had a non-trivial impact on the phase curve parameters, we extended the GP to also include the *Spitzer* phase curve. The amplitude of the pulsations in the *Spitzer* data was also around 100 ppm so this hyperparameter could remain the constant. We decided to use the same kernels in both CHEOPS, TESS, and *Spitzer* to jointly infer the free kernel hyperparameters.

### 3.1.5. CHEOPS

The composite transit, eclipse and thermal phase curve model we use to fit the CHEOPS observations is composed of:

- gravity-darkened transit model;

- planet thermal phase curve model;

- eclipse model;

- stellar pulsation model;

- systematics model.

**Data clipping and systematics model.** For each phase curve, before stitching them together with the other datasets, we sigma-clipped outliers of the centroid position of the target star at $4.5\sigma$ and masked out individual points with anomalously high levels of background or low temperature readings. As discussed in Sect. 2.1, these telescope assembly temperatures have been shown to coincide with a systematic “ramp” feature in the CHEOPS light curves, most likely caused by the change in temperature of the telescope as it adjusts to a new pointing position (Morris et al. 2022). We also masked out the first few hours of data in the second, third and fourth phase curve, as they coincided with a strong, non-linear increase in the telescope temperature. We mask on average 16% of datapoints from each phase curve.

We used a linear detrending model including vectors proportional to the temperature, temperature² and roll angle in our fit to detrend against the correlated noise sources. We chose these basis vectors based on the analysis performed in Morris et al. (2022) which used the BIC minimisation to find the parameters that contributed the most to the model without overfitting.

### 3.1.6. TESS

The TESS phase curves model is a linear combination of:

- gravity-darkened transit model;

- planet thermal phase curve model;

- eclipse model;

- stellar pulsation model.

**Data clipping.** As described in Sect. 2.2, we use the TESS PDCSAP light curves from the SPOC pipeline. Before the fit we also masked out sections of the phase curves that clearly had strong systematics and that were also removed in the analysis performed in Wong et al. (2020). These usually affected the data points shortly before or after a gap in the TESS observations. In total we mask out 5% of the datapoints.

### 3.1.7. Spitzer

The *Spitzer* phase curve model is a linear combination of:

- gravity-darkened transit model;

- planet thermal phase curve model;

- eclipse model;

- stellar pulsation model.

**Pre-fit conditioning and systematics model.** See Sect. 2.3 for a description of the detrending applied to the *Spitzer* dataset. From the BLISS mapping detrending routine, we obtain a detrended lightcurve, which we use for our joint phase curve fit. Observations are also available at 3.6 µm, which will be the subject of a forthcoming paper by Beatty et al. (in prep.).

### 3.1.8. On ellipsoidal variations

To include an ellipsoidal variation model (see, e.g.: Welsh et al. 2010; Gai & Knuth 2018), we initially fit the first two CHEOPS phase curves separately using the Wong et al. (2020) sinusoidal model for the stellar pulsations and ellipsoidal variations, leaving the phase and amplitude of both sinusoids as free parameters. The resulting ellipsoidal phase and amplitude were not consistent between each phase curve. We then fit the combined set of CHEOPS phase curves and the chains converged on solutions with ellipsoidal variation amplitudes of <10 ppm. The theoretical expectation is reported as $44 \pm 6$ ppm in Wong et al. (2020), therefore this fitted amplitude is over $5\sigma$ lower than the theoretical value. We suggest that the variability on the $P/2$ timescale may not be ellipsoidal because it has an unexpected and evolving phase and appears incoherent over 4 orbital timescales (from the first to second phase curve observations).

Carrying out the same test with only the GP stellar pulsation model instead of the pure sinusoids, we found that both models agreed on the amplitude of the stellar pulsations as well as other planetary parameters. As a result of this analysis, we conclude that the GP stellar pulsation model is flexible enough to capture the stellar signals without absorbing the planetary ones.

## 3.2. CHEOPS occultations

As well as the 9 occultation-only CHEOPS observations, 5 eclipses were observed within the 4 CHEOPS phase curves (2 of the 5 eclipses were observed within the same phase curve). To increase our occultation sample size, we added these to the set of occultations, clipping them so that they each had a similar number of CHEOPS orbits to the rest of the occultations. Therefore, there are 13 occultations in total.

Each occultation was then fitted independently with a model containing a linear combination of:

- eclipse model;

- stellar pulsation model;

- systematics model.
We carried out the fit of the occultations with PyMC3 (Salvatièur et al. 2016), which uses a gradient-based Hamiltonian Monte Carlo (HMC) method to integrate for parameter posteriors.

*Data clipping and systematics model.* Similarly to the phase curves, the occultations were sigma-clipped to remove outliers in the target centroid-space and flux values (both by $3\sigma$), and points with anomalously large background or temperature readings were masked out. The data points were also flux sigma-clipped. After clipping, 4–14% of points were removed from each observation.

Since the observation duration for each occultation is short (4–5 CHEOPS orbits), we risk fitting a model that is too complex. To investigate the minimum complexity model needed to reproduce the light curve without over-fitting, we used Leave-One-Out cross-validation (Vehtari et al. 2015) to compare the predictive power for each model containing different basis vectors in the systematics model. The basis vectors we tested included a flux constant, sin (roll angle), cos (roll angle), time, time², temperature (from the *thermFront2* sensor) and background. These have been shown in a previous study (Morris et al. 2022) to be the most influential basis vectors on the light curves.

Overall, the preferred systematic model included every basis vector except the background, however other combinations were preferred by nine of the other occultations. After fitting, we investigated whether there was a trend between the fitted eclipse depth and the number of basis vectors chosen for each occultation. We found no correlation.

### 3.2.1. Stellar pulsation model

As detailed in Sect. 3.1.5, We observed a stellar pulsation signal in the phase curves with an amplitude of around 100 ppm and a timescale of 7.5 h. Assuming it is continuous and also present in the occultation observations, the occultations are taken over a long enough timescale that this signal would vary significantly within a single visit. However as the baseline of the occultations is so short, it is impossible to uniquely infer the phase and amplitude of the pulsations in each observation. Therefore we included a sinusoid in the occultation models with a period fixed at 7.5 h, with the phase of this signal free to vary independently between each occultation. In order to avoid biasing the prior, we implemented a hierarchical Bayesian technique where the prior of the amplitude was set to $N(100, \sigma)$ ppm, and truncated at zero so that the amplitude is always positive, and the standard deviation $\sigma$ became a fitted parameter and was allowed to vary as $U(10, 1000)$ ppm.

### 3.2.2. Eclipse model

We used the `batman` package to create a basic secondary eclipse model for the occultations. The out-of-eclipse observations contain a hint of the shape of the phase curve of KELT-9, which peaks near secondary eclipse. We used the best-fit posteriors from the full phase curve fit to produce a phase curve model from `kelp`. This model is used to scale the out-of-eclipse sections of the basic `batman` model. This model was characterised by the following parameters:

- $t_0$, the time of transit centre*;

- $R_p$, the planetary radius (in units of stellar radii, $R_*$);

- $b$, the impact parameter;

- $P$, the planetary orbital period*;

- $\rho$, the stellar density.

where we fixed starred (*) values and fit for unstarred ones.

## 4. Results

Figure 2 shows the phase folded CHEOPS, TESS, and *Spitzer* phase curves respectively, with the stellar pulsations and other systematic trends removed, and overplotted with the best-fit phase curve, transit and eclipse model. The fitted parameter posteriors are detailed in Table 3 for global variables and Table 4 for bandpass-dependent variables. A clearer view of the transit fits can be seen in Appendix A. The results of the occultation depth analysis are shown in Fig. 8.

In the phase curve analysis, we included both year 1 and year 2 CHEOPS phase curves. We analysed these years of phase curves separately and we saw no evidence of variability from one CHEOPS year to the next. The change in parameter values from using only year 1 to including both years 1 and 2 was, for the majority of parameters, within $1\sigma$ and the rest was within $2\sigma$. This justified our decision to use the same parameters in the fitting procedure for all CHEOPS phase curves.

### 4.1. Thermal map of KELT-9b

As described in Sect. 3.1.2, we use a generalised spherical harmonic temperature map to fit the thermal phase curve variation in the light curves. To a good first-approximation, there are only three free variables that can describe most hot-Jupiter phase curves: $\Delta\phi$ (hotspot offset), $C_{1,1}$ ($C_{ml}$ power coefficient) and $\bar{T}$ (mean background temperature). As each instrument observes in a different bandpass, they are probing different atmospheric depths of KELT-9b and so we used an independent phase curve model for each bandpass. Although Hooton et al. (2022) used spherical harmonic temperature maps in a similar way to fit the phase variations in CHEOPS and *Spitzer* light curves for MASCARA-1b, they used a single set of the same three parameters for bandpasses due to the limited phase coverage of the CHEOPS observations. Figure 4 shows the temperature maps derived from the best-fit phase curve models for each instrument. Similar to other hot Jupiters, the result in all three bandpasses is a hotspot slightly offset to the east. Furthermore, the minimum nightside temperature remains above 2000 K in all bandpasses, further evidence for significant heat-redistribution in the atmosphere. It is also apparent that in the TESS data, there is a smaller temperature contrast between the maximum and minimum temperature than the other datasets.

In Fig. 3, we illustrate the fitted posteriors of the three free variables in these phase curve models. Table 4 details the best-fit values of these parameters along with the other bandpass-dependent parameters. It is clear that these parameters change significantly as we probe different depths of KELT-9b's atmosphere.

The CHEOPS nightside posteriors are considerably broader than those of TESS and *Spitzer*. The TESS data has around 50 days of photometry compared to the ~10 days of CHEOPS, which contributes to narrower TESS posteriors. The CHEOPS nightside flux (in ppm) is also lower in the CHEOPS observations and so is detected at lower significance than the other instruments. Furthermore, *Spitzer* observes a larger planet-to-star contrast than the other two instruments and also other factors such as a larger collecting area and the longer observing durations all contribute to the narrower posteriors. This may also be related to the fact that we detrend the CHEOPS data simultaneously in the fitting procedure but not the TESS and *Spitzer* data. However, we try to account for this by fitting for an additional white noise term to each individual data set, modifying the uncertainties based on the standard deviation of the flux values.
<figure><img src="image_4.png" /><figcaption>(a) CHEOPS phase-folded and detrended phase curves (first 2 panels, 2nd panel is a zoomed-in view of the 1st). 3rd panel shows residuals of the fit.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) TESS phase-folded and detrended phase curves (first 2 panels, 2nd panel is a zoomed-in view of the 1st). 3rd panel shows residuals of the fit.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) *Spitzer* phase-folded and detrended phase curve (first panel). 2nd panel shows residuals of the fit.</figcaption></figure>

**Fig. 2.** CHEOPS, TESS, and *Spitzer* phase-folded and detrended (stellar pulsations and other systematic trends removed) phase curves, overplotted with the best-fit phase curve model, transit model and eclipse model (in red). In black are the binned grey datapoints, with error bars that are smaller than the point size in all panels so they are not visible.
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Posteriori of the phase curve parameters across the three band-passes showing $1\sigma$-confidence contours. Red is CHEOPS, green is TESS and blue is *Spitzer*.</figcaption></figure>

In the second column of Fig. 4, we show the ratio of atomic hydrogen to molecular hydrogen in the atmosphere of KELT-9b, assuming equilibrium chemistry and using the method and equations described in Heng et al. (2016). This supports the theory that the dissociation from atomic to molecular hydrogen occurs near the day-night terminator for the wavelengths observed.

## 4.2. Dayside and nightside brightness temperatures

Although our work reveals 2D information regarding the planet's temperature profile, we believe it is useful to still report the dayside and nightside brightness temperatures that we infer directly from two single ("1D") measurements: the eclipse depth and the flux at half an orbital period away from the centre of the eclipse. Assuming a model spectrum for the star, a black body for the planet and by using each instrument's filter function, these brightness temperatures can be calculated. Although this is not the main focus of this paper, it is useful to report these temperatures due to the wide understanding of this temperature statistic in the community, along with the opportunity for direct comparison to previous work.

We obtain a dayside brightness temperature of CHEOPS, TESS, and *Spitzer* of 4796 ± 46 K, 4643 ± 26 K and 4870<sup>+67</sup><sub>-65</sub> K respectively. We found the nightside brightness temperatures of the three bandpasses (in the same order) to be 3180<sup>+190</sup><sub>-230</sub> K, 3177<sup>+60</sup><sub>-61</sub> K and 2973<sup>+66</sup><sub>-62</sub> K. Figure 5 shows the posteriors of these wavelength-dependent temperatures (derived from the posteriors of the thermal phase curve parameters). It is clear that the highest average dayside temperature is observed in the *Spitzer* bandpass, followed by CHEOPS and then TESS. Generally, the nightside temperatures also have wider posteriors as the lower nightside flux is detected at a lower significance.

Our TESS dayside brightness temperature is consistent with the Wong et al. (2020) value of 4600 ± 100 K and our TESS nightside temperature is also consistent with their nightside temperature of 3040 ± 100 K. This is an encouraging test of our full phase curve model as both of these studies used the same TESS data. For our *Spitzer* dayside temperature, our results are consistent within $2\sigma$ to the dayside brightness temperature

reported in Mansfield et al. (2020), 4566<sup>+140</sup><sub>-136</sub> K. However our nightside temperature is not consistent with their nightside brightness temperature of 2556<sup>+101</sup><sub>-97</sub> K (difference of over $3\sigma$). This could in part be due to the different detrending method applied to *Spitzer* data in this study compared to the Mansfield et al. (2020) study.

Figure 6 shows the nightside temperatures of the most extreme hot Jupiters plotted against their equilibrium temperatures and Fig. 7 shows the dayside temperature of these planets against equilibrium temperature. We define the equilibrium temperature as $T_{eq} = T_* \sqrt{R_* / 2a}$. It has been reported in Keating et al. (2019) that due to the presence of clouds the nightside temperature of most hot Jupiters is around 1100 K, with the caveat that clouds would disperse for hotter planets, and therefore the nightside temperature may increase again proportionally to the amount of stellar irradiation. KELT-9b supports this caveat as it is clearly an exception to the pattern of constant nightside temperature, having a nightside temperature exceeding 2900 K.

## 4.3. Eclipse depths

Figure 8 shows the eclipse depths fitted from the CHEOPS occultation observations. As explained in Sect. 3.2, this dataset includes 9 occultation-only observations and four occultations cut-out from the four CHEOPS phase curves. Together they have a mean eclipse depth of 320 ± 11 ppm. This is consistent within $2.3\sigma$ of the eclipse depth found from the full phase curve fit of the four CHEOPS phase curves which was 367 ± 17 ppm. The mean eclipse depth variation (from the occultations) corresponds to a brightness temperature change of around ±33 K.

To validate our eclipse model uncertainties, we fit the same size error bar to each eclipse depth datapoint at the same time as performing a straight-line fit to the eclipse depths. This fitted error bar amplitude was similar to the average of our occultation fit error bars, therefore we suggest that the error bars are appropriate to justify our results.

As described in Sect. 3.2, fitting the eclipses is particularly challenging due to lack of data out-of-eclipse. This makes finding the correct phase and amplitude of the stellar variations, which we know are present from the CHEOPS, TESS, and *Spitzer* phase curves, very difficult. As the amplitude and half-period of these variations are on the same order as the eclipse depth and duration, this inflates the uncertainty in the eclipse depth measurements, and including them in the eclipse depth fit allows the reported error bars to reflect this uncertainty. This may also be an explanation as to why the mean occultation depth is less than the estimated occultation depth from the full phase curve-only fits, as the stellar pulsations being fited can mimic the dip of the eclipse and produce a good fit with an anomalously low eclipse depth. It is worth noting that this difference is probably not due to the difference in modelling of the out-of-eclipse flux between the occultation-only observations and the occultations within the phase curves as we used a phase curve model in the occultation-only analysis as well (see Sect. 3.2.2). In models without this stellar pulsation model, the eclipse depths we retrieve are very different from one another and from the analysis with the pulsation model, and the errors reported were considerably smaller than the scatter in the depths. In this case it is clear that phase curve observations have been vital for informing the priors and model of the occultation-only observations. In future CHEOPS projects, one must be extremely cautious when working with occultations from variable and pulsating stars. Phase curves are essential in
<figure><img src="image_8.png" /><figcaption>Fig. 4. 2D thermal temperature maps (**left column**) and colour bar) and maps showing ratio of atomic to molecular hydrogen, assuming chemical equilibrium, that would be present given the temperature maps (**right column**) and colour bar) of KELT-9b in each bandpass. The hottest temperature is reached near the substellar point and is just over 5000 K whereas the lowest temperature is reached in the *Spitzer* bandpass at just under 2500 K. On all three maps, the dayside is mostly atomic hydrogen-saturated while across the terminator and over to the nightside, the ratio falls to near zero at the antistellar point.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Posteroirs of the brightness day and nightside temperatures from the different bandpasses.</figcaption></figure>

this case to constrain this stellar source of variability, due to the very limited baseline of occultation observations.

These eclipse depths, especially with the addition of the three occultations observed a year later, suggests a lack of significant variability in the atmosphere of KELT-9b. The variation in eclipse depths is roughly consistent with the error bar of the mean depth. Therefore we set an upper limit of temperature variability of KELT-9b at 1% of the mean brightness temperature (~1σ). This observation is consistent with the lack of variability observed by *Spitzer* for HD 189733 b by Agol et al. (2010), for

HD 189733 b and HD 209458 b by Kilpatrick et al. (2020), and with theoretical expectations by, for example, Showman et al. (2009) and Komacek & Showman (2020).

### 4.4. Albedo and heat redistribution

A thermal phase curve constrains the global dayside and nightside brightness temperatures. If the atmosphere radiates like a perfect blackbody, then the brightness temperature is equal to the real temperature, despite the thermal phase curve being observed only within a limited range of wavelengths. The traditional approach is to then use a 0D “box model” to convert these temperatures into the Bond albedo and heat redistribution efficiency, e.g., Eqs. (4) and (5) of Cowan & Agol (2011). We do not expect KELT-9b to radiate like a perfect blackbody, but our approach for converting the temperature map to fluxes assumes a Planck function (Morris et al. 2022). Generally, the CHEOPS, TESS, and *Spitzer* thermal phase curves are probing different atmospheric layers (across radial distance or atmospheric pressure), which are described by different temperature maps.

We utilise the method described in Morris et al. (2022) to use the entire 2D temperature maps derived from the phase curve fitting to calculate these values. Following this, the Bond albedo is defined as

$$A_B = 1 - \left(\frac{a}{R_*}\right)^2 \frac{\int_{-\pi}^{\pi} \int_{0}^{\pi} F_P(\theta, \phi) \sin\theta \, d\theta \, d\phi}{\pi \sigma T_*^4}, \quad (5)$$
<table><caption>**Table 5.** Fitted eclipse depths and night fluxes derived directly from the phase curve fitting procedure for each of the three different bandpasses.</caption><thead><tr><th></th><th>Eclipse depth (ppm)</th><th>Night flux (ppm)</th><th>Dayside temperature (K)</th><th>Nightside temperature (K)</th></tr></thead><tbody><tr><td>CHEOPS</td><td>367 ± 17</td><td>37<sup>+16</sup><sub>-14</sub></td><td>4796 ± 46</td><td>3180<sup>+190</sup><sub>-230</sub></td></tr><tr><td>TESS</td><td>645 ± 15</td><td>105<sup>+11</sup><sub>-11</sub></td><td>4643 ± 26</td><td>3177<sup>+60</sup><sub>-61</sub></td></tr><tr><td>Spitzer</td><td>3007<sup>+57</sup><sub>-55</sub></td><td>1440<sup>+53</sup><sub>-49</sub></td><td>4870<sup>+67</sup><sub>-65</sub></td><td>2973<sup>+66</sup><sub>-62</sub></td></tr></tbody></table>

**Notes.** The last two columns show the corresponding dayside and nightside brightness temperatures, derived using a PHOENIX model for KELT-9 and the respective bandpass filter functions.

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Nightside integrated temperatures of hot Jupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Dayside integrated temperatures of hot Jupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.</figcaption></figure>

where $\sigma$ is the Stefan-Boltzmann constant and $F_p(\theta, \phi)$ is the flux from the planet derived from the temperature map. We must note here that the Bond albedo is a wavelength-independent quantity, however in this paper we quote three Bond albedos, one derived from each temperature map. We do this by assuming that each temperature element from the temperature map behaves like a blackbody. For each temperature element the flux of that element is estimated to be equal to $F_p(\theta, \phi) = \sigma T_p^4(\theta, \phi)$. This is not a perfect assumption, and so the Bond albedos reported from each temperature map may vary slightly (and can go negative) due to non-blackbody behaviour of the atmosphere. For example, if, at a certain wavelength, we are observing a strong absorption feature, then the estimated bolometric flux and brightness temperature may be less than the actual bolometric flux, as we assume that the small wavelength-band we observe is representative of the entire

spectrum. The Bond albedos reported using the temperature map from each bandpass is reported in Table 6.

Considering the Bond albedo reported from each temperature map, there is strong evidence to suggest it is consistent with zero. This result is expected due to the extreme level of irradiation on the planet, which would contribute to a highly extended atmosphere where light entering from the direction of the star would have a very low probability of escaping before being absorbed. This therefore implies KELT-9b predictably, behaves similar to a blackbody.

The heat redistribution parameter of the atmosphere is defined as the ratio of the nightside flux to the dayside flux and describes the extent to which heat is transported around the planet. As KELT-9b is tidally locked, this parameter is entirely dependent on the dynamics and chemistry of the atmosphere. It is also derived using information from the entire temperature map using the equation from Morris et al. (2022):

$$ \varepsilon = \frac{\int_{\pi/2}^{\pi} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi + \int_{-\pi}^{-\pi/2} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi}{\int_{-\pi/2}^{\pi/2} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi} \quad (6) $$

Figure 9 shows the heat redistribution ($\varepsilon$) of KELT-9b plotted along with other hot-Jupiters. In previous papers (Wallack et al. 2021; Komacek & Showman 2016), it has been shown that for hot Jupiters, the incident stellar flux is the primary decider of the level of $\varepsilon$. However in Bell & Cowan (2018), they predict a rising heat redistribution for UHJ due to the $H_2$ dissociation and recombination increasing the heat transport around the planet. Our work agrees with this theory and in Fig. 9, it appears that $\varepsilon$ does indeed fall with planet equilibrium temperatures up to around 2500 K, however after that the $\varepsilon$ rises again with temperature and KELT-9b's $\varepsilon$ calculated in this paper supports this trend.

## 4.5. Comparison with Spitzer/TESS literature

### 4.5.1. TESS

We find that for the phase curve parameters, our TESS hotspot offset of $-12.32 \pm 0.97^\circ$ is inconsistent with the offset reported in Wong et al. (2020, hereafter W20) ($5.2 \pm 0.9^\circ$) by over $5\sigma$ (N.B. this error is a combination of the uncertainty in this paper's value and W20). However we obtain larger values for the offset in all passbands, and the Spitzer result is consistent with previous analyses. This lends credibility to our TESS measurement. For the other orbital parameters our planetary radius ($R_p = 0.079142 \pm 0.000099 R_\star$) in the TESS bandpass is consistent with W20.

For some of the global parameters, our best-fit impact parameter ($b = 0.195^{+0.016}_{-0.015}$) differs from the W20 value by $1.4\sigma$. Our best-fit semi-major axis ($a = 3.0914^{+0.0090}_{-0.0100} R_\star$) is different from the W20 value by around $4\sigma$ and our best-fit period ($P = 1.48111949 \pm 0.00000034$ days) is different to W20's value
<figure><img src="image_12.png" /></figure>

**Fig. 8.** Fitted eclipse depths from the 9 CHEOPS occultation observations (in blue) and 5 CHEOPS occultations extracted from the phase curves (in red). The eclipse depth derived from the joint CHEOPS phase curve fit is shown in green (367 ± 17 ppm, with $1\sigma$ shaded) and the mean value of these 14 occultations (320 ± 11 ppm) is shown in orange (with $1\sigma$ shaded). These two mean eclipse depths are consistent with each other by just over $2\sigma$. The right-hand y-axis shows how the eclipse depths convert to dayside brightness temperatures, assuming the dayside hemisphere is a blackbody and radiates at a uniform temperature.

<figure><img src="image_13.png" /></figure>

**Fig. 9.** Heat redistribution efficiency of hot Jupiters plotted against their equilibrium temperatures. The coloured points show the KELT-9b efficiency, $\epsilon$, derived in this study.

by $3.5\sigma$. Finally, our best-fit orbital inclination ($i = 86.38^{+0.29}_{-0.30}$°) is just over $1.5\sigma$ away from the inclination reported in W20 and $1.6\sigma$ away from Ahlers et al. (2020).

We find that for the TESS phase curves, the eclipse depth reported in this paper from the phase curve fit ($645 \pm 15$ ppm) is consistent with the eclipse depth found in W20 ($650^{+14}_{-15}$ ppm). As discussed in Sect. 4.2, from our fitted eclipse depth and nightside flux we derived a TESS dayside brightness temperature of $4643 \pm 26$ K which is also consistent with W20 (4600 ± 100 K). Our TESS nightside brightness temperature of $3177^{+60}_{-61}$ K is consistent within $1.2\sigma$ of W20's nightside temperature (3040 ± 100 K). Following from this, in this work we have used the method in Morris et al. (2022) to derive the Bond albedo and heat redistribution efficiency. Our values differ with the Bond albedo in W20 by $2.3\sigma$ and the heat redistribution parameter is consistent within $1\sigma$. The gravity-darkening parameters will be discussed in Sect. 4.5.3.

### 4.5.2. Spitzer

For the orbital and system parameters, the planetary radius and transit duration differ with the same parameters reported in

Mansfield et al. (2020, hereafter M20) by $2.1\sigma$ and $6\sigma$ respectively, indicative of the different reduction methods influencing the derived parameters.

For the phase curve parameters, the Spitzer hotspot offset reported in this paper ($-18.2^{+1.7}_{-1.6}$° eastwards) is consistent with the offset reported in M20 ($-18.7^{+2.1}_{-2.3}$°). However as mentioned in Sect. 4.2, our nightside brightness temperature is inconsistent with the nightside brightness temperature reported in M20 by $3.5\sigma$. The dayside brightness temperature reported in this paper is consistent with M20's value by just under $2\sigma$.

### 4.5.3. Gravity-darkened transits

We find a sky-projected spin orbit angle of $-85.67^{+0.46}_{-0.45}$ degrees which is consistent with the spin orbit angle reported in Ahlers et al. (2020, hereafter A20), and Gaudi et al. (2017). Our value of the gravity-darkening coefficient, $\beta$, ($0.2270^{+0.0077}_{-0.0080}$) is not consistent with A20, but is consistent with the Claret (2016) theoretical value. Our best-fit stellar inclination, $i^* = 47.1 \pm 1.1$° (related to the stellar obliquity $\phi$, specified in A20, by $i^* = 90° - \phi$), and stellar rotation period ($P_{rot} = 18.96 \pm 0.34$ h) are consistent with the value reported in A20.

Using the fitted values for $i^*$, $i$ and $\lambda$ we used Eq. (1) to calculate a value for the true spin orbit angle of $84.36^{+0.37}_{-0.38}$°. This is within $1\sigma$ of the value A20 reported.

## 5. Discussion

### 5.1. Challenges of simulating UHJs

Tidally locked, highly irradiated exoplanets, including UHJs, are complex, 3D objects. General circulation models (GCMs), which are numerical solvers of the 3D fluid equations, have been adapted to study hot Jupiters (see, e.g.: Showman et al. 2009; Rauscher & Menou 2010; Heng et al. 2011; Kataria et al. 2013; Mayne et al. 2014). Recently, GCMs have been used to study UHJs (e.g.: Tan & Komacek 2019).

The higher temperatures of UHJs present additional technical challenges for GCMs. As we have shown in Fig. 4, the atmosphere of KELT-9b transitions from temperatures where it is dominated by atomic hydrogen on its dayside to being dominated by molecular hydrogen on its nightside, verifying
the prediction of Bell & Cowan (2018). Since the specific heat capacity at constant pressure, $c_p$, varies as the reciprocal of the mean molecular mass $m$, i.e., $c_p \propto 1/m$, it implies that $c_p$ changes by a factor of 2 within the atmosphere of an UHJ, unlike for a regular hot Jupiter where it is roughly constant. The transformation from atomic to molecular hydrogen (and back) means that an additional cooling/heating term needs to be inserted into the governing equations (Bell & Cowan 2018; Tan & Komacek 2019). KELT-9b is an extreme example of these processes occurring in UHJs and correctly reproducing the observed dayside-to-nightside flux redistribution will require them to be simulated correctly. The CHEOPS, TESS, and *Spitzer* phase curves presented in the current study, as well as their associated temperature maps, will provide valuable constraints for future GCM studies of KELT-9b.

Another important constraint provided by the current work is that the climate of KELT-9b is somewhat stable: over around 270 orbital periods, the globally averaged dayside temperature varies by less than 50 K (see Sect. 4.3 on results of occultation analysis). These results are consistent with theoretical expectations produced from GCM analysis in works such as Showman et al. (2009) and Komacek & Showman (2020). Simulating variability accurately using GCMs is challenging, as it depends on the choice of governing equations (e.g.: Cho et al. 2008), numerical dissipation (e.g.: Heng et al. 2011; Thrastarson & Cho 2011) and choice of bottom boundary condition (e.g.: Liu & Showman 2013). In particular, our inability to specify numerical dissipation, which is often required to numerically stabilise GCM runs, from first principles (Heng et al. 2011) implies that energy and momentum conservation, and therefore our ability to accurately predict wind speeds and variability, is limited (Goodman 2009).

In these ways, observations of UHJs lead our current ability to simulate them, with KELT-9b presenting the most extreme case study. The empirical constraints derived here therefore provide important checks on future GCMs of UHJs.

## 5.2. Anticipating JWST multi-wavelength phase curves

Formally, the spherical and Bond albedos are monochromatic and bolometric quantities, respectively. The spherical albedo is the monochromatic version of the Bond albedo (monochromatic version of Eq. (5)). One of the key limitations of the current study is that the measured CHEOPS, TESS, and *Spitzer* phase curves are neither monochromatic nor bolometric. Essentially, our ability to extract temperatures and Bond albedos is based on an extrapolation: the assumption that the spectral energy distribution (SED) of KELT-9b follows a perfect blackbody. Generally, a non-blackbody SED sampled in the CHEOPS versus *Spitzer* bandpasses will return different brightness temperatures. These differences translate into differences in the inferred Bond albedos that we report in Table 6.

In the era of the James Webb Space Telescope (JWST), multi-wavelength phase curves, which will be monochromatic to a good approximation, will become the norm. Two questions concerning the data analysis and interpretation of JWST thermal phase curves arise.

### 5.2.1. May we still use box models to extract Bond albedos from JWST thermal phase curves?

Generally, the heating of an atmosphere by starlight occurs in the near-ultraviolet to optical range of wavelengths, which is processed and re-emitted in the infrared range of wavelengths

<table><caption>**Table 6.** Bond albedo and heat redistribution in the three bandpasses.</caption><thead><tr><th>Bandpass</th><th>A<sub>B</sub></th><th>ε</th></tr></thead><tbody><tr><td>CHEOPS</td><td>-0.168<sup>+0.096</sup><sub>-0.092</sub></td><td>0.314<sup>+0.051</sup><sub>-0.058</sub></td></tr><tr><td>TESS</td><td>-0.025 ± 0.032</td><td>0.342<sup>+0.016</sup><sub>-0.017</sub></td></tr><tr><td>*Spitzer*</td><td>-0.109 ± 0.068</td><td>0.257 ± 0.013</td></tr></tbody></table>

as thermal emission. (As already mentioned in the introduction, KELT-9b is a special case where its dayside thermal emission radiates in the optical.) While it is possible to quantify the wavelength-dependent flux of starlight incident upon the top of the atmosphere, it is much more challenging to describe how much starlight penetrates to each atmospheric layer as a function of wavelength. Any such attempt will be model-dependent and probabilistic (in a Bayesian sense).

If one integrates over all wavelengths, the bolometric flux of starlight is simply given by the Stefan-Boltzmann law. If the set of JWST thermal phase curves covers the entire wavelength range of the SED, then one may empirically derive the bolometric thermal flux emitted by the atmosphere. The Bond albedo may then be inferred without any assumption on the nature of the SED of the exoplanetary atmosphere, i.e., it is not necessary to assume a blackbody SED. Essentially, one bypasses the need for a 0D box model.

Alternatively, one could still perform an analysis like we have done in this paper and obtain different Bond albedos for every JWST thermal phase curve, assuming a blackbody SED for each surface element of the planet. However, as previously explained, this assumption will not be perfect for a planet with non-blackbody behaviour. Therefore the former approach, assuming the JWST phase curves cover the entire planetary SED, would be much more accurate to determine the total bolometric flux of the planet, and therefore a more accurate (single value) of the Bond albedo.

### 5.2.2. Is it possible to extract spherical albedos from JWST thermal phase curves?

As mentioned, the spherical albedo is the monochromatic version of the Bond albedo. At first thought, if JWST can produce monochromatic phase curves, it is easy to suggest this could translate into a set of spherical albedos for the planet. However this would be incorrect. In principle, if one knew exactly how much starlight was deposited in a specific atmospheric layer, then one could compare that to the thermal emission from that layer and derive the spherical albedo. However, as previously described this would be a model-dependent exercise that involves some assumption on the model atmosphere in order to perform radiative transfer from the top of the atmosphere to the layer in question.

Specifically, Eq. (5) describes our approach for deriving the Bond albedo from a 2D temperature map. Despite the CHEOPS, TESS, and *Spitzer* thermal phase curves being neither monochromatic nor bolometric, this equation allows one to work with temperatures because of the assumption of a blackbody SED. When JWST multi-wavelength thermal phase curves are available, the numerator of Eq. (5) may be generalised such that $F_p$ is empirically derived from the data with no need to assume a blackbody SED. However, the denominator cannot be generalised in a straightforward way as $\sigma T_*^4$ needs to be replaced by the flux of starlight deposited in the same
atmospheric layer. In principle, it is possible to solve for these
wavelength-dependent fluxes within a holistic framework that
simultaneously interprets phase-dependent emission spectra and
wavelength-dependent phase curves. Such a framework would
have to account for scattered starlight versus thermal emission
as functions of wavelength. However, as explained in the pre-
vious subsection, calculating the Bond albedo would still be
possible from the JWST thermal phase curves alone, provided
they sufficiently spanned the planetary SED. As the Bond albedo
describes the planet's input and output energy bolometrically, it
removes the need to undertake the difficult task of probing the
individual atmospheric layers.

**6. Conclusion**

In this work, we have simultaneously analysed CHEOPS phase curves as well as public phase curves from TESS and Spitzer to infer joint constraints on the phase curve variation, gravity- darkened transits and occultation depth in the three bandpasses (Fig. 2). From this analysis we find the following results:

- We derive 2D temperature maps of the atmosphere at three different depths, and calculate dayside and nightside brightness temperatures of the planet in each bandpass (Fig. 4, Table 5 and Sects. 4.1 and 4.2);

- The day-night heat redistribution of ~0.3 confirms theoretical expectations of enhanced energy transfer to the planetary nightside due to dissociation and recombination of molecular hydrogen in ultra-hot Jupiters (Fig. 9, Table 6 and Sect. 4.4);

- We also find a Bond albedo consistent with zero (Table 6 and Sect. 4.4);

- We also analyse 9 CHEOPS occultations of KELT-9 and find no evidence of variability of the brightness temperature of the planet, excluding variability greater than 1% ($1\sigma$) (Fig. 8 and Sect. 4.3).

Acknowledgements. We gratefully acknowledge the open source software which made this work possible: *astroty (Astropy Collaboration et al. 2013, 2018), batman (Kreidberg 2015), ipython (Pérez & Granger 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013), PyMC3 (Salvatier et al. 2016), pytransit (Parviainen 2015). This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AYA2017-84089. 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. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-190080). 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. S.H. gratefully acknowledges CNES funding through the grant 837319. M.L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. 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 Internacionalizaçã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. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). A.C.C. and T.W. acknowledge support from STFC consolidated grant numbers ST/R000824/1 and

ST/V000861/1, and UKSA grant number ST/R003203/1. L.Bo., G. Br., V.Na.,
I.Pa., G.Pi., R.Ra., G.Sc., V.Si., and T.Zi. acknowledge support from CHEOPS
ASI-INAF agreement no. 2019-29-HH.0. Y.A. and M.J.H. acknowledge the
support of the Swiss National Fund under grant 200020_172746. We acknowl-
edge 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
F.R.S.-FNRS Postdoctoral Researcher. M.F. and C.M.P. gratefully acknowledge
the support of the Swedish National Space Agency (DNR 65/19, 174/18). D.G.
gratefully acknowledges financial support from the CRT foundation under Grant
No. 2018.2323 "Gaseous rocky? Unveiling the nature of small worlds". M.G.
is an F.R.S.-FNRS Senior Research Associate. 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 Ile-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. 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. GyM.Sz. acknowledges the support of the
Hungarian National Research, Development and Innovation Office (NKFIH)
grant K-125015, a PRODEX Experiment Agreement No. 4000137122, 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 F.R.S.-FNRS Research
Associate. N.A.W. acknowledges UKSA grant ST/R004838/1.

**References**

Agol, E., Cowan, N. B., Knutson, H. A., et al. 2010, ApJ, 721, 1861

Ahlers, J. P., Johnson, M. C., Stassun, K. G., et al. 2020, AJ, 160, 4

Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30

Arcangeli, J., Désert, J.-M., Parmentier, V., et al. 2019, A&A, 625, A136

Armstrong, D. J., de Mooij, E., Barstow, J., et al. 2016, Nat. Astron., 1, 0004

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Barnes, J. W. 2009, ApJ, 705, 683

Bell, T. J., & Cowan, N. B. 2018, ApJ, 857, L20

Bello-Arufe, A., Buchhave, L. A., Mendonça, J. M., et al. 2022, A&A, 662, A51

Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109

Blackwell, D. E., & Shallis, M. J. 1977, MNRAS, 180, 177

Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18

Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5

Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157

Borsa, F., Rainer, M., Bonomo, A. S., et al. 2019, A&A, 631, A34

Castelli, F., & Kurucz, R. L. 2003, in IAU Symposium, Modelling of Stellar Atmospheres, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, 210, A20

Cho, J. Y. K., Menou, K., Hansen, B. M. S., & Seager, S. 2008, ApJ, 675, 817

Claret, A. 2016, A&A, 588, A15

Claret, A. 2017, A&A, 600, A30

Claret, A. 2021, VizieR Online Data Catalog: J/other/RNAAS/5

Claret, A., & Bloemen, S. 2011, A&A, 529, A75

Cowan, N. B., & Agol, E. 2011, ApJ, 729, 54

Deline, A., Hooton, M. J., Lendl, M., et al. 2022, A&A, 659, A74

Demory, B.-O., Gillon, M., de Wit, J., et al. 2016a, Nature, 532, 207

Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016b, MNRAS, 455, 2018

Foreman-Mackey, D. 2018, RNAAS, 2, 31
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Gai, A. D., & Knuth, K. H. 2018, ApJ, 853, 49

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514

Gillon, M., Demory, B. O., Benneke, B., et al. 2012, A&A, 539, A28

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456

Goodman, J. 2009, ApJ, 693, 1645

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Heng, K., & Showman, A. P. 2015, Annu. Rev. Earth Planet. Sci., 43, 509

Heng, K., & Workman, J. 2014, ApJS, 213, 27

Heng, K., Menou, K., & Phillipps, P. J. 2011, MNRAS, 413, 2380

Heng, K., Lyons, J. R., & Tsai, S.-M. 2016, ApJ, 816, 96

Hoeijmakers, H. J., Snellen, I. A. G., & van Terwisga, S. E. 2018, A&A, 610, A47

Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165

Hooton, M. J., Watson, C. A., de Mooij, E. J. W., Gibson, N. P., & Kitzmann, D. 2018, ApJ, 869, L25

Hooton, M. J., Hoyer, S., Kitzmann, D., et al. 2022, A&A, 658, A75

Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2016, AJ, 152, 44

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E

Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, ApJ, 767, 76

Keating, D., Cowan, N. B., & Dang, L. 2019, Nat. Astron., 3, 1092

Kilpatrick, B. M., Kataria, T., Lewis, N. K., et al. 2020, AJ, 159, 51

Kipping, D. M. 2013, MNRAS, 435, 2152

Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183

Komacek, T. D., & Showman, A. P. 2016, ApJ, 821, 16

Komacek, T. D., & Showman, A. P. 2020, ApJ, 888, 2

Kreidberg, L. 2015, PASP, 127, 1161

Kreidberg, L., Line, M. R., Parmentier, V., et al. 2018, AJ, 156, 17

Lally, M., & Vanderburg, A. 2022, AJ, 163, 181

Lanotte, A. A., Gillon, M., Demory, B. O., et al. 2014, A&A, 572, A73

Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94

Lightkurve Collaboration (Cardoso, J. V. d. M. a., et al.) 2018, Astrophysics Source Code Library [record ascl:1812.0103]

Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4

Liu, B., & Showman, A. P. 2013, ApJ, 770, 42

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27

Mansfield, M., Bean, J. L., Stevenson, K. B., et al. 2020, ApJ, 888, L15

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77

Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014, A&A, 561, A1

Mendonça, J. M., Malik, M., Demory, B.-O., & Heng, K. 2018, AJ, 155, 150

Morris, B. M., Delrez, L., Brandeker, A., et al. 2021, A&A, 653, A173

Morris, B. M., Heng, K., Jones, K., et al. 2022, A&A, 660, A123

Owens, N., de Mooij, E. J. W., Watson, C. A., & Hooton, M. J. 2021, MNRAS, 503, L38

Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110

Parviainen, H. 2015, MNRAS, 450, 3233

Peixoto, J. P., & Oort, A. H. 1984, Rev. Mod. Phys., 56, 365

Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21

Pino, L., Désert, J.-M., Brogi, M., et al. 2020, ApJ, 894, L27

Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320

Salmon, S. J. A. J., Van Grootel, V., Buldgen, G., Dupret, M. A., & Eggenberger, P. 2021, A&A, 646, A7

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, Astrophysics Source Code Library [record ascl:1610.016]

Schanche, N., Hébrard, G., Collier Cameron, A., et al. 2020, MNRAS, 499, 428

Schwartz, J. C., Kashner, Z., Jovmir, D., & Cowan, N. B. 2017, ApJ, 850, 154

Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83

Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stevenson, K. B., Harrington, J., Fortney, J. J., et al. 2012, ApJ, 754, 136

Sudarsky, D., Burrows, A., & Pinto, P. 2000, ApJ, 538, 885

Tan, X., & Komacek, T. D. 2019, ApJ, 886, 26

Thrastarson, H. T., & Cho, J. Y.-K. 2011, ApJ, 729, 117

Vehtari, A., Gelman, A., & Gabry, J. 2015, arXiv e-prints, [arXiv:1507.04542]

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

von Zeipel, H. 1924, MNRAS, 84, 665

Wallack, N. L., Knutson, H. A., & Deming, D. 2021, AJ, 162, 36

Welsh, W. F., Orosz, J. A., Seager, S., et al. 2010, ApJ, 713, L145

Wong, I., Shporer, A., Kitzmann, D., et al. 2020, AJ, 160, 88

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Wyttenbach, A., Mollière, P., Ehrenreich, D., et al. 2020, A&A, 638, A87

Yan, F., & Henning, T. 2018, Nat. Astron., 2, 714

1 Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland
e-mail: kathryn.jones@unibe.ch

2 Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

3 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK

4 Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

5 Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland

6 Centre Vie dans l'Univers, Faculté des sciences, Université de Genève, Quai Ernest-Ansermet 30, 1211 Genève 4, Switzerland

7 Aix Marseille Univ., CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France

8 Instituto de Astrofísica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

9 Departamento de Física e Astronomia, Faculdade de Ciencias, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal

10 Space Research Institute, Austrian Academy of Sciences, Schmiedelstrasse 6, 8042 Graz, Austria

11 Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK

12 Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

13 Instituto de Astrofísica de Canarias, 38200 La Laguna, Tenerife, Spain

14 INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy

15 Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

16 Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, 08193 Bellaterra, Spain

17 Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

18 Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary

19 Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

20 Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

21 Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden

22 Université de Paris, Institut de physique du globe de Paris, CNRS, 75005 Paris, France

23 Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

24 Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

25 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

26 Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

27 Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden

28 Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, 10125, Torino, Italy

29 University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria

30 Science and Operations Department - Science Division (SCI-SC), Directorate of Science, European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201-AZ Noordwijk, The Netherlands
31 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15–17, Hungary

32 ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

33 IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France

34 Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis bd Arago, 75014 Paris, France

35 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

36 Astrophysics Group, Keele University, Staffordshire ST5 5BG, UK

37 Department of Astrophysics, University of Vienna, Tuerkenschanzstrasse 17, 1180 Vienna, Austria

38 Institute of Optical Sensor Systems, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

39 Dipartimento di Fisica e Astronomia "Galileo Galilei", Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

40 ETH Zurich, Department of Physics, Wolfgang-Pauli-Strasse 2, 8093 Zurich, Switzerland

41 ESTEC, European Space Agency, 2201AZ Noordwijk, NL

42 Zentrum für Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin, Germany

43 Institut für Geologische Wissenschaften, Freie Universität Berlin, 12249 Berlin, Germany

44 ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary

45 MTA-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary

46 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
**Appendix A: Transit best-fits**

Figure A.1 shows the transit models and photometry.

<figure><img src="image_14.png" /><figcaption>Fig. A.1: CHEOPS, TESS, and *Spitzer* phase-folded and detrended transits overplotted with the best-fit transit model (in red), with all other models and systematics removed. In black are the binned grey datapoints, with error bars that are smaller than the point size in all panels so they are not visible.</figcaption></figure>
**Appendix B: Posterior distributions**

Figure B.1 shows the posterior distributions for the phase curve parameters.

<figure><img src="image_15.png" /><figcaption>Fig. B.1: Corner plot showing posteriors of the phase curve parameters in the joint fit. The three contour lines in each subplot refer to the 1-, 2- and 3-sigma contour levels.</figcaption></figure>