TASTE V. A new ground-based investigation of orbital decay in the ultra-hot Jupiter WASP-12b ⋆

The discovery of the first transiting hot Jupiters (HJs), giant planets on orbital periods shorter than P ∼ 10 days, was announced more than 20 years ago. As both ground-and space-based follow-up observations are piling up, we are approaching the temporal baseline required to detect secular variations in their orbital parameters. In particular, several recent studies have focused on constraining the efficiency of the tidal decay mechanism to better understand the evolutionary timescales of HJ migration and engulfment. This can be achieved by measuring a monotonic decrease in orbital period d P / d t < 0 due to mechanical energy being dissipated by tidal friction. WASP-12b was the first HJ for which a tidal decay scenario appeared convincing, even though alternative explanations have been hypothesized. Here we present a new analysis based on 28 unpublished high-precision transit light curves gathered over a 12-yr baseline and combined with all the available archival data, and an updated set of stellar parameters from HARPS-N high-resolution spectra, which are consistent with a main-sequence scenario, close to the hydrogen exhaustion in the core. Our values of d P / d t = − 30 . 72 ± 2 . 67 and Q ′ ∗ = (2 . 13 ± 0 . 18) × 10 5 are statistically consistent with previous studies, and indicate that WASP-12 is undergoing fast tidal dissipation. We additionally report the presence of excess scatter in the timing data and discuss its possible origin.


Introduction
The discovery of Jupiter-sized extrasolar planets orbiting their host stars at unexpectedly close distances ("hot Jupiters"; HJ) posed several new questions about planetary formation and planet migration theories.Over the last two decades, astronomers have taken important steps to understand how HJs originated.As of now, the main process suggested to be responsible for such short orbits is migration from their formation site further out (Dawson & Johnson 2018;Fortney et al. 2021).The two most accredited transport mechanisms are disk-driven migration, due to friction with the gaseous disk around the young star (Lin et al. 1996;Nelson et al. 2000), and high-eccentricity migration (HEM), through planet-planet scattering (Rasio et al. 1996) or secular interactions (Fabrycky & Tremaine 2007;Wu & Lithwick 2011).If we consider the current occurrence estimates of one planet every two stars in the Milky Way (Howard ⋆ Photometric data and dditonal tables are available in electronic form at CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.⋆⋆ E-mail: pietro.leonardi.1@studenti.unipd.it⋆⋆⋆ Email: pietro.leonardi@unitn.it⋆⋆⋆⋆ E-mail: valerio.nascimbeni@inaf.itet al. 2012;Dressing & Charbonneau 2013;Batalha et al. 2013;Silburt et al. 2015), it is clear that only a tiny fraction of the total exoplanet population has been detected and studied.The many planetary systems identified so far range in size from tiny rocky planets to massive gas giants, and showing us that diversity is one of the most important aspects of exoplanet demographics.In recent years a new subclass of HJs has emerged, the so-called "ultra-hot Jupiters" (UHJs), having day-side temperatures higher than 2200 K (Parmentier et al. 2018) and orbital periods between ≲ 1 day up to ∼ 2 days.Given their extreme properties, they are valuable test cases for probing atmospheric chemistry (Stangret et al. 2022), understanding planetary mass loss, and constraining planetary formation and evolution models.
WASP-12b, the subject of this study, is one of such UHJs (T eq ≃ 2600 K).Discovered in 2009 by Hebb et al. (2009), as part of the Wide-Angle Search for Planets (WASP) project, WASP-12b is known to be inflated (R p = 1.825 ± 0.091 R jup ; M p = 1.39±0.12M jup ) and orbiting a an F-type star (T eff = 6250 K) with a small projected rotational velocity (v sin i = 2.2 ± 1.5 km/s, Albrecht et al. 2012) in P ≃ 1.09 days.The peculiar architecture of WASP-12b has inspired several research studies on its atmosphere and planet-star interactions.Strong absorption lines of metals in the near-UV spectra were revealed by the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (Fossati et al. 2010;Haswell et al. 2012;Nichols et al. 2015).Optical transit measurements unveiled that the planet is inflated, and infrared phase curves showed that it is filling ∼ 84% of its Roche lobe resulting in high mass-loss rates of metals and heavy elements (Lai et al. 2010;Li et al. 2010;Bell et al. 2019;Antonetti & Goodman 2022).An escaping exosphere was suggested as the cause of these phenomena (Li et al. 2010;Fossati et al. 2010;Haswell et al. 2012).According to hydrodynamic simulations by Debrecht et al. (2018), a circumstellar disk is forming due to the expanding exospheric gas that overfills the Roche lobe and escapes via the Lagrangian L1 and L2 points.Finally, a new atmospheric model also including the formation of H − through dissociation of H 2 O and H 2 has been proposed by Himes & Harrington (2022) to reconcile the previous findings.
Another notable feature of WASP-12b is its changing orbital period.A single planet on an unperturbed Keplerian orbit should transit at regular intervals, i.e. having a perfectly constant orbital period.If the transits are no longer strictly periodic, among the possible causes for this perturbation there are the presence of additional bodies in the system, tidal forces and relativistic effects (e.g., Agol et al. 2005;Holman & Murray 2005;Antoniciello et al. 2021).On close-in planets, tidal effects can cause orbital circularization and a change in orbital period (Rasio et al. 1996;Levrard et al. 2009).Since the early days of exoplanet research (Rasio et al. 1996;Sasselov 2003), it has already been suggested that short-period planets could be unstable due to tidal dissipation.Detection of orbital decay is now within reach because many transiting exoplanets, predominantly HJs, have been monitored for decades, using both space-and ground-based surveys, as part of long-term transit timing observation campaigns.
The first claim of short-term period variations in WASP-12b was published by Maciejewski et al. (2011) and interpreted as due to the dynamical perturbation from an unseen planetary companion (so-called Transit Time Variations; TTV).This detection prompted other observational campaigns to confirm and assess the cause of this signal.Maciejewski et al. (2016) published additional data and discovered a quadratic trend in the O − C diagram of the mid-transit timing compatible with a uniformly decreasing orbital period.Following studies validated this finding (Patra et al. 2017;Collins et al. 2017;Maciejewski et al. 2018;Baluev et al. 2019), which might be equally explained by 1) an orbital shrinkage related to tidal decay or 2) the Rømer effect from an unseen companion or 3) part of a long-term (≃ 14 yrs) oscillation generated by apsidal orbital precession.The latter hypothesis is plausible if the eccentricity is e ≃ 0.002.However, due to the expected fast rate of tidal circularization, it is unclear how to retain such an eccentricity (Weinberg et al. 2017).Yee et al. (2020) presented compelling evidence against apsidal precession and the Rømer effect as major contributors.The authors demonstrated how the O − C values of the occultations have a decreasing pattern compatible with that of the O − C values of the transits, while one should expect the mid transits and occultations to exhibit an opposite trend when the apparent orbital period variation is due to the apsidal precession.Furthermore, their study revealed no significant acceleration of the barycenter of the system along the line of sight using radial velocity observations, disproving the Roemer effect as the primary cause of the observed decreasing O −C trend.This supported the orbital decay scenario, further confirmed by the two latest studies on this system from the TESS space telescope primary and extended missions: Turner et al. (2021) and Wong et al. (2022).The latter research derived a value for the period change rate dP/dt = −29.81± 0.94 ms/yr by including transit and secondary eclipse light curves.
In this work we present 28 unpublished, high-precision, ground-based light curves of WASP-12b.Additionally, we analyze high-resolution optical spectra from the High Accuracy Radial velocity Planet Searcher in the Northern hemisphere (HARPS-N, Cosentino et al. 2012;2014) at the Telescopio Nazionale Galileo (TNG) in La Palma, to derive the stellar parameters of the host star and improve the calculation of the tidal quality factor.In Section 2 of this paper, we present the new photometric and spectroscopic observations and our data reduction procedures.In section 3, we describe the data analysis.The results of the timing analysis aimed at computing the orbital period change rate and constraining the stellar tidal quality factor Q ′ * are reported in Section 4. Finally, in Section 5 we present a summary of our conclusions.

