Orbital Decay of the Ultra-Hot Jupiter TOI-2109b: Tidal Constraints and Transit-Timing Analysis

The inertial terms in the equation of motion contain the dynamical tide that is forced by the equilibrium tide in convective zones. However, when rotation is present the effect of tides is strongly frequency-dependent, making it a problem harder to solve (e.g., M. Rieutord & L. Valdettaro, 2010).

Tidal interactions between two or more bodies are complex and can affect the equilibrium shape of interacting bodies. Such tidal deformations in a star and planet can be assumed proportional to the their force if the amplitude of the perturbations is small (A. E. H. Love, 1927), and their response to tidal distortion can be quantified via intricate coefficients that, again, strongly depend on the tidal frequency, $\omega$ (e.g., M. Efroimsky, 2012). Those coefficients are the so-called Love numbers, $k_{l}^{m}(\omega)$, for which the imaginary part of the $l=m=2$ component of the tide, $k_{2}^{2}(\omega)$, provides information about the tidal dissipation via IWs⁷⁷7Assuming that the obliquities of the star and the planet are negligible.,

Similar to the formalism of G. I. Ogilvie (2013) that uses an impulsive forcing to calculate IWs, we use here a frequency-averaged formalism to compute the tidal excitation of these waves in convection zones and quantify their typical dissipation. Thus, the previous equation is averaged in $\omega\,\epsilon\,[-2\Omega,2\Omega]$. Prior results have shown a more detailed inclusion of the realistic interior structure of a body (see e.g., F. Gallet et al., 2017; M. Benbakoura et al., 2019). However, using a frequency-averaged dissipation is a coarse estimate that allows us to approach such an intricate problem by providing us a typical quantification of the dissipation produced by IWs. The actual quantification of IW-dissipation⁸⁸8This depends, for example, on the combination of non-linear effects with magnetic fields and/or turbulence in convection zones. This is not addressed here as studying the specific dissipation mechanism of IWs is out of the scope of this work. for a specific tidal frequency in an individual system is subject to large uncertainties (from 2 to 3 orders of magnitude), which is why we use a formalism that is robust and has been widely adopted in previous research.

As shown by some authors (e.g., G. I. Ogilvie, 2013; A. J. Barker, 2020), tidal dissipation due to IWs is significantly more efficient for rapid rotators and will only apply if $|\omega|\leq 2\Omega_{\mathrm{\star}}$ (i.e. while tidal forcing excites IWs), which for circular orbits or small eccentricities will be true only if $P_{\mathrm{orb}}\geq\frac{P_{\mathrm{rot,\star}}}{2}$. Seeking a robust (but straightforward) implementation of tidal dissipation via IWs into the evolution of the TOI-2109 system and using the aforementioned frequency-averaged formalism, we employ a seemingly simple numerical implementation. This is performed with a simplified piece-wise homogeneous two-layer model for both the star (subscript $\star$) and the planet (subscript p). For the star, a fluid-fluid boundary is used between the two layers (S. Mathis, 2015), whereas for the planet it is a solid-fluid boundary (M. Guenel et al., 2014),

where $\epsilon_{\Omega}^{2}\equiv(\Omega/\Omega_{\mathrm{critical}})^{2}$ and $\Omega_{\mathrm{critical}}\equiv(\mathcal{G}M/R^{3})^{1/2}$. Hereafter, $\braket{Q_{\mathrm{\star,\,IW}}^{\prime}}\equiv Q_{\mathrm{\star,\,IW}}^{\prime}$ and $\braket{Q_{\mathrm{p}}^{\prime}}\equiv Q_{\mathrm{p}}^{\prime}$.

Furthermore, A. J. Barker (2020) performed an in-depth study of the dissipation efficiency via IWs for different stellar spectral types, leading to various findings about how these waves are excited in convection zones as well as for which scenarios their dissipation is more dominant. Their modeling, though mathematically speaking does not represent a much higher level of complexity, is more computationally demanding. While it may offer a better representation of stellar physical processes, the simplicity, availability, and widely-adopted use of the two-layer model make it more accessible to a wider range of calculations. However, a detailed comparison of the results achieved for IW-dissipation using a realistic model and a two-layer model showed that the latter underpredicts the dissipation of IWs in F-type stars.

