Looking For Timing Variations in the Transits of 16 Exoplanets

We update the ephemerides of 16 transiting exoplanets using our ground-based observations, new TESS data, and previously published observations including those of amateur astronomers. All these light curves were modeled by making use of a set of quantitative criteria with the EXOFAST code to obtain mid-transit times. We searched for statistically significant secular and/or periodic trends in the mid-transit times. We found that the timing data are well modeled by a linear ephemeris for all systems except for XO-2 b, for which we detect an orbital decay with the rate of -12.95 $\pm$ 1.85 ms/yr that can be confirmed with future observations. We also detect a hint of potential periodic variations in the TTV data of HAT-P-13 b which also requires confirmation with further precise observations.


INTRODUCTION
Since the observations of the first transits in an exoplanet system (Charbonneau et al. 2000), several questions have ‹ E-mail: yalcinkayas@ankara.edu.trarisen regarding their formation, evolution, atmospheric composition, and orbital dynamics.These questions can be further investigated through different observational techniques.For example, radial velocity measurements during transits can be employed to determine the obliquity of a planet's orbit, which in turn can provide important infor-mation for improving theoretical models related to orbital evolution (Mancini et al. 2022).Occultation observations can provide information about the planet's energy budget (Arcangeli et al. 2021), while transmission spectroscopy can reveal its atmospheric composition (Maguire et al. 2022).However, these observations require high levels of precision, which can only be achieved by making use of large groundbased or space-borne telescopes.As the observation time for these instruments is in great demand, accurate predictions of transit and occultation times are crucial.Even small uncertainties in transit times can accumulate over time and require updates to the exoplanet's orbital period and reference mid-transit times (Mallonn et al. 2019).
Tidal interactions can cause the orbit of the planet to shrink (Maciejewski et al. 2016b).The period decrease per year may be much smaller than the uncertainty of the midtransit times, making it difficult to observe.The amplitude of this effect increases over time and may be detected with additional transit (or occultation) observations over a long time range.In addition, transit timing analysis can be used to detect unseen additional bodies in a system that could not be seen with radial velocity (RV) observations due to short phase coverage, stellar activity (Trifonov et al. 2021) or if the host star is too faint for precise RV observations (e.g.Gillon et al. 2017).For eccentric systems, the secular motion of the periastron (i.e.apsidal motion, Giménez & Bastero 1995) is observable with the help of occultation observations (Patra et al. 2017) and could give insights about tidal effects.To identify these effects using the transit timing variation (TTV) technique, it is essential to have transit timing measurements that cover longer time spans and are well sampled.
We selected potential periodic TTV targets depending on the known third bodies in the system or depending on their radial velocity residuals.
In the potential TTV group, there are also orbital decay candidates selected based on their stellar and planetary radii, masses and orbital separation and ages to work on systems with maximum tidal interaction potential.There are unitless metrics that are used to select our candidates.Please see Baştürk et al. (2022).We observed 38 transit of 16 exoplanets (GJ 1214 b, HAT-P-1 b, HAT-P-10 b, HAT-P-13 b, HAT-P-16 b, HAT-P-22 b, HAT-P-30 b, HAT-P-53 b, KELT-3 b, QATAR-2 b, WASP-8 b, WASP-44 b, WASP-50 b, WASP-77A b, WASP-93 b, XO-2 b) that we selected for their potential to display TTVs and/or large shifts in their observed transit timings.
The transit data that we used for timing calculation were obtained from ground-based telescopes and the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al. 2015), and compiled from published observations and open databases1,2 .We performed homogeneous transit timing analyses of these systems and updated their ephemeris information.
This paper is organized as follows.In section 2, we describe the telescopes and the detectors we used for transit observations, data reduction, and photometry procedure as well as light curve selection criteria.TTV analyses and our results are presented in section 3. We discuss our findings in section 4.

Observations
Photometric transit observations were carried out with the T100, T80, ATA50, UT50 and (numbers in the name of telescopes come from primary mirror diameters in cm) CAHA 1.23 m telescopes.Detailed information about the telescopes and their detectors can be found in Baştürk et al. (2022).We also observed a multi-color transit of HAT-P-1 b with the Bonn University Simultaneous Camera (BUSCA) on the CAHA 2.2 m telescope at the Observatory of Calar Alto (Spain).We made use of the well-established defocusing technique (Southworth et al. 2009) in order to increase photometric precision.Exposure times were set to acquire at least " 50 frames per transit.The defocusing amount was determined to keep the detector response within its linearity limits while exposing it for larger durations to increase the Signal-to-Noise ratio (SNR) by reading out from a larger area.In general, we selected the photometric filter that gives the maximum SNR.A detailed log of photometric observations is provided in Table 2.

