# The intermittently-resonant coevolution of migrating planets and their
pulsating stars.

Jared Bryan1 Julien de Wit1 Meng Sun2,3 Zoe Zoë L. de Beurs1 Richard H. D.
Townsend2

###### Abstract

Hot Jupiters are expected to form far from their host star and move toward
close-in, circular orbits via a smooth, monotonic decay due to mild and
constant tidal dissipation[dawson2018origins, eggleton1998equilibrium]. Yet,
three systems have recently been found exhibiting planet-induced stellar
pulsations suggesting unexpectedly strong tidal interactions[de2017planet,
kalman2022gravity, kalman2023discovery]. Here we combine stellar
evolution[paxton2010modules, paxton2013modules, paxton2015modules,
paxton2018modules, paxton2019modules, jermyn2023modules] and tide
models[townsend2013gyre, townsend2018angular, goldstein2020contour,
sun2023tides] to show that dynamical tides raised by eccentric gas giants can
host give rise to chains of resonance locks with multiple modes, enriching the
dynamics seen in single-mode resonance locking of circularized
systems[ma2021orbital]. These series of resonance locks yield orders-of-
magnitude larger changes in eccentricity and harmonic pulsations relative to
those expected from a single episode of resonance locking or nonresonant tidal
interactions. Resonances become more frequent as a star evolves off the main
sequence providing an alternative explanation to the origin of some stellar
pulsators and yield yielding the concept of “dormant migrating giants”. We
also show that evolution Evolution trajectories are characterized by competing
episodes of inward/outward migration and spin-up/-down of the star which are
sensitive to the system parameters, revealing a new challenge in modeling
migration paths and in contextualizing the observed populations of giant
exoplanets and stellar binaries. This sensitivity however offers a new window
to constrain the stellar properties of planetary hosts via tidal
asteroseismology[guo2022new].

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

Department of Astronomy, University of Wisconsin – Madison, Madison, WI, USA;

Center for Interdisciplinary Exploration and Research in Astrophysics,
Department of Physics & Astronomy, Northwestern University, Evanston, IL, USA

##  1 Introduction

Hot Jupiters were the first exoplanets found around main sequence
stars[mayor1995jupiter], but have no analog in our solar system. Understanding
the formation channels of hot Jupiters provides a strong test for planetary
formation theories. Although it has been suggested that Hot Jupiters can form
in situ[batygin2016situ], their formation beyond the ice line is favored,
followed by inward migration. Inward migration proceeds either by interaction
with the protoplanetary disk[goldreich1980disk, lin1986tidal] or by a
combination of eccentricity excitation (e.g. by planet-planet
scattering[rasio1996dynamical] or by cyclic[kozai1962secular,
lidov1962evolution] or chaotic secular interactions[wu2011secular]) followed
by tidal dissipation[barker2009tidal, dawson2018origins] in the
planet[eggleton1998equilibrium] or its host star[jackson2008tidal]. Abundant
evidence for tidal interactions exist for heartbeat stellar binaries in the
form of photometric observations of tidally excited
oscillations[fuller2017accelerated, hambleton2018kic]. High-amplitude stellar
pulsations and rapid orbital evolution have appeared as evidence for
unexpectedly strong tidal interaction in the 2.6 Gyrs-old[bonomo2017gaps]
system hosting the eccentric (e$\sim$0.52) hot Jupiter HAT-P-2b[bakos2007hd,
de2017planet][bakos2007hd, de2017planet, de2023revisiting], with similar high-
amplitude stellar pulsations in the WASP-33[herrero2011wasp,
kalman2022gravity] and HD-31221[kalman2023discovery] systems. Unexpectedly
efficient energy and angular momentum transfer between a hot Jupiter and its
host may provide new insights into their observed population and also yield
frequent planetary engulfment as recently observed in the ZTF SLRN-2020
system[de2023engulf]. Prompted by these systems, we investigate the tidal
evolution of eccentric hot Jupiters.

Tidal migration of hot Jupiters due to stellar tides occurs through damping of
large-scale tidal distortion (the equilibrium tide) by turbulent viscosity in
the convective envelope[zahn1977tidal, hut1981tidal, eggleton1998equilibrium],
as well as by radiative damping of dynamically excited internal gravity waves
(the dynamical tide)[cowling1941non, zahn1977tidal, goodman1998dynamical].
Damping of inertial waves[ogilvie2007tidal, barker2022tidal] and nonlinear
damping due to wave breaking near the stellar center[weinberg2012nonlinear,
weinberg2017tidal] or stellar surface[macleod2022tidal] have also been
considered by previous work.

Most previous work parameterized the stellar response to a tidal perturber in
terms of a constant tidal quality factor[goldreich1966q], or a frequency-
averaged tidal quality factor[ogilvie2004tidal]. Since more sophisticated
models of tidal evolution[witte1999tidal, sun2023tides] demonstrate that
orbital evolution rates are a sensitive function of the forcing frequency and
the stellar structure, these models are insufficient to capture the dynamics
of tidal evolution.

Previous work has also considered tides to be a purely dissipative process.
Some authors have studied inverse tides due to the interaction between tidally
forced oscillations and self-excited stellar oscillations, though only for
circular stellar binaries[willems2003nonadiabatic, fuller2021inverse]. Inverse
tides arising from Doppler shifting of stellar modes due to stellar rotation
have also been studied for solar-type stars[witte2002orbital]. Another thread
of research has dealt with chaotic tidal evolution for high-eccentricity
($e\geq 0.95$) exoplanets where energy can be traded between the planetary
f-mode and the orbital angular momentum[ivanov2004tidal, wu2018diffusive,
veras2019tidal, vick2019chaotic]. The oscillation of specific orbital elements
during single resonances have also been explored, including the orbital
inclination[millholland2019obliquity, alexander2023orbital] through obliquity
tides. In general, these resonances are modeled over short timescales, similar
to applications considering stellar mode amplitude through dynamical resonance
locking[burkart2014dynamical].

Much work on high eccentricity tidal migration considers only the effect of
tides within the planet, neglecting tides in the star[goldreich1966q,
ogilvie2004tidal, vick2019chaotic]. We consider dynamically excited g-modes in
stars. These can be resonantly excited, resulting in orders of magnitude
higher pulsation amplitudes and rates of energy and angular momentum exchange
between the orbit and the star than is possible with equilibrium tides
alone[witte2002orbital, ma2021orbital].

It is also rare to see coupled stellar and orbital evolution, precluding the
existence of resonance locking. Resonance locking occurs when the tidal
forcing frequency and the frequency of a stellar pulsation mode vary in
concert, allowing for resonant interactions to be sustained over longer
timescales than if either of these frequencies were held constant. Resonance
locking has been studied in binary stars[witte1999tidal, fuller2012dynamical,
burkart2012tidal, burkart2013tidal, burkart2014dynamical]. Observational
evidence comes from high-amplitude tidally excited oscillations in heartbeat
stars[fuller2017accelerated, hambleton2018kic]. Recent work examined resonance
locking as a mechanism for rapid orbital migration of a massive planet on a
circular orbit[ma2021orbital]. It has also been proposed to drive orbital
migration of massive exoplanets on eccentric orbits around solar-type
stars[witte2002orbital], though nonlinear tidal dissipation due to wave
breaking in the radiative core limits may prevent the application of resonance
locking in this setting[barker2010internal]to the most massive
exoplanets[barker2010internal, ma2021orbital]. It may also drive the orbital
evolution of Saturn’s and Jupiter’s moons[fuller2016resonance,
lainey2020resonance]. However, resonance locking of massive exoplanets on
eccentric orbits around main sequence stars has received less attention, and
it remains unclear how the enriched tidal forcing spectrum affects their tidal
evolution.