Photometry
The unpublished photometric observations we present were gathered by the TASTE project, a long-term observing program aimed at monitoring transiting planets and creating a library of high-precision light curves to exploit the TTV (Time Transit Variation) technique (Nascimbeni et al. 2011;Granata et al. 2014).The project is mostly based on the Asiago Astrophysical Observatory located at Mount Ekar (elevation: 1366 m) in northern Italy, plus other medium-class telescopes around the globe.The Asiago facilities include the 1.82m Copernico telescope, a Cassegrain telescope equipped with the Asiago Faint Object Spectrograph and Camera (AFOSC) with a 8.8 × 8.8 arcmin 2 field of view, which acquired all but one of our light curves.The AFOSC detector is a back-illuminated 2048 × 2048 E2V pixels CCD.The remaining light curve was obtained at the 67/92 cm Schmidt telescope at the same Observatory.Our observing strategy is based on defocusing the images and applying differential photometry techniques to minimize the impact of systematic errors from instrumental and telluric sources.For all observations, we adopted for consistency the same set of 2 reference stars (TYC 1891-38-1 and TIC 86396443), always imaged over the years in the field of view of the telescope.There are no known variability indicators or higher-than-expected photometric scatter in any of them.A complete description of the TASTE observing strategy and data reduction software is given by Nascimbeni et al. (2011Nascimbeni et al. ( , 2013)).
Our 28 transit light curves were collected over 12 years (2010-2022); a detailed observing log is reported in Table 1, summarizing the dates and other relevant information.All the images were acquired either with a standard Cousins R C filter (λ eff = 672.4nm) or a standard SDSS r ′ (λ eff = 620.4nm), both chosen to maximize the atmospheric extinction and limb darkening effects.The advantage of observing in the R C or r ′ passbands is that the impact of the limb-darkening on the transit profile is reduced with respect to the V passband, thus the timing of the ingress and egress is improved.Moreover, extinction in the Earth's atmosphere is also reduced by observing in those passbands.We applied windowing and 4 × 4 binning to increase the duty cycle of the series as much as possible, while keeping our reference stars within the imaged field.Stars were intentionally defocused to an approximate radius of about 3 arcsec, equivalent to ∼ 12 physical pixels, to avoid saturation and to minimize systematic errors due to pixel inhomogeneities and im-perfect flat-field correction.The exposure time ranged between 2 to 6 s, according to the weather conditions and the amount of defocus applied.For each observation, an out-of-transit part was secured about one hour before ingress and one hour after egress for normalization and detrending purposes.Bias and flat-field frames were collected according to the standard data reduction techniques.We carried out a preliminary selection of our frames, checking for any quality issues (e.g., pixels exceeding the saturation level, strong cosmic-ray hits, etc.).After identifying the problematic frames, we removed the corresponding data points from each light curve.In six of the twenty-eight observations, clouds, high humidity, veils, and other weather conditions prevented a full complete transit from being observed.We did not exclude from our analysis the partial transits when at least the ingress or the egress were well sampled, since most of the timing information lies at those phases where the time derivative of the flux is larger.