Data Reduction
Data reduction (dark, bias and flat correction) and ensemble aperture photometry were performed using the AstroIm-ageJ (hereafter AIJ) (Collins et al. 2017) software package.To increase the precision in photometry, we selected every star similar in brightness to the target in the field as a comparison, AIJ allows the user to visually inspect the relative flux of the target for a combination of different comparison stars.After finding suitable comparison stars, we experimented with different aperture sizes for both the stars and the sky background, AIJ also allows users to visually inspect relative flux change due to different aperture sizes.When selecting comparison stars and aperture sizes, our goal was to minimize the red noise, especially in contact times where the flux change is abrupt.Red noise during ingress and egress can change the mid-transit times dramatically, but may not affect the error bar of individual data points which results in an underestimation of the mid-transit time uncertainty (Pont et al. 2006;Gillon et al. 2006).This could lead to a higher reduced chi-square (χ 2 ν ) for linear ephemeris, which could (incorrectly) be attributed to TTV.In order to avoid that, we detrended relative fluxes by using time-dependent variables such as airmass and target position on the CCD in an interactive manner using AIJ.
For TESS observations, we downloaded the two-minute light curves from Mikulski Archive for Space Telescopes3 (MAST) that are processed by Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016).SPOC generates presearch data conditioning (PDC) light curves and data validation time series (DVT) light curves using simple aperture photometry (SAP).The PDC SAP fluxes are the corrected version of the SAP fluxes from instrumental systematics, outliers and flux contamination from nearby stars.The DVT light curves are created by applying a running median filter to the PDC light curves to remove any long-term systematics and search for transits.We used only the DVT light curves because any signal other than transits will deteriorate the transit profiles, which in turn will increase the uncertainty in the measurement of the midtransit times.For the case of XO-6 b, Ridden-Harper et al.
(2020) has shown that the DVT light curves have least scatter, nevertheless, the transit timings from DVT and PDC light curves are practicaly identical.We have TESS light curves from the SPOC pipeline for all the planets in our sample except for HAT-P-53, which was observed by TESS during Sector 17, but light curves were not produced.Therefore we downloaded the Full Frame Images (FFI) that has 30 minutes cadence from TESScut4 and performed aperture photometry with the lightkurve package (Lightkurve Collaboration et al. 2018) and then detrended the light curve using keplerspline-v25 while ignoring the transit profiles.Final light curves were not suitable for individual modeling due to insufficient sampling so we time-folded the data using a period from our preliminary analysis.We assumed the period of HAT-P-53 to be constant during Sector 17 but this enabled us to measure only a single mid-transit time from TESS observations.We included every TESS light curves until the end of the extended mission 2 (Sector 69) to our analysis.

Light Curve Selection Criteria
The main goal of this work is to search for TTVs in the planetary systems listed in Table 1.This requires precise and accurate mid-transit times measured from high-quality light curves.For this reason, we used the light curve selection criteria given in Baştürk et al. (2022) to select suitable light curves.First, we compiled available transit light curves from literature, open databases of amateur astronomers (Exoplanet Transit Database6 , hereafter ETD and Amateur Exoplanet Archive7 , hereafter AXA) along with our own observations and observations from space telescopes (TESS and Kepler Space Telescope's K2 mission Howell et al. 2014).We did not include light curves that have large gaps inside transit profiles or high-amplitude signatures of correlated noise, especially in the ingress or egress segments.Then we modeled light curves with the exofast (Eastman et al. 2013) (see Section 3 for details) and then calculated photometric noise rate (PNR) (Fulton et al. 2011a) from residuals which indicates white noise.We removed the light curves that have PNR values higher than the transit depth.We binned the residuals between the ingress/egress duration ˘5 minutes with 1-minute steps and calculated the well-known β values as defined in Winn et al. (2008) as a red noise indicator.We removed the light curves with the median β-value larger than 2.5.We also removed the light curves if the transit depth is a 5σ outlier for the given planet.When we visually inspect the removed light curves, we find that this criterion is very useful to detect problematic light curves.Qatar 2 b is an exception because it has K2 light curves with incomparably higher precision than other datasets, affecting the σ-value dramatically.Thus we did not include depth values from K2 for Qatar 2 in the calculations of its 5σ level.

Light Curve Modelling and Measurements of
Mid-Transit Times We followed the same method given in Baştürk et al. (2022) to model the light curves and measure the mid-transit times.Briefly, we used exofast-v1 to model the light curves after converting the observation time to Dynamical Barycentric Julian days (BJD-TDB) and detrending the light curves for the airmass effect that need it with the AIJ.We used our scripts to convert the timings to BJD-TDB and calculated 2 Eliminated because its depth is out of 5σ of the average.
3 Eliminated because its PNR value is larger than its depth.
4 Eliminated because its β-factor is larger than 2.5.
the airmass values by using relevant modules and functions of the astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018) ) library.The centers and widths of the priors were automatically selected from the NASA Exoplanet Archive8 for the atmospheric parameters of the host stars as a Gaussian priors and the orbital periods of the planets as constant values while uniform priors of the limb darkening coefficients were automatically retrieved from Claret & Bloemen (2011) based on the atmospheric parameters of the host stars and the observed passbands.For the passbands that are not available, we choose the passband that has the closest transmission curve (e.g.we choose I for the TESS passband and CoRoT for the clear observations).After selecting the light curves as described in 2.3 and measuring the mid-transit times from the individual transit models using a built-in IDL routine AMOEBA that uses downhill simplex method (Nelder & Mead 1965) to minimise χ 2 , we constructed the TTV diagrams and fitted a linear ephemeris using emcee (Foreman-Mackey et al. 2013) following the recipe given in Baştürk et al. (2022).We discarded the two-tailed 3σ outliers from the linear ephemeris not to bias our final results.These light curves with correlated noise, especially during ingress or egress, may sur-vive the β > 2.5 criteria and result in an inaccurate midtransit time with underestimated error bars.We visually inspected the light curves eliminated based on this criterion and we noticed that this criterion is especially useful for the light curves that come from spectroscopic observations (i.e.white light curve, formed by integrating an observed spectrum over the entire wavelenght coverage) because these light curves usually have very high precision (hence low midtransit time error) but inaccurate mid-transit times due to heavy detrending.We also noticed that this criterion enables us to detect light curves with incorrectly reported time references.We did not apply two-tailed 3σ outlier criteria for the XO-2 system because we detect a statistically significant orbital period decrease.For the Qatar 2 system, we calculated the 3σ value without including K2 light curves but a linear ephemeris was fitted to all data points including K2.

Ephemeris Corrections
For all systems, we fitted independent linear and quadratic ephemeris using the emcee package.We followed the same procedure as described in Baştürk et al. (2022) for selecting random walkers, burn-in period and Markov Chain Monte Carlo (MCMC) steps for convergence.The median values of posterior probability distributions (PPD) of linear elements; slope and y-intercept were added to the reference period (P orb ) and mid-transit time (Tc) respectively.The updated linear ephemerides are listed in Table 3 with their uncertainties calculated from PPD.

