© The Authors 2025

1 Introduction

Most exoplanetary properties, such as temperatures, radii, compositions, and orbital architectures, span a broad range of values that go well beyond those of our Solar System. One of the categories of extra-solar planets that best illustrates this variety are hot Jupiters. These gas giants orbiting close to their stars exhibit extreme conditions in their massive atmospheres with temperatures above ~1000 K. At such small planet-to-star separations, hot Jupiters tend to synchronise their rotation and revolution periods due to immense tidal forces, which results in the stellar insolation always heating the same planetary hemisphere. Tidally locked objects receive large amounts of energy on their permanent dayside, which may lead to high day-to-night contrast and/or strong winds redistributing heat on a planetary scale (Showman & Guillot 2002; Showman et al. 2010; Komacek & Showman 2016; Zhang 2020). The hottest planets in this category are dubbed ultra-hot Jupiters (UHJs), whose atmospheres enter a specific regime where the dayside is dominated by atomic hydrogen (Bell & Cowan 2018), the major source of spectral continuum opacity is induced by hydrogen anions (H⁻) (Arcangeli et al. 2018; Lothringer et al. 2018; Parmentier et al. 2018), and the atmospheric composition of atoms and ions resembles that of a star (Kitzmann et al. 2018). The thermal emission coming from UHJs is so strong that it can be observed from the infrared up to the optical with large occultation signals when the planetary dayside is hidden by the host star. Observations over a large range of wavelengths are thus decisive to get key insights on UHJs.

The UHJ known as WASP-18 b (Hellier et al. 2009) is orbiting an early F6-type star with a very short orbital period of 0.94 day and a separation of 0.02 au (semi-major axis). The planet is a peculiar object of ten Jupiter masses, which is only slightly larger than Jupiter (1.2 R _J ), which makes it denser than the Earth (more than five times denser than Jupiter). With such properties, WASP-18 b resides at the intersection between brown-dwarf and planet populations. It was recently observed by the Characterising Exoplanet Satellite (CHEOPS; Benz et al. 2021) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). These observations continue to contribute to the extensive characterisation of the planet with several instruments, including the Spitzer (Werner et al. 2004) and James Webb space telescopes (JWST; Gardner et al. 2006). Recent works have reported the detection of several species in the atmosphere of WASP-18 b including H⁻, H₂O, OH, and CO (Changeat et al. 2022; Brogi et al. 2023; Yan et al. 2023; Coulombe et al. 2023). Upper limits on the geometric albedo have been computed in the TESS passband (Shporer et al. 2019; Blažek et al. 2022).

In this work, we report the new CHEOPS and TESS phasecurve observations of WASP-18 b. We jointly analysed these data sets with already published TESS observations (Shporer et al. 2019). We also included Spitzer observations: 4 occultations previously analysed in Nymeyer et al. (2011) and 10 occultations observed in 2015 and firstly reported in Deming et al. (2023). We first start by presenting refined stellar properties of WASP-18 in Section 2. Sections 3 and 4 detail the data sets and the modelling framework used in this analysis. We present our results and how we improved the constraints on the planetary parameters (Section 5). Finally, we discuss the implications on the atmospheric properties of the planet in Section 6.

2 Stellar properties of WASP-18

In the framework of this study, we derived the properties of the host star WASP-18 (HD 10069; TOI-185) listed in Table 1, following the methods described below.

We used co-added high-resolution spectra obtained with the HARPS spectrograph (Mayor et al. 2003; Pepe et al. 2004) and analysed them with the spectral synthesis tool SME ¹ (Spectroscopy Made Easy; Valenti & Piskunov 1996; Piskunov & Valenti 2017). This software fits the observations to a set of computed synthetic spectra for a chosen set of parameters based on atomic and molecular line data from the Vienna Atomic Line Database² (VALD; Piskunov et al. 1995; Ryabchikova et al. 2015). We chose the ATLAS12 (Kurucz 2013) stellar atmosphere grid and modelled the stellar effective temperature, T _eff, the surface gravity, log g, abundances, and the projected rotational velocity, v sin i_★ . We used narrow and unblended lines between 6200 Â and 6600 Â to derive abundances and v sin i_★ . The micro- (Bruntt et al. 2008) and macro-turbulent (Doyle et al. 2014) velocities were kept fixed to 1.50 km.s⁻¹ and 5.60km.s⁻¹, respectively. We refer to Persson et al. (2018) for further details on the modelling. Our results suggest that WASP-18 is an early F6V star with T _eff = 6332 ± 60K. All results are listed in Table 1 and are within 1 σ agreement with the Gaia DR3 results and the values listed on the NASA archive³.

We used our stellar spectroscopic parameters to determine the stellar radius of WASP-18 using a MCMC modified infrared flux method (IRFM; Blackwell & Shallis 1977; Schanche et al. 2020). From the stellar properties T _eff, log g, and [Fe/H], we constrained stellar atmospheric models from two catalogues (Kurucz 1993; Castelli & Kurucz 2003). We then constructed spectral energy distributions (SEDs) from which we computed synthetic photometry. We compared the photometry to the broadband observations in the following passbands: Gaia G, G _BP, and G _RP, 2MASS J, H, and K, and WISE W 1 and W2 (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2023). We derived the stellar bolometric flux that we had converted into the stellar angular diameter using the measured effective temperature. The angular diameter and offset-corrected Gaia parallax (Lindegren et al. 2021) were combined to produce the stellar radius. To account and correct for atmospheric model uncertainties, we conducted a Bayesian modelling, averaging of the posterior distributions of the radius produced via this process with both atmospheric catalogues.

The effective temperature, T _eff, metallicity, [Fe/H], and radius, R_★ along with their uncertainties, were used as basic inputs to infer the isochronal mass, M_★ , and age, t_★ , from stellar evolutionary models. We derived two pairs {M_★ , t_★ } from two different approaches. The first pair of mass and age were estimated by applying the isochrone placement algorithm (Bonfanti et al. 2015, 2016) that interpolates the input stellar parameters within pre-computed grids of isochrones and tracks of PARSEC ⁴ v1.2S (Marigo et al. 2017). As outlined in Bonfanti et al. (2016), the isochrone placement code also implements the gyrochrono-logical relation by Barnes (2010) and, thus, we further provided the algorithm with the projected equatorial velocity of the star (v _★ sini_★ ) to benefit from the synergy between isochrone interpolation and gyrochronology (see e.g. Angus et al. 2019), thereby improving the convergence. The second pair of mass and age estimates were derived from the Code Liègeois d’Évolution Stellaire (CLÉS; Scuflaire et al. 2008), which computes the best-fit evolutionary track accounting for the basic input set of stellar parameters following a Levenberg-Marquardt minimisation scheme (Salmon et al. 2021). As detailed in Bonfanti et al. (2021), we checked the mutual consistency between the two respective pairs of outcomes via a χ² -based criterion and then merged the results obtaining $M_{\star }={1.245}_{-0.058}^{+0.057}{M}_{\odot }$ M_⊙ and t _★ = 1.4 ± 0.7 Gyr.

We also note that the star WASP-18 has a very low magnetic activity as indicated by several indices (e.g. Fig. 2 of Fossati et al. 2013). This unusually quiet behaviour and its possible causes have been discussed in several works (e.g. Lanza 2014; Fossati et al. 2015; Lanza & Breton 2024).

3 Observations and data reduction

3.1 CHEOPS observations

We obtained 38 observations of WASP-18 with CHEOPS spanning three observability seasons and a total time range of 2 years and 2 months. These visits were part of the Guaranteed Time Observation programmes of the CHEOPS consortium (CH_PR100012 and CH_PR100016) and covered 25 occultations, 11 transits, and a complete phase-curve. One of the visits, executed on 2022-09-29, was interrupted to perform a collision avoidance manoeuvre and the very short data set (of 1 hour) was not included in the analyses of this work. The visits from programme CH_PR100012 (transits) were observed with an exposure time of 53 sec, and the visits from programme CH_PR100016 (occultations and phase curve) were observed with an exposure time of 50 s. Table A.1 provides a detailed overview of all the CHEOPS observations with their corresponding raw light curves displayed in Fig. A.1.

The photometric light curve was extracted from the raw data using the Data Reduction Pipeline (DRP; Hoyer et al. 2020) that performs circular aperture photometry on images corrected for several effects of instrumental (bias offset, gain, dark current, hot pixels, and flat field) and astrophysical (cosmic rays, stray light, and smearing trails) origins. The DRP provides photometry for a range of aperture radii to make it possible to select the best option. To determine the aperture size with the best photometric precision, we first estimated the white noise level from the raw light curves. For each aperture, we computed the point-to-point difference of the light curves to remove any signal (planetary or red noise) and calculated the standard deviations of the resulting flat data sets using a method robust to outliers. These values were then divided by $\sqrt{2}$ to correct for the point-to-point difference step and obtain a reliable estimate of the white noise level. We then selected the aperture of 24 pixels that consistently gives the smallest noise level over all the 37 visits.

Across all the visits, we identified 400 out of 13 760 points (2.91%) that were flagged by the DRP (e.g. out-of-range temperatures, Earth occultation, and high cosmic ray hits). We also detected 148 outliers (1.11%) by performing 3.5-σ clipping after subtracting the best-fit results of the modelling described in Section 4. Finally, we flagged 407 (3.05%) data points that had background values above 7.8 10⁵ electrons. This threshold was fixed as the point above which the flux does not correlate monotonically with the background level (a similar approach as that shown in Fig. 2 of Deline et al. 2022). The total of discarded points was 509 out of 13 360 (3.81%) since some of the outliers also had high background values.

Due to the orbital configuration of the spacecraft and the requirements on its thermal stability, CHEOPS photometry is generally affected by short interruptions every ~100 min (the orbital period of CHEOPS) mostly due to Earth occultations and also systematic trends induced by background level variations or close-by stars. These systematics usually strongly correlate with the roll angle of the spacecraft, the background level, and the photometric centroid of the target (e.g. Lendl et al. 2020; Delrez et al. 2021; Morris et al. 2021; Barros et al. 2022; Deline et al. 2022; Krenn et al. 2023; Bonfanti et al. 2024; Demangeon et al. 2024). We included the photometric variability caused by systematics in our global model using Gaussian processes (GPs) and polynomial correlation, as detailed in Section 4.1.

3.2 TESS observations

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observed the star WASP-18 in sectors 2, 3, 29, 30, and 69, spanning 5 years from August 2018 to September 2023, and covering up to 115 phase curves of planet b. We downloaded the data sets from the Mikulski Archive for Space Telescopes (MAST⁵). The photometry of sectors 2 and 3 was available in a 2 minute cadence, while sectors 29, 30, and 69 had both 2 minute and 20 second cadences. To optimise the computational cost our our analysis, we only used the 2 minute cadence data for all sectors. TESS photometry was provided by the Science Processing Operations Center (SPOC) in two formats that have undergone different levels of processing. The simple aperture photometry (SAP) is a sum of pixel values within predefined non-circular apertures computed on calibrated images. The pre-search data conditioning SAP (PDCSAP) goes further by correcting the SAP photometry using information from the most common systematic features stored in the so-called cotrending basis vectors (CBVs; Kinemuchi et al. 2012). The resulting PDCSAP photometry is usually much cleaner than SAP, with fewer long-term trends. However, the PDCSAP flux can sometimes feature systematic trends that are not present in the SAP or, in the worst cases, have some low-amplitude astrophysical signals removed. After a careful inspection of the TESS light curves from all sectors, we found that the PDCSAP provides a better photometry with significant corrections of red noise for all sectors, except for sector 3 for which PDCSAP flux has strong systematics that are not present in the SAP flux. We thus selected the cleaner SAP for sector 3 and PDCSAP for the other sectors (2, 29, 30, and 69) as the base TESS flux for this work. We used the QUALITY flag provided for the SPOC pipeline to discard 16 565 out of 93 313 data points (17.75%). Given the focus of our analysis is to measure the phase-curve signal of WASP-18 b, we have not included any complex time-dependent systematic correction (e.g. spline functions or Gaussian processes) as we might lose our ability to accurately retrieve the planetary signal. We thus discarded parts of the light curves with remaining trends or high noise to keep only the cleanest photometry of each sector, as shown in Fig. B.1 and listed in Table B.1. From the best-fit residuals obtained with the model described in Section 4, we performed a 4-σ clipping to remove 28 outliers out of the remaining 67 190 photometric points. The final TESS light curves that we obtained after discarding trends and outliers cover about 101 full orbital period of WASP-18 b.

3.3 Spitzer observations

We recovered the observational data of WASP-18 with the Infrared Array Camera (IRAC) of the Spitzer space telescope (Werner et al. 2004) from the Spitzer Heritage Archive (SHA⁶). We downloaded all available data sets from three General Observing (GO) programmes: 50517 (PI: Harrington), 60185 (PI: Maxted), and 11099 (PI: Kreidberg). Programme 50517 observed two planetary occultations in December 2008, both observed simultaneously with the IRAC channels 1 and 3, and 2 and 4, respectively (Nymeyer et al. 2011). Programme 60185 included two full-orbit observations of WASP-18 b in January and August 2010 with the IRAC channels 1 and 2, respectively (Maxted et al. 2013). Programme 11099 covered ten occultations of planet b with the IRAC channel 2 in September 2015 that were published in Deming et al. (2023). A detailed log of the Spitzer observations is listed in Table C.1.

