# Examining the orbital decay targets KELT-9 b, KELT-16 b, and WASP-4b, and the transit-timing variations of HD 97658 b

J. -V. Harre, A. M. S. Smith, S. C. C. Barros, G. Boué, Sz. Csizmadia, D. Ehrenreich, H. -G. Florén, A. Fortier, P. F. L. Maxted, M. J. Hooton, et al.

► To cite this version:

J. -V. Harre, A. M. S. Smith, S. C. C. Barros, G. Boué, Sz. Csizmadia, et al., Examining the orbital decay targets KELT-9 b, KELT-16 b, and WASP-4b, and the transit-timing variations of HD 97658 b. Astronomy and Astrophysics - A&A, 2023, 669, 10.1051/0004-6361/202244529. insu-04155780

HAL Id: insu-04155780

https://insu.hal.science/insu-04155780

Submitted on 7 Jul 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Examining the orbital decay targets KELT-9 b, KELT-16 b, and WASP-4 b, and the transit-timing variations of HD 97658 b*,**

J.-V. Harre¹, A. M. S. Smith¹, S. C. C. Barros²,³, G. Boué⁴, Sz. Csizmadia¹, D. Ehrenreich⁵,
H.-G. Florén⁶, A. Fortier⁷,⁸, P. F. L. Maxted⁹, M. J. Hooton¹⁰,¹¹, B. Akinsanmi²,⁵, L. M. Serrano¹²,
N. M. Rosário²,³, B.-O. Demory⁸, K. Jones⁸, J. Laskar⁴, V. Adibekyan², Y. Alibert⁷, R. Alonso¹³,¹⁴,
D. R. Anderson¹⁵, G. Anglada¹⁶,¹⁷, J. Asquier¹⁸, T. Bárczy¹⁹, D. Barrado y Navascues²⁰,
W. Baumjohann²¹, M. Beck⁵, T. Beck⁷, W. Benz⁷,⁸, N. Billot⁵, F. Biondi²²,²³, A. Bonfanti²¹, X. Bonfils²⁴,
A. Brandeker⁶, C. Broeg⁷,²⁵, J. Cabrera¹, V. Cessa⁷, S. Charnoz²⁶, A. Collier Cameron²⁷, M. B. Davies²⁸,
M. Deleuil²⁹, L. Delrez³⁰,³¹, O. D. S. Demangeon²,³, A. Erikson¹, L. Fossati²¹, M. Fridlund³²,³³,
D. Gandolfi¹², M. Gillon³⁰, M. Güdel³⁴, C. Hellier⁹, K. Heng⁸,³⁵, S. Hoyer²⁹, K. G. Isaak³⁶,
L. L. Kiss³⁷,³⁸, A. Lecavelier des Etangs³⁹, M. Lendl⁵, C. Lovis⁵, A. Luntzer⁴⁰, D. Magrin²²,
V. Nascimbeni²², G. Olofsson⁶, R. Ottensamer⁴⁰, I. Pagano⁴¹, E. Pallé¹³, C. M. Persson³³, G. Peter⁴²,
G. Piotto²²,⁴³, D. Pollacco³⁵, D. Queloz⁴⁴,¹⁰, R. Ragazzoni²²,⁴³, N. Rando¹⁸, H. Rauer¹,⁴⁵,⁴⁶, I. Ribas¹⁶,¹⁷,
G. R. Ricker⁴⁷,⁴⁸, S. Salmon⁵, N. C. Santos²,³, G. Scandariato⁴¹, S. Seager⁴⁹,⁵⁰,⁵¹, D. Ségransan⁵,
A. E. Simon⁷, S. G. Sousa², M. Steller²¹, Gy. M. Szabó⁵²,⁵³, N. Thomas⁷, S. Udry⁵,
B. Ulmer⁴², V. Van Grootel³¹, N. A. Walton⁵⁴, T. G. Wilson²⁷, J. N. Winn⁵⁵, and B. Wohler⁵⁶,⁵⁷

(Affiliations can be found after the references)

Received 18 July 2022 / Accepted 3 November 2022

ABSTRACT

**Context.** Tidal orbital decay is suspected to occur for hot Jupiters in particular, with the only observationally confirmed case of this being WASP-12 b. By examining this effect, information on the properties of the host star can be obtained using the so-called stellar modified tidal quality factor *Q'*, which describes the efficiency with which the kinetic energy of the planet is dissipated within the star. This can provide information about the interior of the star.

**Aims.** In this study, we aim to improve constraints on the tidal decay of the KELT-9, KELT-16, and WASP-4 systems in order to find evidence for or against the presence of tidal orbital decay. With this, we want to constrain the *Q'* value for each star. In addition, we aim to test the existence of the transit timing variations (TTVs) in the HD 97658 system, which previously favoured a quadratic trend with increasing orbital period.

**Methods.** Making use of newly acquired photometric observations from CHEOPS (CHaracterising Exoplanet Satellite) and TESS (Transiting Exoplanet Survey Satellite), combined with archival transit and occultation data, we use Markov chain Monte Carlo (MCMC) algorithms to fit three models to the data, namely a constant-period model, an orbital-decay model, and an apsidal-precession model.

**Results.** We find that the KELT-9 system is best described by an apsidal-precession model for now, with an orbital decay trend at over 2σ being a possible solution as well. A Keplerian orbit model with a constant orbital period provides the best fit to the transit timings of KELT-16 b because of the scatter and scale of their error bars. The WASP-4 system is best represented by an orbital decay model at a 5σ significance, although apsidal precession cannot be ruled out with the present data. For HD 97658 b, using recently acquired transit observations, we find no conclusive evidence for a previously suspected strong quadratic trend in the data.

**Key words.** planets and satellites: dynamical evolution and stability – planet-star interactions – techniques: photometric

**1. Introduction**

Hot Jupiters, as they are so close to their host stars, are expected
to tidally interact with them. This interaction is expressed as the
planet raising a tidal bulge on the surface of the star due to the
gravitational attraction with the mass of the planet. Vice versa,

the star also raises a bulge on the surface of the planet, but this
latter has a negligible influence on the effect we examine here,
that is, on the tidal decay of the orbit of hot Jupiters. If the stellar
rotation is not synchronized with the planetary orbital period,
the viscosity of the star's plasma leads to a lag between the tide
raised on the surface of the star and the planet-to-star centre,
which leads to a transfer of orbital angular momentum from the
planet to the star if the orbital period of the planet is smaller than
the stellar rotation period. This is described as the equilibrium
tide. This subsequently causes the planet to slowly spiral inwards
and the star to spin up (Counselman 1973; Rasio et al. 1996). An
additional contribution is made from the dynamical tide, which
arises from stellar oscillations (Ogilvie 2014).

<footnote>* Photometry data according to Table 1, as well as the full Tables 4–7 and Table B.1 are only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/669/A124</footnote>

<footnote>** This article uses data from CHEOPS programme CH_PR100012 and CH_PR100013.</footnote>
The sole planet for which tidal orbital decay has been obser-
vationally confirmed is WASP-12 b. Discovered in the Wide
Angle Search for Planets (WASP) project (Pollacco et al. 2006)
by Hebb et al. (2009), WASP-12 b orbits its host star – a late
F-type star – every 1.09 days, has a mass of 1.47 M<sub>J</sub>, and a
radius of 1.90 R<sub>J</sub>. The variations in transit and occultation tim-
ings observed for this planet lead to a decrease in its orbital
period of 29 ms yr<sup>-1</sup>, an orbital decay timescale P/P<sub>0</sub> = 3.25 Myr,
and finally to a modified stellar quality factor of Q'<sub>*</sub> = 1.8 × 10<sup>5</sup>,
as found by Yee et al. (2020). Q'<sub>*</sub> describes the efficiency with
which orbital kinetic energy is dissipated within the star due
to friction; the smaller its value, the stronger the dissipation
(Goldreich & Soter 1966). There may also be a dependence of
Q'<sub>*</sub> on the tidal forcing period P<sub>tide</sub>, as suggested by the results of
Penev et al. (2018).

The observation of hot Jupiters enables us to gain infor-
mation on the planet–star interactions within the system. We
are examining tidal interactions between hot Jupiters and their
host stars by measuring tidal orbital decay and also tidal defor-
mation in the Guaranteed Time Observing (GTO) programme
of the Characterising Exoplanet Satellite (CHEOPS) mission,
under the “Tidal decay (ID 0012)” programme (Benz et al.
2021; Barros et al. 2022). This paper deals with the tidal orbital
decay candidates KELT-9 b, discovered by Gaudí et al. (2017),
KELT-16 b, discovered by Oberst et al. (2017), WASP-4 b, dis-
covered by Wilson et al. (2008), and also the potential transit
timing variation (TTV) system HD 97658 discovered by Howard
et al. (2011). Most of the short-period targets with the highest
expected tidal-decay rates were discovered by ground-based sur-
veys such as that of the Kilodegree Extremely Little Telescope
(KELT) and WASP. The KELT survey (Pepper et al. 2007) was
designed to find close-in orbiting giant planets and concluded in
2020 after 14 yr of observations and more than 20 discovered
exoplanets<fnref n="1" />. The WASP survey (Pollacco et al. 2006) started
in 2004 with SuperWASP-North, and WASP South joining in
2006. This survey discovered more than 180 exoplanets<fnref n="2" />. Our
targets WASP-4 b, KELT-9 b, and KELT-16 b are hot Jupiters
with orbital periods of less than 1.5 days and masses greater than
1.2 times that of Jupiter, making them fall into the category of
potential orbital-decay candidates. The suspected TTV candi-
date HD 97658 b is a super-Earth that orbits its host star every
9.5 days.

Measuring the orbital-decay rate of our candidates allows us
to constrain their *Q*′\* values and to compare them against values
from theory in the literature. Values given in the literature for the
value of the modified stellar quality factor cover a broad range
from 10⁵ up to 10⁸.⁵ (e.g. Meibom & Mathieu 2005; Jackson
et al. 2008; Penev et al. 2012).

Amongst the space telescopes capable of observing exo-
planet transits and occultations photometrically, two are partic-
ularly suited to our needs, namely CHEOPS and the Transiting
Exoplanet Survey Satellite (TESS), while also offering the high
precision necessary for measuring tidal decay. The CHEOPS
space mission (Benz et al. 2021) is an S-class mission in the sci-
ence programme of the European Space Agency (ESA) and was
launched on 18 December 2019. This space telescope is designed
to carry out follow-up observations of planets transiting bright
stars, delivering high-precision photometry to improve system
parameters where applicable. It consists of a single defocused
32 cm telescope, which is used to observe a single target at a
time. We use CHEOPS to acquire new photometric data for our

targets, as it can provide precise timings from its transit and sec-
ondary eclipse observations (Lendl et al. 2020; Borsato et al.
2021; Hooton et al. 2022; Deline et al. 2022; Brandeker et al.
2022; Barros et al. 2022).

In addition to CHEOPS, TESS observations are used to get
precise timings for our analysis. This satellite was launched on
18 April 2018. Unlike CHEOPS, TESS is a survey satellite,
oberving many targets at once during one of its observation sec-
tors. It consists of four cameras, each with a 10.5 cm entrance
pupil diameter lens assembly in front of four CCD detectors
Ricker et al. (2014), also making use of the transit method. There
are more than 220 confirmed exoplanets discovered with TESS,
with almost 5 800 project candidates to date<fnref n="3" /> (Akeson et al. 2013;
Guerrero et al. 2021).

This paper presents observations made to constrain the tidal
decay of WASP-4 b, KELT-9 b, and KELT-16 b and to investi-
gate the apparent TTVs in the HD 97658 b system presented in
Maxted et al. (2022). In Sect. 2, we describe the observations
from the literature and the new observations used for the anal-
ysis described in this work. The methods used in this paper are
described in more detail in Sect. 3, the results from our analysis
are shown and discussed in Sect. 4, and the conclusions can be
found in Sect. 5.

## 2. Observations

2.1. Targets

KELT-9 is a bright, fast-rotating star at the border of the B/A-
types with an effective temperature of around 10 000 K, and
a mass of 1.98±0.02 M<sub>☉</sub>, hosting an inflated hot Jupiter on a
1.48 days polar orbit with a mass of 2.88 ± 0.84 M<sub>J</sub> and a radius
of 1.89<sup>+0.06</sup><sub>-0.05</sub> R<sub>J</sub> (Gaudí et al. 2017; Hoeijmakers et al. 2019).
KELT-9 b is highly irradiated with a dayside temperature of more
than 4000 K. An important feature to be considered in this sys-
tem is gravity darkening caused by the fast rotating nature of the
star. The fast stellar rotation leads to the equatorial radius of the
star being larger as compared to the polar radius due to centrifu-
gal forces in the star. This causes an oblate shape, with higher
surface temperatures at the poles and lower temperatures at the
equator, modifying the shape of transits in light curves (Barnes
2009; Ahlers et al. 2020; Cauley & Ahlers 2022). KELT-9 was
amongst the targets of a recent study by Ivshina & Winn (2022),
who found no evidence for orbital decay in this system.

