# Tidal Circularization of Binaries by Resonance Locking I: The Importance of
the Pre-Main-Sequence

J. J. Zanazzi Canadian Institute for Theoretical Astrophysics, University of
Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada  Yanqin Wu
Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON
M5S 3H4, Canada

###### Abstract

Although tidal dissipation in binary stars has been studied for over a
century, theoretical predictions have yet to match the observed properties of
binary populations. This work quantitatively examines the recent proposal of
tidal circularization by resonance locking, where tidal dissipation arises
from resonances between the star’s natural oscillation frequencies and
harmonics of the orbital frequency, and where resonances are ‘locked’ for an
extended period of time due to concurrent stellar evolution. We focus on tidal
resonances with axi-symmetric gravity-modes, and examine binaries with primary
masses from one to two solar masses. We find that orbital evolution via
resonance locking occurs primarily during the star’s pre-main-sequence phase,
with the main-sequence phase contributing negligibly. Resonance locking,
ignoring non-linearity, can circularize binaries with peri-centre distances
out to $\sim 10$ stellar radii, corresponding to circular periods of $\sim
4-6$ days. However, we find resonantly excited gravity-modes will become
nonlinear in stellar cores, which prevents them from reaching their full,
linear amplitudes. We estimate that such a ‘saturated resonance lock’ reduces
the circularization period by about a third, but resonance locking remains
much more effective than the cumulative actions of equilibrium tides. In a
companion paper, we examine recent binary data to compare against theory.

binary stars— tides – stellar oscillations – Young Stellar objects – stellar
evolution

##  1 Introduction

How tidal dissipation drives orbital and rotational evolution in binary star
systems has been studied for centuries (e.g. Darwin, 1879), but its final
resolution still eludes us (for reviews, see Mazeh, 2008; Ogilvie, 2014).

Tidal circularization is observed in binaries. Short period binaries tend to
be more circular (e.g. Mayor & Mermilliod, 1984; Duquennoy & Mayor, 1991;
Latham et al., 1992, 2002; Mathieu et al., 2004). The most systematic study to
date is that by Meibom & Mathieu (2005), who found that as stars age from Pre-
Main-Sequence (PMS) to late Main-Sequence (MS), the “characteristic
circularization period” ($P_{\rm circ}$, below which orbits are mostly
circular) generally increases with age. While binaries in young open clusters
(less than a few hundred million years) are circularized out to $P_{\rm
circ}\sim 8$ days, older clusters appear to exhibit longer circularization
periods, perhaps out to $P_{\rm circ}\sim 15$ days. Unfortunately, the data
used to establish such a trend are sparse. Each cluster contains only about a
dozen or so spectroscopic binaries, making the measurement of circularization
periods difficult and susceptible to sample variance. In a companion paper to
the current one (Zanazzi & Wu, 2020), we will re-examine the issue of
circularization periods using the much expanded binary samples from APOGEE and
Kepler.

A large number of theoretical studies have been undertaken to explain the
Meibom & Mathieu (2005) results, and all have difficulties accounting for the
long circularization periods in old binaries. In the following, we will review
a few of these theoretical studies that are influential and are relevant for
our study.

Zahn and collaboraters (Zahn, 1966, 1989; Zahn & Bouchet, 1989; Zahn, 2008)
developed the theory of equilibrium tides, where the star is assumed to react
hydrostatically to the binary tidal forcing, and a lag time (or angle) is
introduced to account for the internal dissipation of this response. The
dissipation leads to energy and angular momentum transfer between the binary
orbit and the star. This has been shown to be most effective in stars with
surface convection zones where turbulent eddies can provide kinematic
viscosity.

Zahn & Bouchet (1989) pointed out that, during the pre-main-sequence, this
process is more efficient because stars are more distended and have more
extensive convection zones. They found that binaries can be circularized out
to $P_{\rm circ}\sim 8$ days before they reach the main-sequence. While this
appears to explain the circularization data for young clusters, later studies
have called into question some of their assumptions (see Sec. 6.2 for a
discussion). In addition, Zahn’s equilbrium tide theory has difficulties
explaining the observed spin rates of some binary stars (Meibom et al., 2006;
Marilli et al., 2007; Lurie et al., 2017; Mazeh, 2008), an area not considered
in this work. However, the equilibrium tide theory does enjoy an undisputed
success in the case of evolved stars. Observations show that post-MS binaries
are circularized out to $\sim 100\ {\rm days}$ (e.g. Mayor & Mermilliod, 1984;
Bluhm et al., 2016; Price-Whelan et al., 2020), in accord with predictions
from Zahn’s theory (Verbunt & Phinney, 1995; Price-Whelan & Goodman, 2018).

Another parallel line of investigation is called the theory of dynamical
tides. While the equilibrium tide assumes the star reacts instantaneously to
the forcing, this theory accounts for the tidal excitation of dynamical modes.
In stars with radiative interiors, gravity-modes are the most relevant (e.g
Zahn, 1975; Savonije & Papaloizou, 1983; Terquem et al., 1998; Goodman &
Dickson, 1998), while inertial-modes may also play a role in stars that are
rotating sufficiently fast (e.g. Savonije & Papaloizou, 1997; Ogilvie & Lin,
2007; Ivanov & Papaloizou, 2007). In particular, when the frequency of tidal
forcing lies close to a natural oscillation frequency of the host (a tidal
resonance), the tidal response is strong and its associated dissipation brings
about rapid orbital evolution. However, because the resonance width is
normally very narrow, the binary orbit quickly shrinks through a resonance,
and the latter has little lasting dynamical impact (Terquem et al., 1998;
Goodman & Dickson, 1998).

To address this issue, a variant of the dynamical tide theory, ‘resonance
locking,’ has been developed. Savonije & Papaloizou (1983, 1984); Witte &
Savonije (1999b); Fuller et al. (2016); Fuller (2017) suggested that, when the
mode in resonance increases in frequency (due to evolution in either stellar
structure or stellar spin), in conjunction with the orbital frequency, the
resonance may be maintained for an extended period, during which the binary
continue to “surf” the resonance and experience significant orbital evolution.
It is argued to be important for both early type (Witte & Savonije, 1999b,
2001; Fuller & Lai, 2012a; Burkart et al., 2012) and late type (Witte &
Savonije, 2002; Savonije & Witte, 2002) main-sequence binaries.

A tidal resonance can be expressed as

$k\Omega\approx\sigma=\omega+m\Omega_{\star}\,,$ (1)

where $k\Omega$ is the $k$-order harmonic of the binary orbital frequency, and
$\sigma$ the mode frequency in the inertial frame. The latter is related to
the mode frequency in the star’s rotating frame ($\omega$), the mode’s
azimuthal quantum number ($m$), and the star’s spin rate ($\Omega_{\star}$).
Resonance locking requires that, as the orbit evolves,
$k\dot{\Omega}\approx\dot{\sigma}$. The latter can be caused by stellar
evolution ($\dot{\omega}$), or by changes in the host spin
($\dot{\Omega}_{\star}$), two possibilities addressed respectively by Fuller
et al. (2016); Fuller (2017) and Witte & Savonije (1999b); Fuller & Lai
(2012b).

Despite the promise of resonance locking, it is surprising that no
quantitative works have been performed to study its efficacy in circularizing
solar-type binaries, and to compare the predictions against observations. In
this work and an associated paper, we perform exactly these tasks. We study
resonance locking by tracking how a star’s oscillation frequencies evolve over
its lifetime (similar to Fuller et al., 2016; Fuller, 2017). For the first
time, we follow the evolution both during the pre-main-sequence and the main-
sequence. We find that the pre-main-sequence phase plays a disproportionately
large role, while the main-sequence phase contributes negligibly to tidal
evolution. This is in stark contrast to the trend reported by Meibom & Mathieu
(2005).

This result naturally stimulates us to re-examine the observed data. We are
aided by the large influx of spectrocopic and transiting data (e.g. Milliman
et al., 2014; Van Eylen et al., 2016; Triaud et al., 2017; Windemuth et al.,
2019; Price-Whelan et al., 2020), and we present our results in Paper II
(Zanazzi & Wu, 2020). With these two papers, we hope to establish that: 1)
resonance locking works primarily during the PMS for solar-type binaries, and
2) resonance locking during the PMS may help explain the observed
circularization periods.

While this paper is under review, the referee pointed out an important caveat
in our work. The resonantly excited gravity-modes in our calculations may well
reach such large amplitudes in the stellar radiative cores that the waves may
overturn, a process that has been discussed by Goodman & Dickson (1998);
Barker & Ogilvie (2010, 2011); Su et al. (2020). This is indeed born out by
our estimates. As many of the gravity-modes that can resonantly lock onto the
tidal perturbation are non-linearly damped and saturated at lower amplitudes
compared to our linear results, the efficacy of resonance locking is called
into question. We dedicate a new section to discuss this aspect. The overall
picture remains, at this moment, not as clear as we have hoped.

The structure of this paper is as follows. Section 2 reviews the properties of
gravity-modes (g-modes) in stars, the formalism describing their tidal
excitation. We also describe how we calculate the g-mode properties over the
host star’s lifetime. Section 3 discusses the conditions a binary must satisfy
to maintain resonance lock, as well as the binary’s orbital evolution under
resonance locking. Section 4 presents our results on the tidal evolution of
binary star populations under resonance locking. Section 6 confirms the un-
importance of equilibrium tides, discusses the theoretical uncertainties of
our analysis, and observational predictions. Section 7 summarizes our work.

##  2 Gravity Modes: Properties, Evolution, and Tidal Excitation

Before examining how the evolution of stellar binaries can drive rapid tidal
circularization through resonance locking, we will review the properties of
tidally excited gravity modes (g-modes). Under a tidal forcing
$-{\bm{\nabla}}U({\bm{r}},t)$ produced by a binary companion of mass
$M_{\star}^{\prime}$ orbiting at a semi-major axis $a$ and eccentricity $e$,
the displacement in the host star (with mass $M_{\star}$, radius $R_{\star}$
and rotation frequency $\Omega_{\star}$), ${\bm{\xi}}({\bm{r}},t)$ evolves as

$\frac{\partial^{2}{\bm{\xi}}}{\partial
t^{2}}+\mathcal{R}{\bm{\cdot}}{\bm{\xi}}={{\bm{f}}}-{\bm{\nabla}}U,$ (2)

where $\mathcal{R}$ is an operator acting on ${\bm{\xi}}$ due to restoring
force inside the star, and ${\bm{f}}$ represents any (weak) frictional force.
We restrict ourselves to gravity-modes, where the only contribution to
$\mathcal{R}$ is the buoyancy force. It is reasonable to ignore pressure-modes
(restored by pressure) as their eigenfrequencies lie well above the tidal
forcing frequency. We also assume that the tidal response is linear, that the
host star’s spin axis is aligned with the binary’s orbital plane, and that
tidal dissipation in the companion star is negligible. We further assume the
stellar spin rate is much lower than the g-mode oscillation frequencies
($\Omega_{\star}\ll\omega_{\alpha}$) so we can ignore the Coriolis force. We
also ignore another branch of oscillations, the so-called inertial modes. The
validity of these approximations will be discussed in Section 6.

The above equation simplifies for adiabatic, free oscillations into

${\color[rgb]{0,0,1}-}\omega_{\alpha}^{2}{\bm{\xi}}_{\alpha}+\mathcal{R}{\bm{\cdot}}{\bm{\xi}}_{\alpha}=0\,,$
(3)

where $\omega_{\alpha}$ is the frequency for eigen-mode $\alpha$. We can
express the effect of the weak friction force ${\bm{f}}$ by a damping rate,
using the so-called quasi-adiabatic approximation (see e.g. Friedman & Schutz,
1978a, b; Christensen-Dalsgaard, 2014)

$\gamma_{\alpha}=-\frac{{\rm
i}}{2\omega_{\alpha}}\frac{\langle{\bm{\xi}}_{\alpha}|{\bm{f}}\,{\rm e}^{{\rm
i}\omega_{\alpha}t}\rangle}{\langle{\bm{\xi}}_{\alpha}|{\bm{\xi}}_{\alpha}\rangle}\,.$
(4)

where ${\bm{\xi}}_{\alpha}$ is the adiabatic eigenfunction from equation (3).
Here, $\gamma_{\alpha}>0$ stands for the rate of energy (not amplitude)
dissipation, and $\gamma_{\alpha}\ll\omega_{\alpha}$. The bracket stands for
the volume integral,
$\langle{\bm{u}}|{\bm{v}}\rangle\equiv\int{\bm{u}}^{*}{\bm{\cdot}}{\bm{v}}\rho\text{d}V$,
and the factor of $2$ arises when one converts the energy damping rate to
amplitude damping rate.

Since the tidally forced perturbation is assumed to be linear, we can
decompose the total displacement into a sum over adiabatic eigenfunctions
${\bm{\xi}}_{\alpha}$,

${{\bm{\xi}}({\bm{r}},t)}=\sum_{\alpha}c_{\alpha}(t)\,{\rm e}^{-{\rm
i}\omega_{\alpha}t}\,{{\bm{\xi}}_{\alpha}({\bm{r}})}\,,$ (5)

where the mode amplitude $c_{\alpha}$ varies over a timescale much longer than
the mode frequency. As a result, we can approximately express

$\dot{\bm{\xi}}({\bm{r}},t)\simeq\sum_{\alpha}-{\rm
i}\omega_{\alpha}\,c_{\alpha}(t){\rm e}^{-{\rm
i}\omega_{\alpha}t}{\bm{\xi}}_{\alpha}({\bm{r}})\,.$ (6)

Substituting these above expressions into equation (2), we obtain the
amplitude equation that describes the evolution of $c_{\alpha}(t)$ (e.g.,
Schenk et al., 2002; Burkart et al., 2014; Fuller, 2017)

$\dot{c}_{\alpha}+({\rm
i}\omega_{\alpha}+\gamma_{\alpha})c_{\alpha}=\frac{{\rm
i}}{2\omega_{\alpha}}\frac{\langle{\bm{\xi}}_{\alpha}|-{\bm{\nabla}}U\,{\rm
e}^{{\rm
i}\omega_{\alpha}t}\rangle}{\langle{\bm{\xi}}_{\alpha}|{\bm{\xi}}_{\alpha}\rangle}\,.$
(7)

We review the derivation of equation (7) in Appendix B. In particular, when
the Coriolis force is also present, orthogonality between eigenmodes is called
into question. However, equations (7) and (4) remain valid.

The solution to this in-homogeneous equation describes the evolution of the
tidal response, and through it, the tidal evolution of the binary orbit. In
the following, we consider the quantities $\omega_{\alpha}$, $\gamma_{\alpha}$
and the tidal overlap integral (right-hand-side of eq. [7]), separately. In
particular, since we are interested in tidal evolution that proceeds in the
same timescale as stellar evolution, we pay special attention to how these
quantities change as the star ages.

In the following, we suppress the mode index $\alpha$ for notational
simplicity.

###  2.1 Gravity Mode Frequencies

This section calculates the g-mode eigen-frequencies for our stellar models
(see eq. [3]). Since this exercise is standard, we move most of the relevant
algebra to Appendix A, and only list a few relevant results here. Again, we
ignore the effects of rotation (Coriolis force and centrifugal force), except
the effect of mode frequency shifting by the rotating frame.

Of relevance for later discussions is the structure of a gravity mode: it
propagates (is oscillatory) in the star’s radiative region where its frequency
lies below the local Brunt-Väisälä frequency $N$, and is evanescent in the
convective region. For modes of high radial order ($n\gg 1$), the local
dispersion relation (eq. [A8]) can be integrated across the star to yield a
global dispersion relation

$\omega\simeq\frac{\sqrt{\ell(\ell+1)}}{\pi
n}\int_{0}^{R_{\star}}\frac{N(r)}{r}\text{d}r,$ (8)

where $\omega$ is mode frequency in a non-rotating frame, $n$ is the number of
zero nodes in its radial displacement, and $\ell$ refers to the spherical
degree when the angular dependence is described by the spherical harmonics
function ($Y_{\ell m}$).

To normalize the mode eigenfunction, we adopt (e.g. Weinberg et al., 2012)

$2\langle{\bm{\xi}}|{\bm{\xi}}\rangle\equiv
2\int_{0}^{R_{\star}}\left[\xi_{r}^{2}+\ell(\ell+1)\xi_{\perp}^{2}\right]\rho
r^{2}\text{d}r=1\,.$ (9)

Here, $\xi_{r}$ and $\xi_{\perp}$ are the radial and horizontal displacements,
respectively. Under this normalization, the displacements do not have the
dimension of length, and the mode energy (in the star’s rotating frame) is
$E=\epsilon=\omega^{2}$.

Figure 1:  Internal stratification inside stars, as they age through the pre-
main-sequence and main-sequence phases. Positive (solid) and negative (dashed)
values of the Brunt-Väisälä frequency $N^{2}$ are plotted as functions of
radius $r$, for the stellar masses $M_{\star}$ and ages $t$ indicated. G-modes
propagate where $\omega^{2}\leq N^{2}$.  Figure 2: The radial dependence of a
$\ell=2,n=11$ g-mode’s horizontal displacement ($\xi_{\perp}$) in a
$M_{\star}=1.0\,{\rm M}_{\odot}$ pre-main sequence star (age $1.04\times
10^{7}\,{\rm yrs}$) is shown in blue. The Brunt-Väisälä frequency $N^{2}$ is
shown in red, with positive (negative) values denoted by solid (dashed)
curves. The magnitude scale for the horizontal displacement $\xi_{\perp}$ is
arbitrary. The g-mode propagates where $\omega^{2}\leq N^{2}$, and it becomes
evanescent towards the bottom of the convection zone. Here, the g-mode
frequency is $\omega=0.115\,(G{\rm M}_{\odot}/{\rm
R}_{\odot}^{3})^{1/2}=0.99\,(2\pi/{\rm day})$. The displacement amplitude
remains roughly constant in the convection zone. As the mode propagates into
the stellar inner core, its amplitude rises and its wavelength shortens.

As a star ages and its stratification changes, its g-mode frequencies also
evolve. Figure 1 exhibits the radial profiles for the Brunt-Väisälä frequency
($N^{2}$) for three different stellar masses, at different stages of the
star’s Pre-Main Sequence (PMS) and Main Sequence (MS) evolution. During the
PMS, which lasts from a few to a few tens of Myrs (see Fig. 7), all stellar
types have radiative cores and convective envelopes, with the magnitudes of
$N^{2}$ roughly independent of the star’s mass. The magnitude of $N^{2}$ in
the radiative regions increases with time as the star contracts and its core
becomes more stably stratified. This causes the g-mode frequencies to increase
significantly during the PMS (eq. [8]).