We build on these studies by examining coupled stellar and orbital evolution
for gas giants on eccentricorbits around , spin-aligned orbits around F-type
main sequence stars. We use the MESA stellar evolution code[paxton2010modules,
paxton2013modules, paxton2015modules, paxton2018modules, paxton2019modules,
jermyn2023modules] to simulate the structural evolution of main sequence
stars. We focus on main sequence stars for which radiative damping of the
dynamical tide is the dominant damping mechanism due to the large radiative
envelope (Fig. 1). These stars lack a deep convective envelope, so we neglect
magnetic braking of the stellar rotation[gossage2023magnetic]. We then use the
GYRE stellar pulsation code[townsend2013gyre, townsend2018angular,
goldstein2020contour, sun2023tides] to directly solve for the response of
these stellar models to periodic tidal forcing by an eccentric point-mass
companion. In particular, we use GYRE-tides[sun2023tides] to solve for the
instantaneous rate-change of the orbital eccentricity, orbital semi-major
axis, and stellar spin rate. Beginning with an initial stellar and orbital
configuration similar to that of HAT-P-2b[bonomo2017gaps, stassun2019revised],
we integrate the orbital configuration and stellar spin state forward and
backward in time. In contrast to models of hot Jupiters on circular orbits, we
demonstrate how orbital eccentricity greatly enriches the tidal forcing
spectrum (Fig. 1). In particular, we examine how competing resonances between
different modes can lead to a wide array of tidal evolution behaviors,
including brief resonance sweeps, sustained resonance locking, and wandering
tidal migration.

##  2 Results

###  2.1 Episodic and wandering tidal migration

Periods of slow orbital evolution are punctuated by rapid changes as the tidal
forcing sweeps through resonance with a stellar eigenmode. The equilibrium
tidal response occurs at rates of $\dot{e}\sim 10^{-11}$/yr, $\dot{a}\sim
10^{-12}$ au/yr and acts to circularize and shrink or expand the orbit while
spinning up the star (Fig. 2a,b,e). This diversity in the senses of orbital
evolution is possible due to the eccentric orbit interacting with multiple
modes with different properties (Fig. 2e). However, resonances between the
tidal forcing frequencies and natural stellar pulsation modes result in orders
of magnitude more rapid orbital evolution, up to $\dot{e}\sim 10^{-7}$/yr in
the particular case shown in Fig. 2 and $\dot{e}\sim 10^{-3}$/yr for an
ensemble of stars shown in Fig. 3a. Most of these resonances are exceedingly
brief, lasting $\ll 1$ Myr at their peak. However, when the rate of change of
the frequency of the stellar pulsation mode due to stellar evolution and spin
evolution matches the rate of change of the orbital forcing frequency,
sustained excitation of individual stellar eigenmodes results in orbital
evolution rates of $\dot{e}\sim 10^{-8}-10^{-9}$/yrs over multi-Myr
timescales, thereby driving the bulk of orbital evolution (Fig. 2, Fig. 3b).
This is a process known as resonance locking.

By contrast to planets on circular orbits these resonance locks do not persist
over Gyr timescales[ma2021orbital]. This is because other modes compete with
the resonantly locked mode, eventually pushing it out of resonance and
resulting in relatively brief resonance locks ($\sim 1-5$ Myr). However, the
eccentric orbit enables a rich spectrum of possible resonance locks,
decreasing the time between resonances. In addition, the forest of possible
resonances becomes more dense as the star evolves through the main sequence
and we observe chains of resonance locks over timescales up to $\sim 25$ Myr
(Fig. 2b inset). Resonance locks can be seen as quasi plateaus in the orbital
evolution rate (Fig. 2a,b insets). This episodic tidal migration implies that
observations of orbital evolution rates are generally not reflective of the
long-term rates of orbital migration. Instead, observations of high orbital
evolution rates may be brief episodes in a migration history otherwise
dominated by low tidal interaction strength. Likewise, orbits not currently
observed to be migrating cannot be assumed to be stationary over long
timescales.

###  2.2 Sensitivity of tidal evolution to initial stellar and orbital
configurations.

The tidal evolution trajectory depends on the stellar and orbital
configuration. The array of stellar modes available for tidal interaction
depends sensitively on the stellar structure. In Fig. 3a, we show the rate-
change in eccentricity for a range of stellar models between the zero-age main
sequence and the red-edge of the main sequence for a fixed orbital
configuration. The stellar models range from $M=1.2-1.5M_{\odot}$ with a
constant metallicity of $Z=0.02$. No coherent patterns exist across the HR
diagram due to the superposition of many different stellar modes which can
contribute to the orbital evolution rate (patterns exist when considering one
mode in isolation–see Methods). Without tuning the stellar models to fall
exactly on resonances, we find orbital evolution rates up to $\dot{e}\sim
10^{-3}/yr$, far beyond the equilibrium tidal response of $\dot{e}\sim
10^{-11}$/yr. In simulations with coevolving orbits and stellar structure, we
observe $\dot{e}\sim 10^{-4}$/yr (see Methods).

To further investigate the sensitivity of the tidal evolution of massive
exoplanets, we perform a sensitivity analysis of the tidal evolution
trajectories for uncertainties in the initial orbital and stellar
configuration. To this end, we consider the parameters of the HAT-P-2
system[bonomo2017gaps, stassun2019revised] and sample the $+1\sigma$,
$-2\sigma$, and $+3\sigma$ in $a$, $e$, and $\Omega_{\mathrm{rot}}$, and $\pm
1\sigma$ in the $M_{\odot}$, $Z_{\odot}$$M$, $Z$, and $t_{0}$. We then evolve
the stellar structure and orbital configuration forward and backward in time
and evaluate the similarity of the resulting orbital evolution trajectories.
We find that uncertainties in orbital parameters ($e$ and $a$) are small
enough that orbital trajectories retain similar structure over $\sim 50$ Myr
(Fig. 3b). Similarities in the shape of the orbital trajectories imply
similarities in the set of resonantly excited stellar modes. The exact timing
of these resonances differs on the order of $\sim 10$ Myr, far below the
uncertainty in the stellar age and certainly below the main sequence lifetime
of the star. These small variations in the stellar structure and orbital
configuration after a resonant tidal interaction cascade into larger changes
for the next resonance, and the tidal evolution trajectories soon diverge.
Uncertainties in $t_{0}$, $\Omega_{\mathrm{rot}}$, $M_{\odot}$, and
$Z_{\odot}$ $M$, and $Z$ immediately lead to qualitatively different tidal
evolution behavior, both in the rates of orbital evolution and the specific
timing of resonant interactions. This sensitivity serves as a reminder of
exoplanetary science’s favorite tenet: “Know thy star, know thy planet” and
highlights the challenge to our ability to robustly reconstruct the initial
orbital configuration of hot Jupiters by time-reversing their tidal evolution.
Reducing uncertainties in stellar parameters is thus critical for robust
modeling of exoplanet tidal migration.

##  3 Discussion

###  3.1 Challenges in reversing and propagating tidal evolution trajectories.

Efforts to reconstruct the initial (prior to tidal migration) orbital
configuration of hot Jupiters rely on a model of tidal migration. Simplified
theories based on a constant or frequency-averaged tidal quality factor lead
to a smooth, monotonic evolution toward circular, synchronous, and short
period orbits[jackson2008tidal]. Models based on sophisticated
characterizations of the tidal response have tended to focus on the case of
circular orbits[ma2021orbital], which again gives the impression of smooth,
monotonic evolution of the orbit. We find that orbital evolution is controlled
by excitation of stellar g-modes, which are transiently excited to high
amplitude and lock into resonance due to the coevolution of stellar modes,
stellar spin rate, and orbital forcing frequencies leading to a strong
sensitivity of the overall evolution trajectory to the system parameters.
Small perturbations to the orbital and stellar configuration cause the orbital
trajectories to deviate above $1\sigma$ in the orbital parameters within 1-15
Myr (see methods).