KELT-16 is an F7V star with a brightness of V = 11.7 and an effective temperature of around 6250 K. It has a mass of 1.21<sup>+0.04</sup><sub>-0.05</sub> M<sub>☉</sub>, a radius of 1.36<sup>+0.06</sup><sub>-0.05</sub> R<sub>☉</sub> and hosts a hot Jupiter in a 0.97 day orbit. The planet KELT-16 b has a mass of 2.75<sup>+0.16</sup><sub>-0.15</sub> M<sub>J</sub> and is inflated with a radius of 1.415<sup>+0.084</sup><sub>-0.067</sub> R<sub>J</sub> (Oberst et al. 2017). Its equilibrium temperature is around 2450 K. There have been several studies examining its potential orbital decay, including those of Maciejewski et al. (2018), Patra et al. (2020), Wong et al. (2021), Mancini et al. (2022), and Ivshina & Winn (2022), all of which found no significant period changes and hence no significant signs of orbital decay. A recent analysis by Mancini et al. (2022) found a period derivative of Ṗ = −10.6 ± 13.1 ms yr<sup>-1</sup>, which is consistent with a constant period, and hence no orbital decay. From this, the authors were able to deduce a lower limit on the modified stellar tidal quality factor Q'<sub>*</sub> > (1.9 ± 0.8) × 10<sup>5</sup>.

WASP-4 b is a hot Jupiter with an orbital period of 1.34 days,
a mass of 1.186<sup>+0.090</sup><sub>-0.098</sub> M<sub>J</sub>, a radius of 1.321 ± 0.039 R<sub>J</sub>, and an

<footnote n="1">¹ According to exoplanet.eu as of February 28, 2022.</footnote>

<footnote n="2">² See footnote 1.</footnote>

<footnote n="3">³ According to the NASA Exoplanet Archive, visited on June 30,
2022.</footnote>
<table><caption>**Table 1.** Observation log of new CHEOPS transit observations for all targets.</caption><thead><tr><th>Target</th><th>Visit no.</th><th>Start date (UTC)</th><th>Duration (h)</th><th>No. of data points</th><th>Eff. (%)</th><th>File key</th></tr></thead><tbody><tr><td>KELT-9 b</td><td>1</td><td>2020-08-24 14:10:59</td><td>12.11</td><td>728</td><td>61.1</td><td>CH_PR100013_TG001001_V0200</td></tr><tr><td rowspan="7">KELT-16 b</td><td>1</td><td>2020-08-18 22:23:41</td><td>7.95</td><td>325</td><td>68.0</td><td>CH_PR100012_TG001401_V0200</td></tr><tr><td>2</td><td>2020-09-06 07:44:45</td><td>8.75</td><td>364</td><td>69.2</td><td>CH_PR100012_TG001402_V0200</td></tr><tr><td>3</td><td>2022-07-15 16:23:21</td><td>10.99</td><td>469</td><td>71.0</td><td>CH_PR100012_TG000501_V0200</td></tr><tr><td>4</td><td>2022-08-25 08:26:22</td><td>8.59</td><td>352</td><td>69.0</td><td>CH_PR100012_TG000502_V0200</td></tr><tr><td>5</td><td>2022-09-13 17:32:22</td><td>8.14</td><td>262</td><td>53.6</td><td>CH_PR100012_TG000503_V0200</td></tr><tr><td>6</td><td>2022-09-25 09:27:12</td><td>7.50</td><td>284</td><td>62.9</td><td>CH_PR100012_TG000504_V0200</td></tr><tr><td>7</td><td>2022-10-01 03:59:13</td><td>8.14</td><td>282</td><td>57.6</td><td>CH_PR100012_TG000505_V0200</td></tr><tr><td rowspan="8">WASP-4 b</td><td>1</td><td>2021-07-15 10:19:20</td><td>5.80</td><td>203</td><td>58.1</td><td>CH_PR100012_TG001101_V0200</td></tr><tr><td>2</td><td>2021-08-09 20:22:19</td><td>5.74</td><td>221</td><td>64.0</td><td>CH_PR100012_TG001102_V0200</td></tr><tr><td>3</td><td>2021-08-17 20:21:21</td><td>6.35</td><td>235</td><td>61.5</td><td>CH_PR100012_TG001103_V0200</td></tr><tr><td>4</td><td>2021-08-31 05:32:21</td><td>6.35</td><td>241</td><td>63.1</td><td>CH_PR100012_TG001104_V0200</td></tr><tr><td>5</td><td>2021-09-08 06:20:21</td><td>6.25</td><td>234</td><td>62.2</td><td>CH_PR100012_TG001105_V0200</td></tr><tr><td>6</td><td>2021-09-22 23:48:21</td><td>6.07</td><td>278</td><td>76.1</td><td>CH_PR100012_TG001106_V0200</td></tr><tr><td>7</td><td>2021-10-14 09:58:20</td><td>6.05</td><td>227</td><td>62.3</td><td>CH_PR100012_TG001107_V0200</td></tr><tr><td>8</td><td>2022-08-08 20:15:20</td><td>6.32</td><td>220</td><td>57.9</td><td>CH_PR100012_TG001001_V0200</td></tr><tr><td rowspan="3">HD 97658 b</td><td>1</td><td>2021-03-20 06:54:54</td><td>10.12</td><td>732</td><td>66.2</td><td>CH_PR100012_TG001901_V0200</td></tr><tr><td>2</td><td>2021-03-29 18:28:55</td><td>10.31</td><td>754</td><td>67.0</td><td>CH_PR100012_TG001902_V0200</td></tr><tr><td>3</td><td>2022-01-27 10:16:56</td><td>10.05</td><td>625</td><td>56.9</td><td>CH_PR100012_TG002201_V0200</td></tr></tbody></table>

Notes. “Eff.” stands for efficiency. The number of data points and the efficiency do not take data points flagged by the DRP into account.

equilibrium temperature of around 1700 K. This planet orbits a 5400 K G7V star of V-magnitude 12.5 with a stellar mass of $0.97_{-0.09}^{+0.15} M_{\odot}$ and a radius of $0.89_{-0.06}^{+0.05} R_{\odot}$, and was actually the first exoplanet discovered by the WASP South observatory (Wilson et al. 2008; Bouma et al. 2019). For the WASP-4 system, there are several studies examining potential tidal orbital decay in this system, namely Bouma et al. (2019, 2020), Baluev et al. (2019), Southworth et al. (2019), Baluev et al. (2020), Turner et al. (2022), Maciejewski (2022), and Ivshina & Winn (2022). Of these, Bouma et al. (2020) find the orbital period of WASP-4 b to be decreasing at a rate of $-8.64 \pm 1.26$ ms yr$^{-1}$, but about $-6$ ms yr$^{-1}$ of this to be due to a Doppler effect caused by the movement of the system towards the Earth, as derived from radial velocity measurements. Later that year, Baluev et al. (2020) examined this alleged radial velocity trend, but were not able to recreate it using their new data. However, the authors find a decreasing period at $\dot{P} = -5.4 \pm 1.5$ ms yr$^{-1}$ using previously unpublished data from Hoyer et al. (2013) and Huitson et al. (2017). Turner et al. (2022) find a decay rate of $-7.33 \pm 0.71$ ms yr$^{-1}$ using data from the latest TESS sectors for this target, and suspect an additional planet candidate with an orbital period of around 7000 days and a mass of approximately 5.5 Mj. The newest papers regarding this topic, from Maciejewski (2022) and Ivshina & Winn (2022), find a lower period decay rate of $\dot{P} = -4.8 \pm 1.4$ ms yr$^{-1}$ and $\dot{P} = -5.81 \pm 1.58$ ms yr$^{-1}$ using data covering the same time as Turner et al. (2022), but excluding the Hoyer et al. (2013) and Huitson et al. (2017) light curves for the former, and all available light curves including an additional sector of TESS data for the latter.

The HD 97658 system contains a K1V dwarf star with an effective temperature of around 5200 K, stellar mass of $0.85 \pm 0.08 M_{\odot}$, and a stellar radius of $0.728 \pm 0.008 R_{\odot}$. HD 97658 also contains a dense super-Earth of $7.62 \pm 0.42 M_{\oplus}$ with a radius of $2.293 \pm 0.07 R_{\oplus}$ transiting every 9.49 days (Howard et al. 2011; Ellis et al. 2021; Maxted et al. 2022).

At first, this planet was only detected using radial velocity data from Howard et al. (2011), with Henry et al. (2011) claiming the observation of a transit. The planetary transit was later challenged as Dragomir et al. (2012) were not able to re-observe it using the previously derived parameters a year later. However, using the MOST space telescope (Walker et al. 2003) to observe HD 97658 in 2012 and 2013, Dragomir et al. (2013) were able to observe five transits of HD 97658 b at different times to those claimed by Henry et al. (2011) and confirmed the transiting nature of this planet. The newest study including this system was carried out by Maxted et al. (2022) using archival literature data and also a newly acquired CHEOPS transit observation. Using these data, these latter authors find unexplained apparent TTVs, resulting in a quadratic ephemeris giving a lower BIC value of 37.9 in comparison to 55.5 for the linear model. We examine this system with three additional CHEOPS visits and new TESS Sector 49 data to help us to differentiate between the models.

## 2.2. New CHEOPS observations

For all of our targets, we provide at least one new observation made with the CHEOPS space telescope. In total, we obtained 19 new transit light curves with CHEOPS. Details regarding these can be found in the observation log in Table 1. The ‘OPTIMAL’ aperture was used for every light curve. The exposure times are 36.7 s for KELT-9 b, 60 s for KELT-16 b and WASP-4 b, and 33 s for HD 97658 b. The efficiency describes the percentage of time on the target spent collecting data. With CHEOPS, this is usually less than 100% because its orbital configuration leads to the Earth blocking the line-of-sight from CHEOPS to the target during parts of the satellite’s orbit. All of these light curves can be accessed via the DACE (Data & Analysis Center for Exoplanets) website<fnref n="4" /> of the University of Geneva. Moreover, the light

<footnote n="4">⁴ https://dace.unige.ch/cheopsDatabase/?</footnote>
curves can also be downloaded and processed using PYCHEOPS (Maxted et al. 2022) and the file keys given in Table 1. All of the CHEOPS data were reduced with version 13.1.0 of the data-reduction pipeline (DRP; Hoyer et al. 2020).

## 2.3. New TESS observations

In addition to the new CHEOPS observations, we made use of previously unpublished TESS data for HD 97658 b from Sector 49 (20 s cadence). In all cases in which we made use of data produced by the TESS Science Processing Operations Center (SPOC) at NASA Ames Research Center (Jenkins et al. 2016), except for HD97658b where we use the Simple Aperture Photometry (SAP) flux (Twicken et al. 2010; Morris et al. 2020), we made use of the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP) flux for the light-curve analysis (Stumpe et al. 2012, 2014; Smith et al. 2012). We chose to use the SAP flux for HD 97658 b, because one more transit is contained in this data set, as compared to the PDCSAP flux, where the last transit is also on a slope. All light curves are available at the MAST<fnref n="5" /> portal.

## 2.4. Previously published observations

### 2.4.1. KELT-9b

For the analysis of KELT-9 b, we used the light curves for this target from Gaudí et al. (2017), who supplied them to us in a private communication<fnref n="6" />. In total, there are 23 primary transits and 7 secondary eclipses observed from 2014 to 2015 using different ground-based telescopes and filters; further details can be found in the filenames of the light curves at the ExoFOP website. There are TESS data from Sectors 14, 15, and 41 available for this target at the MAST portal with a cadence of 120 s, which were also used in the study of Ivshina & Winn (2022). In addition to the transit described in Sect. 2.2, several observations of this target were taken with the CHEOPS space telescope, with a total of nine observed occultations and four phase curves. These were observed from July to September 2020 and July to August 2021. In more detail, transits were observed on 1, 2, and 11 September, 2020, and 31 July and 1, 22, and 24 August, 2021 in the phase curves with CHEOPS (Jones et al. 2022). Furthermore, there is one phase curve including one transit and two secondary eclipses available from Spitzer's Infrared Array Camera (IRAC) at 4.5 µm, which were originally published by Mansfield et al. (2020) and subsequently re-reduced by Jones et al. (2022). We use this latter version of the light curve with no further detrending. We homogeneously re-fitted every single one of these light curves.