Spectroscopy
WASP-12 was observed with HARPS-N at TNG at a spectral resolution of 110 000, in the framework of the GAPS project (Covino et al. 2013).A total of 51 HARPS-N spectra, taken from November 2012 to January 2018, were used in this study.An analysis including the first 15 spectra was published by Bonomo et al. (2017) together with the result of the other 44 stars observed in that study.Additionally, we used eight spectra from the same program, which are still unpublished.The main purposes of these first 23 observations were to refine the planetary parameters and to search for additional companions through highprecision RV monitoring (individual RV errors ranging from 2 to 10 m/s).The remaining 28 spectra were obtained during planetary transits and adjacent off-transit windows on the nights 2017-12-23 and 2018-01-14 as part of the GAPS program to study the planetary atmospheres of hot Jupiters (Guilluy et al. 2022).The individual spectra were reduced with the standard Data Reduction Software version 3.8.A coadded spectrum was then created using all the spectra, achieving a mean S/N of ∼ 220 at around 6000 Å.The HARPS-N spectra described above were exploited in this paper to characterise the host star (see Sect. 3.3).For WASP-12, there are no indications of additional companions, and the updated RV analysis yielded a RV semi-amplitude of 219.90±2.2m/s, corresponding to a planetary mass of 1.39±0.12Mjup, with no significant eccentricity (e <0.02).

Light curve extraction
After the images were corrected for bias and flat field, we performed differential aperture photometry and extracted the light curves by employing the STARSKY photometric pipeline, originally designed for the TASTE project (Nascimbeni et al. 2011;Nascimbeni et al. 2013), and based entirely on an empirical approach.This pipeline performs aperture photometry on the target and a set of reference stars, and combine their fluxes to get highprecision differential light curves, where instrumental or telluric systematic errors are mitigated.The error bars on each data point were calculated with the analytic formulae quoted by Nascimbeni et al. (2011) (Eq. 1 and 2) and through standard statistical propagation.
The aperture radius was set to include most of the flux of the target and the reference stars, and to minimize the off-transit scatter of each light curve.The optimal radius fell in the range of 9-23 pixels.All observations utilized the same set of two reference stars as indicated in Sect.2.1 for consistency.Before the extraction, all the light curves were detrended with a linear function of time at the extraction stage, as part of the STARSKY standard processing (Nascimbeni et al. 2011;Nascimbeni et al. 2013).
In order to represent all our time stamps to a single, accurate time standard, we translated them to the mid-exposure instant and converted them to the BJD-TDB standard, following the prescription by Eastman et al. (2010).It is worth mentioning the light curve ID=23 is the only one in our set showing an unexpected feature, a 5-mmag deep, 6-minute long dip right after the end of ingress (see Fig. A.3).Although the overall RMS of this light curve is higher than average due to a high background level (bright moon at ∼ 20 • ), there is no clear correlation between the feature and any of our diagnostic parameters, including guiding drifts, stellar FWHM, transparency, background, etc.The dip is independent of the choice of reference stars, and is even visible in the absolute light curve.We conclude that the feature is probably genuine, and due to the crossing of an active region during the transit.We will come back to this explanation in Section 4.4.