For TOI-2109 where $M_{\mathrm{\star}}\approx 1.4\,\mathrm{M}_{\odot}$ and $t_{\mathrm{\star,age}}$ varies from 1 to 2.5 Gyr, such a comparative analysis reveals that $Q_{\star}^{\prime}$ follows the same qualitative evolution as times goes on, but IW-dissipation is underpredicted in 1–2 order of magnitude in the two-layer model (see fig. 7 in A. J. Barker, 2020). However, it is worth noting that we use the two-layer model not only to quantify IW-dissipation, but also to include the time dependence of $Q_{\star}^{\prime}$ through $\Omega_{\mathrm{\star}}$, which is directly coupled to the evolution equations. Therefore, we account for the quantitative differences of $Q_{\star}^{\prime}$ in the two-layer model by using A. J. Barker (2020) to scale the value of $Q_{\star}^{\prime}$ to the expected level of dissipation in F-type stars such as TOI-2109. By keeping this semi-analytical track of how $Q_{\star}^{\prime}$ changes over time, this approach makes our results easy to be reproduced using the same (or similar) computational methods.

The excitation of IWs by tidal forcing in the convective envelope of TOI-2109 (subsection 3.1) is not the only dissipation mechanism acting in its interior structure. In fact, as IWs in convective envelopes are excited only if $P_{\mathrm{orb}}\geq\frac{P_{\mathrm{rot,\star}}}{2}$, once this condition does not hold true, for instance, when the planet reaches an inner orbit where $P_{\mathrm{orb}}\approx 0.5$ d, the star-planet system will continue to dissipate energy only—at least within this model—through the so-called Internal Gravity Waves (IGWs) produced in the stellar convective/radiative interface. Tidal dissipation via IGWs in TOI-2109 will be significant if IGWs are fully damped, which can happen when IGWs travel and break before reaching another radiative region. In between resonances, travelling waves are the most direct and straightforward way of studying the dissipation of energy via IGWs, thus providing their most efficient tidal dissipation. Moreover, even if resonance exists within tidal forcing, for large enough tidal amplitudes waves can still break and their full damping might still hold true.

Using the formalism presented in A. J. Barker (2020), we study tidal dissipation via IGWs in two scenarios: 1) while $P_{\mathrm{orb}}\geq\frac{P_{\mathrm{rot,\star}}}{2}$ so that IWs and IGWs act together, and 2) when the planet reaches inner positions where IWs are no longer excited (i.e., $P_{\mathrm{orb}}\ngeq\frac{P_{\mathrm{rot,\star}}}{2}$), leaving IGWs as the only active dissipation mechanism for shorter orbital periods. Tidal dissipation due to IGWs in radiative regions strongly depends on the orbital/tidal period. Particularly, after complex numerical calculations coupled to stellar evolution models using MESA, A. J. Barker (2020) found $Q_{\mathrm{\star,\,IGW}}^{\prime}\propto(P_{\mathrm{tide}}/0.5\mathrm{d})^{\frac{8}{3}}$ for main-sequence stars, while $Q_{\mathrm{\star,\,IGW}}^{\prime}\approx 10^{8.5}$ for F-stars of $M_{\mathrm{\star}}\approx 1.45\,\mathrm{M}_{\odot}$. We must, however, mention that the difference between these two cases strongly depends on the parameters of the host star, entailing some caveats that need to be considered within the model adopted for the dissipation of IGWs. Also, the treatment used here neglects the effects of rotation on IGWs. Still, it is important to note that the inclusion of such effects in the tidal frequency $\omega$ can lead to an enhanced efficiency of tidal dissipation via IGWs. Please refer to P. B. Ivanov et al. (2013) for a complete description.

As mentioned before, and further explained in A. J. Barker (2020), tidally excited IGWs will be very efficient if they are in the traveling wave regime due to nonlinear wave breaking in the stellar core. This means that IGWs need to be fully damped, which can happen if radiative diffusion takes place, if large resonant tidal forcing is present (A. J. Barker & G. I. Ogilvie, 2010; A. J. Barker, 2011), and/or if IGWs have large tidal amplitudes so that stratification inside the star can be overturned (see A. J. Barker & G. I. Ogilvie, 2010). Previously, A. J. Barker (2020) found that only the right combination of planetary mass, stellar mass, and stellar age would lead to wave breaking. This condition is enclosed into a mass threshold for the the planet when exceeding the so-called critical mass, $M_{\mathrm{crit}}$, which in turn depends on the interior structure of the star⁹⁹9$M_{\mathrm{crit}}$ varies with $t_{\mathrm{\star,age}}$ and $M_{\mathrm{\star}}$. See fig. 9 in A. J. Barker (2020).. The efficiency of IGWs strongly depends on $M_{\mathrm{crit}}$ and $M_{\mathrm{crit}}$ significantly varies with $t_{\mathrm{\star,age}}$, which for TOI-2109 is highly uncertain (see Table 1). As a result, $t_{\mathrm{\star,age}}$ is the overarching stellar parameter that determines the possible orbital decay rates of TOI-2109b, splitting the results into two different scenarios, as described in Section 4.