### 2.4.2. KELT-16b

For KELT-16 b, we used 19 transit light curves from Oberst et al. (2017) which were obtained using the KELT-North Follow-up Network (KELT-FUN) between May and December 2015. These light curves were provided to us by the authors in a private communication. The observations were made from ten different member observatories of KELT-FUN using different filter sets; details can be found in Table 2 of Oberst et al. (2017). Some of those transits were observed simultaneously using different telescopes or different filters. There are two transits that

were observed with two telescopes or filters each, two that were observed with four telescopes or filters each, and seven that were observed a single time. We also used 11 publicly available transit light curves from Maciejewski et al. (2018). These were acquired between November 2016 and October 2018 using the telescopes and filters stated in Table 1 of their publication. Two more light curves were supplied to us by Patra et al. (2020) in a private communication. These were observed on 10 and 11 June, 2017, with the 1.2 m telescope at the Fred Lawrence Whipple Observatory located in Arizona using images from the KeplerCam detector and a Sloan r'-band filter. We also analysed the public<fnref n="7" /> light curves from the second data release of the ExoClock project (Kokori et al. 2022). This data release contains 32 KELT-16 b transit light curves observed with different telescopes and filters, acquired between July 2018 and November 2020. Details regarding the observation setups can be found at the project homepage under 'ExoClock Observations'. Additionally, we used the 36 public transit light curves from Mancini et al. (2022). These light curves were obtained between June 2016 and June 2021, using the telescopes and filters stated in Table 1 of their publication. In total, 30 planetary transits were observed, with three of them having been observed with two telescopes, and one with three telescopes. TESS data for this target were also obtained from the MAST, including data from Sectors 15 and 41 obtained with a cadence of 120 s and already used in earlier publications. For this target, we also re-fitted all light curves homogeneously.

### 2.4.3. WASP-4b

For WASP-4 b, we made use of the mid-transit times from the homogeneously re-analysed light curves from Baluev et al. (2020). In their re-analysis, the authors process a total of 124 light curves, including those from an earlier paper (Baluev et al. 2019), where they used data from Wilson et al. (2008), Gillon et al. (2009), Sanchis-Ojeda et al. (2011), Nikolov et al. (2012), Petrucci et al. (2013), and amateur observations. Additionally, they used transit light curves from Southworth et al. (2009) which were re-assessed by Southworth et al. (2019) because of potential clock errors, together with new observations. Moreover, Southworth et al. (2019) also obtained previously non-public light curves from Hoyer et al. (2013) and Huitson et al. (2017), with the latter from transmission spectroscopy with GEMINI, offering high precision. Details about the observations can be found in Sect. 2 of the respective publications. The then available TESS Sector 2 light curve was used by these latter authors as well, and also six new amateur observations.

In our analysis, we used the available TESS light curves for this target from Sectors 2 (120 s cadence), 28, and 29 (both 20 s cadence) obtained from the MAST portal, and re-analysed them ourselves. Also, we made use of the publicly available light curves from the ExoClock project again, and re-fit these as well. Details about the observing setups of these can be found on the project homepage. In addition, we obtained the original WASP light curves used in Wilson et al. (2008) from 2006, and also those recorded in the years 2007, 2010, and 2011, and re-fitted them. In addition to the transits, we used four secondary eclipse timings from the literature (Cáceres et al. 2011; Beerer et al. 2011; Zhou et al. 2015), which were observed with ESO's VLT with the Infrared Spectrometer And Array Camera in the *Ks* band, two with the *Spitzer* Space Telescope's IRAC camera (warm Spitzer) at 3.6 µm and 4.5 µm, and one with the

<footnote n="5">⁵ https://mast.stsci.edu/</footnote>

<footnote n="6">⁶ These light curves are now also available at the TESS ExoFOP website.</footnote>

<footnote n="7">⁷ The ExoClock data can be accessed via the project homepage or via at osf.io.</footnote>
Anglo-Australian telescope using the IRIS2 instrument in the Ks-band.

### 2.4.4. HD 97658 b

For the analysis of the HD 97658 system, as Maxted et al. (2022), we also used the published mid-transit time of one Spitzer light curve at 4.5 µm from Van Grootel et al. (2014), and 18 mid-transit times from Guo et al. (2020) using HST/WFC3 spectroscopy, STIS on HST, the Spitzer Space Telescope at 3.6 µm and 4.5 µm, and the MOST Space Telescope in its 0.5 µm band-pass. More details can be found in Sects. 2 and 3 in Guo et al. (2020). Additionally, we also used the single mid-transit time from TESS given in Maxted et al. (2022), but not the time of mid-transit for the then only available CHEOPS light curve from April 2020. This was done to homogeneously re-analyse this single mid-transit time from CHEOPS together with three further light curves recorded by CHEOPS, two of which were observed in March 2021, with the remaining one having been observed in January 2022.

## 3. Transit and secondary eclipse fitting

For the analysis of photometric CHEOPS data, we first use the PYTHON package PYCHEOPS (Maxted et al. 2022). PYCHEOPS is a publicly available<fnref n="8" /> PYTHON module for the analysis of data from the ESA CHEOPS mission. This package can help to deal with the systematic effects present for the photometric data of this space telescope, which are induced by its nadir-locked orbit (Maxted et al. 2022). Moreover, PYCHEOPS offers the ability to detrend the data of other effects, such as correlated noise using Gaussian process regression, fit transit and eclipse models to the light curves to retrieve certain parameters, and more. Fitting can be done using least-squares minimisation and also Markov chain Monte Carlo algorithms. More detailed information can be found in Maxted et al. (2022). In addition to the above and data visualisation, this module also offers a built-in client for data handling, access, and retrieval from the DACE website<fnref n="9" /> (Data & Analysis Center for Exoplanets) hosted by the University of Geneva. General users can access the public data from CHEOPS, with CHEOPS science team members also having the possibility to access proprietary data.

We use this to download the data using the respective file keys. In all cases, we use the 'OPTIMAL' aperture and also use the 'decontaminate' option, which performs a subtraction of the contamination from nearby sources. We plot the light curve using MATPLOTLIB (Hunter 2007) and perform a 2σ outlier clipping using the 'clip_outliers' function from PYCHEOPS. In the case of phase curves, we cut out the individual transits or occultations and separate them. Afterwards, we make use of the PYCHEOPS 'flatten' function with the centre value of the mask being the centre of the transit or occultation from the plot. The mask width value is chosen to fit the transit or occultation width. Next, we save the individual transit or occultation light curves, consisting of time, flux, flux errors, and roll angle values, into two separate files, one including the roll angle information and the other one not. The data handling is carried out using NUMPY (Harris et al. 2020). These two data files will be read in from the Transit and Light Curve Modeller (TLCM, version 97) (Csizmadia 2020). TLCM is a free software tool<fnref n="10" /> that can be

used to analyse, fit, and simulate light curves and radial velocity curves of transiting exoplanets and detached eclipsing binaries. The Mandel & Agol (2002) model is used for the description of transits and occultations, and is enhanced by also accounting for beaming, gravity darkening (see Sect. 4.2.1), reflection, and ellipsoidal effects of both the star and the planet. It also features a wavelet model for red noise based on the model of Carter & Winn (2009). TLCM has the ability to fit photometric and radial velocity data simultaneously, including a simplified Rossiter-McLaughlin effect. A genetic algorithm is used to find the global minimum of the $\chi^2$ or log L values. Subsequently, a simulated annealing algorithm is used to refine the fit, with a Markov chain Monte Carlo algorithm being used for error estimation.

The light curve files are read in by TLCM by entering their names into the configuration file and into a script. The roll angle decorrelation is handled using six parameters to fit the roll angle effect to the data set using the following formula:

$$f_{RA} = p_1 \sin(1 \cdot RA) + p_2 \sin(2 \cdot RA) + p_3 \sin(3 \cdot RA) + p_4 \cos(1 \cdot RA) + p_5 \cos(2 \cdot RA) + p_6 \cos(3 \cdot RA), \quad (1)$$

with $f_{RA}$ being the value needed to add to the corresponding flux value to obtain the decorrelated flux, $p_1$ to $p_6$ being the six fit parameters, and RA being the respective roll angle value in radians. We note that the roll angle decorrelation could also be done in PYCHEOPS.

The necessary fit parameters are the semi-major axis of the planetary orbit, planet-to-star radius ratio, impact parameter, limb-darkening parameters, epoch of the transit, and the orbital period if multiple transits are fitted at the same time. These values, except for the limb-darkening parameters, were taken from TEPCAT (Southworth 2011) as starting points for each fit, except for the transit epoch, which was read off of the transit plots for each transit. After running TLCM, we extracted the transit epochs, as well as their error bars from the results files using the median solution values. In the cases where we created combined models using several transits, for example from TESS, we did a TLCM run while fitting for the parameters as stated above using the epoch of one of the transits near the centre of the data set. After completion of the run, we use the resulting parameter values for the individual transits and fix them, except for the epoch, which is determined visually for each transit as before. As in the case of individual transit fits without a combined model, the resulting epochs are then used in the TTV fits of our three models.

For the occultation fitting, we use the same procedure as before to cut out the occultations from the light curves. Afterwards, we phase-fold the eclipses from one observing season – with the exception of Spitzer data – using the PYASTRONOMY package (Czesla et al. 2019) and fit them using a BATMAN (Kreidberg 2015) box model together with LMFIT (Newville et al. 2014) and EMCEE (Foreman-Mackey et al. 2013) for the Markov chain Monte Carlo (MCMC) fitting. Our box model uses an epoch near the centre of the data set, together with the period, planet radius, semi-major axis, inclination, eccentricity, and uniform limb darkening. The values for these parameters, as before, were taken from TEPCAT. First, we do a preliminary fit to get the combined shape of all occultations in the data set – by leaving all parameters free with reasonable error bars – and the best-fit period from the linear model of the TTV fit for each target. After this, we fix the eclipse shape to the resulting values of the fit and leave only the epoch as a free parameter. Using 5000 steps, 500 steps of burn-in and 200 walkers for the MCMC algorithm of

<footnote n="8">⁸ Available at GitHub.</footnote>

<footnote n="9">⁹ https://dace.unige.ch/</footnote>

<footnote n="10">¹⁰ It is freely available at transits.hu together with a user manual.</footnote>
<table><caption>**Table 2.** Tidal factor *f* in the expression of the period derivative (8) as a function of the orbital obliquity *ψ*, that is the angle between the orbital and stellar rotational plane, measured from the stellar equator, for an equatorial orbit EO (*ψ* = 0), and for a polar orbit PO (*ψ* = 90°).</caption><thead><tr><th>Case</th><th>General expression</th><th>EO</th><th>PO</th></tr></thead><tbody><tr><td>P &lt; P<sub>*</sub></td><td>-<sup>27</sup>&frasl;<sub>2</sub></td><td>-<sup>27</sup>&frasl;<sub>2</sub></td><td>-<sup>27</sup>&frasl;<sub>2</sub></td></tr><tr><td>P = P<sub>*</sub></td><td>-<sup>27</sup>&frasl;<sub>4</sub> + <sup>27</sup>&frasl;<sub>4</sub> cos(ψ) - <sup>27</sup>&frasl;<sub>4</sub> sin<sup>2</sup>(ψ) - <sup>27</sup>&frasl;<sub>8</sub> cos(ψ) sin<sup>2</sup>(ψ) + <sup>27</sup>&frasl;<sub>32</sub> sin<sup>4</sup>(ψ)</td><td>0</td><td>-<sup>405</sup>&frasl;<sub>32</sub></td></tr><tr><td>P<sub>*</sub> &lt; P &lt; 2P<sub>*</sub></td><td><sup>27</sup>&frasl;<sub>2</sub> cos(ψ) - <sup>27</sup>&frasl;<sub>2</sub> sin<sup>2</sup>(ψ) - <sup>27</sup>&frasl;<sub>4</sub> cos(ψ) sin<sup>2</sup>(ψ) + <sup>27</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>+<sup>27</sup>&frasl;<sub>2</sub></td><td>-<sup>189</sup>&frasl;<sub>16</sub></td></tr><tr><td>P = 2P<sub>*</sub></td><td><sup>27</sup>&frasl;<sub>2</sub> cos(ψ) - <sup>27</sup>&frasl;<sub>4</sub> sin<sup>2</sup>(ψ) - <sup>27</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>+<sup>27</sup>&frasl;<sub>2</sub></td><td>-<sup>135</sup>&frasl;<sub>16</sub></td></tr><tr><td>P &gt; 2P<sub>*</sub></td><td><sup>27</sup>&frasl;<sub>2</sub> cos(ψ) + <sup>27</sup>&frasl;<sub>4</sub> cos(ψ) sin<sup>2</sup>(ψ) - <sup>81</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>+<sup>27</sup>&frasl;<sub>2</sub></td><td>-<sup>81</sup>&frasl;<sub>16</sub></td></tr></tbody></table>