We re-extracted the photometry and applied a pre-processing correction that models the IRAC intra-pixel sensitivity (Ingalls et al. 2016) using the bilinearly interpolated subpixel sensitivity (BLISS) mapping (Stevenson et al. 2012a). We also included the possibility for a linear decorrelation as a function of the full-width at half maximum (FWHM) of the pixel response function (PRF). The modelling uncertainties of this correction were added quadratically to the errorbars of the data to ensure consistent error propagation. The final corrected light curves were sampled at cadences of 10.8 and 13.6 seconds (channels 1 and 2 and channels 3 and 4, respectively) for programme 50517, 27.4 seconds for programme 60185, and 129.4 seconds for programme 11099. Demory et al. (2016) and Bourrier et al. (2022) provide a detailed description of this pre-processing step.

When analysing the full-orbit light curves of programme 60185 (PI: Maxted), we noticed abnormal levels of red noise. We made several reduction attempts varying the noise parametrisation in order to mitigate the systematics without success. We also reduced the data with another independent pipeline (Stevenson et al. 2012a,b; Campo et al. 2011; Cubillos et al. 2013, 2014, 2017; Bell et al. 2019), but we ended up not being able to detrend the systematics from the astrophysical signal. The residual red noise was significant and the astrophysical signal depended heavily on the modelling choices without a clear optimal procedure (e.g. which data points or systematics models to include). This resulted in significant inconsistencies for some parameters when comparing these two phase curves with the other light curves from Spitzer, but also CHEOPS and TESS, in particular, for the noise levels, the occultation depths, and the mid-transit times. These inconsistencies were biasing our results when jointly analysing all the light curves. Without being able to explain the source of this red noise and without being able to retrieve reliable results out of these data sets, we chose to discard these two Spitzer phase curves.

We computed a joint fit to the remaining Spitzer data using the modelling of astrophysical signal and systematic noise described in Section 4 and performed a 4.5-sigma clipping to the best-fit residuals. This resulted in flagging and discarding a total of four outliers out of 8519 data points (0.05%).

4 Light-curve analysis

4.1 Normalisation and noise modelling

We adopted a similar normalisation and noise modelling for each of the observations analysed in this work; namely, each CHEOPS visit, each TESS orbit (2 orbits per sector), and each Spitzer visit. At least two parameters were fitted for each observation: the reference flux value used to normalise the photometric flux (see Section 4.2 for details about what a normalised flux of 1 corresponds to) and the noise jitter, σ_w, that is added (or subtracted if σ_w < 0) quadratically to the photometric error bars to account for the under- or over-estimation of the uncertainties of the data points. For each of the CHEOPS and Spitzer visits, we evaluated whether an additional linear slope was necessary to properly model the data. For this, we fit a full model including slopes for all visits to the data and identified every visit with a significant slope (value inconsistent with 0 by > 3σ). All other slopes were fixed to 0 in our final model. This resulted in fitting a slope for CHEOPS visits 3, 9, 31, 32, 33, and 35, and for Spitzer visits 3, 7, 8, 9, and 10.

These corrections are enough for the TESS and Spitzer data sets, as the photometric light curves have been pre-processed with the PDCSAP and the BLISS mapping, respectively.

CHEOPS photometry, however, is not corrected and is strongly affected by systematics, mainly due to background level variations and flux modulation due to the spacecraft rolling around its line of sight, which results in having the field of view rotating around the target star. We included a correction for the induced systematic noise as a function of the background flux and the roll angle of the satellite. We modelled the flux variations as a function of the background with the following formula: $$F=a_{\text{bkg}}{\left[{\log }_{10}\left(B/B_0\right)\right]}^2+b_{\text{bkg}}{\log }_{10}\left(B/B_0\right),$$(1)

where F is the photometric flux, B is the background level, B ₀ is a fixed reference background level, and a _bkg and b _bkg are fitted parameters. We compared the Bayesian information criterion (BIC) of the model with only a _bkg free (b _bkg = 0), only b _bkg free (a _bkg = 0), or both parameters free. The latter option was shown to be significantly favoured (ΔBIC > 70) and we selected it for our flux-background model. We determined the fixed value of B ₀ = 188 425.83 e⁻ from a Gaussian fit to the histogram of all background values, which provides a mean (or median) estimate that is robust to outlier.

The remaining systematic noise in CHEOPS light curves is strongly modulated by the rotation of the field of view, either due to background stars or any diffuse light source (e.g. Earth straylight). These effects on the photometry repeat themselves with each one of CHEOPS’ orbits and can be accurately modelled as a function of the spacecraft roll angle. We adopted a Gaussian process (GP) modelling of the noise as a function of the roll angle, using a Matérn-3/2 kernel from the celerite2 package (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). This means that our GP was fitted to the residual flux phase-folded on the spacecraft roll-angle values, which are ranging roughly from 15 deg to 322 deg. Given the versatility of GPs, we chose to use the same set of hyper-parameters (the standard deviation of the process σ_GP and the correlation scale in roll-angle unit ρ_GP) for all CHEOPS visits, with each visit being modelled independently from the others. This means that every visit has its roll-angle dependency modelled with the same hyper-parameters, but only using the data points of that specific visit. This choice was found to be a good compromise between modelling flexibility and the number of free parameters.

All the steps normalising the flux and modelling the noise were conducted simultaneously with the astrophysical model fit.

4.2 Astrophysical model

4.2.1 Light travel time

We modelled the astrophysical signal of the WASP-18 system by dividing it into two additive contributions: the planetary flux and the stellar flux. Both contributions are accounting for the light travel time (LTT) through the planetary system and synchronised at the time of inferior conjunction, T ₀ (mid-transit time). The corresponding correction is done by converting the observation times, t _obs, into LTT-corrected times, t _corr, using the following equation: $$t_{\text{corr}}=t_{\text{obs}}-\frac{a}{c}[1-\cos \left(2\pi \frac{t_{\text{obs}}-T_0}{P}\right)]\sin \left(i\right),$$(2)

where a is the semi-major axis of the planetary orbit, c is the speed of light, T ₀ is the time of inferior conjunction, and P and i are the orbital period and inclination, respectively. We note that Eq. (2) assumes a circular orbit. We implemented a LTT correction for eccentric orbits but the computational time is much longer and we estimated the model difference to be smaller than ±0.15 ppm (assuming the eccentricity value of 0.0091 from Nymeyer et al. 2011). Therefore, we opted for the circular approximation of the LTT correction. Every time we estimated the model values (i.e. at each step of our parameter space exploration), we first computed the parameter-dependent LTT correction and then calculated our model at times t _corr. The amplitude of the correction for WASP-18 b is of the order of 20 seconds.

Our model fits directly for the orbital period, P, the orbital inclination, i, and the normalised semi-major axis, a/R_★ , where R_★ is the stellar radius. To compute the LTT correction, we included in our model a free parameter R_★ with a normal prior based on the value derived from the IRFM method (see Section 2 and Table 1). This new parameter allows us to compute the LTT, while accounting for uncertainties on the stellar radius.

4.2.2 Planetary flux

We modelled the flux coming from the planet as a combination of the phase curve of the planet and an occultation model when the planet is hidden by the star. The phase curve is computed by a sine function with a possible phase offset with respect to the orbital phase of the planet. This was implemented using the following equation (details provided in Appendix D): $$F_{\text{p}}=\frac{1-\sin (\omega +\nu -\mathrm{\Delta }\varphi )\sin \left(i\right)}{2}(F_{\text{max}}-F_{\text{min}})+F_{\text{min}},$$(3)

where ω is the argument of periastron⁷, ν is the true anomaly, Δφ is the phase offset (e.g. hotspot offset), i is the orbital inclination, and F _min and F _max are the minimum and maximum values of the phase curve as seen from i = 90 deg, respectively. We note that when Δφ = 0, we have F _min and F _max being the nightside and dayside fluxes, respectively. Also, if i ≠ 90 deg, then the value of F _max is greater than the observed occultation depth (even if Δφ = 0).

We computed the planetary flux model by multiplying the phase-curve, F_p , by an occultation model, F _occ, from the batman package (Kreidberg 2015) that has been normalised in such a way that the out-of-occultation value is 1 and the in-occultation value is 0. The occultation model has up to seven free parameters: the mid-transit time, T ₀, the orbital period, P, the planet-to-star radii ratio, R _p /R_★ , the normalised semi-major axis, a/R_★ , the orbital inclination, i, and the eccentricity, e, and the argument of periastron, ω, implemented in the form: {e cos ω, e sin ω}.

4.2.3 Stellar flux

Our model of the stellar flux only includes variability induced by the planet: the ellipsoidal variations, the Doppler boosting and the planetary transit. Other sources of variability, such as granulation or star spots, have not been detected in the data sets and were not modelled. Their small contribution (if any) is thus accounted for by the noise jitter, σ_w, described in Section 4.1. The stellar flux can be expressed with the following equation: $$F_{\star }=(F_{\text{EV}}+F_{\text{DB}})\times F_{\text{tra}},$$(4)

where F _EV is the ellipsoidal variation, F _DB is the Doppler boosting, and F _tra is the flux dip induced by the planetary transit.

The ellipsoidal variations (EVs) are caused by the deformation of the stellar sphere due to the gravitational pull (tidal forces) of the planet. The bulge created on the stellar surface rotates together with the planet, in and out of view of the observer, which in turn induces photometric variability. This effect is well described in Esteves et al. (2013) where they present a model derived from the work on binary stars of Morris (1985) (Eqs. (1)–(3)), which itself is based on the proposed cosine series expansion of Kopal (1959) (Eq. (IV-2-37)). For consistency, Appendix E details the steps to derive the EV model used in this work resulting in the following equations (Eqs. (5)–(10)). Similarly to Eq. (8) of Esteves et al. (2013), we express the ellipsoidal variations of circular orbits using the first three dominating terms of the expansion: $$F_{\text{EV}}=A_{\text{EV}}[\cos \left(2\theta \right)-A_1\sin \left(\theta \right)+A_3\sin \left(3\theta \right)+A_0],$$(5)

where θ = ω + ν with ω the argument of periastron and ν the true anomaly. We chose to define the zero-flux level F _EV = 0 at mid-occultation time (i.e. θ = 3π/2) by adding a constant offset A ₀ = 1 - A ₁ - A ₃ to the EV model. Similarly to Esteves et al. (2013), we can write the coefficients A _EV, A ₁, and A ₂ for circular orbits as follows: $$A_{\text{EV}}={\alpha }_{\text{EV}}\frac{M_p}{M_{\star }}{\left(\frac{R_{\star }}{a}\right)}^3{\sin }^2\left(i\right),$$(6) $$A_1=3{\beta }_{\text{EV}}\frac{R_{\star }}{a}\frac{5{\sin }^2\left(i\right)-4}{\sin \left(i\right)},$$(7) $$A_3=5{\beta }_{\text{EV}}\frac{R_{\star }}{a}\sin \left(i\right),$$(8)

where M _p and M_★ are the masses of the planet and the star, respectively, R_★ is the stellar radius, a is the semi-major axis, and i is the orbital inclination. The coefficients α_EV and β_EV can be expressed as a function of limb-darkening and gravitydarkening coefficients (LDC and GDC, respectively). Esteves et al. (2013) only uses the linear LDC, and we recalculated the expression to account for the quadratic term based on Kopal (1959) and using the limb-darkening law from Manduca et al. (1977) ⁸ implemented in batman for instance. We obtain the following expressions: $${\alpha }_{\text{EV}}=\frac{3}{10}\frac{15+u_1+2u_2}{6-2u_1-u_2}(1+y),$$(9) $${\beta }_{\text{EV}}=\frac{5}{168}\frac{35u_1+22u_2}{15+u_1+2u_2}\frac{2+y}{1+y},$$(10)

where u ₁ and u ₂ are the LDC, and y is the GDC. Note that when u ₂ = 0, we retrieve the equations (12) and (13) of Esteves et al. (2013). These expressions allow us to include the quadratic LDC values fitted to the transit light curve to compute the ellipsoidal variations. Given that a/R_★ and i can be derived from the transit/occultation models, the EV model only introduces two additional free parameters: the amplitude A _EV and the GDC y (via β_EV). We note that not accounting for the secondary terms in the EV model (i.e. A ₁ = A ₃ = 0) results in a difference in EV amplitude of the order of 50 ppm in the CHEOPS passband, which is detectable with the photometric precision (see Fig. 1).

We computed the gravity-darkening coefficient for each passband of our data sets. For Spitzer and TESS, we used directly the values of y tabulated in Claret & Bloemen (2011) and Claret (2017). Based on the stellar properties T _eff, log g, [Fe/H] listed in Table 1, we selected all y values within ±4σ of the stellar parameters, and computed the mean and standard deviation. For CHEOPS, Claret (2021) provides two values, y ₁ and y ₂, that relates to the GDC as y = β₁ y₁ + y ₂, where β₁ is the gravitydarkening exponent (GDE). From the work of Claret (2004), we extracted the values of the GDE using the stellar parameters and their uncertainties and computed β₁ = 0.305 ± 0.004 for WASP-18. Table 2 lists the values of the gravity-darkening coefficients obtained for the different passbands.

We estimated that variations of the GDC of about 0.01 (order of magnitude of our uncertainties) will induce changes of about 0.2 ppm in our EV model. Therefore, for our analysis, we neglected these error bars and fixed the values of the GDC to their mean values for all passbands.

