© The Authors 2026

1 Introduction

Tidal interactions play a crucial role in our understanding of close planetary systems. They drive changes in key orbital parameters (Barker 2020; Ogilvie 2014). Studies have theorised on their role for decades, starting from tight binary star systems (Zahn 1966, 1977, 1989). Tidally driven changes in orbits of exoplanets have been inferred statistically (Schlaufman & Winn 2013; Hamer & Schlaufman 2019) and, more recently, a short-lived optical burst was explained on the basis of a planet inspiralling and finally plunging into its host star (De et al. 2023). However, direct observational confirmation of inspiralling was reached only for WASP-12 b (Maciejewski et al. 2016; Yee et al. 2020; Winn & Stefánsson 2025). This scarcity of detections drives a high demand for more detections of inspiralling planets, to provide better constraints on theoretical models.

Planets and sub-stellar objects on tight orbits interact with their host stars gravitationally, which creates various kinds of responses depending on the interior structure of the host. In convective layers, these interactions are described by equilibrium tides, which are semi-hydrostatic deformations of the star, dissipated through convection. We expect the dominant response in radiative layers to be dynamical tides, especially internal gravity waves (IGW) dissipated by wave breaking (WB). For an overview, we refer to Barker (2020), hereafter B20. In this study, we focus on a scenario with the breaking of the IGW in radiative stellar centres, as this was shown to be the strongest interaction operating in solar-type main sequence stars (B20). Stars of this type are the most commonly encountered hosts of extrasolar planetary systems. The strength of tidal interactions is characterised with the reduced tidal quality factor Q^,, which is defined as the ratio of energy stored within the tidal distortion to the energy dissipated within one tidal period. Therefore, the lower the value of Q^,, the higher the strength of tidal interactions. Alongside Q^,, this theory also provides the critical mass, M_crit, which is the minimal mass of a planet required for the IGW amplitude to overcome stable stratification at the host’s centre and induce wave breaking. If the mass of a planet, M_p, is greater than M_crit in a given system, then wave breaking is predicted to occur. As this model requires the waves to propagate in radiative zones, we assume that it only occurs if the star’s core is radiative and, consequently, no tidal dissipation occurs while it is convective. While equilibrium tides and inertial waves also contribute to the overall tidal dissipation budget (and may dominate in systems where wave breaking does not operate) their effect in the systems analysed here is expected to be significantly weaker than that of wave breaking (B20).

Tidal interaction can be detected in a given planetary system by observing variations in the planet’s orbital parameters on long enough timescales. The best method of testing tidal interactions is transit timing variations (TTV), which rely on detecting an earlier arrival of transits due to a shrinking of the planet’s orbit. There is a significant number of studies trying to find this effect in individual systems (e.g. Southworth et al. 2022; Biagiotti et al. 2024), small samples of systems (e.g. Maciejewski et al. 2022; Harre et al. 2023), as well as rich populations (e.g. Ivshina & Winn 2022). The challenge in detecting orbital shrinking through TTV is acquiring high-quality photometric data spanning many years, even decades. Additionally, several other astrophysical effects can cause a shift in transit central times, producing false positives, such as the Applegate effect (Applegate 1992; Watson & Marsh 2010), the light travel effect (Irwin 1952), apsidal precession (Miralda-Escudé 2002), and perturbations from other objects in the system (Agol et al. 2005). The detection itself can be done by detecting a statistically significant quadratic trend in the shift of central times of transits following a quadratic ephemeris, $$T_{\mathrm{c}}=T_0+P\times E+\frac{1}{2}\frac{\mathrm{d}P}{\mathrm{d}E}\times E^2,$$(1)

where T_c is the central time of the transit, T₀ is the reference time, P is the orbital period of the planet, E is the epoch, which is the number of the transit counted from T₀, and the quadratic term $\frac{\mathrm{d}P}{\mathrm{d}E}$ is the change of orbital period with epoch. The quadratic term can be directly compared with theoretical predictions using the relation presented in Eq. (5) of Maciejewski et al. (2018): $$Q^{\prime }=-\frac{27\pi }{2}\frac{M_{\mathrm{p}}}{M_{\star }}{\left(\frac{a}{R_{\star }}\right)}^{-5}{\left(\frac{\mathrm{d}P}{\mathrm{d}E}\right)}^{-1}P,$$(2)

where $\frac{M_{\mathrm{p}}}{M_{\star }}$ is the planet to star mass ratio and $\frac{a}{R_{\star }}$ is the semi-major axis of the orbit of the planet in host radii.

Additionally, there are theoretical predictions for orbital changes driven by tides to be detectable in radial velocities (Arras et al. 2012; Alvarado et al. 2024). However, no observational confirmation has been reported thus far.

We aim to advance the field by applying detailed theoretical predictions of tidal interactions, based on the IGW dissipation framework assuming total wave breaking, to a large and diverse sample of planetary systems. This approach enables us to identify promising targets for future observational follow-up and to interpret both existing detections and non-detections within a unified theoretical context. For each system studied, we provide specific predictions that can be directly compared with observational results. By analysing a statistically significant sample, we are also able to derive population-level insights into tidal dissipation and orbital evolution. Given the current scarcity of confirmed inspiralling planets, there is a pressing need to identify the most promising candidates for such systems, which we address in this study.

In Section 2, we present the selection process of the studied systems. Section 3 describes the method of calculating the tidal quality factor under the WB regime together with the mass criterion and related parameters, and presents representative examples. In Section 4, we examine the population-level trends, identify the most promising targets for follow-up observations, and discuss individual noteworthy systems. The main conclusions are summarised in Section 5.

2 Sample selection

Our initial sample was drawn from the exoplanet.eu catalogue, downloaded on 23 October 2024, which contained 7341 confirmed planets¹. Systems hosting planets with masses above 0.1 M_J and orbital periods shorter than 10 days (Wright et al. 2012) were selected, yielding 710 hot Jupiters.

From this list, we excluded non-transiting systems, objects with incomplete data, systems with stellar parameters outside the range of our pre-computed stellar evolution model grid (Sect. 3.2), and those with variable host stars. Very young systems were also removed, as the homogeneous isochrone-based approach adopted here does not reliably constrain stellar parameters at early evolutionary stages. The final sample comprises 550 planetary systems, reflecting both data availability and the limits of our stellar model grid.

The catalogue data, supplemented by NASA Exoplanet Archive, were used to extract the orbital period, discovery date, T_dsc, radial velocity amplitude, K, orbital inclination, i, the semi-major axis of the planet in host radii, $\frac{a}{R_{\star }}$, and the planet to star radius ratio, $\frac{R_{\mathrm{p}}}{R_{\star }}$. To ensure reliable estimates of stellar effective temperature, T_eff, surface gravity, log(g), and metallic-ity, [Fe/H], we combined the gathered data with measurements from the SWEET-Cat catalogue (Santos et al. 2013; Sousa et al. 2018, 2021, 2024). The values of log(g) derived directly from SWEET-Cat observations (i.e. excluding those adopted from the literature) were corrected, following Eq. (4) of Sousa et al. (2021), as recommended by the authors.

For some systems, we manually searched the literature to fill data gaps and apply minor corrections. We also adopted Gaia parameters from SWEET-Cat, including parallax Plx and magnitudes in the G, BP, and RP bands, all with reported uncertainties. Missing photometric values were supplemented using Gaia DR3 data retrieved via SIMBAD². All stellar and planetary parameters adopted for our modelling are listed in Table A.1.

3 Data modelling and results

Tidal interaction models can be highly sensitive to small variations in key stellar parameters such as mass and age. Even within the reported uncertainties of these parameters, the corresponding changes in the tidal quality factor induced by wave breaking (denoted as Q′_WB, hereafter) can span up to an order of magnitude. To account for this, our analysis incorporates the uncertainties of the stellar parameters into the predictions of tidal dissipation under the WB regime.

The overall procedure, described in detail in Sects. 3.1-3.5, can be outlined as follows:

1. Redetermination of stellar parameters: spectroscopic and photometric measurements were used as inputs for isochrone fitting to derive updated stellar parameters, such as mass, radius, and age, together with their uncertainties, in a homogeneous and self-consistent manner (Sect. 3.1).

2. Constructing the grid of stellar evolution models: to accelerate the main tidal calculations, we computed a dedicated grid of stellar evolution models covering a broad range of masses and metallicities (Sect. 3.2).