After stars reach the MS, their structures diverge, with the surface
convection zone disappearing in the more massive models ($M_{\star}\gtrsim
1.4\ {\rm M}_{\odot}$). Over much of the MS, the magnitudes of $N^{2}$ in the
radiative regions rise only mildly (if at all), in contrast to the strong
increase in the PMS phase. This leads to a slight increase in g-mode
frequencies.

The radial structure of a representative g-mode is plotted in Figure 2, where
we see that the g-mode has an ‘evanescent tail’ in the surface convection zone
($r\gtrsim 0.6\,{\rm R}_{\odot}$), and its horizontal displacement rapidly
oscillates in the interior. The former is important for coupling to the smooth
tidal potential, while the latter is relevant for the dissipation of these
modes.

###  2.2 Damping Rates

Here, we evaluate how fast the mechanical energy from g-modes can be converted
into heat (mode damping rates). The dissipation causes the tidal displacement
to have a phase lag relative to the tidal forcing, exerting a back-reaction
torque on the perturber, which in turn drives the tidal evolution of the
binary system.

Dissipation of gravity-modes proceeds through three possible channels:
turbulent damping inside the convection zone, radiative damping caused by
photon diffusion, and last, the leakage of mode energy from the upper
atmosphere. The detailed expressions for these rates are enumerated in
Appendix C, where we also obtain order-of-magnitude estimates. Here, we
present these estimates and compare them against numerical results (Fig. 3).

Because a g-mode propagates almost adiabatically across the star, over (under)
densities are hot (cold) spots relative to the stellar background. Radiative
diffusion works to diffuse these temperature inhomogeneities, damping the mode
at a rate

$\gamma_{\rm
diff}\sim\frac{L_{\star}}{M_{\star}R_{\star}^{2}\omega^{2}}\Lambda_{\rm
core},$ (10)

where

$\Lambda_{\rm core}=\frac{M_{\star}R_{\star}^{2}}{\left({\rho
Hr^{4}}\right)_{r=r_{\rm c}}}$ (11)

parameterizes the characteristic mass concentration in the star’s radiative
core (evaluated at the characteristic radius $r_{\rm c}$, see Appendix C.1),
and $L_{\star}$ is the stellar luminosity. The central concentration enters
the estimate because the g-mode energy is mostly removed by radiative
diffusion in the stellar centre, where its wavelength is short and its
displacement is large. The parameter $\Lambda_{\rm core}$ increases with age
as the star becomes more centrally-concentrated. For a solar-mass star,
$\Lambda_{\rm core}\sim 10^{2}-10^{3}$ during the PMS, and $\sim
10^{3}-10^{4}$ during the MS. As a result, the damping rate increases as the
star evolves from the PMS to the MS. Examining the top panel of Figure 3, we
see estimate (10) does well at reproducing the scaling of $\gamma_{\rm diff}$
with frequency, as well as its dependence on stellar age. Higher frequency
modes have lower radial order and longer wavelengths. They suffer
significantly weaker radiative dissipation.

Turbulent damping arises from convective eddies interacting viscously with the
g-mode, causing a dissipation rate of order (see eq. [C20])

$\gamma_{\rm turb}\sim\frac{\ell^{2}}{3n\omega}\left(\frac{z_{\rm
cvz}^{2/3}L_{\star}^{2/3}}{R_{\star}^{2}M_{\rm cvz}^{2/3}}\right)\,,$ (12)

where $z_{\rm cvz}$ and $M_{\rm cvz}$ are the depth and mass of the surface
convection zone, respectively. The bottom panel of Figure 3 compares this
estimate with our numerical results and shows reasonable agreement. Compared
to the radiative damping rates, the turbulent damping rates decrease mildly as
the star ages, and do not depend sensitively on mode frequency.

The g-mode perturbation may propagate in a thin region above the photosphere
if the mode frequency lies below the so-called “acoustic cut-off frequency.”
Our calculation in Appendix C.3 shows this to be negligible compared to other
sources of damping, for g-modes with periods of order a few days.

Comparing the top and bottom panels of. Figure 3, we see radiative diffusion
dominates over turbulent viscosity ($\gamma_{\rm diff}\gtrsim\gamma_{\rm
turb}$) over the PMS and MS evolution of the star, with the exception of modes
of the lowest radial orders ($n\lesssim 10$). For all g-modes of concern, mode
dissipation is weak, relative to the spacing between consecutive g-modes
($\gamma\ll\omega/n$). Therefore, our assumption of individual, discrete modes
is appropriate (Goodman & Dickson, 1998).

Figure 3: G-mode damping rates as functions of mode frequencies ($f=1/P$), in
a sun-like star at the ages as indicated. Here $\ell=2$, and we consider
g-mode frequencies with radial nodes $n=3-200$ during the PMS ($t\lesssim
4\times 10^{7}{\rm yrs}$), and $n=50-500$ during the main-sequence ($t\gtrsim
4\times 10^{7}\,{\rm yrs}$). The top panel shows damping by radiative
diffusion, and the lower panel that by turbulent viscosity. The numerical
results are in solid lines, while our analytical estimates (eqs. [10] & [12])
are in dotted lines. We have adopted $\Lambda_{\rm core}=300$ ($3000$) during
the PMS (MS), respectively, when evaluating equation (10).  Figure 4:
Numerical results (solid lines) and analytic estimates (dotted lines, equation
28) of the tidal overlap integral $I_{\ell m}$ (eq. [19]), as functions of the
g-mode frequency $f=P^{-1}=\omega/2\pi$, for a Sun-like star at different
ages. Here $\ell=2$, and we consider g-modes with radial order $n=3-200$
during the PMS, ($t\lesssim 4\times 10^{7}{\rm yrs}$), and $n=50-500$ during
the MS. ($t\gtrsim 4\times 10^{7}\,{\rm yrs}$). We see that tidal coupling
decreases steeply with radial order, and at the same period, tidal coupling is
the largest during the PMS when the surface convection zone is deeper.

###  2.3 Tidal Excitation

The tidal potential from the binary companion acting on the host star may be
decomposed into spherical harmonics and a Fourier series as (Press &
Teukolsky, 1977)

$\displaystyle U({\bm{r}},t)=\ $
$\displaystyle-\frac{GM_{\star}^{\prime}}{a}\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}\left(\frac{r}{a}\right)^{\ell+1}W_{\ell
m}$ $\displaystyle\times\sum_{k=-\infty}^{\infty}X_{k\ell m}(e){\rm e}^{-{\rm
i}\omega_{km}t}{Y_{\ell m}(\theta,\varphi)}\,,$ (13)

where

$\displaystyle W_{\ell m}=\frac{4\pi}{2\ell+1}Y^{*}_{\ell
m}\left(\frac{\pi}{2},0\right)\,,$ (14) $\displaystyle X_{k\ell
m}(e)=\frac{1}{\pi}\int_{0}^{\pi}\frac{\text{d}E}{(1-e\cos E)^{\ell}}$
$\displaystyle\times\cos\left[k(E-e\sin
E)-2m\tan^{-1}\left(\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}\right)\right]\,,$
(15) $\displaystyle\omega_{km}=k\Omega-m\Omega_{\star}\,,$ (16)

and $\Omega=\sqrt{G(M_{\star}+M_{\star}^{\prime})/a^{3}}$ is the orbital
frequency/mean-motion. From now on, sums over $k$ will implicitly be from
$-\infty$ to $\infty$, $\ell$ from $2$ to $\infty$, and $m$ from $-\ell$ to
$\ell$, unless otherwise noted. We focus on the dominant $\ell=2$ component in
our numerical results, but the following discussion applies to all $\ell$
values.

Using this decomposition of the tidal potential, equation (7) is recast as

$\dot{c}+({\rm i}\omega+\gamma)c={\rm i}\omega\sum_{k\ell m}U_{k\ell m}{\rm
e}^{-{\rm i}\omega_{km}t}\,,$ (17)

where

$\displaystyle U_{k\ell m}$
$\displaystyle=\frac{M_{\star}^{\prime}}{M_{\star}}\left(\frac{R_{\star}}{a}\right)^{\ell+1}\left(\frac{E_{\star}}{\epsilon}\right)W_{\ell
m}I_{\ell m}X_{k\ell m}\,,$ (18) $\displaystyle I_{\ell m}$
$\displaystyle=\frac{\langle{\bm{\xi}}|{\bm{\nabla}}(r^{\ell}Y_{\ell
m})\rangle}{M_{\star}R_{\star}^{\ell}}=\frac{1}{M_{\star}R_{\star}^{\ell}}\int\delta\rho\,r^{\ell+2}\text{d}r\,.$
(19)

Here, $E_{\star}=GM_{\star}^{2}/R_{\star}$, and $\gamma$ is the mode damping
rate (Sec. 2.2).

The in-homogeneous solution to equation (17) is

$c(t)=\sum_{k\ell m}\frac{\omega}{\omega-\omega_{km}-{\rm i}\gamma}U_{k\ell
m}{\rm e}^{-{\rm i}\omega_{km}t}\,.$ (20)

This solution assumes that changes in the frequencies
$(\omega,\omega_{km},\gamma)$ are sufficiently slow that the mode amplitude
$c(t)$ can respond adiabatically. On the other hand, the dissipation needs to
damp the mode amplitude quickly, compared with variations in
$(\omega,\omega_{km},\gamma)$, so that no prior information is retained during
the evolution of $c(t)$. In Appendix D, we further justify the assumption of
adiabaticity for our problem, and argue for adiabadicity conditions contrary
to that suggested in Burkart et al. (2014), who implicitly assumed negligible
damping rates in their derivation. Adiabaticity guarantees the success of
resonance locking, a central topic in this work.

Near a resonance ($\omega\approx\omega_{km}$), one tidal component $(k,m)$
dominates, and solution (7) simplifies to

$c(t)\simeq C(t)\equiv\frac{\omega}{\omega-\omega_{km}-{\rm i}\gamma}U_{k\ell
m}{\rm e}^{-{\rm i}\omega_{km}t}.$ (21)

The maximum amplitude of $C$ occurs at exact resonance ($\omega=\omega_{km}$):

$C^{\rm max}={\rm i}\frac{\omega}{\gamma}U_{k\ell m}{\rm e}^{-{\rm
i}\omega_{km}t}.$ (22)

Each eigen-mode has energy $E_{\rm mode}$ and angular momentum $J_{\rm mode}$
in the inertial frame given by (Friedman & Schutz, 1978a, b; Burkart et al.,
2014; Fuller, 2017)

$E_{\rm mode}=\frac{\sigma\epsilon}{\omega}|c|^{2},\hskip 14.22636ptJ_{\rm
mode}=\frac{m\epsilon}{\omega}|c|^{2},$ (23)

where

$\sigma=\omega+m\Omega_{\star}$ (24)

is the mode frequency in the inertial frame. The mode energy and angular
momentum also attain their maximum values at exact resonance:

$E_{\rm mode}^{\rm max}=\frac{\epsilon\sigma\omega}{\gamma^{2}}|U_{k\ell
m}|^{2}\,,\hskip 14.22636ptJ_{\rm mode}^{\rm
max}=\frac{m\epsilon\omega}{\gamma^{2}}|U_{k\ell m}|^{2}\,.$ (25)

Numerically, the overlap integral $I_{\ell m}$ is delicate to calculate. To
confirm numerical results, we estimate its magnitude here. Because the
Eulerian density perturbation $\delta\rho$ oscillates rapidly in a mode’s
propagating region, most of the contribution to the overlap integral arises
from the evanescent zone, or the surface convection zone:111Here, we consider
low-frequency modes which propagate to the bottom of the convection zone.

$I_{\ell
m}=\frac{1}{M_{\star}R_{\star}^{\ell}}\int\delta\rho\,r^{\ell+2}\text{d}r\sim\frac{1}{M_{\star}}\int\left(\frac{\delta\rho}{\rho}\right)\rho
r^{2}dr\,.$ (26)

The Eulerian density perturbation is as expressed in equation (A4),

$\frac{\delta\rho}{\rho}=\frac{r\omega^{2}}{c_{\rm
s}^{2}}\xi_{\perp}+\frac{N^{2}}{g}\xi_{r}\simeq\frac{r\omega^{2}}{c_{\rm
s}^{2}}\xi_{\perp}\,,$ (27)

where $\xi_{\perp}(\xi_{r})$ is the tangential (radial) displacement (see
Appendix §A) and we have taken $N^{2}\approx 0$ in the convection zone. Since
$c_{s}^{2}\sim gz\propto z$ near the surface, the integrand is the largest at
the radiative-convective boundary, or

$I_{\ell m}\sim\bar{I}\frac{M_{\rm
cvz}}{M_{\star}}\left(\frac{r\omega^{2}\xi_{\perp}}{c_{\rm s}^{2}}\right)_{\rm
cvz}\,,$ (28)

where quantities within the brackets are to be evaluated at the bottom of the
convection zone. Furthermore, we provide an estimate for $\xi_{\perp}$ for our
adopted normalization in equation (A11). We have also inserted a numerical
fudge-factor $\bar{I}$. As Figure 4 shows, our numerical results agree with
this analytical estimate, including its dependencies on stellar structure and
mode frequency, as long as we adopt a constant fudge factor of $\bar{I}\sim
10^{-2}$ ($\bar{I}=0.007$ for this example).

This expression indicates that the overlap is reduced at higher mode order
$n$, and for stars with shallower convection zones. In particular, for stars
without a surface convection zone ($M_{\star}\gtrsim 1.4\,{\rm M}_{\odot}$ on
the MS), the overlap is entirely contributed by the thin region above the
first radial node (not accounted for here) and is very small. This makes
resonant excitation of g-modes difficult for early type main-sequence stars.

##  3 Resonance Locking – linear regime

Now we turn to the central concept in this work, resonance locking between
internal modes and the tidal potential. We review conditions for resonance
locking, and the resulting orbital evolution, building on previous works by
Fuller et al. (2016); Fuller (2017).

###  3.1 Conditions for Stable resonance locking

Figure 5:  Cartoon figure illustrating the conditions for resonance locking.
The timescales for orbital evolution ($t_{\rm orb}=\Omega/\dot{\Omega}$, red
curve) and mode evolution ($t_{\rm evol}=\sigma/\dot{\sigma}$, colored
horizontal lines for different cases) are plotted as functions of the binary
semi-major axis $a$. Under tidal dissipation, the orbit moves left-ward
(orbital frequency increases). Stable resonance locking occurs when: the mode
frequency also increases with time, and the evolution timescale at exact
resonance ($\sigma=k\Omega$, dotted line) satisfies $t_{\rm orb}^{\rm min}\leq
t_{\rm evol}$ . Case (a) is stable, while cases (b) and (c) are unstable.

When a mode is resonantly excited by the tidal potential, it converts orbital
energy into heat and causes orbital decay. The relevant rates of energy and
angular momentum change are (Weinberg et al., 2012; Burkart et al., 2014):

$\displaystyle\dot{E}_{\rm orb}$ $\displaystyle=-\sum_{k\ell
m}2k\Omega\epsilon U_{k\ell m}{\rm Im}\left(c\,{\rm e}^{{\rm
i}\omega_{km}t}\right),$ (29) $\displaystyle\dot{J}_{\rm orb}$
$\displaystyle=\sum_{k\ell m}2m\epsilon U_{k\ell m}{\rm Im}\left(c\,{\rm
e}^{{\rm i}\omega_{km}t}\right),$ (30) $\displaystyle\dot{J}_{\rm spin}$
$\displaystyle=\sum 2\gamma J_{\rm mode}.$ (31)

Here, $E_{\rm orb}=-GM_{\star}M_{\star}^{\prime}/(2a)$ is the orbital energy,
$J_{\rm orb}=-2\sqrt{1-e^{2}}E_{\rm orb}/\Omega$ is the orbital angular
momentum, $J_{\rm spin}=I_{\star}\Omega_{\star}$ is the host star’s spin
angular momentum, with $I_{\star}$ its moment of inertia. ${\rm Im}(X)$
denotes the imagninary part of complex number $X$. The above expressions
contain an implicit summation over all eign-modes, although typically only one
mode is important at any given time.

The most important mode (causing the largest changes in $\dot{E}_{\rm orb}$,
$\dot{J}_{\rm orb}$) is one that lies near resonance,
$\omega\approx\omega_{km}=k\Omega-m\Omega_{*}$, or in the inertial frame
$\sigma=\omega+m\Omega_{*}\approx k\Omega$. Under the action of this mode,
binary orbit evolves over a timescale $t_{\rm
orb}=\Omega/{\dot{\Omega}}=-\frac{2}{3}\,a/{\dot{a}}=-\frac{2}{3}\,[E_{\rm
orb}/(2\gamma E_{\rm mode})]$. Here, $E_{\rm mode}$ is the energy of the
excited mode (eq. [23]). At exact resonance ($\sigma=k\Omega$), the orbital
decay is the fastest, and its associated timescale is at a minimum,

$\displaystyle t_{\rm orb}^{\rm min}=-\frac{E_{\rm orb}}{3\gamma E_{\rm
mode}^{\rm max}}$
$\displaystyle=\frac{\gamma}{6\sigma\omega}\left(\frac{M_{\star}}{M_{\star}^{\prime}}\right)\left(\frac{a}{R_{\star}}\right)^{2\ell+1}\left(\frac{\epsilon}{E_{\star}}\right)W_{\ell
m}^{-2}I_{\ell m}^{-2}X_{k\ell m}^{-2}(e)\,.$ (32)

Resonance locking operates if $\dot{\sigma}\simeq k{\dot{\Omega}}\leq
k{\dot{\Omega}^{\rm max}}$. Let the timescale for mode frequency increase in
the inertial frame be

$t_{\rm evol}\equiv\sigma/{\dot{\sigma}}\,.$ (33)

In general, $\sigma$ is a function of both the stellar structure and the
star’s spin rate.222Even for zonal ($m=0$) modes, the Coriolis force may
affect mode frequencies (see §6). In this work, we ignore the spin
contribution (see discussions in §6), so $t_{\rm evol}$ is set entirely by
stellar evolution.333Tidal heating is insignificant for altering the stellar
structure, for the binary systems we consider. A stable resonance lock then
requires (Fig. 5)

$t_{\rm evol}\geq t_{\rm orb}^{\rm min}\,,$ (34)