Given uncertainties in the age of TOI-2109, the lower-end and central value of the stellar age, $t_{\mathrm{\star,age}}$ (i.e., 1.09 and 1.77 Gyr, see Table 1) do not comply with the requirements to damp IGWs in radiation zones; that is, the planetary companion would have to be much more massive than TOI-2109b, thus indicating that IGWs are not in the wave breaking regime and are instead linear. Our results imply that tidal dissipation due to IGWs in a ‘young’ TOI-2109 will be likely orders of magnitude less efficient at $Q_{\mathrm{\star,\,IGW}}^{\prime}\approx 10^{8.5}$ (cf. A. J. Barker, 2020). That said, most of the dissipation then resulted from the excitation of IWs in the stellar convective envelope, for which we studied the forward orbital evolution of TOI-2109b using the model described in subsection 3.1.

Figure 4 shows the evolution of the semi-major axis for different values of $Q_{\mathrm{\star,\,IW}}^{\prime}$, where tidal dissipation arises mostly from excited IWs. From these tidal simulations we calculate the decay rate of the orbital period as $\dot{P}_{\mathrm{TIDE}}=(-4.186\pm 0.797)\,\mathrm{ms}\,\mathrm{yr}^{-1}$ at $3\sigma$. Using the formulation of P. Goldreich & S. Soter (1966) and $Q_{\star}^{\prime}=3Q_{\star}/(2k_{2,\star})$, this leads to 95% confidence lower limit of the decay timescale, $\tau=\frac{P}{\dot{P}}$, $\tau_{\mathrm{TIDE}}\gtrsim 6$ Myr, and $95\%$ confidence lower limit of $Q_{\mathrm{\star,\,TIDE}}^{\prime}>2.3\times 10^{7}$. The latter is in agreement with the expected values of a type-F star with $1.445\,\mathrm{M}_{\odot}$ (A. J. Barker, 2020), as well as with the expected aspect ratios $\alpha_{\mathrm{\star}}$ and $\beta_{\mathrm{\star}}$ (Figure 3) in the current scenario of stellar tidal dissipation for a ‘young’ TOI-2109 (cf. F. Gallet et al., 2017).

It is particularly interesting that, despite this system is an ideal candidate to study orbital decay in USP Jupiters, the aforementioned decay rate for TOI-2109b might be an indicator of a constant orbital period. Furthermore, a direct comparison of the decay rate calculated here to those predicted by I. Wong et al. (2021) (i.e., 10–740 $\mathrm{ms}\,\mathrm{yr}^{-1}$), indicates that the lower end might still hold while the upper end can be ruled out. As shown later in Section 5, this is further supported by observational constraints using space-based ¹⁰¹⁰10This includes TESS observations from three sectors (see Table 2). and ground-based data of TOI-2109.

As mentioned before, the dissipation of IGWs is driven by the combination of stellar age, stellar mass, and planetary mass. These, together, control the damping of such waves in the stellar radiative regions. In general, the damping of IGWs in young stars requires extremely large planetary masses for IGWs to be in the wave breaking regime. This is why in the previous section the excitation of IGWs was significantly inefficient and most of the tidal dissipation was driven by IWs in the stellar convective envelope. However, it is also worth noting that all stars with $M_{\mathrm{\star}}>0.9\,\mathrm{M}_{\odot}$ will present a steady decrease in $M_{\mathrm{crit}}$ before reaching 10 GYr (A. J. Barker, 2020). This means that for a given stellar mass, IGWs might be fully damped in a later stage of stellar evolution as $M_{\mathrm{crit}}$ decreases to lower values.