This sensitivity to the system parameters highlights a challenge to our
ability to time-reverse and propagate the tidal evolution calculations in
order to deterministically reconstruct the initial orbital configuration of
hot Jupiters, a step that is key to contextualize the observed populations of
hot Jupiters vs warm Jupiters, circularized vs eccentric ones. Further work to
investigate the statistics of tidal evolution trajectories would be needed to
adequately account for and propagate the degree of confidence in the primitive
orbital configurations of particular exoplanets onto their evolution
trajectory.

This sensitivity to initial stellar and orbital configuration is also relevant
to binary population synthesis studies. In binary population synthesis code,
one of the goals is to determine the final mass and orbital distribution of
compact objects, in order to aid in gravitational wave data analysis. To
achieve this, an initial mass, orbital period, and mass ratio distribution are
provided to launch large grids of binary evolution models. Due to
computational constraints, numerical calculations of the tidal dissipation
rate at each stellar evolutionary step, which require a stellar profile, are
often replaced by fitting formulae and semi-analytical equations. In this
work, we demonstrate that tidal evolution in binaries is sensitive to
parameters related to stellar and binary evolution, which implies that
constructing more accurate parameterized equations for fast calculations could
be problematic, particularly for the higher mass-ratio case of stellar
binaries.

###  3.2 The source of sign changes in the orbital evolution rates

The sensitivity of the trajectories also finds its origin in the sign changes
in the orbital evolution rates. This sign of the orbital evolution rates
depends on three quantities: (i) the sign of the Doppler-shifted forcing
frequency, (ii) whether the position of the companion on-net on-average lags
or leads the location of the tidal disturbance for a particular mode, and
(iii) the stability of the excited mode. First, sign changes in the Doppler-
shifted forcing frequency
$\sigma_{m,k}=k\Omega_{\mathrm{orb}}-m\Omega_{\mathrm{rot}}$ determine whether
the excited mode is prograde or retrograde in the corotating frame of the
star. Fig. 2 shows that a star rotating at $1.5\times$ the pseudosynchronous
rotation rate[hut1981tidal, sun2023tides] but otherwise with the orbital
configuration of the HAT-P-2 system[bonomo2017gaps] has $\sigma_{2,k}<0$ for
$k<9$, and thus $\dot{e}>0$ (Fig. 2d). Higher rotation rates lead to
additional retrograde modes, while lower rotation rates lead to only prograde
modes (see Methods). Also visible in Fig. 2d is that $k=1,2$ retain the
circularization sign of $\dot{e}$ ($\dot{e}<0$) despite being retrograde in
the corotating frame. This is because the tidal distortion for $m=2$, $k=1,2$
modes leads the position of the companion in an orbit-averaged sense, while
otherwise for other $k$ it lags. When combined with the negative Doppler-
shifted forcing frequency, this results in modes that act to circularize the
orbit. The final source of sign changes in orbital evolution rates is due to
oscillation mode instability such that the heat-engine of the star supplies
energy to the planet’s orbit, a phenomenon called inverse
tides[willems2003nonadiabatic, fuller2021inverse]. The main sequence stars we
study do not host unstable modes, so this mechanism is not applicable to this
work. Additional discussion of these mechanisms can be found in Sun et al.
(2023).

###  3.3 Insights into stellar properties and orbital architecture from
resonance locking.

The sensitivity of the orbital evolution rates to perturbations in the stellar
and orbital parameters implies in return access to tight observational
constraint on these parameters, also known as tidal
asteroseismology[guo2022new]. The orbital evolution rates across the HR
diagram display no discernible pattern (Fig. 3a), but this complexity is
simply due to the superposition of $\dot{e}$ from every mode. In contrast,
each individual mode is coherent across the HR diagram (see Methods).
Unfortunately, only the net orbital evolution rate is observationally
accessible, but this suggests that measurements of individual mode frequencies
and amplitudes can complement the integrated orbital evolution rates to
provide a strong constraint on stellar structure and orbital configuration.

To highlight this further, we calculate the stellar flux variations due to
tidally excited low-frequency g-modes[townsend2003semi]. We account for band-
limitation of the observer and assume an edge-on viewing geometry. Fig. 4a
shows the stellar flux evolution over 400 Myr, where colored modes are those
that exceed amplitudes of 10 p.p.m., which we take to be the limit of
photometric observability. The amplitudes of these tidally excited g-modes
generally increase with decreasing $k$ due to the concentration of tidal
forcing strength at low frequencies. When the tidal forcing sweeps through
resonance, modes can be excited to high amplitudes (10%), with sustained
excitation of modes up to $k=43$ (present model only account for k$\leq$50) to
amplitudes of 100 p.p.m. during resonance locking (Fig. 4a), consistent with
Ref.[de2017planet].

Since high-frequency modes are often excited to amplitudes an order of
magnitude above the observational limit, wavelength decomposition of the
stellar flux can provide additional information about the nature of the
excited mode and thus can be used as an additional measurement for probing the
stellar interior. Fig. 4b,d,e show the pulsation spectra in six independent
wavelength bands between $0.6-5\mu$m at three different stages in the
pulsation evolution, with Fig. 4c,e,f showing the corresponding synthetic
lightcurves. Fig. 4b,c shows excitation of a single high-frequency mode with
weak wavelength dependence. Fig. 4d,e shows excitation of a single low-
frequency mode with strong wavelength dependence. This is useful since other
sources of power at low orbital harmonics exist, including Doppler beaming and
reflection of light off of the planet[penoyre2019higher]. These models of
wavelength-dependent photometry can be used not only to probe the stellar
interior, but also to refine models of these alternative mechanisms. Fig. 4f,g
shows excitation of two high frequency modes at once, each with a strong
wavelength dependence. This situation is likely less observationally
prevalent, but is nonetheless important to understand since a resonance lock
can be reinforced or broken by excitation of a second mode. Understanding the
dynamics of competing resonance locks thus depends on detailed
characterizations of the photometric observability of multi-mode excitation.

In addition to providing insights into the stellar structure[guo2022new], the
detection of high-frequency pulsations can also be observables for eccentric
companions regardless of their orbital configuration. Such a
strategy[de2017planet] would be uniquely suited to eccentric and inclined
systems[millholland2019obliquity, alexander2023orbital], whose planetary
candidates may later be confirmed via traditional techniques such as radial
velocity, astrometry, and/or direct imaging.

###  3.4 Implications for tidal modification of stellar evolution

Orbital angular momentum is usually assumed to be deposited in the star in a
way that maintains rigid-body rotation. In reality, however, the deposition
will be concentrated where the tidal energy is dissipated–typically, by the
surface layers. Over timescales short compared to stellar evolution, this can
establish a shear layer between the surface and the deeper envelope, as shown
for self-excited oscillations[townsend2018angular]. This can have two
important consequences. First, the shear layer can trigger mixing via the
Kelvin-Helmholtz instability. Second, a rich new phenomenology of
“differential resonances” will be induced, through which only part of a star
participates in a tidal resonance. It is also possible that angular momentum
deposition in the chemically inhomogeneous zone just outside the boundary of
the convective core[townsend2018angular] can lead to mixing of the zone, which
can grow the core size. The mass of the compact object is directly linked to
the core size during the main sequence and later phases of evolution, so
enhanced mixing has important implications for the subsequent evolution and
fate of the star. We note that the only 3 known pulsating exoplanetary
hosts[de2017planet, kalman2022gravity, kalman2023discovery] are found at the
intersection of the instability strip and the end of the main sequence, which
we may warrant future investigations with a focus on linear and non-linear
instabilities (i.e., possible triggering of self-excitation through a
resonance).

###  3.5 Implications for the migration timescales and populations of hot
Jupiters.