EMCEE, we fit for the epoch and its error. The resulting value is then used in the corresponding TTV fit.

# 4. Timing analysis and results

## 4.1. Timing analysis models

### 4.1.1. Keplerian orbit model

To analyse the timing data, we used three models, following Patra et al. (2017). The first model assumes a Keplerian system with a circular orbit, a constant orbital period and hence a linear ephemeris:

$$t_{\text{tra}}(N) = t_0 + N P, \qquad (2)$$

$$t_{\text{occ}}(N) = t_0 + \frac{P}{2} + N P, \qquad (3)$$

with $t_{\text{tra}}$ and $t_{\text{occ}}$ the calculated mid-transit and mid-occultation times, $t_0$ the reference mid-transit time, $N$ the number of orbits from the reference mid-transit time, and $P$ the orbital period of the exoplanet. The parameters $t_0$ and $P$ are fitted in this model.

### 4.1.2. Orbital decay model

The second model is quadratic and assumes a circular orbit with a constant change in the orbital period due to angular momentum transfer from the planet to the star (see e.g. Counselman 1973; Rasio et al. 1996):

$$t_{\text{tra}}(N) = t_0 + N P + \frac{1}{2} \frac{dP}{dN} N^2, \qquad (4)$$

$$t_{\text{occ}}(N) = t_0 + \frac{P}{2} + N P + \frac{1}{2} \frac{dP}{dN} N^2. \qquad (5)$$

In this model, we fit for $t_0$, $P$, and the decay rate $\frac{dP}{dN}$ which is related to the period derivative $\dot{P}$ by $\dot{P} = \frac{dP}{dt} = \frac{1}{P} \frac{dP}{dN}$.

However, if the planet is on an inclined orbit around its host star, this can have an effect on the rate of the tidal decay. According to Kaula's theory of tide (see e.g. Boué & Efroimsky 2019, Eq. (138)), the evolution of the semimajor axis of a circular orbit is, at the quadrupole order $l = 2$,

$$\frac{da}{dt} = -2 a n \frac{M_p}{M_*} \left(\frac{R_*}{a}\right)^5 \sum_{m=0}^{l} \frac{(l-m)!}{(l+m)!} (2 - \delta_{0m}) \cdot \sum_{p=0}^{l} (l-2p) F_{\text{lmpl}}^2(\psi) K_2(\omega_{\text{lmpl}}), \qquad (6)$$

with $\omega_{\text{lmpl}} \approx (l - 2p + q)n - m\dot{\theta}$, where *n* is the orbital mean motion and $\dot{\theta}$ the angular speed of the star. In the semimajor axis variation rate (Eq. (6)), $F_{\text{lmpl}}(\psi)$ are the inclination functions defined by Kaula (1964) and $K_2$ is the quality function of the star. Under the constant-phase-lag model, the expression of $K_2$ in terms of the modified quality factor $Q'_*$ is

$$K_2(\omega) = \frac{3}{Q'_*} \text{Sign}(\omega). \qquad (7)$$

All the tidal frequencies $\omega_{\text{lmpl}}$ have a constant sign except $\omega_{2100} = 2n - \dot{\theta}$ and $\omega_{2200} = 2n - 2\dot{\theta}$. Therefore, in the calculation of $\dot{a}$ (Eq. (6)), five cases have to be considered, namely, $\dot{\theta} < n$, $\dot{\theta} = n$, $n < \dot{\theta} < 2n$, $\dot{\theta} = 2n$, and $\dot{\theta} > 2n$. In terms of the periods $P = 2\pi/n$ and $P_* = 2\pi/\dot{\theta}$, the expression of the orbital period derivative $\dot{P}/P = \frac{3}{2}\dot{a}/a$ is

$$\dot{P} = f \frac{\pi}{Q'_*} \left( \frac{M_p}{M_*} \right) \left( \frac{R_*}{a} \right)^5, \qquad (8)$$

where *f* is a numerical factor given in Table 2. In particular, for a short-period planet, $P < P_*$, we retrieve the formula of Goldreich & Soter (1966):

$$(\dot{P})_{P < P_*} = -\frac{27 \pi}{2 Q'_*} \left(\frac{M_p}{M_*}\right) \left(\frac{R_*}{a}\right)^5. \qquad (9)$$

This enables us to estimate the modified tidal quality factor $Q'_*$ from the measured rate of orbital decay. This rate can also provide us with a lower limit on $Q'_*$, which can be obtained by calculating $Q'_*$ using the 95% confidence lower limit of $\dot{P}$, with the errors being derived from propagating the uncertainties in $M_p/M_*$ and $R_*/a$, as has been done in for example Patra et al. (2020) and Mancini et al. (2022).

### 4.1.3. Apsidal precession model

The third model assumes a non-zero eccentricity *e* and apsidal precession, following the formulations of Giménez & Bastero (1995):

$$t_{\text{tra}}(N) = t_0 + N P_s - \frac{e P_a}{\pi} \cos \omega(N), \qquad (10)$$

$$t_{\text{occ}}(N) = t_0 + \frac{P_a}{2} + N P_s + \frac{e P_a}{\pi} \cos \omega(N), \qquad (11)$$

with $P_s$ the sidereal period, $P_a$ the anomalistic period, and $\omega$ the argument of pericentre. The dependency of $\omega$ on $N$ is as follows:

$$\omega(N) = \omega_0 + \frac{d\omega}{dN} N, \qquad (12)$$
<table><caption>**Table 3.** Tidal factor $f_e$ in the expression of the eccentricity derivative (Eq. (14)) as a function of the orbital obliquity $\psi$, the angle between the orbital and stellar rotational plane, measured from the stellar equator, for an equatorial orbit EO ($\psi = 0$), and for a polar orbit PO ($\psi = 90^\circ$).</caption><thead><tr><th>Case</th><th>General expression</th><th>EO</th><th>PO</th></tr></thead><tbody><tr><td>P &lt; 1/2 P<sub>★</sub></td><td>+ <sup>25</sup>&frasl;<sub>4</sub></td><td>+ <sup>25</sup>&frasl;<sub>4</sub></td><td>+ <sup>25</sup>&frasl;<sub>4</sub></td></tr><tr><td>P = 1/2 P<sub>★</sub></td><td>+ <sup>101</sup>&frasl;<sub>16</sub> + <sup>1</sup>&frasl;<sub>16</sub> cos(ψ) - <sup>1</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) - <sup>1</sup>&frasl;<sub>32</sub> sin<sup>2</sup>(ψ) cos(ψ) - <sup>35</sup>&frasl;<sub>128</sub> sin<sup>4</sup>(ψ)</td><td>+ <sup>51</sup>&frasl;<sub>8</sub></td><td>+ <sup>765</sup>&frasl;<sub>128</sub></td></tr><tr><td>1/2 P<sub>★</sub> &lt; P &lt; P<sub>★</sub></td><td>+ <sup>51</sup>&frasl;<sub>8</sub> + <sup>1</sup>&frasl;<sub>8</sub> cos(ψ) - <sup>1</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) - <sup>1</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) cos(ψ) - <sup>35</sup>&frasl;<sub>64</sub> sin<sup>4</sup>(ψ)</td><td>+ <sup>13</sup>&frasl;<sub>2</sub></td><td>+ <sup>365</sup>&frasl;<sub>64</sub></td></tr><tr><td>P = P<sub>★</sub></td><td>+ <sup>53</sup>&frasl;<sub>8</sub> + <sup>3</sup>&frasl;<sub>8</sub> cos(ψ) - <sup>23</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) - <sup>1</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>37</sup>&frasl;<sub>64</sub> sin<sup>4</sup>(ψ)</td><td>+7</td><td>+ <sup>369</sup>&frasl;<sub>64</sub></td></tr><tr><td>P<sub>★</sub> &lt; P &lt; 3/2 P<sub>★</sub></td><td>+ <sup>55</sup>&frasl;<sub>8</sub> + <sup>5</sup>&frasl;<sub>8</sub> cos(ψ) - <sup>11</sup>&frasl;<sub>4</sub> sin<sup>2</sup>(ψ) - <sup>3</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>109</sup>&frasl;<sub>64</sub> sin<sup>4</sup>(ψ)</td><td>+ <sup>15</sup>&frasl;<sub>2</sub></td><td>+ <sup>373</sup>&frasl;<sub>64</sub></td></tr><tr><td>P = 3/2 P<sub>★</sub></td><td>+ <sup>61</sup>&frasl;<sub>16</sub> - <sup>39</sup>&frasl;<sub>16</sub> cos(ψ) + <sup>5</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) + <sup>43</sup>&frasl;<sub>32</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>169</sup>&frasl;<sub>128</sub> sin<sup>4</sup>(ψ)</td><td>+ <sup>11</sup>&frasl;<sub>8</sub></td><td>+ <sup>697</sup>&frasl;<sub>128</sub></td></tr><tr><td>3/2 P<sub>★</sub> &lt; P &lt; 2P<sub>★</sub></td><td>+ <sup>3</sup>&frasl;<sub>4</sub> - <sup>11</sup>&frasl;<sub>2</sub> cos(ψ) + <sup>27</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) + <sup>23</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>15</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>- <sup>19</sup>&frasl;<sub>4</sub></td><td>+ <sup>81</sup>&frasl;<sub>16</sub></td></tr><tr><td>P = 2P<sub>★</sub></td><td>+ <sup>3</sup>&frasl;<sub>4</sub> - <sup>11</sup>&frasl;<sub>2</sub> cos(ψ) + <sup>29</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) + <sup>25</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>13</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>- <sup>19</sup>&frasl;<sub>4</sub></td><td>+ <sup>83</sup>&frasl;<sub>16</sub></td></tr><tr><td>2P<sub>★</sub> &lt; P &lt; 3P<sub>★</sub></td><td>+ <sup>3</sup>&frasl;<sub>4</sub> - <sup>11</sup>&frasl;<sub>2</sub> cos(ψ) + <sup>31</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) + <sup>27</sup>&frasl;<sub>8</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>11</sup>&frasl;<sub>16</sub> sin<sup>4</sup>(ψ)</td><td>- <sup>19</sup>&frasl;<sub>4</sub></td><td>+ <sup>85</sup>&frasl;<sub>16</sub></td></tr><tr><td>P = 3P<sub>★</sub></td><td>+ <sup>3</sup>&frasl;<sub>4</sub> - <sup>11</sup>&frasl;<sub>2</sub> cos(ψ) + <sup>13</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) + <sup>5</sup>&frasl;<sub>16</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>71</sup>&frasl;<sub>32</sub> sin<sup>4</sup>(ψ)</td><td>- <sup>19</sup>&frasl;<sub>4</sub></td><td>+ <sup>121</sup>&frasl;<sub>32</sub></td></tr><tr><td>P &gt; 3P<sub>★</sub></td><td>+ <sup>3</sup>&frasl;<sub>4</sub> - <sup>11</sup>&frasl;<sub>2</sub> cos(ψ) - <sup>9</sup>&frasl;<sub>4</sub> sin<sup>2</sup>(ψ) - <sup>11</sup>&frasl;<sub>4</sub> sin<sup>2</sup>(ψ) cos(ψ) + <sup>15</sup>&frasl;<sub>4</sub> sin<sup>4</sup>(ψ)</td><td>- <sup>19</sup>&frasl;<sub>4</sub></td><td>+ <sup>9</sup>&frasl;<sub>4</sub></td></tr></tbody></table>

with $\omega_0$ the argument of pericentre at the reference time. The relation of the sidereal period and the anomalistic period can be described as:

$$P_s = P_a \left(1 - \frac{1}{2\pi} \frac{d\omega}{dN}\right). \qquad (13)$$

In this model, we fit for the parameters $t_0$, $P_s$, $e$, $\omega_0$, and $\frac{d\omega}{dN}$. Should there be an eccentricity in the system, we would expect it to be damped by the tidal forces. Both the tides raised on the star and on the synchronous planet tend to circularize the orbit. From Kaula's theory of tides (e.g. Boué & Efroimsky 2019, Eq. (154)), the rate of eccentricity damping, computed at first order in eccentricity, can be put in a form similar to that in Yoder & Peale (1981), namely

$$\frac{de}{dt} = -\frac{1}{3} c (7D + f_e) e, \qquad (14)$$

where the constants $c$ and $D$ are given by