The Doppler boosting (DB) is a combination of two effects related to the motion of the star induced by the gravitational pull of the planet: the first contribution is the Doppler beaming that concentrates the stellar flux in the direction of motion; and the second contribution is the Doppler effect creating a passband-dependent variability due to the blue- and red-shift of the stellar light. The sources and modelling of DB are discussed in more details in the literature (e.g. Hills & Dale 1974; Rybicki & Lightman 1979; Loeb & Gaudi 2003; Zucker et al. 2007; Bloemen et al. 2011; Barclay et al. 2012; Esteves et al. 2013). Overall, the photometric effect of DB is proportional to the radial velocity of the star and can thus be modelled with the following expression: $$F_{\text{DB}}=-A_{\text{DB}}[\cos (\omega +\nu )+e\cos \left(\omega \right)],$$(11)

where A _DB is the amplitude of the Doppler boosting, ω is the argument of periastron, ν is the true anomaly, and e is the eccentricity.

Finally, the stellar flux variability due to ellipsoidal variations and Doppler boosting is further modulated by the partial flux loss when the planet transits its host. We thus multiply the sum of the two effects by a transit model (see Eq. (4)) also from the batman package (Kreidberg 2015) sharing the same parameters as the occultation model (Section 4.2.2), but also including two quadratic limb-darkening parameters, u ₁> and u ₂.

The final model including all astrophysical signals from the planetary system is therefore: $$F_{\text{tot}}=F_p\times F_{\text{occ}}+(F_{\text{EV}}+F_{\text{DB}})\times F_{\text{tra}}.$$(12)

Fig. 1 Top: ellipsoidal variations modelled with and without secondary terms for the best-fit parameters in the CHEOPS passband. The x-axis represents the orbital phase of the planet θ = ν + ω, where ν and ω are the true anomaly and the argument of periastron, respectively. The blue curve (Sine model) is computed as a simple sine wave, i.e. with A 1 = A 3 = 0. The other models include the secondary terms that are computed for several values of the gravity-darkening coefficient y. Bottom: differences between the different EV models and the Sine model. The peak-to-peak difference in amplitude is of the order of 50 ppm for the realistic range of GDC values (0.1 ≤ y ≤ 0.5). The effect of accounting for A 1 and A 3 in the modelling of EV impacts mainly the difference in baseline level between transit and occultation. The baseline flux during transit will be lower than that during occultation (effect similar to a negative nightside flux). If not taken into account, this may lead to overestimating the occultation depth or underestimating the nightside flux from a planet.

4.3 Parameter space exploration and joint fit

The total number of parameters of our model for all passbands can be up to 187, among which 137 are nuisance parameters for noise modelling (84 for CHEOPS, 20 for TESS, and 33 for Spitzer). The 50 remaining parameters define the planetary model described in Section 4.2, with eight related to the planetary orbit (T ₀, P, R _p /R_★ , a/R_★ , i, e cos ω, e sin ω, R_★ ) and seven for each of the six passbands (F _max, F _min, Δφ, A _EV, A _DB, u ₁, u ₂). We did not fit for the passband-dependent variations of R _p /R_★ as the only strong constraints on this parameter came from the transits of CHEOPS and TESS, which have similar passbands hence similar expected planetary radii (negligible difference due to Rayleigh scattering). Moreover, Spitzer occultations would only poorly constrain R_p /R _★ from the durations of the occultations.

Given Spitzer data only covers planetary occultations, we fixed the phase-curve and transit parameters (F _min, A _EV, A _DB, u ₁, u ₂) to 0 as they could not be constrained by the light curves. In addition, we fitted the data with all parameters free and we obtained values for the minimum PC level F _min (CHEOPS and TESS) and the phase offset Δφ consistent with zero (less than 3 σ significance). These parameters were thus fixed to 0 for the final analysis. To further reduce the size of the parameter space and considerably improve the computational time, we assumed the orbit of WASP-18b to be circular (e cos ω = e sin ω = 0). This choice is motivated by the small value of the eccentricity (e ~0.009) reported in previous works (Hellier et al. 2009; Triaud et al. 2010; Nymeyer et al. 2011; Csizmadia et al. 2019). Furthermore, a small eccentricity with the peculiar orbital orientation of WASP-18 b (ω ~ 270 deg) could be explained by radial-velocity signals induced by the bulge of the tidally deformed star (Arras et al. 2012)⁹. In total, we fixed 30 parameters, reducing the total number of free parameters to 157.

We explored the large parameter space with the Nested sampling method (Skilling 2004, 2006) implemented in the package dynesty (Speagle 2020; Koposov et al. 2023). We used a Dynamic Nested Sampler (Higson et al. 2019) with bounding option using multiple ellipsoid bounds (Feroz et al. 2009) and a random slice sampling (Neal 2003; Handley et al. 2015a,b). We started the run with 2000 live points, which corresponds to more than 10 times the number of dimensions and is good compromise between computational time and good convergence. The weight and stopping functions of the dynamic sampling were set to the default values, i.e. a relative fractional importance of 80% on the posterior for the weight, and a stopping criterion of d log Z ≤ 0.01 (with Z being the estimated Bayesian evidence). The run ended with a total number of 333 506 points sampling the parameter space.

5 Refined system properties

5.1 Planetary parameters

The final values of the planetary parameters obtained from our analysis are listed in Table 3. The values of the nuisance parameters fitted for each instrument are listed in appendix (Tables I.1, I.2 and I.3 for CHEOPS, TESS and Spitzer, respectively).

The final detrended light curves are presented in Figs. 2, 3, and 4, for the 36 CHEOPS visits, the 5 TESS sectors, and the 14 Spitzer visits, respectively.

We retrieve a time of inferior conjunction T ₀ consistent at 1.15 σ with the value of Shporer et al. (2019) with a precision improved from 2.2 s to 1.6 s. The uncertainty on T ₀ reaches a minimum of less than 1 s for the date T _0,opt = 2 459 251.66201 BJD_TDB (i.e. on February 6, 2021). The precision of 1.4 ms on the orbital period P is also a significant improvement with respect to previously published results (138 ms for Shporer et al. 2019, 21 ms for Cortés-Zuleta et al. 2020, or 6 ms for Coulombe et al. 2023). The values of P are all consistent within 1 σ. The normalised semi-major axis a/R_★ and the orbital inclination i are both consistent with the values from the literature, with the largest differences from Shporer et al. (2019) of 2.68 σ and 2.15 σ, respectively. This might be explained by the fact that the analysis of Shporer et al. (2019) fitted the data with fixed limb-darkening coefficients and a fixed non-zero eccentricity.

The planet-to-star radii ratio, R _p /R_★ , that we derived from our global fit is consistent with the values from Shporer et al. (2019) (2.24 σ) and Coulombe et al. (2023) (<1 σ). It is, however, significantly smaller than the value from Cortés-Zuleta et al. (2020). We refined the precision down to 0.13% on R _p /R_★ , which converts into 0.65%, or 550 km, on the absolute planetary radius.

In Section 4.3, we explained that the phase offset, Δφ, and the minimum planetary flux, F _min, were fixed to zero as their values were always consistent with zero when let free. The difficulty to constrain Δφ mainly came from the strong degeneracy between a hotspot offset and the effect of Doppler boosting (strong DB can be counterbalanced by having large westward hotspot offset), leading to unrealistically large values for both parameters. We note that fixing Δφ to 0 makes the values of F _max and F _min equivalent to the dayside flux (or occultation depth) and nightside flux, respectively. We were able to fit for F _min in both the CHEOPS and TESS passbands, as we analysed the full phase curves and we derived the following 3 σ upper limits: $F_{\text{min}}^{\text{CHEOPS}}<12.3\text{ppm}$ and $F_{\text{min}}^{\text{TESS}}<18.4\text{ppm}$.

When comparing the values of F _max, A _EV and A _DB for the TESS passband to previously published values, we obtain consistency within ~1σ for all parameters of Shporer et al. (2019) and Coulombe et al. (2023), except for the EV of Shporer et al. (2019) that is marginally consistent at 2.89 σ. The precision on all three parameters is significantly improved. We note that the non-detection of the nightside flux is in agreement with the analysis of Coulombe et al. (2023), and we lowered the 3 σ upper limit by a factor of 2.4 (from 44 ppm to 18.4 ppm).

Following Equations (6) and (9), we could convert the detected EV amplitudes A _EV in both CHEOPS and TESS passbands to planetary masses. We obtained $M_p^{\text{CHEOPS}}={10.03}_{-0.52}^{+0.53}{M}_J$ and $M_p^{\text{TESS}}=8.73\pm 0.45M_J$, which are consistent with each other (1.89 σ) and with the literature. Despite their mutual consistency, the mass difference between the two passbands is probably due to inaccuracies of the GDC estimates (Table 2) between reality and theoretical modelling.

Finally, the dayside flux values, F _max, that we derived for the four Spitzer/IRAC channels are all consistent (<1.1 σ) with the literature (Nymeyer et al. 2011; Maxted et al. 2013; Sheppard et al. 2017; Garhart et al. 2020).

5.2 Orbital eccentricity and argument of periastron

We considered the planetary orbit to be circular for our fit given the small eccentricity reported in the literature and the fact that it may actually be explained by a spurious signal from the bulge of the tidally deformed star (see discussion in Section 4.3). However, Coulombe et al. (2023) recently observed the occultation of WASP-18 b with JWST and derived a mid-occultation time of $T_{\text{occ}}^{\text{JWST}}=2459802.881867\pm 0.000092{\text{BJD}}_{\text{TDB}}$ that is inconsistent with our results. Indeed, the parameter set we compute gives an occultation time of T _occ = 2459802.8826151 ± 0.000015 BJD_TDB after including the LTT of 20.26 ±0.15sec. This means that there is a significant timing inconsistency of ΔT _occ = −64.6 ± 8.1 sec (7.98 σ).

We first discarded the possibility of orbital period variation as the JWST observation occurred in-between the datasets considered in our analysis (nine CHEOPS visits within 2.5 months and TESS sector 69 fewer than 400 days after). We then considered the possibility of having an effect due to the thermal structure of the atmosphere and the fact that different passbands may probe different dayside brightness distribution. As shown in de Wit et al. (2012), this might result in mid-occultation time shifts with respect to expectations from the transit. However, even if we could not reach the required precision in the TESS data, we found a hint of a similar shift ($\mathrm{\Delta }{T}_{\text{occ}}^{\text{TESS}}=-48\pm 32\text{sec}$). In addition, the detection of orbital eccentricity reported from radial velocity (Hellier et al. 2009; Triaud et al. 2010; Nymeyer et al. 2011; Csizmadia et al. 2019) made us consider the eccentric nature of the planetary orbit as the most probable source of the timing inconsistency.

We converted Δ T _occ to orbital phase and obtained an offset of 0.00079 ± 0.00010 corresponding to an occultation phase of 0.49921 ± 0.00010. Due to the high computational cost of running a global fit with an eccentric orbit, we derived constraints on the eccentricity e and the argument of periastron ω indirectly. We considered 3 different sets of prior values on e cos ω and e sin ω from the works of Triaud et al. (2010), Nymeyer et al. (2011) and Csizmadia et al. (2019) that we combined with the likelihood probability of these 2 parameters to match the mid-occultation time $T_{\text{occ}}^{\text{JWST}}$ reported in Coulombe et al. (2023). We explored the 2D-parameter space {e cos ω, e sin ω} with the Markov chain Monte Carlo (MCMC) algorithm emcee (Foreman-Mackey et al. 2013, 2019) implemented in Python. The results from our 3 different approaches (3 different priors) are fully consistent (<1.2 σ) with each other (see posterior distributions in Fig. 5). Table 4 shows the published values of e and ω (used as priors in the form e cos ω and e sin ω) and the results we obtained with the mid-occultation time measured with JWST. Eccentricity values are largely consistent (<0.2 σ) with the published values with a significant improvement in precision for the value of Csizmadia et al. (2019). The same is true for the argument of periastron but with a weaker consistency (≤1.8 σ). We also note that all posterior distributions of ω values converge towards a similar angle of −262 deg despite the 3.2 σ inconsistency between the published values of Nymeyer et al. (2011) and Csizmadia et al. (2019). This highlights the strong constraint provided by the timing inconsistency with $T_{\text{occ}}^{\text{JWST}}$, which translates into unprecedented precision on the value of the argument of periastron (less than ±1.5 deg). As the posterior distribution derived from the published values of Triaud et al. (2010) is the most conservative among our 3 approaches (see blue data in Fig. 5), we selected it for our reference orbital parameters.

Following the confirmation of the eccentric nature of the orbit of WASP-18 b, we computed the circularisation timescale τ_circ of the system using several formulae from the literature (Eq. (1) of Barros et al. 2011; Eqs. (2) or (3) of Adams & Laughlin 2006). We found a 3 σ upper limit of τ_circ ≲ 20 Myr, which is much shorter than the age of the system t _★ = 1.4 ± 0.7 Gyr (see Table 1). The planet should thus have had its orbit circularised and explaining what this is not the case is challenging. A possible explanation could be that an outer companion is slowing down the damping of the eccentricity of WASP-18 b by injecting energy into the planetary orbit inducing (eccentricity excitation; Mardling 2007; Laskar et al. 2012). Pearson (2019) reported a candidate outer planet with an orbital period of 2.16 days and a mass of ~0.2M _J , but there is no firm detection of this object. Another perturbing body that could cause eccentricity excitation is the stellar companion WASP-18 B that is a late-M dwarf at a distance of about 3500 au (Csizmadia et al. 2019). Further observations and analyses of the radial velocity and transit-time variations of the system might help understand the dynamical mechanisms at play and constrain the properties of the potential perturbing body.