Resonance locking may help explain the lack of planets currently observed to
be tidally migrating due to both the enhanced circularization and migration
rates leaving few planets observed in the process of migrating, and the
increased migration efficiency towards the end of the main sequence resulting
in young gas giants with mild migration rates and latent tidal migration
potential (“dormant migrating gas giants”). The latter can be seen in the the
skew of high orbital evolution rates toward the red-edge of the main sequence
(Fig. 3a) and more rapid orbital evolution for older stars (Fig. 3b). Tidal
migration rates increase with stellar age due to the development of a g-mode
cavity at the core-envelope boundary that can strongly couple to the tidal
forcing (Fig. 1c). This leads to an increase in the number of resonances
encountered and thus enhanced orbital migration rates. The increase in stellar
radius as the star ages also acts to increase the strength of the tidal
forcing, accounting for an increase of $\sim 4\times$ over the main sequence
lifetime of a $M=1.35M_{\odot}$ star.

Efficient tidal migration near the end of the main sequence points toward a
tighter limit on the survival of close-in massive planets beyond the main
sequence, with engulfment such as the one recently observed as the ultimate
fate[de2023engulf]. This also suggests the possibility of two populations of
migrating gas giants. A first population yields the hot Jupiters seen around
early main-sequence stars which are unlikely to have formed due to resonance
locking as introduced here. A second population evolves from “dormant
migrating giants” on slowly-evolving eccentric orbits yielding hot Jupiters
once their hosts start evolving off the main-sequence. Such hosts will often
display high-amplitude tidally induced oscillations and support enhanced tidal
migration rates.

##  4 Methods

###  4.1 Stellar Evolution

We use the MESA stellar evolution code[paxton2010modules, paxton2013modules,
paxton2015modules, paxton2018modules, paxton2019modules, jermyn2023modules] to
simulate the main sequence evolution of stars with masses in the range of
$M=1.2M_{\odot}$ to $M=1.5M_{\odot}$, and metallicities in the range
$Z=0.01-0.04$. Stars of this type have a convective core, preventing excited
g-modes from geometrically focusing in the core of the star, which can lead to
nonlinearities from wave breaking[weinberg2012nonlinear] for some massive
exoplanets [ma2021orbital]. We calculate the linear tidal response with GYRE-
tides which excludes inherently nonlinear dissipation mechanisms. Models start
at the zero-age main sequence and are evolved to the terminal-age main
sequence, but the HR diagrams are truncated at the red-edge of the main
sequence for clarityto avoid overlap with the Henyey hook. We evaluate
convective stability with the Ledoux criterion, though semiconvection and
thermohaline mixing are found to have a negligible effect on the stellar
evolution. We include convective overshoot in the exponential scheme. We do
not consider mass loss or stellar wind since these are negligible in the type
of main sequence stars considered in this study. We neglect spin in the
stellar evolution calculations because their evolution does not significantly
change their spin rate, and magnetic braking is negligible for stars with a
radiative envelope. We do, however, consider the evolution of stellar spin in
the orbital evolution calculations. Example MESA inlists are provided in the
supplementary materials.

###  4.2 Tidal Asteroseismology for Orbital Evolution

We consider a binary system consisting of a star with mass $M_{1}$ and radius
$R$, along with a point-mass companion of mass $M_{2}$. The host star spins
with a frequency $\Omega_{\mathrm{rot}}$. We neglect Coriolis and centrifugal
forces because low-mass main-sequence stars rotate well below their critical
rotation rates, which means that that inertial waves cannot be resonantly
excited. We adopt a non-rotating reference frame and assume spin-orbit
alignment. Although spin-orbit misalignment is possible to include in the
linear framework of GYRE-tides, its implementation in these calculations is
beyond the scope of this work. The companion orbits with a frequency
$\Omega_{\mathrm{orb}}$ given by Kepler’s third law

$GM_{1}(1+q)=a^{3}\Omega_{\mathrm{orb}}^{2},$ (1)

where $G$ is the gravitational constant, $q=M_{2}/M_{1}$ is the mass-ratio of
the system, and $a$ is the orbital semi-major axis. We study the tides raised
on the host star due to the tidal portion of the secondary’s gravitational
potential. Using a multipolar expansion in space and Fourier-series expansion
in time, we express the tide-generating potential at position $\mathbf{r}$ and
time $t$ as a superposition of partial tidal potentials

$\displaystyle\Phi_{T}(\mathbf{r},t)=$
$\displaystyle-\varepsilon_{T}\frac{GM}{R}\sum_{l=2}^{\infty}\sum_{m=-l}^{l}\sum_{k=-\infty}^{\infty}\bar{c}_{l,m,k}\left(\frac{r}{R}\right)^{l}\times$
(2) $\displaystyle
Y_{l}^{m}(\theta,\phi)exp\left[-ik\Omega_{\mathrm{orb}}(t-t_{0})\right].$

where $Y_{l}^{m}$ is a spherical harmonic of degree $l$ and order $m$, $t_{0}$
is the time of periastron passage, $\bar{c}_{l,m,k}$ is a tidal expansion
coefficient given by Sun et al. (2023), and $\varepsilon_{T}$ is a
dimensionless parameter describing the strength of the tidal forcing given by

$\varepsilon_{T}=\left(\frac{R}{a}\right)^{3}\left(\frac{M_{2}}{M_{1}}\right).$
(3)

For hot Jupiters orbiting intermediate-mass main-sequence stars,
$\varepsilon_{T}\sim 10^{-5}$.

The tide-generating potential perturbs the spherical symmetry of the star’s
gravitational potential, which in turn perturbs the orbit of the companion
from a pure Keplerian orbit. The magnitude and direction of these
perturbations to the orbital parameters depend on how dissipation of the
tidally-induced pulsations acts to transfer energy and angular momentum
between the star and the orbit. The star’s tidal response consists of the
equilibrium and dynamical tides. We account for damping of the equilibrium
tide by turbulent viscosity in the convective envelope, as well as by
radiative damping of dynamically excited internal gravity waves (the dynamical
tide). We use the GYRE-tides stellar pulsation code[sun2023tides] to directly
solve the linear, non-radial, non-adiabatic stellar oscillation equations with
an inhomogeneous tidal forcing term. Example GYRE-tides inlists are provided
in the supplementary materials.

The changes in orbital elements through these mechanisms are given by

$\displaystyle\left(\frac{de}{dt}\right)_{sec}=4\Omega_{\mathrm{orb}}q\sum_{l,m,k\geq
0}$
$\displaystyle\left(\frac{R}{a}\right)^{l+3}\left(\frac{r_{s}}{R}\right)^{l+1}$
(4)
$\displaystyle\times\kappa_{l,m,k}\text{Im}(\bar{F}_{l,m,k})\bar{G}^{(3)}_{l,m,k}(e)$

where the $sec$ subscript denotes a secular (orbit-averaged) quantity.
$\bar{F}_{l,m,k}$ are dimensionless quantities that measure the response of
the star to the various forcing frequencies $\sigma_{m,k}$. $\kappa_{lmk}$ are
the same as those given by Willems et al. (2010). $\bar{G}^{(3)}_{lmk}$ is
given by Sun et al. (2023), along with $\bar{G}^{(3)}_{lmk}$ corresponding to
and differs only by sign convention from the coefficients derived by Smeyers
et al. (1998) and Willems et al. (2003). Similar $\bar{G}$ coefficients exist
for changes in the semi-major axis, stellar rotational angular momentum, and
argument of periastron. We account for the effect of rotation only in its
Doppler-shift of the co-rotating frame forcing frequency.

