The tidal deformation and atmosphere of WASP-12 b from its phase curve ⋆ , ⋆⋆

Context. Ultra-hot Jupiters present a unique opportunity to understand the physics and chemistry of planets, their atmospheres, and interiors at extreme conditions. WASP-12b stands out as an archetype of this class of exoplanets, with a close-in orbit around its star that results in intense stellar irradiation and tidal e ff ects. Aims. The goals are to measure the planet’s tidal deformation, atmospheric properties, and also to refine its orbital decay rate. Methods. We performed comprehensive analyses of the transits, occultations, and phase curves of WASP-12b by combining new CHEOPS observations with previous TESS and Spitzer data. The planet was modeled as a triaxial ellipsoid parameterized by the second-order fluid Love number of the planet, h 2 , which quantifies its radial deformation and provides insight into the interior structure. Results. We measured the tidal deformation of WASP-12b and estimated a Love number of h 2 = 1 . 55 + 0 . 45 − 0 . 49 (at 3.2 σ ) from its phase curve. We measured occultation depths of 333 ± 24ppm and 493 ± 29ppm in the CHEOPS and TESS bands, respectively, while the nightside fluxes are consistent with zero, and also marginal eastward phase o ff sets. Our modeling of the dayside emission spectrum indicates that CHEOPS and TESS probe similar pressure levels in the atmosphere at a temperature of ∼ 2900K. We also estimated low geometric albedos of A g = 0 . 086 ± 0 . 017 and A g = 0 . 01 ± 0 . 023 in the CHEOPS and TESS passbands, respectively, suggesting the absence of reflective clouds in the high-temperature dayside of the planet. The CHEOPS occultations do not show strong evidence for variability in the dayside atmosphere of the planet at the median occultation depth precision of 120ppm attained. Finally, combining the new CHEOPS timings with previous measurements refines the precision of the orbital decay rate by 12% to a value of –30.23 ± 0.82 ms / yr, resulting in a modified stellar tidal quality factor of Q ′ ⋆ = 1 . 70 ± 0 . 14 × 10 5 . Conclusions. WASP-12b becomes the second exoplanet, after WASP-103b, for which the Love number has been measured from the e ff ect of tidal deformation in the light curve. However, constraining the core mass fraction of the planet requires measuring h 2 with a higher precision. This can be achieved with high signal-to-noise observations with JWST since the phase curve amplitude, and consequently the induced tidal deformation e ff ect, is higher in the infrared.


Introduction
Ultra-hot Jupiters (UHJs) orbit very close to their host stars and are subjected to immense tidal forces and irradiation which im-⋆ Based on data from CHEOPS guaranteed time observations (GTO) with Program IDs: CH_PR100013, CH_PR100016, and CH_PR330093 ⋆⋆ The CHEOPS photometric time-series data used in this paper are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgibin/qcat?J/A+A/ pact the orbital, atmospheric, and geometric characteristics of the planets.Depending on the stellar properties, the strong tidal interaction may cause the orbits of the planets to become circular and coplanar (zero eccentricity and obliquity), while their rotation rates and orbital periods may become synchronized (tidal locking; Hut 1980).These effects can impact the atmospheric circulation and also result in tidal deformation of the planet's shape in response to the perturbing force (Correia & Rodríguez 2013;Correia et al. 2014).Another consequence of the tidal in-teraction is the shrinkage of the planetary orbit (tidal decay) due to loss of angular momentum to the star if the planetary orbital period is not synchronized with the stellar rotation period.
The intense irradiation received by UHJs results in extremely high dayside atmospheric temperatures such that molecular water (H 2 O) is thermally unstable in favor of atomic hydrogen (Woitke et al. 2018).The hot daysides of UHJs generally favor the atomic form rather than the molecular form (e.g., Si/SiO or Mg/MgH).Other metal elements may appear in their ionized form, such as Na + , K + , but also in the less abundant Ti + and Al + (Helling et al. 2019(Helling et al. , 2021)).Therefore, UHJs are unique laboratories to study the physics and chemistry of planets at extreme conditions, shedding valuable insights into their orbital evolution, atmospheres, and interiors.WASP-12 b stands out as one of the most irradiated UHJs, with an orbital period of only 1.09 days around a G0 star of T eff =6300 K (Hebb et al. 2008).Separated from the star by less than 3 stellar radii, the planet is exposed to strong tidal forces and irradiation, which makes it an attractive target for characterization via transit, eclipse, and phase-curve observations.WASP-12 b is indeed one of the most extensively studied exoplanets, with numerous observations ranging from ultraviolet to infrared wavelengths, revealing remarkable properties of the planet.The planet is inflated with a large radius of 1.9R Jup likely due to high stellar irradiation.Evidence suggests that it is undergoing atmospheric mass loss with gas outflowing toward the star (Fossati et al. 2010).The proximity of its orbit to the Roche limit of its star makes it one of the few exoplanets where the tidal distortion of the planet can be probed from its light curve (Correia 2014;Akinsanmi et al. 2019).It is also the only exoplanet that has been observationally confirmed to be spiraling into the star due to tidal orbital decay (Yee et al. 2019).Secondary eclipse observations of WASP-12 b have resulted in discrepant eclipse depth measurements at various passbands which could be indicative of variability in the dayside atmosphere (Hooton et al. 2019).Previous work on WASP-12 b found evidence for water absorption in its terminator (Stevenson et al. 2014;Kreidberg 2015), whereas the dayside spectrum showed no signs of water (Swain et al. 2013), supporting previous gas-phase modeling results (e.g., Helling et al. 2019).
In this paper, we report on the transit, occultation, and phase curve observations of WASP-12 b by the CHaracterizing ExOplanet Satellite (CHEOPS; Benz et al. 2021).Previous CHEOPS observations (e.g., Lendl et al. 2020;Barros et al. 2022;Deline et al. 2022;Hooton et al. 2022;Ehrenreich et al. 2023) have shown its remarkable photometric capability in characterizing exoplanets.Here, we analyze the CHEOPS observations alongside archival data of WASP-12 to characterize the shape, orbit, and atmosphere of the planet.In Section 2.1, we described the observations obtained by the different instruments, and also the pre-processing of the datasets.We also derived the stellar properties and summarized the theoretical background on planetary tidal deformation.The analyses of the datasets are detailed in Section 3. In Section 4, we use the ellipsoidal planet model to probe the tidal deformation of the planet in the joint phase curve and transit observations.We characterize the atmosphere of the planet in Section 5: modeling its emission spectrum and probing for variability in its dayside atmosphere.In Section 6, we use our CHEOPS transit timing measurements along with published ones to refine the tidal decay timescale of the planet and the stellar tidal dissipation.Finally, we summarize the main results of our work in Section 7.

Observations and system properties
2.1.Observations WASP-12 has been observed by several space-and ground-based instruments, capturing the transit, eclipse, and also the phase curve of the planet, WASP-12 b.In this paper, we analyze new observations from CHEOPS in addition to existing space-based observations from TESS and Spitzer.Although HST/WFC3 observations of WASP-12 b are also available, we excluded them due to poor transit ingress and/or egress coverage that makes it challenging to probe the planetary deformation.