Light curve modeling
We analyzed the whole set of twenty-eight transit light curves simultaneously using the package PyORBIT1 (Malavolta et al. 2016;2018) developed for modelling planetary transits and radial velocities.In the fitting procedure, the central time of transit T 0 was left free to vary for every transit.As each data set includes a single transit, a uniform prior on each central time of transit T 0 is automatically generated by taking as boundaries the first and last epoch of the corresponding data set.In those cases where only a partial light curve has been observed and the T 0 falls outside the observed window (namely, ID# 10,14,15,16,27), the uniform prior was set manually via visual inspection.
We imposed a Gaussian prior on the orbital period (P) and the stellar radius, based on the values from Bonomo et al. (2017).We adopted a quadratic law to model WASP-12's limb darkening (LD) effect, employing the limb darkening parameterization (q 1 , q 2 ) introduced by Kipping (2013).The limb darkening coefficients are associated with two sets of Gaussian priors, one for each filter, centred on the prediction of the PHOENIX atmospheric models (Husser et al. 2013) with an added conservative uncertainty of 0.10.
We assumed a circular orbit.The effect of non-zero eccentricity on the light curve implies a difference in the ingress and egress times.This difference is 10 −2 × e for a close-in planet (Winn 2010).Since our eccentricity is < 0.02 (See Sect.3.3) the effect on the light curve is negligible.Additionally, for a short period planet orbiting a star as old as WASP-12 (See Table 5), we can consider that the orbit is circularized (Nagasawa et al. 2008, Fig. 4).
For each data set, we included a jitter term to account for possible underestimation of measurement errors and short-term stellar activity white noise.We additionally used a quadratic baseline (three extra free parameters for each light curve) for transits ID# 3, 4, 5, 14, 15, 21, and 22, i.e., for those where the corresponding BIC value indicated a significant improvement in fit as a result of the inclusion of the detrending baseline.Finally, because no contaminants fall inside the photometric aperture adopted, the dilution factor is negligible, and is it not included in the fit.Notes.The columns show: a numerical identifier, the evening date of the observation, the telescope used, the number N of valid scientific frames acquired during the observation, the photometric filter, the average exposure time in seconds, the RMS of the residuals and the approximate phase coverage.
The free parameters in the fit accounted for: three parameters for the transit shape, four for the LD parameters, 28 for the transit times, 28 for the jitter and 15 (3 × 15) for the detrending parameters.This amounts to a total of 78 free parameters.
All the transit models were computed with the popular package batman (Kreidberg 2015), with an exposure time of 5 s and an oversampling factor of 1.The differential evolution algorithm PyDE 2 was used to conduct global parameter optimization.The output parameters served as the starting point for the Bayesian analysis, which was carried out using the emcee (Foreman-Mackey et al. 2013) package, a Markov Chain Monte Carlo algorithm.We ran an auto-correlation analysis on the chains: if the chains were longer than 100 times the estimated auto-correlation time and this estimate varied by less than 1%, the chains were deemed converged.We cautiously set the burn-in value to a number greater than the previously stated convergence point, and we used a precautionary value of 1000 for the thinning factor.We ran the sampler for 75 000 steps, with 680 walkers and a burn-in cut of 25 000 steps.

Stellar parameters
The coadded spectrum was analyzed as in Biazzo et al. (2022) to derive the effective temperature T eff , the surface gravity log g, the microturbulence velocity ξ, the iron abundance [Fe/H], and the rotational velocity v sin i ⋆ of WASP-12.For T eff , log g, ξ, and [Fe/H] we applied a method based on equivalent widths (EWs) of 81 Fe i and 11 Fe ii lines taken from Biazzo et al. (2015) and Biazzo et al. (2022) and the spectral analysis package MOOG (Sneden 1973;version 2017).We then adopted the Castelli & Kurucz (2003) grid of model atmospheres with solarscaled chemical composition and new opacities.T eff and ξ were derived by imposing that the abundance of the Fe i lines is not dependent on the line excitation potentials and the reduced equivalent widths (i.e.EW/λ), respectively, while log g was obtained by imposing the ionization equilibrium condition between the abundances of Fe i and Fe ii lines.
The v sin i ⋆ was measured with the same MOOG code and applying the spectral synthesis technique of three regions around 5400 , 6200 and 6700 Å.We adopted the grid of model atmosphere mentioned above and, after fixing the macroturbulence velocity to the value of 5.5 km/s from the relationship by Doyle et al. (2014), we find an upper limit for the stellar v sin i ⋆ of 1.9 km/s.From the same coadded spectrum, we also derived the abundance of the lithium line at ∼6707.8 Å log A(Li), after measuring the Li EW (39.5 ± 1.5 Å) and considering our stellar parameters previously derived together with the non-LTE corrections by Lind et al. (2009).We therefore obtained 2.55 ± 0.05 dex as NLTE lithium abundance.Lithium equivalent width and elemental abundance values for a star with effective temperature like WASP-12 are intermediate between those of clusters of 2 Gyr, such as NGC752, and ∼4 Gyr, such as M67 (Sestito & Randich 2005;Jeffries et al. 2023).Final atmospheric stellar parameters, together with iron and lithium elemental abundances are listed in Table 4, where uncertainties were computed as in Biazzo et al. (2022).
In order to determine the mass, radius, density and age, we followed the same approach as described in Sect.3.2 of Lacedelli et al. (2021).We used the code isochrones (Morton 2015), with posterior sampling through the code MultiNest (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019)).We provide as input the photometry from the Two Micron All Sky Survey 2MASS, WISE, plus the parallax of the target from the Gaia eDR3 (See Table 5 for the adopted stellar parameters).Analysis has been performed using two evolutionary stellar models, the MESA Isochrones & Stellar Tracks (MIST) (Dotter 2016;Choi et al. 2016;Paxton et al. 2011) and the Dartmouth Stellar Evolution Database (Dotter et al. 2008).
We derived from the mean and standard deviation of all posterior samplings, M ⋆ = 1.325 +0.026 −0.018 M ⊙ , R ⋆ = 1.690 +0.019 −0.018 R ⊙ .The stellar density ρ ⋆ = 0.27±0.010ρ ⊙ was derived directly from the posterior distributions of M ⋆ and R ⋆ .These stellar parameters, together with the age, are listed in Table 4.