Fig. 2 CHEOPS phase-folded phase curve. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. From top to bottom are show the full phase curve, a zoomed-in version to highlight the phase-curve signal, and the residuals in ppm.

Fig. 3 TESS phase-folded phase curve. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. From top to bottom are show the full phase curve, a zoomed-in version to highlight the phase-curve signal, and the residuals in ppm.

Fig. 4 Spitzer phase-folded occultations. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. The 4 IRAC channels 1, 2, 3, 4, centred at 3.6 μm, 4.5 μm, 5.8 μm, 8.0 μm, respectively, are shown at the top, with their corresponding residuals in ppm at the bottom.

Fig. 5 Correlation plot of the posterior distribution of the eccentricity e and argument of periastron ω of the orbit of WASP-18 b that match the mid-occultation time $T_{\text{occ}}^{\text{JWST}}$ from Coulombe et al. (2023). We used our posterior distribution on planetary parameters from Table 3, and Gaussian priors on e cos ω and e sin ω from Triaud et al. (2010) in blue, Nymeyer et al. (2011) in orange, and Csizmadia et al. (2019) in green. Overall, the 3 posterior distributions are all consistent with each other within <1.2 σ. This plot made use of the package corner.py (Foreman-Mackey 2016).

5.3 Occultation depth variability

We checked for variability of the occultation depth in each instrument passband by fixing all parameters to their best-fit values from our global fit (see Table 3), letting free only the occultation depth (via the dayside flux F _max), the noise jitter, σ_w, and the normalisation flux, f₀ . The 3D parameter space was explored with the MCMC algorithm emcee (Foreman-Mackey et al. 2013, 2019).

We fitted each CHEOPS visit individually and extracted the posterior values of $F_{\text{max}}^{\text{CHEOPS}}$ for each epoch. The dayside flux values are represented in Fig. 6. Naturally, the visits not covering any occultation event are poorly constrained, but we obtain a precision of the order of 30 ppm for each of the other ones. We do not detect any strong sign of variability with discrepancies being below 3 σ. To quantify the dayside flux scattering, we calculated the multiplicative factor to be applied to the errorbar to find a distribution consistent with a Normal distribution. For the CHEOPS’ dayside flux, we found a multiplicative factor of 1.05, which is close to 1 and thus showing consistency with Gaussian statistics. We computed the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) of the series but could not identify any significant periodicity.

We split the TESS sectors into 110 short light curves with a duration of one planetary orbit. Given the TESS cadence, we had about 680 exposures for a full coverage of a planetary orbit, and we discarded 12 out of the 110 data sets that had fewer than 400 points (poor phase-curve coverage). For each of the 98 remaining light curves, we applied the same procedure as for the CHEOPS visits and obtained the $F_{\text{max}}^{\text{TESS}}$ values represented in Fig. 7. The scatter is larger than what is allowed from the uncertainties with 17 data points at >3 σ (all ≲5 σ but 1 at 9.5 σ). The errorbar multiplicative factor is this time of 2.21, further confirming the large inconsistency with Gaussian statistics. We investigated the outliers’ light curves individually and found local variability likely of instrumental origin. No specific periodicity nor dominant timescale could be identified as a significant signal in the Lomb-Scargle periodogram. We note that this variability could also be of stellar origin and induced by supergranulation (e.g. HAT-P-7 b; Lally & Vanderburg 2022).

We selected the Spitzer occultations of the IRAC Channel 2 and fitted them individually similar to what was done for CHEOPS and TESS. The dayside flux values we recover were all consistent with each other (see Fig. 8). The errorbar multiplicative factor is 0.88 for the Spitzer F _max values, hence showing clear consistency with Normal distribution statistics. We compared our individual values to the ones published in Deming et al. (2023) and found most of them being consistent (<3 σ) with only 2 out of 11 showing a significant discrepancy (4.1 σ and 3.2σ for AOR keys 53516800 and 53515520, respectively).

The analysis of the occultation depth variability of all three passbands did not reveal any significant signal on the short timescale. Thanks to the very long periods covered by our data sets (3 years for CHEOPS, 5 years for TESS, and 6.7 years for Spitzer), our results demonstrate long-term stability of the dayside of WASP-18 b.

Fig. 6 Dayside flux observed in the CHEOPS passband over the three seasons of observations (all three panels). The observations covering an occultation are represented in blue (including a phase curve), and the ones only covering transits in pink. The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{CHEOPS}}={211.9}_{-7.5}^{+7.6}\text{ppm}$ from Table 3.

Fig. 7 Dayside flux observed in the TESS passband over the 5 sectors: sectors 2 and 3 on the left, sectors 29 and 30 in the center, and sector 69 on the right. The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{TESS}}=340.4\pm 7.0\text{ppm}$ from Table 3.

Fig. 8 Dayside flux observed in the Spitzer/IRAC/Channel 2 passband for programmes 50517 (left panel) and 11099 (right panel). The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{SPITZER 2}}=3935\pm 25\text{ppm}$ from Table 3.

6 Atmospheric characterisation

6.1 Planet thermal emission

We discussed in Section 5.1 the upper limits that we obtained on the minimum PC flux values, F _min, from both CHEOPS and TESS passbands. Given that we fixed the phase offset Δφ to 0, F _min actually corresponds to the nightside flux. Assuming a black-body emission for the nightside of WASP-18 b and using synthetic stellar spectra from the PHOENIX library (Husser et al. 2013), we could convert upper limits from flux to temperature. We selected the PHOENIX spectrum with the closest properties to our stellar parameters: T _eff = 6300 K, logg = 4.5 and [Fe/H] = 0.0. We modelled the thermal emission of the planet following relationship: $$\frac{F_{\text{night}}}{F_{\star }}={\left(\frac{R_p}{R_{\star }}\right)}^2\frac{{\int }_{\lambda =0}^{+\mathrm{\infty }}\mathcal{S}(\lambda ,T_{\text{night}}){\mathcal{T}}_{\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}}\left(\lambda \right)\lambda d\lambda }{{\int }_{\lambda =0}^{+\mathrm{\infty }}\mathcal{S}(\lambda ,T_{\star }){\mathcal{T}}_{\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}}\left(\lambda \right)\lambda d\lambda }$$(13)

where R _p and R_★ are the planetary and stellar radii, respectively, S(λ, T _night) and S(λ, T_★ ) are the flux emission spectra of the planet and the star, respectively, and ${\mathcal{T}}_{\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}}\left(\lambda \right)$ is the instrumental passband. We downloaded the CHEOPS passband from the CHEOPS Archive Browser¹⁰, and the TESS passband from the TESS Science Support Center¹¹. From $F_{\text{min}}^{\text{CHEOPS}}<12.3\text{ppm}$ and $F_{\text{min}}^{\text{TESS}}<18.4\text{ppm}$, we could derive the nightside black-body temperatures, $T_{\text{night}}^{\text{CHEOPS}}<2090K$ and $T_{\text{night}}^{\text{TESS}}<1970K$, for the two instruments, leading to an overall limit of T _night < 1970 K.

We performed a similar analysis for the dayside emission of WASP-18 b using Eq. (13) and replacing F _night by F _day = F _max, and T _night by the so-called brightness temperature T _b. We used F _max from all passbands but Spitzer/IRAC Channels 3 and 4, as their wavelength coverage goes beyond the limit of 5.5μm of the PHOENIX spectra. We obtained the following brightness temperatures: $T_{\text{b}}^{\text{CHEOPS}}=3027\pm 17K$, $T_{\text{b}}^{\text{TESS}}=2956\pm 10K$, $T_{\text{b}}^{\text{SPITZER 1}}=2957\pm 58K$, and $T_{\text{b}}^{\text{SPITZER 2}}=3132\pm 13K$. The uncertainties are expected to be underestimated given the fact that we did not propagate the uncertainties of the stellar properties by using a single PHOENIX spectrum. Also, we might be slightly biasing our results by assuming a black-body emission for the planet. However, these results already show some interesting outcomes as the TESS and Spitzer/IRAC/Channel 1 passband agree on a brightness temperature of 2960 K, whereas CHEOPS and Spitzer/IRAC/Channel 2 would be capturing an excess of flux (converted into a higher temperature here). The excess in the CHEOPS passband might actually indicate the presence of reflected light as we probe visible wavelengths. As for Spitzer/IRAC/Channel 2, the larger occultation depth could be explained by the presence of CO (Brogi et al. 2023; Yan et al. 2023) that have emission lines at 4.5 μm. Using the framework of Cowan & Agol (2011) and following the methodology detailed in Section 6.2 of Deline et al. (2022), we estimated the values of the Bond albedo A _B and the heat redistribution efficiency ε (more details in Appendix F). As usual and expected for UHJs, we obtained very low values of both A_B and ε implying that the planetary atmosphere absorbs most of the incoming energy and redistributes it inefficiently. The narrowest constraints come from the TESS passbands with A _B ≤ 0.13 and ε ≤ 0.21. Interestingly, the maximum dayside temperature that this analytical approach allows is about 3060 K for the extreme case where A_B = ε = 0. This strengthens the hint for an excess of flux, especially in the passband of Spitzer/IRAC/Channel 2. Given the several assumptions made for this analysis of the dayside emission, we consider the aforementioned results as indicative. In the following Sections 6.2 and 6.3, we analysed the dayside emission of WASP-18 b including the recently published JWST values (Coulombe et al. 2023) to complete our 0.6-to-8.0 μm wavelength coverage of occultation depth measurements (see Fig. 9). We used two different approaches to characterise the planetary atmosphere: General Circulation Modelling (GCM) and atmospheric retrieval.

6.2 General circulation model

We performed forward modelling of the atmosphere of WASP-18 b with the general circulation model (GCM) ExoRad, using the version with full radiative transfer (expeRT/MITgcm; Carone et al. 2020; Schneider et al. 2022)¹². We generated synthetic emission spectra of the planetary dayside using the parameters listed in Table 3 as inputs, and normalised the outcomes with the PHOENIX stellar spectra (Husser et al. 2013) up to 5.5 μm, assuming T _eff = 6300 K, and a matching black-body SED at longer wavelengths. To calculate the radiative forcing in the model, we assumed solar metallicity and equilibrium chemistry. We included correlated-k tabulated opacities combined to 11 spectral bins, corresponding to the S1 resolution as specified in Schneider et al. (2022) for the following species: H₂O (from ExoMol¹³; Tennyson et al. 2016, 2020), Na (Allard et al. 2019), K (Allard et al. 2019), CO₂, CH₄, NH₃, CO, H₂S, HCN, SiO, PH₃ and FeH, as well as H⁻ absorption suitable for an ionised atmosphere.

We generated several GCM simulations to derive dayside spectra. We first tested runs without TiO and VO opacities with strong magnetic drag τ_fric and no drag. Here, we found that already the inclusion of magnetic drag alone yielded a vast improvement in the data fit. We then included TiO and VO opacities (McKemmish et al. 2016, 2019). TiO and VO produce an upper atmosphere thermal inversion that impacts in particular the dayside emission in the optical as observed with CHEOPS and TESS. We then tested various values for magnetic drag and including TiO/VO (see Table 5). Magnetic drag is treated during the simulation by applying uniform friction to the horizontal wind field (u, v) via: $$\frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{t}}=-\frac{u}{{\tau }_{\text{fric}}},$$(14) $$\frac{\mathrm{d}\mathrm{v}}{\mathrm{d}\mathrm{t}}=-\frac{v}{{\tau }_{\text{fric}}},$$(15)

where u and v are the zonal and meridional wind fields, respectively.

The inclusion of magnetic drag is motivated by the substantial degree of ionisation of the dayside atmosphere resulting in the presence of several ions such as Na⁺, K⁺, Ca⁺, Fe⁺, Al⁺ and Ti⁺ (Helling et al. 2019, 2021). The partially ionised flow induces magnetic coupling that in turn affects the winds as already pointed out by Perna et al. (2010); Rodríguez-Barrera et al. (2015). The necessity to account for the impact of magnetic fields in UHJs such as WASP-18 b has been further confirmed several times (Wardenier et al. 2023; Beltz et al. 2022; Demangeon et al. 2024). Recent JWST observations also support the impact of magnetic drag (e.g. WASP-18 b in Coulombe et al. 2023). Table 5 clearly shows that a combination of drag and TiO/VO is needed and that the fit to the data improves with decreasing τ_fric = 10⁶, 10⁵, 10⁴ s. Coulombe et al. (2023) found with the GCM SPARC/MITgcm a lower value of τ_fric = 10³ s. These authors tested, however, only two scenarios with weak τ_fric = 10⁶ s and strong drag τ_fric = 10³ s, respectively. They further did not perform a sensitivity study with the GCM with respect to drag and TiO, as we did in this paper. We thus note that the choice of τ_fric is currently still open for debate, in particular in the uniform drag assumption as implemented here (Perna et al. 2010; Tan & Komacek 2019; Coulombe et al. 2023; Beltz et al. 2022). We found with expeRT/MITcgm that τ_fric = 10⁴ s already effectively disrupts superrotation at the dayside and yields agreement with the JWST emission spectrum. We further note that τ_fric = 10⁴ s is currently the smallest drag time scale that was tested in expeRT/MITcgm. Smaller τ_fric time scales as well as other implementations of magnetic field coupling are currently under investigation. In any case, only full 3D GCMs can assess the impact of magnetic drag on the climate state that shapes horizontal heat transfer and thus the dayside emission. The resulting strong horizontal temperature gradient over the dayside, Δ T > 1000 K, also ensures that the TiO and VO abundances are not constant across the dayside.