3. Modelling tidal dissipation: for each system, 100 random realisations of the stellar parameters were drawn from the posterior distributions of the isochrone fits. For each realisation, we computed Q′_WB, M_crit, and the orbital decay rate, Ṗ. These results were then used to derive the median values and uncertainties of these quantities within a 1 σ age range, as well as the probability of wave breaking f_crit, inferred from the M_crit distribution (Sect. 3.3).

4. Transit timing diagnostics: from the predicted orbital decay rates, we derived key timing diagnostics, including the epoch when the cumulative transit-time shift reaches 5 minutes, to aid in assessing observational detectability (Sect. 3.4).

5. Orbital evolution: the long-term orbital evolution driven by tidal dissipation through wave breaking was simulated to predict the fate of each planet and its destruction timescale (Sect. 3.5).

3.1 Stellar parameter fitting

We used the isochrones (Morton 2015) Python package to translate the observed stellar quantities in our sample into the fundamental physical parameters of the host stars, with the results shown in Table A.2. The input parameters were the effective temperature, T_eff, surface gravity, log(g), metallicity, [Fe/H], Gaia parallax, Plx, and the magnitudes in Gaia G, BP, and RP bands. To accelerate convergence, we adopted flat priors on distance between 0 and 10 kpc and interstellar reddening between 0 and 2 magnitudes in the G band. Bounds on equivalent evolutionary points (EEPs, Dotter 2016) were set between 203 to 1700 to exclude pre-zero-age main sequence (pre-ZAMS) and late-stage configurations.

The fits were performed using nested sampling implemented through MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2019), run with the default settings of the isochrones package. The best-fitting stellar parameters and their uncertainties were derived from the posterior distributions produced by isochrones. To propagate the uncertainties in stellar mass and metallicity into our tidal dissipation predictions, we drew 100 unique realisations of the stellar parameters from the posterior distributions for each system and used them in the subsequent analysis.

3.2 Grid of stellar evolution models

The wave-breaking formalism introduced in B20 requires a onedimensional model of the stellar interior. For this purpose, we used the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). Using MESA offers two advantages: it enables direct comparison with the results of B20, who also employed MESA for his stellar models, and ensures consistency with the stellar parameter fitting (Sect. 3.1), which is based on the MESA isochrones and stellar tracks (MIST, Choi et al. 2016).

To compute the wave-breaking predictions, we used MESA models corresponding to the appropriate stellar mass and metal-licity. Since these parameters recur across many systems in our sample, we constructed a grid of models instead of computing individual ones for each system’s realisation. These models were used by selecting the closest one to the target stellar parameters, without interpolation. To achieve sufficient accuracy, the grid was made dense, with steps of 0.01 M_⊙ in stellar mass and 0.01 dex in metallicity. It spans from 0.7 to 1.5 M_⊙ and from −0.4 to +0.6 dex in [Fe/H], yielding 8000 MESA evolutionary models and covering the bulk of the host stars in our sample. This grid-based approach enables the efficient modelling of a large number of configurations, while minimising computational cost.

Our additional tests confirmed that using grid-based approximations introduces no measurable difference compared to dedicated MESA models.

To maintain compatibility between our grid and the MIST evolutionary tracks, we adopted several input settings from the MIST configuration³. We made a compromise between reproducing MIST data accurately and ensuring numerical efficiency. The final set of parameters used for the MESA models is given in Appendix C.

3.3 Tidal dissipation modelling

To calculate Q′_WB and M_crit, we followed the formalism of Eqs. (41)-(43) and (45)-(48) from B20. We assumed planets orbit their host stars on circular, aligned orbits and that the stellar rotation is slow. Under these assumptions, only the tidal components with a mode degree of l = 2, azimuthal order of m = 2, and tidal frequency of ω = 4π/P need to be considered. For each stellar model, we extracted the interior radial profiles of density, pressure, mass, and Brunt-Väisälä frequency from a single MESA profile. The boundary between the radiative and convective zones was determined from the Brunt-Väisälä frequency profile. Local derivatives at this interface were then computed and used in Eqs. (42) and (43) of B20. The Eulerian gravitational potential perturbation Φ′ from Eq. (15) of B20 was computed from Eq. (12) from B20, using Eqs. (13) and (14) of B20 as boundary conditions. The equilibrium displacement field ξ_e,r is solved using Eq. (15) of B20, with a tidal potential of Ψ = Ar² (Barker 2020), where A is the tidal amplitude calculated using Eqs. (45) and (46) of B20, and r is the distance from the star’s centre. Then, ξ_e,r is used, from the right side of Eq. (3) of Goodman & Dickson (1998), to compute ξ_d,r, which allows us to calculate its local derivative at the zones interface. These quantities allowed us to solve Eqs. (42) and (43) of B20 and then evaluate Q′_WB using Eq. (41) of B20.

The critical mass, M_crit, was derived from Eqs. (45)-(48) of B20. We calculated the tidal amplitude A using Eqs. (45) and (46) of B20. Following Eq. (48) of B20, we fit a linear function to the Brunt-Väisälä frequency profile near the stellar centre to determine the stratification coefficient, C. These allowed us to solve Eq. (47) of B20 yielding $A_{nl}^2$, which indicates whether the tidal amplitude, A, is high enough to cause wave breaking. Rearranging this relation together with Eq. (46) of B20 gave an expression for M_crit, representing the minimum planetary mass required to trigger wave breaking.

To ensure consistency with fitted values of host mass, M_*, we re-determined the planet’s mass M_p based on the radial velocity amplitude, K, using the expression (based on Eq. (5) from Maciejewski et al. 2022): $$M_{\mathrm{p}}=4.917\times {10}^{-3}K{P}^{1/3}{M}_{\star }^{2/3}{\sin }^{-1}i,$$(3)

where K is given in m s⁻¹, M_* is in solar masses, M_⊙, and the outcome is in Jupiter masses, M_J. For 31 systems lacking either K or i, the planet masses from exoplanet.eu were adopted.

This procedure yielded Q′_WB and M_crit for a single MESA profile and, by repetition, their evolution along a complete stellar track. We applied this algorithm to each of the 100 random realisations of stellar parameters obtained from the isochrones posterior distributions. For each realisation, we selected the corresponding MESA model from our grid by rounding its M_* and [Fe/H] to two decimal places. The few realisations falling outside the grid were omitted.

The MESA models evolve with adaptive time-steps, which can introduce differences between tracks. To homogenise the output, we interpolated each model over an evenly spaced grid of 1000 points spanning the ±1σ age range of the system. For each of the 100 realisations, Q′_WB values were stored if wave breaking was expected. From this distribution, the median and 15.9 and 84.1 percentile values were adopted as the central value of Q′_WB and its 1σ uncertainties.

The same procedure was applied to M_crit To quantify the likelihood of wave breaking, we defined the fraction, $$f_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\frac{N(M_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}<M_{\mathrm{p}})}{N_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}},$$(4)

where N(M_crit < M_p) is the number of interpolated points for which wave breaking is expected, and N_total is the total number of sampled points, including those when M_crit is undefined due to a convective stellar core. The resulting f_crit, expressed as a percentage, represents the overall probability of wave breaking in a given system within its 1 σ age interval.

This approach naturally accounts for the three principal configurations of the stellar interior that determine the occurrence of wave breaking. When the stellar core is convective, wave breaking cannot occur. If the core is radiative but M_crit exceeds M_p, the internal gravity waves do not reach sufficient amplitude to break. Conversely, when the core is radiative and M_crit falls below M_p, wave breaking is expected to take place.

For systems in which planets are insufficiently massive to trigger wave breaking during radiative-core phases (i.e. f_crit = 0), we nevertheless computed the corresponding Q′_WB values following the same procedure but without imposing the M_crit criterion. As discussed by Barker (2020), such values represent the formal wave-breaking dissipation efficiency and remain useful for comparison with observational constraints, even if wave breaking is not expected to occur under the current system parameters.

We also calculated the value of the change of orbital period Ṗ using $$\dot{P}=-\frac{27}{2}\frac{\pi }{Q_{\mathrm{W}\mathrm{B}}^{\prime }}\frac{M_{\mathrm{P}}}{M_{\star }}{\left(\frac{a}{R_{\star }}\right)}^{-5}.$$(5)