Timing Analysis
We searched for the shrinking of the orbital period To constrain the orbital decay effects.Similar to Patra et al. (2020), Yee et al. (2020) and Turner et al. (2021) we implement two models to fit the timing data.The first one is the linear ephemeris model, which assumes that the orbital period is constant: where E, T 0 and P are respectively the transit epoch, the reference mid-transit time (corresponding to E = 0) and the orbital period.The second model is a quadratic one, which assumes that the orbital period is varying uniformly over time: where dP/dE is the change in the orbital period between consecutive transits, which is assumed constant.The best fitting model parameters were found by performing a DE-MCMC analysis (100 000 steps with 1000 burn-in steps) of all the timings in Table 2.We adopted as the new zeroth epoch the transit time closest to the average of the available timings, to minimize the correlation between the ephemeris parameters, as done for instance by  6).
To further test the orbital decay hypothesis, we combine our homogeneous ground-based data set with all the available midtransit times (356) compiled by Bai et al. 2022.The data encompass timings from ground-based facilities as well as TESS results and data by amateur astronomers collected by ETD (Poddaný et al. 2010).All the literature timings and the filters used are listed in the additonal tables at CDS.The resulting parameters of the quadratic fit, performed as in the former case, are:  1.39 ± 0.08 We used these parameters to compute the O − C diagram shown in Fig. 2. We point out that the planetary parameters reported in Table 3 were derived based on Asiago photometry only, and are therefore independent from any literature work.
As already shown by several previous studies (Patra et al. 2017;Yee et al. 2020;Turner et al. 2021) the quadratic model provides a much better fit with respect to the linear ephemeris.This is also confirmed in our goodness-of-fit statistic, by the quadratic model (25 degrees of freedom, d.o.f) chi-square (χ 2 ) 770.87 against the linear model (26 d.o.f) value of 2463.66 and by the Bayesian Information Criterion (BIC; Schwarz 1978).The orbital decay model is favored with a value of ∆(BIC) = 1693 against the linear model.

Tidal Orbital Decay
The measured fast rate of orbital decay for WASP-12b has been attributed to tidal dissipation inside its host star.A short-period planetary system experiences a tidal orbital decay if the host star References.
(1) Two Micron All Sky Survey (2MASS; Cutri et al. 2003;Skrutskie et al. 2006); (2) Gaia eDR3 (Gaia Collaboration 2020) (3) Tycho-2 Catalogue (Høg et al. 2000) (4) Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) rotates slower than the orbital period of the planet.In this scenario, the orbital angular momentum of the planet is transferred to the spin of the star because the tidal bulge lags behind the position of the planet in its orbit.The orbital energy of the planet, deposited in the tidal bulge within the host star, is dissipated by tides, causing the star to spin up (Zahn 1977;Hut 1981;Eggleton et al. 1998).The transfer of angular momentum from the planetary orbit to the stellar spin causes the orbital separation between the planet and the host star to shrink -along with orbital circularization-leading ultimately to the engulfment of the planet after the stellar Roche radius limit is reached (Levrard et al. 2009).Orbital decay in a close-in, short-period exoplanet can be detected via long-monitoring transit timing variations ob-  1 and  2).The reference ephemeris here is the linear part of Eq. 3, while the additional quadratic term (based on dP/dE) is plotted as a black solid line.An excess scatter in the data points with respect to the quadratic model can be seen in the figure.This is addressed in sec.4.4.