Comparing the outcome with the measured occultation depths from CHEOPS, TESS, JWST and Spitzer reveals agreement when both TiO and VO and magnetic drag are used (see Fig. 9 and Table 5). Magnetic drag reduces circulation efficiency, resulting in a hotter dayside such that, even without TiO and VO, the GCM provides a relatively close match to the data, with a reduced χ² (${\chi }_r^2$) of 8.5¹⁴. The agreement with the data is nonetheless further improved when including TiO and VO, even with reduced drag. Our best GCM match to the data is obtained when both TiO/VO and a strong drag are present (${\chi }_r^2=2.6$). Therefore, both the magnetic drag and the presence of a strong, local gas-phase absorber in the optical like TiO and VO known from cool star atmosphere modelling (e.g. Gustafsson et al. 2008; Van Eck et al. 2017) are needed to match the data from the optical to the IR wavelengths. Models without TiO and VO are not shown for clarity but underpredict the emitted planetary flux. We particularly note that the GCM output qualitatively matches with the atmospheric retrieval models (see Section 6.3), including the flattening of the spectrum between 1 and 2 micron due to H₂O dissociation. In cooler planets, this wavelength region is dominated by water bands. We note, however, that while our GCM with TiO/VO and magnetic drag shows an upper atmosphere temperature inversion, it is not sufficiently hot for CO emission. Retrieval models show that temperatures higher than 3500 K are required to trigger CO emission (see Fig. 10).

Fig. 9 Occultation depths as a function of wavelength. The black points with errorbars represent the measurement from this work (CHEOPS, TESS and Spitzer; 2 leftmost and 4 rightmost points) and from Coulombe et al. (2023) with JWST (0.8-3 μm). The GCM simulations including TiO and VO (Section 6.2) with and without magnetic drag are shown in orange and light brown, respectively. The orange-filled diamonds mark the passband-integrated GCM values with magnetic drag (τfric = 104 s). The retrieval runs (Section 6.3) are shown in blue, pink, purple and green (same colours as in Fig. 10 and Table 6) depending on the data points included in the fit. An inset zoomed-in view of the CHEOPS-to-1.75 μm range is shown in the lower right corner for convenience.

6.3 Atmospheric retrieval

We performed atmospheric-retrieval analyses to characterise the thermal emission properties of WASP-18 b, hence, inferring its thermal and composition structure. We employed the open-source PYRAT BAY package (Cubillos & Blecic 2021) to perform the atmospheric modelling, spectral synthesis, and Bayesian posterior sampling. Since the model only accounts for thermal emission, we primarily constrained the retrievals using the infrared JWST and Spitzer occultations, while running additional fits with and without the optical CHEOPS and TESS occultations to assess how stellar reflected light affected the retrieval results.

The atmospheric model consists of a set of parameterized temperature and composition profiles as a function of pressure (81 layers between 100 and 10⁻⁹ bar). We adopted the Madhusudhan & Seager (2009) model to parameterise the temperature profile. The atmosphere of the UHJ WASP-18 b is expected to exceed 2000 K, with the hottest part of the dayside that could reach 3000 K. Since disequilibrium-chemistry processes such as photochemistry and transport-induced quenching become less and less important with increasing effective temperature, at the extreme temperatures seen on WASP-18 b, the chemical reaction rates are fast enough to overcome the effect of disequilibrium chemistry (Kopparapu et al. 2012; Moses 2014; Madhusudhan et al. 2016; Venot et al. 2018). If disequilibrium chemistry occurs at all, it would occur at high altitudes above the pressures probed by the observations presented in this work, and thus the modeled emission spectra of planets such as WASP-18 b would not be significantly impacted (Shulyak et al. 2020)¹⁵. Thermochemical equilibrium is therefore a reasonable assumption to model the dayside atmospheric composition of WASP-18 b and its corresponding emission spectra. The chemical network consists of 45 neutral and charged species, accounting for the main H, He, C, O, N, Na, K, S, Si, Fe, Ti, and V-bearing species expected to form in the atmosphere. We defined three retrieval parameters to sample the composition parameter space: the oxygen elemental abundance relative to solar values, [O/H]; the carbon abundance relative to solar values, [C/H]; anda catchall parameter to scale all other metal abundances relative to solar values, [M/H]. Finally, we computed the height of the pressure layers assuming hydrostatic equilibrium. Table 6 lists the retrieval parameters, priors, and posterior values.

For the spectral synthesis, we considered line-sampled opacities from the dominant molecular species from the HITEMP (CH₄, CO, and CO₂; Rothman et al. 2010) and ExoMol databases (H₂O, HCN, NH₃, TiO, VO, and C₂H₂; Tennyson et al. 2016, 2020). Prior to the retrievals, we tabulated the opacity line lists into a cross-section grid of temperatures, pressures, and wavenumbers (at a constant resolving power of R = 20 000). We pre-processed the larger line lists with the repack package to extract the dominant line transitions (Cubillos 2017). Additionally, we considered opacities from the Na and K resonant lines (Burrows et al. 2000), collision-induced absorption for H₂-H₂ (Borysow et al. 2001; Borysow 2002; Jørgensen et al. 2000) and H₂-He (Borysow et al. 1988; Borysow & Frommhold 1989; Borysow et al. 1989), Rayleigh scattering for H, H₂, and He (Kurucz 1970), and H⁻ free-free and bound-free (John 1988). We adopted a PHOENIX stellar model as the flux spectrum for WASP-18. We performed the Bayesian posterior sampling using the Nested Sampling algorithm (Skilling 2004, 2006) as implemented in the MultiNest package (Feroz et al. 2009; Buchner et al. 2014), employing 1500 live points.

To assess the impact of each observation in constraining the atmospheric properties, we performed four retrieval analyses considering different sets of observing constraints. First, we performed retrievals using only the JWST occultations. Next, we added the Spitzer occultations to test the role of infrared absorbers. Finally, we included the TESS and CHEOPS occultations to test the impact of optical absorbers and reflected light. As expected, the JWST observations dominate the retrieval results, resulting in a consistent fit to the H₂O bands between 1.0 and 2.5 μm. However, we found that retrievals with and without the Spitzer constraints lead to widely different temperature and composition outcomes (see Figs. 10 and 11, and Table 6). The driving factor is the high signal-to-noise IRAC2 occultation depth at 4.5 μm. The much larger brightness temperature at 4.5 μm than that of the JWST occultations requires an atmosphere with significant CO/CO₂ absorption (as highlighted by the diagonally hatched area in left panel of Fig. 10), which in turn leads to super-solar carbon and oxygen abundances of 3-5 times solar. We also noted that incorporating the Spitzer occultations to the retrieval constraints results in steeper temperature gradients. This adjustment is expected given the stark brightness-temperature difference between the JWST and the Spitzer 4.5-μm observations (Fig. 10, right panel). In contrast, the retrieval of the JWST data alone returns an atmosphere with solar-to-subsolar abundances, well in agreement with the previous analysis by Coulombe et al. (2023).

Finally, we compared the retrieved thermal-emission spectra with the measurements in the CHEOPS and TESS passbands to look for the presence of reflected light coming from the planet. The left panel in Fig. 10 displays the retrieval models integrated over all instrumental passbands as coloured diamonds. Table 7 lists the thermal flux values retrieved for each one of the four retrieval models. We detail below the comparison within the CHEOPS and TESS passbands in the cases for which the Spitzer data were included in the retrieval fit (pink, purple and green). First, when only JWST and Spitzer data points were included in the retrieval fit (pink model), we found that the measured occultation depth in the TESS passband was significantly below the retrieval value (at 4.1 σ). Therefore, the inclusion of the TESS data point in the retrieval runs is necessary to properly model the data by mitigating the amount of thermal emission in the TESS passband. With the TESS data point included, we compared two retrieval runs with and without fitting the CHEOPS data point (green and purple models, respectively). When discarding CHEOPS from the fit, the retrieval model showed a smaller flux difference in the TESS passband (−11.1 ± 8.2 ppm). If we convert the flux difference into a 3-σ upper limit on the reflected light in the TESS passband, we get $F_{\text{refl}}^{\text{TESS}}<13.5$ ppm, which corresponds to a geometric albedo¹⁶ of $A_g^{\text{TESS}}<0.017$.

The passband-integrated value of this model for CHEOPS gives a thermal flux of ${190.7}_{-3.5}^{+3.3}$ ppm, which is smaller than the measured dayside flux by 2.6 σ. Explaining this difference with reflected light means $F_{\text{refl}}^{\text{CHEOPS}}=21.2\pm 8.2$ ppm and a corresponding geometric albedo of $A_g^{\text{CHEOPS}}=0.027\pm 0.011$. In our last run, we included all data points in the fit (green model) and analysed the effect on the limits of reflected light in both TESS and CHEOPS passbands. As listed in Table 7, we obtained a larger flux difference for TESS and smaller for CHEOPS, at −2.2 σ and +2.0 σ, respectively. These correspond to limits on geometric albedos of $A_g^{\text{TESS}}<0.009$ and $A_g^{\text{CHEOPS}}=0.020\pm 0.010$.

Our retrievals indicate that a significant fraction of the short-wavelength thermal emission originates from TiO/VO (Figure 10, right panel). Although the high planetary mass of WASP-18 b can potentially enhance the efficiency of cold-trapping - hence depleting heavier metals from the atmosphere (Beatty et al. 2017) -, our retrieved thermal profiles lie well above the thermal stability curves for Ti-bearing and other refractory species. This suggests that the dayside of WASP-18 b remains sufficiently hot to inhibit vertical cold-trapping (Spiegel et al. 2009), and thus allow for the presence of heavy metals like TiO, VO, and alkali metals in gas form in its atmosphere. While day-night cold-trapping remains a plausible scenario, additional phase-curve observations of the nightside emission are needed to evaluate its potential (Parmentier et al. 2013). However, Beatty et al. (2017) argued that observations of planets with daysides of 3000 K and hotter like this planet appear to escape nightside trapping. Figure 10 also elucidates why this would be the case. The dayside is too hot for TiO thermal stability even down to 100 bar. For such a high gravity and fast rotating planet, superrotation is expected to extend well down to 100 bar and deeper into the planet (Carone et al. 2020). Thus, an efficient horizontal heat and material transport is expected below 10 bar. Even if TiO would condense out (i.e. be thermally stable) on the nightside, it would evaporate again at depth and be transported towards the dayside, where it could again dominate the gas absorption in the optical wavelength range in agreement with the results in this work. This is also in line with the study of Helling et al. (2019), where they infer the absence of clouds on the dayside of WASP-18 b from cloud formation models.

We conclude that both the scenario including TESS only and the scenario including CHEOPS+TESS strongly constrain the amount of reflected light in the TESS passbands to extremely low values (< 0.017). This is not unexpected as many atmospheres of UHJs have been reported to be very dark (e.g. Wong et al. 2021). In the shorter wavelength CHEOPS band, we find tentative (2.6 σ) evidence for reflected light, indicative of stronger atmospheric scattering than in the TESS passband. Even so, it’s amplitude remains low, constrained to a geometric albedo below 0.059. The overall reflectivity remains nonetheless very low as well for CHEOPS, in line with expectations for UHJs. As the geometric albedo is wavelength-dependent, the difference in passbands between CHEOPS and TESS prevents us from performing a quantitative comparison between the two instruments. Fig. 12 shows how the reflected light detection of WASP-18 b computed in this work aligns well with previously published $A_g^{\text{CHEOPS}}$ values for other hot exoplanets (Lendl et al. 2020; Hooton et al. 2022; Deline et al. 2022; Brandeker et al. 2022; Scandariato et al. 2022; Parviainen et al. 2022; Demory et al. 2023; Krenn et al. 2023; Pagano et al. 2024; Singh et al. 2024; Demangeon et al. 2024; Akinsanmi et al. 2024; Scandariato et al. 2024). We also included measurements from the Kepler instrument (Koch et al. 1998; Borucki et al. 2010) as its passband probes essentially the same wavelength range as CHEOPS (Deline et al. 2020). Geometric albedos values in the Kepler passband have been reported in several publications (Esteves et al. 2015; Angerhausen et al. 2015; Heng et al. 2021; Morris et al. 2024). We selected values from the works of Esteves et al. (2015) and Morris et al. (2024) as they performed self-consistent analysis accounting for the contribution of thermal emission in the planetary dayside flux¹⁷.

From our analysis, we derived a strong upper limit on the contribution of the reflected light of WASP-18 b in the TESS passband. Despite its significant overlap with the CHEOPS passband (Deline et al. 2020), the latter has a much less stringent constraint with a hint of detection in CHEOPS at 2.6 σ. We explored the scattering properties of the atmosphere that could produce reflected light signals consistent with our computed geometric albedo values in both passbands. Since the dayside of WASP-18 b is likely too hot for condensates to be formed, the only reflective component of the atmosphere would be Rayleigh scattering by molecular and atomic hydrogen as well as helium. We estimated the geometric albedo of the dayside based on the retrieval results shown on the left-hand side of Fig. 11 and the corresponding median temperature profile depicted in Fig. 10. Given the temperature profile and the stated element abundances, we calculated the chemical composition (assuming chemical equilibrium) and the atmospheric absorption coefficients for the chemical species described above and Rayleigh scattering contributions by atomic and molecular hydrogen as well as helium. For the pressure, we chose 0.1 bar, which is the region where the contribution function peak, as shown in Fig. 10. Based on the computed absorption and scattering coefficients we then obtained the wavelength-dependent geometric albedo A _g following the analytic solution derived by Heng et al. (2021). The corresponding stellar and planet fluxes were then finally integrated over the CHEOPS and TESS passbands to estimate the A _g-contribution of the corresponding measurements.