where both the magnitudes and the signs of the timescales matter. Here,
$t_{\rm orb}^{\rm min}$ depends on both mode properties and binary properties.
It is shorter for modes that are more strongly coupled to the tidal potential,
are more weakly damped, and for binaries that have shorter orbital periods,
and higher eccentricities.

In addition, stable locking also requires $\gamma$ to be large enough for the
mode amplitude $c$ to evolve adiabatically (eq. [21]):

$\gamma^{-1}\ll t_{\rm orb}^{\rm min},t_{\rm evol}\,.$ (35)

Given our calculated mode damping rates (Fig. 3), this requirement is always
satisfied.

###  3.2 Orbital Evolution with resonance locking

Figure 6:  Evolution of binary eccentricity during inward migration for
different resonant locks to tidal frequency components
$\omega_{km}=k\Omega-m\Omega_{\star}$, using equation (41). Different lines
denote different initial eccentricities, with $e(0)=0.2$ (solid), $e(0)=0.5$
(dashed), and $e(0)=0.8$ (dotted).

As orbital energy is dissipated in the host star, the binary orbit shrinks.
Meanwhile, the total angular momentum is conserved. This section derives the
resultant evolution in eccentricity $e$ and semi-major axis $a$ during a
resonance lock.

Conservation of angular momentum requires $\dot{J}_{\rm mode}=-\dot{J}_{\rm
orb}-\dot{J}_{\rm spin}$. It may be shown using equations (22), (30), and (31)
that $\dot{J}_{\rm orb}\simeq-\dot{J}_{\rm spin}$ when
$\omega_{km}\simeq\omega$, so all binary orbital angular momentum is
transferred to the host star’s spin angular momentum near resonance
($\dot{J}_{\rm mode}\simeq 0$ when $\omega\simeq\omega_{km}$). Because $E_{\rm
mode}$ is related to $J_{\rm mode}$ via $J_{\rm mode}=-\frac{m}{\sigma}E_{\rm
mode}$, the evolution of orbital energy and angular momentum is related as
(Fuller, 2017)

$\dot{J}_{\rm orb}=-\frac{m}{\sigma}\dot{E}_{\rm orb}.$ (36)

While resonance locking operates,

$\displaystyle k\dot{\Omega}$
$\displaystyle\simeq\dot{\sigma}=\frac{\partial\sigma}{\partial\omega}\dot{\omega}+\frac{\partial\sigma}{\partial\Omega_{\star}}\dot{\Omega}_{\star}$
$\displaystyle=\frac{\omega}{t_{\rm
mode}}+mB\frac{\Omega_{\star}}{t_{\star}},$ (37)

where $t_{\star}=\Omega_{\star}/\dot{\Omega}_{\star}$ is the host star’s spin
up/down timescale, $t_{\rm mode}=\omega/\dot{\omega}$ is the mode frequency
evolution timescale (in the star’s rotating frame), while the parameter

$B=\frac{1}{m}\frac{\partial\sigma}{\partial\Omega_{\star}}$ (38)

takes into account the dependence of the inertial frame frequency $\sigma$ on
$\Omega_{\star}$ (see, e.g. Vick et al., 2019).

Since ${\dot{\Omega}}/{\Omega}={3\dot{E}_{\rm orb}}/({2E_{\rm orb}})$, the
binary orbital energy evolves as

$\frac{\dot{E}_{\rm orb}}{E_{\rm orb}}=\frac{2}{3}\left[\frac{1}{t_{\rm
mode}}+\frac{m}{k}\left(\frac{B}{t_{\star}}-\frac{1}{t_{\rm
mode}}\right)\frac{\Omega_{\star}}{\Omega}\right].$ (39)

To obtain an evolution equation for the binary eccentricity, we differentiate
$J_{\rm orb}=-2\sqrt{1-e^{2}}E_{\rm orb}/\Omega$ to obtain

$\frac{\text{d}e}{\text{d}t}=-\frac{1-e^{2}}{2e}\left(1-\frac{m}{k\sqrt{1-e^{2}}}\right)\frac{\dot{E}_{\rm
orb}}{E_{\rm orb}}.$ (40)

Expressing ${\dot{a}}/a=-{\dot{E}_{\rm orb}}/{E_{\rm orb}},$ we integrate
equation (40) to obtain the circularization track for binary evolution via
resonance locking:

$a\left(\sqrt{1-e^{2}}-\frac{m}{k}\right)^{2}=\text{constant}.$ (41)

The ratio $m/k$ is related to the ratio between the angular momentum and
energy a given mode extracts from the binary’s orbit. Notice equation (41) is
independent of $\Omega_{\star}$. For zonal modes ($m=0$), the mode angular
momentum is zero ($J_{\rm mode}=0$), and equation (41) yields
$a(1-e^{2})\propto J_{\rm orb}^{2}=\text{constant}$ for all $k$.

Some example tracks for different $(k,m)$ are plotted in Figure 6. We see that
when $m=0$, $e$ always decreases with time, while $m\neq 0$ evolution
frequently leads to more eccentric orbits (as is noted by, e.g. Witte &
Savonije, 1999b; Fuller, 2017).

Fortunately, as we argue in Appendix E, the relevant resonance lock is with
the zonal modes, while locking with the sectoral modes ($m\neq 0$) drives
negligible orbital evolution. This is related to the fact that the star’s
moment of inertia ($I_{\star}$) is much smaller than the binary’s moment of
inertia ($\mu a^{2}$, where
$\mu=M_{\star}M_{\star}^{\prime}/(M_{\star}+M_{\star}^{\prime})$ is the
reduced mass of the binary). The conservation of angular momentum then demands
that the spin evolution proceeds at a much higher rate than the orbital
evolution, $\dot{\Omega}_{\star}\sim(\mu
a^{2}/I_{\star})\dot{\Omega}\gg\dot{\Omega}$. Meanwhile, resonance lock
requires

$\dot{\omega}\simeq
k\dot{\Omega}-m\dot{\Omega}_{\star}\simeq-m\dot{\Omega}_{\star}\,,$ (42)

or the sectoral mode frequency in the star’s frame has to evolve rapidly. This
is difficult to satisfy (App. E.1), unless $k\gg m$.444Except for the highly
eccentric systems like the ‘heart-beat’ stars, one usually finds $k\not\gg m$.
We also consider when the binary spin has been nearly tidally synchronized
with the binary orbit, for arbitrary eccentricity (App. E.2). Because near
tidal synchronization of the stellar spin happens when $k\Omega\approx
m\Omega_{\star}$, resonant locks occur only on high-order g-modes
($\omega\approx k\Omega-m\Omega_{\star}\approx 0$ implies $n$ large, see eq.
[A9]). Hence, we find resonance locking is only possible for much higher-order
g-modes when compared with zonal modes (locks occur when $\omega\ll
k\Omega/m$). So from now on, we focus exclusively on zonal modes, where tidal
dissipation always lead to orbital circularization.

##  4 Application to binaries – linear study

This section applies the formalism developed above to test the effectiveness
of resonance locking with zonal g-modes ($m=0$) in circularizing binary
systems. We apply the Modules for Experiments in Stellar Astrophysics (MESA)
(Paxton et al., 2011, 2013, 2015, 2018, 2019) code to model the stellar
evolution. We use the pre_ms_to_wd test case within MESA to track stars with
masses from $M_{\star}=1.0,1.5$, to $2.0\,{\rm M}_{\odot}$, and from the
beginning of the PMS to near the end of the MS.

First, tidal circularization of stellar binaries by resonance locking usually
requires mode frequencies to increase with time. This is possible for gravity-
modes during both the PMS and the MS, with the PMS stars showing much more
marked changes. This condition, however, precludes inertial modes as good
candidates for resonance locking.

Second, resonance locking is more feasible for modes that are strongly coupled
to the tidal potential, and are weakly damped. Both of these are optimized
during the PMS, where the stars are less heavily stratified and g-modes of the
requisite periods have much lower radial orders.

Third, as we have argued above, sectoral modes ($m\neq 0$) are likely un-
important for the long-term evolution. So we can focus only on the $m=0$ tidal
potential. On the other hand, the relevant orbital harmonics ($k$) run from
$1$ to some small integers, with the maximum value depending on the orbital
eccentricity.

###  4.1 Regions of Stable resonance locking

Figure 7:  G-mode periods as functions of stellar age, for primary stars with
masses $M_{\star}$ as indicated. Stellar evolution generally leads to decrease
in g-mode periods, facilitating resonance locking with tidal forcing. The line
colors denote stable (blue, $t_{\rm orb}^{\min}<t_{\rm evol}$) and unstable
(red, $t_{\rm orb}^{\rm min}>t_{\rm evol}$) resonance locking, with the tide
excited by a low-mass companion ($0.3\,M_{\star}$) at the same orbital period
(component $(k,\ell,m)=(1,2,0)$) and at a modest eccentricity ($e=0.4$). We
consider g-modes with radial node orders $n=3-200$ during the PMS (to the left
the black dashed line), and $n=50-500$ during the MS (to the right). Resonance
locking occurs pervasively during the PMS, but is more sporadic during the MS.

We first numerically evaluate the regions of parameter space where stable
locking with g-modes can occur. Figure 7 depicts the evolution of g-mode
periods as the stars age, where we mark out modes that can stably lock with
the $k=1$ tidal potential, excited by a lower mass companion
($0.3\,M_{\star}$) at a moderate eccentricity ($e=0.4$). For this potential
(the dominant tidal term), locking requires $P=P_{\rm orb}=2\pi/\omega$.

To ascertain these numerical results, we estimate the boundary for resonance
locking, using the Solar model at $10^{7}$ yrs. To compare with Figure 7, we
only consider the $k=1$ orbital harmonics, and approximate $X_{120}(e)\simeq
3e/2$, valid for $e\ll 1$ (Weinberg et al., 2012). Binaries that allow
resonance locking satisfy equation (34) with the torque time taken from
equation (32). Using estimates in equations (8), (10), (28), and (A11), with
relevant parameters taken from our MESA model and $t_{\rm evol}$ set to
$10^{7}$ yrs, we find that such binaries satisfy

$P_{\rm orb}\lesssim P_{\rm lock}\sim 6\,{\rm
days}\,\,\left(\frac{e}{{0.4}}\right)^{6/31}\,.$ (43)

This is in reasonable agreement with results shown in the top panel of Fig. 7
when $e\ll 1$.

So our numerical calculations show resonance locking to be pervasive during
the PMS. Solar-type binaries with periods as long as $6$ days (at $e=0.4$) can
find stable g-mode locking, while this period increases to $\sim 10$ days for
A-type binaries ($2\,{\rm M}_{\odot}$). This is due to the better tidal
coupling and weaker mode damping for PMS stars – their g-mode at a given
frequency has a lower radial order. For stars we study, one expects
significant orbital evolution over the PMS phase.

In contrast, resonance locking during the MS is much less extensive. Figure 7
shows that only solar-type binaries can sustain significant resonance locking.
At the same orbital period, the g-modes of relevance are much higher in radial
order in the MS compared to that in PMS, leading to weaker tidal coupling.
This is even more exasperated for stars above the so-called ‘Kraft break’
($M_{\star}\gtrsim 1.4\,{\rm M}_{\odot}$, Kraft, 1967). These stars have no
surface convection zones, eliminating the evanescent tails through which
g-modes couple effectively to the tidal potential. Resonance locking, at least
with the dominant $\\{k,\ell\\}=\\{1,2\\}$ potential, is impossible for these
stars.

###  4.2 Orbital Evolution of Binaries

Figure 8:  Orbital evolution of three solar-type binaries during the PMS, as a
result of resonance locking. Starting from initial periods $P_{\rm
orb}(0)=4.0,6.6,8.5$ days, and initial eccentricities $e(0)=0.2,0.4,0.6$, the
orbital evolution proceeds as the solid curves. The symbols represent g-mode
periods (or their $k$ harmonics) during the same period, with the
corresponding ${n,k}$ values. Resonance locking requires that
$k\Omega=\omega$, or $P_{\rm orb}=kP$. Orbital evolution follows the mode
evolution while resonance locking is effective. The $n=9$ g-mode is
responsible for circularizing two of the systems, forced by two different
orbital harmonics ($k=1$ and $2$), while the middle binary system locks onto
two different g-modes in turn. After a few million years, resonance locking no
longer operates and the binary evolution is stalled.

We now use the above results to calculate a binary’s orbital evolution,
adopting a simplified approach. Instead of directly integrating equations
(39)-(40) to calculate the binary evolution, we simply look for g-modes with
frequencies $\omega_{\alpha}\simeq k\Omega$ at some point in the host star’s
evolution. We use equation (34) to determine if the $\\{k,\ell,m\\}$ component
of the tidal potential can lock onto mode $\alpha$. If so, we assume that the
orbital period evolves in sync with the mode evolution
($\Omega\simeq\omega_{\alpha}/k$), and the orbit circularizes at a constant
angular momentum (eq. [41]). If resonance locking is not possible, the orbit
does not evolve, while the star continues to age.

Figure 8 displays a few examples of this scheme. Three solar-type PMS
binaries, starting at 4.0, 6.6 and 8.5 day orbital periods and different
eccentricities, are found to stably lock onto a few g-modes of $n\sim 10$. The
tidal potential that is reponsible for locking has harmonics value $k=1$ or
$2$. As the stars contract along their PMS tracks, the g-mode periods shorten,
forcing the orbital period to shorten in step. The associated energy
dissipation brings about orbital circularization, which in turn reduces the
tidal forcing over time and increases the timescale for orbital evolution (eq.
[32]). For the cases considered here, resonance locking becomes ineffective
after only a few million years, but by then, significant orbital
circularization (and shrinking) has occurred. It is interesting to notice that
the binary that starts with the widest orbit ($P_{\rm orb}=8.5$ days) actually
progresses further than the next widest one ($P_{\rm orb}=6.6$ days).

Figure 8 also shows that multiple g-modes and orbital harmonics can work
together to circularize a binary. For the binary with an initial orbital
period of 6.6 days, after locking between the $n=17$ g-mode and the $k=1$
harmonic stops, the $\\{n,k\\}=\\{8,2\\}$ mode “picks-up” the orbit, and
causes a brief episode of additional circularization. There appears little
shortage of appropriate g-modes for resonance locking, at least during the PMS
stage.

Figure 9:  Evolutionary tracks (in lines) and final states (as magenta points)
of binary circularization by resonance locking, over the host star’s lifetime.
Starting from an initially uniform grid in $\log P_{\rm orb}-e$, cyan lines
trace binary evolution during the pre-main sequence, while magenta lines that
during the main-sequence – the latter are barely visible in most cases. Here,
we have included the $(\ell,m)=(2,0)$ tidal potential with orbital harmonics
$k=1,2,3,5,7,10,13,17,22,27,33$ (integers roughly uniformly spaced in
logarithmic values). Different panels are for the three host star masses
$M_{\star}$ as indicated. The dashed curves crudely describe the upper bounds
of binary orbits, as a result of resonance locking. They are orbits at
constant periaps ($a(1-e)$), with $P_{\rm circ}$ indicating the period where
$e=0$. The maximum period out to which linear resonance locking operates (eq.
[43]) is displayed by the dot-dashed blue line, while estimated range for
saturated resonance locking is shown by the black dotted line (eq. [46]).
Figure 10:  Same as the top panel of Figure 9 (host star mass
$M_{\star}=1.0\,{\rm M}_{\odot}$), except we integrate binary circularization
tracks starting at $t=4\times 10^{7}\,{\rm yrs}$ (i.e., zero-age-main-
sequence). While some circularization occurs, many binary systems remain
eccentric above the dashed curve, the rough boundary of circularization due to
PMS resonance locking.

In the following, we report the results of orbital evolution for a population
of binaries. Figure 9 displays the binary circularization tracks and final
orbits, starting from a regular grid of (logarithmic) orbital periods and
eccentricities, for three different host star masses $M_{\star}$.

It appears that resonance locking can effectively circularize all binaries
with periapses smaller than some critical values, corresponding to circular
orbits that range from $\sim 4$ days for solar-type to $\sim 6$ days for
A-type binaries.555In Zanazzi & Wu (2020), we provide a more rigorous measure
for the circularization periods. Binaries with periaps smaller than these
values should be circular. Even very wide orbits can circularize, as long as
they are sufficiently eccentric (e.g., $P_{\rm orb}=100$ days requires
$e\gtrsim 0.8$). This is due to the effective resonance locking between their
high-order tidal harmonics (large $k$) and medium order g-modes ($n\sim 10$).

We find that almost the entirety of orbital evolution in Figure 9 is solely
due to resonance locking during the PMS phase. Little to no resonance locking
is found to occur during the MS (with some minor exceptions in solar-type
binaries, see top panel of Fig. 9). To illustrate this issue more clearly, we
repeat our calculations for solar-type binaries, but without the benefits of
the PMS phase. This is shown in Figure 10. Many eccentric binaries with short
orbital periods ($P_{\rm orb}\lesssim 5\,{\rm days}$) remain eccentric, some
highly so.

The inefficiency of tidal circularization during the MS stage can be
attributed to two factors. First, the stronger stable stratification during
the MS (Fig. 1) shifts g-mode frequencies upward. So the resonant g-modes have
in general much higher radial orders. This leads to both weaker tidal coupling
(Fig. 4), and stronger damping (Fig. 3). Second, the MS stage is largely
static, compared to the PMS stage. G-mode periods vary only mildly, if at all
(Fig. 7). Many binaries could not find suitable g-modes to lock onto. Both
factors are more severe for more massive stars that have fully radiative
envelopes.

##  5 The Role of non-linearity

The calculations presented above assume that the gravity-waves remain linear
in its propagation and damping. However, we are prompted by the referee to
consider the role of non-linearity in resonance locking: if the tidal
displacement becomes non-linear, the wave can overturn the stratification and
dissipate energy rapidly. This physics complicates the resonant locking
picture substantially.

In stars with radiative interiors, as gravity-waves propagate towards the
centre, their wavelengths shorten and their amplitudes rise. This is most
extreme near the inner turning point (where $\omega\approx N$), where the
background stratification can be overturned by the large amplitude wave, akin
to the physics of Kelvin-Helmholtz instability (e.g. Goodman & Dickson, 1998;
Barker & Ogilvie, 2010, 2011; Barker, 2011; Weinberg et al., 2012). For stars
like the Sun, Goodman & Dickson (1998); Barker & Ogilvie (2010) showed that
gravity-waves excited by stellar or planetary companions can find themselves
to be nonlinear on their first travel inward to the core (we call this non-
resonant excitation). The non-linearity is stronger for longer period
companions, because they excite lower frequency waves, which have shorter
wavelengths, and can penetrate into more central regions. Given this, it is
natural to ask whether similar physics operates for pre-main-sequence
binaries, and if resonantly excited modes, with their large amplitudes, can in
fact exist.