While eq. 2 and eq. 4 provide the exact expressions for the orbital evolution
rate, computational considerations require us to truncate the summations in
$l$ and $k$. Since the orbital evolution rate scales with $(R/a)^{l+3}$, where
$(R/a)\approx 0.1$, the strength of tidal forcing falls off rapidly above
$l=2$. We restrict $m$ by noting that $l=2$, $m=\pm 1$ modes are not excited
by the tidal potential, and we find that orbital evolution rates due to $l=2$,
$m=-2$ modes are consistently several orders of magnitude smaller than for
$m=0,2$ modes. We calculate only the $l=2$, $m=0,2$ modes[witte2002orbital].
The orbital evolution rates also decrease for sufficiently high $k$, and we
truncate the summation at $k=50$.

###  4.3 Benchmarking accelerated tidal migration

In order to evaluate whether these tidal processes can contribute to rapid
circularization and inward migration, we compare the orbital evolution
produced by GYRE-tides and by a constant-Q model of stellar and planetary
tides from Jackson et al. (2008). In this framework[goldreich1966q,
jackson2008tidal], the rate-change in eccentricity is given by

$\frac{de}{dt}=-e\left[\frac{63}{4}(GM_{1}^{3})^{1/2}\frac{R_{2}^{5}}{Q_{p}M_{2}}+\frac{171}{16}(G/M_{1})^{1/2}\frac{R_{1}^{5}M_{2}}{Q_{\ast}}\right]a^{-13/2},$
(5)

where $R_{1}$ is the stellar radius, $R_{2}$ is the planetary radius, and
$Q_{\ast}$ and $Q_{p}$ are quality factors for the star and the planet,
respectively. A similar expression exists for $da/dt$. Jackson et al. (2008)
estimate $Q_{\ast}=10^{5.5}$ and $Q_{p}=10^{6.5}$, though they consider a wide
range of possible values from $10^{4}-10^{8}$ in each parameter. Extended data
Fig. 2 compares the variability in tidal evolution trajectories for the GYRE-
tides model to the constant-Q model. Parameters in the constant-Q model are
perturbed from $-3\sigma$ to $+3\sigma$ in each parameter, while the
perturbations to the GYRE-tides model are given in the legend to the right of
extended data Fig. 2. While the perturbations to parameters in the GYRE-tides
model are generally smaller (e.g. perturbations of $1\sigma$ in $M_{\odot}$ in
the GYRE-tides model vs. perturbations of $3\sigma$ in the constant-Q model),
the resulting variability in tidal evolution trajectories is much higher for
the GYRE-tides model. This is especially true for well-constrained parameters
such as $a$, where no deviation of the tidal evolution trajectories is visible
in the constant-Q case, while significant divergence in tidal trajectories
emerges after only $\sim 50$ Myr for the GYRE-tides model. Variations in $e$
in the constant-Q model results result in parallel tidal evolution
trajectories that maintain the initial uncertainty in $e$. The GYRE-tides
model, on the other hand, increases the difference between tidal evolution
trajectories as time elapses. Of particular interest are the variations in the
stellar mass. The constant-Q model vastly underestimates the variation in
tidal evolution trajectories due to variations in stellar mass. Perturbations
of $\pm 3\sigma$ in stellar mass result in smaller variations than
perturbations of $\pm 1\sigma$ in the GYRE-tides model.

We also see that the tidal evolution trajectories are most sensitive to the
choice of stellar quality factor, with the planetary quality factor playing a
less significant role. The wide range of orbital evolution rates resulting
from varying $Q_{\ast}$ over four orders of magnitude provides a convenient
way to benchmark whether a tidal evolution trajectory displays accelerated
tidal migration relative to the baseline $Q_{\ast}=10^{5.5}$.
$Q_{\ast}=10^{5.5}$ results in $\dot{e}$ on the order of $10^{-10}$/yr over
the studied period.

The background $\dot{e}$ is on the order of $10^{-11}$/yr for the GYRE-tides
model, so without resonance the circularization timescale should be an order
of magnitude larger (extended data Fig. 1, extended data Fig. 3, extended data
Fig. 4). However, brief resonances sweep through $\dot{e}\sim 10^{-5}$/yr and
resonance locks sustain rates of $10^{-8}-10^{-9}$$10^{-8}$/yr, which
compensates for the lower background orbital evolution rates. In cases where
few resonances are encountered, such as for early young main sequence stars
($t_{0}-1\sigma$), the equivalent $Q_{\ast}$ is as high as $\sim 10^{7}$. For
cases where many resonances are encountered, such as for late old main
sequence stars ($t_{0}+1\sigma$), the equivalent $Q_{\ast}$ is closer to $\sim
10^{5}$. Similarly, perturbing the stellar mass by $+1\sigma$ results in an
effective $Q_{\ast}$ of $10^{4.5}$. Increasing the mass of the star while
holding the age constant has the same effect as advancing the age while
holding the mass constant since higher mass stars have short main sequence
lifetimes. This emphasizes that tides can be an effective driver of
accelerated orbital evolution in late old main sequence stars, but are
unlikely to accelerate orbital migration in early young main sequence stars.

##  5 Tidal asteroseismology for photometric observations

###  4.1 Limits of orbital trajectory predictability

The deviation of orbital evolution trajectories due to small perturbations in
stellar properties is of fundamental importance to understanding whether tidal
evolution calculations can be used to deterministically predict the future
state of exoplanet systems or reconstruct their formation state. We
characterize the predictability of tidal evolution trajectories through a
series of fixed-orbit calculations similar to those in fig. 3a, but we seek to
understand whether features of the net orbital evolution rate are preserved
under small perturbations to the stellar parameters. We simulate the stellar
evolution of 15 stars, with perturbations in mass and metallicity ranging from
0.1 to 5% with the baseline star from the main text ($M=1.36M_{\odot}$,
$Z=0.02262$) used as a starting point. For each of these stars, we calculate
the orbital evolution rates for a fixed-orbit tidal perturber for 8 different
spin-orbit states; for similarity to fig. 3, we focus on an orbit with
$e=0.575$, $a=0.153$ au, and a pseudosynchronous stellar spin rate. Individual
resonances can be easily matched between the reference star and those with
small mass perturbations (say, 0.1%) (extended data fig. 5). Larger mass
perturbations distort the shape and location of these resonances, and only the
general increase in density of resonances toward the end of the main sequence
remains visible at perturbations of $5\%$. Perturbations in metallicity result
in smaller deviations in the shape of orbital evolution rates, though the
location of individual resonances in time is nonlinear with increases in $Z$,
another representation of the zig-zag pattern in extended data fig. 12b.

We now calculate how this variation in orbital evolution rates affects
predictability of the orbital migration trajectories; we quantify the
predictability of an orbital trajectory via the timescale over which its
deviation from a reference trajectory grows to observable levels (say,
$1\sigma$ in e or a). We estimate this timescale by calculating the
instantaneous deviation in the orbital evolution rates. For eccentricity, this
is

$\DIFadd{\Delta\dot{e}=\dot{e}(M_{i},Z_{i},t)-\dot{e}(M_{0},Z_{0},t),}$ (6)

where $M_{0}$, $Z_{0}$ are the reference mass and metallicity and $M_{i}$,
$Z_{i}$ are the perturbed mass and metallicity, $t$ is a common stellar age.
We define a timescale over which the orbital trajectories remain
observationally indistinguishable,

$\DIFadd{T_{e}=\frac{\Delta\dot{e}}{\dot{e}}.}$ (7)

To make this calculation in practice, we begin with time series of $\dot{e}$
for a fixed spin-orbit configuration for a range of stellar mass and
metallicity perturbations (extended data fig. 5, extended data fig. 7). We
aggregate these $\dot{e}$ values into a distribution (extended data fig. 6a),
revealing that the frequency-rate distribution is distributed according to a
power law with slope $-1/2$. This long-tailed distribution of orbital
evolution rates indicates the importance of rare episodes of high orbital
evolution rate in controlling the orbital migration even before considering
resonance locking. The other implication of this distribution is that there is
no well-defined average orbital evolution rate.