Transit Timing Analyses
In order to detect potential secular changes in the orbital periods, we fitted quadratic functions to the TTVs of all planets using the method described in section 3.2.We compared the quadratic ephemeris with the linear to detect any significant secular change for the planets in our sample.
In Table 4, we report the Akaike Information Criterion differences (∆AIC) and Bayesian Information Criterion differences (∆BIC) values between linear and quadratic ephemeris and the rate of the secular period change calculated from the coefficient of the second-degree term of the quadratic ephemeris.We only consider ∆BIC > 10 as suggested strong evidence by Raftery (1995) to favour quadratic over the linear ephemerides.
After correcting the ephemerides (displayed in Figure 1) using the linear coefficients, we performed a frequency analysis to search for potential periodic variations that can be caused by orbital perturbers or the apsidal motion of the planets.We used the astropy's Lomb-Scargle function (VanderPlas 2018) to find possible frequencies and their False Alarm Probabilities (FAP).

GJ 1214 System
GJ 1214 b is a sub-Neptune planet (Mp = 6.55 M C , Rp = 2.678 R C ) that orbits an M dwarf star.It has a very high Transmission Spectroscopy Metric (TSM) (Kempton et al. 2018), making it one of the most favorable sub-Neptune planets for atmospheric studies (Charbonneau et al. 2009).Additional bodies in the system have been searched for using the radial velocity (RV) method with 165 RV points spanning 10 years (Cloutier et al. 2021), as well as with the transit method using a continuous observing run for " 21 days from the Spitzer Space Telescope (Gillon et al. 2014).Follow-up transit observations have been performed multiple times to investigate TTVs or the atmospheric properties of the planet (Kundurthy et al. 2011;de Mooij et al. 2012;Harpsøe et al. 2013;Narita et al. 2013;Cáceres et al. 2014;Nascimbeni et al. 2015;Parviainen et al. 2015;Rackham et al. 2017;Angerhausen et al. 2017;Mallonn et al. 2018;Orell-Miquel et al. 2022;Spake et al. 2022;Lampón et al. 2023;Gao et al. 2023).We selected the planet for its potential to display TTVs as well as updating its ephemeris for future observation plans especially to understand its atmospheric properties.
We analyzed a total of 48 light curves, including 6 from the Exoplanet Transit Database (ETD), 37 from the literature, and 5 from our observations.However, 5 of the light curves did not meet our selection criteria and were eliminated (as explained in Section 2.3).The total data span 10 years of observations, but there was a 5-year gap in the TTV diagram.After analyzing the TTV diagram, we did not detect any significant period change in the GJ 1214 system.

HAT-P-1 System
HAT-P-1 b is a warm Jupiter with low-density (Mp = 0.53 MJ, Rp = 1.36 RJ) orbiting a G0 V type star discovered by Bakos et al. (2007).The host star is part of a wide binary system with a companion (HAT-P-1A) of similar effective temperature, making it an excellent comparison star for atmospheric observations in high angular resolution.The planet has a relatively high TSM, which makes it a favorable object for atmospheric studies using groundbased and space-borne telescopes (e.g.Montalto et al. 2015;Wakeford et al. 2013).Bakos et al. (2007) suggested a small eccentricity that could be attributed to perturbations by an outer companion, which could be discovered by RV or TTV observations.With follow-up RV observations, Ment et al. (2018) rejected the eccentric orbit and Johnson et al. (2008) found that the spin of the orbit of HAT-P-1 b is aligned with the stellar rotation axis.Winn et al. (2007) and Johnson et al. (2008) found no significant TTVs in the system.
Here we analyzed 26 transit light curves, 3 of which were eliminated, to update the ephemeris of HAT-P-1 b.We found no statistically significant periodic or parabolic change in the period analysis.We updated the ephemeris of transit which can be very useful for future atmospheric observations.