Stellar tidal quality factor
The tidal quality factor Q is defined as the ratio between the total energy of the tide and the total energy dissipated in one tidal period (see, e.g.Eq. 2.19 of Zahn 2008).When dealing with tidal dissipation, we pair Q with the Love number k 2 , which takes into account the internal density stratification of the object.Directly measuring the decrease of the orbital period can be used to estimate such a so-called modified stellar tidal quality factor Q ′ * ≡ 3Q * /2k 2 (Matsumura et al. 2010;Hoyer et al. 2016;Penev et al. 2018;Turner et al. 2022).This dimensionless parameter measures the efficiency of the dissipation of the kinetic energy of the tides inside the star.We can deduce its value using the constant-phase-lag (CPL) model of Goldreich & Soter (1966): where M p /M ⋆ is the planet-to-star mass ratio, R ⋆ is the stellar radius, a is the semi-major axis and dP/dt is the period derivative.In the literature, there is a substantial disagreement between the value of Q ′ * determined from WASP-12's studies and the observed or theoretically predicted ranges of values.Our derived value is at the lower limit of the observed ranges as derived from population studies of hot Jupiters (10 5.5 -10 6.5 ; Jackson et al. 2008;Husnoo et al. 2012;Bonomo et al. 2017;Barker 2020) and binary star systems (10 5 -10 7 ; Meibom & Mathieu 2005;Ogilvie & Lin 2007;Lanza et al. 2011;Meibom et al. 2015).Collier Cameron & Jardine (2018) used a form of Bayesian hierarchical inference to numerically estimate the modified stellar quality factor for hot Jupiters.Their analysis found mean values of Q ′ * ∼ 10 7.3 and 10 8.3 , for dynamical and equilibrium tide regimes, respectively.The authors expected a departure from a linear ephemeris for WASP-12b of less than 1 s, on a 20-yr baseline.We can see from Fig. 1 and Fig. 2 that our observations greatly exceed that prediction showing a variation of minutes on a 12-yr baseline.Hamer & Schlaufman (2019) and Barker (2020) argued that the findings of population-wide studies of hot Jupiters should be taken as approximate estimates considering that the value of Q ′ * depends on the physical and mechanical properties (orbital period, stellar interior structure and rotation period, planet and stellar mass) of the individual star-planet systems.In this regard, Rosário et al. (2022) studied WASP-18b, a massive hot Jupiter (∼11 M j ) on a 0.94 days period orbiting a star with a spectral type and mass similar to those of WAPS-12, finding however a lower limit for the tidal quality factor of two orders of magnitude higher than that of WASP-12b.
Furthermore, based on the current understanding of tidal dissipation processes, the physical mechanism driving this efficient dissipation inside WASP-12 is yet to be understood.The tide in the star is frequently divided into two components: an equilibrium tide and a dynamical tide (Zahn 1977).The damping of the equilibrium tide in the convective layers is too weak to explain the observed strong dissipation.Linear and non-linear wave-breaking of the dynamical tide (g-modes) near the radiative core of the star would lead to a sufficiently enhanced dissipation (Q ′ * ∼ 10 5 ) only if the star has entered its subgiant stage of evolution (Weinberg et al. 2017).
The presence and size of the radiative core, which can be derived from the internal structure of the star, is fundamental to understand the tidal decay mechanism.If WASP-12 is still on the main sequence, it would possess a convective envelope and a small convective core.However if it has evolved to become a subgiant it would hold a radiative core, allowing gravity waves to deposit their angular momentum at the radiative-convective boundaries by breaking, resulting in a fast tidal decay (Weinberg et al. 2017;Weinberg et al. 2024).Our precise observational constraints on the stellar parameters, such as mass, effective temperature, age and metallicity, are however consistent with the main-sequence structure if compared to the theoretical models of Weinberg et al. (2017) and Bailey & Goodman (2019).In order to further investigate the evolutionary phase of WASP-12 we used the MESA Isochrone & Stellar Tracks (MIST) stellar evolution models, to extract the theoretical isochrone equivalent to our derived value of [Fe/H] and stellar age (see Table. 4).We placed the star in the log(T eff )-log(L/L ⊙ ) diagram (H-R diagram), plotted in Fig B .1, where the isochrone is displayed as a solid black line.We can notice that the star is located in an ambiguous position, between the main sequence and the turn-off, but far away from the subgiant branch.This could tell us, as suggested by Efroimsky & Makarov (2022), that the star has entered the post-turn-off stage where hydrogen fusion in the core has ceased,ultimately leading to the core contraction and the burning of hydrogen in a shell around the core.
Recent studies proposed alternative hypotheses in which the observed tidal orbital decay can be explained by obliquity tides (Millholland & Laughlin 2018), planetary eccentricity tides (Efroimsky & Makarov 2022), or a spin-orbit coupling mediated by a non-axisymmetric gravitational quadrupole moment inside the planet (Lanza 2020).Although significant improvements have been made in our observations of the system, there is still a need for more precise modeling to address the numerous unresolved questions and test both stellar evolution and tidal theories.A promising approach to determine whether WASP-12 is a main-sequence or a sub-giant star by measuring the age of the star, the mixing length parameter α, the size of the convective core or the presence of a core rotation, is asteroseismology, as suggested by previous studies (Deheuvels et al. 2016;Weinberg et al. 2017;Bailey & Goodman 2019;Fellay et al. 2023).

The search for orbital decay in other systems
Among the population of hot Jupiters, WASP-12b is the first case where orbital decay has been detected.Several studies over the years investigated the decaying orbit of hot Jupiters, which possess orbital features similar to the target of our study (Patra et al. 2020;Ivshina & Winn 2022;Hagey et al. 2022).Harre et al. (2023) used CHEOPS and TESS high-precision photometric data to investigate the orbital decay of WASP-4 b, KELT-9 b and KELT-16 b.They found, in agreement with Turner et al. (2022) amd Bouma et al. (2019), a changing orbit only for WASP-4b, but even if the orbital decay model is preferred, apsidal precession cannot be ruled out with the available data.Additional transit and occultations are needed to clearly identify the cause of the period variation.Vissapragada et al. (2022) presented strong evidence for the tidal decay of Kepler-1658 b, a close-in giant planet orbiting a sub-giant star, using a thirteen-year observational baseline.Using RVs, the authors ruled out apsidal precession on theoretical grounds and line-of-sight acceleration.The computed value of the modified tidal quality factors (Q ′ * = 2.50 +0.85 −0.62 × 10 4 ) fully agrees with what should result from the dissipation of inertial waves in the convective zone.The recently discovered hot-Jupiter TOI-2109b (Wong et al. 2021) presents orbital and stellar characteristics that could lead to a stronger tidal dissipation than WASP-12b.Still new observations with a long observational baseline are needed to detect a potential orbital decay.
According to expected tidal decay rates and host star properties, many other planets should exhibit significant rates of orbital decay, e.g., WASP-18b (Collier Cameron & Jardine 2018), WASP-19b (Rosário et al. 2022), WASP-43b (Patra et al. 2020), WASP-103b (Barros et al. 2022), WASP-161b (Yang & Chary 2022), KELT-16b (Harre et al. 2023), HATS-18b (Southworth et al. 2022).So far, not enough evidence has been found from the available data to confirm a decaying orbit for any of these systems.Several effects can produce an O − C variation over a long baseline, thus to have a successful validation of a planet tidal orbital inspiral and discard any other different possible cause it is necessary a decade-long observational baseline.Many of these targets will be re-observed by TESS, CHEOPS, PLATO (Nascimbeni et al. 2022) and Ariel (Borsato et al. 2022) in the near future and could be used in a future collective study.However, as advocated by Hamer & Schlaufman (2019), we additionally suggest that the focus of orbital decay studies should be on planets orbiting host stars at the end of their main-sequence lifetime.