###  5.1 Magnitudes of non-linearity

Non-linearity for non-resonantly excited waves is insignificant for pre-main-
sequence stars. Inserting values relevant for the solar model at 10 Myrs into
equation (14) of Goodman & Dickson (1998) (see also Ogilvie & Lin 2007; Barker
& Ogilvie 2010), we find that at the lower turning point ($r=r_{m}$), a
$({k,\ell})=(1,2)$ gravity-wave reaches a dimensionless amplitude

$|k_{r}\xi_{r}|_{r=r_{\rm m}}^{\rm non-res}\sim 6\times
10^{-4}\left(\frac{e}{{0.4}}\right)\,\left(\frac{P_{\rm orb}}{1\,{\rm
day}}\right)^{1/6}\,,$ (44)

or, non-resonant gravity-waves are linear for systems of interest (eq. [43]).
This non-linearity is about $10^{6}$ weaker than that in a Solar model (eq.
[34] in that work). This mostly results from the reduced stratification in PMS
cores, which moves the inner turning point outward (eq. [F2]).

In contrast, resonantly excited gravity-modes can reach much higher
amplitudes. In Appendix F.1, we estimate for the same gravity-wave but now
resonantly excited to find (eq. [F5])

$|k_{r}\xi_{r}|_{r=r_{\rm m}}^{\rm max}\sim 1.4\times
10^{4}\left(\frac{e}{0.4}\right)\left(\frac{1\,{\rm day}}{P_{\rm
orb}}\right)^{7/2}.$ (45)

For all systems of interest (eq. [43]), non-linearity is severe.

So in conclusion, during the PMS phase, tidally excited gravity-waves can
still form discrete modes (not mere running waves), but they could not reach
the amplitudes expected from (linear) resonant driving due to nonlinear
breaking. The efficacy of resonant locking is severely compromised by the
presence of non-linearity.

###  5.2 Saturated Resonant Lock?

What is the consequence of this compromise? While resonance lock still
operates, the amplitude a mode can reach is now limited to a lower value than
our linear estimates, leading to a reduced rate of tidal dissipation. This
leads us to consider a new scenario, saturated resonant lock.

We discuss this scenario in detail in Appendix F.2 . Here is a brief re-cap.
Consider a gravity-mode growing in amplitude due to resonant tidal forcing. It
reaches the threshold amplitude for breaking, in a time shorter than it takes
to reach the linear estimate. At this point, its energy may quickly be
converted into heat in the stellar core. Dividing this energy by the growth
time, we obtain an effective rate of energy dissipation. We then compare the
rate of orbital evolution with the rate of stellar evolution to find a new
condition for resonance lock. Inserting values appropriate for the Sun at 10
Myrs, we find that resonance locking is successful (satisfying eq. [34]) for
binaries with

$P_{\rm orb}\lesssim P_{\rm sat}\sim 3.6\,{\rm
days}\,\left(\frac{e}{{0.4}}\right)^{6/41}\,.$ (46)

This is $\sim 35\%$ shorter than the estimate based on linear resonance lock
(eq. 43), but still much beyond what equilibrium tides can accomplish (see
below).

Whether resonance lock in the presence of nonlinear saturation operates the
way we speculate here is unclear. A more careful look, involving perhaps
hydrodynamical simulations, is needed.

###  5.3 Nonlinear Wave Breaking during the main-sequence

Although the non-linearity of non-resonant dynamical tides is small during the
PMS phase, this is not so during the main-sequence, and Goodman & Dickson
(1998) argued that non-resonant dynamical tides should circularize two solar-
type binaries with ages $t$ out to orbital periods (see also Barker & Ogilvie,
2010; Barker, 2020)

$P_{\rm orb}\approx 6.4\left(\frac{t}{{10^{10}{\rm yrs}}}\right)^{1/7}\,{\rm
days}\,.$ (47)

Combining our picture of PMS resonant locking with the Goodman & Dickson
(1998) picture, we predict that solar-type binaries should be circularized out
to $\sim 3$ days by the end of the PMS, and exhibit a slow rise of the
circularization period out to $\sim 6$ days, as binary ages approach 10 Gyrs.
Stars above the Kraft break ($M_{\star}\gtrsim 1.4\,{\rm M}_{\odot}$), in
contrast, lack radiative interiors and should not display additional
circularization with age. Both these predictions can be tested.

##  6 Discussions

In the following, we first compare resonance locking with other proposed tidal
mechanisms. We then discuss uncertainties in our calculations and possible
improvements.

###  6.1 The Equivalent $Q$-value for Resonance Locking

Figure 11: Tidal circularization by resonance locking in solar-type binaries
(magenta dots, from the top panel of Fig. 9), are compared with
circularization expected from the equilibrium tide (dotted line), and that
obtained using a parametrized tidal-Q model (colored curves standing for
different $Q^{\prime}$ values). Tidal circularization by the equilibrium tide
has comparable contribution from both the PMS and the MS, and we have also
adopted for it an eccentricity dependence as that in equation (48). The
tidal-Q model is evaluated for the PMS phase, setting $R_{\star}=2\,{\rm
R}_{\odot}$. The curve with an effective $Q^{\prime}=(3/2)k_{2,\star}Q\sim
3\times 10^{6}$, corresponding to a circularization time that is comparable to
the PMS lifetime ($\tau_{\rm circ}\sim 4\times 10^{7}$ yrs), provides a good
resemblance to the results of resonance locking.

A constant time-lag model is often invoked to parametrize tidal dissipation.
This assumes the tidal bulge lags behind the perturber by an amount $\Delta
t_{\star}$ that is independent of the perturbing tidal frequency. This model
is convenient to integrate, explaining its outsized popularity, but does not
correspond to any particular physical tidal theory. This model is often also
parameterized by the so-called tidal Q-factor, where $Q=\Delta
t_{\star}\Omega$, or the modified tidal Q-factor,
$Q^{\prime}=\frac{3}{2}k_{2,\star}Q$, where $k_{2,\star}$ is the star’s tidal
Love number. Here, we infer the equivalent value of $Q^{\prime}$ for our
resonance locking model.

In this model, the circularization time (cf eq. 50) is expressed as (Hut,
1981)

$\tau_{\rm
circ}=\frac{Q^{\prime}}{11\Omega}\left(\frac{M_{\star}}{M_{\star}^{\prime}}\right)\left(\frac{a}{R_{\star}}\right)^{5}\left[\frac{F_{e}(e)G(e)}{F(e)}-\frac{18}{11}G_{e}(e)\right]^{-1}$
(48)

where the functions $F(e),G(e),F_{e}(e),G_{e}(e)$ (as defined in Hut, 1981)
account for the rise of tidal forcing with orbital eccentricity. We adopt a
Love number $k_{2,\star}\simeq 0.2$, the value for a fully convective star
(e.g., Chandrasekhar, 1939; Batygin & Adams, 2013), for the PMS phase.

The comparison in Fig. 11 shows that, to reproduce the amount of tidal
circularization achieved by resonance locking during the PMS, we require a
value of $Q\sim 3Q^{\prime}\sim 10^{7}$, equivalent to $\tau_{\rm circ}\sim
4\times 10^{7}$ yrs, the PMS lifetime for a solar-type star. We do not provide
an equivalent $Q$ value for the MS phase, because there is little
circularization then (Fig. 10). Moreover, nonlinear effects (§5) increases the
effective $Q$ value over that inferred here.

We caution that the tidal circularization by resonance locking is not properly
described by a constant-$Q$ (constant time lag) model. Rather it is an on-off
process that depends on whether the locking condition can be satisfied or not.

###  6.2 The un-importance of Equilibrium Tides

Another physically motivated model for tidal dissipation is the equilibrium
tide model (for reviews, see Zahn, 2008; Ogilvie, 2014). Instead of
considering dynamical response of the star, one assumes that the tidal
frequency is so low that the star adjusts instantaneously to the tidal
forcing. This describes the tidal response well when the binary lies far away
from any dynamical resonance, with the resulting displacements (equilibrium
tides) varying over global scales within the star, similar to the fundamental-
modes (see, e.g. Terquem et al., 1998) This minimizes their dissipation via
radiative diffusion. In stars with surface convection zones, turbulent
viscosity is the dominant dissipation for these waves. Dissipation produces a
phase lag between the tidal forcing and the response, which then leads to
tidal circularization and synchronization.

Here, we calculate the impacts of equilibrium tides, both during the MS and
the PMS phases. This calculation shows that the equilibrium tide is
ineffective in both stages, compared to resonance locking.

In the limit of low eccentricity, the equilibrium tidal theory of Zahn (1977,
1978, 1989) results in the orbital evolution

$\displaystyle\frac{1}{a}\frac{\text{d}a}{\text{d}t}$
$\displaystyle=-\frac{2}{\tau_{\rm circ}}e^{2},$ (49)
$\displaystyle\frac{1}{e}\frac{\text{d}e}{\text{d}t}$
$\displaystyle=-\frac{1}{\tau_{\rm circ}},$ (50)

where the circularization time $\tau_{\rm circ}$ is the critical parameter
that determines the efficiency of eccentricity damping (eq. [21] of Zahn
1989). The derivation of equations (49) and (50) assumes $e\ll 1$.
Circularization is expected to be faster for high eccentricity systems (Hut,
1981; Vick & Lai, 2020). The timescale $\tau_{\rm circ}$ is recast by Goodman
& Oh (1997) as their equation (45). We re-evaluate that equation for a binary
with mass ratio $q=0.3$. Considering only tides in the primary star, this is

$\tau_{\rm
circ}^{-1}=1684\,\left(\frac{{R_{\star}}}{a}\right)^{8}\,\frac{1}{{M_{\star}R_{\star}^{8}}}\,{\int\rho\nu_{T}r^{8}\text{d}r}\,,$
(51)

where $\nu_{T}$ is the turbulent viscosity appropriate for the equilibrium
tide. Over most of the convection zone, the tidal period ($P_{\rm tide}$) is
short compared to the convection turn-over time. Adopting a reduction in
turbulent viscosity that scales as the ratio of the two timescales to the
squared power (Goldreich & Nicholson, 1977; Goodman & Oh, 1997; Ogilvie &
Lesur, 2012; Duguid et al., 2020; Vidal & Barker, 2020), Goodman & Oh (1997)
expressed, for a convection zone that can be approximated as an $n=1.5$
polytrope,

$\rho\nu_{T}=\frac{9}{{128}}\left(\frac{P_{\rm orb}}{2\pi
t_{f}}\right)^{2}\frac{{M_{\star}}}{{4\pi R_{\star}t_{f}}}x^{-3}(1-x)^{-1}\,,$
(52)

where $x=r/R_{\star}$, and we have taken $P_{\rm tide}=P_{\rm orb}/2$. Here,
$t_{f}=(M_{\star}R_{\star}^{2}/L_{\star})^{1/3}$ is an estimate for the global
convection turn-over time. For the Sun, $t_{f}\sim 0.45$ yrs.

Assuming a thin surface convection zone, we can evaluate the integral in
equation (51) and write

$\tau_{\rm circ}^{-1}\approx
0.05\left(\frac{{R_{\star}}}{a}\right)^{8}\,\frac{P_{\rm
orb}^{2}}{t_{f}^{3}}\,.$ (53)

In the following, we enhance the dissipation rate by a factor of $10$ to
account for the possibility of a fully convective star, and for our perhaps
overly severe reduction of the turbulent viscosity (Fig. 2 of Goodman & Oh,
1997). Expressing the stellar dynamical time as $t_{\rm
dyn}=(GM_{\star}/R_{\star}^{3})^{-1/2}$, the final expression is then,

$\displaystyle\tau_{\rm circ}$ $\displaystyle\approx$ $\displaystyle 1.5\times
10^{-4}\,t_{f}^{3}\,t_{\rm dyn}^{-16/3}\,P_{\rm orb}^{10/3}\,$ (54)
$\displaystyle=$ $\displaystyle 7\times 10^{12}\,{\rm yrs}\left(\frac{P_{\rm
orb}}{10\,{\rm days}}\right)^{10/3}$
$\displaystyle\times\left(\frac{M_{\star}}{{\rm
M}_{\odot}}\right)^{11/3}\,\left(\frac{R_{\star}}{{\rm
R}_{\odot}}\right)^{-6}\,\left({L_{\star}}\over{{\rm
L}_{\odot}}\right)^{-1}\,.$

To estimate the maximum orbital period that can become circularized over the
main-sequence, we follow Goodman & Oh (1997) by setting $\tau_{\rm circ}\sim
1/3\times t_{\rm MS}$, the main-sequence lifetime. For the Sun, this yields a
very tight orbit, $P_{\rm orb}\sim 1.3$ days. This is about twice shorter than
the $2.5$ days estimate in Goodman & Oh (1997), because of all the reduction
factors we have included above.

As stars are puffier during their PMS phase, action during that stage can be
more significant (Zahn, 1989). PMS evolution of the Sun (e.g., Gough, 1980),
can be approximately fitted as

$\displaystyle R_{\star}(t)$ $\displaystyle\approx$ $\displaystyle
1.26\,R_{\odot}\left({t\over{10^{7}{\rm yrs}}}\right)^{-0.3}\,,$
$\displaystyle L_{\star}(t)$ $\displaystyle\approx$ $\displaystyle
0.4\,L_{\odot}\left({t\over{10^{7}{\rm yrs}}}\right)^{-0.6}\,,$ (55)

valid from $\sim 3\times 10^{5}$ yrs to $\sim 3\times 10^{7}$ yrs. To compare
the relative importance of the PMS action, we integrate equation (50) to
obtain the relative contributions to circularization

$\displaystyle{{{\Delta e}|_{\rm PMS}}\over{\Delta e|_{\rm MS}}}$
$\displaystyle\approx{{1.1\times 10^{7}{\rm yrs}}\over{t_{\rm
MS}}}\,\left({{t_{0}}\over{10^{7}{\rm yrs}}}\right)^{-1.4}$ $\displaystyle\sim
0.6\left({{t_{0}}\over{10^{5}{\rm yrs}}}\right)^{-1.4}\,,$ (56)

where $t_{0}$ is the start of our PMS stage, and we have evaluated the
expression for a solar-type star. So while the rate of tidal damping is
enhanced during the PMS, its short duration limits its contribution to the
overall effect.666Compared to the PMS, a star on the Hayashi track is even
fluffier. And it has a higher luminosity powered by accretion. The impact of
this stage is unclear. This conclusion (also see Goodman & Oh, 1997) differs
from that obtained in Zahn (1989); Zahn & Bouchet (1989), mainly because we
have adopted the sharper viscosity reduction, as is advocated by Goldreich &
Nicholson (1977); Goodman & Oh (1997).

A graphical comparison of the two models, the equilibrium tide model and the
resonance lock model, is presented in Fig. 11. The equilibrium tide is un-
important at all times.

###  6.3 Examining the Assumptions

In this work, we have adopted a simple approach to study resonance locking. We
first calculate the properties of the gravity-modes in a given star (mass,
age), then evaluate whether any component (quantum number $\\{k,\ell,m=0\\}$)
of the tidal force can successfully lock on to a single g-mode (eq. [34]). If
so, the orbital frequency $\Omega(t)$ is assumed to evolve in step with the
frequency of this mode $\omega(t)$, $\Omega\simeq\omega/k$, as stellar
evolution modifies the mode frequency. The orbit is assumed to circularize at
a constant angular momentum (eq. [41]). We assume the orbital evolution stalls
when there is no appropriate resonance locking. We have justified the last
point by arguing that the equilibrium tide is ineffective. Dynamical tides
without resonance locking has also previously been shown to be irrelevant for
solar-type stars (Terquem et al., 1998; Goodman & Dickson, 1998).

In these calculations, we have only considered locking with zonal ($m=0$)
g-modes, and have ignored the effects of Coriolis force and nonlinear wave
interactions. We have adopted a form for the turbulent viscosity that is
strongly reduced for slow tides. Here, we discuss these simplifications and
assumptions in turn.

First, we are justified to consider locking with only one mode at any given
time. G-modes are sparse, and the resonance width
($|\delta\omega|\sim|\omega-\omega_{km}|\sim|\gamma|$) is much smaller than
the frequency difference between two consecutive g-modes
($|\Delta\omega|\sim\omega/n$ ; see Fig. 3).

We consider only the zonal modes because of the mis-match in orbital and spin
angular momenta (§ 3.2). The adoption of zonal modes means we can forgo the
spin evolution.777In fact, our approach is fully self-consistent when the star
is not spinning, since zonal modes then do not carry any angular momentum.
This simplifies the evolution significantly, and it ensures that the resonance
locking is stable – Burkart et al. (2014) showed that, if both energy and
angular momentum are transferred between the binary orbit and the stellar spin
(e.g., by $m\neq 0$ modes or the equilibrium tide), it may be difficult to
maintain resonance locking. We further justify our choice of zonal modes by
investigating the inefficiency of sectoral ($m\neq 0$) modes in Appendix E.

However, our approach assumes that the frequencies of zonal modes are
independent of the spin rate. This is only correct if
$\Omega_{\star}\ll\omega$. Unfortunately, gravity-modes of interest have
periods up to a few days (Fig. 7), comparable in value to the spin periods of
T Tauri stars (e.g. Bouvier, 2013), or of stars that are spin synchronized
with the orbit. For modes in these stars, the Coriolis force provides a
substantial restoring force, and their frequencies (and structure) become
dependent on the spin rate (for an example, see Fig. 2 of Bildsten et al.,
1996), in addition to their dependencies on stellar structure. This could
drive the system into or away from resonance locking (e.g. Lainey et al.,
2020).

We have opted not to include the Coriolis effect in this work, for two reasons
that are both to do with our inadequacies. First, including it will make
g-modes much more analytically un-tractable – while one can use the
‘traditional approximation’ to calculate g-modes for fully stratified stars
(when $\Omega_{*}\ll N$, Chapman & Lindzen, 1970; Bildsten et al., 1996), our
stars typically have convective outer shells. Buoyancy-dominated g-modes
become inertial-gravity modes and there is currently no good technique to
tackle them (but see attempts in e.g. Lai, 1997; Witte & Savonije, 1999a;
Savonije & Witte, 2002; Ogilvie & Lin, 2004; Lai & Wu, 2006; Ogilvie & Lin,
2007; Wu, 2005a, b; Ogilvie, 2013; Lin & Ogilvie, 2017). Second, as the mode
frequency depends on the spin rate, we have to account for all physical
processes that change the spin, including magnetic braking, stellar
contraction, tidal torques from both the equilibrium tides and other dynamical
tides. We also have to consider the re-distribution of angular momentum inside
stars. This is beyond our current capability. The neglect of Coriolis force
remains one of our most serious caveats.