$$c = \frac{27}{2Q_{\star}} \frac{M_p}{M_{\star}} \left(\frac{R_{\star}}{a}\right)^5 n \quad \text{and} \quad D = \frac{k_2}{3} \left(\frac{R_p}{R_{\star}}\right)^5 \left(\frac{M_{\star}}{M_p}\right)^2 \frac{Q'_{\star}}{Q_p}. \qquad (15)$$

In these formulae, $k_2$ is the planet's second Love number, $Q_p$ its tidal quality factor, and $R_p$ its radius. The parameter $f_e$, provided in Table 3, represents the contribution of the tides raised on the star to the circularization.

## 4.2. Results

### 4.2.1. KELT-9 b

**Transit fitting.** For the transit light curves publicly available from Gaudí et al. (2017), we fitted each transit individually with TLCM. The reason for this is that the observations were made using different telescopes and filters, making it too computationally expensive to create a combined model from all light curves at once and achieve convergence, because of all of the extra parameters associated with the implementation of a different noise model per instrument and also different limb-darkening parameters per filter. On the other hand, for the TESS

and CHEOPS data, a combined model was created by jointly fitting both data sets with TLCM to improve the transit fits especially for the CHEOPS observations, because some of those lack parts of the ingress or egress, leading to uncertain fits. This is countered by fixing the transit shape from the combined model in the individual transit fits. We do not use priors for these, because our approach leads to the same mid-transit times as when using priors (within 1.5 s), but gives less strict error bars by about 15% to 20%. For this system, CHEOPS and TESS have relatively close theoretical limb-darkening parameters with $a_{TESS} = 0.1712$, $b_{TESS} = 0.2399$, $a_{CHEOPS} = 0.2571$, and $b_{CHEOPS} = 0.3264$ (Claret 2018, 2021). The validity of the combination of the two data sets to create a combined model is given by the comparison of the resulting transit times of the TESS and CHEOPS combined model with those of the combined TESS model from all TESS transit observations for this target, with a maximum difference in the mid-transit times of 3.6 s and most of them being smaller than 2 s. An additional feature of the transit light curves for this target is gravity darkening, which can distort the transit shape, but account for this requires extra modelling (see e.g. Ahlers et al. (2020)). However, gravity darkening is implemented in TLCM, and so we accounted for this in the creation of the combined model from CHEOPS and TESS and fixed it afterwards when fitting for the individual mid-transit times.

TLCM parameterises the gravity darkening by fitting for two angles, the inclination of the stellar rotation axis, and the angle between it and celestial north, $\Omega_*$ which is related to the sky-projected obliquity, $\lambda$ by $\Omega_* = 90^\circ - \lambda$ (Lendl et al. 2020). The gravity-darkening coefficient, $\beta$, is calculated from the relation between stellar oblateness and $\beta$ proposed by Espinosa Lara & Rieutord (2011); for KELT-9, $\beta = 0.22$. The previously unpublished transit observed with CHEOPS is shown in Fig. 1 together with the already published transits, and the resulting mid-transit times are given in Table 4.

**Secondary eclipse fitting.** The occultation observations for this target yield a total of eight occultation times: one from the ground-based observations of Gaudí et al. (2017) with a relatively large error bar, two from Spitzer, three from TESS, and two from CHEOPS. The nature of the Spitzer data for this
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Transits of KELT-9 b observed with CHEOPS. The x-axis shows the orbital phase, with the y-axis showing the normalised flux from the star. There is an offset of 0.007 in terms of normalised flux between the individual transits. The black dots represent the data, which were corrected for the roll angle of the satellite using PYCHEOPS and were fitted with TLCM, leading to the transit models (red lines). The lowest transit resulted from the single transit observation listed in Table 1; the rest were observed in phase curves (Jones et al. 2022).</figcaption></figure>

<table><caption>**Table 4.** KELT-9 b mid-transit times and errors.</caption><thead><tr><th>Mid-transit time<br>(BJD<sub>TDB</sub>)</th><th>Error<br>(days)</th><th>Epoch</th><th>Source</th></tr></thead><tbody><tr><td>2456873.51672</td><td>0.00119</td><td>-870</td><td>Gaudi et al. (2017)</td></tr><tr><td>2456876.47963</td><td>0.00117</td><td>-868</td><td>Gaudi et al. (2017)</td></tr><tr><td>2456876.47981</td><td>0.00062</td><td>-868</td><td>Gaudi et al. (2017)</td></tr><tr><td>2456895.73193</td><td>0.00134</td><td>-855</td><td>Gaudi et al. (2017)</td></tr><tr><td>2456895.73416</td><td>0.00151</td><td>-855</td><td>Gaudi et al. (2017)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** Epoch given from the middle of the data set. The ‘source’ column describes the source of the respective light curve. Only a portion of the data is shown here, the full table is available at the CDS.

target allows the individual secondary eclipses of this target to be fitted with a similar precision to the transit fits, which is not the case for TESS and CHEOPS, here only the combination of multiple occultation observations allows a similar precision. Although the precision of the CHEOPS secondary eclipse light curves would be sufficient for individual occultation fits, the missing data points in between CHEOPS observations – when the target cannot be observed due to the Earth being in the line of sight of the target – make it hard to get precise timings, especially if parts of the ingress or egress are missing.

When mid-occultation times are derived from multiple occultation light curves, the reference time of conjunction is chosen to be as close to the middle of the data set as possible. This is done because the occultation observations are relatively evenly spaced in each data set. The resulting mid-occultation times can be found in Table 5.

<table><caption>**Table 5.** KELT-9 b mid-occultation times and errors.</caption><thead><tr><th>Mid-occultation time<br>(BJD<sub>TDB</sub>)</th><th>Error<br>(days)</th><th>Epoch</th><th>Source</th></tr></thead><tbody><tr><td>2457117.15824</td><td>0.00138</td><td>-705.5</td><td>Gaudi et al. (2017)</td></tr><tr><td>2458414.62225</td><td>0.00031</td><td>170.5</td><td>Mansfield et al. (2020)</td></tr><tr><td>2458416.10255</td><td>0.00028</td><td>171.5</td><td>Mansfield et al. (2020)</td></tr><tr><td>2458694.55498</td><td>0.00044</td><td>359.5</td><td>TESS</td></tr><tr><td>2458725.65775</td><td>0.00061</td><td>380.5</td><td>TESS</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** The epochs are given relative to the mid-time of the data set. The ‘source’ column denotes the source of the respective light curves. Only a portion of the data is shown here, the full table is available at the CDS.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** O–C plot showing the deviations in transit time from the best-fit linear ephemeris (grey dashed line) for KELT-9 b. The transit number is shown on the x-axis. The y-axis shows the difference in observed and calculated mid-transit time. The orange line shows the best orbital decay fit to the KELT-9 b TTV data, with the red line showing the best apsidal precession fit. CHEOPS data are highlighted in black.</figcaption></figure>

**Timing analysis.** In combination, the transit and occultation data yield 76 times of mid-transit or mid-occultation, with the timings from the ground-based observations of Gaudi et al. (2017) having relatively large error bars in comparison to the those of CHEOPS, Spitzer, and TESS. These data points are fitted using the three models described in Sect. 4.1 using EMCEE. First of all, we gather the necessary planetary and stellar parameters from TEPCAT. Then, the three models are fitted to the transit and occultation timing data using an MCMC algorithm with 30 000 steps, a burn-in period of 10 000 steps, and 500 walkers. The data points are weighted according to their timing errors. The resulting fits are shown in Figs. 2 and 3, and the fit parameters for each model are available in Table B.2.

The BIC values indicate the preference of the apsidal precession model with $BIC_p = 90.5$ over the orbital decay ($BIC_d = 95.3$) and the Keplerian orbit model ($BIC_k = 103.7$), taking into account the number of free parameters of each fit. Furthermore, the $d\omega/dN$ value is in agreement with a constant period only at over $5\sigma$, and the resulting orbital period deviates from the other two at around $3\sigma$. The best-fit eccentricity value agrees with a circular orbit only at $4\sigma$. Furthermore, Stephan et al. (2022) recently found a nodal precession trend using radial velocity measurements, which could also have an effect on the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** KELT-9 b occultations from various sources, as stated in the legend. The x-axis shows the epoch of the occultation, and the y-axis shows the difference in observed and calculated mid-occultation time, assuming a linear ephemeris. The best fits to the constant period, orbital decay, and apsidal precession models are shown as the grey dashed, orange solid, and red solid lines.</figcaption></figure>

transit timings, albeit generally a much smaller signal in comparison to apsidal precession for typical hot Jupiters. Transit duration variations are a better indicator for nodal precession (Miralda-Escudé 2002; Ragozzine & Wolf 2009; Damiani & Lanza 2011).

Moreover, the fits to the Spitzer data yield small error bars and are also in agreement with a constant period, which has an influence especially on the apsidal precession fit. The early data do not help in constraining the fit to either of the models since the scatter and also the error bars are large in comparison to the newer data sets. The occultation data from the various sources is also inconclusive in differentiating the three models because the scatter of the data points is too large (see Fig. 3).

Further high-precision observations of this target will be necessary in the future to constrain the models further and to differentiate between them. KELT-9 is observed in TESS Sector 55 in camera 3, from 5 August to 1 September, 2022, which could enable better fits and might replace the apsidal precession model as the best model, as this is heavily influenced by the Spitzer transit data point.

Using the results from the orbital decay fit, we can calculate $\dot{P} = (-24.42 \pm 10.66)$ ms yr$^{-1}$, assuming the observed trend is true and only caused by this effect. Calculating the 95% confidence lower limit on the orbital decay timescale $\tau = \frac{P}{|P|}$ yields $\tau > (2.8 \pm 0.7)$ Myr.

Due to the polar orbit of this planet, we have to pay attention to the tidal factor $f$ from Eq. (8) in Table 2 when determining $Q'_{*}$. The stellar rotation period of KELT-9 is $P_{*} = (18.96 \pm 0.34)$ h (Jones et al. 2022), with the orbital period of KELT-9 b being nearly double the rotation period of the star, with a value of around $P = 1.48$ days (Cauley & Ahlers 2022). This leads to a smaller absolute value of the tidal factor of $f = -135/16$ instead of $-27/2$. With this, we can also give a 95% confidence lower limit on $Q'_{*}$ of $Q'_{*} > (9.0 \pm 3.7) \times 10^{4}$, with the best-fit parameter being $Q'_{*,best} = (1.7 \pm 0.7) \times 10^{5}$. This would place KELT-9 near the lower edge of the theoretical predictions, which range from $10^{5}$ to $10^{8.5}$ in the literature (Meibom & Mathieu 2005; Jackson et al. 2008; Hansen 2010; Husnoo et al. 2012; Penev et al. 2012, 2018; Bonomo et al. 2017; Patel & Penev 2022), if the decay trend is found to be genuine. However, some of these studies assumed $Q_{*}$ to be a universal constant. Still, this value is on a similar level to the tidal quality factor of WASP-12 with $Q'_{*} = 1.75 \times 10^{5}$.

Nevertheless, as the apsidal-precession model gives the best fit for now and we get a non-zero eccentricity ($e = 0.00122$) from

it, we can compute the expected tidal circularisation timescale $\tau_{e} = \frac{e}{|P|}$. Using Eq. (14) with $f_{e} = \frac{83}{16}$ from Table 3 with the corresponding system parameters for KELT-9, Jupiter's value for the second Love number ($k_{2} = 0.565$) from the Juno mission (Durante et al. 2020), an estimate for its tidal quality factor ($Q_{p} = 5 \times 10^{5}$, e.g. Goldreich & Soter (1966); Wu (2005)), and the stellar modified quality factor from our tidal decay analysis yields

$$ \tau_{e} = 0.78 \text{ Myr.} \tag{16} $$

This estimate indicates that the circularisation of the orbit should be completed before the planet is tidally disrupted by the star, should orbital decay also be happening. This is especially the case considering that the semi-major axis would be shrinking during that time as well, leading to an acceleration of the circularisation process. Concerning the origin of the eccentricity, it is possible that it is justified by the migration history of the planet, for example if the planet migrated inwards due to high-eccentricity migration or planet-planet scattering. Depending on the final eccentricity after the migration, it is plausible that the eccentricity could be maintained until now, especially considering the relatively small value provided by the apsidal precession fit. Another possibility is that the eccentricity could be excited by an unseen third body in the system.

### 4.2.2. KELT-16 b

**Transit fitting.** For KELT-16 b, we fitted a total of 151 transits. All of these, except those of CHEOPS, were modelled individually with TLCM, as many different telescopes and filters were used for the ground-based observations. Even though there are many transit observations, they offer less precision in comparison to KELT-9 b. The CHEOPS light curves suffer from gaps in between observations, especially if they are located near the ingress or egress phase of the transit (see Fig. 4), meaning the resulting mid-transit times cannot be as precise as the data would allow (Barros et al. 2013; Borsato et al. 2021). Because of this, a combined model was created and applied to the final transit fitting of all seven CHEOPS light curves. The previously unpublished observed transits from CHEOPS are shown in Fig. 4 and the obtained mid-transit times can be found in Table 6.

**Timing analysis.** Fitting the obtained times of mid-transit using the three models as described in Sect. 4.2.1, we obtain Fig. 5. The fit parameters of the three models can be found in Table B.2.

As we do not have observations of secondary eclipses, and those observed by TESS suffer from a low signal-to-noise ratio due to the faintness of the target, we cannot obtain accurate mid-occultation times. This limits us in examining the differences in tidal orbital decay and apsidal precession. However, the Keplerian orbit model is the favoured model in any case, having the lowest BIC value at $BIC_{k} = 292.9$, with the orbital decay model leading to $BIC_{d} = 297.6$, and the apsidal precession model yielding $BIC_{p} = 311.0$. As with KELT-9, the orbital periods resulting from the fits agree with each other well within $1\sigma$. In addition, the $dP/dN$ of the orbital decay model also agrees with a constant period within $1\sigma$, with the resulting eccentricity $e$ of the apsidal precession model agreeing with a circular orbit at $2\sigma$. The $d\omega/dN$ parameter agrees with a Keplerian orbit at only just over $4\sigma$. Comparing our orbital decay fit results to earlier results from the literature, we find that our finding of $dP/dN = (-2.73 \pm 3.77) \times 10^{-10}$ days per orbit
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** KELT-16 b transits observed with CHEOPS. The *x*-axis shows the orbital phase of the transits and the *y*-axis shows the normalised flux. There is an offset of 0.015 in terms of normalised flux between the individual transits. The data points (black dots) were corrected for the roll angle of the satellite using PYCHEOPS and fitted with TLCM. The resulting transit models are indicated by the red lines.</figcaption></figure>

from our independent analysis agrees well with the values of $dP/dN = (-0.1 \pm 1.4) \times 10^{-9}$ days per orbit from Maciejewski et al. (2018), $dP/dN = (-0.6 \pm 1.4) \times 10^{-9}$ days per orbit from Patra et al. (2020), and $dP/dN = (-3.2 \pm 4.0) \times 10^{-10}$ days per orbit from Mancini et al. (2022). These literature results also all agree with a constant orbital period. In general, it can be said that for this object there is an overall insufficiency in precision to find an effect that is as small as the one we are looking for, leading to scatter in the data and relatively large error bars. A longer baseline of observations is needed, with the current baseline only spanning about 6.5 yr. However, this will be extended with TESS in Sector 55, where KELT-16 will be located in the field of view of camera 2. Nevertheless, we can use our tidal decay fit result to constrain the modified stellar tidal quality factor, assuming that the tidal decay trend is genuine. Our calculated lower limit of $Q'_* > (2.1 \pm 0.9) \times 10^5$ for this is close to that of KELT-9, and with that also close to the actual $Q'_*$ value for WASP-12, meaning that we expect a slower orbital decay for KELT-16 b than for WASP-12b. Moreover, our lower limit value is in agreement with that from Mancini et al. (2022) whose value is $(1.9 \pm 0.8) \times 10^5$. Using our best-fit value of $\dot{P} = (-8.94 \pm 12.35) \text{ ms yr}^{-1}$ to calculate the orbital decay timescale $\tau$, we obtain a 95% confidence lower limit of $\tau > (2.5 \pm 0.9) \text{ Myr}$.

### 4.2.3. WASP-4 b

**Transit fitting.** We have a total of 172 transit timing data points for this target, many of which offer high precision. To test the validity of the timings obtained by Baluev et al. (2020), we compared them against those from the analysis of Southworth et al. (2019), those from their own earlier paper (Baluev et al. 2019), and ours. Particularly, the timings with very high precision were re-analysed and verified. Still, there is the possibility that the errors on the measured flux may be underestimated in

<table><caption>**Table 6.** KELT-16 b mid-transit times and errors.</caption><thead><tr><th>Mid-transit time<br>(BJDTDB)</th><th>Error<br>(days)</th><th>Epoch</th><th>Source</th></tr></thead><tbody><tr><td>2457165.85120</td><td>0.00076</td><td>-1177</td><td>Oberst et al. (2017)</td></tr><tr><td>2457166.82162</td><td>0.00076</td><td>-1176</td><td>Oberst et al. (2017)</td></tr><tr><td>2457166.82490</td><td>0.00114</td><td>-1176</td><td>Oberst et al. (2017)</td></tr><tr><td>2457168.75747</td><td>0.00132</td><td>-1174</td><td>Oberst et al. (2017)</td></tr><tr><td>2457196.86000</td><td>0.00132</td><td>-1145</td><td>Oberst et al. (2017)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** Epochs are relative to the middle of the data set. The ‘Source’ column describes the source of the respective light curve. Only a portion of the data is shown here; the full table is available at the CDS.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** O–C plot showing the deviations in transit time from the best-fit linear ephemeris (grey dashed line) for KELT-16 b. The transit number is shown on the *x*-axis, and the timing difference on the *y*-axis. The orange line shows the best quadratic fit to the KELT-16 b TTV data, with the red line showing the best apsidal precession fit. CHEOPS data are highlighted in black. The data from the literature are described in more detail in Sect. 2.4.2.</figcaption></figure>

some of those cases. However, they were recorded with high-precision ground-based instruments, like those of the GEMINI observatory or ESO’s VLT. All but six of their 29 timings from the light curves of Gillon et al. (2009), Sanchis-Ojeda et al. (2011), and Southworth et al. (2019) agree with our results within $3\sigma$, with a median value of $0.69\sigma$. The ExoClock data are scattered and in some cases still only have relatively small error bars, which could be caused by underestimated uncertainties in the original light curves. The transits observed with CHEOPS for this target are shown in Fig. 6, and the literature and newly fitted transit times can be found in Table 7.

**Timing analysis.** Fitting the obtained timings as described in Sect. 4.2.1 using the MCMC algorithms yields the models and parameters in Figs. 7, 8, and in Table B.2.

From the final fit parameters, we can see that the orbital decay model is the preferred model with the lowest BIC value of $BIC_d = 244.63$. The apsidal precession and Keplerian orbit models lead to higher BIC values of $BIC_p = 252.31$ and $BIC_k = 339.52$, respectively. There is no clear favoured model from the occultation measurements due to the scatter of the data and their error bars. There are also only four secondary eclipse timings in the literature and the eclipses are not visible with TESS, as Turner et al. (2022) also found. This hinders differentiation between the models using only the available occultations.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Transits of WASP-4 b observed with CHEOPS. The x-axis shows the orbital phase and the y-axis the normalised flux. There is an offset of 0.03 in terms of normalised flux between the individual transits. The data points (black dots) were corrected for the roll angle of the satellite using PYCHEOPS and fitted with TLCM.</figcaption></figure>

<table><caption>**Table 7.** WASP-4 b mid-transit times and errors.</caption><thead><tr><th>Mid-transit time (BJD<sub>TDB</sub>) (BJD<sub>TDB</sub>)</th><th>Error (days) (days)</th><th>Epoch</th><th>Source</th></tr></thead><tbody><tr><td>2453960.43148</td><td>0.00182</td><td>-2071</td><td>WASP</td></tr><tr><td>2454361.90048</td><td>0.00188</td><td>-1771</td><td>WASP</td></tr><tr><td>2454396.69616</td><td>0.00008</td><td>-1745</td><td>Baluev et al. (2020)</td></tr><tr><td>2454697.79815</td><td>0.00006</td><td>-1520</td><td>Baluev et al. (2020)</td></tr><tr><td>2454697.79830</td><td>0.00013</td><td>-1520</td><td>Baluev et al. (2020)</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** The epochs are given in relation to the middle of the data set. The ‘Source’ column describes the source of the respective light curve or the source of the timing, which is described in more detail in Sect. 2. Only a portion of the data is shown here; the full table is available at the CDS.

For this target, the best-fit orbital period values agree with each other within $1\sigma$. The best-fit $dP/dN$ value of the orbital decay model is non-zero at a $5\sigma$ level, disagreeing with a linear ephemeris and making this trend significant. Nevertheless, the eccentricity value of the apsidal precession model is also in disagreement with a Keplerian orbit at more than $3\sigma$. This, in combination with the position of the CHEOPS data points in the O–C plot in Fig. 7, indicates that the system is not Keplerian. More observations in the future are needed to help to differentiate between the two preferred models, and TESS could contribute to this with a planned observation of this target during Sector 69 of its second extended mission in camera 2 from August to September 2023.

Assuming that the fitted orbital decay trend is genuine, we find that our derived decay rate of $dP/dN = (-2.62 \pm 0.49) \times 10^{-10}$ days per orbit agrees with the values previously given in

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** O–C plot showing the deviations in transit time from the best-fit linear ephemeris (grey dashed line) for WASP-4 b. The transit number is shown on the x-axis, with the difference in timing being shown on the y-axis. The orange line shows the quadratic fit to the WASP-4 b transit timing variation data, and the red line shows the apsidal precession fit. CHEOPS data are highlighted in black.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Occultation timings for WASP-4 b from the literature, with the three lines representing our three models according to the legend. The x-axis shows the epoch relative to the middle of the data set, and the y-axis shows the difference in measured and calculated (linear) timing.</figcaption></figure>

the literature, which range from $(-1.0 \pm 1.0) \times 10^{10}$ days per orbit in the study of Baluev et al. (2019), which was adjusted to $(-1.7 \pm 0.5) \times 10^{-10}$ days per orbit after re-assessing older data and adding new data (Baluev et al. 2020), to $(-5.4 \pm 0.5) \times 10^{-10}$ days per orbit (Bouma et al. 2019), with other literature values lying between them (Southworth et al. 2019; Bouma et al. 2020; Turner et al. 2022). The most recent value from Maciejewski (2022) with $dP/dN = (-2.1 \pm 0.6) \times 10^{-10}$ is in good agreement with our analysis. Under the same assumption, we can also compare our lower limit on $Q'_*$ of $Q'_* > (4.1 \pm 0.9) \times 10^4$ and best-fit value of $Q'_{*,best} = (5.7 \pm 1.0) \times 10^4$ with the literature values and notice that it agrees with the above-mentioned publications, and that it is significantly lower than the theoretical estimates stated in Sect. 4.2.1. If we calculate the potential orbital decay timescale of WASP-4 b, we get:

$$ \tau = 18.73 \pm 3.48 \text{ Myr}, $$

which leads to a 95% confidence lower limit of $\tau_{low} = (13.7 \pm 1.8) \text{ Myr}$. Comparing this to the WASP-12 system, we notice that the orbital decay timescale is much shorter there with $\tau = 3.25_{-0.21}^{+0.24} \text{ Myr}$, even though the tidal quality factor of WASP-12 is higher with $1.75_{-0.11}^{+0.13} \times 10^5$ as compared to WASP-4 with $(5.6 \pm 1.0) \times 10^4$, if we assume that orbital decay is actually happening. However, this can be explained by the combination of
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** New transit observations of HD 97658 b with CHEOPS. The *x*-axis shows the orbital phase, and the *y*-axis shows the normalised flux. There is an offset of 0.0015 in terms of normalised flux between the individual transits. The data points (black dots) were corrected for the roll angle of the satellite using PYCHEOPS and fitted with TLCM. The resulting transit models are shown with red lines for each transit.</figcaption></figure>

the system parameters of these two systems, namely *M*<sub>P</sub>, *M*<sub>*</sub>, *R*<sub>*</sub> and *a*, and their influence in Eq. (9).

### 4.2.4. HD 97658 b

**Transit fitting.** We fitted the four available CHEOPS transits using TLCM by first creating a combined model and then fitting the transits individually using the parameters from the combined model, which is useful due to the gaps in the light curve. For the TESS Sector 49 data, we fitted the data of the whole sector together using the SAP light curve in order to get a more precise timing than what would be possible using individual fits, and because it was already done this way for the earlier TESS observation. The SAP light curve was chosen over the PDCSAP light curve because the latter excludes one of the transits near the centre of the light curve, and the last transit is on a slope, which it is not in the SAP light curve. The SAP light curve was reasonably flat in each transit cut-out and thus required no further detrending. Nevertheless, we still fitted the transits using the PDCSAP flux as well, which resulted in similar mid-transit times (< 0.12 σ difference) and also error bars in comparison to the individual transit fits using the SAP flux. The remaining mid-transit times were taken from Maxted et al. (2022) to have a direct comparison. The newly observed CHEOPS transits are shown in Fig. 9, the fit to the new TESS observation in Fig. A.1, and all obtained mid-transit times can be found in Table B.1.

**Timing analysis.** The timing analysis for this target was done as described in Sect. 4.2.1, with the exception of an apsidal precession fit. In addition to the new quadratic and linear fits for this target, we also show the fits of Guo et al. (2020) and Maxted et al. (2022) in Fig. 10. Our fit parameters can be found in Table B.2.

Comparing our best-fit results for a linear and a quadratic model to the now available data with the linear fit of Guo et al. (2020) and the linear and quadratic fits of Maxted et al. (2022), we find that the orbital period found in the former publication

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** O–C plot showing the deviations in transit time from the best-fit linear ephemeris from Guo et al. (2020) (grey dashed line) for HD 97658 b to have a direct comparison to the Maxted et al. (2022) best-fit linear and quadratic models (red dashed and green solid lines). Our best-fit models are shown as an orange dashed line for the linear model, and as a purple solid line for the quadratic model. The transit number is shown on the *x*-axis, and the difference in timing is shown on the *y*-axis. CHEOPS data are highlighted in black.</figcaption></figure>

of $P_G = (9.489295 \pm 0.000005)$ days still agrees with our linear best-fit period of $P = (9.48930773 \pm 0.00000149)$ days at just over $3\sigma$. The updated orbital period of Maxted et al. (2022) with a value of $P_M = (9.4893072 \pm 0.0000025)$ days improves the earlier result and agrees with our findings at well within $1\sigma$. Comparing the quadratic fits, which was their preferred model by means of BIC, these latter authors find an orbital period of $P_{M,q} = (9.4892968 \pm 0.0000038)$ days, and period change per orbit of $dP/dN_M = (1.46 \pm 0.48) \times 10^{-7}$ days per orbit, whereas we find $P_q = (9.48929994 \pm 0.00000824)$ days and $dP/dN = (0.42 \pm 0.44) \times 10^{-7}$ using the now available data. The quadratic model shows weak evidence against the linear model with $\Delta\text{BIC} = 58.13 - 55.05 = 3.08$, but also agrees with a linear trend within $1\sigma$, rendering this as inconclusive. The scatter in the early data contributes to this significantly. Looking at the Spitzer transits in Fig. 9 in the publication of Guo et al. (2020), one notice the relatively shallow transit depth of the transits in relation to the noise of the data. Comparing this figure with the resulting uncertainties on the transit timings indicates that the uncertainties of the transit fits could be underestimated. This could explain the scatter of this data set in relation to their error bars and our best-fit models. This suspicion is supported by the fact that Van Grootel et al. (2014) also analysed one Spitzer transit recorded with Spitzer’s IRAC at 4.5 µm, with their analysis yielding approximately double the uncertainty on the timing. However, improvements to the Spitzer data reduction pipeline or improved fitting methods might also lead to more precise transit fits. An important thing to mention is that there are only a few outliers, and that most of the Spitzer mid-transit times analysed by Guo et al. (2020) are in agreement with both of the linear models.

In conclusion, by taking new observations of the transits of HD 97658 b, we find a weaker quadratic trend than Maxted et al. (2022), and so are unable to conclusively distinguish between a linear and a quadratic trend.
## 5. Conclusions

We homogeneously analysed new CHEOPS and TESS photometric data and re-analysed archival photometric data for KELT-9b, KELT-16b, WASP-4b, and HD 97658b using state-of-the-art software, namely TLCM and PYCHEOPS, and developed our own MCMC routines for the analysis of this data. Transit fitting was done using TLCM, occultation fitting using BATMAN in combination with an MCMC approach, and the TTVs for each system were fitted using MCMC algorithms for three different models. The first of these models is a constant-period model assuming a Keplerian orbit, the second is an orbital-decay model accounting for a changing orbital period, and the third is an apsidal-precession model assuming that the orbit is slightly eccentric, leading to a precessing orbit that can mimic orbital decay for short baselines but can be distinguished using secondary eclipses or with long baselines.

For the targets KELT-9b and KELT-16b, we did a homogeneous re-analysis of every available transit from the literature and added new observations to these. For KELT-9b, we find that the timing deviations are best described by an apsidal precession model with an eccentricity of $e = 0.00122 \pm 0.00025$; however, we note that these could already be different with new observations, for example those of TESS. Nevertheless, in general, apsidal precession does not rule out tidal decay and vice versa, and a combination of both is likely. The KELT-16 system is best described with a Keplerian orbit model as the scatter and error bars are both relatively large for the observations of this system. Nevertheless, there will also soon be new TESS observations of this system available, which could favour one of the models.

For WASP-4b, we re-analysed data from a recent publication (Turner et al. 2022), as well as TESS data and added new CHEOPS transit observations to the data set. The remaining data were taken from a recent study where the authors also re-analysed earlier transit data in a homogeneous way (Baluev et al. 2020). We find the orbital decay model to describe the data the best and find a similar significant trend at $5\sigma$ as earlier evaluations of this system (Bouma et al. 2020; Baluev et al. 2020; Turner et al. 2022; Ivshina & Winn 2022), but apsidal precession cannot yet be ruled out. The TESS observations of this target towards the end of 2023 could already rule out one of the two preferred models.

In the case of HD 97658b, we re-analyse the only available CHEOPS transit at the time of the study (Maxted et al. 2022) and combine the analysis with three recently acquired transits of this target from CHEOPS. Moreover, we also analyse new TESS data from Sector 49, including three transits of this target. Adding these new transit timings to the data set of Maxted et al. (2022), we find a weaker quadratic trend than these authors, which only deviates from a linear trend at $1\sigma$, leaving our findings inconclusive.

**Acknowledgements.** We thank the anonymous referee for their helpful comments. The CHEOPS photometric data underlying this article is publicly accessible using the CHEOPS archive browser (https://cheops-archive.astro.unige.ch/archive browser/), at the DACE website in the CHEOPS database, or via PYCHEOPS using the file keys provided in Table 1. The TESS data used in this publication is publicly available at the MAST archive. The timing data for the various planets and the photometric CHEOPS data is available with the online version of this article and at the CDS. CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved.

Their flexibility and willingness to explore new approaches were essential to the success of this mission. JVH acknowledges the support of the DFG priority programme SPP 1992 "Exploring the Diversity of Extrasolar Planets (SM 4862-1)". S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project FOUR ACES. P.M. acknowledges support from STFC research grant number ST/M001040/1. L.M.S. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 'Gaseous or rocky? Unveiling the nature of small worlds'. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-190080). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. Y.A. and M.J.H. acknowledge the support of the Swiss National Fund under grant 200020_172746. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia María de Maeztu-Centro de Astrobiología (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. X.B., S.C., D.G., M.F. and J.L. acknowledge their role as ESA-appointed CHEOPS science team members. A.Br. was supported by the SNSA. A.C.C. acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L.D. is an F.R.S.-FNRS Postdoctoral Researcher. This work was supported by FCT - Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987, O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. M.F. and C.M.P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). DG gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 "Gaseous or rocky? Unveiling the nature of small worlds". M.G. is an F.R.S.-FNRS Senior Research Associate. SH gratefully acknowledges CNES funding through the grant 837319. K.G.I. is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. M.L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. LBo, GBr, VNa, IPa, GPi, RRa, GSc, VSi, and TZi acknowledge support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). IRI acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme. S.G.S. acknowledge support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). GyMSz acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a a PRODEX Experiment Agreement No. 4000137122, the Lendület LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. V.V.G. is an F.R.S-FNRS Research Associate. N.A.W. acknowledges UKSA grant ST/R004838/1. A.C.C. and T.W. acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA's Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Programme through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products.
References