Excess scatter in the timing data
An evident feature of our data is the presence of an excess scatter in the O − C plot (Fig. 1): the reduced χ 2 of the fit on the Asiago data set is χ 2 r = 5.52 (χ 2 ∼ 138 for 25 d.o.f), implying that the error bars are about ∼ 2.3 times too small with respect to a well-behaved Gaussian distribution having the same variance.This cannot be attributed to issues in the absolute calibration of our time stamps, or to systematic errors left uncorrected by our photometric pipeline, as other TASTE studies carried out with the same setup and data analysis technique demonstrated a timing accuracy at the 1 s level or better (Brown-Sevilla et al. 2021).Such an excess is not unique to our unpublished data set: the best fit to the overall O − C diagram (Fig. 2) yields χ 2 r ∼ 2, and other independent data sets reach χ 2 r significantly larger than one (e.g., χ 2 r ∼ 2.1 for Baluev et al. 2019 and χ 2 r ∼ 8 for ETD).We investigated the impact of correlated noise on our results by analyzing the residuals of our light curve, i.e. after the best-fit transit model has been subtracted.This was done by calculating the β factor after binning the time series at different time scales (Gillon et al. 2006;Winn et al. 2007a).When correlated noise is present, the binned rms σ is larger by a factor of β with respect to the ideal white-noise case, where the rms is supposed to decrease by the square root of the number of points per bin ( √ N).At the time scale corresponding to the ingress/egress duration (∼23 min for WASP-12b), i.e., the most relevant time scale for the transit time estimation, β usually lies in the 1.5-2 range for wellbehaved ground-based photometry (Winn et al. 2007b).Among our 28 light curves, only four of them (#2, 7, 20, 28) show an unusually large amount of red noise (β ≫ 2), but #20 is the only one also being a significant outlier in the O − C diagram.
As a further check, we binned each complete light curve at the time scale of the ingress duration (0.0159 days), then redid the global fitting process on complete transits in the exact same way as explained in Section 3.2.The resulting χ 2 r is 3.2 (χ 2 = 64.8 with 20 d. o. f.), that is, still considerably larger than one.Even though the error bars increased, the sign of the most significant outliers is unchanged.We also emphasize that our best-fit dP/dt in the binned case is −29 ± 4 msec/year, fully in agreement with the previous value of −30.7 ± 2.7 msec/year obtained without binning the light curves (Eq.3).
These exercises demonstrate that, while a few light curves show a higher-than-usual amount of correlated noise, there is no obvious link with the excess scatter we found in the measured transit times.We also emphasize that the presence of partial transits in our data set cannot explain said excess: the reduced χ 2 for just our 22 complete transits is χ 2 r = 5.78 (χ 2 = 110 for 19 d. o. f.), almost identical to the χ 2 r = 5.52 value for the full distribution (χ 2 = 137 for 25 d.o. f.).
It is well known in the literature that comparing timing data extracted from ground-based transit light curves can be prone to systematic errors due to a combination of telluric/instrumental red noise and different reduction techniques (Baluev et al. 2019).On the other hand, excess scatter is sometimes noticed also when dealing with space-based data (as in Turner et al. 2022;Harre et al. 2023 in the case of WASP-4), and it is usually attributed to the presence of stellar activity.Little is known about the stellar activity properties of WASP-12b, since its anomalously low R HK value could be due to the presence of extra absorption from escaped circumstellar gas (Fossati et al. 2010) and no coherent rotational modulation can be detected from its light curve.Yet, significant bumps are visible in the residuals of the TESS data (Wong et al. 2022), and at least one of our light curves (#23) shows an in-transit anomaly which is difficult to explain with systematic errors (Section 3), suggesting that in some cases inhomogeneities may arise on the stellar photosphere of WASP-12 but do not result in a detectable rotational modulation because of the high inclination of the stellar axis with respect to the observer, also suggest by the unusually low v sin i.If the pole-on hypothesis is true, the timing scatter we see could be the only observable hint so far that WASP-12 is actually a significantly active star.
Intriguingly, even though running a GLS periodogram on our O − C residuals does not result in any significant frequency, the excess scatter appears to be a time-dependent phenomenon.Transits gathered during some observing seasons (such as 2011-2012 and 2019-2020) line up within seconds from the average quadratic ephemeris and return χ 2 r ∼ 1 even on the Asiago subset, while others (such as 2020-2021) appears much more scattered on different and independent data sets.
Forthcoming additional photometry from TESS (scheduled on two sectors -71 and 72-in Fall 2023) will possibly shed some more light on this behaviour.Meanwhile, we emphasize that, despite the scatter excess, our homogeneous data set, published here for the first time, represents a fully independent confirmation of the tidal decay process in WASP-12b, at an estimated rate perfectly consistent with those derived by previous studies.

Summary and Conclusions
In this paper we inspected the tidal orbital decay of the Ultrahot Jupiter WASP-12b, using 28 unpublished, high-precision, ground-based light curves collected over twelve years, from 2010 to 2022, as part of the TASTE program.We applied a DE-MCMC approach, for a linear (constant period) and quadratic model (orbital decay), to the timing of our homogeneously collected observations confirming the decaying orbit of WASP-12b.
Our new planetary and transit parameters are in good agreement with previous studies.However, combining all the timing data from literature we found some discrepancies in the O − C plot.The timing data shows an excess of scatter that cannot be attributed to errors in the photometric pipeline or the calibration of the time stamps.This raises several questions which will require further study to resolve.
We used HARPS-N observations to investigate the evolutionary stage of the host star, which is fundamental to explain the high rate of tidal dissipation.The applied spectroscopic observations allowed us to derive precise stellar parameters of the host star, however, we did not arrive at a firm conclusion regarding the evolutionary stage.The stellar parameters enable us to place the star in the main sequence close to the turn-off point and away from the subgiant branch.However, the fast tidal dissipation requires a subgiant scenario.
: TASTE V. A new ground-based investigation of orbital decay in the ultra-hot Jupiter WASP-12b Dragomir et al. (2011) andMallonn et al. (2019).The results of the quadratic model are:T 0 = 2457280.093510± 0.000067 (BJD-TDB), corresponding O − C ("Observed minus Calculated") diagram is shown in Fig.1as a function of the epoch E. The derivative of the orbital period dP/dE is derived from the quadratic coefficient of the best fitting parabola in Eq. 2, and then translated into dP/dt.Its negative sign indicates a decrease in the orbital period with a rate of 0.03071 seconds per year.The corresponding orbital decay timescale would be τ = −(dP/dt)/P = 3.07 ± 0.26 Myr, a value consistent with the results ofTurner et al. (2021) andWong et al. (2022) (See Table

Fig. 1 .
Fig. 1.Observed minus calculated (O − C) diagram of WASP-12b for the unpublished transits collected at Asiago (Tables1 and 2).The reference ephemeris here is the linear part of Eq. 3, while the additional quadratic term (based on dP/dE) is plotted as a black solid line.An excess scatter in the data points with respect to the quadratic model can be seen in the figure.This is addressed in sec.4.4.

Fig. 2 .
Fig. 2. Observed minus calculated (O − C) diagram of WASP-12b for all the available transit timing data.The reference ephemeris here is the linear part of Eq. 4, while the additional quadratic term (based on dP/dE) is plotted as a black solid line.The Asiago data are plotted as black and red points, while the remaining points come from the literature: Bai et al. (2022) (red points, including TESS), Wong et al. (2022) (light green), Baluev et al. (2019) (light blue), other literature works (available at CDS, pink).Amateur timings from ETD are plotted in grey.
By taking our derived values of M ⋆ , dP/dt, R ⋆ , a, see Sect.3.3 and Sect.3.4, and the planetary mass of Chakrabarty & Sengupta (2019), we obtained a modified tidal quality factor of Q ′ * = (2.13 ± 0.18) × 10 5 .(6) The planet mass is assumed to be constant, which we know from Lai et al. (2010) not to be the case.This value is slightly higher, but still consistent, compared to what is derived by Turner et al. (2021) and Wong et al. (2022).

Table 1 .
Observing log of our unpublished transit light curves of WASP-12b.

Table 2
shows all the 28 transit central times obtained from the PyORBIT fit.The resulting best-fit parameters and the priors are shown in Table3.The light curves and the corresponding PyORBIT best-fitting models are plotted in Figures A.1-A.4.

Table 2 .
New best-fit transit times of WASP-12b.
(Eastman et al. 2010)Ds match those assigned on Table1.The transit times are given in the BJD-TDB standard(Eastman et al. 2010); the third column reports the associated 1-σ error.

Table 3 .
WASP-12b orbital parameters derived from PyORBIT Markov Chain Monte-Carlo analysis.The listed best-fit values and uncertainties are the medians and 15.865th-84.135thpercentiles of the posterior distributions, respectively.

Table 4 .
Stellar parameters of WASP-12 derived from the analysis of the HARPS-N spectra and from stellar evolutionary tracks.