Damping of g-modes directly affects resonance locking. If the modes are more
strongly damped, the tidal torque will be smaller and resonance locks may
break. We have experimented with enhancing the damping by a factor of $10$,
and find that the circularization period moves inward by about a day.

We have briefly considered the role of nonlinear wave breaking on our results.
In stars with radiative interiors, mode interactions occur predominantly near
the inner cores where the gravity-wave wavelengths are short and the
amplitudes are large. We find that, even in the PMS phase, resonantly excited
gravity-modes invariably reach nonlinear amplitudes, though they remain
discrete modes. The presence of non-linearity reduces the efficacy of
resonance locking. Our crude estimate reduces the critical binary period by
$\sim 35\%$ (or separation by $\sim 23\%$), compared to when non-linearity is
ignored. A more thorough investigation is needed.

Concurring with Goodman & Oh (1997), we find that the effects of equilibrium
tide to be negligible in Sun-like stars, during both the PMS and MS phases.
This conclusion depends on the turbulent viscosity form that we have adopted
(Goldreich & Nicholson, 1977; Goodman & Oh, 1997). If the form advocated by
Zahn (1966) is adopted instead, with a less severe reduction for the cases
where convective turn-over is much slower than the tide, Zahn & Bouchet (1989)
found that binaries out to $P_{\rm orb}\approx 8$ days can be circularized
during the PMS phase. This would then make the equilibrium tide comparable in
importance to resonance locking. If so, an intricate inter-play between spin
and orbital evolution may result (e.g. Witte & Savonije, 2002; Savonije &
Witte, 2002; Fuller & Lai, 2012a; Burkart et al., 2012). As the equilibrium
tide drives the star towards spin synchronization (Hut, 1981), resonance
locking of g-modes may be modified or disrupted as the frequencies of these
modes are altered by the spin. This is made more complicated by the fact that
the equilibrium tide deposits its energy and therefore angular momentum
primarily in the outer convective shell. In contrast, much of g-mode
properties (frequency, structure and damping) are determined in the stellar
interior. It is yet unclear how the two regions communicate.

There are a few situations where dynamical processes different what we study
here may enter the play. For more massive stars without surface convection
zones, tidal excitation and dissipation of gravity-waves behave very
differently. For instance, the absence of an evanescent tail weakens the tidal
coupling of g-modes. Moreover, dissipation of gravity waves may be dominated
by the so-called “critical layers” (Su et al., 2020). In another scenario,
relevant for highly eccentric stellar binaries ($e\gtrsim 0.9$), f-modes may
be stochastically excited at each periastron passage (e.g. Press & Teukolsky,
1977; Mardling, 1995a, b; Vick & Lai, 2018; Wu, 2018). The diffusive growth of
f-modes can efficiently extract orbital energy and circularizes the binary’s
orbit rapidly. Wu (2018) estimated an effective $Q\sim 1$ during this phase.

###  6.4 Observational Comparison

We argue in this work that most of the orbital evolution has finished by the
time stars enter the main-sequence. Moreover, we find that the circularization
periods only extend to $\sim 4$ days or shorter for solar-type binaries. This
runs counter to the results by Meibom & Mathieu (2005); Meibom et al. (2006),
where they showed that the longest circular orbits appear to be wider in more
evolved stars and can extend to $\sim 15$ days. In Paper II, we will present a
comprehensive analysis of current data, and compare it against our theory. We
suggest that the data are more compatible with shorter circularization
periods. A brief comparison is shown in Figure 12, showing observed binaries,
alongside the prediction from linear resonance lock during the PMS. Non-
linearity during the PMS and the MS will modify these predictions and are not
yet included in this prediction.

Figure 12: A preview of observation vs. theory. More details see Paper II.
Here, the colored dots represent observed binaries that contain G-star
primaries, with green circles being spectroscopic binaries from Price-Whelan
et al. (2020), and blue triangles eclipsing binaries from Windemuth et al.
(2019). The upside-down cyan triangles mark eclipsing binaries that display
eclipse-timing variations and so may be in triple systems. We do not display
data from Price-Whelan et al. (2020) when $P_{\rm orb}<5\ {\rm day}$, because
the cadence of the observations likely lead to an over-estimate of $e$ at
short $P_{\rm orb}$ values (see Paper II for further discussion). The two
theory lines are taken from Figure 11, where the line for resonance locking
has not accounted for non-linearity.

Observant readers may notice that Figure 12 seems to show an excess of
circular orbits out to $\sim$10 days. In paper II, we will show this excess is
statisitically significant, and is not explainable within our framework of
resonance locking. Additional mechanisms, such as formation of short-period
nearly-circular binaries via disk migration (e.g. Kratter et al., 2010; Moe &
Kratter, 2018; Tokovinin & Moe, 2020), may need to be invoked.

Another way to test resonance locking during the pre-main sequence is to
observe resonance locking in action. Pulsations of resonantly forced gravity-
modes may become observable in pre-main-sequence binaries. Resonance locking
theory can provide clear predictions on their amplitudes and phases (Burkart
et al., 2012; Fuller & Lai, 2012b; Fuller, 2017; Fuller et al., 2017).

We have skirted the issue of spin evolution in this work, by invoking tidal
forcing of zonal g-modes. As such, we can not make predictions on stellar
spin. Observationally, measurements of $P_{\rm rot}$ for short-period binaries
show that most binaries with $P_{\rm orb}\lesssim 10\,{\rm days}$ have nearly
synchronous rotation periods, but with some notable exceptions (e.g. Meibom et
al., 2006; Marilli et al., 2007; Lurie et al., 2017). It is possible that
resonance locking is instrumental in the spin evolution, but further work is
required to address this issue.

##  7 Conclusions

We have investigated the process of resonance locking to circularize binary
systems. Tidal excitation of dynamical tides transfers orbital energy into
mode energy, which is subsequently deposited in the stellar centre as heat.
Where resonance between the tidal potential and an internal mode can be
sustained, the orbit evolves over the same timescale as stellar evolution. We
find this to be the case for pre-main-sequence binaries out to a few days.
Main-sequence binaries, in contrast, rarely enjoy the benefit of resonance
locking, and we find negligible tidal evolution during this phase of the
stellar evolution.

Stars both above and below the Kraft break harbour surface convection zones
during the PMS. Our calculations show that they can both be circularized to
similar distances. Ignoring non-linearity, we find that as these binaries exit
the PMS, their orbits should have evolved to fall under a curve of $a(1-e)\geq
a_{\rm crit}$, where $a_{\rm crit}$ corresponds to a period of $4-6$ days
(Fig. 9). Expressed in units of MS stellar radii, we have

$a_{\rm crit}/R_{\star}\approx 10\,.$ (57)

In contrast, equilibrium tide can only circularize these stars out to $\sim
1.3$ day periods, both during the PMS and the MS, or $a_{\rm
crit}/R_{\star}\sim 5.5$.

Non-linearity in the stellar cores, however, is substantial for all models of
interest. Our preliminary attempt to estimate its impact finds that resonance
locking should still operate during the PMS, but with a lower efficiency. This
process which we call ‘saturated resonance locking’ may circularize binaries
out to $a_{\rm crit}/R_{\star}\sim 8$ (equivalent to $\sim 35\%$ reduction in
period). More work is warranted.

These circularization limits seem to fall far short of what have been reported
in the literature ($\sim 10-18$ days for old stellar clusters and halo stars,
Meibom & Mathieu, 2005; Meibom et al., 2006). However, as we will argue in a
companion paper, re-examining these observational claims lead to very
different conclusions.

Despite our progress, there remains a number of theoretical uncertainties in
the physics of resonance locking. We have shown that zonal gravity-modes
(azimuthal number $m=0$) play the most important role in resonance locking.
However, the frequencies and structure of these modes are affected by stellar
spin, when the latter is sufficiently fast. So the spin evolution, ignored in
this paper, may upset resonance locking in ways not accounted for here.
Nonlinear mode interactions, which may enhance mode dissipation, could also
affect our conclusions.

We thank our anonymous referee for insightful comments. We also thank Eugene
Chiang, Jim Fuller, Kaitlin Kratter, Dong Lai, Maxwell Moe, Gordon Ogilvie,
and Michelle Vick for helpful conversations. YW acknowledges NSERC for
research funding.

##  Appendix A Gravity-modes: Notations and Properties

This section reviews the main properties of g-mode oscillations (see also e.g.
Unno et al., 1979; Christensen-Dalsgaard, 2014, for reviews). Because g-modes
rapidly propigate across the star before energy can be diffused, they are
approximatly adabatic, or $\Delta p=c_{\rm s}^{2}\Delta\rho$. Here $\Delta p$
and $\Delta\rho$ are the Lagrangian perturbed pressure and density, while
$c_{\rm s}^{2}=(\text{d}p/\text{d}\rho)_{S}$ is the adiabatic sound-speed
squared. We cut off the outer radius of the star where non-adiabatic effects
are important, and set $R_{\star}$ to be 97% of the outer radius MESA
calculates ($R_{\star}=0.97\,R_{\star,{\rm MESA}}$).

To calculate g-modes in the star, we perturb the star’s equilibrium state
$X(r)$ by Eularian perturbations $\delta X(r)Y_{\ell m}(\theta,\varphi){\rm
e}^{-{\rm i}\omega t}$, where $Y_{\ell m}$ are spherical harmonics. We
decompose the angular dependence of the displacement vector
${\bm{\xi}}({\bm{r}},t)$ into vector spherical harmonics:

${\bm{\xi}}({\bm{r}},t)=\xi_{r}(r){\bm{\hat{r}}}Y_{\ell m}(\theta,\varphi){\rm
e}^{-{\rm i}\omega t}+\xi_{\perp}(r)\,\left[r{\bm{\nabla}}Y_{\ell
m}(\theta,\varphi)\right]{\rm e}^{-{\rm i}\omega t}.$ (A1)

Linearizing the continuity and momentum equations gives a coupled set of
ordinary differential equations for a g-mode’s displacement:

$\displaystyle\frac{\text{d}}{\text{d}r}(r^{2}\xi_{r})$
$\displaystyle=\frac{g}{c_{\rm s}^{2}}(r^{2}\xi_{r})+\frac{r^{2}}{c_{\rm
s}^{2}}\left(\frac{L_{\ell}^{2}}{\omega^{2}}-1\right)\frac{\delta p}{\rho},$
(A2) $\displaystyle\frac{\text{d}}{\text{d}r}\left(\frac{\delta
p}{\rho}\right)$
$\displaystyle=\frac{(\omega^{2}-N^{2})}{r^{2}}(r^{2}\xi_{r})+\frac{N^{2}}{g}\frac{\delta
p}{\rho},$ (A3)

where $\rho$ is the star’s density, while $L_{\ell}^{2}$ and $N^{2}$ are the
squared Lamb and Brunt-Väisälä frequencies, respectively. Since we will be
concerned primarily with high-order g-modes, we use the Cowling approximation
(Cowling, 1941). The transverse displacement $\xi_{\perp}$ and Eularian
density perturbation $\delta\rho$ are given by

$\xi_{\perp}(r)=\frac{1}{\omega^{2}r}\frac{\delta p}{\rho},\hskip
14.22636pt\delta\rho(r)=\frac{\delta p}{c_{\rm s}^{2}}+\frac{\rho
N^{2}}{g}\xi_{r}.$ (A4)

To solve equations (A2)-(A3), we impose a few boundary conditions. Regularity
at the origin ($\text{d}[r^{2}\xi_{r}]/\text{d}r|_{r=0}$ and $\text{d}[\delta
p/\rho]/\text{d}r|_{r=0}$ finite) requires

$\lim_{r\to 0}(r^{2}\xi_{r})=\frac{\ell}{\omega^{2}}Y_{0}r^{\ell+1},\hskip
14.22636pt\lim_{r\to 0}\left(\frac{\delta p}{\rho}\right)=Y_{0}r^{\ell},$ (A5)

where $Y_{0}$ is determined by the mode normailization. A free boundary a the
star’s surface requires $\Delta p$ to vanish, or

$\left(\frac{\delta p}{\rho}-g\xi_{r}\right)_{r=R_{\star}}=0.$ (A6)

We use a shooting algorithm to calculate the displacement amplitudes
throughout the star (Press et al., 2002).

The dispersion relation of g-modes gives insight into their behavior. Assuming
$\delta p,\xi_{r}\propto{\rm e}^{{\rm i}k_{r}r}$ with $k_{r}\gg r^{-1}$, one
can show using equations (A2)-(A3) that the radial wave-number $k_{r}$ is

$k_{r}^{2}=\frac{(L_{\ell}^{2}-\omega^{2})(N^{2}-\omega^{2})}{c_{\rm
s}^{2}\omega^{2}}.$ (A7)

Low-frequency g-modes propagate where $\omega^{2}<L_{\ell}^{2},N^{2}$, and are
always evanescent in the star’s convective regions ($N^{2}\leq 0$). Since over
most of the star’s radiative region ($N^{2}>0$), g-mode frequencies satisfy
$\omega^{2}\ll L_{\ell}^{2},N^{2}$, we will frequently use the approximate
expression

$k_{r}\simeq{N\over{\omega}}k_{\perp}\simeq\frac{\sqrt{\ell(\ell+1)}N}{r\omega},$
(A8)

where $k_{\perp}=\sqrt{\ell(\ell+1)}/r$ is the azimuthal wavevector. This
dispersion relation can in turn be used to solve for the g-mode’s
eigenfrequency $\omega$ (e.g. Lai, 1994; Fuller & Lai, 2012b):

$\omega\simeq\frac{\sqrt{\ell(\ell+1)}}{\pi
n}\int_{0}^{R_{\star}}\frac{N(r)}{r}\text{d}r,$ (A9)

where $n$ is the number of radial displacement nodes (number of locations
$\xi_{r}[r]=0$), while the integral is taken over the star’s radiative zone.
We have implicitly assumed the radial node-number is large ($n\gg 1$). For the
Solar model, a 10-day period $\ell=2$ g-mode has $n\sim 1250$.

We assume a conventional normalization (e.g. Unno et al., 1979; Lai, 1994),
where the mode integral satisfies

$2\langle{\bm{\xi}}|{\bm{\xi}}\rangle=2\int_{0}^{R_{\star}}\left[\xi_{r}^{2}+\ell(\ell+1)\xi_{\perp}^{2}\right]\rho
r^{2}\text{d}r=1.$ (A10)

Here,
$\langle{\bm{u}}|{\bm{v}}\rangle\equiv\int{\bm{u}}^{*}{\bm{\cdot}}{\bm{v}}\rho\text{d}V$,
where ${\bm{u}}^{*}$ is the complex conjugate of ${\bm{u}}$. Although we will
use normailization (A10) throughout this work, because different works
frequently use different mode normailizations (e.g. Weinberg et al., 2012;
Burkart et al., 2012), it is of use to write quantities in terms of
$\epsilon=2\omega^{2}\langle{\bm{\xi}}|{\bm{\xi}}\rangle$. The quantity
$\epsilon$ is related to the mode’s energy, and with our normalization
$\epsilon=\omega^{2}$.

For our order-of-magnitude calculations to follow, we estimate amplitude for
the normalized eigen-function at the bottom of the star’s convection zone. For
a Brunt-Väisälä that rises with depth $z$ roughly as $N^{2}\sim g/z$, the
gravity-modes propagate below a depth $z_{1}$, where the upper turning point
$z_{1}\sim\omega^{2}/(gk_{h}^{2})$. For high order modes
($n\geq{\sqrt{{R_{\star}}/z_{\rm cvz}}\sim{\rm few}}$ for the sun), this
turning point lies above the bottom of the convection zone $z_{\rm cvz}$, and
$z_{1}\sim z_{\rm cvz}$. We consider this case here.

In the energy integral (eq. [A10]), every node contributes equally to the mode
energy. So integration over the first node should be $1/n$. Combine this with
$z_{1}\sim z_{\rm cvz}$, the horizontal displacement (the dominant component)
has a normalized magnitude of

$\left|{\tilde{\xi}_{\perp}}\right|^{2}_{z=z_{1}}\sim{1\over{n\,\ell^{2}\,M_{\rm
cvz}}}\,$ (A11)

where $M_{\rm cvz}$ is the mass contained within the convection zone. The odd
unit for the displacement is a result of our normalization convention. For
modes with $z_{1}\gtrsim z_{\rm cvz}$ (including cases where the star has no
surface convection zone), one should substitute $M_{\rm cvz}$ by $M_{z_{1}}$,
the mass above the first turning point. We do not consider modes in stars
without surface convection zones.

Inside the propagating cavity, according to the WKB analysis, the total energy
flux carried by the mode is constant across all radii,

$r^{2}\,\rho\xi_{\perp}^{2}\,v_{\rm group,r}\simeq{\rm constant}\,,$ (A12)

where the radial group speed $v_{\rm group,r}={{\partial\omega}\over{\partial
k_{r}}}\simeq\omega/k_{r}$. Using expression (A8), we obtain the following
radial dependence for the envelope of the pulsation amplitude (the WKB
envelope),

$\left|{\bar{\xi}}_{\perp}(r)\right|^{2}\sim{{R_{\star}^{2}\rho(z_{1})\left|{\bar{\xi}_{\perp}}(z=z_{1})\right|^{2}\,z_{1}}\over{r^{2}\rho(r)/k_{r}(r)}}\sim{{k_{r}(r)}\over{n\ell^{2}\,r^{2}\rho(r)}}\simeq{1\over{\omega
n\ell}}\,{N(r)\over{\rho r^{3}}}\,.$ (A13)

This is also consistent with the statement that every radial node contributes
the same amount of energy. Moreover, the above expression gives
$(k_{h}\xi_{\perp})^{2}\propto N/(\rho r^{5})$. This exhibits a steep rise
toward the stellar centre.

##  Appendix B Derivation of mode amplitude evolution and damping rate

In this appendix, we review the derivation of equation (D2) derived by Schenk
et al. (2002), and perturbative calculation for the damping rate (4). We do so
to ensure the perturbative expression (4) holds when the Corioulus force is
included. We begin with the tidal oscillations in the star’s rotating frame:

$\ddot{\bm{\xi}}+\mathcal{C}{\bm{\cdot}}\dot{\bm{\xi}}+\mathcal{R}{\bm{\cdot}}{\bm{\xi}}={\bm{f}}_{\rm
ext},$ (B1)

where ${\bm{f}}_{\rm ext}({\bm{r}},t)$ is an arbitrary external force, while

$\mathcal{C}{\bm{\cdot}}\dot{\bm{\xi}}=-2{\bm{\Omega}}_{\star}{\bm{\times}}\dot{\bm{\xi}}$
(B2)

is an operator for the Corioulus force, with $\mathcal{R}$ defined as before.
When $\Omega_{\star}\neq 0$, eigenfunctions ${\bm{\xi}}_{\alpha}({\bm{r}}){\rm
e}^{-{\rm i}\lambda_{\alpha}t}$ are not orthogonal, or
$\langle{\bm{\xi}}_{\alpha}|{\bm{\xi}}_{\beta}\rangle\neq 0$ when
$\alpha\neq\beta$. They do however satisfy

$\langle{\bm{\xi}}_{\alpha}|{\rm
i}\mathcal{C}{\bm{\cdot}}{\bm{\xi}}\rangle+(\lambda_{\alpha}+\lambda_{\beta})\langle{\bm{\xi}}_{\alpha}|{\bm{\xi}}_{\beta}\rangle=0$
(B3)

when $\alpha\neq\beta$. Because of this, it is non-trivial to determine
${\bm{\xi}}({\bm{r}},t)$ in terms of the star’s natural oscillation modes
${\bm{\xi}}_{\alpha}{\rm e}^{-{\rm i}\lambda_{\alpha}t}$. To find a
decomposition, we re-write equation (B1) in a higher-dimensional vector space.
Defining a new vector

${\bm{\zeta}}({\bm{r}},t)=\left[\begin{array}[]{c}{\bm{\xi}}\\\
\dot{\bm{\xi}}+\frac{1}{2}\mathcal{C}{\bm{\cdot}}{\bm{\xi}}\end{array}\right],$
(B4)

equation (B1) may be re-written as

$\dot{\bm{\zeta}}+\mathcal{M}{\bm{\cdot}}{\bm{\zeta}}={\bm{F}}({\bm{r}},t),$
(B5)

where

$\mathcal{M}=\left[\begin{array}[]{cc}-\frac{1}{2}\mathcal{C}&1\\\
-\mathcal{R}+\frac{1}{4}\mathcal{C}^{2}&-\frac{1}{2}\mathcal{C}\end{array}\right],\hskip
14.22636pt{\bm{F}}({\bm{r}},t)=\left[\begin{array}[]{c}0\\\ {\bm{f}}_{\rm
ext}({\bm{r}},t)\end{array}\right].$ (B6)

We then decompose ${\bm{\zeta}}({\bm{r}},t)$ as

$\zeta({\bm{r}},t)=\sum_{A}c_{A}(t){\bm{\zeta}}_{A}({\bm{r}}){\rm e}^{-{\rm
i}\lambda_{A}t},$ (B7)

where ${\bm{\zeta}}_{A}({\bm{r}})$ satisfies

$(\mathcal{M}+{\rm i}\lambda_{A}){\bm{\cdot}}{\bm{\zeta}}_{A}=0.$ (B8)

As in Schenk et al. (2002), we use the subscript $A$ instead of $\alpha$,
because the basis for ${\bm{\zeta}}$ has a higher dimension than ${\bm{\xi}}$.
Each left eigenvector ${\bm{\zeta}}_{A}$ has a right eigenvector
${\bm{\chi}}_{A}$ which satisfies

$(\mathcal{M}^{\dagger}-{\rm i}\lambda_{A}^{*}){\bm{\cdot}}{\bm{\chi}}_{A}=0,$
(B9)

where $\mathcal{M}^{\dagger}$ is the Hermitian conjugate of $\mathcal{M}$. It
may be shown (Schenk et al., 2002) ${\bm{\chi}}_{A}$ can be re-written as

${\bm{\chi}}_{A}=\left[\begin{array}[]{c}{\rm
i}\lambda_{A}^{*}{\bm{\tau}}_{A}-\frac{1}{2}\mathcal{C}{\bm{\cdot}}{\bm{\tau}}_{A}\\\
{\bm{\tau}}_{A}\end{array}\right],$ (B10)

where ${\bm{\tau}}_{A}$ may be written in terms of ${\bm{\xi}}_{A}$ using

${\bm{\tau}}_{A}=-\frac{1}{b_{A}^{*}}{\bm{\xi}}_{A},$ (B11)

where the constant $b_{A}$ has yet to be determined. Notice because

$\langle{\bm{\chi}}_{A}|{\bm{\zeta}}_{B}\rangle=\frac{1}{b_{A}}\big{[}\langle{\bm{\xi}}_{A}|{\rm
i}\mathcal{C}{\bm{\cdot}}{\bm{\xi}}_{B}\rangle-{\rm
i}(\lambda_{A}+\lambda_{B})\langle{\bm{\xi}}_{A}|{\bm{\xi}}_{B}\rangle\big{]},$
(B12)

${\bm{\chi}}_{A}$ is orthogonal to ${\bm{\zeta}}_{B}$. We choose $b_{A}$ so
${\bm{\chi}}_{A}$ is orthonormal to ${\bm{\zeta}}_{B}$
($\langle{\bm{\chi}}_{A}|{\bm{\zeta}}_{B}\rangle=\delta_{AB}$):

$b_{A}=\langle{\bm{\xi}}_{A}|{\rm
i}\mathcal{C}{\bm{\cdot}}{\bm{\xi}}_{A}\rangle+2\lambda_{A}\langle{\bm{\xi}}_{A}|{\bm{\xi}}_{A}\rangle.$
(B13)

Inserting equation (B7) into (B5), and using the orthonormality of
${\bm{\zeta}}_{A}$ with respect to ${\bm{\chi}}_{A}$, we obtain equation (7).

Our calculations include weak dissipative forces $\mathcal{R}^{(1)}$ in
$\mathcal{R}$, which we assume are much smaller in magnitude than
$\mathcal{R}^{(0)}=\mathcal{R}-\mathcal{R}^{(1)}$. To calculate the effects of
the dissipative forces, we decompose
$\mathcal{M}=\mathcal{M}^{(0)}+\mathcal{M}^{(1)}$, where

$\mathcal{M}^{(0)}=\left[\begin{array}[]{cc}-\frac{1}{2}\mathcal{C}&1\\\
-\mathcal{R}^{(0)}+\frac{1}{4}\mathcal{C}^{2}&-\frac{1}{2}\mathcal{C}\end{array}\right],\hskip
14.22636pt\mathcal{M}^{(1)}=\left[\begin{array}[]{cc}0&0\\\
-\mathcal{R}^{(1)}&0\end{array}\right].$ (B14)

Expanding
${\bm{\zeta}}_{A}={\bm{\zeta}}_{A}^{(0)}+{\bm{\zeta}}_{A}^{(1)}+\dots$,
$\lambda_{A}=\lambda_{A}^{(0)}+\lambda_{A}^{(1)}+\dots$, and matching terms of
like order, we find at order (1)

$(\mathcal{M}^{(0)}+{\rm
i}\lambda_{A}^{(0)}){\bm{\cdot}}{\bm{\zeta}}_{A}^{(1)}=-(\mathcal{M}^{(1)}+{\rm
i}\lambda_{A}^{(1)}){\bm{\cdot}}{\bm{\zeta}}_{A}^{(0)}.$ (B15)

Because both ${\bm{\zeta}}_{A}^{(1)}$ and $\lambda_{A}^{(1)}$ are undetermined
in equation (B15), we may choose a solution ${\bm{\zeta}}_{A}^{(1)}$ which
satisfies

$\langle{\bm{\chi}}_{A}^{(0)}|\mathcal{M}^{(0)}{\bm{\cdot}}{\bm{\zeta}}^{(1)}\rangle+{\rm
i}\lambda^{(0)}\langle{\bm{\chi}}_{A}^{(0)}|{\bm{\zeta}}^{(1)}\rangle=0.$
(B16)

Hence $\lambda_{A}^{(1)}$ is given by

$\lambda_{A}^{(1)}=-{\rm
i}\langle{\bm{\chi}}_{A}^{(0)}|\mathcal{M}^{(1)}{\bm{\cdot}}{\bm{\zeta}}_{A}^{(0)}\rangle=\frac{1}{b_{A}}\langle{\bm{\xi}}_{A}^{(0)}|\mathcal{R}^{(1)}{\bm{\cdot}}{\bm{\xi}}_{A}^{(0)}\rangle.$
(B17)

Dropping the zero superscript for ${\bm{\xi}}_{A}$, letting
$\omega_{\alpha}=\lambda_{A}^{(0)}$, ${\rm
i}\gamma_{\alpha}=\lambda_{A}^{(1)}$, and ${\bm{f}}_{\rm
diss}=\mathcal{R}^{(1)}{\bm{\cdot}}{\bm{\xi}}_{A}^{(0)}$, we obtain equation
(4) when $\mathcal{C}=0$.

##  Appendix C Damping rates for gravity Modes

In this appendix, we estimate the damping rates $\gamma$ (eq. [4]) for g-mode
oscillations due to three different drag forces ${{\bm{f}}}$: radiative
diffusion in the star’s radiative regions, turbulent diffusivity in the star’s
convective region, and wave leakage from the upper atmosphere. The damping
rates we are after are the time-averaged rates of energy dissipation

$\gamma\equiv-{{\langle{\dot{E}}\rangle}\over E}=-{{\omega}\over{2\pi}}{1\over
E}\oint\text{d}t{{\text{d}E}\over{\text{d}t}}\,,$ (C1)

where our convention is that positive rate for damping. In the following, we
present both formal integrals for these rates, as well as their order-of-
magnitude estimates. This latter exercise is necessary as it acts both as a
sanity check for our numerical results (presented in Fig. 3), and it guides
our understanding on how these rates depends on global properties (e.g.,
stellar luminosity, mass, mode frequency, etc). Moreover, we take the solar
model as a reference when estimating the rates: this assumes a surface
convection zone of depth $z_{\rm cvz}$ and mass $M_{\rm cvz}$. We also discuss
cases where a surface convection zone is absent.

###  C.1 Radiative Diffusion Damping

The radiative diffusion arises when heat is lost from the pulsational motion
via photo diffusion. This could occur anywhere inside the star, even in the
convective region, as long as radiative diffusion is carrying some flux. In
the following, we focus on the radiative zone and follow Goldreich & Wu (1999)
for much of the discussions. We will frequently use the fact that radiative
diffusion damps g-modes primarily in the radiative zone where $k_{r}r\gg 1$.

Radiative diffusion causes departure from adiabaticity and removes mode energy
at a rate of $T\text{d}s/\text{d}t$, where $s$ is entropy and is related to
the flux variation by the equation of energy conservation

$T{{\text{d}\Delta
s}\over{\text{d}t}}=-{1\over\rho}{\text{d}\over{\text{d}r}}\Delta F_{\rm
rad}\simeq-{{F_{\rm rad}}\over{\rho}}{\text{d}\over{\text{d}r}}\left({{\Delta
F_{\rm rad}}\over{F_{\rm rad}}}\right)\,,$ (C2)

where the flux variation can be related to other pulsational quantities by the
equation of radiative diffusion and is (see Goldreich & Wu, 1999)

${{\Delta F_{\rm rad}}\over{F_{\rm
rad}}}=-(1+\kappa_{\rho}){{\Delta\rho}\over\rho}+(4-\kappa_{T}){{\Delta
T}\over T}-{{\text{d}\xi_{r}}\over{\text{d}r}}+\left({{\text{d}\ln
T}\over{\text{d}r}}\right)^{-1}{\text{d}\over{\text{d}r}}\left({{\Delta
T}\over{T}}\right)\simeq\left(\frac{\text{d}\ln
T}{\text{d}r}\right)^{-1}\frac{\text{d}}{\text{d}r}\left(\frac{\Delta
T}{T}\right)-\frac{\text{d}\xi_{r}}{\text{d}r}\,.$ (C3)

Here $\kappa_{\rho}=\text{d}\ln\kappa/\text{d}\ln\rho$ and
$\kappa_{T}=\text{d}\ln\kappa/\text{d}\ln T$ are the opacity dependence on
density and temperature, neglected from here on out. Further using
$\Delta\rho/\rho=-{\bm{\nabla}}{\bm{\cdot}}{\bm{\xi}}\simeq-\xi_{r}/H$
everywhere for g-modes (where $H=c_{\rm s}^{2}/g$ is the scaleheight, valid
when $\omega^{2}\ll L_{\ell}^{2}$),

$\frac{\text{d}\ln T}{\text{d}r}=-\frac{\text{d}\ln
T}{\text{d}\ln\rho}\frac{g}{c_{\rm e}^{2}},$ (C4)

where $c_{\rm e}^{2}=\text{d}p/\text{d}\rho$ is the equilibrium sound speed,
and the quasi-adiabatic approximation
($\Delta\rho/\rho=[\text{d}\ln\rho/\text{d}\ln T]_{S}[\Delta T/T]$), the flux
perturbation becomes

$\frac{\Delta F_{\rm rad}}{F_{\rm rad}}\simeq{\rm
i}k_{r}\left[\left(\frac{\text{d}\ln\rho}{\text{d}\ln
T}\right)_{S}\frac{c_{\rm s}^{2}}{g}-\left(\frac{\text{d}\ln\rho}{\text{d}\ln
T}\right)\frac{c_{\rm e}^{2}}{g}\right]\frac{\Delta T}{T}.$ (C5)

In addition, we make the crude approximation

$\left(\frac{\text{d}\ln\rho}{\text{d}\ln T}\right)_{S}\frac{c_{\rm
s}^{2}}{g}-\left(\frac{\text{d}\ln\rho}{\text{d}\ln T}\right)\frac{c_{\rm
e}^{2}}{g}\approx\left(\frac{\text{d}\ln\rho}{\text{d}\ln
T}\right)_{S}\frac{c_{\rm s}^{2}}{g}.$ (C6)

Although this approximation will overestimate the magnitude of $\gamma_{\rm
diff}$, it simplifies its numerical calculation, and still gives damping rates
within a factor of $\sim 2-3$ of more detailed calculations (e.g. Fuller,
2017, Fig. 2). Inserting equation (C2) into (C1), we obtain

$\gamma_{\rm diff}=\frac{\omega}{2\pi\epsilon}\oint\text{d}t\int_{\rm
rad}\left(\frac{\Delta T}{T}\right)\frac{\text{d}\Delta F_{\rm
rad}}{\text{d}r}r^{2}\text{d}r\simeq\frac{1}{\epsilon}\int_{\rm rad}F_{\rm
rad}(\Gamma_{3}-1)Hk_{r}^{2}\left|\frac{\Delta\rho}{\rho}\right|^{2}r^{2}\text{d}r,$
(C7)

where the adiabatic exponent $\Gamma_{3}-1=(\text{d}\ln
T/\text{d}\ln\rho)_{S}$. Equation (C7) is the expression we use to numerically
calculate $\gamma_{\rm diff}$ for our stellar models.

To estimate the value of $\gamma_{\rm diff}$ analytically, we use
$|\Delta\rho/\rho|\simeq|\xi_{r}/H|\simeq|(k_{\perp}/k_{r})(\xi_{\perp}/H)|$
and $r^{2}F_{\rm rad}\sim L_{\star}$ in the outer regions of the star’s
radiative core ($r\gtrsim 0.2\,{\rm R}_{\odot}$ for PMS and MS models). Then
$\gamma_{\rm diff}$ becomes

$\gamma_{\rm diff}\sim\frac{L_{\star}}{\epsilon}\int_{\rm
rad}\frac{k_{\perp}^{2}\xi_{\perp}^{2}}{H}\text{d}r.$ (C8)

Because $r^{2}F_{\rm rad}\approx{\rm constant}$ only applies to the outer
regions of the radiative core, and $r^{2}F_{\rm rad}\to 0$ as the core’s
center is approached, the integrand remains finite even as
$k_{\perp}^{2}\xi_{\perp}^{2}/H$ diverges as $r\to 0$. To evaluate this
integral, we use our estimate for the WKB envelope of $|\xi_{\perp}|^{2}$ in
the star’s radiative zone (eq. [A13]),
$k_{r}\simeq(N/\omega)k_{\perp}\sim(N/\omega)(\ell/r)$, and the fact every
radial node of $\xi_{\perp}$ ha s a width $\Delta r\sim
k_{r}^{-1}\sim\text{d}r$ to obtain

$\gamma_{\rm
diff}\sim\frac{L_{\star}}{\epsilon}\left[n\left(\frac{r\omega}{\ell
N}\right)\frac{1}{H}\left(\frac{\ell^{2}}{r^{2}}\right)\left(\frac{1}{\omega
n\ell}\frac{N}{\rho r^{3}}\right)\right]_{\rm
rad}\sim\frac{L_{\star}}{M_{\star}R_{\star}^{2}\omega^{2}}\Lambda_{\rm core},$
(C9)

where the quantity

$\Lambda_{\rm core}=\left(\frac{M_{\star}R_{\star}^{2}}{\rho
Hr^{4}}\right)_{r=r_{\rm c}}$ (C10)

measures the mass concentration in the star’s radiative core at the
characteristic radius $r_{\rm c}$, and has values of order $\Lambda_{\rm
core}\sim 10^{2}-10^{3}$ ($10^{3}-10^{4}$) during the PMS (MS), and increases
as the star evolves and becomes more centrally concentrated. The value of
$\Lambda_{\rm core}$ diverges near the center, but does not overwhelm the
contribution to the integral for $\gamma_{\rm diff}$ is because $r^{2}F_{\rm
rad}$ vanishes there.

Evaluating this for a 10 day period g-mode in the sun, we obtain $\gamma_{\rm
rad}\sim 2\times 10^{-3}\ {\rm day}^{-1}\sim 2\times 10^{-8}\ {\rm s}^{-1}$.
This agrees with the numerical calculation in Table 1 of Terquem et al. (1998)
and the illustrative estimate in equation (21) of Goodman & Dickson (1998).

We ignore radiative diffusion from inside the convection zone. This is
justified because radiation carries only a tiny fraction of the total flux
there. Our estimated damping rate (eq. [C9]) does not apply to stars without a
surface convection zone, since most of the diffusion occurs in the deep
interior.

