TESS Revisits WASP-12: Updated Orbital Decay Rate and Constraints on Atmospheric Variability

After observing WASP-12 in the second year of the primary mission, the Transiting Exoplanet Survey Satellite (TESS) revisited the system in late 2021 during its extended mission. In this paper, we incorporate the new TESS photometry into a reanalysis of the transits, secondary eclipses, and phase curve. We also present a new $K_s$-band occultation observation of WASP-12b obtained with the Palomar/Wide-field Infrared Camera instrument. The latest TESS photometry spans three consecutive months, quadrupling the total length of the TESS WASP-12 light curve and extending the overall time baseline by almost two years. Based on the full set of available transit and occultation timings, we find that the orbital period is shrinking at a rate of $-29.81 \pm 0.94$ ms yr$^{-1}$. The additional data also increase the measurement precision of the transit depth, orbital parameters, and phase-curve amplitudes. We obtain a secondary eclipse depth of $466 \pm 35$ ppm, a $2\sigma$ upper limit on the nightside brightness of 70 ppm, and a marginal $6\overset{\circ}{.}2 \pm 2\overset{\circ}{.}8$ eastward offset between the dayside hotspot and the substellar point. The voluminous TESS dataset allows us to assess the level of atmospheric variability on timescales of days, months, and years. We do not detect any statistically significant modulations in the secondary eclipse depth or day-night brightness contrast. Likewise, our measured $K_s$-band occultation depth of $2810 \pm 390$ ppm is consistent with most $\sim$2.2 $\mu$m observations in the literature.


INTRODUCTION
Since its discovery in 2009 (Hebb et al. 2009), WASP-12b has become one of the most intensively studied exoplanets. This archetypal ultra-hot Jupiter, with a mass of 1.5 M Jup and an inflated radius of 1.9 R Jup (Collins et al. 2017), has a 1.09 day orbit around a 6150 K F-type star (e.g., Stassun et al. 2019). The intense stellar irradiation and the high dayside equilibrium temperature of ≈2600 K have made WASP-12b a prime target for atmospheric study. Extensive spectroscopic and photometric transit and secondary eclipse measurements have been obtained at infrared wavelengths (Campo et al. 2011;Croll et al. 2011;Cowan et al. 2012;Crossfield et al. 2012;Zhao et al. 2012;Sing et al. 2013;Stevenson et al. 2014a,b;Croll et al. 2015;Kreidberg et al. 2015;Bell et al. 2019;Arcangeli et al. 2021) and visible wavelengths (López-Morales et al. 2010;Copperwheat et al. 2013;Föhring et al. 2013;Sing et al. 2013;Bell et al. 2017;Hooton et al. 2019;von Essen et al. 2019). Some highlights from these observations include the detection of water vapor in the upper atmosphere along the day-night terminator  and the very low dayside reflectivity (Bell et al. 2017). In addition, full-orbit infrared phase curves from the Spitzer Space Telescope have provided a glimpse of the planet's global temperature distribution (Cowan et al. 2012;Bell et al. 2019).
WASP-12b is notable in several respects. First, the planet is undergoing atmospheric mass loss. Repeated observations * NASA Postdoctoral Program Fellow of the planet's transit in the near-ultraviolet and in H-alpha have revealed a strongly asymmetric signal with an early ingress, suggestive of a distended circumplanetary disk of stripped material (Fossati et al. 2010;Haswell et al. 2012;Nichols et al. 2015;Jensen et al. 2018). This scenario is supported by the recent analysis of the Spitzer 4.5 µm phase curve, which displays a significant modulation with a period equal to half of the orbital period that could be explained by CO emission from escaping gas flowing from the planet into the host star (Bell et al. 2019).
Another peculiarity of WASP-12b is the apparent variation in its dayside atmospheric brightness. Mutually discrepant secondary eclipse depths have been reported at various optical wavelengths (López-Morales et al. 2010;Föhring et al. 2013;Hooton et al. 2019;von Essen et al. 2019), as well as in the 2 µm region (Crossfield et al. 2012;Zhao et al. 2012;Croll et al. 2015). In light of other purported detections of atmospheric variability in hot Jupiters (e.g., Armstrong et al. 2016;Jackson et al. 2019; but also see Lally & Vanderburg 2022), the question of whether these differences are real or are instead caused by instrumental effects and/or observing conditions has come to the fore. Meanwhile, recent theoretical work on the atmospheric dynamics of hot Jupiters has shown that hydrodynamic instabilities and other transient atmospheric processes may lead to orbit-to-orbit evolution in the global average temperature, the effects of which may be detectable with current and near-future instruments (e.g., Tan & Komacek 2019;Komacek & Showman 2020;Tan & Showman 2020).
Lastly, WASP-12b is the only exoplanet for which there is clear evidence for a gradually shrinking orbital period. Evidence for a change in period was first reported by Maciejewski et al. (2016). Subsequent transit and occultation observations solidified the detection and ruled out apsidal precession and the Rømer effect as explanations for the changing period (Patra et al. 2017;Yee et al. 2020). The best explanation appears to be tidal orbital decay, which has motivated theorists to try to understand the observed rate of decay in terms of tidal dissipation mechanisms within the star (e.g., Weinberg et al. 2017;Bailey & Goodman 2019;Ma & Fuller 2021) From 2019 December to 2020 January, the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observed WASP-12 as part of its ongoing mission to search the entire sky for transiting exoplanets. The light curves provided additional transit timings, which were used to refine the estimate of the planet's orbital decay rate (Turner et al. 2021). Meanwhile, the planetary phase curve was analyzed by Owens et al. (2021) and Wong et al. (2021a), yielding a robust secondary eclipse measurement and broad constraints on the day-night temperature contrast. No evidence for significant orbit-to-orbit variability in the dayside brightness was found in the TESS photometry. TESS has since revisited the WASP-12 system during its ongoing extended mission, providing an additional two and a half months of nearly continuous photometric observations.
In this paper, we present a follow-up study of WASP-12b's orbital decay, phase curve, and potential atmospheric variability, incorporating the latest TESS light curves from 2021 as well as a new ground-based K s -band observation of an occultation. The relevant datasets are described in Section 2. An overview of our light-curve modeling techniques and our results is provided in Section 3. In Section 4, we use our measurements to update the orbital ephemeris and derive constraints on the level of atmospheric variability. Finally, the main findings of our work are summarized in Section 5.