The results of the tidal dissipation modelling for all systems are summarised in Table A.3. Additionally, the results are presented for each planet in a dedicated multi-panel figure, such as Fig. 1, shown as an illustrative example. Each figure is divided into four columns. The first column shows a histogram of the Q′_WB values within the 1 σ age interval, colour-coded as follows: black for convective phases (wave-breaking cannot occur), grey for radiative phases without wave breaking (M_crit is greater than M_p), and blue for radiative phases where wave breaking occurs. The orange solid and dotted lines indicate the median and ±1σ ranges of Q′_WB, respectively. The second column displays the evolution of Q′_WB, with the bottom panel showing the full evolutionary track and the top panel zooming in around the 1 σ age interval (white background). Each of the 100 realisations is shown with the same colour scheme. The final Q′_WB value with its uncertainties is given at the top. The third column presents the evolution of M_crit, plotted analogously to Q′_WB, but with red lines marking radiative phases where M_crit < M_p. The corresponding probability of wave breaking, f_crit, is given above the panel, and the planetary mass M_p is shown as a purple dashed line. The fourth column illustrates the evolution of P, using green lines to denote radiative phases where wave breaking is expected, with the final P value and its uncertainties provided at the top.

Fig. 1 Predictions of tidal interactions in the WASP-4 system, shown as an illustrative example of our modelling outcome. The first column shows the histogram of Q′WB values. The second column shows the evolution of Q′WB, with the top panel zoomed in around the 1σ age estimate. Each line corresponds to one realisation of stellar parameters. Black lines indicate phases when the host’s core is convective, grey lines mark phases when the host’s core is radiative but wave breaking is not expected, and blue lines denote phases when the host’s core is radiative and wave breaking occurs. The estimated Q′WB with uncertainties is shown above the column. The third column shows the evolution of Mcrit, with red lines indicating phases when the host’s core is radiative and wave breaking occurs. The estimated probability of wave breaking is shown above the column. The fourth column presents the evolution of Ṗ, with green lines marking phases when the host s core is radiative and wave breaking occurs. The estimated value of Ṗ with uncertainties is displayed above the column.

3.4 Transit timing diagnostics

In addition to the parameters described above, we calculated the quantity T_shift, which approximates the cumulative shift in the central transit time since the system’s discovery. It is derived using $$T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}=\frac{27}{4}\frac{\pi }{Q_{\mathrm{W}\mathrm{B}}^{\prime }}\frac{M_{\mathrm{p}}}{M_{\star }}{\left(\frac{a}{R_{\star }}\right)}^{-5}{P}^{-1}{(T_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{e}}-T_{\mathrm{d}\mathrm{s}\mathrm{c}})}^2,$$(6)

where T_date is the date of new observations and T_dsc is the date of the earliest available data. In the Table A.3, we provide T_shift values calculated for the year 2025.

We also estimated the epoch at which the cumulative timing shift reaches 5 minutes, using $$T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}=\sqrt{\frac{5\mathrm{m}\mathrm{i}\mathrm{n}}{T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}}}(2025-T_{\mathrm{d}\mathrm{s}\mathrm{c}})+T_{\mathrm{d}\mathrm{s}\mathrm{c}},$$(7)

where $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$ is the cumulative transit timing shift calculated for the year 2025. The $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ estimate is particularly useful for recently discovered systems, for which even strong tidal interactions may not yet have produced an observable shift, due to the short timespan of available observations. The threshold of 5 minutes was adopted following the observed TTV of WASP-12 b (see Fig. 2 in Maciejewski et al. 2016), where a cumulative shift of approximately this magnitude was measured.

The final estimates of $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$ and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$, together with their uncertainties, were derived in the same manner as Q′_WB. The resulting values for all analysed systems are listed in Table A.3. We emphasise that the detectability of TTV varies between systems, depending on factors such as stellar brightness, transit duration and depth, photometric stability of the host, and the precision of available light curves. We therefore leave the assessment of detectability to the reader, providing $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$ and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ as reference values.

3.5 Orbital evolution

The orbital evolution of each system was computed by integrating the time-dependent Ṗ, recalculated at each step together with Q′_WB and M_crit from the evolving stellar and orbital parameters. The procedure follows the method outlined in Sect. 3.3, extended to account for the coupled evolution of the star and orbit.

The integration time step was dynamically adjusted to ensure numerical stability while minimising computational cost. The integration was carried out until the planet’s orbital semi-major axis decreased below the stellar radius - at which point the planet was considered engulfed - or the planet’s disruption at Roche limit (Faber et al. 2005; Ford & Rasio 2006), which was calculated based on the equation presented in Bonomo et al. (2017): $$a_{\mathrm{R}\mathrm{o}\mathrm{c}\mathrm{h}\mathrm{e}}=2.16R_{\mathrm{p}}{\left(\frac{M_{\star }}{M_{\mathrm{p}}}\right)}^{1/3}.$$(8)

The semi-major axis was calculated from P assuming a circular orbit, as described in Sect. 3.3. The integration was also terminated if the system’s age exceeds the upper age limit of our MESA grid (18 Gyr).

If the orbital evolution led to planetary destruction, either by engulfment or by disruption at the Roche limit, the host star’s age and EEP at that moment were recorded as the age and EEP of destruction. After completing the integrations for all 100 realisations, the resulting distributions of ages and EEPs of destruction were used to compute the median and lσ confidence intervals for both quantities. These results are listed in Table A.3 and illustrated for a representative case in Fig. 2.

Fig. 2 Orbital evolution of the planet MASCARA-3 A b, shown as an illustrative example of our modelling outcome. Each line in the panels corresponds to one of the 100 realisations used for this system. Colours indicate evolutionary regimes: red, brown, and blue show Mcrit, P, and Q′WB, respectively, during phases with active wave breaking. Grey denotes radiative cores where the wave-breaking criterion is not fulfilled. Black marks convective-core phases. The black dotted and dash-dotted vertical lines indicate the positions of the intermediateage main sequence (IAMS; EEP 353) and terminal-age main sequence (TAMS; EEP 454), respectively. The planet’s mass is shown as a purple dashed line in the top two panels. Orange solid and dashed lines trace the tracks ending in the median and 1σ ranges of the age or EEP of destruction. Orange circles mark the starting points of these tracks, while orange crosses denote the moment of planetary destruction.

4 Discussion

To place our results in a broader context, we examined the population-level behaviour of tidal dissipation driven by the WB mechanism across our sample, identified the emerging trends, and discussed their physical origins. We also assessed the implications for orbital evolution and the prospects for observational detection. While the full matrix of marginal and joint posterior distributions for all key parameters is presented in Fig. B.1 (Appendix B), the discussion below is illustrated with appropriately selected diagrams in Figs. 3-8.

4.1 Population statistics

With the full set of tidal-evolution calculations in hand, we examined the global behaviour of the sample. Figure 3 presents the joint distribution of the host-star evolutionary phase, expressed through EEP, and the present-day orbital period P, accompanied by marginal histograms. The use of EEP enables a uniform comparison across stars with different masses and metallicities, providing a physically consistent measure of evolutionary stage. The colour scale encodes the probability of wave breaking f_crit. Grey points correspond to systems with radiative cores for which wave breaking is not expected because M_crit > M_p throughout the 1σ age range; the plotted values therefore represent formal Q′_WB estimates. Black points mark systems hosted by stars maintaining a convective core over this interval and therefore cannot undergo wave breaking. The purple symbol marks the WASP-12 subgiant scenario (Sect. 4.5), whose artificially enhanced dissipation is consistent with its well-documented rapid orbital decay. In our framework, it has f_crit = 100% and serves as a reference case for a system firmly in the wave-breaking regime.

A strong dependence on the stellar evolutionary state is evident. Pre-IAMS systems (squares) show virtually no wave breaking, consistent with their structural inability to support gravity waves of sufficient amplitude. The probability of wave breaking rises markedly between the IAMS and TAMS (circles), and remains high for stars beyond the TAMS (triangles), provided that the planet survives to such phases. This behaviour is visible not only in the distribution of data points but also in the EEP histogram, where the fraction of systems with high f_crit increases progressively with increasing EEP. The orbital-period histogram, in turn, reveals the well-known period-distribution structure of close-in giants, with a pronounced pile-up at P ≈ 34 d (Udry et al. 2003; Yee & Winn 2023), common to all evolutionary groups represented in the diagram.