Ahlers, J. P., Johnson, M. C., Stassun, K. G., et al. 2020, AJ, 160, 4

Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989

Baluev, R. V., Sokov, E. N., Jones, H. R. A., et al. 2019, MNRAS, 490, 1294

Baluev, R. V., Sokov, E. N., Hoyer, S., et al. 2020, MNRAS, 496, L11

Barnes, J. W. 2009, ApJ, 705, 683

Barros, S. C. C., Boué, G., Gibson, N. P., et al. 2013, MNRAS, 430, 3032

Barros, S. C. C., Akinsanmi, B., Boué, G., et al. 2022, A&A, 657, A52

Beerer, I. M., Knutson, H. A., Burrows, A., et al. 2011, ApJ, 727, 23

Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109

Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107

Borsato, L., Piotto, G., Gandolfi, D., et al. 2021, MNRAS, 506, 3810

Boué, G., & Efroimsky, M. 2019, Celest. Mech. Dyn. Astron., 131, 30

Bouma, L. G., Winn, J. N., Baxter, C., et al. 2019, AJ, 157, 217

Bouma, L. G., Winn, J. N., Howard, A. W., et al. 2020, ApJ, 893, L29

Brandeker, A., Heng, K., Lendl, M., et al. 2022, A&A, 659, A4

Cáceres, C., Ivanov, V. D., Minniti, D., et al. 2011, A&A, 530, A5