CHEOPS
We obtained 47 visits of WASP-12 spanning 3 observation seasons between 2022-11-02 and 2022-12-24 using CHEOPS as part of the Guaranteed Time Observations (See observation log in Table A.1).The visits, identified by unique file keys, consist of 21 transits, 25 occultations, and half of a phase curve of WASP-12 b, all taken at an exposure time of 60 s.The phase curve visit lasted for 24 hrs, starting before an occultation and ending after transit.The visit durations of the transits and occultations range between 7.1 -12.6 hrs, capturing significant baselines before and after the transits/occultations (see Figs. B.1 and B.2). Due to the short orbital period of the planet, the transit and occultation visits also fortuitously combine to construct a phase curve for the planet.In some cases, consecutive observations of transits and occultations result in half or full phase curves.
Due to the low-Earth orbit of CHEOPS, its line of sight is often interrupted by Earth occultations of the target or spacecraft passages through the South Atlantic Anomaly resulting in data gaps.For our observations of WASP-12, this resulted in light curve efficiencies between 49 and 63%.CHEOPS data are automatically processed by the official Data Reduction Pipeline (DRP version 13;Hoyer et al. 2020) which performs aperture photometry after calibration of the images and correcting for instrumental and environmental effects.The DRP provides light curves extracted with different aperture radii.Point-spread function (PSF) photometry can also be extracted using the PIPE 1 package developed specifically for CHEOPS data.
In our analyses of WASP-12, we used the PSF-extracted light curves which are less sensitive to contamination from background stars (see e.g., Morris et al. 2020;Brandeker et al. 2022;Delrez et al. 2023).The resulting light curves have a lower scatter than the DRP apertures in all visits.Data points that were flagged to have poor photometry (e.g., due to cosmic ray hits or bad pixels) were discarded.We further removed points with high background (BG > 3 times the median background level) where the correlation with the flux becomes nonlinear due to scattered light from the moon or the Earth's limb.Finally, a 15point moving median filter was used to eliminate points >5 times the median absolute deviation (MAD).In total, 987 points were discarded corresponding to 6.2% of the data points across all visits.
The Nadir-locked orientation of CHEOPS, as it orbits around the Earth, causes its field of view to rotate around the target.Combined with the irregular shape of the CHEOPS PSF, this leads to time-variable flux contamination in the aperture that is correlated with the spacecraft's roll-angle.Thus, it is usually necessary to decorrelate against the roll-angle when analyzing CHEOPS data (e.g., Lendl et al. 2020;Morris et al. 2021;Barros et al. 2022).Spacecraft pointing jitter can also result in flux trends that can be accounted for by decorrelating against the X and Y centroid positions of the target PSF on the CCD.Furthermore, CHEOPS observations can feature ramp effects at the beginning of each visit caused by the thermal settling of the telescope as it adjusts to a new target position.The ramp is accounted for by decorrelating the flux against the deviation of telescope tube temperature from the median value ∆T tube (see e.g., Morris et al. 2021;Deline et al. 2022).Data points (∼45) in the first orbit of visits 44, 45, and 47 were removed as they featured a strong, nonlinear increase in the telescope temperature.
We model the systematic trends in each CHEOPS visit using spline decorrelations2 against the roll-angle, background flux, telescope tube temperature, and the X and Y centroid positions.The spline fit is performed simultaneously with the fit of the astrophysical model ( §3.1) and involves successively fitting splines to the residuals of the astrophysical model and then evaluating the likelihood of the joint model.First, we performed a 2D spline fit of the residual against BG and ∆T tube .The resulting residual is then used for another 2D spline fit against the X and Y centroid positions.Since the flux trends with these variables are approximately linear, the 2D spline functions are defined with a single degree and knot in each dimension.Finally, we model the rollangle trend with a 1D cubic spline fit with knots every 18°.

TESS
TESS observed WASP-12 with 2-minute cadence in sectors 20, 43, 44, and 45 with a span of almost 2 years between December 2019 and December 2021.Across the four sectors, TESS observed 74 transits and occultations.Details of the TESS dataset are given in Table A.2.We utilized the Pre-search Data Conditioning Single Aperture Photometry (PDCSAP) light curve data produced by the Science Processing Operations Center (SPOC) pipeline which has been corrected for known instrumental systematics and contamination (Stumpe et al. 2012;Smith et al. 2012).The light curves from these sectors were recently published by Wong et al. (2022) where the transits, occultations, phase curve, and transit timings for WASP-12b were analyzed.We also analyzed these light curves to complement our CHEOPS observations.
The lightcurves were downloaded from the Mikulski Archive for Space Telescope (MAST) archive using the lightkurve python package (Lightkurve Collaboration et al. 2018).Data points flagged by the SPOC pipeline were removed after which a 15-point moving median filter was used to remove points >5×MAD.We separated the light curves into segments by the time of the momentum dumps and removed any strong flux ramps at the beginning of data segments.In fitting the astrophysical model, we simultaneously account for long-term temporal trends in each data segment using a cubic spline with knots every 3 days so as to preserve the phase variation within an orbital period.

Spitzer
We also analyzed archival Spitzer data in the 3.6µm and 4.5µm channels of the InfraRed Array Camera (IRAC) that have already been published (Cowan et al. 2012;Bell et al. 2019).These observations consist of 2 phase curves in each channel acquired in 2010 (PID 70060, PI P. Machalek) and 2013 (PID 90186, PI K. Todorov).We downloaded the data from the Spitzer Heritage Archive3 .The reduction and analysis of these datasets are similar to Demory et al. (2016b), where we modeled the IRAC intra-pixel sensitivity (Ingalls et al. 2012) using a modified implementation of the BLISS (BiLinearly-Interpolated Sub-pixel Sensitivity) mapping algorithm (Stevenson et al. 2012).In addition to the BLISS mapping (BM), our baseline model includes a linear function of the Point Response Function's (PRF) FWHM along the x and y axes, which significantly reduces the level of correlated noise as shown in previous studies (e.g., Demory et al. 2016a,b;Mendonça et al. 2018;Barros et al. 2022;Jones et al. 2022).This baseline model (BM + PRF FWHM) does not include time-dependent parameters.We implemented this instrumental model in a Markov Chain Monte Carlo (MCMC) framework already presented in the literature (Gillon et al. 2012).We included all data described in the paragraph above in the same fit.We ran two chains of 200,000 steps each to determine baseline corrected light curves at 3.6 and 4.5 µm that were used subsequently.
Previous analyses of these datasets by Cowan et al. (2012) and Bell et al. (2019) reported an anomalous phase modulation (with a periodicity of half the planet's orbital period) in the 4.5µm data but not at 3.6µm.If the anomalous modulation is due to the tidal deformation of the planet, Cowan et al. (2012) estimated that the substellar axis would have to be at least 1.5 times longer than the polar axis.However, the effect of such a large deformation is not supported by the transit.Bell et al. (2019) instead proposed that the anomalous signal at 4.5µm may be due to heated CO emission from gas outflowing from the planet.Bell et al. (2019) also found that the measured hotspot offset of the 3.6µm phase curves significantly changed from eastward in 2010 to westward in 2013.As the unexpected phase curve features in these datasets make them difficult to combine, we chose to use only the transit regions of the Spitzer data in our analysis (0.2 d before and after mid-transit).Similar to the CHEOPS and TESS datasets, we remove outlier points >5×MAD with a 15point moving median filter.Details of the Spitzer datasets are given in Table A.3.

Stellar parameters
To facilitate our analysis of the observations, we refined the stellar parameters of WASP-12 (V = 11.5) as shown in Table 1.The spectroscopic stellar parameters (T eff , log g, microturbulence, [Fe/H]) were originally taken from a previous version of SWEET-Cat (Santos et al. 2013;Sousa et al. 2018).The spectroscopic parameters for WASP-12 were estimated with the ARES+MOOG methodology where we used the latest version of ARES 4 (Sousa et al. 2007, 2015) to consistently measure the equivalent widths (EW) of selected iron lines on the spectrum of WASP-12.The list of iron lines is the same as the one presented in Sousa et al. (2008).In this analysis, we used a minimization process to find the ionization and excitation equilibrium to converge on the best set of spectroscopic parameters.This process makes use of the ATLAS grid of stellar model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973).More recently the same methodology was applied on a combined HARPS-N spectrum where we derived completely consistent spectroscopic stellar parameters (T eff = 6301 ± 64 K, log g = 4.26 ± 0.10 dex, and [Fe/H] = 0.18 ± 0.04 dex; Sousa et al. 2021).We also derived a more accurate trigonometric surface gravity using recent GAIA data following the same procedure as described in Sousa et al. (2021) which provided a consistent value when compared with the spectroscopic surface gravity (4.23 ± 0.01 dex).
We determined the radius of WASP-12 using the infrared flux method (IRFM) with an MCMC approach (Blackwell & Shallis 1977;Schanche et al. 2020).We conducted a comparison between observed and synthetic broadband photometry to determine the stellar effective temperature and angular diameter that is converted to the stellar radius with knowledge of the target's parallax.We constructed the spectral energy distributions (SEDs) of WASP-12 using the results of our spectral analysis above and the ATLAS stellar atmosphere models.We then produce synthetic photometry in the Gaia, 2MASS, and WISE passbands that are compared to Gaia G, G BP , and G RP , 2MASS J, H, and K, and WISE W1 and W2 broadband fluxes from the most recent data releases (Skrutskie et al. 2006;Wright et al. 2010;Gaia Collaboration et al. 2022).Using the offset-corrected Gaia DR3 parallax, we obtained the stellar radius that is reported in Table 1.
We used T eff , [Fe/H], and R ⋆ along with their uncertainties to determine the stellar mass M ⋆ and age t ⋆ by employing two different stellar evolution models.In fact, a first pair of mass and age values (M ⋆,1 , t ⋆,1 ) were computed by the isochrone placement algorithm (Bonfanti et al. 2015(Bonfanti et al. , 2016) that interpolates the input values within precomputed grids of PARSEC5 v1.2S (Marigo et al. 2017) isochrones and tracks.A second pair (M ⋆,2 , t ⋆,2 ), instead, was computed via the CLES (Code Liègeois d'Évolution Stellaire Scuflaire et al. 2008) code, which generates the best-fit evolutionary track according to the provided input and the Levenberg-Marquadt minimization scheme (Salmon et al. 2021).We finally merged the two pairs of outcomes after successfully checking their mutual consistency through the χ 2 -based criterion as described in Bonfanti et al. (2021) and we obtained M ⋆ = 1.422 +0.077 −0.069 M ⊙ and t ⋆ = 2.3 ± 0.5 Gyr, consistent with literature values.

Stellar limb darkening
Accurate modeling of stellar limb darkening (LD) is important when analyzing exoplanetary transits in order to obtain unbiased transit parameters, and also to measure higher-order effects such as tidal deformation (Espinoza & Jordán 2015;Akinsanmi et al. 2019).The theoretical LD profile of different stars can be obtained from stellar atmosphere models that compute the stellar intensities as a function of the foreshortening angle µ measured from the limb to the center of the stellar disk.The most widely used stellar libraries for this purpose are the PHOENIX (Husser et al. 2013) and ATLAS (Kurucz 1993) stellar models.Since they are both theoretical models, they may not always provide an accurate representation of the actual stellar intensity profile (see e.g., Espinoza & Jordán 2015;Patel & Espinoza 2022).Indeed, both libraries predict slightly different intensity profiles for the same star in the same passband, making it difficult to select one over the other.However, since these libraries represent our current knowledge of stellar atmospheres, they may still be used to put useful priors on the stellar limb darkening profiles.Obtaining priors from these models can also be beneficial to multipassband observations since the derived priors will ensure that the LD profile in all passbands relates to the same star.
First, we compute stellar intensity profiles for each passband using the LDCU python package 6 which queries both stellar libraries using the stellar parameters (T eff , log g, and [Fe/H]) given in Table 1.Following Claret & Bloemen (2011), the obtained intensity profile computed from each library, and in each passband, is represented by 100 interpolated points (evenly spaced in µ).The uncertainties in the stellar parameters are accounted for by generating 10000 intensity profile samples using random stellar parameter sets drawn from their normal distribution.The median and standard deviation of the profiles are then computed.Thus, for each passband, we have a median intensity profile from each stellar library with error bars at each µ point.
Different parametric limb darkening laws can be adopted to approximate the model intensity profile to derive limb darkening coefficients (LDCs) that can be used in transit analyses.In this work, we adopted the power-2 LD law (Hestroffer 1997) parameterized by the LDCs, c and α.The power-2 law has been shown to be superior to other two-parameter laws in modeling intensity profiles generated by stellar atmosphere models (Morello et al. 2017;Claret & Southworth 2022).The power-2 LD law is only surpassed by the four-parameter law which is difficult to use in transit model fitting due to the higher number of parameters and the strong correlations between them, which can lead to nonphysical intensity profiles.
Similar to Barros et al. (2022), we leveraged the computed model intensity profiles from the two libraries to derive priors on the LDCs.We obtain LDCs in each passband by fitting the preferred power-2 law to the corresponding combined PHOENIX and ATLAS model intensity profiles.The 1σ uncertainties of the obtained LDCs are inflated such that the allowed parameter space encompasses both intensity profiles and associated 1σ uncertainties.This approach is illustrated in Fig. B.3a and the derived LDC priors for the passbands are reported in Table A .4.Using the derived LDCs as priors allows the transit fit to determine the best-fit limb darkening profile without being too restricted to the predictions of either library.Alternative approaches to modeling limb darkening using the intensity profiles are presented in Section 4.2.

The planet: Tidal deformation
WASP-12 b orbits so close to its host star that it is predicted to be one of the most tidally deformed planets (Akinsanmi et al. 2019;Hellard et al. 2019;Berardo & De Wit 2022).The deformation of a planet in response to perturbing forces depends on its interior structure and can be quantified by the second-degree Love number for radial displacement h 2 (Love 1911;Kellermann et al. 2018).For a fluid planet, h 2 is related to the tidal Love number k 2 (as h 2 = 1 + k 2 ) which measures the distribution of mass within the planet.Therefore, measurement of h 2 from the detection of tidal deformation allows constraining the planet's interior structure (Kramm et al. 2011(Kramm et al. , 2012)).The relationship between the tidal response of a planet and its core mass has been investigated in previous studies (e.g., Batygin et al. 2009;Ragozzine & Wolf 2009;Kramm et al. 2011) showing that an incompressible fluid planet with a homogeneous interior mass distribution will have the highest h 2 value of 2.5.However, the value of h 2 decreases as a planet becomes more centrally condensed (more mass at the core), since the presence of a massive core reduces the response of a planet to the perturbing potential (Leconte et al. 2011b).For instance, the lower measured Love number of 1.39 for Saturn (Lainey et al. 2017) reflects a core mass fraction higher than that of Jupiter which has a Love number of 1.565 (Durante et al. 2020).
The shape of a tidally deformed planet can be described by a triaxial ellipsoid (with axes r 1 , r 2 , r 3 ) where r 1 is the planet radius oriented along the star-planet (substellar) axis, r 2 is along the orbital direction (dawn-dusk axis), and r 3 is along the polar axis.The volumetric radius of the ellipsoid is given by R v = (r 1 r 2 r 3 ) 1/3 .According to the analytical shape model formulated in Correia (2014), the axes of the ellipsoid follow the relation 7 : where q is an asymmetry parameter defined as Therefore the axes r 1 , r 2 , r 3 of the deformed planet depend on the Love number h 2 , the proximity of the planetary orbit to the star a, the planet-to-star mass ratio Q M = M p /M ⋆ , and the volumetric radius of the planet R v .
The nonspherical shape of the ellipsoidal planet causes the projected cross-sectional area to vary as it rotates with orbital phase.The projected area as a function of orbital phase angle (ϕ = 2π phase) is given (e.g., by Leconte et al. 2011a) as where i is the orbital inclination of the planet.The projected area can be calculated using Eqs.( 1) and ( 2), given the values of h 2 , R v , and Q M .The ellipsoid projects its maximum area at quadrature, while the minimum area is projected at mid-transit (and mid-occultation).The effective planetary radius at mid-transit (ϕ = 0) can be obtained from Eq. (3) as R eff p = A ϕ=0 /π.The varying planet projection can lead to anomalies in the transit light curve (Correia 2014;Akinsanmi et al. 2019) and also in the shape of the full-orbit phase curve (Leconte et al. 2011b;Akinsanmi et al. 2023) when compared to a spherical planet.Measuring the deviation of the planet's shape from sphericity in high-precision light curves can thus provide a measurement of the Love number.

Light curve analysis
In this section, we describe the analytical light curve model that we use to model the phase curves, transits, and occultations of WASP-12 b.We further describe our fitting methodology for each analysis.

Light curve model
We adopt an analytical light curve model composed of the transit (F tra ) and occultation (F occ ) signals, the phase variation signal by the planet (F p : due to reflection and thermal emission from the atmosphere), and also the phase variation signal by the star (F ⋆ : due to ellipsoidal distortion of the star by the planet and Doppler beaming of the stellar light).The total phase curve model is given as a function of the orbital phase angle as: Given the expected deformed shape of WASP-12 b, an adequate light curve model should account for the deformation in the relevant component signals (Akinsanmi et al. 2023).We describe the components of Eq. ( 4) in the following sections.

Stellar phase variation model
The flux from the star F ⋆ varies as a function of phase as such that F ⋆ is unity at mid-transit and mid-eclipse.
where Q M is a again the planet-to-star mass ratio, a/R ⋆ is the semi-major axis scaled by the stellar radius, i is the orbital inclination, K RV is the radial velocity (RV) semi-amplitude, and c is the speed of light.The coefficient α EV depends on the linear limb darkening coefficient u, and gravity darkening coefficient g as while the coefficient α DB depends on the stellar flux F λ and passband transmission T λ at wavelength λ as

Transit and occultation models
The transit, F tra , and occultation, F occ , signals are generated using the ellc transit tool (Maxted 2016) which allows modeling the planet shape as a sphere or as an ellipsoid parameterized by the Love number as implemented in Akinsanmi et al. (2019).The ellipsoidal planet model parameters are the same as the usual spherical planet model except that the spherical planet radius R p is replaced by the volumetric radius R v8 , and the addition of h 2 and Q M .

Planetary phase variation model
Although the flux from the planet (F p ) is composed of both reflected light and thermal emission from the atmosphere, the degeneracy between both components makes it challenging to model them simultaneously (see e.g., Lendl et al. 2020;Deline et al. 2022;Parviainen et al. 2022).Following recent optical phase curve studies (e.g., Shporer et al. 2019;Wong et al. 2021;Daylan et al. 2021), we model the total planetary phase variation using a sinusoidal function given by: where F max and F min are the maximum and minimum planet fluxes respectively, and δ is the hotspot offset (positive value means an eastward offset).The planet's dayside flux F d (i.e., occultation depth) and nightside flux F n are derived as the value of F p (ϕ) at ϕ = π and ϕ = 0 respectively9 .The semi-amplitude of the atmospheric phase variation A atm is (F max − F min )/2.Following the recommendation of Akinsanmi et al. (2023), we account for the planetary tidal deformation in the outof-transit phases by multiplying the planet's phase variation (Eq.10) by the normalized phase-dependent projected area of the ellipsoid (Eq. 3) to have This implies that the phase variation of the deformed planet depends on h 2 , R v , and Q M .Bell et al. ( 2019) employed a similar approach to model the tidal deformation of WASP-12 b in an attempt to explain its anomalous phase curve shape in the Spitzer 4.5µm band.However, the planet's shape was modeled as a biaxial ellipsoid (instead of the triaxial model used here).They found that even if tidal deformation might contribute to the phase curve, it is not sufficient to explain the observed anomaly present only in the 4.5µm phase curve.Similarly, in the analysis of the HST and Spitzer phase curves of WASP-103 b, Kreidberg et al. ( 2018) accounted for the planetary deformation by including the normalized phase-dependent projected area of the ellipsoid.However, instead of fitting the shape of the planet, it was fixed based on the tabulated predictions in Leconte et al. (2011b).

Light travel time
We corrected for light-travel time across the planetary system by converting the observation times (t obs ) into a reference time (t ref ).The reference time accounts for the projected distance between the current position of the planet along the orbit and its position at inferior conjunction.This is given for a circular orbit as (Deline et al. 2022)

Fitting process
To obtain results regarding different planetary properties, we perform different model fits to the datasets: -We analyzed all datasets by performing a global fit using the ellipsoidal and spherical planet models and comparing the results ( §3.3).This analysis allows to constrain the shape and atmospheric properties of the planet.-We analyzed the CHEOPS transit observations individually ( §3.4) to derive transit timings for orbital decay analysis .-We analyzed the CHEOPS occultation observations individually ( §3.5) to measure the occultation depth of each visit and probe for potential atmospheric variability.
In these fits, the astrophysical and systematic trends are modeled simultaneously.In cases where we perform model comparison, we sample the parameter space using the nested sampling algorithm, dynesty (Speagle 2019) which provides posteriors of the fit and also the Bayesian evidence for each model.We used 1000 live points to explore the parameter space until the estimated log-evidence was smaller than 0.1.In other fits, we use the affine-invariant MCMC ensemble sampler implemented in emcee (Foreman-Mackey et al. 2013).The results of these analyses are discussed in Sections 4, 5 and 6.

Phase curve analysis
We used dynesty to fit the light curve model (Eq.4) to all the datasets simultaneously, considering both a spherical and an ellipsoidal planet.We fit the CHEOPS and TESS observations using the full phase curve model, and the Spitzer observations using only the transit component.We fit for the orbital parameters: planetary period P, mid-transit time T 0 , planet-star radius ratio R p /R ⋆ , scaled semi-major axis a/R ⋆ and impact parameter b.We also fit for the phase curve parameters: F d , F n , δ, A EV , and A DB .Fitting for the total EV amplitude, A EV , ensures that we consider all possible combinations of the component limb and gravity darkening coefficients (in Eqs.6 and 8) and not just their theoretically estimated values.
As seen in Table 2, we assumed wide uniform priors on all parameters except for the LDCs which have Gaussian priors as derived in Section 2.2.2, and A DB for which we used Gaussian priors centered on theoretical values calculated from Eq. ( 9) for each passband with a standard deviation of 0.2.The prior on F n includes negative values to ensure Gaussian posteriors and avoid overestimating the nightside flux.For the fit using the ellipsoidal planet model, we adopt a uniform prior for h 2 and a normal prior for the log of the mass ratio (log Q M ) based on literature values.
We performed the analysis using a two-step fitting process in which an initial fit of the data from each instrument is performed to estimate the noise properties of each visit/sector (e.g., Lendl et al. 2017;Wong et al. 2021;Demory et al. 2023).We estimated the amplitude of additional white noise β w , which is calculated as the ratio of the residual RMS to the mean photometric uncertainty.We also calculated the amplitude of the red noise β r by taking the average of the ratio of the binned residuals at several timescales to the expected Gaussian 1/ √ n noise scaling at that timescale (e.g., Winn et al. 2008).We then rescaled the photometric uncertainties of the light curve by multiplying them by β w β r .This method allows us to propagate the extra noise contributions to the best-fit parameters.Fitting instead for a jitter parameter to account for extra white noise in each visit/sector of the instruments would introduce 54 additional parameters to the global fit of the datasets, making the convergence time of the fit much longer.We chose not to model the red noise using Gaussian Processes since this might be capable of absorbing the subtle signal of tidal deformation in the light curve.
After determining the β w β r for all datasets, the global results are obtained from the final joint fit of the datasets.The orbital parameters (P, T 0 , b, a/R ⋆ ) are common between the datasets while other passband-dependent parameters including R v /R ⋆ are different between the datasets.We also fit the models using fixed quadratic ephemeris determined from our orbital decay result in Section 6, and find that the derived parameters are consistent within 1σ.The results of the spherical and ellipsoidal fits are given in Table 2 while the phase-folded data and best-fit models are shown in Fig. 1 and Fig. 2.
The results of the phase curve fit are discussed in Sections 4.1 and 5.1.

CHEOPS transit analysis
With the aim of deriving precise transit times, we individually analyzed the 22 CHEOPS transit observations using the F tra model in Eq. ( 4).Since the phase variation signals (F p and F ⋆ ) impact the transit baseline of each visit, we first subtracted out the bestfit phase variation signals determined from the phase curve fit (in § 3.3).We propagate the uncertainties of the F p and F ⋆ models by quadratically adding their standard deviations to the flux error bars.Each transit was allowed to have a unique mid-transit time and systematics model, but the shape parameters (i.e., R p /R ⋆ , b, a/R ⋆ ) were constrained using Gaussian priors based on the parameter posteriors from the spherical planet phase curve fit.We sample the parameter space using emcee.A.5 with a mean precision of 27.5 s.

CHEOPS occultation analysis
With the aim of measuring the individual occultation depths, we individually analyzed the 26 CHEOPS occultation observations using the F occ model in Eq. ( 4).Similar to the transit analyses, we also subtracted out the best-fit phase variation signals (F p and F ⋆ ) determined from the phase curve fit and propagated their uncertainties to the flux error bars.Here, we fit for the occultation depth but fixed the orbital and shape parameters (i.e., P, R p /R ⋆ , b, a/R ⋆ ) to the best-fit CHEOPS values obtained from the phase curve fit.However, since precise timing measurements cannot be obtained from the individual occultations due to the low signal-to-noise, we use Gaussian prior on T 0 based on the posterior from the phase curve fit.

Measuring deformation from phase curve
From the joint fit of the CHEOPS, TESS, and Spitzer datasets, we compare the results of the ellipsoidal and spherical planet models.Table 2 reports the median posterior and 1σ uncertainties of the parameters of both models.The posterior probability distributions of some relevant parameters are also shown in Fig. B.5 where we see greater uncertainties in the determination of the planet-star radius ratios for the ellipsoidal model due to strong correlations with the Love number.
From the ellipsoidal planet phase curve fit, we measured a Love number of 1.55 +0.45  −0.49corresponding to a 3.16σ detection.This is the first measurement of the Love number of a planet from the analysis of its full-orbit phase curve.Previous work by Barros et al. (2022) obtained a 3σ measurement of the Love number of WASP-103 b from the analysis of transit-only observations from different space telescopes.Since planetary tidal deformation signal in the phase curve is correlated with the stellar ellipsoidal variation signal, we additionally follow the strategy of Barros et al. (2022) to measure the deformation in the transit-only regions where the stellar ellipsoidal variation is insignificant.We measured h 2 =1.56 +0.47  −0.52 in agreement with the phase curve derived value, but at a slightly reduced significance of ∼3σ.This indicates that the detection of deformation is robust against the modeling of ellipsoidal variation and that the inclusion of out-of-transit data in our phase curve fit slightly enhances the detection of deformation (see Akinsanmi et al. 2023).Both detections of tidal deformation from its induced subtle effects on light curves have been facilitated by CHEOPS, allowing to reach the 3σ measurement significance.We further assess the significance of our detection by computing the Bayes factor from the log-evidence of each model obtained from the dynesty fit.The Bayes factor, B ES is computed as the exponent of the difference in log-evidence between the ellipsoidal and spherical models.We obtained a Bayes factor of 6.7 in positive favor of the ellipsoidal planet model (Kass & Raftery 1995).
Figure 1 shows the best-fit phase curve models to the CHEOPS and TESS light curves and the contribution of the component signals.We see that the major impact of tidal deformation on the total phase curve is the modification of the planet's atmospheric phase variation F p (green curves).Compared to the spherical planet case, F p for the deformed planet features additional flux contribution between transit and eclipse due to the larger and varying projected size of the ellipsoid at these phases.The contribution of tidal deformation to the total ellipsoidal planet phase curve (yellow curve at the bottom of panel b) peaks between quadrature and eclipse (at phases 0.3 and 0.69).The spherical planet model fit attempts to account for the deformation contribution by increasing the amplitude of the stellar ellipsoidal variation F sph EV (dashed black curve).However, since the peaks of F EV occur at quadratures and are always spaced by 0.5 phases, it is unable to completely absorb the deformation signal (which has shorter peak spacing) even if it is shifted in phase.The bottom two panels of Fig. 1 show the residuals from the fits.In each passband, the difference between the ellipsoidal and spherical planet models is overplotted on the spherical planet fit residuals, highlighting the deformation signature -that is, the remaining signal due to deformation that cannot be accounted for by a spherical planet model.The deformation signature is concentrated within the in-transit phases while the out-of-transit phases show only slight variations.This is due to the low amplitude of the planetary phase variation in these optical bands (a few U(0, 0.  Figure 2 shows the transit phases of the fit for all the datasets.We see that the amplitude and shape of the deformation signature are passband dependent, respectively, due to the varying size of the planet and also different limb darkening compensation at each passband.The ratios of the ellipsoidal planet axes are also passband-dependent since the asymmetry between the axes depends on the radius at that passband (Eq.2).Using Eqs. ( 1) and (2) with the derived values from the ellipsoidal planet model fit, we calculate the axis ratios, r 1 :r 2 :r 3 , in the different passbands and report them in Table 2.As opposed to the large r 1 /r 3 ratio of 1.5 estimated in (Cowan et al. 2012) to explain the Spitzer 4.5µm phase curve anomaly, our fit to the datasets results in a more physical ratio of 1.16 in this passband.Therefore, the source of the Spitzer 4.5µm remains unexplained and requires further investigation.
Comparing the posterior parameters between the ellipsoidal and spherical planet model fits, we find that most parameters differ only by ≲2σ due to the limited precision of h 2 obtained from these datasets.Figure 3 compares the physical radius and mean density of the planet in the different passbands, calculated using the posteriors from the fits and stellar parameters in Table 1.The results show that the assumption of sphericity leads to an underestimation of the volumetric radius of WASP-12 b by 4.0-5.2%and an overestimation of the planetary density by ≳12%, in line with theoretical expectations (Leconte et al. 2011b;Burton et al. 2014;Correia 2014).The biases in radius and density are respec-tively only ∼2σ and ∼1σ significant due to the limited precision of the derived h 2 from these datasets and the planetary mass from radial velocity (RV).Therefore, an increase in photometric and RV precisions will lead to more significant parameter biases and the derived properties of the planet if sphericity is assumed (Berardo & De Wit 2022).
For both deformed and spherical model phase curve fits, the measured stellar ellipsoidal variation semi-amplitude A EV is larger in the CHEOPS band compared to TESS.This is due to the passband-dependent limb and gravity darkening parameters in α EV in Eq. ( 8).The measured A EV in both bands can be used to independently estimate the mass of WASP-12 b using Eqs.( 6) and ( 8), and we expect to obtain the same mass estimate in both passbands.First, the linear limb-darkening parameter for each passband was estimated using LDCU, while the gravity darkening coefficient was estimated from the tables in Claret (2021) and Claret (2017) for the CHEOPS and TESS bands, respectively.We use the stellar mass from Table 1 and orbital parameters from the joint fit in Table 2.For the spherical planet model, we derive a planetary mass of 2.34 ± 0.48 M Jup from CHEOPS parameters and 2.09 ± 0.49 M Jup from TESS which are consistent with one another within 1σ.These masses are higher but consistent with the value of 1.47 ± 0.073 M Jup derived from radial velocity (RV; Collins et al. 2017) at 1.2-1.8σ.For the ellipsoidal planet model, we derive 1.89 ± 0.53 M Jup for CHEOPS and 1.49 ± 0.58 M Jup for TESS which are more consistent with the RV value at <0.8σ.Although the discrepancy between the masses from EV and RV in the spherical model case is not significant due to the large uncertainties on EV masses, we see a slight indication that accounting for deformation reconciles both estimates better.
The derived h 2 for WASP-12 b is unexpectedly close to the value of Jupiter despite its higher insolation and mass.This was also the case for WASP-103 b with the measured h 2 of 1.585 (Barros et al. 2022).Theoretical models predict a decrease in the Love number with mass above ∼1 M Jup since more massive objects are more compressible and therefore tend to become more centrally condensed (Leconte et al. 2011b).Lower h 2 values are also expected at higher equilibrium temperatures due to the lower density of the planetary envelope compared to the core (Kramm et al. 2012;Wahl et al. 2021).Nonetheless, our derived h 2 value is consistent with the theoretical models (e.g., Leconte et al. 2011b;Wahl et al. 2021) predicting a maximum h 2 of 1.6 for tidally locked hot Jupiters.It has been suggested that the ellipsoidal planet shape model might not sufficiently account for nonlinear tidal response of the planet causing it to systematically overestimate the Love number (Wahl et al. 2021).However, transit models that account for these are not currently available and will require more parameters to model the planetary shape which will be strongly correlated.
Following the procedure of Buhler et al. (2016), we attempt to estimate the core mass fraction of the WASP-12 b using the measured h 2 value.We found that the 3σ measurement of the Love number does not provide valuable constraints as the result remain consistent with a core mass fraction of 0 and 1.Indeed, Akinsanmi et al. (2023) showed that valuable constraints on the core mass fraction require h 2 precisions higher than 4σ.They also showed that a single JWST phase curve of a target such as WASP-12 b is capable of attaining 17σ measurement of h 2 due to the large phase curve amplitude in the infrared, the reduced effect of limb darkening, and the unparalled precision of the instrument.

Impact of limb darkening
As most of the optical band deformation signature is concentrated within the transit phases, the detection is sensitive to the modeling of the limb darkening profile.In the fit for tidal deformation, we used the method described in Section 2.2.2 to derive priors on the LDCs based on model intensity profiles from the two spectral synthesis libraries (we call this our fiducial analysis).We also explore two alternative methods to leveraging the generated model intensity profiles, which involve simultaneously fitting the model intensity profiles and the transit observations in order to find the best-fit LD profile.
• Alternative-1: This approach merges the PHOENIX and AT-LAS model intensity profiles to create a new joint intensity profile whose 1σ uncertainty at each µ encompasses the 1σ uncertainty of the individual model profiles (see Fig. B.3b).The joint intensity profile is fitted with the power-2 LD law simultaneously with the transit observations so that a joint likelihood of the suggested limb darkening profile and transit model is obtained (for each passband) at each iteration in our sampling process.This approach allows the joint model intensity profile and transit observations to determine the best-fit limb darkening profile.The parameters of the LD law here have uninformative wide priors since their values are constrained during the joint fit.
• Alternative-2: In this approach, a uniform bound is defined that spans the median of both model intensity profiles (see Fig. B.3c).During the likelihood estimation, each µ point in the suggested limb darkening profile that is anywhere within the defined bounds is assigned zero likelihood.Profile points outside the bound are still acceptable but have lower likelihood values (Gaussian) depending on their distance from the bounds.This approach is agnostic to the eventual choice of limb darkening profile as long as its points are within or close to the bounds defined by both theoretical stellar intensity profiles.
We perform the ellipsoidal planet model fits to the data again using the aforementioned limb darkening alternative approaches and compare the results with our fiducial analysis.The result is given in Table 3 where a consistent value of h 2 is derived with the 3 methods, indicating that our fit is robust against the modeling of limb darkening.

Atmospheric characterization
Phase curves provide a wealth of information to characterize the atmosphere of a planet such as the day and nightside temperatures, the longitudinal temperature and reflectivity map, and also the efficiency of heat transport among others (see e.g., Cowan & Agol 2011;Shporer 2017;Parmentier & Crossfield 2018).In this section, we infer the atmospheric properties of WASP-12 b from our CHEOPS and TESS phase curve analyses ( §3.3).Unless otherwise stated, the discussion below is based mostly on the phase curve fit using the ellipsoidal planet model.However, we still report the derived parameters for both planet models in Table 2 and we found them to be in agreement within 1σ.

Phase curve constraints
The results of the phase curve fit for the ellipsoidal and spherical planet model fits (Table 2) reveal a larger planet-to-star radius ratio in the CHEOPS band compared to the TESS band, indicating stronger atmospheric opacity in the bluer CHEOPS band.
We also measure a significantly lower dayside flux (occultation depth) in CHEOPS compared to TESS.The nightside fluxes in both passbands are consistent with zero at <1σ (2σ upper limits of ∼50 ppm).The best-fit phase curves show only marginally significant eastward phase offsets of 9 ± 5°and 6 ± 3°in the CHEOPS and TESS bands, respectively.The low nightside flux and small phase offset in both bands indicate low-efficiency daynight heat redistribution at the atmospheric layers probed by the instruments.This is consistent with the theoretical and observed trend of decreasing phase offset with increasing temperature for ultra-hot Jupiters (Parmentier et al. 2016;Komacek & Showman 2016;Komacek et al. 2017) possibly due to their short radiative timescales (Perna et al. 2012).

Atmospheric modeling
We model the atmosphere of WASP-12 b by performing 1D retrievals on the emission spectrum and computing the forward global circulation model.This will facilitate the proper interpretation of the phase curve, thereby enabling the determination of the relative contribution of both reflection and thermal emission to the observed occultation depths.

1D retrieval
We model the emission spectra of WASP-12 b using the opensource Pyrat Bay framework (Cubillos & Blecic 2021) to constrain its atmospheric properties based on the tabulated occultation depth measurements in Hooton et al. (2019, and references therein).This dataset consists of several ground-based, HST/WFC3, and Spitzer measurements spanning 0.5-8µm.Since reflection is not accounted for in Pyrat Bay, our model only considered the infrared observation (λ > 1.0 µm) to avoid interference from the reflected flux at shorter wavelengths.The retrieved thermal emission spectrum can then be computed including shorter wavelengths to estimate the thermal contribution in the CHEOPS and TESS passbands.
We model the atmosphere of WASP-12 b between 10 2 -10 −9 bar adopting the parametric temperature-pressure (T-P) prescription of Guillot (2010).The parametric model depends on the irradiation temperature T irr , the mean thermal opacity log κ ′ , and the ratio of the visible to thermal opacities log γ.For this analysis, we modeled the composition under thermochemical equilibrium, considering the most relevant neutral and ionic species expected for hot-Jupiter atmospheres.The chemistry was parameterized by the abundance of carbon [C/H], oxygen [O/H], and all other metals [M/H] relative to solar-abundance values.The radiative transfer calculation considered HITEMP and Exo-Mol opacity line lists for H 2 O, CO, CO 2 , CH 4 , HCN, NH 3 , TiO, VO (Rothman et al. 2010;Li et al. 2015;Hargreaves et al. 2020;Harris et al. 2006Harris et al. , 2008;;Polyansky et al. 2018;Coles et al. 2019;Yurchenko 2015;McKemmish et al. 2016McKemmish et al. , 2019)), which were preprocessed using the repack algorithm to extract the dominant line transitions (Cubillos 2017).Additional opacities include the Na and K resonant lines (Burrows et al. 2000), collision-induced absorption from H 2 -H 2 (Borysow et al. 2001;Borysow 2002), Rayleigh opacity from H 2 , H, and He (Kurucz 1970), and H − free-free and bound-free opacity (John 1988).For the stellar spectrum, we used a synthetic PHOENIX spectrum (Husser et al. 2013) according to the stellar properties (Table 1).Finally, the atmospheric Bayesian retrieval employed a differential-evolution MCMC algorithm implemented in Cubillos et al. (2017) to construct posterior distributions for the atmospheric parameters.We tested four retrieval scenarios: with and without the Spitzer 5.8 and 8.0µm data points; and also with and without TiO/VO absorption (since these species may or may not be present in the atmosphere).We found statistically consistent results between all scenarios, thus in the following, we report the results of the retrieval including TiO/VO and all Spitzer observations.
Figure 4 shows the retrieved T-P profile and thermal spectrum, which we extended over the CHEOPS and TESS passbands.The retrieval spectral fit is mostly driven by the HST/WFC3 and the Spitzer 3.6 and 4.5 µm observations.The model fits most of the occultation depths relatively well, although several mutually inconsistent measurements lead to a reduced χ 2 of 2.9.The observations at wavelengths λ < 2 µm are well fit by a ∼3000 K blackbody model, indicating that these bands probe a nearly isothermal region of the atmosphere.The depths at the Spitzer 3.6 µm and 4.5 µm bands imply lower brightness temperatures (∼2600 K) than in the HST/WFC3 band, possibly due to absorption by carbon-bearing species.These measurements drive the retrieved model towards a non-inverted T-P profile since these bands probe the upper atmosphere (see contribution function and T-P profile in the right panel of Fig. 4).Previous work by Stevenson et al. (2014a) and Oreshenko et al. (2017) similarly finds a noninverted T-P profile.The retrieval model underestimates the Spitzer 5.8 µm and 8.0 µm which would require higher brightness temperatures of ∼3400 K, and thus are inconsistent with the rest of the observations.These Spitzer points have larger uncertainties and are therefore not very constraining in the fit.
In all scenarios, we find insignificant contributions from Nbearing species, which are not expected to be prominent at hightemperature atmospheres in thermochemical equilibrium.Likewise, TiO and VO do not seem to have an impact on the retrieved spectrum, probably because the bluer end of the spectrum (where these molecules absorb the most) is nearly isothermal, which decreases the amplitude of spectral features.On the other hand, Cand O-bearing species (CO, CO2, CH4, H 2 O) are more prominent.Figure B.4 shows the retrieval posterior distribution, we find super-solar abundances for both carbon and oxygen while keeping mostly C/O > 1 ratios.For other metals, we find subsolar abundances.These results are consistent with previous retrieval analyses of the WASP-12 b emission spectrum (e.g., Oreshenko et al. 2017;Himes & Harrington 2022).

3D global circulation model (GCM)
We also compare the measured occultation depths with the output of a forward 3D GCM that self-consistently calculates the thermal emission from the planet.Here, we use expeRT/MITgcm as introduced by Carone et al. (2020); Schneider et al. (2022).The expeRT/MITgcm uses a pre-calculated grid of correlated-k opacities, where we employ here the S1 spectral resolution as described in Schneider et al. (2022) and opacity sources for Na, K, CH 4 , H 2 O, CO 2 , CO, H 2 S, HCN, SiO, PH 3 and FeH with collision-induced absorption from H 2 , He and H − and Rayleigh scattering of H 2 and He.We omit VO and TiO because our retrievals showed no evidence of a thermal inversion in the upper atmosphere, which would indicate the presence of these molecules.In post-processing, we calculate the thermal emission of the planet for different orbital phases (including the dayside emission) using a wavelength resolution R=100, where we reduce the spatial resolution of the dayside to 15 degrees in longitude and latitude.During the GCM simulation as well as during post-processing, the abundance of molecules is computed using equilibrium chemistry assuming solar ([Fe/H]=0) metallicity.Thus, we include the effect of thermal ionization of H 2 and H 2 O on our predicted phase curves and dayside emission spectrum.
The dayside emission of the GCM is overplotted in Fig. 4 where we see a good agreement with the retrieval model.In particular, the GCM model correctly predicts that water features are muted in the HST/WFC3 wavelength range (between 1.1 and 1.6µm) due to a combination of thermal ionization and thus dissociation of H 2 O over large parts of the dayside and electron opacities that suppress molecular features in this wavelength range.For larger wavelengths, where Spitzer (and now JWST) is sensitive, molecular features are present again.Although the GCM dayside emission yields smaller amplitude features compared to the best retrieval model, it is still qualitatively in agreement with the available data.
Using the phase-resolved thermal emission spectrum from the GCM, we integrate the flux with the CHEOPS and TESS response functions to obtain the thermal phase variation of the planet in both bands.We overplot the thermal GCM along with the fitted planetary models in Fig. 1.For CHEOPS, the lower amplitude of the thermal GCM compared to the fitted F p models indicates the need for additional flux from reflection.For TESS, the computed thermal GCM is a good representation of the observed atmospheric phase variation, although the maximum flux and hotspot offset are slightly overestimated; as often seen when comparing GCMs to observations (Parmentier & Crossfield 2018;Zhang et al. 2018).For ultra-hot Jupiters like WASP-12b, it also has been recognized that their daysides are so hot that the gas is dominated by atoms and ions and no longer by molecules like H 2 O, TiO, and VO.The locally high gas temperatures, therefore, lead to a moderate ionization of the atmosphere (Arcangeli et al. 2018;Parmentier et al. 2018).While expeRT/MITgcm takes the H 2 O thermal stability into account by calculating a chemical equilibrium model and adjusting abundances accordingly, we did not take into account that this partially ionized flow may couple to the exoplanet's magnetic field.Which processes dominate a magnetic coupling will depend on the planet's magnetic field and the local degree of ionization.A partially ionized gas may couple to the magnetic field and some species may, hence, be transported with the flow (Helling et al. 2021).Such interactions may modify the wind flow to dampen the hot spot offset and change the maximum flux (Rogers & Komacek 2014;Beltz et al. 2022), in particular in the infrared (May et al. 2022).In this work, we opted to not add drag to the GCM to mimic the unknown magnetic field in WASP-12 b and minimize the number of uncertainties.

Albedos and dayside/nightside temperatures
The observed dayside flux in each passband is composed of the thermal emission from the planet and the reflected light by the atmosphere in that passband.The observed flux is thus given (e.g., Esteves et al. 2013;Shporer 2017) as: For a planet with radius R p and semi-major axis a, the reflective contribution is determined by the geometric albedo A g , while the thermal contribution is determined by the emission spectra of the planet and the star F λ (T p ) and F λ (T eff ), respectively.Therefore, interpreting the measured dayside fluxes requires determining the relative contributions of reflection and thermal emission in the CHEOPS and TESS bands.
As seen in the T-P profile in Fig. 4, CHEOPS and TESS probe similar pressure levels in the atmosphere consistent with blackbody dayside temperatures of 2915±25 K and 2821±25 K, respectively.The average measured dayside temperature across both bands is 2868 ± 17 K, in agreement with the effective dayside temperature of 2864 ± 15 K derived in Schwartz et al. (2017).Using the retrieved thermal emission spectra, we calculated the thermal contribution in the CHEOPS and TESS passbands as 205±10 and 480±19 ppm, respectively.Subtracting the thermal contribution from the observed occultation depths, we used Eq. ( 13) to estimate the geometric albedo in the passbands.We find A g = 0.083±0.015 in the CHEOPS band and A g = 0.010±0.023 in the TESS band, indicating that the atmosphere has non-negligible reflectivity in the CHEOPS band but much lower reflectivity in the TESS band where A g is consistent with zero with 2σ upper limit of 0.06.The derived A g values are consistent with the results from the shorter wavelengths of HST/STIS where Bell et al. (2017) found A g < 0.064.
The derived low geometric albedo of WASP-12 b follows the trend of low reflectivity (Ag ≲ 0.2) of UHJs in the optical to near-infrared transition bands (Mallonn et al. 2019), which is supported by the difficulty in forming condensates at such high temperatures (Parmentier et al. 2018;Wakeford et al. 2017).Perhaps unlikely, but the higher albedo in the CHEOPS band may be due to high-temperature condensates (e.g., silicates, Al 2 O 3 , CaTiO 3 ) on the western terminator that have been transported from the nightside.Indeed, Sing et al. (2013) found that the bestfit model to the HST/STIS, HST/WFC3, and Spitzer transmission spectrum of WASP-12 b was Mie scattering by Al 2 O 3 haze but Bell et al. (2017) ruled out this scenario based on low dayside reflectivity measured in the HST/STIS band.They instead favored thermal emission and Rayleigh scattering from atomic hydrogen and helium.
We used Eq. ( 13) to convert the measured nightside fluxes to nightside brightness temperatures and obtained 2σ upper limits of T night ≃ 2000 K in both passbands (Table 2).This implies a large day-night temperature contrast (>45%) as expected for UHJs due to relatively inefficient heat redistribution (Perna et al. 2012;Komacek et al. 2017).The derived high nightside temperature limit is consistent with the expectation of increasing values as a function of stellar irradiation (Keating et al. 2019) as night clouds disperse for highly irradiated planets.

Dayside atmospheric variability
Several authors have reported hints of possible time variability in the dayside atmospheric brightness of WASP-12 b due to discrepant secondary eclipse depth measurements obtained at various wavelength bands (see Fig. 4).For example, analyses of two i ′ -band observations from different ground-based telescopes resulted in an eclipse depth difference of more than 2σ (Hooton et al. 2019)  month of each other revealed eclipse depths that are discrepant with 4.5σ significance (von Essen et al. 2019).There is also a 2.5σ discrepancy between published z ′ -band depth measurements (López-Morales et al. 2010;Föhring et al. 2013) and similar level of disagreement in measurements taken around ∼2µm (see Hooton et al. 2019), and also at Spitzer 3.6µm and 4.5µm (e.g., Cowan et al. 2012;Stevenson et al. 2014b).Given that the reported eclipse depths have come from different instruments and authors, it is possible that some of the observed disagreements could be due to instrument systematics, observing conditions, telluric contamination in ground-based observations, or even in the analysis of the data.However, Komacek & Showman (2019) showed that the hydrodynamic instabilities in hot Jupiter atmospheres can impact the thermal emission leading to variability at the 2% level in eclipse depth measurements.WASP-12 b can particularly show variability due to magnetohydrodynamic effects in the partially ionized atmosphere (Rogers & Komacek 2014).It remains unclear whether the observed disagreements are astrophysical or due to systematics.
Contrary to these results, the analysis of 4 sectors of TESS data by Wong et al. (2022) found no strong evidence of variability between the individual eclipse measurements, although the largest discrepancy between any two measurements is 3.1σ.Given the low signal-to-noise of the individual TESS eclipses, they were only able to place 2σ upper limits of 450 ppm and 80 ppm on orbit-to-orbit and month-long variability, respectively.
We further investigate this potential variability with the CHEOPS occultation observations spanning 3 seasons.The analyses of the 26 occultations have been described in Section 3.5.Figure 5 shows the derived occultation depths and their 1σ uncertainties ordered chronologically.The depths derived from the joint fit of occultations with each season of observation are also shown.The individual occultation depths agree with each other within 1σ and also with the joint fit of each season better than 1.2σ.The largest discrepancy between any two depth measurements is 1.8σ while the joint fits for the seasons are in agreement within 1σ.We search for periodicity in the occultation depths using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) and found the highest peak periodicity of 8.3 d to be nonsignificant with a false alarm probability of ∼20%.Therefore, we do not find strong evidence for variability in the dayside atmosphere of WASP-12 b at the median depth precision of 120 ppm attained by CHEOPS.

Tidal decay
Previous analyses of the transit times of WASP-12 b have clearly shown that the orbit of the planet is decaying due to tidal interaction with the star, causing the planet to lose angular momentum to the star (Maciejewski et al. 2016;Yee et al. 2019;Turner et al. 2021;Wong et al. 2022).The estimate of the planet's decay rate was possible due to the long baseline of available timing measurements.Since our CHEOPS observations further extend the time baseline, we perform a fit to refine the ephemeris and decay rate of the WASP-12 b by combining our transit and occultation timing measurements given in Table A.5.We combined these CHEOPS transit timings with prior transit and occultation timings compiled by Yee et al. (2019) from various authors and TESS timings derived in Wong et al. (2022).
We model the orbital decay of the planet using a quadratic ephemeris model which gives the transit and occultation times (e.g., Turner et al. 2021) as: where T 0 is the reference transit time closest to the middle of the entire time series, E is the transit epoch, and dP/dE is the orbital decay rate from which the time decay rate Ṗ = 1 P dP dE and the decay timescale τ = P/ Ṗ can be derived.We used emcee to simultaneously fit Eq. ( 14) to the transit and occultation times with T 0 , P, and dP/dE as free parameters.The result from the model fit is given as: The decay rate Ṗ derived from our fit including new timing measurements by CHEOPS remains consistent with the estimate of −29.81 ± 0.94 ms/yr obtained in Wong et al. (2022) but improves the precision by 12%.Our revised orbital decay timescale τ is slightly shorter but 13% more precise than the value of 3.16 ± 0.1 Myr derived in Wong et al. (2022).Figure 6 shows the transit and occultation timing deviations of WASP-12 b and the orbital decay model fit after subtraction of the bestfit linear ephemeris model.
The rate at which the planet's orbital energy (E p ) and angular momentum (L p ) are being lost to the star can be calculated (e.g., Yee et al. 2019) as: The energy is then dissipated inside the star as the tidal oscillations are converted into heat.The efficiency of tidal dissipation is quantified by the modified tidal quality factor of the star Q ′ ⋆ which can we derive from Ṗ using the constant lag model of Goldreich & Soter (1966) as: Population studies of stars show that the value of Q ′ ⋆ ranges between 10 5 -10 6.5 for hot Jupiter systems and may extend up to 10 7 for binary star systems (Jackson et al. 2008;Ogilvie 2014).Lower values of Q ′ ⋆ imply more efficient tidal dissipation.Our derived value in agreement with the value of Q ′ ⋆ = (1.75±0.12)× 10 5 derived in Yee et al. (2019).The derived value is on the low end of the range of Q ′ ⋆ values and implies a much higher dissipation rate for WASP-12 than is expected for a main-sequence star (Turner et al. 2021;Yee et al. 2019).Since the efficiency of tidal dissipation depends also on the structure and evolutionary state of the star, one explanation would be that WASP-12 is actually a subgiant star capable of such efficient dissipation due to nonlinear breaking of gravity waves close to the center of the star (Weinberg et al. 2017).However, further stellar modeling reports that the observed characteristics of WASP-12 are consistent with a main-sequence star rather than a subgiant (Bailey & Goodman 2019).Our modeling also confirms a young stellar age of 2.3 Gyr ( § 2.2.1).Tidal decay has not been confirmed for any other exoplanet despite several ultra-hot Jupiter systems having similar planetary and orbital characteristics as WASP-12.Recently, Vissapragada et al. (2022) reported evidence for the tidal decay of Kepler-1658 b orbiting an evolved star and derived a decay rate of 131 +20 −22 ms/yr.However, further observation of this system will be needed to improve the precision of the orbital decay rate.As this is the first reported case of planetary inspiral around an evolved star, its confirmation will give support to the theoretical expectation of planetary engulfment with stellar evolution, which results in the observed dearth of hot Jupiters around evolved stars (Grunblatt et al. 2022).Harre et al. (2023) also reported a 5σ significant measurement of orbital decay for WASP-4 b but found that more observations are required to differentiate between the tidal decay and apsidal precession scenarios.

Conclusion
In this work, we have presented a detailed analysis of WASP-12 b using observations by the CHEOPS spacecraft alongside publicly available data from TESS and Spitzer.We leverage these datasets to put constraints on the shape, atmosphere, and orbital characteristics of the planet.We summarize the main results below: -The large number of CHEOPS visits allowed us to construct a phase curve which we analyzed together with TESS phase curve and Spitzer transits to constrain the tidal deformation and atmospheric properties of the planet.From our global fit to the datasets, we measured a Love number of h 2 = 1.55 +0.45 −0.49corresponding to a 3.16σ detection.This measurement makes WASP-12 b the second planet, after WASP-103 b, where tidal deformation has been significantly detected from the light curve.The 3σ Love number measurements for these planets are still consistent with those of Jupiter despite the strong irradiation of these planets.There is a need to improve the precision of Love number measurements in order to perform comparative interior structure analysis between these highly irradiated, tidally locked planets and the cooler Jupiter.Phase curve observations of such planets with JWST will provide the precisions needed to measure the Love number at ∼12σ significance, which will in turn better constrain their core mass fractions and interior structures (Akinsanmi et al. 2023).
-We significantly measure occultation depths of 333 ± 24 ppm and 493 ± 29 ppm in the CHEOPS and TESS bands, respectively.The nightside fluxes are consistent with zero at both bands.We compared our phase curve with the output of 3D-GCM and found close agreement between them.We measured marginal eastward phase offset of in both bands.We also detected stellar ellipsoidal variation in both passbands, which leads to mass estimates that marginally differ from the RV-derived mass at the 1.2-1.8σlevel if a spherical planet shape is assumed.However, when we account for tidal deformation, the mass estimates agree with the RV value at < 1σ indicating a preference for the ellipsoidal model.-We model the emission spectrum of WASP-12 b using published occultation depth measurements spanning 0.5-8µm to derive the thermal profile of the planet.Our best-fit model indicates that CHEOPS and TESS are probing the same pressure level in the atmosphere of WASP-12 b with an average brightness temperature of ∼2868 K.We found no evidence of temperature inversion in the atmosphere in agreement with the conclusion from Madhusudhan et al. (2011).We additionally estimate the geometric albedo of the planet and find the planet to be more reflective in the CHEOPS band with A g = 0.086 ± 0.017 than in TESS with A g = 0.01 ± 0.023.-Our analysis of the CHEOPS occultations did not show strong evidence of variability in the dayside atmosphere of the planet at the median occultation depth precision of 120 ppm attained by CHEOPS.-Our analysis of the tidal decay of the planet using the new CHEOPS observations refines the orbital decay rate of the planet to −30.13 ± 0.82 ms/yr corresponding to a precision improvement of 12% compared to the latest estimates from Wong et al. (2021).

Figure
Figure B.1 shows the best-fit transit+systematics model overplotted on each transit light curve.It also shows the bestfit transit model overplotted on each detrended light curve and annotated with the obtained timing precision.Finally, the residual of each light curve fit is shown with the measured RMS.The mid-transit times obtained for the transits and their 1σ uncertainties are given in TableA.5 with a mean precision of 27.5 s.
Figure B.2 shows the best-fit occultation+systematics model overplotted on each occultation light curve.It also shows the best-fit occultation model overplotted on each detrended light curve and lastly, the residuals of each light curve fit and its RMS.We further jointly analyzed the occultation observations within each of the three seasons, freely fitting for the midoccultation time and occultation depth.The epoch of the occultation time for each season was set to that of the occultation closest to the center of the time series.The resulting timings for each of the three seasons (labeled S1-S3) are listed in Table A.5. Section 5.4 discusses the measured occultation depths with respect to previously reported hints of atmospheric variability in dayside of WASP-12 b.

Fig. 1 :Fig. 2 :
Fig. 1: CHEOPS (left) and TESS (right) phase curves of WASP-12 b.Panels (a): The total phase curve models overplotted on the detrended and phase-folded data.Panels (b): Zoom of panel (a)showing the different components of the total phase curve.The solid blue curve is the total phase curve model for an ellipsoidal planet, while the red dashed curve is for a spherical planet.The solid and dashed green curves represent the best-fit planetary phase variation component (F p ) for the deformed and spherical planets, respectively.The cyan curve is the computed thermal-only 3D GCM for the planet (see §5.2.2) which shows for CHEOPS that some reflection from the atmosphere is required to reach the amplitude of the green curves.The black curves at the bottom of this panel show the best-fit stellar ellipsoidal variation (F EV ) from the deformed (solid) and spherical (dashed) planet model fits.The yellow curve represents the contribution of tidal deformation to the total phase curve, with circles denoting the peaks.Panels (c) & (d): The residuals of the phase curve fits using the ellipsoidal (blue) and spherical (red) planet models.The difference between the models (ellipsoidal-spherical) representing the deformation signature in each passband is overplotted on the spherical model residuals.