We also aggregate the $\Delta\dot{e}$ values for each perturbed stellar model
into a distribution (extended data fig. 6c). The deviation in orbital
evolution rates, $\Delta\dot{e}$, grows with mass perturbations, but is more
stable for metallicity perturbtaions (extended data fig. 8c). To estimate
$T_{e}$, we sample the distribution of $\Delta\dot{e}$ for orbital evolution
rates, which we integrate to produce an ensemble of orbital trajectories. For
each member of the ensemble, we record the time needed to accumulate a change
of at least $\sigma_{e}=0.011$, the uncertainty in $e$ for HAT-P-2b as
reported in [de2023revisiting]. These trajectories do not possess the same
intermittent migration structure due to the removal of time-dependence when
constructing the distribution, but the ensemble allows us to estimate a
distribution of $T_{e}$. The distribution of $T_{e}$ shows that larger mass
perturbations require less time to accumulate observable deviations in orbital
evolution rates, i.e. larger mass differences accumulate differences faster
(extended data fig. 6e). For the spin-orbit state and stars considered in this
example, eccentricity differences become observable over $\sim 5-15$ Myr
timescales.

We also calculate the distribution of $\dot{e}$ and $\dot{a}$ for the live-
orbit simulations with a coevolving star and orbit (extended data fig. 9).
Coevolution of the stellar and orbital configuration modifies the distribution
such that it has a well-defined peak. The live-orbit distribution is depleted
in the lowest $\dot{e}$ and doesn’t reach the highest $\dot{e}$ achieved in
the fixed-orbit simulation. We also calculated $T_{e}$ and $T_{a}$ for the
live-orbit case (extended data fig. 10). We did this by calculating the time
needed to accumulate a $1\sigma$ deviation in e or a from every point in the
orbital trajectories. That is to say, we begin with the initial orbital
configuration at $t_{0}$ and find the first $t$ such that $t>t_{0}$ and
$e(t)-e(t_{0})>e$. We see that the coevolution of the star and orbit enables a
much broader range of $T_{e}$, with the most likely $T_{e}$ around 1 Myr,
though several rare cases reach as low as 1 kyr or as high as 1 Gyr. This
indicates that models of tidal evolution and backpropagation are robust in
comparison to observational constraints over timescales of typically 1 Myr.

###  4.2 Tidal asteroseismology for photometric observations

We calculate the stellar flux variations due to tidally excited low-frequency
g-modes[townsend2003semi]. Stellar flux variations are due to a combination of
radial displacements at the stellar surface, and Lagrangian perturbations to
the radiative luminosity[burkart2012tidal, fuller2012dynamical]. Rather than
assuming black body radiation, we use the CAP18 photospheric
grids[prieto2018collection] to estimate the stellar flux variations between
wavelengths of $0.6-5.0\mu$m at each spherical and orbital harmonic. Since
these photospheric models are not calculated for stars of precisely the mass,
metallicity, age, etc. that we study, we use the MSG code[townsend2023msg] to
interpolate the stellar spectra. Finally, we calculate the observable flux
variations by averaging the luminosity variations over the visible disk. We
assume an edge-on viewing geometry, but do not include transit effects.

Figure 4a shows the evolution of the observed flux variation over the full
$0.6-5.0\mu$m wavelength range. The tidal forcing is strongest for low
frequencies, leading to persistent power on the $k<13$ orbital harmonics at
the $10^{-5}$ level. During resonance between the tidal forcing and a stellar
eigenmode, additional higher frequency modes can be excited to observable
levels (Fig. 4b,c), low-frequency modes can be excited to higher amplitudes
(Fig. 4d,e), and multiple modes can be excited at once (Fig. 4f, g). Also
visible in Fig. 4a is complex resonance locking dyanmics. For example,
excitation of the $k=32$ mode is sustained between $2586-2588$ Myr at an
amplitude of $\sim 10^{-4}$ (Fig. 4a,b,c). $k=4,14,16,30$ modes are resonantly
excited during this time, but the resulting orbital evolution is not enough to
break the resonance lock with the $k=32$ mode. Excitation of the $k=22$ mode
at $2588$ Myr eventually breaks the resonance lock, and the mode amplitude
oscillates over the next $8$ Myr until another resonance lock with the $k=30$
mode is established. Although individual modes are typically locked into
resonance for $<5$ Myr, multiple resonance locks can be chained together
resulting in sustained photometric observability over timescales of $\sim 50$
Myr.

In addition to structure in the dynamics of resonance locking through time,
the tidal response of the star has interesting structure across the HR
diagram. For example, although the net orbital evolution rate has no clear
pattern across the HR diagram (Fig. 3a), this is simply due to the
superposition of many different modes, each with a coherent and unique
fingerprint across the HR diagram (extended data Fig. 311). The orbital
evolution due to a particular mode is not observable since the orbit changes
due to a summation over all modes. However, photometric observations enable
mode identification, and individual modes also have a coherent photometric
signature across the HR diagram (extended data Fig. 412). This tidal
asteroseismology can thus be used to place strong constraints on stellar
structure. For example, if both a $k=10$ and a $k=30$ mode were identified
(extended data Fig. 412), the space of stellar structures can be strongly
constrained to a few specific points on the HR diagram. That is to say,
multiple mode excitation constrains the possible stellar model to the
intersection of the sets of resonances in extended data Fig. 4a and 4b12a and
12b.

###  4.3 Numerical integration for orbital trajectories.

In order to convert the estimates of orbital evolution rates (e.g. $\dot{e}$)
into trajectories in the orbital parameters, we solve a system of differential
equations

$\frac{dy^{m}_{n}}{dt}=f_{n}(t,y)$ (8)

iteratively for the spin-orbit state

$y^{m}_{n}(t+dt_{n})$ (9)

where $dt_{n}$ is the $n^{th}$ timestep and $y^{m}_{n}$ stands in for the
$m^{th}$ spin-orbit parameter (eccentricity, semi-major axis, or stellar
rotation rate). We update the orbital elements and stellar spin frequency at
each time step with the Runge-Kutta-Fehlberg method[press1992adaptive], a
fourth-order accurate method embedded in a fifth-order method used for error
estimation. Error estimation allows for adaptive timestepping depending on the
current dynamics of the system, i.e. the timestep is reduced as the tidal
forcing sweeps through resonances with the stellar eigenmodes. With this
method, the solution is a weighted average of six increments to the spin-orbit
parameters

$y^{m}_{n+1}=y^{m}_{n}+\sum_{i=1}^{6}c_{i}k_{i}+\mathcal{O}({dt^{6}}).$ (10)

where the intermediate increments are given by

$k_{i}=dt_{n}\frac{dy^{m}}{dt}\left(t_{n}+a_{i}dt_{n};y^{m}_{n}+\sum_{j<i}b_{ij}k_{i}\right).\\\
$ (11)

The embedded fourth-order formula gives an alternate estimate of the updated
system state

$y^{m\ast}_{n+1}=y^{m}_{n}+\sum_{i=1}^{6}c^{\ast}_{i}k_{i}+\mathcal{O}({dt^{5}}).$
(12)

The parameters $a_{i}$, $b_{ij}$, $c_{i}$, and $c_{i}^{\ast}$ can be found in
a Butcher tableau for the Runge-Kutta-Fehlberg method[press1992adaptive]. The
error estimate is given by

$\Delta^{m}_{n+1}=y^{m}_{n+1}-y^{m\ast}_{n+1}=\sum_{i=1}^{6}(c_{i}-c^{\ast}_{i})k_{i}.$
(13)

For each orbital element, the error between the 4th and 5th order solutions
suggests a new timestep, given by

$dt^{m}_{n+1}=\begin{cases}S\cdot
dt_{n}\left(\frac{\Delta_{max}}{\Delta^{m}_{n+1}}\right)^{1/5}&(\Delta_{max}\geq\Delta^{m}_{n+1})\\\
S\cdot
dt_{n}\left(\frac{\Delta_{max}}{\Delta^{m}_{n+1}}\right)^{1/4}&(\Delta_{max}\leq\Delta^{m}_{n+1})\end{cases}$
(14)