Interestingly, the systems with the very shortest orbital periods - where the strongest tidal interactions might be expected -generally exhibit low f_crit. Among planets with P < 1.0 d, the highest wave-breaking probability is f_crit = 23%, clearly indicating that a short orbital period is not a sufficient condition for efficient wave breaking. This absence of systems undergoing wave-breaking extends down to P < 1.5 d, as shown in the P histogram. The only case in this regime that attains f_crit = 100% is the WASP-12 subgiant scenario.

This behaviour reflects a strong evolutionary filtering effect: very short-period planets are rarely found around evolved hosts. For P < 1.5 d, the WASP-12 subgiant scenario is by far the most evolved case ($\mathrm{E}\mathrm{E}\mathrm{P}={426.3}_{-1.6}^{+1.7}$), whereas no other system in this group exceeds EEP = 400. This finding suggests that wave breaking efficiently removes short-period planets once their hosts evolve beyond the IAMS. The grey lines show evolution of radii for 0.8-1.5 M_⊙ hosts, transformed into the engulfment orbital periods. They typically do not exceed 0.4 d between the IAMS and TAMS, showing that the lack of short-period planets after the IAMS cannot be explained by stellar evolution alone.

The same pattern is visible in Fig. 4, which shows the joint distribution of systems in the Q′_WB-|Ṗ| plane. As expected, the strongest orbital decay rates correspond to the lowest predicted values of Q′_WB, reflecting the tight theoretical coupling between these two quantities. However, most of the systems that occupy this high-|Ṗ|-low-Q′_WB region have hosts that are still residing prior to the IAMS and the M_crit condition has not been fulfilled for them (or is met only with very low probability) across their 1σ age range. In such systems, wave breaking is therefore not expected, and consequently, no rapid orbital decay driven by this mechanism should occur - despite their apparently favourable locations in the diagram. These cases include WASP-19, WASP-43, and HATS-18, which are highlighted in Fig. 4 and discussed further in Sect. 4.3. The WASP-12 subgiant scenario remains the principal exception, standing out as the only system that combines low Q′_WB, large |P|, and a high probability of wave breaking.

The Q′_WB histogram displays a single broad peak centred at ∼ 3 × 10⁶, a behaviour shared across systems independent of their wave-breaking probability. In terms of orbital decay, these values correspond to |P| ~ 10⁻³−10⁻² ms yr⁻¹. Systems with low f_crit predominantly occupy the lower end of this range, while those with high f_crit cluster toward the upper end.

Figure 5 shows the distribution of the sample systems in the EEP-host-mass plane. As evident from both the main panel and the mass-side histogram, two regimes emerge, separated at M_* ≈ 1.2 M_⊙. Below this threshold, hosts begin their evolution with radiative cores, and wave breaking typically becomes efficient only near or after the IAMS. At earlier phases, WB is suppressed because the M_crit condition is not satisfied. Above M_* ≈ 1.2 M_⊙, stars start out with convective cores and can enter the wave-breaking regime only at more advanced stages, typically between the IAMS and TAMS rather than near the IAMS itself. Both populations converge around EEP ≈400, where f_crit becomes moderate or high across the full mass range.

A notable feature is the near-complete absence of systems with M_* > 1.2 M_⊙ between EEP = 400 and the TAMS. We identify this gap as only apparent, caused by the rapid evolution of such stars. In this parameter range, the timestep per EEP decreases by nearly two orders of magnitude, greatly reducing the probability of detecting a system within that narrow evolutionary interval. This finer resolution in EEP also manifests as larger uncertainties in EEP for systems in this region.

The lack of post-IAMS planets around the lowest-mass hosts (M_* < 0.8 M_⊙) can instead be attributed to their longer evolutionary timescales. Stars in this mass range evolve so slowly that even the oldest systems in this group, with ages comparable to the age of the Universe, have not yet reached the IAMS. Consequently, the absence of post-IAMS planets in this region does not reflect a physical clearing mechanism, but rather the limited evolutionary progress such stars can achieve within the Hubble time.

Fig. 3 Joint plot of EEP and orbital period, P, for the studied systems. The colour of the points, from blue to red, corresponds to the probability of wave breaking, fcrit, as indicated by the colour bar in the top-right corner. Grey points denote systems with host stars with radiative cores and the wave breaking criterion being unfulfilled, while black points indicate systems with convective core stars. The purple point is the WASP-12 subgiant scenario, as described in Sect. 4.5. The shape of the points indicates the evolutionary stage of the host: squares for systems before the IAMS, circles for those between the IAMS and TAMS, and triangles for those beyond the TAMS. The dotted and dash-dotted vertical lines mark IAMS and TAMS evolutionary phases, respectively. The colour scale in the histograms corresponds to the colour of the points. The grey lines show the engulfment period for stars with masses between 0.8 and 1.5 M⊙.

Fig. 4 Relation between the tidal quality factor Q′WB on the orbital period derivative, P, due to wave breaking for the studied systems. Colours and shapes of the points follow the same scheme as in Fig. 3.

Fig. 5 Joint plot of the host mass and EEP for the studied systems, with symbols and colours defined as in Fig. 3.

4.2 Population-level orbital evolution

The evolution of the orbital period for all systems in our sample is shown in Fig. 6. Each coloured line traces the orbital trajectory of an individual planet, with the colour indicating its eventual fate. Olive lines correspond to planets that are ultimately engulfed by their host stars, while cyan lines mark systems in which planets first reach the Roche limit and are disrupted. In total, 235 planets (≈42.7%) are destroyed between the IAMS and TAMS, and a comparable number, 227 (≈41.3%), are destroyed after the TAMS through one of these two channels. No planets are destroyed before their hosts reach the IAMS, consistent with the general inefficiency of tidal dissipation at earlier evolutionary phases. A further 89 systems (≈16.2%) do not reach destruction within the evolutionary range covered by our MESA grid; their tracks are shown in brown.

Among the systems that are destroyed in our simulations, 79 planets (≈14.4% of the entire sample) undergo direct engulf-ment, while 383 (≈69.6%) are tidally disrupted at the Roche limit. Disruption is particularly dominant in the pre-TAMS regime: 219 planets are disrupted before the TAMS, whereas only 16 are engulfed. Between the TAMS and EEP ≈500, engulf-ment becomes increasingly important, with 36 planets engulfed and 164 disrupted. Beyond EEP ≈500, engulfment becomes the dominant channel, reflecting the rapidly rising stellar radii in post-main sequence evolution. In this regime, tidal decay driven by wave breaking still contributes to orbital shrinkage, but stellar expansion increasingly governs the final outcome. The eventual fate of planets from this group is therefore set by the combined action of tidal dissipation and the accelerated growth of the host-star envelope.

Orbital decay driven by wave breaking becomes extremely rapid for systems with P ≲ 2 d whose hosts have evolved beyond the IAMS. As shown in Fig. 6, systems in this regime undergo orbital collapse on timescales far shorter than stellar evolution. This is consistent with the trend noted in Sect. 4.1, where no very-short-period planets are found around hosts beyond EEP ≈400, with the sole exception of the WASP-12 subgiant scenario. Thus, the short-period planets observed today most likely orbit hosts that remain before EEP ≈350, a phase in which wave breaking is not yet operating.

The efficient destruction at short orbital periods is further clarified by examining a relation between the EEP of destruction EEP_dest and the present-day orbital period, shown in Fig. 7. This diagnostic directly links current system architecture to its expected evolutionary fate and provides a population-level view complementary to the phase-dependent orbital tracks shown in Fig. 6.

Figure 7 reveals a clear trend: systems with longer presentday orbital periods are destroyed at systematically later evolutionary stages, with the shortest-period planets being removed shortly after their hosts cross the IAMS. The side histogram along P shows the marginal distributions for planets destroyed before (green) and after (pink) the TAMS. Both distributions drop sharply at their point of intersection, near P ≈ 3.5 d, marking a well-defined division between the two sub-populations.