###  C.2 Turbulent Damping

Our treatment of the force acting on g-modes from turbulent eddies in the
star’s convective zone is more heuristic. We approximate the viscous force
(which acts on the gradient of perturbed velocity) as

${\bm{f}}_{\rm turb}=-\frac{{\rm i}\omega}{r^{2}\rho}\frac{\partial}{\partial
r}\left(r^{2}\rho\nu_{\rm turb}\frac{\partial{\bm{\xi}}}{\partial
r}\right)\,,$ (C11)

where $\nu_{\rm turb}$ is the effective turbulent viscosity. To account for
the two extremes, when the turn-over time is much shorter or much longer than
the mode period, we adopt the following prescription for $\nu_{\rm turb}$
(Goodman & Oh, 1997),

$\nu_{\rm turb}=\frac{H_{\rm p}v_{\rm cv}}{1+(\omega t_{\rm cv}/2\pi)^{s}}.$
(C12)

Here $H_{\rm p}=(\text{d}\ln p/\text{d}r)^{-1}\sim z$ is the local pressure
scale-height, and $z$ the depth from the surface. $v_{\rm cv}$ is the
convective velocity, and $t_{\rm cv}=H_{\rm p}/v_{\rm cv}$ is the eddy
turnover time. In the limit of slow convection, Goldreich & Nicholson (1977)
have argued that a Kolmogorov turbulent spectrum gives rise to $s=2$. In
particular, for eddies of a range of sizes, the most relevant ones for
turbulent viscosity are those that have turn-over time comparable to the tidal
period, $t_{\lambda}=\lambda/v_{\lambda}$ ($t_{\lambda}\omega\sim 2\pi$),
where $\lambda$ and $v_{\lambda}$ are their length-scales and velocities,
respectively. This, also advocated later by Goodman & Oh (1997), differs from
the prescription by Zahn (1966), where he set $s=1$. Direct numerical
simulations have yet to conclusively settle this issue (e.g. Penev et al.,
2007; Ogilvie & Lesur, 2012; Duguid et al., 2020; Vidal & Barker, 2020).

We insert this into equation (4), integrate by parts, and obtain an expression
for the turbulent damping rate,

$\gamma_{{\rm turb}}=-\frac{\omega^{2}}{\epsilon}\int_{\rm cv}\nu_{\rm
turb}\left|\frac{\partial{\bm{\xi}}}{\partial r}\right|^{2}\rho
r^{2}\text{d}r\,,$ (C13)

where $\epsilon$ is the mode energy,
$\epsilon=2\omega^{2}\langle{\bm{\xi}}|{\bm{\xi}}\rangle=\omega^{2}$, for our
adopted normalization.

To evaluate this integral at the order-of-magnitude level, we first estimate
the magnitude of $\xi_{r}$ and $\xi_{\perp}$ in the convective zone. Since
$N^{2}\approx 0$, $H\ll r$, and $\omega^{2}\ll L_{\ell}^{2}$, the ODEs for
$\xi_{r}$ and $\xi_{\perp}$ are

$\frac{\text{d}\xi_{r}}{\text{d}r}\simeq\frac{\xi_{r}}{H}+k_{\perp}^{2}r\xi_{\perp},\hskip
14.22636pt\frac{\text{d}\xi_{\perp}}{\text{d}r}\simeq\frac{\xi_{r}}{r}-\frac{\xi_{\perp}}{r}.$
(C14)

Because $H\sim z=R_{\star}-r\ll r$ in the convective zone, equation (C14) has
the approximate solution

$\xi_{r}\approx-\frac{1}{2}k_{\perp}^{2}rz\xi_{\perp}+\frac{A}{z^{2}},$ (C15)

where the constant $A$ is determined by regularity conditions. Regularity as
$z\to 0$ requires $A=0$, so an order-of-magnitude estimate for $\xi_{r}$ in
the convective zone is

$\xi_{r}\sim-k_{\perp}^{2}rH\xi_{\perp}.$ (C16)

Notice equations (C14) and (C16) imply $|\xi_{r}|\lesssim|\xi_{\perp}|$ but
$|\text{d}\xi_{r}/\text{d}r|\gtrsim|\text{d}\xi_{\perp}/\text{d}r|$. Hence in
the convective zone, $|\partial{\bm{\xi}}/\partial r|\sim
k_{\perp}^{2}r\xi_{\perp}$.

In the deeper part of the convection zone, $\omega t_{\rm cv}\gtrsim 2\pi$,
and

$\rho\nu_{\rm turb}\sim{{\rho v_{\rm cv}^{3}}\over{\omega^{2}z}}\sim{{F_{\rm
cv}}\over{\omega^{2}z}}\,,$ (C17)

where $F_{\rm cv}$ is the convective flux and is of order the total flux,
$F_{\rm cv}=L_{\star}/4\pi R^{2}$. In contrast, near the surface, $\omega
t_{\rm cv}\ll 1$, and

$\rho\nu_{\rm turb}\sim\rho v_{\rm cv}z\sim\rho^{2/3}zF_{\rm cv}^{1/3}\,.$
(C18)

As a result, the dominant contribution to viscous damping arises where $\omega
t_{\rm cv}\sim 1$, (see also Wu & Lithwick, 2019), at the depth

$z_{\nu}\sim\frac{v_{\rm cv}}{\omega}\sim\frac{1}{\omega}\left(\frac{F_{\rm
cv}}{\rho_{\rm cv}}\right)^{1/3},$ (C19)

where $\rho_{\rm cv}\sim M_{\rm cvz}/(R_{\star}^{2}z_{\rm cvz})$. Hence, the
turbulent damping rate is of order

$\displaystyle\gamma_{\rm turb}$
$\displaystyle\sim-\frac{\omega^{2}}{\epsilon}\int_{\rm cvz}\left(\frac{F_{\rm
cv}}{\omega^{2}z_{\nu}}\right)(k_{\perp}^{2}r\xi_{\perp})^{2}r^{2}\text{d}r\sim\bar{\gamma}_{\rm
turb}\frac{\ell^{2}\omega}{n}\left(\frac{z_{\rm
cvz}^{2/3}L_{\star}^{2/3}}{\omega^{2}R_{\star}^{2}M_{\rm cvz}^{2/3}}\right)$
$\displaystyle\sim\bar{\gamma}_{\rm turb}\
\frac{\omega}{n}\times\frac{1}{(\tau_{\rm th}\omega_{\rm
dyn})^{2/3}}\times\left(\frac{z_{\rm
cvz}}{z_{1}}\right)\times\left(\frac{z_{\rm cvz}}{R_{\star}}\right)^{1/3},$
(C20)

where $\tau_{\rm th}\sim M_{\rm cvz}z_{\rm cvz}g/L_{\star}$ is the thermal
time at the base of the convection zone, and $\bar{\gamma}_{\rm turb}$ is an
order-unity coefficient used to reduce the magnitude of $\gamma_{\rm turb}$
(since we evaluated $\rho\nu_{\rm turb}$ at it’s maximum value). We find the
reduction in the turbulent damping estimate (C20) of $\bar{\gamma}_{\rm
turb}\sim 1/3$ better matches the more detailed calculations using equation
(C13).

For the Solar-model, $\gamma_{\rm turb}\sim 6\times 10^{-4}(\omega/n)\sim
4\times 10^{-7}\,{\rm day}^{-1}\sim 5\times 10^{-12}\,{\rm s}^{-1}$ for a
$\ell=2$ mode at $10$ day period, assuming $n\sim 1000$. Pre-main-sequence
stars, with their deeper convection zones (and slower convection), experience
stronger damping. In contrast, stars with shallow convection zones can’t
remove mode energy as efficiently. Lastly, given the dispersion relation for
g-modes ($\omega\propto l/n$), the damping rate is independent of mode radial
order, $\gamma_{\rm turb}\propto\ell$.

Luan et al. (2018) got a different scaling with mode number ($\gamma_{\rm
diff}\propto n$, see their eq. [A13]), because the turbulent dissipation rate
$\rho\nu_{\rm turb}$ was evaluated at the base of the convection zone ($z_{\rm
cvz}$), rather than where turbulence peaked within the convection zone
($z_{\nu}$). The numerical calculations in this work support our estimate for
$\gamma_{\rm turb}$.

###  C.3 Wave Leakage

Even for stars with a surface convection zone, the atmosphere above the
photosphere is stratified. This creates a g-mode cavity through which the
internal g-modes that are trapped below the surface convection zone can tunnel
through and escape. We proceed to estimate the associated wave damping rate.
Some of these calculations were also detailed in Wu & Lithwick (2019) in the
context of f-modes.

In this region, the stratification $N\sim\sqrt{g/H}$, where $H$ is the
photospheric scale height, and the sound speed $c_{s}\sim\sqrt{gH}\sim NH$.
The dispersion relation is

$k_{z}={1\over 2}{{\Gamma_{1}}\over H}\pm{1\over
2}\sqrt{\left({{\Gamma_{1}}\over{H}}\right)^{2}-4\left[k_{\perp}^{2}\left({{N^{2}}\over{\omega^{2}}}-1\right)+{{\omega^{2}}\over{c_{s}^{2}}}\right]}\,.$
(C21)

So propagative waves have either frequencies higher than the acoustic cut-off
frequency $\omega_{\rm ac}={\Gamma_{1}c_{s}}/(2H)$, or frequencies lower than

$\omega_{\rm ac,2}=\omega_{\rm ac}{H\over
R_{\star}}=\Gamma_{1}{{c_{s}}\over{2R_{\star}}}\,,$ (C22)

Low frequency gravity-waves with $\omega<\omega_{\rm ac,2}\sim 2\pi/(5\,{\rm
days})$ (evaluated with solar values, $H=20$km and $T=5800$K) can leak out,
with a simplified dispersion relation of the form

$k_{z}=k_{\perp}{N\over\omega}\approx{{\ell}\over{H}}\left({{\omega_{\rm
ac,2}\over{\omega}}}\right)\,.$ (C23)

The rate of energy leakage is determined by the flux of kinetic energy
(dominated by the horizontal velocity) at the photosphere

${\dot{E}}_{\rm
leak}=R_{\star}^{2}\times\left[\rho\omega^{2}\ell(\ell+1)\left|{\xi_{\perp}}\right|^{2}\times
v_{\rm group,z}\right]_{\rm ph}\,,$ (C24)

where $V_{\rm group,z}\approx\omega/k_{z}$ is the wave’s vertical group
velocity. This leads us to an estimate for the damping rate

$\gamma_{\rm leak}\equiv{{{\dot{E}}_{\rm
leak}}\over{E}}=R_{\star}^{2}\rho_{\rm
ph}\left|{\tilde{\xi}}_{\perp}\right|^{2}{\omega\over{k_{z}}}\approx{\omega\over{n\ell}}\times\left({{M_{\rm
ph}}\over{M_{\rm cvz}}}\right)\times\left({\omega\over{\omega_{\rm
ac,2}}}\right)\,,$ (C25)

where the photospheric mass $M_{\rm ph}\approx R_{\star}^{2}\rho_{\rm ph}H$,
and we have used equation (C23) above, and equation (A11) for the normalized
horizontal displacement, assuming that there is little gradient in the
horizontal displacement across the convection zone.

For the Solar model, $M_{\rm ph}/M_{\rm cvz}\sim 10^{-9}$, and $\gamma_{\rm
leak}\sim 10^{-9}(\omega/n\ell)$ for a mode near $\omega_{\rm ac,2}$, and is
much below the convective damping estimated above. However, leaking becomes
increasingly more important as the convection zone thins. In the limit of a
vanishing convection zone, g-modes with $\omega<\omega_{ac,2}$ are no longer
trapped inside the star and are instead running waves with $\gamma_{\rm
leak}\sim\omega/n$.

##  Appendix D Adiabadic Evolution of The Amplitude Equation

This appendix derives the regime of validity for the mode amplitude $c(t)$ to
be given by equation (21). We start with the evolution equation for $c$
excited by the tidal component $\\{k,\ell,m\\}$:

$\dot{c}+({\rm i}\omega+\gamma)c={\rm i}\omega U_{k\ell m}{\rm e}^{-{\rm
i}\omega_{km}t}.$ (D1)

Equation (D1) has the general solution (e.g. Fuller & Lai, 2012a; Vick et al.,
2019)

$c(t)=c(0){\rm e}^{-\int_{0}^{t}({\rm i}\omega+\gamma)\text{d}t^{\prime}}+{\rm
e}^{-\int_{0}^{t}({\rm i}\omega+\gamma)\text{d}t^{\prime}}\int_{0}^{t}{\rm
i}\omega U_{k\ell m}{\rm e}^{\int_{0}^{t^{\prime}}\left[{\rm
i}(\omega-\omega_{km})+\gamma\right]\text{d}t^{\prime\prime}}\text{d}t^{\prime}$
(D2)

The first and second terms in equation (D2) are homogeneous and inhomogeneous
solutions to equation (D1), respectively. To show how equation (D2) reduces to
(21), we first notice for the systems of interest, their damping timescales
are much shorter than their ages ($\gamma^{-1}\ll t$), which implies the
homogeneous solution quickly damps by age $t$ to give

$c(t)\simeq\int_{0}^{t}{\rm i}\omega U_{k\ell m}{\rm
e}^{-\int_{t^{\prime}}^{t}({\rm i}\omega+\gamma)\text{d}t^{\prime\prime}}{\rm
e}^{-{\rm
i}\int_{0}^{t^{\prime}}\omega_{km}\text{d}t^{\prime\prime}}\text{d}t^{\prime}.$
(D3)

We further note that binaries typically have mode damping times much shorter
than mode evolution times ($\gamma^{-1}\ll t_{\rm evol}$), which implies all
expressions outside the damped exponential

$\exp\left(-\int_{t^{\prime}}^{t}\gamma\text{d}t^{\prime\prime}\right)\simeq\exp\left[-\gamma(t)(t-t^{\prime})\right]$
(D4)

contribute negligibly to the integral from $0$ to $t$ in (D3) when
$t^{\prime}\ll t$. Also assuming $\omega,\omega_{km}\gg t_{\rm evol}^{-1}$, we
obtain equation (21). The key condition for solution (21) to hold is the mode
damping rate be sufficiently high ($\gamma\gg t^{-1},t_{\rm evol}^{-1}$).

Burkart et al. (2014) derived a different condition for adiabadicity:

$\gamma^{2}\gtrsim\omega/t_{\rm evol}.$ (D5)

We briefly repeat the derivation of equation (D5). When $|c|\ll 1$,
$\omega=\omega_{km}$, $\omega\gg t_{\rm evol}^{-1}$, and crucially
$\gamma^{-1}\gg t$ (where $t$ is the age of the star), equation (D2) reduces
to

$c(t)\simeq{\rm i}\omega tU_{k\ell m}{\rm e}^{-{\rm i}\omega t}.$ (D6)

The time derivative of equation (D6) is assumed to be the maximum evolution
rate a mode near resonance could evolve as:

$\left.\frac{\text{d}|c|}{\text{d}t}\right|_{\rm max}=\omega U_{k\ell m}.$
(D7)

When equation (D7) is less than

$\left.\frac{\text{d}|c|}{\text{d}t}\right|_{\rm
ev}\sim\frac{\omega^{2}}{(\omega-\omega_{km})^{2}-\gamma^{2}}\frac{U_{k\ell
m}}{t_{\rm evol}},$ (D8)

which is roughly the time derivative of the in-homogeneous mode amplitude
(21), the locked mode can no longer keep up with the systems evolution, and
equation (21) no longer holds. However, a crucial step in deriving solution
(D6) for $c(t)$ at resonance is the assumption $\gamma^{-1}\gg t$. Since most
astrophysical systems satisfy $\gamma^{-1}\ll t$, equation (D5) is usually an
invalid condition for adiabadicity.

Physically, the value of $\text{d}|c|/\text{d}t|_{\rm max}$ given by equation
(D7) is roughly the rate a mode out of equilibrium “catches up” to its
equilibrium amplitude (21) when mode damping $\gamma$ is negligible. However,
if the mode damping rate is sufficiently high, any memory of a previous mode
amplitude is erased, and remains approximately equal to its adiabatic
amplitude (21).

##  Appendix E resonance locking from sectoral Modes

This section looks at resonance locking from stellar evolution due to sectoral
modes ($m\neq 0$), which exert a torque on the host star driving the evolution
of the stellar spin frequency $\Omega_{\star}$. As was stated in Section 3.2,
these $m\neq 0$ modes are either very difficult for binaries to lock onto, or
will drive negligible inward migration/circularization compared to zonal
modes. This Appendix discusses why these modes are unstable to resonance
locking.

###  E.1 Stable resonance locks with evolving stellar spin

When a non-zero torque is exerted on the host star, both $\Omega$ and
$\Omega_{\star}$ evolve as the binary tidally migrates inward. The tidal
torque may be written as (Burkart et al., 2014; Fuller, 2017)

$\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
tide}=-\frac{m}{k\Omega I_{\star}}\dot{E}_{\rm
orb}=\frac{m\mu_{\star}^{\prime}}{2k\kappa_{\star}}\left(\frac{a}{R_{\star}}\right)^{2}\Omega\frac{\dot{E}_{\rm
orb}}{E_{\rm orb}},$ (E1)

where $\kappa_{\star}=I_{\star}/(M_{\star}R_{\star}^{2})$ and
$\mu_{\star}^{\prime}=M_{\star}^{\prime}/(M_{\star}+M_{\star}^{\prime})=(1+q)^{-1}$.
Assuming tidal torques dominate the spin evolution
($\Omega_{\star}/t_{\star}\simeq\text{d}\Omega_{\star}/\text{d}t|_{\rm tide}$,
see eqs. [E9]-[E11] below), equation (39) may be solved for the binary orbital
evolution during resonance locks:

$\frac{\dot{E}_{\rm orb}}{E_{\rm orb}}=\frac{2}{3t_{\rm
mode}}\left(1-\frac{m\Omega_{\star}}{k\Omega}\right)\Bigg{/}\left[1-\frac{m^{2}\mu_{\star}^{\prime}B}{3k^{2}\kappa_{\star}}\left(\frac{a}{R_{\star}}\right)^{2}\right].$
(E2)

Because the orbital evolution at resonance is (using eq. [29] with $\sum
c\simeq C^{\rm max}$)

$\frac{\dot{E}_{{\rm orb};k\ell m}^{\rm max}}{E_{\rm
orb}}=\frac{4k^{2}\epsilon{\rm Re}(U_{k\ell m})}{\mu_{\star}^{\prime}\gamma
a^{2}M_{\star}}\left(1-\frac{m\Omega_{\star}}{k\Omega}\right),$ (E3)