where $S$ is a safety factor we set to $0.9$. The new timestep is chosen to be
the most conservative of the proposed timesteps

$dt_{n+1}=min\\{dt_{n+1}^{m},dt_{max}\\}_{m}$ (15)

where $dt_{max}$ is an additional safeguard to ensure we do not miss brief
resonances due to extending time steps between resonances. Timesteps near
$dt_{max}$ generally indicate slow orbital evolution, and this parameter is
empirically set to $dt=100$ kyr in order to avoid aliasing the stellar
evolution. Smaller timesteps indicate rapidly changing orbital evolution
rates, where the timestep is reduced to maintain errors in the orbital
parameters below some fractional error threshold; we set this threshold to
$10^{-9}$, i.e. the update in accurate to one part in a billion.

###  4.4 Stellar model interpolation

Full coupling of stellar and orbital evolution requires the ability to solve
for the tidal response of a stellar model at arbitrary times. This means that
we either need to fully couple the MESA and GYRE-tides codes, or else devise
some interpolation scheme whereby the tidal response of the interpolated
stellar model closely matches that of the stellar model produced by the
stellar evolution code. We prefer the second approach since it permits
decoupling of the stellar evolution and orbital evolution segments of the
calculation. This interpolation scheme must adapt to changes in radius and
radial sampling of the star while preserving spatial resolution around key
features of the stellar interior such as the boundary between the convective
core and radiative envelope and near the stellar surface.

We adopt a simple linear interpolation between pairs of points in two stellar
models. This requires a sufficiently small timestep in the stellar model such
that perturbations in stellar structure are locally linear; we use a maximum
timestep of $10$ kyr in the stellar models. Interpolation errors can introduce
false peaks in the orbital evolution rates, leading to misinterpretation as
being due to resonances (Extended data Fig. 5a13a,b,c). However, interpolation
errors do not become significant until $dt>200$ kyr (Extended data Fig.
5d13d,e,f), with timesteps up to $1$ Myr retaining the same general shape of
the orbital evolution rates (Extended data Fig. 5a13a,b,c).

When the number of radial samples is equal between stellar models, the
interpolation is trivial since the radius provides a natural ordering and thus
natural coupling. When the number of radial samples differs, we use an optimal
transport based method for first finding the optimal matching between points
between the two curves, and then perform the interpolation with the coupled
points[chewi2021fast]. To apply optimal transport, we reinterpret the
cumulative number of points as a function of radius as a cumulative
distribution function. The optimal transport map between the two radial sample
profiles is then the map which matches quantiles. This allows single points to
match with multiple other points (i.e. allows mass splitting), which accounts
for the difference in sampling. We approximate one stellar model at the
sampling of the other by barycentric projection of the optimal transport map.
This produces radial profiles with the same sampling, admitting a natural
pairing of points, and thus allows simple linear interpolation to produce a
stellar model at any arbitrary time.

##  5 Captions of Figures

Figure 1. The rich tidal forcing spectrum of eccentric orbits Snapshots of the
stellar structure structure of an $M=1.36M_{\odot}$, $Z=0.02262$ star as
represented by a propagation diagram at three times: (a) Zero-age main
sequence, (b) midway through the main sequence, and (c) Red-edge of the main
sequence. The Blue line is the Brunt-Väisälä frequency $N$ and the red line is
the $l=2$ Lamb frequency $S_{2}$. g-modes can propagate in the blue region,
where $\sigma$ is below both $N$ and $S_{2}$. The orbital frequency is plotted
as a line near the bottom, with harmonics $0<k<50$ of the orbital frequency
shown adjacent to (a). Circular orbits excite only the $k=2$ harmonic, while
eccentric orbits spread energy across all $k$.

Figure 2. Episodic and wandering tidal migration by alternating resonance
locks. (Top) Rate-change of the orbital eccentricity induced by excitation of
quadrupolar (l=2) modes with (a) $m=0$ and (b) $m=2$ by a tidal perturber with
evolving orbital configuration. The host star has $M=1.36M_{\odot}$ and
$Z=0.02262$. The baseline $\Omega_{\mathrm{rot}}$ is set to
1.5$\Omega_{\mathrm{ps}}$, the crossover point for the tidal torque in Sun et
al. (2023). Red and blue points correspond to positive and negative $de/dt$,
respectively. Inset in (a) shows $5$ Myr of tidal evolution centered on a
resonance lock with the $m=0$, $k=12$ mode, and the inset in (b) shows $25$
Myr of tidal evolution with an isolated resonance lock with the $m=2$, $k=30$
mode followed by alternating resonance locks with the $m=2$, $k=32$, $22$,
$15$, $30$, and $19$ modes. (Bottom) Rate-change of the orbital eccentricity
decomposed into orbital harmonics $k=\sigma_{m,k}/\Omega_{\mathrm{orb}}$ for
(c) $m=0$ and (d) m=2. (e) Crossplot of the rate-change of the orbital
elements showing the diversity of different senses of orbital evolution.

Figure 3. Sensitivity of orbital evolution to stellar and orbital initial
conditions. (a) Rate-change of eccentricity across the HR diagram from the
zero-age main sequence to the red-edge of the main sequence for a fixed
orbital configuration, corresponding to the initial orbital state of the black
curve in Fig. 4b3b, for a fixed $Z=0.02$ and a range of masses. The colormap
is truncated at $\dot{e}\sim 10^{-10}/yr$ below at $\dot{e}\approx
10^{-10}/yr$ to emphasize the resonances. (b) Forward and time-reversed
orbital evolution trajectories with perturbations in each of the orbital and
stellar parameters, with colors and perturbation sizes given by the legend on
the right. The “x” ending each trajectory marks the terminal age main
sequence. The insets in (b) show a zoom of tidal evolution over the final 150
Myr of the main sequence.

Figure 4. Observability of Tidally Excited Oscillations (a) Evolution of the
spectrum of relative flux variations over wavelengths of $0.6-5$ $\mu$m.
(b,d,f) Stellar pulsation spectra at three times decomposed into $6$
wavelength bins between $0.6-5$ $\mu$m. Lines of different wavelength are all
on integer $k$, but are offset for visibility. (c,e,g) Wavelength-dependent
synthetic lightcurves corresponding to the pulsation spectra in (b,d,f). The
pulsating star has $M=1.36M_{\odot}$ and $Z=0.02262$.

J.B. and Z.L.d.B. acknowledge the National Science Foundation for supporting
this work through the Graduate Research Fellowship program under Grant No.
1745302 and the MIT Presidential Fellowship. The simulations presented in this
paper were performed on the Engaging cluster at the MGHPCC facility.

J.d.W. designed the study. J.B. developed the computational framework for the
study with notable support from M.S. and R.H.D.T regarding the stellar and
orbital modeling and from J.d.W. and Z.L.d.B. regarding exoplanetary context
and observables. All authors contributed to the manuscript writing, which was
led by J.B.

The authors declare that they have no competing financial interests.

We make use of the orbital configuration of HAT-P-2b and the stellar
properties of HAT-P-2[bonomo2017gaps, stassun2019revised]. The stellar tidal
response was converted to synthetic photometry using the CAP18 photospheric
grids[prieto2018collection].

This work makes use of the following publicly available codes:
MESA[paxton2010modules, paxton2013modules, paxton2015modules,
paxton2018modules, paxton2019modules, jermyn2023modules] for stellar
evolution, GYRE[townsend2013gyre, townsend2018angular, goldstein2020contour,
sun2023tides] for calculating stellar pulsation properties and orbital
evolution rates, and MSG[townsend2023msg] for converting the stellar tidal
response into synthetic photometric observables. Tools for coupling these
codes for live-planet simulations will be made publicly available at
https://github.com/jaredbryan881/orbev.