Fig. 3 :
Fig. 3: Physical radius and density for WASP-12 b derived from the spherical and ellipsoidal planet model fits.

Fig. 4 :
Fig. 4: Eclipse spectrum of WASP-12 b.Left: Observed occultation depths of WASP-12 b (gray points) with the CHEOPS and TESS occultation depths are shown as colored diamonds.The best-fit thermal retrieval model and its 1σ uncertainties are overplotted in green while the forward thermal 3D GCM spectrum is plotted in violet.The gray dashed curves show the blackbody model spectra for brightness temperatures of 2600 K, 3000 K, and 3400 K. Right: T-P profile for the retrieval showing a gradual temperature decrease with decreasing pressure depth between 10 −1 -10 −4 bar indicative of no thermal inversion in the atmosphere of WASP-12 b.The contribution function at each passband is also shown indicating that TESS and CHEOPS essentially probe the same pressure levels within the atmosphere.

Fig. 5 :
Fig.5: Occultation depth measurements of individual CHEOPS visits in the three observation seasons.The horizontal lines and shaded regions represent the measured depth from jointly fitting all the visits in each season.The depths measured between seasons are consistent within 1σ.The dotted sinusoidal signal is the best-fit periodic signal (P=8.3) .to the depth measurements across all seasons.

Fig. 6 :
Fig. 6: Deviation in the transit (top) and occultation (bottom) timings of WASP-12b compared to the best-fit linear ephemeris model.The black and green points are published measurements taken from Yee et al. (2019) and Wong et al. (2022), respectively, while the red points are the new CHEOPS timing measurements.

Fig
Fig. B.2: CHEOPS occultation light curves of WASP-12b labeled according to their visit numbers.Left: The best-fit occultation and systematics model is overplotted on the data.Middle: Systematics detrended flux with occultation model overplotted.Right: Residuals after subtraction of best-fit occultation and systematics model.The 30-min bins and the root-mean-square (rms) of each visit are also shown.

Table 1 :
Properties of the star WASP-12 system

Table 2 :
Result of the fit to the CHEOPS and TESS phase curves.We note that N(µ, σ) represents a Gaussian prior with mean µ and standard deviation σ while U(a, b) is a uniform prior between a and b.

Table 3 :
Sensitivity of derived Love number on different approaches of modeling limb darkening.

Table A .
1: CHEOPS observation log of WASP-12.The visit types are either occultation (occ), transit (tra), or phase-curve (PC).β w β r gives the white and red noise correction factor to the flux uncertainties of each visit.

Table A .
4: Derived power-2 LDC priors in the different passbands and the posterior from the ellipsoidal and spherical planet model fits.Article number, page 22 of 29 Akinsanmi et al.: The tidal deformation and atmosphere of WASP-12b Table A.5: Derived mid-transit times for the individual CHEOPS transit observations of WASP-12 b.