The wavelength-dependent albedo and the passbands are shown in Fig. 13. The atmospheric species that decisively affects the visible wavelength region for a hot planet such as WASP-18 b is the hydrogen anion H⁻. Its very strong continuum absorption tends to produce black-body like spectra with only very little contributions by scattering (see Kitzmann et al. 2018, for example). As the yellow curve in Fig. 13 suggests, the geometric albedo of WASP-18 b would be quite low. The albedo still increases towards lower wavelengths due to the λ⁻⁴-dependence of Rayleigh scattering, such that the contributions by scattering are higher at the blue end of the CHEOPS passband than in the one of TESS. However, with a passband-integrated geometric albedo for CHEOPS of 0.0013, the estimated value remained below the retrieved value of $A_g^{\text{CHEOPS}}=0.027\pm 0.011$.

Since the most important species affecting the shortwave albedo is H⁻ (see Fig. H.2), we performed two additional calculations that are depicted in Fig. 13. In the first we artificially removed H⁻ (blue line). This leads to a strong increase in the geometric albedo, producing a CHEOPS Ag value of 0.09, while the corresponding value for TESS would be 0.017. In the second additional calculation, the chemical-equilibrium abundance of H⁻ was reduced by 98% (orange curve), which yields $A_g^{\text{CHEOPS}}=0.028$ and $A_g^{\text{TESS}}=0.004$, more comparable to the retrieval results.

Besides H⁻, other short-wave absorbers, such as TiO and VO or the alkali metals Na and K can also impact the amount of scattered light and, thus, the geometric albedo. For the case of no H⁻ absorption, we also show the resulting geometric albedos when removing these species as opacity sources in Fig. H.2. As the results suggest, the impact of TiO and VO on the geometric albedo is negligible for the atmosphere of WASP-18 b. The alkali metals, however, clearly have a stronger impact. When removing them, in addition to H⁻, we obtain the analytical limit of 0.75 for the geometric albedo in case of Rayleigh scattering.

However, the impact of these species can only be noticed when H⁻ is not present. In Fig. H.1 of Appendix H we show the impact of the additional shortwave absorbers on A_g as a function of H⁻. These results suggest that as long H⁻ is not removed entirely, the effect of all of these additional species on the geometric albedo is minimal.

A H⁻ abundance lower than the one expected from chemical equilibrium can clearly be caused by non-equilibrium effects, including, for example, photochemistry. Furthermore, our estimated value of A_g is based on a one-dimensional description. In the real three-dimensional atmosphere, the abundance of H⁻ should be expected to change considerably across the visible dayside since its abundance is very strongly temperature-dependent.

Fig. 10 Brightness-temperature spectra derived from the atmospheric retrieval analyses. The circle markers with error bars show the observed occultation depths with uncertainties. The two leftmost circles (red) correspond to CHEOPS and TESS passbands, while the four rightmost black points cover the Spitzer/IRAC channels. The JWST spectroscopic observations are shown in black between 0.8 and 3.0 μm and have been binned down for better visualisation. The coloured solid curves with shaded areas denote the retrieved median spectra and span of the 1σ credible interval for the four model fits (see legend). The diamond markers show the model spectra integrated over the photometric bands. The grey shaded areas highlight the contributions of some relevant species to thermal emission: alkali/H− and TiO/VO in the blue end, and CO/CO2 in the infrared. Right panel: the median (solid curves) and 1σ credible interval (shaded area) of the temperature-profile posterior distributions (same colour coding as in the left panel). The curves on the left edge show the normalized contribution functions that indicate the pressures mainly probed by the observations according to each retrieval analysis. The dashed black line shows the thermal stability curve for TiO2. Other thermal stability curves for species of interest, such as alkali metals (Na2S or KCl), are situated at even lower temperatures around ~1000 K.

Fig. 11 Atmospheric retrieval compositions (each quadrant shows one of the four retrieval runs, labelled in bold text). Top: pairwise and marginal posterior distributions for the abundance free parameters. The quoted values on top of the marginal histograms denote the median and 1σ credible interval for each parameter (span of the central 68% percentiles, also denoted as the shaded areas). The black diagonal dashed lines denote the parameter space values where both parameter are scaled by the same amount (relative to solar abundances). Bottom: volume mixing ratios (VMR) for each retrieval corresponding to the posterior distributions shown above. The shaded areas denote the 1σ credible-interval span in VMR for selected species (see labels at the right of the bottom right quadrant). The solid curves on the left edge of the VMR panels show the contribution functions that indicate the pressures probed by the observations.

Fig. 12 Geometric albedo of several hot and ultra-hot exoplanets as function of their estimated dayside temperature. The detection values (blue) or upper limits (green) are measured in the CHEOPS passband ($=A_g^{\text{CHEOPS}}$) and reported in the following works: Lendl et al. (2020), Hooton et al. (2022), Deline et al. (2022), Brandeker et al. (2022), Scandariato et al. (2022), Parviainen et al. (2022), Demory et al. (2023), Krenn et al. (2023), Pagano et al. (2024), Singh et al. (2024), Demangeon et al. (2024), Akinsanmi et al. (2024), and Scandariato et al. (2024). Measurements in the Kepler passband from Esteves et al. (2015) and Morris et al. (2024) are shown in grey.

Fig. 13 Estimated geometric albedo of WASP-18 b based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The green line shows the geometric albedo with an H− abundance derived by assuming chemical equilibrium (extremely low Ag values, reaching 0.03 at short wavelengths). For the purple line, this abundance is artificially reduced by 98% (yielding Ag values comparable to our results in the CHEOPS and TESS passbands), while for the yellow case the H− continuum is entirely removed as an opacity source. We further explored the effect of removing TiO and VO (red) as well as Na and K (blue). The impact of the former is negligible and the red curve lies hidden behind the yellow model. The latter reveals however more significant scattering spectral features, especially in the CHEOPS and TESS passbands. The strong increase of Ag towards lower wavelengths is caused by the λ−4-dependence of Rayleigh scattering.

7 Conclusion

We performed a joint analysis of the phase-curve and occultation observations of the planet WASP-18 b, covering six passbands from the visible to the mid-infrared with the space-based instruments CHEOPS, TESS, and the four IRAC channels of Spitzer. We included several signals related to the planet in our modelling; namely, the transit, the occultation, the phase signal, the ellipsoidal variations, the Doppler boosting, and the light travel time. We derived new ephemerides with precisions of 1 second and 1.4 millisecond for the time of inferior conjunction and the orbital period, respectively. We measured a planetary radius of R_p = 1.1926 ± 0.0077 R_J , reaching a precision of 0.65% (equivalent to 550 km).

We interpreted a timing inconsistency with the time of superior conjunction from the JWST observation (Coulombe et al. 2023) as a consequence of the slight orbital eccentricity of the planet. Using priors on the eccentricity and the argument of periastron from radial-velocity measurements (Triaud et al. 2010; Nymeyer et al. 2011; Csizmadia et al. 2019), we substantially improved the constraints on the argument of periastron with a value of $\omega ={261.9}_{-1.4}^{+1.3}\mathrm{deg}$ for an eccentricity of e = 0.00852 ± 0.00091. Despite the relative youthfulness of the WASP-18 system, having such a short-period massive planet on an eccentric orbit is peculiar and might require further study to understand the formation and stability of the planetary orbit.

We investigated the occultation depth values across the 15 years covered by our data sets and we were not able to detect either short-term nor long-term variability.

From the upper limits on the nightside flux in the CHEOPS and TESS passbands, we derived a maximum nightside temperature estimate of T _night < 2000 K. Modelling the dayside emission with general circulation models (GCMs) showed that the presence of magnetic-driven friction and super-solar metallicity are necessary to explain our observations. We also fit the occultation depths, including the data from JWST (Coulombe et al. 2023), with an atmospheric retrieval model. We were able to show that new Spitzer data at 4.5 μm reveal super-solar carbon and oxygen abundances, with significant CO/CO₂ absorption, and a steep inverted temperature gradient throughout the planetary atmosphere. Overall, we found that the retrievals provided a good fit to the observations. However, the difference in brightness temperature inferred from the CHEOPS and TESS occultation depths was not captured by any of the models. This could be indicative of missing physics in the retrievals, such as the omission of reflected stellar light or atmospheric variations across latitude and longitude. We note that while GCMs do help in investigating the 3D effects of the physics, they are not strict fits to the data. Ultimately, future phase-curve or dayside-mapping observations will be critical to better constrain the spatial variation of the physical properties of WASP-18 b, allowing us to fine-tune our GCMs and construct more sophisticated retrieval models.

Finally, we explored the potential contribution of reflected light to the dayside brightness in the CHEOPS and TESS passbands. We strongly constrained the geometric albedo in TESS with a 3-σ upper limit of $A_g^{\text{TESS}}<0.017$. In parallel, the CHEOPS passband showed a hint of flux excess at 2.6σ that corresponds to a geometric albedo $A_g^{\text{CHEOPS}}=0.027\pm 0.011$. This result suggests the presence of a very weak reflectivity of the planetary atmosphere. We modelled the scattering properties of WASP-18 b to match our derived A_g values. We quantitatively reproduced this behaviour as the result of Rayleigh scattering in the presence of continuum absorption by H⁻. The abundance of H⁻ needs to deviate from chemical equilibrium to match the geometric albedo in the CHEOPS and TESS passbands, which can be caused by non-equilibrium effects (e.g. photochemistry) or significant temperature variations across the planetary dayside.

Data availability

Raw and detrended light curves are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5 or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A150

Appendix A CHEOPS observations

Fig. A.1 CHEOPS raw light curves. The 37 CHEOPS visits are displayed in coloured data points from top to bottom with their corresponding models shown as black lines. The planet model at the top (continuous black line) is the final model with only astrophysical signals (see Eq. 12).

Appendix B TESS observations

Fig. B.1 TESS raw light curves extracted by the SPOC. The coloured points, blue and green, represent the photometric data of the two orbits of each sector. The black points are the data binned at a cadence of 30 min. The shaded areas highlight the photometric ranges that were discarded due to the presence of remaining systematic trends. The times in BJDTDB of the different cuts are listed chronologically in Table B.1. The light curve of sector 3 corresponds to the SAP flux (second panel) and the other sectors (2, 29, 30 and 69) are all PDCSAP fluxes.

Appendix C Spitzer observations

Appendix D Phase-curve model

We provide details about the expression of the phase-curve model used in this work and reported in Eq 3.

We start by defining the phase-curve signal of the planet as a function of the planetary phase angle α: $$F_{\text{p}}=a_{\text{PC}}\cos \left(\alpha \right)+b_{\text{PC}},$$(D.1)

where a _PC and b _PC are the phase-curve semi-amplitude and a constant offset, respectively. We defined those two parameters in such a way that the phase-curve signal F _p reaches a maximum Fmax when α = 0 deg and a minimum F _min when α = 180 deg, which corresponds to rewriting the previous equation as follows: $$F_{\text{p}}=\frac{1+\cos \left(\alpha \right)}{2}(F_{\text{max}}-F_{\text{min}})+F_{\text{min}}.$$(D.2)

We then used the expression of the phase angle α as a function of the angular position of the planet θ and the orbital inclination i: $$\cos \left(\alpha \right)=-\sin \left(\theta \right)\sin \left(i\right).$$(D.3)

Here we have defined θ in such a way that the inferior conjunction occurs at θ = 90 deg. We note that when the planetary orbit is not viewed edge-on (i.e. i ≠ 90 deg), one has 0deg < α < 180 deg and the phase-curve signal F _p never reaches the values F _min and F _max.

The angular position of the planet can be written θ = ω + ν, with ω and ν being the argument of periastron and the true anomaly, respectively. In order to account for potential phase shift between the superior conjunction (θ = 270 deg) and the phase-curve maximum value (e.g. due to a hotspot offset), we introduce the term Δφ in our expression of θ: $$\theta =\omega +\nu -\mathrm{\Delta }\varphi .$$(D.4)

We chose to subtract the phase shift value, Δφ, to ensure that a positive Δφ induces a positive time delay, which means that the phase-curve maximum will occur after the superior conjunction (e.g. westward hotspot offset).

We obtain the phase-curve model reported in Eq. 3 by merging Eqs. D.4 and D.3 into Eq. D.2.

Appendix E Ellipsoidal variations model

The ellipsoidal variations correspond to the photometric variability created by the rotation of a distorted star, which deformation is caused by the tidal forces of a close-by massive planet.

In the work of Kopal (1959), the light variation due to stellar distortion is derived in Section IV.2 of chapter IV as a series function of tesseral harmonics. considering the tidal lag negligible, the tesseral harmonics reduce to Legendre polynomials and the light variation can be simplified (Eq. IV-2-37 of Kopal 1959) as shown below up to the third harmonic: $$\begin{array}{rl}\delta \mathfrak{L}=& X_2^{(k)}(1+\frac{{\beta }_2}{4})[\frac{v_1^{(2)}}{3}{P}_2\left(n_0\right)-w_1^{(2)}{P}_2\left(l_0\right)]\\ & -X_3^{(k)}(1+\frac{{\beta }_3}{10})w_1^{(3)}{P}_3\left(l_0\right)\\ & +\dots ,\end{array}$$(E.1)

where:

• $X_j^{(k)}$ are functions of the limb-darkening coefficients (k = 2 for the linear law, k = 3 for the quadratic law) detailed in Eqs. IV-2-38 and IV-2-39 of Kopal (1959),

• β _j = [1 + η _j (R_★ )] y from Eq. IV-2-31 of Kopal (1959), where ${\eta }_j\left(R_{\star }\right)\in [j-2;j+1]$ depends on the density distribution within the star (Eqs. II-2-5 and II-2-6 of Kopal 1959 with a ₁ ≡ R_★ ), and y is the gravity-darkening coefficient,

• P_n (x) is the n ^th Legendre polynomials, i.e. $P_2\left(x\right)=\frac{1}{2}(3x^2-1)$ and $P_3\left(x\right)=\frac{1}{2}(5x^3-3x)$,

• n ₀ = cos(i) and l ₀ = cos(ψ) sin(i), where i is the orbital inclination (from the plane of the sky) and ψ is the true anomaly reckoned from the moment of inferior conjunction¹⁸, i.e. ψ = ω + ν - π/2 and l ₀ = sin(ω + ν) sin(i ),

• $v_1^{(2)}$ refers to the rotational distortion and will therefore be ignored (this term is time-independent),

• $w_1^{(j)}$ refers to the tidal distortion, i.e. ellipsoidal variations, and $w_1^{(j)}={\mathrm{\Delta }}_j\frac{M_p}{M_{\star }}{\left(\frac{R_{\star }}{a}\frac{1+e\cos \left(\nu \right)}{1-e^2}\right)}^{j+1}$ from Eqs. IV-2-33 and IV-2-56 of Kopal (1959) (with $a_1\equiv R_{\star }$, and $R\equiv a\frac{1-e^2}{1+e\cos \left(\nu \right)}$) where ${\mathrm{\Delta }}_j=\frac{2j+1}{j+{\eta }_j\left(R_{\star }\right)}$ (Eq. II-1-27 of Kopal 1959).

The two limits reported on η _j (R_★ ) depend on the density distribution within the star (Eqs. II-2-4 to II-2-8 of Kopal 1959). The lower limit j - 2 assumes a constant density distribution throughout the stellar interior. The upper limit j + 1 represents the case where all the mass is concentrated at the center of the star. We follow Morris (1985) and thus Esteves et al. (2013), assuming the latter case to simplify the equations and obtain ${\mathrm{\Delta }}_j=1$ and ${\beta }_j=(j+2)y$.

When discarding the rotational distortion ($v_1^{(2)}=0$) and assuming a point-mass stellar density distribution (η _j (R_★ ) = j + 1), we can rewrite Eq. E.1, describing the light variation due to stellar distortion, as follows: $$\begin{array}{rl}\delta \mathfrak{L}=& X_2^{(k)}(1+y)\frac{M_p}{M_{\star }}{\left[\frac{R_{\star }}{r\left(\theta \right)}\right]}^3\frac{1-3{\sin }^2\left(\theta \right){\sin }^2\left(i\right)}{2}\\ & X_3^{(k)}(1+\frac{y}{2})\frac{M_p}{M_{\star }}{\left[\frac{R_{\star }}{r\left(\theta \right)}\right]}^4\frac{3\sin \left(\theta \right)\sin \left(i\right)-5{\sin }^3\left(\theta \right){\sin }^3\left(i\right)}{2}\\ & +\dots \end{array}$$(E.2)

Here, $r\left(\theta \right)=\frac{a(1-e^2)}{1+e\cos (\theta -\omega )}$ and θ = ω + ν = ψ + π/2.

We can further simplify the expression of the light variation in the case of a circular orbit (e = 0) where one has r(θ) = a. By discarding all but the time-dependent terms (i.e. the terms depending on θ), we can obtain from Eq. E.2 the following expression of $\delta \mathfrak{L}$: $$\begin{gathered} \begin{array}{r}\delta \mathfrak{L}=\frac{M_p}{M_{\star }}{\left(\frac{R_{\star }}{a}\right)}^3[\frac{3}{4}{X}_2^{(k)}(1+y){\sin }^2\left(i\right)\cos \left(2\theta \right) \\ +\frac{3}{4}{X}_3^{(k)}(2+y)\frac{R_{\star }}{a}(1-\frac{5}{4}{\sin }^2\left(i\right))\sin \left(i\right)\sin \left(\theta \right) \\ +\frac{5}{16}{X}_3^{(k)}(2+y)\frac{R_{\star }}{a}{\sin }^3\left(i\right)\sin \left(3\theta \right) \\ +\dots ]\end{array} \end{gathered}$$(E.3)

This expression is the one reported in Eqs. 1–3 of Morris (1985) where the higher order terms (represented by the dots) have been neglected¹⁹.

We can finally write the expression of the light variation of Eq. E.3 in a more compact way: $$\delta \mathfrak{L}=A_{\text{EV}}[\cos \left(2\theta \right)-A_1\sin \left(\theta \right)+A_3\sin \left(3\theta \right)+\dots ],$$(E.4)

where $$A_{\text{EV}}=\frac{3}{4}{X}_2^{(k)}(1+y)\frac{M_p}{M_{\star }}{\left(\frac{R_{\star }}{a}\right)}^3{\sin }^2\left(i\right),$$(E.5) $$A_1=\frac{1}{4}\frac{X_3^{(k)}}{X_2^{(k)}}\frac{2+y}{1+y}\frac{R_{\star }}{a}\frac{5{\sin }^2\left(i\right)-4}{\sin \left(i\right)},$$(E.6) $$A_3=\frac{5}{12}\frac{X_3^{(k)}}{X_2^{(k)}}\frac{2+y}{1+y}\frac{R_{\star }}{a}\sin \left(i\right).$$(E.7)

Here, we retrieve the expressions of the Eqs. 5, 6, 7 and 8 of Section 4.2.3.

Based on the Eqs. IV-2-15 and IV-2-39 of Kopal (1959), we can write the functions of the limb-darkening coefficients $X_j^{(k)}$ for the quadratic limb-darkening law, i.e. for k = 3. For this, we also need to take into account that quadratic LD law used in Kopal (1959) is different from the commonly used law reported in Manduca et al. (1977) and implemented in the Python code batman (Kreidberg 2015). Indeed, the LD law from Manduca et al. (1977) is $\mathcal{I}\left(\mu \right)/{\mathcal{I}}_0=1-u_1(1-\mu )-u_2{(1-\mu )}^2$, whereas the one from Kopal (1959) is $\mathcal{I}\left(\mu \right)/{\mathcal{I}}_0=1-{\upsilon }_1(1-\mu )-{\upsilon }_2(1-{\mu }^2)$. By applying the transformations ${\upsilon }_1=u_1+2u_2$ and υ₂ = -υ ₂ = in Eqs. IV-2-15 and IV-2-39 of Kopal (1959), we obtain: $$X_2^{(3)}=\frac{2}{5}\frac{15+u_1+2u_2}{6-2u_1-u_2}$$(E.8) $$X_3^{(3)}=\frac{1}{7}\frac{35u_1+22u_2}{6-2u_1-u_2}$$(E.9)

from which the expressions of α_EV and β_EV of Section 4.2.3 are derived (see Eqs. 9 and 10). Note that when u ₂ = υ₂ = 0, we get the expression for the linear limb-darkening law reported in Morris (1985) and Esteves et al. (2013).

Appendix F Bond albedo and heat redistribution efficiency

Here, we briefly explain an analytical method detailed in Deline et al. (2022) to estimate the Bond albedo A_B and the heat redistribution ε.

We define the average effective temperature of a planet ${\tilde{T}}_p$ by: $${\tilde{T}}_p^4=\frac{1}{4\pi }{\iint }_{\text{planet}}{T}_p^{4}d\mathrm{\Omega },$$(F.1)

where T _p is the local effective temperature in the planet atmosphere and dΩ describes a surface element of the atmosphere.

Assuming the planet is emitting as a black body at thermal equilibrium (i.e. absorbed flux equals emitted flux), we can write the following relationship with the average effective temperature of the planet: $${\tilde{T}}_p^4=\frac{1-A_B}{4{\sigma }_{\text{SB}}}{\left(\frac{R_{\star }}{a}\right)}^2{\int }_{\lambda =0}^{+\mathrm{\infty }}\mathcal{S}(\lambda ,T_{\star })d\lambda ,$$(F.2)

where A _B is the Bond albedo, σ_SB is the Stefan-Boltzmann constant, R_★ is the stellar radius, a is the semi-major axis, and S(λ, T_★ ) is the flux emission spectrum of the star as a function of its temperature, $T_{\star }$, and the wavelength, λ.

The parametrisation of Cowan & Agol (2011) allows one to estimate the dayside and nightside effective temperatures as functions of the heat redistribution efficiency ε and ${\tilde{T}}_p$: $$T_{\text{day}}={\left(\frac{8-5\epsilon }{3}\right)}^{\frac{1}{4}}{\tilde{T}}_p,$$(F.3) $$T_{\text{night}}={\epsilon }^{\frac{1}{4}}{\tilde{T}}_p.$$(F.4)

Therefore, from these equations, we can compute a relationship between A _B and ε for a given dayside or nightside temperature, and constrain the range of possible values for these two parameters. Usually, the measured and constrained value is the dayside temperature, T _day, from the occultation depth, and we derive the range of values $A_B=f(\epsilon ,T_{\text{day}})$. In the framework of Cowan & Agol (2011), the maximum dayside temperature is reached when A _B = ε = 0, and one obtains $$T_{\text{day, max}}^4=\frac{2}{3{\sigma }_{\text{SB}}}{\left(\frac{R_{\star }}{a}\right)}^2{\int }_{\lambda =0}^{+\mathrm{\infty }}\mathcal{S}(\lambda ,T_{\star })d\lambda .$$(F.5)

Appendix G Detailed Exorad description

The Exorad GCM used the MITgcm dynamical core to solve the hydrostatic primitive equations (e.g. Showman et al. 2009) on a rotating sphere in an Arakawa C-type cubed-sphere grid that spans in the horizontal plane 128 x 64 cells in longitude and latitude, respectively. In vertical direction it follows the formalism established in Carone et al. (2020): in the vertical direction, 41 logarithmically spaced grid cells between 10⁻⁵ bar and 100 bar are combined with six linearly spaced grid cells between 100 bar and 700 bar, resulting in 47 vertical cells.

Following Showman et al. (2009), a fourth-order Shapiro filter with dampening time scale τ _shap = 25 s is used to remove small grid scale noise after each time step from the horizontal velocity fields²⁰. The GCM is stabilised against non-physical gravity wave reflection on top of the modelling domain by implementing a sponge layer between 10⁻⁴ and 10⁻⁵. The zonal horizontal velocity u is damped by a Rayleigh friction term towards its longitudinally averaged mean ū via: $$\frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{t}}=-k(u-\overline{u})$$(G.1)

where t is the time and k is the strength of the Rayleigh friction applied within the sponge layer, depending on pressure p by: $$k=k_{\mathrm{t}\mathrm{o}\mathrm{p}}\cdot \mathrm{m}\mathrm{a}\mathrm{x}{[0,1-{\left(\frac{\mathrm{p}}{{\mathrm{p}}_{\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{g}\mathrm{e}}}\right)}^2]}^2$$(G.2)

The control parameters p _sponge and k _top determine the position and strength of the applied Rayleigh friction in the sponge layer. In this paper, the default values k _top = 20 days⁻¹ and 4 _sponge = 10⁻⁴ bar are used.

To stabilise the deep layers, basal drag is applied to the zonal u and meridional wind velocity v in pressure layers deeper than 400 bar via: $$\frac{\mathrm{d}\mathrm{u}}{\mathrm{d}\mathrm{t}}=-k_{\mathrm{d}\mathrm{e}\mathrm{e}\mathrm{p}}\cdot u,$$(G.3) $$\frac{\mathrm{d}\mathrm{v}}{\mathrm{d}\mathrm{t}}=-k_{\mathrm{d}\mathrm{e}\mathrm{e}\mathrm{p}}\cdot v,$$(G.4)

where the control parameter k _deep is defined as $$k_{\mathrm{d}\mathrm{e}\mathrm{e}\mathrm{p}}=k_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}}\cdot \mathrm{m}\mathrm{a}\mathrm{x}[0,\frac{\mathrm{p}-490\mathrm{b}\mathrm{a}\mathrm{r}}{700\mathrm{b}\mathrm{a}\mathrm{r}-490\mathrm{b}\mathrm{a}\mathrm{r}}]$$(G.5)

with k _bottom = 1 day⁻¹.

The model was run with a dynamical timestep Δ _t = 25 s for 1000 days simulation time to ensure that the temperature structure does not evolve anymore in the modelling domain· Fluxes are recalculated every fourth dynamical timestep.

Performance tests for the sponge layer and basal drag can be found in Carone et al. (2020). Full radiative transfer was established with the expeRT/MITgcm-branch of the ExoRad framework. A more detailed description of the radiative transfer implementation and performance tests can be found in Schneider et al. (2022).

Appendix H Additional geometric albedo calculations