The marginal distribution along EEP_dest exhibits three prominent peaks, each corresponding to a distinct physical regime. The first peak, at EEP_dest ≈ 405, traces short-period systems destroyed before the TAMS. Nearly all planets in this group (≈88 %) have P ≲ 3.5 d, consistent with the sharp boundary seen in the period histogram. Systems with P ≳ 3.5 d typically reach the TAMS and may survive through it or be destroyed around that phase, producing the second peak around EEP_dest ≈ 455. At wider separations (P ≳ 5-6 d), planets begin to avoid rapid tidal inspiral and instead follow an engulfment sequence in which destruction is dictated primarily by stellar expansion. This behaviour gives rise to the third peak, located near EEP_dest ≈ 525, beyond which engulfment dominates and the EEP of destruction increases systematically with orbital period.

Predictions arising from our framework may manifest in population-level trends by modulating the occurrence rate of hot Jupiters across different evolutionary phases of their host stars. The sub-sample of systems whose hosts remain before the IAMS can reasonably be regarded as unaffected by the wave-breaking mechanism, thereby providing a suitable reference population against which to assess relative occurrence rates. Within this group, comprising 227 systems, we examined the distribution of EEP_dest and identified 115 systems predicted to be destroyed before the TAMS, 36 between the TAMS and EEP = 465 (the boundary between subgiant and red giant stars defined by Bryant & Van Eylen 2025), and a further 36 disrupted by EEP = 530 during the early red giant phase.

These numbers were translated into relative occurrence rates across these three evolutionary intervals. Taking the subgiant occurrence rate as the reference level, we find that the occurrence rate for main sequence hosts is higher by a factor of 1.80 ± 0.13, whereas that for early red giant hosts is lower by a factor of 1.6 ± 0.3. Bryant & Van Eylen (2025) reported a decline in the occurrence rate of giant planets with advancing stellar evolution, comparing post-main sequence hosts with literature determinations for the main sequence. A qualitative comparison with their Fig. 10 indicates that our inferred trend is broadly consistent with previous findings. We note, however, that occurrence-rate estimates for main sequence stars exhibit substantial scatter, and our relative occurrence rate for early red giants is higher, potentially hinting at a modest underestimation of the efficiency of planet destruction in this evolutionary stage. This could be explained by other tidal interactions, including convective damping of equilibrium tides, which becomes more vigorous in early red giant stars (Duguid et al. 2020).

Fig. 6 Same as Fig. 3, but including the evolutionary tracks of planetary orbits for all systems in our sample. The coloured solid lines trace the orbital evolution driven by wave breaking: olive lines indicate systems ending in engulfment, cyan lines show those disrupted at the Roche limit, and brown lines denote tracks that terminate early due to the limits of the MESA grid.

Fig. 7 Joint plot of EEP of destruction and the present-day orbital period. The point types and colour coding follow those in Fig. 3. The right-hand panel shows two separate histograms: one for systems destroyed before the TAMS (green), one for systems destroyed after the TAMS (pink). The overlapping region of the distributions is shown with semi-transparent shading. A clear separation between the two populations emerges at P ≈ 3.5 d.

Fig. 8 Predicted year at which the cumulative shift in central transit times reaches 5 minutes, plotted as a function of the absolute value of the orbital period change |P| due to wave breaking. The horizontal black line indicates the year 2035.

4.3 Implications for follow-up observations

We assessed which of the studied systems offer the most promising opportunities for detecting orbital decay through TTVs. To this end, we analysed the parameters $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$ and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$, which, together with the predicted f_crlt values, provide a useful basis for planning and prioritising follow-up campaigns. In the most favourable cases, $f_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=100\mathrm{\%}$ and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ is close to the current epoch, or lies within a timescale accessible to feasible observations. For systems benefiting from higher photometric precision or long-term monitoring, smaller T_shift values may also yield detectable trends.

The calculated $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$ values depend on the available observational baseline for each system and were derived using Eq. (6) as (2025 − T_dsc). The practical ability to detect shifts in transit times may, however, vary depending on the quality, quantity, and temporal coverage of the available data, and is therefore beyond the scope of this work. For uniformity, we adopted the discovery year (i.e. the year the planet was announced) as a reference epoch across the entire sample. In specific cases, earlier photometric data suitable for TTV analysis may exist, leading to slightly higher effective T_shift values, particularly for recently published systems. Conversely, the effective baseline may be shorter if the earliest observations are of insufficient quality or unavailable. Since T_shift scales quadratically with the duration of the observational baseline, adjustments to account for the true first-observation epoch are straightforward to apply.

Having established the framework for assessing detectability through T_shift, we identified systems in which a measurable shift in transit central times could occur in the near future. Figure 8 presents a diagnostic diagram of $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ against |Ṗ|, restricted to the coming decades. This plot provides a concise overview of the expected detectability of tidal orbital decay driven by wave breaking. Systems located towards the lower-right corner of the diagram, characterised by large |Ṗ| and small $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$, represent the candidates which can produce measurable TTVs within the next decade, provided wave breaking is active. For systems with f _crit = 0% the presented values of T_shift are calculated using the formal Q′_WB estimate. We stress, however, that if wave breaking is indeed inactive, no measurable 5-minute transit-time shift is expected on practical observational timescales.

Among the systems shown in Fig. 8, we identified nine planets with the highest |Ṗ| values and the earliest predicted $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ before 2035: HAT-P-23 b, HATS-18 b, HIP 65 A b, NGTS-10 A b, TOI-1937 b, WASP-4 b, WASP-19 A b, WASP-43 b, and WASP-103 b. The individual characteristics and current observational constraints for these planets are discussed below. Their relevant parameters, including $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ and f_crit, are summarised in the upper part of Table 1.

HAT-P-23 b. This system hosts a hot Jupiter of 2.1 M_J orbiting a G-type star every 1.21 d (Bakos et al. 2011). It was investigated using the transit timing method by Maciejewski et al. (2022), who combined over a decade of photometric observations. Their analysis revealed no statistically significant trend in transit times, resulting in a lower limit on Q′ of (2.76 ± 0.21) × 10⁶. Our modelling yields an effective value of $Q_{\mathrm{W}\mathrm{B}}^{\prime }={6.7}_{-1.2}^{+1.8}\times {10}^5$, but predicts f_crit = 0%, as the host star remains before the IAMS at EEP ≈330, below the evolutionary stage at which efficient wave breaking typically develops. Wave breaking is therefore not expected for the present stellar model. If it were operating at the level implied by Q′_WB, the resulting orbital decay would produce a readily detectable TTV signal. The absence of a measurable P is thus fully consistent with the star s evolutionary state and the tidal regime inferred from our models.

HATS-18 b. This system consists of a G-type host star and a hot Jupiter with a mass of approximately 2 M_J on a 0.84-day orbit, and was originally highlighted as a promising target for investigating tidal interactions (Penev et al. 2016). Transit-timing studies were conducted by Patra et al. (2020), Southworth et al. (2022), and Maciejewski et al. (2024). The first lacked sufficient observational coverage for meaningful constraints, while the latter two found no detectable trend in transit times, placing lower limits on $Q_{\star }^{\prime }of{1.29}_{-0.11}^{+0.12}\times {10}^5$ and ${3.5}_{-0.7}^{+0.5}\times {10}^5$, respectively.

These limits exceed our predicted value of ${7.13}_{-0.68}^{+0.69}\times {10}^4$, indicating that wave breaking is not currently active in this system. If it were operating at the level implied by Q′_WB, a 5-minute shift in the transit central times would have been expected by 2022, which would be readily detectable in TTV monitoring.

Our predicted probability of wave breaking, f_crit = 20.6%, naturally explains this non-detection and suggests that HATS-18 more likely lies towards the younger side of its 1σ age interval, between 2.2 and 4.3 Gyr. These findings imply that alternative tidal mechanisms should be considered as potential contributors to the observed spin-up of the host (Penev et al. 2016, 2018). Improved constraints on the stellar parameters would help to clarify the origin of this spin-up and the role of tides in the system.

HIP 65 A b. This system consists of a G-type star orbited by a hot Jupiter of 3.2 M_J on a short 0.98-day orbit (Nielsen et al. 2020). Recent studies by Maciejewski et al. (2024) and Griffiths et al. (2026) analysed over a decade-long transit timing observations, finding no significant deviation from a constant orbital period and establishing a lower limit on the stellar modified tidal quality factor of $Q_{\star }^{\prime }>{7.6}_{-0.7}^{+0.8}\times {10}^4$ (Maciejewski et al. 2024). Our modelling predicts that the host star, still before the IAMS (EEP ≈ 340), remains below the evolutionary stage at which internal gravity waves can reach amplitudes sufficient for breaking. If, however, wave breaking were to occur, it would imply an effective Q′_WB ≈ 1.4 × 10⁵, which remains above the threshold detectable with current transit timing precision.