In light of the above, we notice that if the stellar age of TOI-2109 is closer to the upper-end of the 1-sigma range (i.e., $t_{\mathrm{\star,age}}=2.65$ Gyr, see Table 1), IGWs might break in the stellar radiation region as $M_{\mathrm{crit}}$ would be between 1 and 10 Jupiter masses. Following A. J. Barker (2020), we notice that a planetary mass of $\sim 5\,\mathrm{M}_{\mathrm{J}}$ would make IGWs that are launched from the radiative/convective regions of the star go into the wave breaking regime and subsequently be fully damped. An ‘old’ TOI-2109 would follow $Q_{\mathrm{\star,\,IGW}}^{\prime}\propto(P_{\mathrm{tide}}/0.5\mathrm{d})^{\frac{8}{3}}$, presenting a more efficient tidal dissipation due to IGWs. In Figure 5, we show the semi-major axis as a function of time and $Q_{\mathrm{\star}}^{\prime}$. The latter is now a significant combination of $Q_{\mathrm{\star,\,IW}}^{\prime}$ and $Q_{\mathrm{\star,\,IGW}}^{\prime}$, noting that during half of the orbital evolution both dissipation mechanisms (IWs and IGWs) acted together.

Once $P_{\mathrm{orb}}\ngeq\frac{P_{\mathrm{rot,\star}}}{2}$ (i.e., no excitation of IWs), IGWs acted alone and dominated the angular momentum transport in the star. This behavior was expected because IGWs, if efficiently damped, produce the strongest tidal dissipation rates for short orbital periods (A. J. Barker, 2020). From the tidal simulations presented in Figure 5, we calculated the orbital decay rate as $\dot{P}_{\mathrm{TIDE}}=(-1107\pm 648)\,\mathrm{ms}\,\mathrm{yr}^{-1}$, resulting in a decay timescale $\tau_{\mathrm{TIDE}}\gtrsim 17$ kyr for the $Q_{\mathrm{\star}}^{\prime}$ that were studied. The previous decay rate puts a lower limit of $Q_{\mathrm{\star,\,TIDE}}>8.7\times 10^{4}$, which for the adopted values of $\alpha_{\mathrm{\star}}$ and $\beta_{\mathrm{\star}}$ is close to the lowest-limit of $Q_{\mathrm{\star}}^{\prime}$ in F-stars, ranging from $10^{5}$ to $10^{7}$ (G. I. Ogilvie & D. N. C. Lin, 2007; A. F. Lanza et al., 2011).

Distinct values of $C$ in $Q_{\mathrm{\star,\,IGW}}^{\prime}=C\,(P_{\mathrm{tide}}/0.5\mathrm{d})^{\frac{8}{3}}$ could change the previous lower limits, which corresponded to a $C=3$. A. J. Barker (2020) found that $C$ can vary from 1 to 3 for main-sequence stars with 0.5–1.1$\,\mathrm{M}_{\odot}$, thus expecting $C\gtrsim 3$ for TOI-2109. To quantify these changes, we increased $C$ by an inflated factor of two (as $C$ would just be slightly larger), which led to a decrease of $\sim 34\%$ in $\dot{P}_{\mathrm{TIDE}}$, and an increase of $\sim 45\%$ in $\tau_{\mathrm{TIDE}}$ and $Q_{\mathrm{\star,\,TIDE}}^{\prime}$. As studied in the next Section, however, the results from this scenario—regardless of $C$—are still vastly inconsistent with the constraints from transit-timing observations reported to date. Hence, seeking conservative values in this study case, we kept $C=3$ in all our simulations for an ‘old’ TOI-2109.

In comparison to USP-J systems such as NGTS-10b and WASP-19b that could spin up the rotational rate of their host star by about 30 per cent (see J. A. Alvarado-Montes et al., 2021), we found that the rotational rate of TOI-2109 was not significantly spun up during the shrinking of the orbit of its planetary companion. This may be a consequence of the moment of inertia of TOI-2109 (i.e. $I_{\mathrm{\star}}\propto M_{\mathrm{\star}}R_{\mathrm{\star}}^{2}$), which is significantly higher than that of the previously mentioned systems. All our dynamical simulations resulted in final stellar rotation periods longer than the saturation rotation period of TOI-2109 (i.e. $\bar{P}_{\mathrm{\star}}=0.559$ d), using the saturation rate $\bar{\Omega}_{\mathrm{\star}}$ provided by A. Collier Cameron & J. Li (1994). We conclude that during the orbital evolution of TOI-2109b, the host star’s magnetic activity remains a function of the stellar rotational rate and a constant stellar activity level is not achieved.