TESS Light Curves
TESS first observed WASP-12 (TIC 86396382) during Sector 20 of the primary mission, which lasted from UT 2019 December 24 to 2020 January 21. The WASP-12 data were obtained with Camera 1 in the two-minute cadence mode. The data were downlinked from the spacecraft and passed through the Science Processing Operations Center (SPOC) pipeline hosted at the NASA Ames Research Center (Jenkins et al. 2016). The pipeline determined the optimal photometric aperture and extracted the light curve, which was deposited in the Mikulski Archive for Space Telescopes (MAST).
Analyses of the Sector 20 light curves of WASP-12 have been published by several teams. Turner et al. (2021) analyzed the transits and secondary eclipses and used the measured timings to calculate a refined orbital decay rate for WASP-12b. The full-orbit phase curve of the WASP-12 system was analyzed independently by Owens et al. (2021) in a standalone paper and by Wong et al. (2021a) as part of a systematic phase-curve study of short-period exoplanet systems observed during the second year of the TESS mission.
TESS revisited the northern ecliptic hemisphere beginning in mid-2021. Due to the fortuitous positioning of the fields of view, WASP-12 was observed in three consecutive sectors: Sector 43 (2021 September 16 to October 12), Sector 44 (2021 October 12 to November 6), and Sector 45 (2021 November 6 to December 2). During these three roughly month-long observations, the target was located on Cameras 4, 3, and 1, respectively. The combined TESS dataset for WASP-12 now spans almost two years.
In the analysis of the TESS WASP-12 photometry presented in this paper, we utilized the pre-search data conditioning (PDC) light curves, which were corrected for common-mode instrumental systematic trends by the SPOC pipeline using empirically derived cotrending basis vectors Stumpe et al. 2012Stumpe et al. , 2014. Our experience with TESS light-curve fitting has demonstrated that the PDC light curves typically exhibit lower levels of timecorrelated noise when compared to the raw simple aperture photometry (SAP) light curves, while preserving genuine astrophysical variability, including the orbital phase curve (see, for example, Wong et al. 2020b,d).
The preprocessing steps that we completed prior to fitting were largely identical to the methodology described in our previous TESS light-curve papers (e.g., Wong et al. 2020cWong et al. , 2021a. First, we removed all flagged points and trimmed 3σ outliers from the transit-masked photometry using a 16point-wide moving median filter. Next, we separated the light curves into data segments that are separated by the scheduled momentum dumps and data downlink interruptions. Lastly, we searched for significant flux ramps that often occur at the beginning or end of data segments and removed them from the time series. The list of TESS light-curve segments is provided in Table 1, including the number of datapoints, time ranges, and any trimmed portions. Figure 1 shows the PDC light curves from the four sectors of the TESS observations, with the momentum dumps and trimmed ramps marked.

Palomar/WIRC Secondary Eclipse
We observed a full secondary eclipse of WASP-12b on UT 2021 February 18 with the Wide-field Infrared Camera (WIRC; Wilson et al. 2003) on the Hale 200 Telescope at Palomar Observatory, California, USA. Images were taken in the K s band (λ eff = 2.13 µm) using a special beam-shaping diffuser for improved duty cycle and guiding stability; see Stefansson et al. (2017) and Vissapragada et al. (2020) for details concerning the observing setup. Following an initial two-point dither near the target to construct a background frame and mitigate short-cadence detector systematics, we obtained 178 exposures with a total integration time of 42 s (14 coadds × 3 s).
After dark-correcting, flat-fielding, and backgroundsubtracting each image using the methods described by Vissapragada et al. (2020), we extracted the photometry of WASP-12 and 10 nearby companion stars using circular apertures ranging in diameter from 10 to 25 pixels (in 1 pixel steps). For our subsequent analysis, we chose a 17 pixel aperture, which was found to minimize the per-point scatter in the resultant WASP-12 light curve.

DATA ANALYSIS
For our analysis of the TESS photometry, we used the Python-based ExoTEP pipeline (e.g., Benneke et al. 2019;Wong et al. 2020a). In this section, we provide a brief a The numbers indicate the TESS sector, spacecraft orbit (two per sector), and segment number, respectively. b Number of datapoints in each data segment before and after removing flagged points, outliers, and flux ramps. c Start and end times of each data segment (BJD TDB − 2,458,000). d Order of the polynomial detrending model used in the full-orbit phase-curve fitting. overview of our fitting methodology for the transits, secondary eclipses, and full-orbit phase curves. We refer the reader to previous studies, such as Shporer et al. (2019) and Wong et al. (2020c), for more detailed descriptions of the light-curve modeling and error analysis.

Transits
We considered all 74 full-transit events in the trimmed TESS light curves that were not interrupted by a momentum dump: 19 in Sector 20, 18 in Sector 43, 18 in Sector 44, and 19 in Sector 45. These transits are labeled sequentially within each sector in Figure 1. To construct each individual transit time series, we selected all datapoints with timestamps within 0.1P of the predicted mid-transit time, where P is the orbital period. There were about 160 selected points on average. This procedure ensured a sufficient out-of-transit baseline, while simultaneously minimizing the effect of the orbital phase-curve variation on the fitted transit depth.
Within the ExoTEP pipeline, transit light-curve modeling is handled by batman (Kreidberg 2015), which is based on the analytical formulation of Mandel & Agol (2002). For every transit, the unconstrained free parameters were the midtransit time T c , planet-star radius ratio R p /R * , impact parameter b, and scaled semimajor axis a/R * . The orbital period P was held fixed at the extrapolated value based on the quadratic transit ephemeris reported by Turner et al. (2021). We assumed a quadratic limb-darkening law with coefficients constrained by Gaussian priors: the median values (u 1 = 0.24, u 2 = 0.31) were calculated by interpolating the coefficients tabulated by Claret (2017) based on the measured properties of the host star WASP-12 (from Version 8 of the TESS Input Catalog; Stassun et al. 2019), and the Gaussian widths were set to 0.05 to account for the uncertainties in the stellar parameters.
Parameter estimation was carried out using the affineinvariant Markov Chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013). In addition to the free parameters listed above, we included a per-point uncertainty parameter σ that was allowed to vary freely to ensure that the resultant reduced χ 2 was at unity. For normalization and systematics detrending, we multiplied the transit model by either a constant multiplicative constant or a linear function of time. We utilized the Bayesian information criterion (BIC) to determine the optimal systematics model on a transit-by-transit basis; the linear trend was preferred in only one case -the second transit of Sector 43 (S43-2).
The results of our individual transit fits are listed in Table A1 in the Appendix. The measured transit depths show a very high level of mutual consistency, with the largest pairwise discrepancy being 2.3σ. Likewise, the transit-shape parameters b and a/R * show little variation across the 74 transits, with the spread in values lying within 2.4 and 2.5σ, respectively. The median 1σ timing uncertainty for the TESS transits is 42 s. Compilation plots of the normalized transit light curves and best-fit models are provided in Figures A1-A4.