NGTS-10 A b. This system hosts a massive hot Jupiter (M_p = 2.16 M_J) on one of the tightest orbits in our sample (P = 0.77 d) around a K-type star (McCormac et al. 2020). The host star is still before the IAMS and our modelling predicts no wave breaking. If wave breaking were somehow active, it would correspond to an effective $Q_{\mathrm{W}\mathrm{B}}^{\prime }={9.0}_{-0.6}^{+0.8}\times {10}^4$, placing the system among the few in our sample with Q′_WB below 10⁵ (see Fig. 4). In that scenario, the cumulative 5-minute shift in transit central times would occur as early as 2027, making the signal readily detectable through a dedicated TTV campaign. Although Ivshina & Winn (2022) included NGTS-10 A b in their survey, their analysis was based on only a single light curve. More recently, Griffiths et al. (2026) derived a lower limit of $Q^{\prime }>{5.4}_{-2.0}^{+3.0}\times {10}^4$, closely consistent with our prediction. Continued monitoring of this system could therefore provide a stringent constraint on the onset of tidal dissipation as the star approaches the IAMS.

TOI-1937 b. This system hosts a G-type star orbited by a massive hot Jupiter (M_p ≈ 2 M_J) on a short-period orbit (P = 0.95 d) (Yee et al. 2023). It belongs to the small group of systems in our sample with P < 1 d and a non-zero probability of wave breaking. The predicted likelihood remains modest, with f_crit = 23.1%, and the expected cumulative 5-minute shift in transit central times is projected to occur around 2030. Despite a relatively large predicted period change of |Ṗ| = 21 ms yr⁻¹, no such shift has yet been observed. A recent TTV analysis by Jankowski et al. (2025) reported an upper limit of |P| < 90 ms yr⁻¹, which does not rule out ongoing wave breaking. Continued high-precision photometric monitoring will therefore be crucial to more tightly constrain tidal dissipation in this system.

WASP-4 b. One of the earliest exoplanets discovered by the WASP programme, this hot Jupiter (M_p ≈ 1.2 M_J) orbits its G-type host star with a period of 1.34 d (Wilson et al. 2008). Deviations from a linear ephemeris were first reported by Bouma et al. (2019), prompting various interpretations, including acceleration towards the Earth (Bouma et al. 2020), the presence of an additional planet (Turner et al. 2022), and the Rømer effect (Harre & Smith 2023). More recently, Baştürk et al. (2025) proposed that orbital decay driven by tidal dissipation through wave breaking could explain the observed period shortening. Their theoretical estimate of Q′ ≈ (2-5) × 10⁵ aligns with our predicted $Q_{\mathrm{W}\mathrm{B}}^{\prime }={2.950}_{-0.043}^{+0.061}\times {10}^5$, although both remain higher than the observationally derived Q′ ≈ 8 × 10⁴. This moderate discrepancy has been attributed to a possible tension between the stellar interior structure and standard isochronal models, assuming that wave breaking were the sole mechanism driving inspiral.

More recent analysis by Winn & Stefánsson (2025) showed that the observed timing variations are more plausibly explained by the light travel-time effect induced by an additional planet in the system. Our models also indicate a relatively low probability of wave breaking (≈28%), implying that tidal dissipation through this mechanism is unlikely to be active at the current evolutionary stage of the host. Taken together, these results disfavour wave breaking as the driver of the observed period change in WASP-4 b.

WASP-19 A b. With one of the tightest orbits among known hot Jupiters (P = 0.78 d), WASP-19 A b orbits a G-type star (Hebb et al. 2010). This system has been extensively studied in the context of tidal interactions, yielding varied results (Ivshina & Winn 2022; Petrucci et al. 2020; Rosàrio et al. 2022; Adams et al. 2024; Maciejewski et al. 2024). Our models predict that wave breaking is not occurring, in agreement with the nondetections reported in the literature. The observational lower limits on Q′ from these studies - 2.51 × 10⁶ (Adams et al. 2024), 1.26 × 10⁶ (Rosàrio et al. 2022), 1.23 × 10⁶ (Petrucci et al. 2020), and 4.79 × 10⁶ (Maciejewski et al. 2024) - exceed our predicted Q′_WB = (6.92 ± 0.03) × 10⁴ by more than an order of magnitude, implying that wave breaking would already have produced a clearly detectable signal. Indeed, with $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}=2014.9\pm 0.2$, a measurable shortening of the orbital period would have been observed by now if wave breaking were active. A tentative detection of orbital decay reported by Ivshina & Winn (2022) (Ṗ = −3.54 ± 1.18 ms yr⁻¹) is an order of magnitude smaller than our predicted P = −38 ± 2 ms yr⁻¹ and, therefore, cannot be attributed to wave breaking.

WASP-43 b. This system consists of a K-type star hosting a 2 M_J hot Jupiter on a 0.81 d orbit (Hellier et al. 2011). Due to its extremely short period, WASP-43 b has been extensively studied in the context of tidal interactions. Several analyses have found no evidence of orbital decay, placing lower limits on the modified tidal quality factor of Q′ > 3.9 × 10⁵ (Alvarado et al. 2024) and Q′ > 1.15 × 10⁶ (Maciejewski et al. 2024). A combined orbital decay and apsidal motion scenario was proposed by Bernabò et al. (2025), though its interpretation remains uncertain. Our models predict that wave breaking is not currently operating in this system. If it were active, the expected |Ṗ| = 14 ms yr⁻¹ would produce a cumulative 10-minute shift in transit timing by 2025 - readily detectable with current observational capabilities. Any reported timing variations are therefore unlikely to be caused by wave breaking.

WASP-103 b. This ultra-hot Jupiter orbits an F-type host star (Gillon et al. 2014) and was identified by Barker (2020) as a particularly interesting system for testing tidal dissipation models. That study predicted that wave breaking should not occur owing to the presence of a convective core, but noted that the system’s relatively low Q′_WB could nonetheless be observationally constrained. Our models agree, yielding only a 0.5% probability of wave breaking, consistent with convective suppression. If the stellar core were radiative, a 5-minute cumulative shift in transit times would be expected around 2032. Recent analyses have found no statistically significant departure from a linear ephemeris (Alvarado et al. 2024; Maciejewski et al. 2022; Baştürk et al. 2022). The most stringent constraint to date, Q′ > 1.18 × 10⁶ (Alvarado et al. 2024), approaches the threshold required to test our predicted value of $Q_{\mathrm{W}\mathrm{B}}^{\prime }={1.31}_{-0.37}^{+0.62}\times {10}^6$.

Despite having the highest predicted strength of tidal interactions in our sample, all of these systems exhibit low probabilities of wave breaking, with f_crit values below 30% or consistent with zero. This suggests that wave breaking is unlikely to be the dominant tidal dissipation mechanism currently operating in any of these systems. The absence of significant P detections in recent transit-timing studies is therefore fully consistent with our theoretical predictions.

Beyond the systems discussed above, our models also identify a broader group of targets with weaker observational prospects. These planets are characterised by smaller absolute values of P and larger $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$, implying that any measurable transit-timing shifts driven by tidal dissipation would accumulate only over longer timescales. Among them, five systems stand out in Fig. 8 - HATS-9 b, TOI-1855 b, TOI-2337 b, WASP-135 b, and WASP-163 b - which exhibit f_crit values above ≈60% and predicted $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ occurring before 2060. These systems, together with their key parameters (f_crit, $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$, and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{2025}$), are listed in the lower part of Table 1. They represent the most promising candidates for detecting orbital decay due to wave breaking on decadal timescales.