Fig. H.1 Estimated geometric albedo of WASP-18 b based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The top panel shows the geometric albedo with the H− abundance derived by assuming chemical equilibrium and for several species removed as opacity source (note that all curves overlap as H− dominates and other opacity sources are negligible). In the middle panel, H− abundance is artificially reduced by 98%, for which A g values are comparable to our results in the CHEOPS and TESS passbands (i.e. the red and yellow curves overlap, but the effect of removing Na and K starts to appear around 600 nm). In the lower panel, the H− continuum is entirely removed, and the removal of alkali (Na and K) reveals even more strongly the scattering features, while the effect of TiO and VO remains negligible (overlap of red and yellow models). The strong increase of A g towards lower wavelengths is caused by the λ−4-dependence of Rayleigh scattering. In the last case, when Na and K are removed as well, the geometric albedo converges towards its analytical limit of 0.75 for Rayleigh scattering.

Fig. H.2 Estimated geometric albedo of WASP-18 b as a function of the abundance of H−. The A g values are computed based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The H− abundance that best matches our measurements reported in Section 6.3 lies around 2% of the expected chemical equilibrium value (100% on the x-axis), i.e. a reduction of 98%.

Appendix I Joint-fit values of the systematic parameters

Appendix J Acknowledgements

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). A.De. acknowledges the financial support of the SNSF.

P.E.C. is funded by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship, program J4595-N.

LCa and CHe acknowledge financial support from the Österreichische Akademie 1158 der Wissenschaften and from the European Union H2020-MSCA-ITN-2019 1159 under Grant Agreement no. 860470 (CHAMELEON). Calculations were performed using supercomputer resources provided by the Vienna Scientific Cluster (VSC).

B.-O. D. acknowledges support from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00046.

ML acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576.

ABr was supported by the SNSA.

MNG is the ESA CHEOPS Project Scientist and Mission Representative, and as such also responsible for the Guest Observers (GO) Programme. MNG does not relay proprietary information between the GO and Guaranteed Time Observation (GTO) Programmes, and does not decide on the definition and target selection of the GTO Programme.

S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004.

GBr, VSi, LBo, VNa, IPa, GPi, RRa, and GSc acknowledge support from CHEOPS ASI-INAF agreement n. 2019-29-HH.0.

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.

S.G.S. acknowledge support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC).

The Portuguese team thanks the Portuguese Space Agency for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142255.

TWi acknowledges support from the UKSA and the University of Warwick.

This project has received funding from the Swiss National Science Foundation for project 200021_200726. It has also been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grant 51NF40_205606. The authors acknowledge the financial support of the SNSF.

YAl acknowledges support from the Swiss National Science Foundation (SNSF) under grant 200020_192038.

We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2019-107061GB-C61, PID2019-107061GB-C66, PID2021-125627OB-C31, and PID2021-125627OB-C32, from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), from the Centre of Excellence “María de Maeztu” award to the Institut de Ciències de l’Espai (CEX2020-001058-M), and from the Generalitat de Catalunya/CERCA programme.

DBa, EPa, and IRi acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2019-107061GB-C61, PID2019-107061GB-C66, PID2021-125627OB-C31, and PID2021-125627OB-C32, from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), from the Centre of Excellence “María de Maeztu” award to the Institut de Ciències de l’Espai (CEX2020-001058-M), and from the Generalitat de Catalunya/CERCA programme.

CBr and ASi acknowledge support from the Swiss Space Office through the ESA PRODEX program.

ACC acknowledges support from STFC consolidated grant number ST/V000861/1, and UKSA grant number ST/X002217/1.

ACMC acknowledges support from the FCT, Portugal, through the CFisUC projects UIDB/04564/2020 and UIDP/04564/2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, respectively.

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. thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142531.

This work was supported by FCT - Fundaçâo para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 through the research grants UIDB/04434/2020, UIDP/04434/2020, 2022.06962.PTDC.

O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT.

MF and CMP gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18).

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

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

K.W.F.L. was supported by Deutsche Forschungsgemeinschaft grants RA714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets.

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.

PM acknowledges support from STFC research grant number ST/R000638/1.

This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127).

NCSa acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

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

JV acknowledges support from the Swiss National Science Foundation (SNSF) under grant PZ00P2_208945.

NAW acknowledges UKSA grant ST/R004838/1.

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. CHEOPS data analysed in this article will be made available in the CHEOPS mission archive (https://cheops.unige.ch/archive_browser/).

This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).

Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

We thank both referees for their insightful comments that helped improve the quality of this work.

References

All Tables

All Figures

Fig. 1 Top: ellipsoidal variations modelled with and without secondary terms for the best-fit parameters in the CHEOPS passband. The x-axis represents the orbital phase of the planet θ = ν + ω, where ν and ω are the true anomaly and the argument of periastron, respectively. The blue curve (Sine model) is computed as a simple sine wave, i.e. with A 1 = A 3 = 0. The other models include the secondary terms that are computed for several values of the gravity-darkening coefficient y. Bottom: differences between the different EV models and the Sine model. The peak-to-peak difference in amplitude is of the order of 50 ppm for the realistic range of GDC values (0.1 ≤ y ≤ 0.5). The effect of accounting for A 1 and A 3 in the modelling of EV impacts mainly the difference in baseline level between transit and occultation. The baseline flux during transit will be lower than that during occultation (effect similar to a negative nightside flux). If not taken into account, this may lead to overestimating the occultation depth or underestimating the nightside flux from a planet.

In the text

Fig. 2 CHEOPS phase-folded phase curve. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. From top to bottom are show the full phase curve, a zoomed-in version to highlight the phase-curve signal, and the residuals in ppm.

In the text

Fig. 3 TESS phase-folded phase curve. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. From top to bottom are show the full phase curve, a zoomed-in version to highlight the phase-curve signal, and the residuals in ppm.

In the text

Fig. 4 Spitzer phase-folded occultations. Data points are shown in blue, and their 10-min binned counterparts are shown in black with error bars. The mean and 3 σ uncertainty of the model are shown in orange (line and shaded area, respectively) and have been computed from a set of 2000 models randomly drawn from the posterior distribution. The 4 IRAC channels 1, 2, 3, 4, centred at 3.6 μm, 4.5 μm, 5.8 μm, 8.0 μm, respectively, are shown at the top, with their corresponding residuals in ppm at the bottom.

In the text

Fig. 5 Correlation plot of the posterior distribution of the eccentricity e and argument of periastron ω of the orbit of WASP-18 b that match the mid-occultation time $T_{\text{occ}}^{\text{JWST}}$ from Coulombe et al. (2023). We used our posterior distribution on planetary parameters from Table 3, and Gaussian priors on e cos ω and e sin ω from Triaud et al. (2010) in blue, Nymeyer et al. (2011) in orange, and Csizmadia et al. (2019) in green. Overall, the 3 posterior distributions are all consistent with each other within <1.2 σ. This plot made use of the package corner.py (Foreman-Mackey 2016).

In the text

	Fig. 6 Dayside flux observed in the CHEOPS passband over the three seasons of observations (all three panels). The observations covering an occultation are represented in blue (including a phase curve), and the ones only covering transits in pink. The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{CHEOPS}}={211.9}_{-7.5}^{+7.6}\text{ppm}$ from Table 3.
In the text

	Fig. 7 Dayside flux observed in the TESS passband over the 5 sectors: sectors 2 and 3 on the left, sectors 29 and 30 in the center, and sector 69 on the right. The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{TESS}}=340.4\pm 7.0\text{ppm}$ from Table 3.
In the text

	Fig. 8 Dayside flux observed in the Spitzer/IRAC/Channel 2 passband for programmes 50517 (left panel) and 11099 (right panel). The horizontal black line and shaded areas mark the value $F_{\text{max}}^{\text{SPITZER 2}}=3935\pm 25\text{ppm}$ from Table 3.
In the text

Fig. 9 Occultation depths as a function of wavelength. The black points with errorbars represent the measurement from this work (CHEOPS, TESS and Spitzer; 2 leftmost and 4 rightmost points) and from Coulombe et al. (2023) with JWST (0.8-3 μm). The GCM simulations including TiO and VO (Section 6.2) with and without magnetic drag are shown in orange and light brown, respectively. The orange-filled diamonds mark the passband-integrated GCM values with magnetic drag (τfric = 104 s). The retrieval runs (Section 6.3) are shown in blue, pink, purple and green (same colours as in Fig. 10 and Table 6) depending on the data points included in the fit. An inset zoomed-in view of the CHEOPS-to-1.75 μm range is shown in the lower right corner for convenience.

In the text

Fig. 10 Brightness-temperature spectra derived from the atmospheric retrieval analyses. The circle markers with error bars show the observed occultation depths with uncertainties. The two leftmost circles (red) correspond to CHEOPS and TESS passbands, while the four rightmost black points cover the Spitzer/IRAC channels. The JWST spectroscopic observations are shown in black between 0.8 and 3.0 μm and have been binned down for better visualisation. The coloured solid curves with shaded areas denote the retrieved median spectra and span of the 1σ credible interval for the four model fits (see legend). The diamond markers show the model spectra integrated over the photometric bands. The grey shaded areas highlight the contributions of some relevant species to thermal emission: alkali/H− and TiO/VO in the blue end, and CO/CO2 in the infrared. Right panel: the median (solid curves) and 1σ credible interval (shaded area) of the temperature-profile posterior distributions (same colour coding as in the left panel). The curves on the left edge show the normalized contribution functions that indicate the pressures mainly probed by the observations according to each retrieval analysis. The dashed black line shows the thermal stability curve for TiO2. Other thermal stability curves for species of interest, such as alkali metals (Na2S or KCl), are situated at even lower temperatures around ~1000 K.

In the text

Fig. 11 Atmospheric retrieval compositions (each quadrant shows one of the four retrieval runs, labelled in bold text). Top: pairwise and marginal posterior distributions for the abundance free parameters. The quoted values on top of the marginal histograms denote the median and 1σ credible interval for each parameter (span of the central 68% percentiles, also denoted as the shaded areas). The black diagonal dashed lines denote the parameter space values where both parameter are scaled by the same amount (relative to solar abundances). Bottom: volume mixing ratios (VMR) for each retrieval corresponding to the posterior distributions shown above. The shaded areas denote the 1σ credible-interval span in VMR for selected species (see labels at the right of the bottom right quadrant). The solid curves on the left edge of the VMR panels show the contribution functions that indicate the pressures probed by the observations.

In the text

Fig. 12 Geometric albedo of several hot and ultra-hot exoplanets as function of their estimated dayside temperature. The detection values (blue) or upper limits (green) are measured in the CHEOPS passband ($=A_g^{\text{CHEOPS}}$) and reported in the following works: Lendl et al. (2020), Hooton et al. (2022), Deline et al. (2022), Brandeker et al. (2022), Scandariato et al. (2022), Parviainen et al. (2022), Demory et al. (2023), Krenn et al. (2023), Pagano et al. (2024), Singh et al. (2024), Demangeon et al. (2024), Akinsanmi et al. (2024), and Scandariato et al. (2024). Measurements in the Kepler passband from Esteves et al. (2015) and Morris et al. (2024) are shown in grey.

In the text

Fig. 13 Estimated geometric albedo of WASP-18 b based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The green line shows the geometric albedo with an H− abundance derived by assuming chemical equilibrium (extremely low Ag values, reaching 0.03 at short wavelengths). For the purple line, this abundance is artificially reduced by 98% (yielding Ag values comparable to our results in the CHEOPS and TESS passbands), while for the yellow case the H− continuum is entirely removed as an opacity source. We further explored the effect of removing TiO and VO (red) as well as Na and K (blue). The impact of the former is negligible and the red curve lies hidden behind the yellow model. The latter reveals however more significant scattering spectral features, especially in the CHEOPS and TESS passbands. The strong increase of Ag towards lower wavelengths is caused by the λ−4-dependence of Rayleigh scattering.

In the text

	Fig. A.1 CHEOPS raw light curves. The 37 CHEOPS visits are displayed in coloured data points from top to bottom with their corresponding models shown as black lines. The planet model at the top (continuous black line) is the final model with only astrophysical signals (see Eq. 12).
In the text

Fig. B.1 TESS raw light curves extracted by the SPOC. The coloured points, blue and green, represent the photometric data of the two orbits of each sector. The black points are the data binned at a cadence of 30 min. The shaded areas highlight the photometric ranges that were discarded due to the presence of remaining systematic trends. The times in BJDTDB of the different cuts are listed chronologically in Table B.1. The light curve of sector 3 corresponds to the SAP flux (second panel) and the other sectors (2, 29, 30 and 69) are all PDCSAP fluxes.

In the text

Fig. H.1 Estimated geometric albedo of WASP-18 b based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The top panel shows the geometric albedo with the H− abundance derived by assuming chemical equilibrium and for several species removed as opacity source (note that all curves overlap as H− dominates and other opacity sources are negligible). In the middle panel, H− abundance is artificially reduced by 98%, for which A g values are comparable to our results in the CHEOPS and TESS passbands (i.e. the red and yellow curves overlap, but the effect of removing Na and K starts to appear around 600 nm). In the lower panel, the H− continuum is entirely removed, and the removal of alkali (Na and K) reveals even more strongly the scattering features, while the effect of TiO and VO remains negligible (overlap of red and yellow models). The strong increase of A g towards lower wavelengths is caused by the λ−4-dependence of Rayleigh scattering. In the last case, when Na and K are removed as well, the geometric albedo converges towards its analytical limit of 0.75 for Rayleigh scattering.

In the text

	Fig. H.2 Estimated geometric albedo of WASP-18 b as a function of the abundance of H−. The A g values are computed based on the retrieval results in comparison to the passbands of CHEOPS and TESS. The H− abundance that best matches our measurements reported in Section 6.3 lies around 2% of the expected chemical equilibrium value (100% on the x-axis), i.e. a reduction of 98%.
In the text