Secondary Eclipses
The individual secondary eclipse events are not detected with a sufficiently high signal-to-noise ratio to permit reliable timing measurements. Instead, we combined all full and uninterrupted occultation events within a TESS sector and calculated a single effective mid-occultation time. In an analogous manner to our transit analysis, we selected points within 10% of an orbit from the predicted mid-occultation time. After median-normalizing each secondary eclipse event, we phase-folded the datapoints. For each sector, we held P, R p /R * , b, and a/R * fixed to the best-fit values determined The red dots indicate portions of the photometry that were trimmed prior to fitting due to short-term flux ramps. The vertical green lines denote the scheduled momentum dumps, and each data segment is labeled with the corresponding sector, orbit, and segment numbers. The primary transits that were included in the transit-timing analysis are labeled sequentially in blue.
from the corresponding individual sector full-orbit phasecurve fit (Section 3.3, Table 2).
We fit for the mid-occultation time T occ and secondary eclipse depth D d , along with the uniform per-point uncertainty σ. The epoch of T occ for each sector was set to the secondary eclipse event closest to the midpoint of the time series. The best-fit occultation timings for the four TESS sectors are (in BJD TDB ) T occ,S20 = 2,458,855.5543 ± 0.0029, T occ,S43 = 2,459,485.3066 ± 0.0021, T occ,S44 = 2,459,511.4976 ± 0.0026, T occ,S45 = 2,459,537.6908 ± 0.0020.
(1) The corresponding measured secondary eclipse depths are 534±54, 497±49, 450±51, and 498±46 ppm, respectively. Later in this paper, we will present improved measurements of the occultation depths derived from the full-orbit phasecurve fits. The improvement comes from a more robust treatment of the time-correlated noise and self-consistently accounting for the out-of-eclipse phase-curve modulations as well as the uncertainties in the orbital parameters. We fit the Palomar/WIRC K s -band secondary eclipse light curve with the exoplanet toolkit (Foreman-Mackey et al. 2021). For detrending vectors, we used the extracted flux arrays from the 10 companion stars, airmass, background flux, and the relative centroid offset of WASP-12. These vectors were combined in linear combination with a constant multiplicative offset, and the coefficients were included in the MCMC fit as free parameters along with the mid-occultation time and eclipse depth. We applied Gaussian priors to the remaining orbital parameters based on the measured values from the previously published Sector 20 phase-curve fit in Wong et al. (2021a).
We obtained a mid-occultation time of T occ,WIRC = 2,459,263.74611 ± 0.00090 BJD TDB (2) and an eclipse depth of 2810 ± 390 ppm. The precision of this timing measurement is significantly higher than that of the TESS occultation timings. The normalized and detrended K s -band secondary eclipse light curve is plotted in Figure 2, along with the best-fit occultation model and 1σ confidence region.