HATS-9 b, TOI-1855 b, TOI-2337 b, WASP-135 b, and WASP-163 b. These systems host massive hot Jupiters on shortperiod orbits (P < 3 d) around G-type dwarfs or slightly evolved stars, whose internal structures are favourable for efficient tidal dissipation. HATS-9 b stands out with f_crit = 100% and a predicted $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ of ~ 2040. TOI-1855 b has f_crit ≈ 60% and is expected to reach a 5-minute shift around 2055. TOI-2337 b exhibits f_crit ≈ 96%, though with larger uncertainties, placing its $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ at ${2054}_{-17}^{+10}$. WASP-135 b shows f_crit ≈ 8/% and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}$ just before 2050. WASP-163 b, with a f_crit = 100%, reaches a shift of 5 minutes on a similar timescale, albeit with higher uncertainty. Although no deviations from a linear ephemeris have yet been reported for these systems (Ivshina & Winn 2022), their predicted decay rates of |Ṗ| ≈ 2-5 ms yr⁻¹ suggest that sustained, high-precision transit monitoring in the coming decades could reveal measurable orbital decay consistent with tidal dissipation through wave breaking.

Future space-based missions will offer new opportunities to test these predictions. Several systems from our sample lie within the planned PLATO (Rauer et al. 2014) fields and will benefit from precise asteroseismic characterisation of their host stars, thereby reducing the current uncertainties in stellar parameters and evolutionary stage. Moreover, departures from linear ephemerides may be detected well before the epoch at which the cumulative timing shift reaches 5 minutes, owing to the exceptional photometric precision of space-borne facilities.

One particularly interesting target is NGTS-3 A b, which has a 52.6% probability of wave breaking and a relatively low predicted Q′_WB ≃ (5.5 ± 1.2) × 10⁵. The expected epoch for a 5-minute timing shift is ${2062}_{-4.9}^{+4.0}$, but the host’s poorly constrained age (spanning about 8 Gyr at 1σ) leaves open the possibility of stronger tidal dissipation if the star turns out to be more evolved than currently inferred. High-precision, long-term monitoring of this and similar systems with PLATO and other photometric facilities will therefore be essential for probing tidal decay on decadal timescales and for testing the predictions of the wave-breaking framework.

4.4 Systems with reported orbital decay

Our results also provide a theoretical framework to interpret several systems in which orbital decay has been claimed or investigated through transit timing analyses. In several such cases, the reported departures from linear ephemerides cannot be reproduced by the internal structure of the host stars, which in our models prevents wave breaking from occurring.

A number of systems - KELT-1 b (Baştürk et al. 2023), KELT-16 b (Mancini et al. 2022), WASP-99 b (Ivshina & Winn 2022; Sodickson & Grunblatt 2025), HAT-P-32 b (Hagey et al. 2022), XO-3 b (Ivshina & Winn 2022), and WASP-32 b (Sun et al. 2023) - are predicted to have stars with fully convective cores in their current stage of evolution, thereby effectively suppressing the wave-breaking mechanism. Similarly, systems such as HAT-P-19 b, TrES-1 A b, and TrES-2 b, for which tentative orbital decay was reported by Hagey et al. (2022), exhibit f_crit = 0% in our models, confirming that tidal dissipation driven by wave breaking cannot account for the observed timing trends.

Two additional systems, HAT-P-51 b and HAT-P-53 b, have recently been reported by Yeh et al. (2024) to show apparent orbital decay. Our models predict wave-breaking probabilities of f_crit = 9.2% and 94%, respectively. However, the reported decay rates - P = −134 and −82 ms yr⁻¹ (derived from Table 2 of Yeh et al. 2024) - exceed our theoretical expectations of 8.2 × 10⁻⁴ and 0.2 ms yr⁻¹ by several orders of magnitude. If these departures from linear ephemerides are confirmed, they cannot be attributed to tidal dissipation via wave breaking. Alternative mechanisms, such as observational systematics or dynamical perturbations, should be considered.

4.5 Testing a low-mass host-star scenario for WASP-12

The WASP-12 system has long served as a benchmark for studies of tidal interactions in close-in giant planets. Its 1.47 M_J hot Jupiter orbits an F-type host star every 1.09 days (Hebb et al. 2009), and WASP-12b was the first exoplanet for which tidally driven orbital decay was detected (Maciejewski et al. 2016). As a result, the system provides a crucial test case for theories of tidal dissipation.

A persistent discrepancy exists between the stellar parameters inferred from standard isochrone fitting and those implied from tidal evolution constraints (B20). Spectroscopic parameters favour a 1.3-1.4 M_⊙ main sequence host star (Weinberg et al. 2017; Bailey & Goodman 2019). In this mass range, the star possesses a convective core, which suppresses wave breaking and leads to tidal dissipation efficiencies far too weak - by nearly an order of magnitude - to reproduce the empirically inferred Q′ ≈ 2 × 10⁵. In contrast, the observed decay rate can be matched if the host is instead a lower-mass (∼1.2 M_⊙) subgiant with core convection already ceased at this evolutionary phase (Bailey & Goodman 2019).

Our modelling leads to the same conclusion. Using the spectroscopic parameters from the SWEET-Cat catalogue, the isochrone fit yields a mass of 1.343 ± 0.008 M_⊙ and EEP of ${366.0}_{-4.6}^{+1.9}$, placing this post-IAMS star firmly in the regime with a convective core, where wave breaking cannot occur. Consequently, in the absence of wave breaking, the tidal dissipation is insufficient to account for the observed decay of WASP-12b.

To test whether our wave-breaking model can reproduce the observed decay for the subgiant stellar structure, we followed the approach of B20 and constructed an additional, deliberately constrained scenario, referred to as WASP-12 subgiant. Using the same spectroscopic inputs as in our default analysis, we restricted the isochrone fit to the mass interval 1.18-1.22 M_⊙ to enforce a radiative-core configuration. This artificially narrowed model was then processed through our pipeline to compute Q′_WB, f_crit, and the corresponding orbital decay parameters.

Being located at EEP = 426.3 ± 1.7, the star is already significantly evolved (though still prior to the TAMS) and attains f_crit = 100%, yielding Q′_WB = 1.65 ± 0.50 × 10⁵, $\dot{P}=-{35}_{-8}^{+16}$ ms yr⁻¹, and $T_{\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{f}\mathrm{t}}^{5\mathrm{m}\mathrm{i}\mathrm{n}}={2015.1}_{-1.2}^{+1.0}$. These values are in close agreement with observational estimates - Q′ = 1.8 × 10⁵, P = −29 ± 2 ms yr⁻¹ (Yee et al. 2020) - and with the reported 5-minute timing shift detected in 2016 (Maciejewski et al. 2016). Taken together, these results demonstrate that the low-mass, radiative-core scenario successfully reproduces all the key observational characteristics of WASP-12 b’s infall. This suggests that wave breaking can plausibly account for the observed orbital decay provided that the host star is indeed less massive and more evolved than inferred from standard isochronal fits.

The source of the tension between spectroscopic and tidal-inferred stellar properties for WASP-12 remains unidentified (Bailey & Goodman 2019) and lies beyond the scope of this paper. These uncertainties may reflect limitations in current stellar modelling for WASP-12 and highlight the need for caution when interpreting stellar parameters derived solely from isochrone fits.

5 Conclusions

We developed a framework for modelling the inspiral of hot Jupiters driven by the WB mechanism, following the formalism of B20. This approach yields quantitative predictions for presentday orbital-period decay and enables us to project the subsequent dynamical evolution of planetary systems. While WB provides the most efficient dissipation of dynamical tides, complementary mechanisms, including weakly non-linear wave-wave interactions in radiative interiors (Weinberg et al. 2024) and inertial waves in convective envelopes (B20), remain necessary to achieve a comprehensive understanding of tidal evolution.

Applying this pipeline to a sample of 550 hot-Jupiter systems, we homogeneously redetermined the fundamental stellar properties, particularly mass and age, which critically influence tidal dissipation. Uncertainties in these parameters were fully propagated through the calculation of all derived quantities, including the tidal quality factor. Each system was assigned a probability that the WB mechanism is currently operating.

Our results indicate that WB becomes effective in stars with M_* < 1.2 M_⊙ predominantly after they cross the IAMS, while in more massive stars it begins during the evolutionary interval between the IAMS and the TAMS. Therefore, the pre-IAMS population of hot Jupiters is unperturbed by wave breaking. Within this framework, 43% of hot Jupiters (primarily those with orbital periods shorter than ∼3.5 d) are predicted to inspiral during the main sequence. A further 41% inspiral during post-main sequence evolution as their hosts become subgiants or early red giants, while the remaining 16% survive throughout the evolutionary stages covered by our models, although their eventual engulfment at later phases is inevitable.