Correspondence and requests for materials should be addressed to J.B. (email:
<EMAIL_ADDRESS>or J.d.W. (email: jdewit@mit.edu).

## References

##  6 Tables and Figures for Methods

Extended Data Figure 1. Decomposition of tidal response of a slowly rotating
star into different modes. (Top) Rate-change of the orbital eccentricity
induced by excitation of quadrupolar (l=2) modes with $m=0,2$ in a host star
with $M=1.36M_{\odot}$ and $Z=0.02262$ by a tidal perturber with evolving
orbital configuration. The tidal migration in this case begins at
$t_{0}+1\sigma$. Red and blue points correspond to positive and negative
de/dt, respectively. (Bottom) Rate-change of the orbital eccentricity
decomposed into orbital harmonics $k=\sigma_{m,k}/\Omega_{\mathrm{orb}}$. The
baseline $\Omega_{\mathrm{rot}}$ is set to the value estimated for HAT-P-2 by
Bonomo et al. (2017). Colors correspond to those in Fig. 2.

Extended Data Figure 2 — 2. Accelerated tidal evolution by intermittent
resonance locking. Comparison between GYRE-tides and constant-Q tidal
evolution models. Forward and time-reversed orbital evolution trajectories
with perturbations in each of the orbital and stellar parameters. The baseline
$\Omega_{\mathrm{rot}}$ is set to the value estimated for HAT-P-2 by
Ref.[bonomo2017gaps]. The solid black line corresponds to the baseline GYRE-
tides model. Colored lines correspond to those in Fig. 4b3b. The dotted black
line corresponds to the best-fit constant-Q model from Jackson et al. (2008).
Gray lines are evenly spaced between $\pm 3\sigma$ in the given parameter
except for $Q_{\ast}$ and $Q_{p}$, who span the full range considered by
Jackson et al. (2008).

Extended Data Figure 3. Comparison of equilibrium and dynamical tidal
migration. Comparison of GYRE-tides models considering only equilibrium tides
to GYRE-tides models considering the full tidal response. The solid black line
corresponds to the baseline GYRE-tides model with the full tidal response. The
dotted black line corresponds to the GYRE-tides model with only the
equilibrium tidal response. Colored lines correspond to those in Fig. 3b, but
for only the equilibrium tidal response. The baseline $\Omega_{\mathrm{rot}}$
is set to the value estimated for HAT-P-2 by Ref.[bonomo2017gaps].

Extended Data Figure 4. Comparison of equilibrium and dynamical tidal
migration. Comparison of GYRE-tides models considering only equilibrium tides
to GYRE-tides models considering the full tidal response. The solid black line
corresponds to the baseline GYRE-tides model with the full tidal response. The
dotted black line corresponds to the GYRE-tides model with only the
equilibrium tidal response. Colored lines correspond to those in Fig. 3b, but
for only the equilibrium tidal response. The layout follows extended data Fig.
3. The baseline $\Omega_{\mathrm{rot}}$ is set to 1.5$\Omega_{\mathrm{ps}}$,
the crossover point for the tidal torque in Ref.[sun2023tides].

Extended Data Figure 5. Coherence of net orbital evolution rates for small
stellar perturbations. Rate-change of orbital eccentricity from 2.6 Gyr to the
end of the main sequence for a $Z=0.02262$ star with a range of mass
perturbations, using a base mass of $M=1.36M_{\odot}$. The orbital evolution
rate is calculated with a fixed-orbit tidal perturber with $e=0.575$,
$a=0.153$, and a pseudo synchronous stellar rotation rate.

Extended Data Figure 6. Limits of tidal evolution predictability for small
stellar perturbations. (a) Distribution of $\dot{e}$ and (b) distribution of
$\dot{a}$. (c) Distribution of $\dot{e}$ and (d) distribution of $\dot{a}$.
(e) Distribution of time needed to accumulate a $1\sigma$ deviation in $e$ and
(f) distribution of time needed to accumulate a $1\sigma$ deviation in $a$.
The different colors of distributions are given by the colorbar and denote the
perturbations in stellar mass.

Extended Data Figure 7. Coherence of net orbital evolution rates for small
stellar perturbations. Rate-change of orbital eccentricity from 2.6 Gyr to the
end of the main sequence for a $M=1.36M_{\odot}$ star with a range of
metallicity perturbations, using a base metallicity of $Z=0.02262$. The
orbital evolution rate is calculated with a fixed-orbit tidal perturber with
$e=0.575$, $a=0.153$, and a pseudo synchronous stellar rotation rate.

Extended Data Figure 8. Limits of tidal evolution predictability for small
stellar perturbations. (a) Distribution of $\dot{e}$ and (b) distribution of
$\dot{a}$. (c) Distribution of $\dot{e}$ and (d) distribution of $\dot{a}$.
(e) Distribution of time needed to accumulate a $1\sigma$ deviation in $e$ and
(f) distribution of time needed to accumulate a $1\sigma$ deviation in $a$.
The different colors of distributions are given by the colorbar and denote the
perturbations in stellar metallicity.

Extended Data Figure 9. Statistics of orbital evolution rates for
perturbations in orbital and stellar perturbations (left column) Distribution
of $\dot{e}$ and (right column) distribution of $\dot{a}$ corresponding to the
orbital evolution trajectories of fig. 3b in the manuscript. The panels follow
fig. 3b with trajectories differing by perturbations to the initial (a-b)
eccentricity, (c-d) semi-major axis, (e-f) stellar rotation rate, (g-h)
stellar age, (i-j) stellar mass, and (k-l) stellar metallicity.

Extended Data Figure 10. Limits of tidal evolution predictability within
observational bounds (left column) Distribution of time needed to accumulate a
$1\sigma$ deviation in e (right column) and distribution of time needed to
accumulate a $1\sigma$ deviation in a corresponding to the orbital evolution
trajectories of fig. 3b in the manuscript. The panels follow fig. 3b with
trajectories differing by perturbations to the initial (a-b) eccentricity,
(c-d) semi-major axis, (e-f) stellar rotation rate, (g-h) stellar age, (i-j)
stellar mass, and (k-l) stellar metallicity.

Extended Data Figure 3. 11. Coherence of single mode damping rates across
stellar models. HR diagram showing the contribution to $de/dt$ of the $m=2$,
$k=10$ orbital harmonic excited by tidal forcing with a fixed orbital
configuration and a fixed $Z=0.02$ for a range of masses.

Extended Data Figure 4. 12. Coherence of individual modes across stellar
models. HR diagrams of the amplitude of the (a) $k=10$ and (b) $k=30$ orbital
harmonics excited by tidal forcing of a host star with $M=1.36M_{\odot}$ and
$Z=0.02262$ with a fixed orbital configuration. Amplitudes are for an edge-on
view angle.

Extended Data Figure 5. 13. MESA model interpolation error. (a-b) Orbital
evolution rates, plotted in absolute value and on a log-scale, for
eccentricity $e$, semi-major axis $a$, and stellar angular momentum $J$. The
orbital configuration is fixed and the star evolves. The ground truth (black)
is given by the orbital evolution rates evaluated directly on the stellar
profiles produced by MESA ($dt_{gyre}=dt_{mesa}=10^{4}$ yr). Colored lines
correspond to orbital evolution rates evaluated on interpolated stellar models
(i.e. $dt_{gyre}<dt_{mesa}$). (d-f) L2 norm between the (log absolute value)
orbital evolution rates obtained in the ground truth model and for sparser
models. $dt_{mesa}$ is given on the x-axis and $dt_{gyre}$ is fixed at
$dt=10^{4}$ yr. All calculations use the reference properties of the HAT-P-2
system[bonomo2017gaps].