Full-orbit Phase-curve Fits
The phase curve consists of the transits, secondary eclipses, and the synchronous brightness modulation of the entire system. The flux variability contains contributions from the changing viewing geometry of the planet's longitudinal atmospheric brightness modulation as it orbits around the host star, as well as the ellipsoidal distortion of both orbiting companions due to the mutual gravitational interaction and the Doppler boosting of the star's spectrum (see the review by Shporer 2017). In tandem with previous TESS phase-curve studies, we used the following simple sinusoidal light-curve model (e.g., Shporer et al. 2019;Wong et al. 2020bWong et al. ,c,d, 2021a: Here, λ t and λ e denote the unit-normalized transit and secondary eclipse light curves. The flux from the planet, F p , is parameterized by the average flux levelf p , the semiamplitude of the atmospheric brightness modulation A atm , and the phase shift in the planet's phase curve δ. The secondary eclipse depth (i.e., dayside hemispheric brightness) is derived from the aforementioned three parameters via D d = f p − A atm cos(π + δ). Similarly, the nightside flux is given by D n =f p −A atm cos δ. The stellar flux F * includes two harmonic amplitudes, A ellip and A Dopp , which correspond to the ellipsoidal distortion and Doppler boosting signals (e.g., Loeb & Gaudi 2003;Faigler & Mazeh 2011, 2015Shporer 2017).
The PDC light curves are affected by low-level residual systematics in the form of long-term temporal trends. To model these systematics, we multiplied the astrophysical phase-curve model by polynomial functions of time, with a separate systematics function defined for each data segment. We carried out preliminary phase-curve fits to each data segment and determined the optimal polynomial order using the BIC. These selected orders are listed in Table 1 and are predominantly 0 or 1.
Each sector's light curve was analyzed separately using a two-step fitting process. All of the orbital parameters (T c , P, b, and a/R * ) and most transit/eclipse/phase-curve parameters (R p /R * ,f p , A atm , δ, and A ellip ) were allowed to vary freely. We accounted for the light-travel time across the system when modeling the secondary eclipses by adding 2a/c ≈ 23 s to the calculated time of superior conjunction. Our model neglects the gradual decrease in the orbital period over the timescale of a TESS sector, because the expected transittiming deviation is only about 6 ms. The limb-darkening coefficients were constrained by the same Gaussian priors as in our transit fits (u 1 = 0.24 ± 0.05, u 2 = 0.31 ± 0.05). Meanwhile, given the small amplitude of the Doppler boosting signal and its degeneracy with the planet's phase-curve offset δ, we followed the technique used by Wong et al. (2020c) and Wong et al. (2021a) and applied a Gaussian prior to A Dopp . The median and width of this prior was derived from the theoretical formulation described by Shporer (2017) and the system parameters published by Stassun et al. (2019) In the first MCMC fit, we allowed the per-point uncertainty σ to vary. The best-fit σ values for the four TESS sectors were 1880, 1790, 1820, and 1780 ppm, respectively. Upon examination of the residuals, we found discernible timecorrelated noise (i.e., red noise). To propagate this red-noise contribution to the best-fit phase-curve parameters and uncertainties, we computed the scatter of the residuals binned at intervals ranging from 20 minutes to 8 hr and calculated the average fractional deviation β from the expected 1/ √ n trend for pure white noise (see Wong et al. 2020c for a full description of this methodology). We then multiplied the best-fit σ  (2017): u 1 = 0.24 ± 0.05, u 2 = 0.31 ± 0.05. c D d and Dn denote the dayside and nightside fluxes, respectively. The dayside flux is equivalent to the secondary eclipse depth. from the first fit by β to obtain the red-noise-enhanced perpoint uncertainty σ r ≡ βσ. For the second phase-curve fit, we fixed the flux uncertainties to σ r and reran the MCMC analysis to arrive at the final set of results. Across the four TESS sectors, β ranged from 1.01 to 1.23, illustrating the small and variable contribution of red noise throughout the light curves.
The results of our individual sector phase-curve fits are listed in Table 2. There is a high level of consistency in the measured astrophysical parameters across the four sectors. The planet-star radius ratios agree at better than the 1.1σ level, and the transit-shape parameters (b, a/R * ) are mutually consistent to well within 1σ. Meanwhile, the secondary eclipse depths D d and atmospheric brightness modulation amplitudes A atm lie within 1.0 and 2.2σ, respectively. Likewise, the planet's phase-curve offsets δ are in excellent agreement, and the nightside fluxes are all consistent with zero. The star's ellipsoidal distortion signal amplitudes show a 1.3σ spread. All in all, there is no sign of statistically significant alterations in the shape of the WASP-12 phase curve throughout the TESS observations. Constraints on atmospheric variability are discussed in more detail in Section 4.2.
Having established that the transits, secondary eclipses, and synchronous photometric modulations of the WASP-12 system are comparable across the four TESS sectors, we proceeded with a joint phase-curve fit to the combined TESS light curve. Just as in our individual sector phase-curve analyses, we assumed the orbital period to be constant throughout the time interval of the TESS observations. Choosing the reference epoch to be the transit event nearest to the midpoint of the full TESS time series, the accumulated transit-time deviation due to the shrinking period is expected to be less than 15 s at the first and last epochs. A deviation of 15 s is shorter than both the 120 s cadence of the TESS observations and the average precision with which the individual mid-transit times are measured (42 s). Prior to fitting, we divided away the best-fit systematics models from the light curve (calculated in the individual sector phase-curve fits) and fixed the per-point uncertainties to the respective red-noise-enhanced σ r values that we computed previously.
The results of our combined phase-curve analysis are listed in the rightmost column of Table 2. The phase-folded light curve and best-fit model are shown in Figure 3. The fitted values yield a more than twofold improvement in precision when compared to the previously published Sector 20 phasecurve fit results in Wong et al. (2021a). The transit-shape parameters were tightly constrained: b = 0.344 ± 0.036, a/R * = 3.061 ± 0.034, and i = 83 • .54 ± 0 • .74. The planetstar radius ratio of R p /R * = 0.11600 ± 0.00045 is among the most precise in the literature. Using the stellar radius R * = 1.75 ± 0.09 R published in Stassun et al. (2019), we find R p = 1.98 ± 0.10 R Jup and a = 0.0249 ± 0.0013 au.
The TESS-band secondary eclipse depth of WASP-12b is 466 ± 35 ppm, and the nightside flux is consistent with zero at 1.6σ (2σ upper limit of 70 ppm). The day-night brightness contrast has a peak-to-peak amplitude of 524 ± 30 ppm. There is a small, marginally significant 6 • .2 ± 2 • .8 eastward shift in the location of maximum dayside brightness. Given the low optical geometric albedo of WASP-12b (Bell et al. 2017;Wong et al. 2021a), the planetary phase curve observed in the TESS bandpass is dominated by the thermal emission from the atmosphere. It follows that the measured phase offset suggests the presence of superrotating equatorial winds that advect heat toward the colder nightside.
The star's ellipsoidal distortion signal has a semiamplitude of 65 ± 14 ppm. Using the theoretical relationship between the ellipsoidal distortion amplitude and the planetstar mass ratio (e.g., Shporer 2017), along with the measured system parameters and coefficients for limb-darkening and gravity-darkening interpolated from the tables in Claret (2017), we arrived at an independent estimate of WASP-12b's mass: 1.96 ± 0.57 M Jup . This value is consistent with the planet mass presented in the discovery paper (1.41 ± 0.10 M Jup ; Hebb et al. 2009) to within 1σ. The close agreement indicates that the photometric modulation at the first  . Top: the combined four-sector TESS phase curve, phasefolded and binned in 10-minute intervals. The red curve shows the best-fit full-orbit phase-curve model. Middle: a zoomed-in view of the phase-curve variability. The solid and dashed blue curves denote the individual components of the phase curve that correspond to atmospheric brightness modulation and ellipsoidal distortion, respectively. Bottom: the residuals from the best-fit model. harmonic of the orbital frequency is adequately modeled by the classical formalism of stellar tidal distortion.
A recent radial velocity (RV) analysis of the WASP-12 system demonstrated that the apparent signature of a small orbital eccentricity could be explained instead by the flow velocity of the stellar atmosphere due to the tidal bulge raised by the orbiting planet (Maciejewski et al. 2020b). The measured tidal RV amplitude was 7.3 m s −1 , which is expected to produce a photometric modulation with a semi-amplitude of roughly 80 ppm, in good agreement with our measured ellipsoidal distortion signal.

Orbital Decay
Previous transit-timing analyses (Patra et al. 2017;Maciejewski et al. 2018;Yee et al. 2020;Turner et al. 2021) uncovered strong evidence for a shrinking orbital period. Having obtained a large number of new transit timings from TESS Sectors 43-45 and extended the overall time baseline by almost 2 years, we carried out an updated ephemeris fit to refine the orbital decay-rate estimate.
We combined the TESS transit timings listed in Table A1 and the five mid-occultation timings presented in Equations (1) and (2) with the full list of prior transit and occultation times compiled by Yee et al. (2020). 1 There are a total of 213 transit timings and 24 occultation timings. The transitand occultation-time models for orbital decay are The variable dP/dE represents the decay rate in units of days per orbit, T 0 is the zeroth epoch transit time, and E is the transit epoch. Before fitting the measured timings with emcee, we subtracted the light-travel time of 23 s from the mid-occultation times. We designated the transit event closest to the weighted average of timings as the zeroth epoch in our fit. The results of our MCMC analysis are T 0 = 2,457,103.283654 ± 0.000032 BJD TDB , P = 1.091419370 ± 0.000000020 days, dP dE = (−1.031 ± 0.033) × 10 −9 days orbit −1 , dP dt = −29.81 ± 0.94 ms yr −1 .
The updated orbital decay timescale is τ = P/|Ṗ| = 3.16 ± 0.10 Myr, which lies in between the published estimates of Yee et al. (2020) and Turner et al. (2021). Figure 4 shows the timing deviations of the WASP-12b transits and secondary eclipses relative to a constant-period ephemeris model, with the best-fit orbital decay model and corresponding 1σ confidence regions plotted in blue. The fit quality is excellent, with a reduced χ 2 of 1.05 and 234 degrees of freedom (237 datapoints and 3 free parameters). The additional time observations from the TESS extended mission improved the precision of the measured orbital decay rate by almost 30% when compared to the estimate by Turner et al. (2021), which incorporated timings from Sector 20 only. The magnitude of our derived decay rate (29.81 ± 0.94 ms yr −1 ) is 1.5σ larger than the previously published value (27.3 ± 1.4 ms yr −1 ). 2 To assess the influence of the occultation timings on the measured orbital decay rate, we carried out an analogous fit without the occultation timings and obtained 30.01 ± 0.96 ms yr −1 . This indicates that while the occultations are crucial for distinguishing between the orbital decay and apsidal precession scenarios (Yee et al. 2020;1 The previous transit and secondary eclipse timings were contributed by  (Table A1 and Section 3.2). The red datapoint is based on the Ks-band Palomar/WIRC secondary eclipse observation. The green curves and shaded areas indicate the best-fit orbital decay ephemeris model and corresponding 1σ confidence regions. Turner et al. 2021), the timings do not have a substantial impact on the decay-rate estimate. The measured orbital decay rate dP/dt can be used to constrain the host star's modified 3 tidal quality factor Q * , which in the model of Goldreich & Soter (1966) is related to the period derivative via Using the planet-star mass ratio q ≡ M p /M * from Collins et al. (2017) and a/R * from our combined phase-curve fit (Table 2), we arrived at Q * = (1.50 ± 0.11) × 10 5 . As has been discussed in the previous orbital decay studies (Yee et al. 2020;Turner et al. 2021), the measured Q * value of WASP-12 lies at the lower end of most other estimates in the literature based on theoretical calculations or less direct empirical methods. The range of tidal quality factors inferred from analyzing the eccentricity distribution of main-sequence stellar binaries is 10 5 -10 7 (e.g., Ogilvie & Lin 2007;Lanza et al. 2011;Meibom et al. 2015), while analogous population studies of hot-Jupiter hosts point toward Q * values in the range of 10 5.5 -10 6.5 (e.g., Jackson et al. 2008;Husnoo et al. 2012;Barker 2020). Meanwhile, the theoretically expected values of Q * are generally higher (10 7 -10 10 ; see Ogilvie 2014 and references therein). Weinberg et al. (2017) proposed that WASP-12 is a subgiant star for which more rapid tidal dissipation might be expected, due to breaking of internal gravity waves near the star's core. This hypothesis was investigated further by Bailey & Goodman (2019), who agreed with the theoretical premise but found that the observed characteristics of WASP-12 made it more likely to be a main-sequence star. Barker (2020) also calculated tidal dissipation rates for WASP-12 and demonstrated that the theoretical and observed decay rates could only be brought into agreement if the star is a subgiant. Taking a different track, Millholland & Laughlin (2018) proposed that the orbital decay is driven by tidal dissipation within the planet, rather than the star; this would be due to a hypothetical nonzero planetary obliquity that is maintained by a nearby planetary companion. This intriguing scenario was investigated further by Su & Lai (2021), who concluded that it required too much fine tuning to be plausible.
In short, there is not yet a completely satisfactory explanation for the rapid orbital decay of WASP-12b. Some priorities for future progress in this realm are better characterization of the host star and its evolutionary state, as well as more sensitive searches for nearby planetary companions.
Orbital decay has not yet been detected for any other planet, at least not with nearly the same high confidence level as it has been detected for WASP-12b. However, looking at the wider population of hot Jupiters, there are several other planets with characteristics similar to those of WASP-12b that might be expected to show comparable rates of or-bital decay. In particular, if we assume Equation (9) is correct and that Q * is the same for all hot-Jupiter systems, then WASP-103b, KELT-16b, WASP-18b, and the recently discovered ultra-hot Jupiter TOI-2109b have larger predicted |dP/dt| values than WASP-12b (see Figure 14 in Wong et al. 2021b). Notably, the hosts of all five aforementioned planets are F-type stars. The transits of WASP-18b have been monitored for over 15 years, with no sign yet of a decaying orbital period; this nondetection constrains WASP-18's tidal quality factor to Q * > 3.9 × 10 6 at 2σ (Maciejewski et al. 2020a). Meanwhile, a recent analysis of WASP-103b transit timings spanning seven years yielded a 3σ lower limit of Q * > 1.6 × 10 6 ( Barros et al. 2022). Perhaps there is a diversity in the dissipative properties of planet-hosting stars, or a very sensitive dependence on the forcing frequency and internal structure of the star.
TOI-2109b is a particularly promising candidate for probing tidal orbital decay. Due to its ultra-short 0.67 day orbital period and its large mass of roughly 5 M Jup , this extreme hot Jupiter has a predicted orbital decay rate that is more than an order of magnitude faster than that of WASP-12b (assuming the same Q * ). Follow-up transit observations of the TOI-2109 system over the coming years should soon allow us to determine if the host star has a Q * value comparable to that of WASP-12, or a higher Q * similar to WASP-18.

Constraints on Atmospheric Variability
As was discussed in the Introduction, several prior secondary eclipse observations of WASP-12b have yielded discrepant depth measurements, raising the possibility of timevarying behavior of the planet's atmosphere. However, most of these previous results were obtained using different ground-based facilities, and as such it is difficult to ascertain the extent to which variable observing conditions and photometric extraction methodologies may have caused the formal uncertainties in the occultation depths to be underestimated. Wong et al. (2021a) fit the individual WASP-12b secondary eclipses from the TESS Sector 20 light curve and found that the measured depths show no sign of variability to within the measurement precision. We expanded this analysis to all four TESS sectors. After dividing out the best-fit systematics model from each sector's light curve, we fixed the mideclipse times based on the best-fit decaying transit ephemeris from Equation (8). The parameters b, a/R * , and R p /R * were fixed to the best-fit values computed from our combined phase-curve fit (Table 2), and the per-point photometric uncertainties were set to the respective red-noise-enhanced σ r values. No additional systematics modeling was included. Figure 5 shows the depths measured from the 76 secondary eclipse light curves, arranged sequentially by sector. The occultation depths derived from the individual sector phasecurve fits are also shown. Due to the low signal-to-noise ratio of the eclipse light curves, the precision of the individual depth measurements is poor. Nevertheless, we can still place broad constraints on any orbit-to-orbit atmospheric variability that may be present. The majority of individual eclipse depths agree with the corresponding sector-wide values at better than the 1σ level. The largest discrepancy between a pair of eclipse depths is 3.1σ. The overall scatter in the measurements (223 ppm) is smaller than the median uncertainty (242 ppm), indicating that the individual eclipse depths are consistent with being drawn from the same dis-tribution. We calculated the Lomb-Scargle periodogram for the eclipse measurements from the three consecutive sectors (43-45) and detected no significant periodicities. We are able to place a 2σ upper limit of roughly 450 ppm on the orbitto-orbit variability of WASP-12's dayside brightness in the TESS bandpass.
We now compare this constraint to the previously reported instances of variable secondary eclipse depths. The most significant discrepancy in the literature lies between two Vband observations in von Essen et al. (2019): 280 ± 100 and 1160 ± 170 ppm. This 4.5σ difference greatly exceeds our upper limit; moreover, given the lower planet-star contrast ratio in the V -band, such a discrepancy in secondary eclipse depth would translate to an even larger absolute change in the TESS-band dayside brightness. It is important to mention that a third V -band eclipse light curve from the same instrument yielded a nondetection (−30 ± 170 ppm). Therefore, barring the serendipitous observation of a hitherto unexplained astrophysical event, we suggest that uncorrected instrumental systematics and/or observing conditions led to significant biases in one or both of the reported V -band secondary eclipse depths in von Essen et al. (2019).
The remaining instances of short-wavelength (<1 µm) secondary eclipse depth discrepancies in the literature are smaller in magnitude. Hooton et al. (2019) measured i -band depths of 440 ± 210 and 970 ± 140 mmag using two different telescopes. The former value is in good agreement with our TESS-band sector-wide depth measurements. The 2.1σ difference between the two depths is only marginally larger than our 2σ upper limit. Lastly, the three z -band eclipse depths published in López-Morales et al. (2010) and Föhring et al. (2013) are 680 ± 210, 1000 ± 230, and 1300 ± 130 ppm, with a 2.5σ mutual discrepancy between the largest and smallest values. In light of the relatively low statistical significance of these mismatches, as well as the lack of uniformity in instrumentation and possible i -and z -band telluric contamination in ground-based observations, we should be cautious about interpreting these measurements as evidence for atmospheric variability.
The TESS light curves provide much tighter constraints on longer-term variability of WASP-12b's atmosphere. The phase-curve results that we obtained for the individual sectors offer four independent snapshots of the planet's properties, each averaged over a roughly one-month period. As illustrated in Figure 5, the best-fit secondary eclipse depths for the four sectors are in excellent agreement, with a standard deviation of 37 ppm and a median uncertainty of 80 ppm. Moreover, there is no evidence of a systematic shift in the size of the occultation across the roughly 600 days that elapsed between Sector 20 and Sectors 43-45. We place a 2σ upper limit of ≈80 ppm on the TESS-band dayside atmospheric variability of WASP-12b across both month-long and year-long timescales. In other words, the planet's dayside brightness is invariable to within ≈17%. Likewise, there are no significant variations in the day-night temperature contrast, with the atmospheric brightness amplitudes (A atm ) from the four TESS sectors consistent to within 2.2σ (see Table 2).

Sector 20
Sector 43 Sector 44 Sector 45 Figure 5. Plot of the individual secondary eclipse depths measured from the TESS light curve (black). The datapoints are arranged sequentially by epoch, with occultations that coincided with momentum dumps removed. The colored lines and shaded areas indicate the median eclipse depths and 1σ confidence regions that we derived from full-orbit phase-curve fits for each TESS sector (see Table 2). The scatter in the individual eclipse depths is comparable to the typical uncertainties, and the depths measured for each sector are mutually consistent at better than the 1σ level. Zhao et al. (2012) obtained two secondary eclipse measurements in the K s band with the TIFKAM instrument on the Hiltner Telescope: 2810 ± 850 and 3160 ± 950 ppm, with a joint depth of 2990 ± 640 ppm. Narrowband 2.312 µm observations at the Subaru Observatory yielded a notably deeper occultation depth of 4500 ± 600 ppm (Crossfield et al. 2012). Croll et al. (2015) used the Canada-France-Hawaii Telescope/Wide-field InfraRed Camera (CFHT/WIRFCam) instrument on three nights and measured K s -band secondary eclipse depths of 2840 ± 200, 2890 ± 180, and 2590 ± 420 ppm; the same team also observed the system in the K CONT band (∼2.22 µm) and calculated a depth of 2640 ± 450 ppm.
The new K s -band measurement is consistent with all prior values, except for the significantly larger depth reported by Crossfield et al. (2012).

SUMMARY
In this paper, we analyzed the full TESS light curve of WASP-12, including three sectors' worth of photometry obtained during the spacecraft's ongoing extended mission. We also presented a new K s -band secondary eclipse observation taken with the Palomar/WIRC instrument. Below, we reprise the main results of our work.
1. From our combined fit to all four sectors of TESS data (Section 3.3; Table 2), we significantly improved the precision of the measured phase-curve properties. We obtained a secondary eclipse depth of 466 ± 35 ppm, a peak-to-peak atmospheric brightness amplitude of 524 ± 30 ppm, and a nightside flux that is consistent with zero. A marginal phase-curve offset of 6 • .2 ± 2 • .8 was also detected. The observed ellipsoidal distortion amplitude of 65 ± 14 ppm is consistent with theoretical predictions based on the measured planet-star mass ratio. 2. We fit for 74 individual transit timings (Section 3.1; Table A1) and five mid-occultation times (four from the TESS light curves and one from the Palomar/WIRC observation; Section 3.2). The updated orbital period is 1.091419370 ± 0.000000020 days, with a decay rate of −29.81 ± 0.94 ms yr −1 .
3. From the sector-wide phase-curve fits, we found close agreement in the measured secondary eclipse depths, planetary phase-curve amplitudes, and phase-curve offsets. We placed a tight 2σ upper limit of ≈80 ppm (≈17%) on the dayside atmospheric brightness variability in the TESS bandpass occurring on month-long and year-long timescales.
4. The set of individual TESS-band secondary eclipse depths shows no signs of variability ( Figure 5). We placed a broad 2σ upper limit of 450 ppm on the orbitto-orbit variability of WASP-12b's dayside brightness.

5.
Our new K s -band occultation depth measurement (2810 ± 390 ppm) agrees with all previously published values in the literature, except for the anomalously high value from Crossfield et al. (2012), which we consider to be an outlier.
Funding for the TESS mission is provided by NASA's Science Mission directorate. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. We thank the Palomar Observatory team, particularly Paul Nied, Carolyn Heffner, and Tom Barlow, for enabling the K s -band secondary eclipse observations and facilitating remote operations on the Hale 200 Telescope. I.W. is supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space  Table A1 lists the results from our individual TESS transit light-curve fits. Each transit is labeled by the sector and sequential transit number; Figure 1 shows the location of the corresponding transit events. The second through fifth columns provide the median and 1σ values for the mid-transit time T c , planet-star radius ratio R p /R * , impact parameter b, and scaled orbital semimajor axis a/R * . The fitted transit depths and transit-shape parameter values are mutually consistent to within 2.  Figure A2. Same as Figure A1, but for the Sector 43 transits.  Figure A3. Same as Figure A1, but for the Sector 44 transits.  Figure A4. Same as Figure A1, but for the Sector 45 transits.  Table A1 continued