Tidal disruption at the Roche limit emerges as the dominant destruction channel for most systems, with direct engulfment becoming the primary pathway only at more advanced stellar evolutionary stages, particularly for planets with orbital periods exceeding roughly 5-6 d. Our predictions for main sequence destruction provide a physical interpretation for the observed trend that hot Jupiters preferentially orbit younger stars (Schlaufman & Winn 2013; Hamer & Schlaufman 2019). The distribution of predicted tidal quality factors peaks between 10⁶ and 10⁷, in agreement with these population-level inferences.

Systems with the shortest orbital periods, which could, in principle, experience the strongest tidal dissipation (and therefore the fastest inspiral rates) are unlikely to undergo WB, leaving their planets on dynamically stable orbits. Instead, the most rapidly inspiralling systems with a high WB probability are predicted to display measurable orbital-period shortening within the next few decades. Although such predictions offer limited prospects for immediate observational confirmation (unless additional WASP-12-like systems are identified), the demographic imprint of the WB mechanism on hot-Jupiter occurrence rates is expected to be observable, with the first indications already emerging in current occurrence-rate estimates.

Data availability

Tables A.1, A.2, and A.3 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A60. The result figures for all systems studied in this paper are publicly available in the Zenodo repository https://doi.org/10.5281/zenodo.18772779.

Acknowledgements

We are in debt to Adrian Barker for valuable discussions on tidal dissipation calculations during the early stages of this research. JG acknowledges the financial support from the National Science Centre, Poland, through grant no. 2025/57/N/ST9/02205. This research has made use of the SIMBAD database (Wenger et al. 2000), operated at CDS, Strasbourg, France, and of data and tools provided by the Exoplanet Encyclopaedia (exoplanet.eu). We also made use of the NumPy (Harris et al. 2020) and SciPy (Virtanen et al. 2020) scientific libraries.

References

Appendix A Input parameters, isochrone fits, and tidal modelling results

Appendix B Correlation plots for the modelled sample

Fig. B.1 Mosaic of selected parameters for the studied systems, illustrating correlations between stellar, planetary, and tidal properties derived from our modelling. The dotted and dash-dotted vertical lines on EEP and EEPdest plots mark IAMS and TAMS evolutionary phases, respectively.

Appendix C MESA settings

We used MESA version r23.05.1. We used a custom pro-file_columns.list file to reduce calculation time and disk space usage, in which we only saved mass, logR, logRho and logP. Our inlist_project file was as follows:

&star_job
    create_pre_main_sequence_model = .true.
    set_uniform_initial_composition = .true.
    initial_h1 = H1
    initial_h2 = H2
    initial_he3 = He3
    initial_he4 = He4
    initial_zfracs = 6
    change_net = .true.
    new_net_name = 'pp_and_cno_extras.net'
    change_initial_net = .true.
&kap
    use_Type2_opacities = .false.
    kap_file_prefix = 'a09'
    kap_lowT_prefix = 'lowT_fa05_a09p'
    kap_CO_prefix = 'a09_co'
&controls
    initial_mass = M
    initial_z = Z
    MLT_option = 'Henyey'
    delta_lg_XH_cntr_max = −1
    mixing_length_alpha = 1.82
    mesh_delta_coeff = 0.7
    max_age = 1.8e10
    history_interval = 1
    profile_interval = 1
    max_num_profile_models = 500
    max_years_for_timestep = 2.5d8
    when_to_stop_rtol = 1d-6
    when_to_stop_atol = 1d-6
    max_model_number = 500
    stop_at_phase_He_Burn = .true.

M and Z are mass and metallicity of the given model in the grid, with Z calculated from [Fe/H] using: Z = Z_proto × 10^FeH with the base of Z_proto = 0.0142 taken from MIST. H1, H2, He3, He4 are also calculated for each model, changing with metallic-ity (based on Eqs. 1-3 from Choi et al. 2016). Some models had errors, which we fixed by increasing the mesh grid by lowering the mesh_delta_coeff to 0.5. For especially stubborn cases, we increased the grid further with max_dq = 1d-3.

All Tables

All Figures

Fig. 1 Predictions of tidal interactions in the WASP-4 system, shown as an illustrative example of our modelling outcome. The first column shows the histogram of Q′WB values. The second column shows the evolution of Q′WB, with the top panel zoomed in around the 1σ age estimate. Each line corresponds to one realisation of stellar parameters. Black lines indicate phases when the host’s core is convective, grey lines mark phases when the host’s core is radiative but wave breaking is not expected, and blue lines denote phases when the host’s core is radiative and wave breaking occurs. The estimated Q′WB with uncertainties is shown above the column. The third column shows the evolution of Mcrit, with red lines indicating phases when the host’s core is radiative and wave breaking occurs. The estimated probability of wave breaking is shown above the column. The fourth column presents the evolution of Ṗ, with green lines marking phases when the host s core is radiative and wave breaking occurs. The estimated value of Ṗ with uncertainties is displayed above the column.

In the text

Fig. 2 Orbital evolution of the planet MASCARA-3 A b, shown as an illustrative example of our modelling outcome. Each line in the panels corresponds to one of the 100 realisations used for this system. Colours indicate evolutionary regimes: red, brown, and blue show Mcrit, P, and Q′WB, respectively, during phases with active wave breaking. Grey denotes radiative cores where the wave-breaking criterion is not fulfilled. Black marks convective-core phases. The black dotted and dash-dotted vertical lines indicate the positions of the intermediateage main sequence (IAMS; EEP 353) and terminal-age main sequence (TAMS; EEP 454), respectively. The planet’s mass is shown as a purple dashed line in the top two panels. Orange solid and dashed lines trace the tracks ending in the median and 1σ ranges of the age or EEP of destruction. Orange circles mark the starting points of these tracks, while orange crosses denote the moment of planetary destruction.

In the text

Fig. 3 Joint plot of EEP and orbital period, P, for the studied systems. The colour of the points, from blue to red, corresponds to the probability of wave breaking, fcrit, as indicated by the colour bar in the top-right corner. Grey points denote systems with host stars with radiative cores and the wave breaking criterion being unfulfilled, while black points indicate systems with convective core stars. The purple point is the WASP-12 subgiant scenario, as described in Sect. 4.5. The shape of the points indicates the evolutionary stage of the host: squares for systems before the IAMS, circles for those between the IAMS and TAMS, and triangles for those beyond the TAMS. The dotted and dash-dotted vertical lines mark IAMS and TAMS evolutionary phases, respectively. The colour scale in the histograms corresponds to the colour of the points. The grey lines show the engulfment period for stars with masses between 0.8 and 1.5 M⊙.

In the text

	Fig. 4 Relation between the tidal quality factor Q′WB on the orbital period derivative, P, due to wave breaking for the studied systems. Colours and shapes of the points follow the same scheme as in Fig. 3.
In the text

	Fig. 5 Joint plot of the host mass and EEP for the studied systems, with symbols and colours defined as in Fig. 3.
In the text

	Fig. 6 Same as Fig. 3, but including the evolutionary tracks of planetary orbits for all systems in our sample. The coloured solid lines trace the orbital evolution driven by wave breaking: olive lines indicate systems ending in engulfment, cyan lines show those disrupted at the Roche limit, and brown lines denote tracks that terminate early due to the limits of the MESA grid.
In the text

Fig. 7 Joint plot of EEP of destruction and the present-day orbital period. The point types and colour coding follow those in Fig. 3. The right-hand panel shows two separate histograms: one for systems destroyed before the TAMS (green), one for systems destroyed after the TAMS (pink). The overlapping region of the distributions is shown with semi-transparent shading. A clear separation between the two populations emerges at P ≈ 3.5 d.

In the text

	Fig. 8 Predicted year at which the cumulative shift in central transit times reaches 5 minutes, plotted as a function of the absolute value of the orbital period change |P| due to wave breaking. The horizontal black line indicates the year 2035.
In the text

	Fig. B.1 Mosaic of selected parameters for the studied systems, illustrating correlations between stellar, planetary, and tidal properties derived from our modelling. The dotted and dash-dotted vertical lines on EEP and EEPdest plots mark IAMS and TAMS evolutionary phases, respectively.
In the text