HAT-P-10/WASP-11 System
HAT-P-10 b is a low-mass, hot Jupiter that was independently discovered by Bakos et al. (2009a) and West et al. (2009).Follow-up radial velocity (RV) observations by Knutson et al. (2014) revealed a linear trend that suggested the presence of a stellar-mass companion.Adaptive-optic (AO) observations by Ngo et al. (2015) revealed the existence of a 0.36 Md companion at a distance of 42 AU ("0.235 2 ), which can explain the RV trend.Ngo et al. (2016) showed that the companion can not cause Kozai-Lidov migration of the planet, and the eccentricity of the planet is consistent with zero as expected.The Rossiter-Mclaughlin (RM) observations by Mancini et al. (2015) indicate that the system is aligned, and this alignment has a primordial origin rather than being due to tidal interactions, owing to the relatively long distance between the star and the planet.Therefore, we do not expect to observe orbital decay in this system.Wang et al. (2014) investigated the TTVs to detect any outer companion with the light-time effect (LiTE), but found the orbital period of HAT-P-10 b to be constant.We included this system in our study for the same reasons and studied its TTV diagram with more data spanning a longer baseline.
We conducted an analysis of 29 transit light curves, consisting of 16 from ETD, 7 from literature, 4 from TESS, and 2 from our own observations.However, we excluded 4 of them and ultimately derived a TTV diagram from 25 midtransit times that were evenly distributed across a span of 13 years.Our analysis of the TTV diagram did not reveal any significant periodic changes or deviations from a constant period.

HAT-P-13 System
HAT-P-13 b is a warm Jupiter discovered by Bakos et al. (2009b), revolving around a Solar-like, metal rich (T eff " 5653 K, rFe{Hs " 0.41) and slightly evolved star.The system consists of at least another planet, HAT-P-13 c, highly eccentric (e " 0.691), long period (Pc " 446.27 days), massive (Mp sin i " 15.2 MJ) outer companion discovered with RV observations.The presence of another outer companion is suggested by the linear trend of RV residuals, as noted by Winn et al. (2010a) and Knutson et al. (2014).HAT-P-13 was suggested to have a cooler companion (T eff " 3900 K; Piskorz et al. 2015) blending its lines in its infrared spectrum.However, the AO observations do not reveal a companion (Ngo et al. 2015), within the limits of the study given in their Figure 4 making the system worthwhile for TTV investigations.We used 28 light curves to construct the TTV diagram of HAT-P-13 b, after eliminating seven of them.We found that the period of HAT-P-13 b deviates from a constant period.The frequency analysis revealed a peak at 479.52 days with a FAP of 0.0007 and the full TTV amplitude is " 321 seconds (see Fig 2).Assuming that planet c is the perturber, and the system is coplanar, the TTV amplitude caused by a planet c should be approximately 40 seconds, as previously calculated by Bakos et al. (2009b).The RM observations by Winn et al. (2010a) revealed that the orbit of HAT-P-13 b is aligned, which supports the coplanar scenario.However, the transit of HAT-P-13 c has not been observed in longterm observations by Fulton et al. (2011b) and Szabó et al. (2010).Therefore, we conducted a preliminary Newtonian orbital analysis to fit the RVs and TTVs and found that the inclination of the putative planet c must be " 2 ˝, and its mass should be " 0.4 Md to cause a 321-second TTV.If this is the case, we would expect the impact parameter, b, to vary over time, making HAT-P-13 b's orbit misaligned.Some of the transit light curves of HAT-P-13 b exhibit modulations that can be attributed to star spots.This makes it challenging to accurately measure the mid-transit times, which could introduce fallacious TTVs.Additional observations, including upcoming TESS data and new ground-based observations, are needed to determine the true ephemeris of HAT-P-13 b.We suggest that these light curves require special treatment, such as Gaussian Process (Yalçınkaya et al. 2021) or spot modeling (Mancini et al. 2017) for better accuracy.

HAT-P-16 System
HAT-P-16 b is a dense, (Mp " 4.193 MJ, Rp " 1.289 RJ) hot Jupiter (P = 2.775960 days) orbiting an F8 dwarf, discovered by Buchhave et al. (2010).The planet was found to have a small but statistically significant eccentricity based on its RV observations (Buchhave et al. 2010;Bonomo et al. 2017) and its projected spin-orbit angle suggests that it is aligned HK index as -4.863, which indicates low magnetic activity (e.g.Noyes et al. 1984), andCiceri et al. (2013) found no starspot induced anomalies in the transit light curves, which is also indicative of low magnetic activity.Hence, the system should not be very young.Using the log pR 1 HK q´stellar rotation period (Prot) calibration from Suárez Mascareño et al. (2015), we found Prot " 22.5 days and vrot " 2.8 km s ´1, meaning that the stellar inclination (I‹) is consistent with 90 ˝within uncertainties (Vsini = 3.5˘0.5,Buchhave et al. 2010).This result, combined with the RM values, suggests that the orbit of HAT-P-16 b is well-aligned.Winn et al. (2010b) speculates that hot Jupiter systems may have primordial misaligned orbits, but the tidal dissipation in the convective zones of their host stars can lead to spin-orbit alignment.Considering the relatively high effective temperature of HAT-P-16, the star should have a thin convective zone.The T eff cut-off at which the star will have a negligible convective mass was determined at 6250 K by Pinsonneault et al. (2001), while HAT-P-16's T eff is 6158 K.Then, it should take a few Gyr for HAT-P-16 to diminish the primordial obliquity.On the other hand, using Eq. ( 2) from Adams & Laughlin (2006b) and the limits for Qp between 10 5 and 10 6 as given by them, the tidal circularization timescale is only 400 Myr even if the Qp is taken to be 10 6 .Assuming the system is at least a few Gyrs old based on its magnetic activity, the non-zero eccentricity may have been caused by an outer companion (Adams & Laughlin 2006a), which may have led to Kozai-Lidov oscillations.Sada & Ramón-Fox (2016) searched for TTVs, but did not detect a definitive signal because there were too few observations.Sun et al. ( 2023) detected orbital decay with ∆BIC = 167 and apsidal motion with ∆BIC = 317.We included light curves from several works (Aladag et al. 2021;Buchhave et al. 2010;Sada & Ramón-Fox 2016;Ciceri et al. 2013;Pearson et al. 2014), adding up to a total of 62 light curves, nine of which were eliminated, hence we were able to form a TTV diagram covering the widest time range available for analysis.The recently published TESS sector ruled out the orbital decay suggested by Sun et al. (2023).We also did not detect any significant cyclic TTV, as suggested by Sun et al. (2023), that can be caused by the apsidal motion of the orbit.Although we did not detect any significant cyclic or parabolic changes, we updated the ephemeris for future observations.

HAT-P-22 System
HAT-P-22 b is a relatively dense (Mp " 2.147 MJ, Rp " 1.080 RJ), slightly eccentric (e " 0.0064 `0.0080 ´0.0046 , Knutson et al. 2014), probably aligned (true spin-orbit angle, Ψ " 25 ˝˘18 ˝, Mancini et al. 2018) hot Jupiter discovered by Bakos et al. (2011).Linear trend in the radial velocity residuals has been detected by Knutson et al. (2014), they suggest that this acceleration is an evidence of a presence of at least one additional body in the system.Later on, Piskorz et al. (2015) detected a spectroscopic companion with an effective temperature of 4000 K.However, this companion could not be seen in the AO observations (Ngo et al. 2015), based on which Piskorz et al. (2015) calculated the mass of the potential companion to be " 660 MJ with a maximum separation of 33 AU.If this companion is responsible for the RV trend, then it should have a face on orbit (e.g. the inclination of the companion's orbit must be close to 0 ˝).The companion also could have separation larger than 33 AU but observed at the time when the angular separation is low, which explains the non-AO detection.Based on the mass ratio of the host star and the companion, it is possible for companion to excite the Kozai-Lidov oscillations for HAT-P-22 b from 33 AU distance (see Fig 5 .in Ngo et al. 2015).Moreover, the small eccentricity could be a hint for such oscillation.HAT-P-22 b is also one of the most favorable exoplanet for atmospheric characterization with TSM = 582.We included HAT-P-22 b to our list to attempt to detect TTV and/or update the ephemeris for future observations.
Ground-based photometric follow-up transit observations have been carried out by Hinse et al. (2015) and Wang et al. (2021).We used all available observations from the literature, ETD, TESS and our observations, which passed our criteria (six of them were eliminated), to form a TTV diagram of 19 data points spanning a baseline of 13 years.We did not detect any parabolic or periodic changes, and we updated the ephemeris of the exoplanet HAT-P-22 b as a result.

HAT-P-30/WASP-51 System
HAT-P-30 b is a hot Jupiter (Mp " 0.711 MJ, Rp " 1.340 RJ), independently discovered by Johnson et al. (2011) andEnoch et al. (2011).RV observations show that the planet has a highly oblique (i.e.misaligned) orbit but no potential perturbing companion has been detected in the system through spectral (Piskorz et al. 2015) or AO observations (Ngo et al. 2015).Enoch et al. ( 2011) detected a strong Lithium absorption line indicating the system is young (< 1 Gyr).Therefore, it is possible that the planet has not had enough time to damp its obliquity with tidal dissipation (Winn et al. 2010a).Bai et al. (2022) detected TTV for HAT-P-30 b that could be caused by apsidal precession or additional perturbing body.We selected this system to investigate the findings of Bai et al. (2022) with the new transit observations.We analyzed a total of 48 light curves (6 of them were eliminated), including three TESS sectors and followup observations from the literature (Wang et al. 2021; Maciejewski et al. 2016a) and ETD observations to form the TTV diagram spanning more than 12 years.We were able to accurately update the ephemeris thanks to the multisector TESS observations.However, in contrast to Bai et al. (2022), we did not find any statistically significant TTVs.

HAT-P-53 System
HAT-P-53 b is a hot Jupiter (Mp " 1.484 MJ, Rp " 1.318 RJ ) that orbits a Sun-like star (Hartman et al. 2015).Unfortunately, RV follow-up observations are not sufficiently precise to measure the orbital eccentricity, RM effect or to search for additional bodies even unless they are too massive though Hartman et al. (2015) denoted the system's RV can be precisely measured despite the relatively faint host star thanks to its slow rotation and low surface temperature.Although the star rotates slowly, the planet moves rapidly in its orbit.Because of that, HAT-P-53 system is suggested to be a good example for tidal spin-up by Gallet (2020).As the angular momentum transferred from planet's orbit to star, star will rotate faster while planet's orbit shrinks.This effect could be observable with the TTV method if the observed time range is long enough as the amplitude of this effect increases within time.We selected this system to analyze its TTV to attempt to detect such variation.
Photometric transit follow-up observations have been carried out by Kjurkchieva et al. (2018) and Wang et al. (2021).The system has only 3 transit observations that cover the full transit and does not have 2-minute TESS observations.Nevertheless, we were able to update the ephemeris of the system using the combination of ETD, literature and 30-min TESS data for future observations.We did not detect any deviations from the linear ephemeris in the system in our analysis.

KELT-3 System
KELT-3 b is a hot Jupiter (Mp " 1.477 MJ, Rp " 1.345 RJ) orbiting a bright, (V " 9.8 magnitude) late F star discovered by Pepper et al. (2013).A faint nearby star at 3.74 arcseconds angular distance was detected from direct imaging observations (Wöllert & Brandner 2015;Pepper et al. 2013).Gaia revealed that this neighbor is actually bound to the system at a linear distance of " 800 au.Relatively high surface temperature and brightness of the host star make the KELT-3 system an excellent candidate for probing the atmosphere of its planet in shorter wavelengths with transit observations (e.g.Cauley et al. 2017;Corrales et al. 2021) or with the occultation observations in longer wavelengths (Emission Spectroscopy Metric, ESM = 170, Kempton et al. 2018).Despite having a bright host star, the system does not have many follow-up observations in the literature.In fact, our observation is the only one that covers a full transit.Mallonn et al. (2019) updated the ephemeris using the transit observations from ETD and observation from Pepper et al. (2013).Wang et al. ( 2021) observed two transits but they were not able to cover the full duration.
We have refined the ephemeris of KELT-3 b using two sectors of TESS observations and our observation.Based on our criteria, we eliminated the ETD and previous observations from Pepper et al. (2013).However, the TESS observations are relatively precise, thanks to the bright host star, enabling us to correct the ephemeris for future observations of this bright system.

Qatar-2 System
Qatar-2 b is a short period (P " 1.337 d) hot Jupiter (Mp " 2.487 MJ, Rp " 1.144 RJ) discovered by Bryan et al. (2012).Some of the transit light curves show star-spot occultations by the planet's disk.Mancini et al. (2014)  The short period and relatively high mass ratio of Qatar-2 b make it a potential target to observe an orbital decay (Dai et al. 2017 and references therein), which manifests itself as a parabolic change in the TTV diagram.The amplitude of this effect increases over time, making it detectable with ground-based observations.We observed the transit of Qatar-2 b after "1250 epochs later from Kepler observations but we did not detect statistically significant parabolic change.We also did not detect statistically significant periodic TTV in the system.Although the TTV diagram of Qatar-2 b is not well-sampled, the updated ephemeris pre-cision is the highest among the exoplanets in this study (except WASP-50 b), thanks to the ultra-precise Kepler data.

WASP-8 System
WASP-8 b is a warm Jupiter (P " 8.1587 d, Mp " 2.244 MJ, Rp " 1.038 RJ), orbiting a bright (V " 9.87 magnitude) solar-like star, discovered by Queloz et al. (2010).The planet is very interesting due to its eccentric (e " 0.3044; Knutson et al. 2014) and misaligned and retrograde orbit (λ " ´143 ˝; Bourrier et al. 2017).In the discovery paper, the radial velocity residuals show a linear drift, potentially caused by a companion.The system consists of a physically bound faint M-dwarf, located " 4.5 2 (" 440 au) away from WASP-8A (Ngo et al. 2015).Follow-up RV observations by Knutson et al. (2014) revealed that only a part of the observed slope in RV residuals can be due to the presence of WASP-8 B. Instead, another planet, WASP-8 c (Pc " 4323 d, Mc sin ic " 9.45 MJ) was found to be responsible for the RV variation.The only photometric follow-up observations were carried out by Borsato et al. (2021) with The CHaracterising ExOPlanet Satellite (CHEOPS, Benz et al. 2021) to improve the precision of ephemeris.
WASP-8 b is a promising TTV candidate and has a very high TSM (421), suitable for atmospheric observations.However, its equatorial position (δ " ´35 ˝) and long period (hence long transit duration) make it difficult to observe its full transits.We updated the ephemeris of WASP-8 b with three sectors of TESS observations and two light curves from previously published observations.This updated ephemeris will be useful for future ground and space-based observations of the system.

WASP-44 System
WASP-44 b is a hot Jupiter (Mp " 0.889 MJ, Rp " 1.14RJ) discovered by Anderson et al. (2012).Mancini et al. (2013) found that the radius of the planet is smaller by 10% than first measured and there is no extreme radius variation in the optical wavelengths from multi-band photometry.However, Turner et al. (2016) reported the radius of the planet is 1.4σ larger in the near-ultraviolet.After "6.5 years, Addison et al. (2019) observed a transit and updated the ephemeris.Similiar to the HAT-P-30 system, WASP-44 system is also exposed to tidal spin-up (Gallet 2020).Similiar study has been carried out by Brown (2014) and isochrone age was found to be significantly older than gyrochronological age.The angular momentum transfer from planet's orbit to rotation of the star could manifest itself as orbital decay in TTV diagram.We selected this system to attempt to detect such variation.
We analyzed the follow-up transit observations mentioned above, along with our own observations, TESS and ETD data, to update the ephemeris.However, we did not find any evident periodic or secular TTVs in its timing data.

WASP-50 System
WASP-50 b is a hot Jupiter (Mp " 1.468 MJ, Rp " 1.53 RJ), discovered by Gillon et al. (2011) revolving on a circular orbit (Bonomo et al. 2017).Follow-up photometric transit ob-servations were carried out by Tregloan-Reed & Southworth (2013), Sada et al. (2012), andSada (2018) to update the ephemeris or increase the precision of its transit parameters.Gillon et al. (2011) measured the rotation period of WASP-50 from two seasons of WASP photometry as 16.3 ˘0.5 days; however, Canto Martins et al. (2020) found this value to be only 5.488 days from TESS sector-4 light curve.We performed a preliminary analysis of the sector 31 PDCSAP FLUX of TESS and confirmed the finding by Canto Martins et al. (2020).The measured log R 1 HK and Prot values by Gillon et al. (2011) are in excellent agreement with each other comparing to the empirical values calculated using the Prot -log R 1 HK relation presented by Suárez Mascareño et al. (2015).Moreover, a rotation rate of 5.488 days indicates a very young age (" 80 Myr; Barnes 2007); however, according to the lithium abundance, the system should be at least 0.6 ˘0.2 Gyr old (Gillon et al. 2011).If the true rotation period is 5.488 days, then the lack of lithium suggests that the star could be a good example of tidal spin-up (e.g.Gallet 2020).Tejada Arevalo et al. (2021) have suggested that even after orbital circularization, the planet's orbit may shrink by transferring angular momentum to its host star and causing its rotation rate to increase.We selected this system to observe such effect via TTV method as it should manifest itself as orbital decay.
We analyzed 45 light curves (9 of them were eliminated) spanning 10 years of best observations available.We did not detect any parabolic TTV but the frequency analysis peaked at 34.45 days with a false alarm probability of 2 per cent.The amplitude of this periodic variation is 57 seconds, which is compatible with our average mid-transit uncertainty.As a result, our findings are inconclusive.Further precise observations are required to confirm this hint of a periodic TTV.

WASP-77 System
WASP-77A b is a short period (P " 1.36 days) hot Jupiter (Mp " 1.76 MJ, Rp " 1.21 RJ) revolving around a G8 V type, bright (V " 10.12 magnitude), wide binary with the component WASP-77B at a projected angular distance of " 3.5 arcseconds (Maxted et al. 2013).Photometric followup transit observations that confirmed the transit parameters of the discovery paper were carried out by Turner et al. (2016) and Cortés-Zuleta et al. (2020).The planet has relatively high TSM and ESM (ESM = 333, TSM = 770) and its wide companion (WASP-77B) can be used as a comparison star, making it favorable for atmospheric observations via transmission or emission spectroscopy from the ground (Line et al. 2021;Reggiani et al. 2022) or space-borne observations (Mansfield et al. 2022).Gallet (2020) suggested the host star might have been affected by tidal-spin-up by its planet WASP-77 b.
Cortés-Zuleta et al. (2020) performed a TTV analysis for WASP-77 b in a similar way within this work.We added additional transit light curves from TESS sector-31, our observations and the newly available light curves from ETD.As a result, we were able to update the ephemeris with increased precision, thanks to transit light curves covering a longer baseline.As in Cortés-Zuleta et al. ( 2020), we found no significant secular or periodic TTVs.

WASP-93 System
WASP-93 b is a hot Jupiter (Mp " 1.47 MJ, Rp " 1.597 RJ) orbiting a fast-rotating (v sin i " 37 ˘3 km s ´1) F4 V star discovered by Hay et al. (2016).RM observations were attempted twice by Hay et al. (2016); however, the first observation was unable to cover the transit due to ephemeris uncertainty, and the combination of the first and second attempts resulted in inconclusive results due to insufficient RV precision.Gajdoš et al. (2019) searched for TTVs in the system using only ETD observations and they did not observe any significant deviation from linear ephemeris.Although WASP-93 b has relatively high TSM and ESM, TESS observations do not show significant phase modulations or an occultation signal (Wong et al. 2021).
TESS observed WASP-93 during sectors 17, 57, and 58 but, unfortunately, the object was too close to the edge of the camera in sector-58 observations, hence we were unable to use it.We observed 7 transits of WASP-93 b and used available observations in the literature and ETD to update its ephemeris.The RM observations by Hay et al. (2016) show the hint of a retrograde orbit if the transit ephemeris arrived earlier by "35 min.The timing difference between our ephemeris and the ephemeris from Hay et al. ( 2016) is only ´46 ˘48 seconds at the time of the second RM observations.Therefore, we rule out the early transit-retrograde orbit scenario even if the two independent RM observations agree with each other (see Figure 8 in Hay et al. 2016).We also did not detect any secular or periodic TTV signal.

XO-2N System
XO-2N b is a hot Jupiter (Mp " 0.57 MJ, Rp " 0.98 RJ) orbiting around a metal rich, rM{Hs " 0.44 ˘0.02 dex, wide binary component XO-2N in an aligned orbit (Narita et al. 2011) discovered by Burke et al. (2007).A binary companion, XO-2S, which is also a metal-rich star, resides " 30 2 away from XO-2N and has at least two planets discovered with RV observations (Desidera et al. 2014).The visual binary components have similar effective temperatures, thus XO-2S is a great comparison star for transmission spectroscopy of XO-2N b (e.g.Sing et al. 2012;Crouzet et al. 2012).With additional RV observations, Knutson et al. (2014) detected that the radial velocity residuals show a linear trend within time, possibly caused by an outer companion.Later on, Damasso et al. (2015) revealed that the linear RV residuals are actually only a part of the curve in RV, possibly caused by an outer companion XO-2N c or by stellar activity.
We obtained the transit light curves from several works (Burke et al. 2007;Fernandez et al. 2009;Kundurthy et al. 2011;Damasso et al. 2015;Maciejewski et al. 2018;Wang et al. 2021), three sectors of observation (20, 47, 60) from TESS, 13 from amateur astronomers and our three new observations to form the TTV diagram with a total of 42 light curves, spanning almost 16 years.Our analyses indicated that the parabolic ephemeris fit the data better than linear with the ∆BIC of 42.98 and ∆AIC of 45.15, suggesting the orbit of XO-2 b is decaying with a rate of ´12.95 ˘1.85 ms yr ´1 (Figure 3).The parabolic ephemeris is as follows: Although the parabolic ephemeris is statistically significant, it does not agree well with the latest TESS observations (sectors 47 and 60).TESS will observe the XO-2N system during cycle 6 of its mission.However, the wide binary component is only " 30 2 away, and a TESS pixel has a 21 2 field of view.This may add extra red and white noise because the light from XO-2S blends into the aperture selected for XO-2N.Therefore, ground-based observations could be a better option for confirming or discarding the parabolic trend.

DISCUSSION
We constructed the TTV diagrams of 16 exoplanets consisting of the most precise and complete light curves with the longest time span for each of the planets in our sample (Figure 1).This allowed us to increase the precision of the orbital period and the accuracy of the ephemeris information for future follow-up observations.Based on the ephemeris information given in Table 3, the uncertainty on the predicted ephemeris will provide transit timings with a precision below 5 minutes with an average of 1.8 minutes until 2070 for all the systems except for XO-2N and HAT-P-13 where we detect deviation from the linear ephemeris.
We detect a decrease in the orbital period of XO-2N b may have been caused by several events if it is real.As discussed by Vissapragada et al. (2022), such a decrease can be observed if the system is accelerating towards us, making the transits observed earlier than expected.In that case, the radial velocity residuals should have a slope of ´0.05 m s ´1 d ´1 ( 9 ν = c 9 P {P , Vissapragada et al. 2022) but Damasso et al. (2015) reported this value as `0.0017 m s ´1 d ´1.Damasso et al. (2015) also reported the eccentricity of the planet is consistent with zero.Therefore we did not consider a scenario based on precession.The RV residuals show a parabolic variation that can be caused by a long-period outer companion.This companion may be causing LiTE and changing the orbital period of XO-2 b within a longer time interval, thus the parabolic TTV could be a part of this periodic variation.We neglect this scenario also because Damasso et al. (2015) showed that the parabolic RV variation is due to magnetic activity.Even if the magnetic activity was not the reason for the RV residuals, the phases of LiTE and long-term parabolic RV change do not match.
We detect a cyclic variation in the period of HAT-P-13 b which has a semi-amplitude of 160.8 seconds and 479.52 day periodicity.This variation could be caused by the known outer companion of the system, HAT-P-13 c (P " 445.82 ˘0.11 days, Mp sin i " 14.61 `0.46 ´0.48 MJ, Knutson et al. 2014).However, in order to cause such a high-amplitude TTV, the orbital inclination of HAT-P-13 c needs to be " 2 ˝, which translates into a mass of " 0.4 Md using the Mp sin i value.Piskorz et al. (2015) detected another star in the spectrum of HAT-P-13 that has an effective temperature of 3900 `300 ´350 (Mcompanion " 0.6 `0.086 ´0.179 Md).Piskorz et al. (2015) discussed if HAT-P-13 d, detected by linear drift in the RV residuals, is that spectral companion.But in that case, the inclination of the planet d should be " 5 ˝, which could lead to Kozai-Lidov oscillations and make the orbit of HAT-P-13 b misaligned Winn et al. (2010a).Our TTV analysis suggests that the observed spectral companion could be HAT-P-13 c instead.This could explain the non-AO detection by Ngo et al. (2015) due to short angular separation.However, such a close companion would have catastrophic effects on the system stability aside from making the orbit of HAT-P-13 b misaligned.
Nevertheless, in order to assess the timescale of any potential orbital instability of the planet c with an orbital inclination of 2 ˝, we run an N-body simulation by using Rebound code (Rein & Liu 2012), with IAS15 (Rein & Spiegel 2015) integrator.We used the masses and orbital parameters of the planets b and c derived from Winn et al. (2010a).We stress that we did not include the potential additional planet (planet d in Knutson et al. (2014)), since the orbital parameters and mass of this putative object are currently unknown.Results of our simulation show that the system becomes unstable on the order of 30,000 yrs.Therefore, we think that the amplitude of the detected TTV signal in our analysis might have been overestimated due to the high scatter caused by correlated noise in the transit light curves.
We do not find any statistically significant periodicities in our timing analysis of any other system.However, we simulated the WASP-8 system in a same way as we did for HAT-P-13 and we found that the TTV due to WASP-8 c should be observable.We assumed the transiting planet WASP-4 b is coplanar with the RV planet WASP-8 c.Using the absolute and orbital parameters from Queloz et al. (2010) and Knutson et al. (2014), we found that the full TTV amplitude of WASP-8 b should be "40 seconds due to LiTE and "2 seconds due to gravitational interaction.The typical midtransit measurement error of TESS data for this system is about 20 seconds, so the TTV of this system can be detected with light curves that has similar precision as TESS.We could not detect this signal due to poor phase coverage.

Discussion on Tidal Quality Factors
Although we do not have any observational data on the rotation rates of the exoplanets in our sample, the rotational periods of their hosts are all longer compared to their orbital periods.The energy raised in those tidal interactions can be dissipated in the convective envelope of their host stars, transferring angular momentum from the planet to the star that would cause the star to spin-up while the planet to fallin (Counselman 1973;Rasio et al. 1996;Essick & Weinberg 2016;Harre et al. 2023;Weinberg et al. 2023).In addition to this equilibrium tide, dynamical tide excites internal gravity waves, which dissipate the energy through also the secondary waves it generates via wave-wave interactions (Barker & Ogilvie 2010;Ivanov et al. 2013;Barker 2020).This latter mechanism is especially dominant in a system with a solar-type host and a hot Jupiter-type planet (Essick & Weinberg 2016;Weinberg et al. 2023).
If we assume the planetary mass to be constant, then the rate of change in the orbital period can be related to the so-called reduced tidal quality factor of the host star by the constant phase lag model of Goldreich & Soter (1966) defined as where Mp is the mass of the planet, M‹ and R‹ are the mass and the radius of the host star and a is the semi-major axis of the planet's orbit.The rate of orbital decay, 9 P, is derived from the timing analysis, which is twice the value of the quadratic coefficient of the best-fitting parabola.If that best-fitting model is not found to be statistically superior to the linear model, then an orbital decay cannot be argued and the quadratic coefficient can only be used to derive a lower limit for the reduced tidal quality factor, which is the case for all systems in our sample except XO-2N.We derived these limits based on the fundamental parameters of the objects in our sample, which we provide in Table 4 together with ∆AIC and ∆BIC values, indicating the statistical significance of the quadratic model in each of the cases.Positive values for both statistics hint that the quadratic model should be favored.
XO-2N b is not one of the prime candidates for orbital decay due to tidal interactions with its host star because it is not a particularly big planet.Although its solar-like host star should have a convective envelope to dissipate the tidal energy, its reduced tidal quality factor should be larger than 550 as derived from the rate of observed period decrease when it is compared with other host stars with similar spectral types and evolutionary history.Infall time calculated from the same rate is P { 9 P « 17.45 Myr; which is too short compared to age of the system.TESS observations do not follow the orbital decay model found to be statistically superior to the linear model.However, the precision of TESS observations of XO-2N is lower than that achieved for the stars with similar magnitudes.This is because the wide binary component XO-2S is only " 30 2 away from the XO-2N, blending its flux in the TESS aperture, adding extra noise and diluting the transit depth, resulting in high scatter on the TTV diagram.Because of this reason, we encourage observations of the transit of XO-2N b with ground-based telescopes in high angular resolution in the future to confirm the orbital decay scenario, which we find unlikely at the moment.

Figure 1 .
Figure 1.Linear residuals of TTV diagrams for all the planets in our sample based on observations from open databases (green), our observations (red), TESS observations (magenta), Kepler observations (yellow) and light curves published in the literature (blue).

Figure 2 .
Figure 2. Top: Lomb-Scargle periodogram of TTV of HAT-P-13 b.The horizontal dotted lines correspond to false alarm probabilities of 0.1, 0.01, 0.001, from bottom to top.Bottom: phase folded TTV to the frequency with the highest power.
observed consecutive transits of QATAR-2 b with ground-based telescopes and they concluded that the planet's orbital plane is aligned with the stellar rotation by tracking the change in position of one star-spot.Later on, Dai et al. (2017) and Močnik et al. (2017) used Kepler observations and also found that the planet's orbit is aligned.The findings were further confirmed by Esposito et al. (2017) based on radial velocity observations during a transit, which revealed a symmetric Rossiter-McLaughlin effect in the prograde direction.

Figure 3 .
Figure 3. Quadratic TTV model of XO-2N b (black dashed line with 3 σ uncertainty shown in grey shaded region).The linear term was subtracted for display purposes.

Table 1 .
Fundamental stellar and planetary properties and the number of light curves analyzed for each planetary system in our sample.

Table 2 .
The log of photometric observations performed for this study.The dates of the light curves that are eliminated and hence not used in the TTV diagrams are marked and the reasons for their elimination are given in the footnotes.

Table 4 .
Lower limits for the Reduced Tidal Quality Factors (Q 1 ‹ ) for the host stars in our sample at 95% confidence level.∆AIC and ∆BIC values were used for comparisons between linear and quadratic models in this order.The quadratic model is favored in the cases where ∆BIC > 10.This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI).Funding for the Kepler mission is provided by the NASA Science Mission Directorate.STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement.