© K. Jones et al. 2022

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, different fractions of the dayside of the exoplanet will be visible to us at different times in its orbit. There are three main observables 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 single 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 ~2500K. 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 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 ± 180 K (Gaudi et al. 2017), It has a near-polar, 1.48-day orbit around a rapidly rotating A0/B9 star (T_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σ) 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 <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 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¹. 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 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_BP, and G_RP, 2MASS J, H, and K, and WISE W1 and W2 (Gaia Collaboration 2021; Skrutskie et al. 2006; Wright et al. 2010) 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_* = 2.379 ± 0.038 R_⊙.

The stellar radius R* together with the effective temperature T_eff and the metallicity [Fe/H] constitute the input set to determine the isochronal age t_* and mass M_*. To make the derivation process more robust we employed two different stellar evolutionary models, namely PARSEC v1.2S² (Marigo et al. 2017) and CLES (Code Liègeois d’Evolution 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_* = 0.3 ± 0.1 Gyr and M_⊙.

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 know exoplanet with a dayside temperature of ~4600K 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 80000 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_*);

• 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_rot, the rotation period of the star;

• T_pole, 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_pole = 10170K, e = 0 and ω = 90°. We also fixed R_* = 2.379 R_⊙ (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 (Gaudi 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: (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): (2)

where θ and ϕ are latitude and longitude. T is a background temperature upon which 2D perturbations exist and are quantified by the “alphabet” or basis functions h_mℓ as defined by Eq. (258) of Heng & Workman (2014). These functions are also dependent on the dimensionless parameters α and ω_drag. 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 ω_drag 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 ω_drag ≳ 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;

• ω_drag, dimensionless drag frequency;

• C_mℓ, the power of individual harmonic modes;

• ℓ_max, describes the highest spherical harmonic mode in the model;

• planetary parameters including a, R_p, T_pole;

• telescope bandpass transmission efficiency;

• 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: (3)

where I, the intensity, is defined by (4)

where F_λ is the instrument-specific filter response function and B_λ(T(θ, ϕ)) 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, ω_drag = 4.5 and ℓ_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_mℓ = 0 apart from C₁₁ 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 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.

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.

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σ 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 ± 6 ppm in Wong et al. (2020), therefore this fitted amplitude is over 5σ 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 (Salvatier 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σ), 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 occulations. 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, σ) ppm, and truncated at zero so that the amplitude is always positive, and the standard deviation σ 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₀, 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*;

• ρ, 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σ and the rest was within 2σ. 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: ∆ϕ (hotspot offset), C_1,1 (C_mℓ power coefficient) and 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.

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.

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.

Fig. 3 Posteriors of the phase curve parameters across the three band-passes showing 1σ-confidence contours. Red is CHEOPS, green is TESS and blue is Spitzer.

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 (“ID”) 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 K respectively. We found the nightside brightness temperatures of the three bandpasses (in the same order) to be K, K and 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σ to the dayside brightness temperature reported in Mansfield et al. (2020), K. However our nightside temperature is not consistent with their nightside brightness temperature of K (difference of over 3σ). 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 . 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 HOOK, 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σ 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 fitted 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 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).

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.

Fig. 5 Posteriors of the brightness day and nightside temperatures from the different bandpasses.

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

where σ is the Stefan-Boltzmann constant and F_p(θ, ϕ) 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 . 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): (6)

Figure 9 shows the heat redistribution (ε) 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 ε. However in Bell & Cowan (2018), they predict a rising heat redistribution for UHJ due to the H₂ dissociation and recombination increasing the heat transport around the planet. Our work agrees with this theory and in Fig. 9, it appears that ε does indeed fall with planet equilibrium temperatures up to around 2500 K, however after that the ε rises again with temperature and KELT-9b’s ε calculated in this paper supports this trend.

Fig. 6 Nightside integrated temperatures of hot lupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.

Fig. 7 Dayside integrated temperatures of hot lupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.

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 ± 0.97° is inconsistent with the offset reported in Wong et al. (2020, hereafter W20) (5.2 ± 0.9°) by over 5σ (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 ± 0.000099 R_*) in the TESS bandpass is consistent with W20.

For some of the global parameters, our best-fit impact parameter differs from the W20 value by 1.4σ. Our best-fit semi-major axis is different from the W20 value by around 4σ and our best-fit period (P = 1.48111949 ± 0.00000034 days) is different to W20’s value by 3.5σ. Finally, our best-fit orbital inclination is just over 1.5σ away from the inclination reported in W20 and 1.6σ 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 ± 15 ppm) is consistent with the eclipse depth found in W20 . As discussed in Sect. 4.2, from our fitted eclipse depth and nightside flux we derived a TESS dayside brightness temperature of 4643 ± 26 K which is also consistent with W20 (4600 ± 100 K). Our TESS nightside brightness temperature of K is consistent within 1.2σ 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.3cr and the heat redistribution parameter is consistent within 1σ. The gravity-darkening parameters will be discussed in Sect. 4.5.3.

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σ shaded) and the mean value of these 14 occultations (320 ± 11 ppm) is shown in orange (with 1σ shaded). These two mean eclipse depths are consistent with each other by just over 2σ. 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.

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

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σ and 6σ respectively, indicative of the different reduction methods influencing the derived parameters.

For the phase curve parameters, the Spitzer hotspot offset reported in this paper ( eastwards) is consistent the offset reported in M20 . However as mentioned in Sect. 4.2, our nightside brightness temperature is inconsistent with the nightside brightness temperature reported in M20 by 3.5σ. The dayside brightness temperature reported in this paper is consistent with M20’s value by just under 2σ.

4.5.3 Gravity-darkened transits

We find a sky-projected spin orbit angle of 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, β, is not consistent with A20, but is consistent with the Claret (2016) theoretical value. Our best-fit stellar inclination, i*= 47.1 ± 1.1° (related to the stellar obliquity ϕ, specified in A20, by i* = 90° − ϕ), and stellar rotation period (P_rot = 18.96 ± 0.34 h) are consistent with the value reported in A20.

Using the fitted values for i*, i and λ we used Eq. (1) to calculate a value for the true spin orbit angle of . This is within ler 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 ∝ 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 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 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 previous 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σ) (Fig. 8 and Sect. 4.3).

Acknowledgements

We gratefully acknowledge the open source software which made this work possible: astropy (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 Counci (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 Internacionalizacão by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987, O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. 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 acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiología (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract no. 4000124370. S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. X.B., S.C., D.G., M.F. and J.L. acknowledge their role as ESA-appointed CHEOPS science team members. A.Br was supported by the SNSA. A.C.C. acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L.D. is an 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 “Gaseousor 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 Region 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 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.

Appendix A Transit best-fits

Figure A.1 shows the transit models and photometry.

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.

Appendix B Posterior distributions

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

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.

References

All Tables

All Figures

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.

In the text

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

	Fig. 3 Posteriors of the phase curve parameters across the three band-passes showing 1σ-confidence contours. Red is CHEOPS, green is TESS and blue is Spitzer.
In the text

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.

In the text

	Fig. 5 Posteriors of the brightness day and nightside temperatures from the different bandpasses.
In the text

	Fig. 6 Nightside integrated temperatures of hot lupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.
In the text

	Fig. 7 Dayside integrated temperatures of hot lupiters plotted against their equilibrium temperatures (Keating et al. 2019). The coloured points show the KELT-9b temperatures derived in this study.
In the text

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σ shaded) and the mean value of these 14 occultations (320 ± 11 ppm) is shown in orange (with 1σ shaded). These two mean eclipse depths are consistent with each other by just over 2σ. 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.

In the text

	Fig. 9 Heat redistribution efficiency of hot Jupiters plotted against their equilibrium temperatures. The coloured points show the KELT-9b efficiency, ε, derived in this study.
In the text

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

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