TOI-2109 or TIC 392476080, as listed in the TESS Input Catalog (TIC, K. G. Stassun et al., 2018) has been observed by TESS in: 2020 (May 13 to June 8, Sector 25), 2022 (May 18 to June 13, Sector 52), and 2024 (May 21 to Jun 18, Sector 79). All TESS data was downloaded from the Mikulski Archive for Space Telescopes (MAST¹¹¹¹11https://mast.stsci.edu/) using Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration et al., 2018). In Table 2 we list the TESS data we processed in this work, indicating the source, cadence, and sector. The simple aperture photometry (SAP) light-curve files¹²¹²12These consist of the flux extracted from the optimized aperture on the full-frame images (FFI). were previously extracted via the Science Processing and Operations Center (SPOC) as part of the TESS Light Curves from Full Frame Images (TESS-SPOC) High Level Science Products project (D. A. Caldwell et al., 2020).

We process the SAP light curves ourselves and measure the mid-transit time of each observation, using the formalism and code for transit-timing analysis presented in E. S. Ivshina & J. N. Winn (2022)—hereafter IV22—which has been employed to calculate transit times in several other systems. We input the ephemeris of TOI-2109b from I. Wong et al. (2021) into the transit-timing code of IV22¹³¹³13This code was originally designed for data with a nominal cadence of 2, 10, and 30 minutes. The formalism, however, can be easily adapted to be used with a cadence of 20 seconds (i.e., for the data from sector 79)., which analyses individual transit intervals by only selecting transit data with optimal time coverage and removing any transit interval with a lack of data points. Occasional outliers were removed by following the Median Absolute Deviation (MAD) of residuals used in IV22, which is also employed in all the photometric data to identify and reject transits of unusually poor quality. We present the flux time series of TOI-2109 from TESS Sector 52, obtained using IV22 and with all selected transits for transit-timing analysis marked in red (Figure 6).

To remove photometric trends such as stellar variability, IV22 fits a model to the out-of-transit data and uses it to normalize ("rectify") the entire light curve. The individual transits are then phase-folded (Figure 7) to extract key parameters that describe the geometry of the transits, including $R_{\mathrm{p}}/R_{\mathrm{\star}}$, $a/R_{\star}$, impact parameter ($b$), and limb-darkening law. Once reliable values for these parameters were obtained from the best-fit transit model of the phase-folded light curve, they were used to fit each individual transit and determine the corresponding mid-transit time. In Figure 8, we show the selected transits of TOI-2109b for Sector 52 along with the best-fit transit model (in red) and residuals. Each subplot in this figure is titled with the number of the orbit, as calculated from the zero epoch in I. Wong et al. (2021). To estimate the mid-transit times (and uncertainties) of the selected transits, we used a Markov Chain Monte Carlo (MCMC) method. We applied each step of the previous process¹⁴¹⁴14For a more detailed explanation, see section 3 (TESS Data Analysis) in E. S. Ivshina & J. N. Winn (2022). to the TESS data in Table 2, measuring 31, 21, and 22 reference mid-transit times for Sectors 25, 52, and 79, respectively. These will be refined later in subsection 5.3 for transit-timing analysis.

From 2021 to 2023, TOI-2109 was observed 31 times by CHEOPS, 26 of which contain transits with significant coverage. The data is fully accessible, for example, via the CHEOPS archive¹⁵¹⁵15https://cheops-archive.astro.unige.ch/archive_browser/ through the Data and Analysis Center for Exoplanets (DACE) website¹⁶¹⁶16https://dace.unige.ch/cheopsDatabase/? of the University of Geneva. However, we did not process the data ourselves but instead adopted the results of J. V. Harre et al. (2024)¹⁷¹⁷17Available at https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/692/A254, who calculated the 26 reference mid-transit times in the CHEOPS data. They extracted these values from their best-fit model after preparing and reducing the data with two different photometry pipelines: the CHEOPS Data Reduction Pipeline (DRP; S. Hoyer et al., 2020) and the PSF Imagette Photometric Extraction (PIPE; see e.g., A. Brandeker et al., 2022, and references therein). The resulting mid-transit values from both photometry pipelines are systematically within 1-sigma of each other’s uncertainties. Since both sets of values are not substantially different from each other, we chose to work with the DRP mid-transit values (see table C.1 in J. V. Harre et al., 2024).

Long-term variations in the orbital period of TOI-2109b can only be detected by having transit data from a well-sampled, long-time baseline. As mentioned before, ever since TOI-2109b was discovered back in 2020, TESS has only observed it three times. Therefore, to have a four-year baseline that is better sampled to do transit-timing analysis, we use archival data obtained from The Exoplanet Follow-up Observing Program (ExoFOP). These include the 20 partial-transit observations of I. Wong et al. (2021) collected by multiple ground-based telescopes: the Fred Lawrence Whipple Observatory (FLWO, 1.2 m) on Mt. Hopkins, Arizona, USA; the University of Louisville Manner Telescope (ULMT, 0.61 m) at Mt. Lemmon, Arizona, USA; Las Cumbres Observatory Global Telescope network (LCOGT, 0.4–1.0 m; T. M. Brown et al., 2013); the Maury Lewin Astronomical Observatory (MLO, 0.36 m) in California, USA; the Telescopio Carlos Sánchez (TCS, 1.5 m) at Teide Observatory, Spain; the Wild Boar Remote Observatory (WBRO, 0.24 m) near Florence, Italy; the Grand-Pra (GdP, 0.4 m) Observatory in Switzerland; and the Faulkes Telescope North (FTN, 2.0 m) at Haleakala Observatory on Maui, Hawai’i. The ground-based data of TOI-2109 covers from 2020 to 2023, with a cadence as low as 3 s and using different filters, as shown in Table 3.

We independently analyzed each dataset at their corresponding epoch, using the PyTransit package (H. Parviainen, 2015) for the light-curve modeling. For each dataset, a K. Mandel & E. Agol (2002) model was fitted to the phase-folded light curve. The parameters that we used for this modeling included planet-to-star radius ratio $R_{\mathrm{p}}/R_{\mathrm{\star}}$, stellar density in g cm^-3, impact parameter $b$, orbital period $P$, and zero epoch $T_{0}$. To fit the eccentricity and argument of periastron, PyTransit used the parametrization $\sqrt{e}\,{\rm cos}\omega$ and $\sqrt{e}\,{\rm sin}\omega$. We followed a similar process to that of subsection 5.1, that is, we fitted first the phase-folded light curve to extract reliable values of the parameters mentioned before and then use these values to fit each individual transit, thus recovering its reference mid-transit time. We fully process all the ground-based data ourselves, recovering 20 reference mid-transit times.

Data from the All-Sky Automated Survey for Supernovae (ASAS-SN) were available for TOI-2109, which spans for over 10 years. However, these observations have significant variations in cadence and most of them do not even complete a partial transit. For the purpose of calculating mid-transit times, at least a partial transit is needed to identify the mid-point and thus recover the transit time. For this reason, we decided to neglect ASAS-SN data and only use data from other dedicated follow-up missions, as previously mentioned. Moreover, SuperWASP-N (Wide Angle Search for Planets; D. L. Pollacco et al., 2006) observed TOI-2109 in several fields, from 2006 to 2011. We do not process this data ourselves but instead use the only reference mid-transit time recovered by J. V. Harre et al. (2024), as most of the data was filtered out due to poor transit-timing quality caused by bad transit coverage and/or high scatter. This data point, however, is significant and provides a longer baseline that can allow us to further constrain the orbital decay rate of TOI-2109b.

Transit-timing variations (TTVs) can be used as possible indicators of orbital decay. TTVs can be recovered by comparing the reference mid-transit times with calculated, refined values that are obtained from a fitting process to the transit-timing data. For TOI-2109b, we use the reference mid-transit times previously found (subsections 5.1 and 5.2) into a quadratic ephemeris model described as follows

where $T_{0}$, $P_{0}$, $N$, and $T_{N}$ is the reference mid-transit time, the initial orbital period of the planet, the number of orbits since the reference transit, and the calculated mid-transit time, respectively. The constant rate of change of the orbital period, $\dot{P}$, in Equation (20) implies that there is a small difference between the orbital period of two successive transits. In this model, we fit the free parameters $T_{0}$, $P_{0}$, and $\dot{P}$, taking initial guesses for $T_{0}$ and $P_{0}$ from I. Wong et al. (2021).

To fit the previous model to our reference transit-timing data¹⁸¹⁸18This can be found at https://github.com/JAAlvarado-Montes/TOI-2109b, we use the software package PdotQuest provided by W. Wang et al. (2024). By inputting the reference mid-transit times of TOI-2109b and their uncertainties into PdotQuest, we calculated each new mid-transit time and associated uncertainty, thus recovering $\dot{P}$ in the process. Internally, PdotQuest uses the expression $\frac{{\rm d}P}{{\rm d}t}=\frac{1}{P}\frac{{\rm d}P}{{\rm d}N}$ to adopt the quadratic form of equations (19) and (20) as in G. Maciejewski et al. (2021).

It is worth noting that, despite TESS provides the majority of mid-transit times, the photometric quality of TESS data hinders the transit-timing analysis. An extreme example of this is the data from TESS Sector 25 where some of the reference mid-transit times had uncertainties up to $\sim$ 10 min. For a better fit to the transit-timing data, any mid-transit points with $\Delta T_{0}\gtrsim 2$ min were clipped out before using PdotQuest. This allowed us to use a wider variety of data sources while only keeping the data with a higher timing precision for a more robust fit. This process filtered out 8 mid-transit times from the ground-based data and 6 from TESS data (5 from sector 25 and 1 from sector 79). To avoid including data with high scatter, PdotQuest uses an iterative data-clipping technique to remove outliers in the TTVs (see section 3 in W. Wang et al., 2024), also called $n$-sigma rejection areas. Such outliers are represented as gray dots in Figure 9, which shows TTVs of TOI-2109b for a $1\sigma$ (upper panel) and $3\sigma$ (lower panel) rejection area, respectively. While the $3\sigma$ fit used 102 mid-transit times out of the 117 that were provided—including several high-scatter data points—the results from this fit are completely consistent with the $1\sigma$ calculation that only uses the 53 more precise mid-transit times.

The complete results from our transit-timing analysis are reported in Table 4. We can see how, for both sigma-rejection areas, our calculated orbital decay rate is only compatible with the lower limit of the range 10–740 $\mathrm{ms}\,\mathrm{yr}^{-1}$ first reported by I. Wong et al. (2021), ruling out any $\dot{P}\gtrsim 10\,\mathrm{ms}\,\mathrm{yr}^{-1}$. Also, J. V. Harre et al. (2024) calculated an orbital decay rate of just $\dot{P}\sim-1\,\mathrm{ms}\,\mathrm{yr}^{-1}$ by using ground- and space-based data. Our best-fitting value of $\dot{P}$ is consistent with these results and is based on one additional sector of TESS data (i.e., sector 79) that was unavailable in previous works (I. Wong et al., 2021; J. V. Harre et al., 2024). With a a cadence of 20 seconds, this dataset provided a well-mapped transit timeseries with an improved precision on the mid-transit times as compared to those from previous TESS sectors (i.e., 25 and 52). These 22 new data points from TESS sector 79 were included in the TTV analysis of TOI-2109b and allowed us to further constrain its orbital decay rate, yielding a value of $\dot{P}_{\mathrm{OBSV}}=(-2.616\pm 1.285)\,\mathrm{ms}\,\mathrm{yr}^{-1}$ at $3\sigma$. This allowed us to put a lower limit $Q_{\mathrm{\star,\,OBSV}}^{\prime}>3.7\times 10^{7}$, which is in agreement with those of F-stars that range from $10^{5}$ to $10^{7}$.

TTVs in a compact, ultra-short-period planetary system such as TOI-2109b can arise from multiple physical mechanisms, including gravitational interactions with additional planetary companions (J. V. Harre et al., 2024), general relativistic effects, and deviations from spherical symmetry (oblateness) of both the planet and the host star (see, e.g., M. Sucerquia & N. Cuello, 2025, and references therein). To quantify these effects, we performed numerical simulations using the REBOUNDx package (H. Rein & S. F. Liu, 2012; D. Tamayo et al., 2020), systematically evaluating each mechanism’s contribution to the TTV signal. First, we simulated gravitational perturbations induced by a proposed planetary companion with a mass of $0.2\mathrm{M}_{\mathrm{J}}$, a period of $P_{\mathrm{c}}\sim 1.125$ days, and a moderate eccentricity ($e=0.05$). For TOI-2109b, we adopted $e=0.04$, consistent with the parameter space explored by J. V. Harre et al. (2024). The resulting TTVs exhibit a clear, short-term periodic signature with amplitudes of approximately 1–2 minutes, emphasizing the companion’s significant gravitational influence. Importantly, this model does not include any tidal orbital decay; instead, it illustrates how precession and dynamical perturbations alone can generate a quasi-secular TTV signal that may superficially resemble decay, highlighting the need for caution in its interpretation.

To evaluate the impact of rotationally induced oblateness, we calculated the gravitational quadrupole moment $J_{2}$ for TOI-2109b assuming synchronous rotation (rotational period equal to orbital period, $P_{\mathrm{rot}}=P_{\mathrm{orb}}$). This yielded $J_{2}\approx 3.119\times 10^{-4}$, indicative of significant oblateness comparable to Saturn. Additionally, we included stellar oblateness characterized by a similar quadrupole moment calculation (assuming a stellar rotation period of about 1.14 days and a Love number $k_{2}=0.03$), and general relativistic corrections appropriate for the system’s compactness and eccentricity.

We performed these simulations using the fast Wisdom-Holman symplectic integrator WHFast (J. Wisdom & M. Holman, 1991; H. Rein & D. Tamayo, 2015), with a timestep of $1\times 10^{-5}$ years. This integrator is particularly suited for compact systems, offering high numerical stability and accuracy over long integration times. Transit times in our simulations were precisely determined by detecting intersections of the planet with the $y=0$ plane for $x>0$, refined with a bisection algorithm to achieve a temporal accuracy of $10^{-8}$ years.

The resulting simulated transit times were analyzed by removing a linear trend through least-squares fitting, thus isolating the variations attributable to the studied dynamical mechanisms. All physical effects were consistently incorporated via REBOUNDx’s specialized modules. General relativity was included using the gr_potential module, while the oblateness of both planet and star were implemented using the gravitational_harmonics module. For the planet, synchronous rotation dictated the above-mentioned $J_{2}$ and $J_{4}=-J_{2}/5$ values.

To account for the stellar quadrupole moment, we adopted a rotational period of 1.4 days and a stellar Love number $k_{2}=0.03$. The corresponding $J_{2}^{\star}$ was computed following

and the octupole moment was set as $J_{4}^{\star}=-J_{2}^{\star}/5$, consistent with the quadrupolar potential expansion. This value captures the contribution of stellar oblateness to the long-term orbital precession of TOI-2109b.

Although Equations (21) and (22) share a similar functional form—both scaling as $\Omega^{2}R^{3}/(GM)$—they originate from different physical considerations. Equation (21) is derived from the hydrostatic equilibrium of a rotating fluid planet with a simplified internal structure (e.g., a polytrope of index $n=1$), yielding a direct expression for $J_{2}$ without requiring additional structural parameters (e.g., C. D. Murray & S. F. Dermott, 2000). In contrast, Equation (22) incorporates the stellar Love number $k_{2}$, which quantifies the star’s response to rotational deformation and encapsulates details of its internal density distribution and stratification (e.g., Z. Kopal, 1960; V. N. Zharkov & V. P. Trubitsyn, 1978). The use of $k_{2}$ makes the stellar expression more general and applicable to objects with complex internal profiles. Despite their formal similarity, the coefficients and underlying assumptions of each equation reflect the differing physical regimes of gaseous giant planets and rotating stars.

Figure 10 shows the resulting TTVs under different dynamical scenarios. The green dots represent TTVs caused solely by gravitational interactions with the additional planetary companion. The orange dots incorporate both stellar and planetary oblateness effects along with general relativistic precession, showing a slower, secular trend with smaller amplitudes. Finally, the blue dots integrate all these mechanisms, revealing a signal composed of both periodic, short-term variations from the planetary companion and longer-term secular drift from oblateness and relativistic effects.

Our results reveal distinct dynamical signatures in the TTV signal: short-term variability is primarily driven by the gravitational perturbations of an external planetary companion, while stellar and planetary oblateness, together with general relativistic precession, induce subtle but coherent secular trends. This decomposition is essential for accurately interpreting observed TTVs and disentangling the underlying physical mechanisms shaping the orbital architecture of compact systems like TOI-2109b—particularly in the context of long-duration, high-precision monitoring.