resonance locking is unstable when the denominator of equation (E2) is
negative, since the resonant mode locations evolve $E_{\rm orb}$ in the
opposite direction as tidal angular momentum and energy dissipation (eq. [E3];
see Fig. 5, case $t_{\rm evol}$ [c]). Here ${\rm Re}(X)$ denotes the real part
of $X$.

Stable resonance locks require the denominator of equation (E2) to be
positive. Solving for $k=k_{\rm c}$ when the denominator of equation (E2) is
zero gives

$k_{\rm c}=\pm
m\sqrt{\frac{B\mu_{\star}^{\prime}}{3\kappa_{\star}}}\left(\frac{a}{R_{\star}}\right).$
(E4)

The critical orbital harmonic $k_{\rm c}$ was also derived by Fuller & Lai
(2012a); Burkart et al. (2014), where orbital harmonics $|k|<|k_{\rm c}|$
cannot stably resonance lock. Since $k_{\rm c}\gg 1$ for typical binary star
parameters, $m\neq 0$ tidal components lock less easily than $m=0$ tidal
components.

###  E.2 Stable resonance locks with near-synchronous stellar spin

Instead of assuming $\Omega_{\star}$ evolves with the locked mode as in the
previous section, this section considers a situation where $\Omega_{\star}$
evolves to an equilibrium state over a timescale much faster than the binary
orbital evolution (since $I_{\star}\ll M_{\star}a^{2}$ typically). The stellar
spin $\Omega_{\star}(t)$ for our model (eq. [E5]) is then set by the balance
of tidal torques, stellar winds and contraction/expansion:

$\frac{\text{d}\Omega_{\star}}{\text{d}t}=\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
tide}+\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
wind}+\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm inert}=0.$
(E5)

The first, second, and third terms in equation (E5) represent stellar spin
up/down from the tidal torque from mode $\alpha$ and potential component
$\\{k,\ell,m\\}$, stellar winds, and stellar contraction/expansion,
respectively. Equations (39) and (E1) may be shown to reduce to (when
$t_{\star}\to\infty$)

$\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
tide}=\frac{m\mu_{\star}^{\prime}}{3k\kappa_{\star}}\left(\frac{a}{R_{\star}}\right)^{2}\left(1-\frac{m\Omega}{k\Omega_{\star}}\right)\frac{\Omega}{t_{\rm
mode}}$ (E6)

We parameterize stellar wind spin-down with (Skumanich, 1972)

$\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
wind}=-\alpha_{\star}\Omega_{\star}^{3},$ (E7)

where $\alpha=1.5\times 10^{-14}$ for a solar-type star (e.g. Verbunt & Zwaan,
1981; Barker & Ogilvie, 2009). Spin up (down) from stellar contraction
(expansion) is

$\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
inert}=-\frac{\Omega_{\star}}{\tau_{\star}},$ (E8)

where $\tau_{\star}=I_{\star}(\text{d}I_{\star}/\text{d}t)^{-1}$ is the spin
up/down timescale from changes in $I_{\star}$. The $\dot{\Omega}_{\star}$
contributions from stellar winds and contraction/expansion may be normalized
to the tidal torque acting on the host star:

$\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
wind}=\frac{m\mu_{\star}^{\prime}}{3k\kappa_{\star}}\left(\frac{a}{R_{\star}}\right)^{2}\frac{\Omega}{t_{\rm
mode}}\Gamma_{km}^{\rm w}\bar{\Omega}_{\star}^{3},\hskip
14.22636pt\left.\frac{\text{d}\Omega_{\star}}{\text{d}t}\right|_{\rm
inert}=\frac{m\mu_{\star}^{\prime}}{3k\kappa_{\star}}\left(\frac{a}{R_{\star}}\right)^{2}\frac{\Omega}{t_{\rm
mode}}\Gamma_{km}^{\star}\bar{\Omega}_{\star},$ (E9)

where $\bar{\Omega}_{\star}=\Omega_{\star}/\Omega$,

$\displaystyle\Gamma_{km}^{\rm w}$
$\displaystyle=\frac{3k\kappa_{\star}}{m\mu_{\star}^{\prime}}\left(\frac{R_{\star}}{a}\right)^{2}t_{\rm
mode}\alpha_{\star}\Omega^{2}=0.0395\frac{k}{m}\left(\frac{\kappa_{\star}}{0.05}\right)\left(\frac{0.3}{\mu_{\star}^{\prime}}\right)\left(\frac{R_{\star}}{a}\right)^{2}\left(\frac{t_{\rm
mode}}{10^{6}\,{\rm yrs}}\right)\left(\frac{\alpha_{\star}}{1.5\times
10^{-14}\,{\rm yrs}}\right)\left(\frac{1\,{\rm day}}{P_{\rm orb}}\right)^{2},$
(E10) $\displaystyle\Gamma_{km}^{\star}$
$\displaystyle=\frac{3k\kappa_{\star}}{m\mu_{\star}^{\prime}}\left(\frac{R_{\star}}{a}\right)^{2}\frac{t_{\rm
mode}}{\tau_{\star}}=0.5\frac{k}{m}\left(\frac{\kappa_{\star}}{0.05}\right)\left(\frac{0.3}{\mu_{\star}^{\prime}}\right)\left(\frac{R_{\star}}{a}\right)^{2}\frac{t_{\rm
mode}}{\tau_{\star}}.$ (E11)

Since $a>R_{\star}$ and $\tau_{\star}\gtrsim t_{\rm evol}$, we typically have
$\Gamma_{km}^{\rm w},\Gamma_{km}^{\star}\ll 1$, so tidal torques dominate the
$\Omega_{\star}$ evolution.

Equations (E9)-(E11) may be re-arranged to give an equation for the
equilibrium $\Omega_{\star}$:

$1-\frac{m}{k}\bar{\Omega}_{\star}-\Gamma_{km}^{\rm
w}\bar{\Omega}_{\star}^{3}-\Gamma_{km}^{\star}\bar{\Omega}_{\star}=0.$ (E12)

The general analytic solution to the cubic equation (E12) is complicated, but
since $\Gamma_{km}^{\rm w},\Gamma_{km}^{\star}\ll 1$ for typical binary star
parameters, defining

$\Delta\bar{\Omega}_{km}\equiv 1-\frac{m}{k}\bar{\Omega}_{\star},$ (E13)

and assuming $\Delta\bar{\Omega}_{km}\ll 1$, we get the analytic approximate
solution

$\Delta\bar{\Omega}_{km}\simeq\frac{k}{m}\left(\Gamma_{km}^{\star}+\frac{k^{2}}{m^{2}}\Gamma_{km}^{\rm
w}\right).$ (E14)

Thus the semi-major axis evolution may be calculated using

$\Omega(a)=\frac{\omega}{k\Delta\bar{\Omega}_{km}}\simeq\frac{m\omega}{k}\left(\Gamma_{km}^{\star}+\frac{k^{2}}{m^{2}}\Gamma_{km}^{\rm
w}\right)^{-1}.$ (E15)

When $m\neq 0$, resonance locking is only possible for high-order modes
($\omega\ll k\Omega/m$), which will have higher damping rates and lower
overlap integrals than $m=0$ modes (since $\omega\simeq k\Omega$; see Fig. 3).
Therefore, $m\neq 0$ resonance locks with nearly synchronous spins are
significantly less stable than $m=0$ resonance locks, contributing far less to
the circularization of binary star populations.

##  Appendix F Non-linearity and dissipation of resonantly excited gravity-
modes

###  F.1 Estimate of Core Non-linearity

This section estimates the degree of non-linearity for resonantly driven
g-modes, as they propagate into the stably stratified core of a young star.
For definiteness, we focus on a 10 Myr old model of the Sun, with
$R_{\star}\approx 1.2R_{\odot}$, and a core density $\rho_{\rm c}\approx
10\,{\rm g}/{\rm cm}^{3}$. Estimates using a 3 Myr model lead to similar
conclusions.

Near the centre, the stratification runs as (Fig. 1)

$N(r)\approx
2\left(\frac{r}{R_{\star}}\right)\left(\frac{GM_{\star}}{R_{\star}^{3}}\right)^{1/2}\,,$
(F1)

so the lower resonant cavity of a g-mode, where $\omega=N(r)$, is

$r_{\rm m}\approx
0.5\left(\frac{\omega^{2}R_{\star}^{3}}{GM_{\star}}\right)^{1/2}\approx
0.075\,k\left(\frac{1\,{\rm day}}{P_{\rm orb}}\right){\rm R}_{\odot}\,,$ (F2)

or, almost $100$ times further out than that in a 5 Gyrs old Sun (eq. [32] in
Goodman & Dickson 1998).

Using the WKB envelope for the horizontal displacement (eq. [A13]), the
maximal non-linearity of a resonantly locked g-mode ($\omega\simeq k\Omega$,
and $\Delta=|\omega-\omega_{km}|/\omega\leq\gamma/\omega$) is reached near
$r_{\rm m}$,

$|k_{r}\xi_{r}|_{r=r_{\rm m}}^{\rm max}\simeq|k_{\perp}\xi_{\perp}|_{r=r_{\rm
m}}^{\rm max}=C^{\rm max}|k_{\perp}\tilde{\xi}_{\perp}|_{r=r_{\rm
m}}\sim\left(\frac{\omega}{\gamma}\right)U_{k\ell m}\frac{\ell}{r_{\rm
m}}\left(\frac{1}{n\ell}\frac{1}{\rho_{\rm c}r_{\rm m}^{3}}\right)^{1/2}.$
(F3)

The dispersion relation for a $\\{\ell,m\\}=\\{2,0\\}$ mode in such a model is
(eq. [A9]), $n\approx 10\,k^{-1}({P_{\rm orb}}/{1\,{\rm day}})$. We also
obtain the following numerical values from our calculations (see Figs. 3 \- 4,
eqs. [10] & [28])

$\displaystyle\frac{\gamma}{\omega}$ $\displaystyle\approx$
$\displaystyle\frac{\gamma_{\rm diff}}{\omega}\approx 1.2\times 10^{-8}\
k^{-3}\left(\frac{P_{\rm orb}}{1\,{\rm day}}\right)^{3},$ $\displaystyle
I_{20}$ $\displaystyle\approx$ $\displaystyle 5\times 10^{-5}\
k^{5/2}\left(\frac{1\,{\rm day}}{P_{\rm orb}}\right)^{5/2}{\rm
M}_{\odot}^{-1/2}{\rm R}_{\odot}^{-1}$ $\displaystyle\epsilon$
$\displaystyle\approx$ $\displaystyle 0.013\ k^{2}\left(\frac{1\,{\rm
day}}{P_{\rm orb}}\right)^{2}\frac{G{\rm M}_{\odot}}{{\rm R}_{\odot}^{3}},$
$\displaystyle\frac{a}{R_{\star}}$ $\displaystyle\approx$ $\displaystyle
3.9\left(\frac{P_{\rm orb}}{1\,{\rm day}}\right)^{2/3}$ $\displaystyle
U_{k\ell m}$ $\displaystyle\approx$ $\displaystyle 1.3\times 10^{-5}\
X_{k20}(e)k^{1/2}\left(\frac{1\,{\rm day}}{P_{\rm orb}}\right)^{5/2}{\rm
M}_{\odot}^{1/2}{\rm R}_{\odot}.$ (F4)

Putting everything together, we find the dimensionless non-linearity parameter
to be

$|k_{r}\xi_{r}|_{r=r_{\rm m}}^{\rm max}\sim 2.4\times 10^{5}\
X_{k20}(e)k^{3/2}\left(\frac{1\,{\rm day}}{P_{\rm orb}}\right)^{7/2}.$ (F5)

For our discussion, we set the threshold for strong non-linearity to be
$|k_{r}\xi_{r}|_{r=r_{m}}=10$, on account of the fact that WKB estimate fails
near the turning point for which the above estimate is made (Goodman &
Dickson, 1998). For $k=1$ and $e\ll 1$, $X_{120}\simeq 3e/2$, we find that
non-linearity is important if the binary period

$P_{\rm orb}\lesssim P_{\rm nl}\sim 15.5\,{\rm
days}\,\left({e\over{0.4}}\right)^{2/7}\,.$ (F6)

This is longer than the limit for resonance locking (eq. [43]), indicating
that for the Sun at 10 Myrs old, the efficacy of resonance locking with
gravity-modes is severely compromised by non-linearity.

###  F.2 Mode amplitude at nonlinear Breaking And Effective Dissipation

Here, we estimate the amplitude of a g-mode when wave-breaking becomes
important ($C=C_{\rm sat}$). Scaling from equation (F5), an estimate for
$C=C^{\rm max}$, we find a threshold amplitude
($|k_{r}\xi_{r}(r=r_{m})|\approx 10$) to be

$\frac{|C^{\rm sat}|}{|C^{\rm max}|}\sim 4.2\times
10^{-5}\,X_{k20}^{-1}(e)k^{-3/2}\left(\frac{P_{\rm orb}}{1\,{\rm
day}}\right)^{7/2}\,.$ (F7)

We speculate that the following dynamical picture determines the effective
dissipation. Without nonlinear breaking, the g-mode is pumped by tidal forcing
to reach its maximum amplitude ($C^{\rm max}$) after a time $\Delta
t\sim\gamma_{\rm diff}^{-1}$. With its damping dominated by linear
dissipation, the rate of energy dissipation is ${\dot{E}}^{\rm max}_{\rm
orb}=\epsilon|C^{\rm max}|^{2}/\Delta t$.

However, if non-linearity sets in earlier at a smaller amplitude $C^{\rm
sat}$, the effective damping is reduced. The wave reaches this amplitude
within a shorter time, of order $\Delta t^{\prime}\sim\gamma_{\rm
diff}^{-1}|C^{\rm sat}|/|C^{\rm max}|$. And within a few radial propagation
times ($t_{\rm group}\sim R_{\star}[\partial\omega/\partial k_{r}]^{-1}\sim
n/\omega\ll\Delta t^{\prime}$), most wave flux carried to the centre is
converted to heat by overturning. This is consistent with results from Barker
& Ogilvie (2010); Su et al. (2020) where they found that wave breaking reduces
the amplitude of the outgoing wave to be a fraction of the ingoing wave (of
order $20\%-50\%$). The effective rate of energy loss is then $\dot{E}_{\rm
orb}^{\rm sat}\sim\epsilon|C^{\rm sat}|^{2}/\Delta t^{\prime}$ and is reduced
from the linear form as

$\frac{\dot{E}_{\rm orb}^{\rm sat}}{\dot{E}_{\rm orb}^{\rm
max}}\approx\frac{|C^{\rm sat}|}{|C^{\rm max}|}\sim 4.2\times
10^{-5}\,X_{k20}^{-1}(e)k^{-3/2}\left(\frac{P_{\rm orb}}{1\,{\rm
day}}\right)^{7/2}.$ (F8)

Following the procedure that leads to equation (43), we estimate a new binary
period out to which (saturated) resonance locking can be effective in
circularizing the binary. Setting again $t_{\rm evol}\sim 10^{7}\ {\rm yrs}$,
but now increasing $t_{\rm orb}$ in equation (32) by the inverse of the above
ratio, we find that $t_{\rm orb}\leq t_{\rm evol}$ when

$P_{\rm orb}\leq P_{\rm sat}\sim 3.6\,{\rm
days}\,\left({e\over{0.4}}\right)^{6/41}\,.$ (F9)

Compared to the criterion for linear resonance locking to be effective (eq.
[43]), saturated resonance locking can only circularize binaries out to a
shorter orbital period.

## References

  * Barker (2020) Barker, A. J. 2020, MNRAS, 498, 2270. doi:10.1093/mnras/staa2405
  * Barker (2011) Barker, A. J. 2011, MNRAS, 414, 1365
  * Barker & Ogilvie (2011) Barker, A. J., & Ogilvie, G. I. 2011, MNRAS, 417, 745
  * Barker & Ogilvie (2010) Barker, A. J., & Ogilvie, G. I. 2010, MNRAS, 404, 1849
  * Barker & Ogilvie (2009) Barker, A. J., & Ogilvie, G. I. 2009, MNRAS, 395, 2268
  * Batygin & Adams (2013) Batygin, K. & Adams, F. C. 2013, ApJ, 778, 169
  * Bildsten et al. (1996) Bildsten, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, 460, 827
  * Bluhm et al. (2016) Bluhm, P., Jones, M. I., Vanzi, L., et al. 2016, A&A, 593, A133
  * Bouvier (2013) Bouvier, J. 2013, EAS Publications Series, 143
  * Burkart et al. (2014) Burkart, J., Quataert, E., & Arras, P. 2014, MNRAS, 443, 2957
  * Burkart et al. (2012) Burkart, J., Quataert, E., Arras, P., et al. 2012, MNRAS, 421, 983
  * Chandrasekhar (1939) Chandrasekhar, S. 1939, Chicago, Ill., The University of Chicago press [1939]
  * Chapman & Lindzen (1970) Chapman, S., & Lindzen, R. 1970, Dordrecht: Reidel
  * Christensen-Dalsgaard (2014) Christensen-Dalsgaard, J. 2014, Lecture Notes on Stellar Oscillations, Fifth Edition, http://astro.phys.au.dk/jcd/oscilnotes/
Lecture_Notes_on_Stellar_Oscillations.pdf

  * Cowling (1941) Cowling, T. G. 1941, MNRAS, 101, 367
  * Darwin (1879) Darwin, G. H. 1879, Phil. Trans. Roy. Soc., 170, 1
  * Duquennoy & Mayor (1991) Duquennoy, A., & Mayor, M. 1991, A&A, 500, 337
  * Duguid et al. (2020) Duguid, C. D., Barker, A. J., & Jones, C. A. 2020, MNRAS, 491, 923
  * Friedman & Schutz (1978a) Friedman, J. L., & Schutz, B. F. 1978a, ApJ, 221, 937
  * Friedman & Schutz (1978b) Friedman, J. L., & Schutz, B. F. 1978b, ApJ, 222, 281
  * Fuller et al. (2017) Fuller J., Hambleton K., Shporer A., Isaacson H., Thompson S., 2017, MNRAS, 472, L25
  * Fuller (2017) Fuller, J. 2017, MNRAS, 472, 1538
  * Fuller et al. (2016) Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867
  * Fuller & Lai (2012a) Fuller, J., & Lai, D. 2012a, MNRAS, 420, 3126