Carter, J. A., & Winn, J. N. 2009, ApJ, 704, 51

Cauley, P. W., & Ahlers, J. P. 2022, AJ, 163, 122

Claret, A. 2018, A&A, 618, A20

Claret, A. 2021, RNAAS, 5, 13

Counselman, C. C., I. 1973, ApJ, 180, 307

Csizmadia, S. 2020, MNRAS, 496, 4442

Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, Astrophysics Source Code Library [record ascl:1906.010]

Damiani, C., & Lanza, A. F. 2011, A&A, 535, A116

Deline, A., Hooton, M. J., Lendl, M., et al. 2022, A&A, 659, A74

Dragomir, D., Matthews, J. M., Howard, A. W., et al. 2012, ApJ, 759, L41

Dragomir, D., Matthews, J. M., Eastman, J. D., et al. 2013, ApJ, 772, L2

Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572

Ellis, T. G., Boyajian, T., von Braun, K., et al. 2021, AJ, 162, 118

Lara, E., F., & Rieutord, M. 2011, A&A, 533, A43

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaudí, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514

Gillon, M., Smalley, B., Hebb, L., et al. 2009, A&A, 496, 259

Giménez, A., & Bastero, M. 1995, Ap&SS, 226, 99

Goldreich, P., & Soter, S. 1966, Icarus, 5, 375

Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39

Guo, X., Crossfield, I. J. M., Dragomir, D., et al. 2020, AJ, 159, 239

Hansen, B. M. S. 2010, ApJ, 723, 285

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hebb, L., Collier-Cameron, A., Loeillet, B., et al. 2009, ApJ, 693, 1920

Henry, G. W., Howard, A. W., Marcy, G. W., Fischer, D. A., & Johnson, J. A. 2011, ArXiv e-prints, [arXiv:1109.2549]

Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165

Hooton, M. J., Hoyer, S., Kitzmann, D., et al. 2022, A&A, 658, A75

Howard, A. W., Johnson, J. A., Marcy, G. W., et al. 2011, ApJ, 730, 10

Hoyer, S., López-Morales, M., Rojo, P., et al. 2013, MNRAS, 434, 46

Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24

Huitson, C. M., Désert, J. M., Bean, J. L., et al. 2017, AJ, 154, 95

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Husnoo, N., Pont, F., Mazeh, T., et al. 2012, MNRAS, 422, 3151

Ivshina, E. S., & Winn, J. N. 2022, ApJS, 259, 62

Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE 9913, 99133E

Jones, K. D., Morris, B. M., Demory, B. O., et al. 2022, A&A, 666, A118

Kaula, W. M. 1964, Rev. Geophys. Space Phys., 2, 661

Kokori, A., Tsiaras, A., Edwards, B., et al. 2022, ApJS, 258, 40

Kreidberg, L. 2015, PASP, 127, 1161

Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94

Maciejewski, G. 2022, Acta Astron., 72, 1

Maciejewski, G., Fernández, M., Aceituno, F., et al. 2018, Acta Astron., 68, 371

Mancini, L., Southworth, J., Naponiello, L., et al. 2022, MNRAS, 509, 1447

Mandel, K., & Agol, E. 2002, ApJ, 580, L171

Mansfield, M., Bean, J. L., Stevenson, K. B., et al. 2020, ApJ, 888, L15

Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77

Meibom, S., & Mathieu, R. D. 2005, ApJ, 620, 970

Miralda-Escudé, J. 2002, ApJ, 564, 1019

Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2020, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-003

Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, https://doi.org/10.5281/zenodo.11813

Nikolov, N., Henning, T., Koppenhoefer, J., et al. 2012, A&A, 539, A159

Oberst, T. E., Rodriguez, J. E., Colón, K. D., et al. 2017, AJ, 153, 97

Ogilvie, G. I. 2014, ARA&A, 52, 171

Patel, R., & Penev, K. 2022, MNRAS, 512, 3651

Patra, K. C., Winn, J. N., Holman, M. J., et al. 2017, AJ, 154, 4

Patra, K. C., Winn, J. N., Holman, M. J., et al. 2020, AJ, 159, 150

Penev, K., Jackson, B., Spada, F., & Thom, N. 2012, ApJ, 751, 96

Penev, K., Bouma, L. G., Winn, J. N., & Hartman, J. D. 2018, AJ, 155, 165

Pepper, J., Pogge, R. W., DePoy, D. L., et al. 2007, PASP, 119, 923

Petrucci, R., Jofré, E., Schwartz, M., et al. 2013, ApJ, 779, L23

Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407

Ragozzine, D., & Wolf, A. S. 2009, ApJ, 698, 1778

Rasio, F. A., Tout, C. A., Lubow, S. H., & Livio, M. 1996, ApJ, 470, 1187

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Ser., 9143, 914320

Sanchis-Ojeda, R., Winn, J. N., Holman, M. J., et al. 2011, ApJ, 733, 127

Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000

Southworth, J. 2011, MNRAS, 417, 2166

Southworth, J., Hinse, T. C., Burgdorf, M. J., et al. 2009, MNRAS, 399, 287

Southworth, J., Dominik, M., Jørgensen, U. G., et al. 2019, MNRAS, 490, 4230

Stephan, A. P., Wang, J., Cauley, P. W., et al. 2022, ApJ, 931, 111

Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985

Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100

Turner, J. D., Flagg, L., Ridden-Harper, A., & Jayawardhana, R. 2022, AJ, 163, 281

Twicken, J. D., Clarke, B. D., Bryson, S. T., et al. 2010, Proc. SPIE, 7740, 774023

Van Grootel, V., Gillon, M., Valencia, D., et al. 2014, ApJ, 786, 2

Walker, G., Matthews, J., Kuschnig, R., et al. 2003, PASP, 115, 1023

Wilson, D. M., Gillon, M., Hellier, C., et al. 2008, ApJ, 675, L113

Wong, I., Kitzmann, D., Shporer, A., et al. 2021, AJ, 162, 127

Wu, Y. 2005, ApJ, 635, 688

Yee, S. W., Winn, J. N., Knutson, H. A., et al. 2020, ApJ, 888, L5

Yoder, C. F., & Peale, S. J. 1981, Icarus, 47, 1

Zhou, G., Bayliss, D. D. R., Kedziora-Chudczer, L., et al. 2015, MNRAS, 454, 3002

<footnote n="1">¹ Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany
e-mail: jan-vincent.harre@dlr.de</footnote>

<footnote n="2">² Instituto de Astrofísica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="3">³ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="4">⁴ IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote n="5">⁵ Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="6">⁶ Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden</footnote>

<footnote n="7">⁷ Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland</footnote>

<footnote n="8">⁸ Center for Space and Habitability, University of Bern, Gesellschaftstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="9">⁹ Astrophysics Group, Keele University, Staffordshire, ST5 5BG, UK</footnote>

<footnote n="10">¹⁰ Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="11">¹¹ Physikalisches Institut, University of Bern, Gesellschaftstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="12">¹² Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, I-10125, Torino, Italy</footnote>

<footnote n="13">¹³ Instituto de Astrofísica de Canarias, 38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="14">¹⁴ Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="15">¹⁵ Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK</footnote>

<footnote n="16">¹⁶ Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, 08193 Bellaterra, Spain</footnote>

<footnote n="17">¹⁷ Institut d'Étudis Espaciaux de Catalunya (IEEC), 08034 Barcelona, Spain</footnote>
18 ESTEC, European Space Agency, 2201AZ, Noordwijk, NL

19 Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary

20 Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

21 Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria

22 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

23 Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstr. 1, 85748 Garching, Germany

24 Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

25 Center for Space and Habitability, Gesellschaftstrasse 6, 3012 Bern, Switzerland

26 Université de Paris, Institut de physique du globe de Paris, CNRS, 75005 Paris, France

27 Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK

28 Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

29 Aix Marseille Univ, CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France

30 Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

31 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

32 Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

33 Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden

34 University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria

35 Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

36 Science and Operations Department - Science Division (SCI-SC), Directorate of Science, European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201-AZ Noordwijk, The Netherlands

37 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15-17, Hungary

38 ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

39 Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis blvd. Arago, 75014 Paris, France

40 Department of Astrophysics, University of Vienna, Tuerkenschanzstrasse 17, 1180 Vienna, Austria

41 INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy

42 Institute of Optical Sensor Systems, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

43 Dipartimento di Fisica e Astronomia "Galileo Galilei", Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

44 ETH Zurich, Department of Physics, Wolfgang-Pauli-Strasse 2, 8093 Zurich, Switzerland

45 Zentrum für Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin, Germany

46 Institut für Geologische Wissenschaften, Freie Universität Berlin, 12249 Berlin, Germany

47 MIT Kavli Institute for Astrophysics and Space Research, 70 Vassar St, Cambridge, MA 02139, USA

48 MIT Physics Department, 182 Memorial Dr, Cambridge, MA 02142, USA

49 Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

50 Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

51 Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

52 ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary

53 MTA-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary

54 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

55 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

56 SETI Institute, Mountain View, CA 94043, USA

57 NASA Ames Research Center, Moffett Field, CA 94035, USA
## Appendix A: Additional figure

<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Previously unpublished phase-folded TESS transit data for HD 97658 b (black dots). The red line corresponds to the model resulting from the fits to the three available transits in the SAP flux with TLCM.</figcaption></figure>

## Appendix B: Additional tables

<table><caption>**Table B.1.** HD 97658 b mid-transit times and errors.</caption><thead><tr><th>Mid-transit time<br>[BJD<sub>TDB</sub>]</th><th>Error<br>[d]</th><th>Epoch</th><th>Source</th></tr></thead><tbody><tr><td>6361.80500</td><td>0.00330</td><td>0</td><td>MOST (Guo et al. 2020)</td></tr><tr><td>6523.12544</td><td>0.00061</td><td>17</td><td>Spitzer (Van Grootel et al. 2014)</td></tr><tr><td>6646.48290</td><td>0.00110</td><td>30</td><td>HST/WFC3 (Guo et al. 2020)</td></tr><tr><td>6665.46210</td><td>0.00120</td><td>32</td><td>HST/WFC3 (Guo et al. 2020)</td></tr><tr><td>6712.90960</td><td>0.00240</td><td>37</td><td>MOST (Guo et al. 2020)</td></tr><tr><td>6864.73938</td><td>0.00027</td><td>53</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>6883.71892</td><td>0.00027</td><td>55</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>6893.20742</td><td>0.00026</td><td>56</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7082.99324</td><td>0.00032</td><td>76</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7092.48243</td><td>0.00034</td><td>77</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7101.97164</td><td>0.00034</td><td>78</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7234.82179</td><td>0.00026</td><td>92</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7244.31157</td><td>0.00022</td><td>93</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7253.80001</td><td>0.00045</td><td>94</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7472.05445</td><td>0.00025</td><td>117</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7481.54364</td><td>0.00031</td><td>118</td><td>Spitzer (Guo et al. 2020)</td></tr><tr><td>7491.03120</td><td>0.00110</td><td>119</td><td>HST/WFC3 (Guo et al. 2020)</td></tr><tr><td>7785.20210</td><td>0.00110</td><td>150</td><td>HST/WFC3 (Guo et al. 2020)</td></tr><tr><td>8904.94070</td><td>0.00100</td><td>268</td><td>TESS (Maxted et al. 2022)</td></tr><tr><td>8961.87669</td><td>0.00028</td><td>274</td><td>CHEOPS</td></tr><tr><td>9294.00198</td><td>0.00047</td><td>309</td><td>CHEOPS</td></tr><tr><td>9303.49113</td><td>0.00041</td><td>310</td><td>CHEOPS</td></tr><tr><td>9607.14869</td><td>0.00051</td><td>342</td><td>CHEOPS</td></tr><tr><td>9654.59570</td><td>0.00084</td><td>347</td><td>TESS</td></tr></tbody></table>

**Notes.** The respective epochs are given relative to the beginning of the data set. The ‘Source’ column describes the source of the respective light curve or timing. This is described in more detail in Section 2. The timings without references were obtained by the present authors. This table is available at the CDS.
<table><caption>**Table B.2.** Fit results of all models for KELT-9 b, KELT-16 b, WASP-4 b, and HD 97658 b.</caption><thead><tr><th>Planet</th><th>Model</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="12">KELT-9 b</td><td rowspan="3">Constant period</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458162.09128 ± 0.00009</td></tr><tr><td>P [d]</td><td>1.48111916 ± 0.00000013</td></tr><tr><td>BIC</td><td>103.73</td></tr><tr><td rowspan="4">Orbital decay</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458162.09150 ± 0.00013</td></tr><tr><td>P [d]</td><td>1.48111926 ± 0.00000013</td></tr><tr><td>dP/dN [d/orbit]</td><td>(-11.41 ± 4.98) × 10<sup>-10</sup></td></tr><tr><td>BIC</td><td>95.28</td></tr><tr><td rowspan="6">Apsidal precession</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458162.09119 ± 0.00009</td></tr><tr><td>P<sub>s</sub> [d]</td><td>1.48111973 ± 0.00000017</td></tr><tr><td>e</td><td>0.00122 ± 0.00025</td></tr><tr><td>ω<sub>0</sub> [rad]</td><td>4.426 ± 0.254</td></tr><tr><td>dω/dN [rad/orbit]</td><td>0.001624 ± 0.000300</td></tr><tr><td>BIC</td><td>90.45</td></tr><tr><td rowspan="12">KELT-16 b</td><td rowspan="3">Constant period</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458509.84580 ± 0.00010</td></tr><tr><td>P</td><td>0.96899284 ± 0.00000014</td></tr><tr><td>BIC</td><td>292.90</td></tr><tr><td rowspan="4">Orbital decay</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458509.84586 ± 0.00013</td></tr><tr><td>P [d]</td><td>0.96899283 ± 0.00000014</td></tr><tr><td>dP/dN [d/orbit]</td><td>(-2.73 ± 3.77) × 10<sup>-10</sup></td></tr><tr><td>BIC</td><td>297.63</td></tr><tr><td rowspan="6">Apsidal precession</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2458509.84576 ± 0.00013</td></tr><tr><td>P<sub>s</sub> [d]</td><td>0.96899281 ± 0.00000023</td></tr><tr><td>e</td><td>0.00065 ± 0.00047</td></tr><tr><td>ω<sub>0</sub> [rad]</td><td>2.838 ± 1.685</td></tr><tr><td>dω/dN [rad/orbit]</td><td>0.001441 ± 0.000335</td></tr><tr><td>BIC</td><td>310.95</td></tr><tr><td rowspan="12">WASP-4 b</td><td rowspan="3">Constant period</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2456880.45366 ± 0.00003</td></tr><tr><td>P [d]</td><td>1.33823132 ± 0.00000002</td></tr><tr><td>BIC</td><td>351.03</td></tr><tr><td rowspan="4">Orbital decay</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2456880.45387 ± 0.00005</td></tr><tr><td>P [d]</td><td>1.33823133 ± 0.00000002</td></tr><tr><td>dP/dN [d/orbit]</td><td>(-2.62 ± 0.49) × 10<sup>-10</sup></td></tr><tr><td>BIC</td><td>245.34</td></tr><tr><td rowspan="6">Apsidal precession</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2456880.45359 ± 0.00005</td></tr><tr><td>P<sub>s</sub> [d]</td><td>1.33823140 ± 0.00000008</td></tr><tr><td>e</td><td>0.00079 ± 0.00020</td></tr><tr><td>ω<sub>0</sub> [rad]</td><td>3.492 ± 0.310</td></tr><tr><td>dω/dN [rad/orbit]</td><td>0.001089 ± 0.000103</td></tr><tr><td>BIC</td><td>252.25</td></tr><tr><td rowspan="12">HD 97658 b</td><td rowspan="3">Constant period</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2456361.80580 ± 0.00024</td></tr><tr><td>P [d]</td><td>9.48930773 ± 0.00000149</td></tr><tr><td>BIC</td><td>58.13</td></tr><tr><td rowspan="4">Quadratic fit</td><td>t<sub>0</sub> [BJD<sub>TDB</sub>]</td><td>2456361.80625 ± 0.00053</td></tr><tr><td>P [d]</td><td>9.48929994 ± 0.00000824</td></tr><tr><td>dP/dN [d/orbit]</td><td>(4.20 ± 4.37) × 10<sup>-8</sup></td></tr><tr><td>BIC</td><td>55.05</td></tr></tbody></table>

**Notes.** BIC stands for Bayesian information criterion; the lower the value, the better the fit, but the number of parameters is also taken into account.