
\section{Introduction} \label{s:intro}

The GYRE\ stellar oscillation code
\citep{Townsend:2013,Townsend:2018,Goldstein:2020} is a open-source
software instrument that solves the linear, non-radial, non-adiabatic
oscillation equations for an input stellar model. Released in 2013, it
has been productively used to study of heat-driven oscillations in
$\gamma$ Doradus and $\delta$ Scuti pulsators
\citep[e.g.,][]{Van-Reeth:2022,Murphy:2022}, slowly-pulsating B stars
\citep[e.g.,][]{Michielsen:2021}, variable sub-dwarf B stars
\citep[e.g.][]{Silvotti:2022}, pulsating pre-main sequence stars
\citep[e.g.,][]{Steindl:2021}, DBV white dwarfs
\citep[e.g.,][]{Chidester:2021}, and hypothetical `dark' stars
\citep{Rindler-Daller:2021}; to explore stochastically excited
oscillations in solar-like, subgiant and red-giant stars
\citep[e.g.,][]{Nsamba:2021,LiTanda:2020,LiTanda:2022}; and to model
oscillations of uncertain origin in red supergiant stars
\citep{Goldberg:2020}, post-outburst recurrent novae \citep{Wolf:2018}
and even gas-giant planets \citep{Mankovich:2019}.

This paper describes new functionality in GYRE\ for modeling static
and dynamic tides in binary\footnote{While our narrative focuses on
binary \emph{star} systems, it remains equally applicable to
star-planet systems.} systems. The equations governing small tidal
perturbations to one component of a binary are the linear
oscillation equations with extra terms representing the gravitational
forcing by the companion. Release 7.0 of GYRE\ implements these terms
and the supporting infrastructure necessary to solve the tidal
equations.

The view of astrophysical tides through the lens of forced
oscillations was pioneered in a pair of seminal papers by
\citet{Zahn:1970,Zahn:1975}. These papers also introduce complementary
approaches to solving the tidal equations, that we dub `modal
decomposition' (MD) and `direct solution' (DS). In MD the tidal
perturbations are decomposed as a superposition of the star's
free-oscillation modes, with weights determined from overlap integrals
between the mode eigenfunctions and the tidal force field. In DS the
two-point boundary value problem (BVP) posed by the tidal equations is
solved directly using a standard approach such as shooting or
relaxation. Examples of MD are given by \cite{Kumar:1995},
\citet{Lai:1997}, \citet{Schenk:2001}, \citet{Arras:2003} and \citet{Fuller:2012}; and of DS by
\citet{Savonije:1983,Savonije:1984}, \citet{Pfahl:2008}, and
\citet{Valsecchi:2013}. The study by \citet[][hereafter
  B12]{Burkart:2012} is noteworthy in that it adopts \emph{both}
approaches, although no direct comparison is made between them (a
lacuna that appears to extend into the wider literature).

The following section lays out the theoretical foundations for our
treatment of tides. Section~\ref{s:implement} describes the
modifications to GYRE\ to implement this formalism via a DS
methodology, and Section~\ref{s:calcs} presents illustrative calculations
focused on selected problems in the published
literature. Section~\ref{s:discuss} summarizes the paper, discusses
potential applications of the new GYRE\ functionality, and outlines
future improvements.


\section{Theoretical Formalism} \label{s:formalism}

Rather an exhaustive derivation of all equations, we opt to focus on
the key expressions that support our narrative and define the choices
(e.g., normalizations, sign conventions) dictated by the existing
numerical framework of GYRE. For more-detailed exposition,
we refer the reader to the papers by \citet{Polfliet:1990},
\citet{Smeyers:1991}, \citet{Smeyers:1998}, and
\citet{Willems:2003,Willems:2010}.

\subsection{Binary Configuration} \label{s:formalism-config}

Consider a binary system comprising a primary star of mass $M$
and photospheric radius $R$, together with a secondary star of
mass $q M$. To model the tides raised on the primary by the
secondary, we adopt a non-rotating reference frame with the primary's
center-of-mass fixed at the origin, the orbit of the secondary lying
in the Cartesian $xy$-plane, and the line of apsides coinciding with
the $x$-axis. The position vector of the secondary at time $t$ is then
given by
\begin{equation} \label{e:vsec}
  \mathbf{r}_{2}(t) = \frac{a(1 - e^{2})}{1 + e \cos\upsilon} \left( \mathbf{e}_{x} \cos
  \upsilon + \mathbf{e}_{y} \sin \upsilon \right),
\end{equation}
where $\mathbf{e}_{x}$ and $\mathbf{e}_{y}$ are the unit vectors along the $x$ and $y$
axes, respectively, $a$ is the orbital semi-major axis, and $e$ the
eccentricity. The true anomaly $\upsilon$ is linked to $t$ via Kepler's
equation
\begin{equation} \label{e:kepler}
  \mathcal{E} - e \sin\mathcal{E} = \mathcal{M}
\end{equation}
and the auxiliary relations
\begin{gather}
  \mathcal{M} = \Omega_{\rm orb} (t - t_{0}), \\
  \tan \frac{\upsilon}{2} = \sqrt{\frac{1-e}{1+e}} \tan \frac{\mathcal{E}}{2},
\end{gather}
where $t_{0}$ is a time of periastron passage, and $\mathcal{M}$ and
$\mathcal{E}$ are the mean and eccentric anomalies, respectively.  The
orbital angular frequency $\Omega_{\rm orb}$ is given by Kepler's third law,
\begin{equation}
  G M (1 + q) = a^{3} \Omega_{\rm orb}^{2},
\end{equation}
with $G$ the gravitational constant.

\subsection{Tidal Potential} \label{s:formalism-pot}

Tides are raised on the primary star by the forces arising from the
gravitational potential $\Phi_{2}$ of the secondary, which at position
vector $\mathbf{r}$ and time $t$ is
\begin{equation}
  \Phi_{2}(\mathbf{r};t) = - \frac{q G M}{|\mathbf{r}_{2} - \mathbf{r}|}.
\end{equation}
Using a multipolar expansion in space and a Fourier-series expansion
in time, this expression can be recast as
\begin{multline} \label{e:pot-expand}
  \Phi_{2}(\mathbf{r};t) =
  - \frac{q G M}{r_{2}}
  - \frac{q G M}{r_{2}^{2}} r \sin\vartheta \cos(\varphi - \upsilon) + \mbox{} \\
  \pottide(\mathbf{r};t)
  \end{multline}
where $(r,\vartheta,\varphi)$ are the spherical-polar radius,
colatitude and azimuth coordinates corresponding to $\mathbf{r}$. The first
term on the right-hand side is constant and therefore does not
generate a force. The second term produces a spatially uniform force
directed from the primary star toward the secondary, and precisely
cancels the inertial force arising from the orbital motion of the
primary center-of-mass about the system center-of-mass. The third term
represents the tidal part of the secondary potential, and is expressed
as a superposition
\begin{equation} \label{e:pottide}
  \pottide(\mathbf{r};t) =
  \sum_{\ell=2}^{\infty} \sum_{m=-\ell}^{\ell} \sum_{k=-\infty}^{\infty}
  \pottide[\ell,m,k] (\mathbf{r};t)
\end{equation}
of partial tidal potentials defined by
\begin{multline} \label{e:pottide-part}
  \pottide[\ell,m,k](\mathbf{r};t) =
  - \varepsilon_{\rm T} \,
  \frac{G M}{R} \,
  \cbar[\ell,m,k]
  \left( \frac{r}{R} \right)^{\ell} \times \mbox{} \\
  Y^{m}_{\ell}(\vartheta,\varphi) \,
  \ee{- \mathrm{i} k \mathcal{M}}.
\end{multline}
Here,
\begin{equation}
  \varepsilon_{\rm T} \equiv
  \left( \frac{R}{a} \right)^{3} \, q =
  \frac{\Omega_{\rm orb}^{2} R^{3}}{G M} \, \frac{q}{1 + q}
\end{equation}
is a dimensionless parameter that quantifies the overall strength of
the tidal forcing, $\cbar[\ell,m,k]$ is an tidal expansion coefficient
(Appendix~\ref{a:exp-coeff}), and $Y^{m}_{\ell}$ is a spherical
harmonic (Appendix~\ref{a:sph-harm}).

\subsection{Tidal Response} \label{s:formalism-resp}

Appendix~\ref{a:lin-eqns} summarizes the set of linearized equations
governing the response of the primary star to the tidal potential
$\pottide$. Based on the form~(\ref{e:pottide-part}) of the partial
potentials, solutions to these equations take the form
\begin{multline} \label{e:sol-xi}
  \boldsymbol{\xi}(\mathbf{r};t) =
  \sum_{\ell,m,k}
  \left[
    \txir[\ell,m,k](r) \, \mathbf{e}_{r} + \phantom{\left( \mathbf{e}_{\vartheta} \, \pderiv{}{\vartheta} \right)}
  \right. \\
  \left.
    \txih[\ell,m,k](r) \left( \mathbf{e}_{\vartheta} \, \pderiv{}{\vartheta} + \frac{\mathbf{e}_{\varphi}}{\sin\vartheta} \, \pderiv{}{\varphi}\right)
  \right] \,
  \angtime[\ell,m,k](\vartheta,\varphi;t)
\end{multline}
for the displacement perturbation vector $\boldsymbol{\xi}$, and
\begin{equation} \label{e:sol-f}
  f'(\mathbf{r};t) =
  \sum_{\ell,m,k}
  \tf[\ell,m,k]'(r) \, \angtime(\vartheta,\varphi;t)
\end{equation}
for the Eulerian ($f'$) perturbation to a scalar variable $f$ (the
corresponding Lagrangian perturbation $\delta f$ follows from
equation~\ref{e:lag-eul}). In these expressions, the notation
$\sum_{\ell,m,k}$ abbreviates the triple sum appearing in
equation~(\ref{e:pottide}), while $\mathbf{e}_{r}$, $\mathbf{e}_{\vartheta}$ and $\mathbf{e}_{\varphi}$ are the
unit basis vectors in the radial, polar, and azimuthal directions,
respectively. The functions
\begin{equation} \label{e:angtime}
  \angtime[\ell,m,k](\vartheta,\varphi;t) \equiv Y^{m}_{\ell}(\vartheta,\varphi) \ee{-\mathrm{i} k \mathcal{M}}
\end{equation}
describe the angular and time dependence of the response, while the
functions with tilde accents ($\txir,\txih,\tf'$) encapsulate the
radial dependence. The latter are found as solutions to a system of
tidal equations summarized in
Appendix~\ref{a:tidal-eqns}. Importantly, the set of radial functions
for a given combination of indices $\{\ell,m,k\}$ can be determined
independently of any other combination.

\section{Implementation in GYRE} \label{s:implement}

To implement the tidal equations~(\ref{e:tidal-mass}--\ref{e:tidal-diff}) in GYRE, which follows a DS methodology, we first
transform to a dimensionless independent variable $x \equiv r/R$
and a set of dimensionless dependent variables
\begin{equation}
  \begin{aligned}
    y_{1} &\equiv x^{2 - \ell} \frac{\txir}{r}, \\
    y_{2} &\equiv x^{2 - \ell} \frac{\tP'}{\rho g r}, \\
    y_{3} &\equiv x^{2 - \ell} \frac{\tpottot'}{gr}, \\
    y_{4} &\equiv x^{2 - \ell} \frac{1}{g} \deriv{\tpottot'}{r}, \\
    y_{5} &\equiv x^{2 - \ell} \frac{\delta \tS}{c_{P}}, \\
    y_{6} &\equiv x^{-1 - \ell} \frac{\delta \tLrad}{L}.
  \end{aligned}
\end{equation}
Here $g \equiv \sderiv{\Phi}{r}$ is the scalar gravity, and the other symbols are defined in Appendix~\ref{a:lin-eqns} (for notational simplicity we neglect the $\ell,m,k$ subscripts from perturbed quantities). With these transformations, the differential equations and boundary conditions
governing $y_{1},\ldots,y_{6}$ can be written in a form almost
identical to the linear, non-radial, non-adiabatic free oscillation equations detailed in
Appendix B2 of \citet{Townsend:2018}. The only differences are that
the interpretation of the $y_{3}$ and $y_{4}$ variables is altered;
the outer boundary condition governing the gravitational potential
acquires an inhomogeneous term on the right-hand side
\begin{equation} \label{e:inhom-bc}
 U y_{1} +  (\ell + 1) y_{3} + y_{4} = (2 \ell + 1) y_{\rm T},
\end{equation}
where $U$ is the usual homology invariant \citep[e.g.,][]{Kippenhahn:2013} and
\begin{equation} \label{e:y-T}
  y_{\rm T} \equiv x^{2 - \ell} \frac{\tpottide[\ell,m,k]}{gr};
\end{equation}
and the dimensionless oscillation frequency is defined by
\begin{equation} \label{e:omega}
  \omega \equiv \sqrt{\frac{R^{3}}{G M}} \sigma_{m,k},
\end{equation}
where $\sigma_{m,k}$ is the rotating-frame frequency defined in
equation~(\ref{e:sigmk}).

GYRE\ uses a multiple shooting algorithm to discretize BVPs on a
radial grid $x = x_{1}, \ldots, x_{N}$. For a system of of $n$
differential equations and $n$ boundary conditions (in the present
case, $n=6$) this leads to a corresponding system of linear algebraic
equations with the form
\begin{equation} \label{e:lin-sys}
  \mathsf{\mathbf{S}} \mathbf{u} = \mathbf{b},
\end{equation}
where $\mathsf{\mathbf{S}} \in \mathbb{C}^{N n \times N n}$ is a block staircase
matrix and the solution vector $\mathbf{u} \in \mathbb{C}^{N n}$ contains the
dependent variables $y_{1},\ldots,y_{n}$ evaluated at successive grid
points. For free-oscillation problems the right-hand side vector $\mathbf{b}
\in \mathbb{C}^{N n}$ is identically zero, and so the linear
system~(\ref{e:lin-sys}) is homogeneous; non-trivial solutions exist
only when the determinant of $\mathsf{\mathbf{S}}$ vanishes, a condition that
ultimately determines the dimensionless eigenfrequencies $\omega$ of
the star. However, for the tidal problem considered here, $\mathbf{b}$ is
non-zero due to the inhomogeneous term appearing on the right-hand
side of the boundary condition~(\ref{e:inhom-bc}). The linear system
can then be solved for any choice of $\omega$ --- that is, the
dimensionless frequency is an \emph{input} rather than an output.

The reordered flow of execution ($\omega$ as input rather than output)
motivates the decision to provide two separate executables in release
7.0 of GYRE: \texttt{gyre}\ for modeling free stellar oscillations (the same
as in previous releases), and \texttt{gyre\_tides}\ for modeling stellar tides
(new to this release). A full description of these programs, including
their input parameters and output data, is provided on the
GYRE\ documentation
site\footnote{\url{https://gyre.readthedocs.io/en/stable/}}


\section{Example Calculations} \label{s:calcs}

\subsection{Surface Perturbations in a KOI-54 Model} \label{s:calcs-surf}

\begin{table}
  \begin{centering}
  \begin{tabular}{cccccc}
    Model  & $M$ (\si{M_{\odot}}) & $R$ (\si{R_{\odot}}) & $T_{\rm eff}$ (\si{\kelvin}) & $Z_{\rm i}$ & $X_{\rm c}$ \\ \hline
    KOI-54 & 2.32 & 2.19 & \num{9400}  & 0.0328 & 0.487 \\
    B-star & 5.00 & 2.80 & \num{16700} & 0.0200 & 0.658 \\
  \end{tabular}
  \caption{Fundamental parameters for the two stellar models
    discussed in the text: mass $M$, radius $R$, effective
    temperature $T_{\rm eff}$, initial metal mass fraction $Z_{\rm i}$ and center hydrogen mass fraction $X_{\rm c}$.} \label{t:models}
  \end{centering}
\end{table}

\begin{figure}
  \includegraphics{fig_prop_diag_koi54}
  \caption{Propagation diagram for the KOI-54 model discussed in the
    text, plotting the Brunt-V\"ais\"al\"a\ ($N$) and Lamb ($S_{\ell}$) frequencies as a
    function of pressure coordinate for $\ell=2$ modes. The hatched
    regions indicate the acoustic wave ($\sigma > N,S_{\ell}$) and
    internal gravity wave q($\sigma < N,S_{\ell}$) propagation
    regions. The sloped dashed line shows the frequency corresponding
to the local thermal timescale $t_{\rm thm}$ (equation~\ref{e:t-thm}), while the upper (lower) horizontal dashed line corresponds to a frequency $\sigma/\Omega_{\rm orb} = 30$ (15).}
  \label{f:prop-diag-koi54}
\end{figure}

KOI-54 is a highly-eccentric ($e \approx 0.8$), near face-on binary
system comprising a pair of A-type stars. Discovered by
\citet{Welsh:2011} in Kepler observations, it serves as the archetype
of the heartbeat class of periodic variables
\citep{Thompson:2012}. B12 and \citet{Fuller:2012} each present
initial attempts to model the light curve of KOI-54 in terms of
contributions from the equilibrium tide, dynamical tides and stellar
irradiation. Here we undertake a calculation to reproduce Fig.~6 of
B12, which illustrates how a model for the KOI-54 primary responds to
forcing by a single partial tidal potential.

Guided by the parameters given in Table~1 of B12, we use release
r22.11.1 of the MESA\ software instrument
\citep{Paxton:2011,Paxton:2013,Paxton:2015,Paxton:2018,Paxton:2019,Jermyn:2022}
to evolve a \SI{2.32}{M_{\odot}} model from the zero-age main sequence
(ZAMS) until its photospheric radius has expanded to
\SI{2.19}{R_{\odot}}. Input files for this and subsequent
MESA\ calculations are available through Zenodo at
\dataset[10.5281/zenodo.7489814]{\doi{10.5281/zenodo.7489814}}. The
growth of the core is followed using the convective premixing
algorithm \citep{Paxton:2019}, but rotation is neglected. We avoid any
smoothing of the Brunt-V\"ais\"al\"a\ frequency profile, as this can introduce small
but consequential departures from mass conservation. Fundamental
parameters for the model are summarized in Table~\ref{t:models}, and
its propagation diagram is plotted in Fig.~\ref{f:prop-diag-koi54}.

\begin{figure*}
  \includegraphics{fig_scan_sigma}
  \caption{Plot of the radial displacement perturbation $\txir$ and
    the Lagrangian radial flux perturbation $\delta \tFradr$ at the
    surface of the KOI-54 model, as a function of forcing frequency
    $\sigma$. Selected peaks are labeled by the resonant mode's
    classification within the Eckart-Osacki-Scuflaire scheme.}
  \label{f:scan-sigma}
\end{figure*}

Figure 6 of B12 plots the modulus of $\txir$ and $\delta \tFradr$ as a
function of $\sigma/\Omega_{\rm orb}$, for a single partial potential with
$\ell=2$. The choices of $m$ and $k$ are left undetermined because B12
neglect rotation in their figure, and treat $\sigma$ as a free
parameter rather than being constrained by
equation~(\ref{e:sigmk}). They normalize the strength of the partial
potential in a manner equivalent to setting $y_{\rm T} = 1$ in
equation~(\ref{e:inhom-bc}). We use \texttt{gyre\_tides}\ to repeat these steps
for our KOI-54 model, plotting the results in Fig.~\ref{f:scan-sigma}.

This figure shows qualitative agreement with Fig.~6 of B12. For
$\sigma/\Omega_{\rm orb} \gtrsim 30$ the surface perturbations exhibit distinct
peaks, corresponding to resonances with the star's $\ell=2$
free-oscillation modes. Selected peaks are labeled with the resonant
mode's classification within the Eckart-Osaki-Scuflaire scheme
\citep[e.g.,][]{Unno:1989}. For $\sigma/\Omega_{\rm orb} \lesssim 30$ the peaks
merge together and dissolve, because the periods of the resonant modes
become appreciably shorter than the local thermal timescale
\begin{equation} \label{e:t-thm}
  t_{\rm thm}(r) \equiv \int_{r}^{R} \frac{4 \pi r^{2} \rho T c_{P}}{L} \diff{r}
\end{equation}
in the outer parts of the main mode-trapping cavity (see
Fig.~\ref{f:prop-diag-koi54}), resulting in significant
radiative-diffusion damping that broadens and suppresses the
resonances. Eventually, for $\sigma/\Omega_{\rm orb} \lesssim 15$ the surface
perturbations reach the limits $\txir/R \rightarrow 1$ and $\delta
\tFradr/F_{{\rm rad},r} \rightarrow 4$ corresponding to the $\ell=2$
equilibrium tide (see Section~6.2 of B12 for a discussion of these
limits).

On closer inspection, some differences between the two figures are
apparent. While the KOI-54 stellar models in B12 and the present work
are similar (in particular, sharing the same $M$ and $R$),
they are not identical; therefore, the locations and heights of the
resonance peaks are not the same in each figure. More significantly,
over the range $5 \lesssim \sigma/\Omega_{\rm orb} \lesssim 30$ the flux
perturbation behaves much more smoothly in Fig.~\ref{f:scan-sigma}
than in Fig.~6 of B12. The reason for this discrepancy is not obvious,
but we speculate that it may be linked to differences in the
near-surface convection of the models. The propagation diagram for our
KOI-54 model (Fig.~\ref{f:prop-diag-koi54}) shows a He\textsc{ii}
convection zone at $\log(P/\si{\dyne\per\cm\squared}) \approx 5.5$,
but this zone is absent in B12's model (cf. their Fig.~1).


\subsection{Direct Solution versus Mode Decomposition in a KOI-54 Model} \label{s:calcs-ds-vs-md}

\begin{figure*}
  \includegraphics{fig_scan_sigma_zoom}
  \caption{Zoom-in on the four labeled resonance peaks in
    Fig.~\ref{f:scan-sigma}, plotting the complex amplitude (upper
    panels; rescaled to have a maximum value of unity) and phase
    (lower panels) of $\txir/R$ as a function of normalized detuning
    parameter $\epsilon$ (equation~\ref{e:detune}). Separate curves are
    shown for the DS (\texttt{gyre\_tides}) and MD approaches.}
  \label{f:scan-sigma-zoom}
\end{figure*}

\begin{figure}
  \includegraphics{fig_eigfuncs}
  \caption{Real (upper panel) and imaginary (lower panel) parts of the
    wavefunctions associated with the radial displacement perturbation
    $\txir$ at the peak of the $\gmode{35}$ resonance, plotted as a
    function of pressure coordinate. The separate curves correspond to
    the DS (\texttt{gyre\_tides}) and MD approaches.}
  \label{f:eigfuncs}
\end{figure}

As a validation of the results presented in Fig.~\ref{f:scan-sigma},
we re-calculate the surface perturbations using the MD approach. We
follow the formalism laid out in Section 3.2 of B12, although adopting
non-adiabatic eigenfrequencies $\sigma_{n,\ell}$ and damping rates $\gamma_{n,\ell}$
provided by \texttt{gyre}\ to evaluate the Lorentzian factor
$\Delta_{n,\ell,m,k}$ (their equation 13; here, $n$ is the mode radial
order). To evaluate the overlap integrals $Q_{n,\ell}$ that weight the
contribution of each free-oscillation mode in the MD superposition, we
use the third expression of equation~(9) in B12; we find that the the
first expression yields unreliable values when $|n| \gtrsim 20$,
because the integrand is highly oscillatory and suffers from
significant cancellation.

Figure~\ref{f:scan-sigma-zoom} zooms in on the four labeled
resonance peaks from Fig.~\ref{f:scan-sigma}, plotting the
complex amplitude and phase of $\txir$ as a function of normalized
detuning parameter
\begin{equation} \label{e:detune}
  \bar{\epsilon} \equiv \frac{\sigma - \sigma_{n,\ell}}{|\gamma_{n,\ell}|},
\end{equation}
for the two approaches.  DS (i.e., \texttt{gyre\_tides}) and MD agree at higher
forcing frequencies (right-hand panels), but show mismatches toward
lower frequencies (left-hand panels).

To delve further into these discrepancies, Fig.~\ref{f:eigfuncs} plots
the $\txir$ wavefunction evaluated at the peak ($\bar{\detune}=0$) of the
$\gmode{35}$ resonance. The imaginary part of the wavefunction is
spatially oscillatory because it is dominated by the dynamical tide,
comprising the resonantly forced oscillation mode. The real part is
non-oscillatory and corresponds to the equilibrium tide, comprising
the superposition of the other, non-resonant $\ell=2$ modes.

The discrepancies in the wavefunction are restricted to the outer
parts of the stellar envelope. For $\operatorname{Re}(\txir)$, the DS and MD
curves begin to diverge at $\log (P/\si{\dyne\per\cm\squared}) \approx
8$ (corresponding to $r/R \gtrsim 0.96$), while for $\operatorname{Im}(\txir)$
the divergence begins further out at $\log
(P/\si{\dyne\per\cm\squared}) \approx 6$ ($r/R \approx 0.99$). We
hypothesize that these divergences arise because $\sigma \ll
2\pi/t_{\rm thm}$ in these superficial layers, leading to significant
non-adiabaticity that MD is unable to correctly reproduce (see Section
6.2 of B12; also, Section 6 of \citealp{Fuller:2017}).


\subsection{Circularization in a B-star Model} \label{s:calcs-circ}

\begin{figure}
  \includegraphics{fig_prop_diag_5msun}
  \caption{Propagation diagram for the B-star model introduced in the
    text (cf. Fig.~\ref{f:prop-diag-koi54}). The dip in the
    Brunt-V\"ais\"al\"a\ frequency around $\log (P/\si{\dyne\per\cm\squared})
    \approx 7.5$ is caused by the iron opacity bump responsible for
    the overstability of the $\gmode{9}-\gmode{17}$ modes.}
  \label{f:prop-diag-5msun}
\end{figure}

\begin{figure*}
  \includegraphics{fig_secular}
  \caption{The circularization timescale $t_{e}$ plotted as a
    function of orbital period $P_{\rm orb}$, for the binary system with the
    B-star primary. Resonances are labeled beneath with the mode
    classification, and in parentheses the harmonic index $k$ and the
    sense of propagation in the co-rotating frame (p=prograde,
    r=retrograde). If the resonance is with an overstable mode, then
    the label border is bolded.}
  \label{f:secular}
\end{figure*}

\begin{figure}
  \includegraphics{fig_eigfreqs}
  \caption{Damping rates $\gamma_{n,\ell}$ plotted against
    eigenfrequencies $\sigma_{n,\ell}$ for $\gmode{1}$--$\gmode{20}$ modes of the
    B-star model, as calculated using \texttt{gyre}\ and as tabulated
    by W03. Open symbols indicate stable ($\gamma_{n,\ell} > 0$) modes, and
    filled symbols overstable ($\gamma_{n,\ell} < 0$) modes.}
  \label{f:eigfreqs}
\end{figure}

\citet[][hereafter W03]{Willems:2003} explore the secular orbital
changes due to tides in a binary system comprising a \SI{5}{M_{\odot}}
B-star primary and a \SI{1.4}{M_{\odot}} secondary on an $e=0.5$
orbit. They adopt the MD approach, but include only a single term at a
time in the modal superposition. The calculations illustrated in
Fig.~3 of W03 are used by \citet{Valsecchi:2013} to benchmark their
CAFein\ code, motivating us to do likewise here with \texttt{gyre\_tides}.

We use MESA\ as before to evolve a \SI{5}{M_{\odot}} model from the ZAMS
until its photospheric radius has grown to match the \SI{2.80}{R_{\odot}}
of W03's model. Fundamental parameters for this model are
summarized in Table~\ref{t:models}, and its propagation diagram is
plotted in Fig.~\ref{f:prop-diag-5msun}. Then, we apply
\texttt{gyre\_tides}\ to evaluate the star's response to the dominant
contributions in the tidal potential~(\ref{e:pottide}), comprising the
$\ell=2$ terms with $|m| = 2$ and $m k > 0$. We further restrict the
summation over $k$ to terms with a magnitude at least $10^{-12}$ times
that of the largest term.  As in W03, we assume the stellar angular
rotation frequency is equal to the periastron angular velocity of the
secondary,
\begin{equation}
  \Omega_{\rm peri} = \Omega_{\rm orb} \, \sqrt{\frac{(1+e)}{(1-e)^{3}}}.
\end{equation}
Based on this configuration, Fig.~\ref{f:secular} plots the circularization
timescale
\begin{equation}
  t_{e} \equiv - \left[ \frac{1}{e} \left( \deriv{e}{t} \right) \right]^{-1}_{\rm sec}
\end{equation}
as a function of orbital period $P_{\rm orb} \equiv 2\pi/\Omega_{\rm orb}$ over a pair
of intervals\footnote{These correspond to the short- and long-period
limits in Fig.~3 of W03.}. The secular rate-of-change of eccentricity
is evaluated via
\begin{multline} \label{e:secular-edot}
  \left( \deriv{e}{t} \right)_{\rm sec} = 4 \,
  \Omega_{\rm orb} \,
  q
  \sum_{\ell,m,k \geq 0}
  \left( \frac{R}{a} \right)^{\ell+3} \,
  \left( \frac{r_{\rm s}}{R} \right)^{\ell+1} \\
  \mbox{} \times \kap[\ell,m,k] \,
  \operatorname{Im}(\Fbar[\ell,m,k]) \,
  \Gbar[\ell,m,k]{3};
\end{multline}
this comes from equation~(55) of \citet{Willems:2010}, with
$F_{\ell,m,k}$ replaced by $\Fbar[\ell,m,k] \equiv F_{\ell,m,-k}$ and
$G^{(3)}_{\ell,m,k}$ by $\Gbar[\ell,m,k]{3} \equiv G^{(3)}_{\ell,m,-k}$ to account for
the differing sign convention in the assumed time dependence of
partial tides (see equation~\ref{e:angtime}). Note that the summation
is now restricted to $k \geq 0$.

The $|t_{e}|$ data plotted in the figure vary smoothly with $P_{\rm orb}$
between a series of sharp extrema. The maxima (more correctly,
singularities) arise when $(\sderiv{e}{t})_{\rm sec}$ passes through
zero. The minima arise from resonances with the star's $\ell=2$
free-oscillation modes, similar to the peaks seen in
Fig.~\ref{f:scan-sigma}. However, a key difference here is that the
star is being forced with a superposition~(\ref{e:pottide}) of partial
tidal potentials, rather than a single one as before. The criterion
for resonance $\sigma_{m,k} \approx \sigma_{n,\ell}$ can be satisfied for many
different values of $k$, leading to multiple resonances with the same
mode. This can be seen in the figure; for instance, the upper panel
shows resonances between the $f$ mode and the partial tidal potentials with
$k=31,\ldots,35$.

In the vicinity of some of the resonances shown in the figure, $t_{e}
< 0$ (blue) indicates that the tide acts to increase the eccentricity
of the orbit, driving it further away from circular. This behavior is
an instance of the `inverse tides' phenomenon discussed by
\citet{Fuller:2021}, and occurs when the summation in
equation~(\ref{e:secular-edot}) is dominated by a single,
\emph{positive} term. There are four distinct configurations that lead
to this outcome:
\begin{enumerate}[label=\Roman*.]
\item $\Gbar[\ell,m,k]{3} > 0$ and $\operatorname{Im}(\Fbar[\ell,m,k]) > 0$, the latter because
  \begin{enumerate}
  \item the resonant mode is prograde in the co-rotating frame ($\sigma_{n,\ell}/m > 0$)
    and stable ($\gamma_{n,\ell} > 0$); or \label{l:pos-pro-s}
  \item the resonant mode is retrograde in the co-rotating frame
    ($\sigma_{n,\ell}/m < 0$) and overstable ($\gamma_{n,\ell} < 0$). \label{l:pos-ret-o}
  \end{enumerate}
\item $\Gbar[\ell,m,k]{3} < 0$ and $\operatorname{Im}(\Fbar[\ell,m,k]) < 0$, the latter because
  \begin{enumerate}
  \item the resonant mode is prograde in the co-rotating frame ($\sigma_{n,\ell}/m > 0$)
    and overstable ($\gamma_{n,\ell} < 0$); or \label{l:neg-pro-o}
  \item the resonant mode is retrograde in the co-rotating frame
    ($\sigma_{n,\ell}/m < 0$) and stable ($\gamma_{n,\ell} > 0$). \label{l:neg-ret-s}

  \end{enumerate}
\end{enumerate}
All of the $t_{e} < 0$ resonances seen in Fig.~\ref{f:secular} are
instances of cases~(\ref{l:pos-ret-o}) or~(\ref{l:neg-pro-o}), and
therefore involve overstable modes. The overstability is caused by the
iron-bump opacity mechanism responsible for the slowly pulsating B
(SPB) stars \citep[e.g.,][]{Gautschy:1993,Dziembowski:1993}; in the
B-star model considered here, which falls well inside the SPB
instability strip \citep[e.g.,][]{Pamyatnykh:1999,Paxton:2015}, this
mechanism excites the $\ell=2$ $\gmode{9}$--$\gmode{17}$ modes.

Comparing Fig.~\ref{f:secular} against Fig.~3 of W03 reveals some
important differences. The latter shows numerous gaps and
discontinuities in $t_{e}$, that appear to arise because W03 only
allow a given mode to contribute toward the MD superposition when its
detuning parameter (equation~\ref{e:detune}) satisfies $0.1\,\varepsilon_{\rm T}
\leq |\bar{\detune} \gamma_{n,\ell}/ \sigma_{n,\ell}| \leq 10\,\varepsilon_{\rm T}$. The lower bound on
$\bar{\detune}$ means that the central parts of each resonance are omitted,
and so Fig.~3 of W03 does not fully reveal how small $t_{e}$ can
become when close to a resonance.

At the short-period limit of the range shown in the figures, there is
also disagreement between the typical magnitude of $t_{e}$ between
the resonances; Fig.~\ref{f:secular} shows an inter-resonance
$|t_{e}| \approx \SI{e9}{\year}$ , whereas for W03 it is 2--4 orders
of magnitude shorter. This is likely a consequence of W03
over-estimating the damping for the $\gmode{1}$--$\gmode{4}$ modes, which
dominate the tidal response at short orbital
periods. Figure~\ref{f:eigfreqs} plots eigenfrequencies $\sigma_{n,\ell}$ and
damping rates $\gamma_{n,\ell}$ for the $\gmode{1}$--$\gmode{20}$ modes, as calculated
using \texttt{gyre}\ and as obtained from Table~1 of W03. The data for the
$\gmode{5}$--$\gmode{16}$ modes are in reasonable agreement, especially given
that the underlying stellar models do not have the exact same internal
structure. However, the W03 damping rates for the $\gmode{1}$--$\gmode{4}$ modes
are four-to-five orders of magnitude larger than the \texttt{gyre}\ ones.


\subsection{Pseudo-Synchronization in a KOI-54 Model} \label{s:calcs-sync}

\begin{figure}
  \includegraphics{fig_torque}
  \caption{Secular torque $\torque_{\rm sec}$ plotted against stellar
    angular rotation frequency $\Omega_{\rm rot}$ for the KOI-54 primary
    model. The dashed line shows the equilibrium-tide torque extracted
    from Fig. 4 of B12.}
  \label{f:torque}
\end{figure}

B12 explore how tides can modify the rotation of the primary star in
the KOI-54 system, by evaluating the secular tidal torque $\torque_{\rm sec}$ on
the primary star as a function of the star's rotation rate. Using an
MD approach (see their Appendix C), they consider contributions toward
the torque from the $\ell=2$ partial tidal potentials. Their Fig. 4
shows a smoothly varying $\torque_{\rm sec}$ punctuated by many narrow peaks due
to modal resonances. The smooth part corresponds to the torque from
the equilibrium tide, and passes through zero at the
pseudo-synchronous angular frequency
\begin{equation}
  \Omega_{\rm ps} \equiv \Omega_{\rm orb} \frac{1 + (15/2)e^{2} + (45/8)e^{4} + (5/16)e^{6}}{[1 + 3e^{2} + (3/8)e^4](1 - e)^{3/2}}
\end{equation}
first derived by \citet{Hut:1981} in his theory of tides in the weak
friction limit. For the orbital parameters of KOI-54, the
pseudo-synchronous frequency is $\Omega_{\rm ps}/(2\pi) = \SI{0.395}{\per\day}$.

To repeat this calculation, we use \texttt{gyre\_tides}\ to evaluate the
response of the KOI-54 model (Section~\ref{s:calcs-surf}) to the
$\ell=2$ terms in the tidal potential~(Equation~\ref{e:pottide}) for
\num{25000} rotation frequencies in the interval $\SI{0}{\per\day}
\leq \Omega_{\rm rot}/(2\pi) \leq \SI{1}{\per\day}$. As in the preceding section,
we restrict the summation over $k$ to terms with a magnitude at least
$10^{-12}$ times that of the largest term. Then, we evaluate the tidal
torque via
\begin{multline}
  \torque_{\rm sec} = 
  4 \Omega_{\rm orb} \sqrt{\frac{G M^{3} q^{2}}{1 + q}} \, q \, a^{1/2}
  \sum_{\ell,m,k \geq 0}
  \left( \frac{R}{a} \right)^{\ell+3} \,
  \left( \frac{r_{\rm s}}{R} \right)^{\ell+1} \\
  \times \kap[\ell,m,k] \,
  \operatorname{Im}(\Fbar[\ell,m,k]) \,
  \Gbar[\ell,m,k]{4};
\end{multline}
this comes from equation~(63) of \citet{Willems:2010}, with
$G^{(4)}_{\ell,m,k}$ replaced by $\Gbar[\ell,m,k]{4} \equiv
G^{(4)}_{\ell,m,-k}$.

Fig.~\ref{f:torque} plots $\torque_{\rm sec}$ as a function of $\Omega_{\rm rot}$, using
discrete points rather than a continuous line because we are
significantly undersampling the dense forest of resonances. Also
plotted for comparison is the smooth (equilibrium tide) part of the
torque extracted from Fig.~4 of B12.  Clearly, there are some
significant differences between the two figures. Ours shows a torque
that's generally positive torque at small rotation frequencies, and
negative at high frequencies; however, the switch-over point is not
nearly as sharply defined as in B12, and occurs at a frequency
$\Omega_{\rm rot}/(2\pi) \approx \SI{0.5}{\per\day}$ around 25\% higher than
$\Omega_{\rm ps}$. Moreover, on either side of the switch-over, the lower
envelope of our torque values is around two orders of magnitude
smaller than the B12 curve.

Exploratory calculations indicate that these differences are not a
result of inaccurate overlap integrals (as was the case in
Section~\ref{s:calcs-ds-vs-md}), but rather due to a genuine
incompatibility between the DS and MD approaches. While a detailed
investigation of this problem is beyond the scope of the present
paper, we believe the fault lies with MD's over-estimation of
radiative dissipation for the equilibrium tide. If this hypothesis is correct, an immediate corollary is that pseudo-synchronization as envisaged by \citet{Hut:1981} does not operate for stars in which radiative dissipation dominates the tidal damping.


\section{Summary and Discussion} \label{s:discuss}

To briefly summarize the preceding sections: we establish the
theoretical foundations for our tides treatment
(Section~\ref{s:formalism}), describe modifications to GYRE\ to
implement tides (Section~\ref{s:implement}), and then apply the new
\texttt{gyre\_tides}\ executable to selected problems (Section~\ref{s:calcs}).
These example calculations uncover disagreements between the DS and MS
approaches, arising for a variety of reasons --- from numerical
inaccuracies in overlap integrals (Section~\ref{s:calcs-ds-vs-md}),
through to what appears to be an incompatibility between the
approaches. We look forward to future opportunities to investigate
these disagreements.

We also plan a number of enhancements to \texttt{gyre\_tides}\ that will extend
its capabilities. Key milestones include adding the ability to model
systems with spin-orbit misalignments \citep[e.g., following the
  formalism described by][]{Fuller:2017}; the incorporation of
additional damping mechanisms beyond radiative diffusion \citep[for
  instance, turbulent viscosity within convection zones; see
][]{Willems:2010}; and treating the effects of the Coriolis force,
which was neglected in deriving the linearized equations
(Appendix~\ref{a:lin-eqns}). Partial treatment of the Coriolis force
is already included in the main \texttt{gyre}\ executable, via the
traditional approximation of rotation (TAR; see, e.g.,
\citealp{Bildsten:1996,Lee:1997,Townsend:2003}). However, in its
current form this implementation cannot be used in \texttt{gyre\_tides},
because the angular operator appearing in the linearized continuity
equation within the TAR \citep[see equation~6 of][]{Bildsten:1996}
does not commute with the angular part of the Laplacian operator
appearing in the linearized Poisson equation~(\ref{e:lin-poisson}).

GYRE\ is not the first software package that adopts the DS approach
to model stellar tides; \citet{Pfahl:2008} and \citet{Valsecchi:2013}
describe functionally similar codes. Although the former authors' code
has never been made publicly available, the latters' CAFein\ code is
accessible on
GitHub\footnote{https://github.com/FrancescaV/CAFein}. After fixing a
number of bugs in CAFein\ (for instance, relating to misinterpreting
cell-centered quantities in MESA\ models as face-centered), we have
undertaken exploratory calculations comparing it against \texttt{gyre\_tides},
and find the two codes are in general agreement. Given that
CAFein\ is unmaintained, we decided that more-detailed comparison
would not be a worthwhile exercise.

It is our hope that \texttt{gyre\_tides}\ will provide a standardized and
well-supported community tool for simulating tides of spherical stars
within the linear limit. Specific areas where we anticipate productive
applications include modeling the many heartbeat systems discovered by
Kepler \citep[e.g.,][]{Thompson:2012}, TESS
\citep[e.g.,][]{Kolaczek:2021} and OGLE \citep{Wrona:2022};
investigating why most of these systems rotate faster than the
pseudo-synchronous rate \citep{Zimmerman:2017}; and exploring orbital
and rotational evolution in more-general star-star and star-planet
systems. These latter activities will initially be restricted to cases
where radiative diffusion dominates the tidal damping \citep[e.g., the
  $\gamma$ Doradus stars considered by ][]{LiGang:2020}; but with the
planned addition of other damping mechanisms, they can be extended
more broadly.


\section*{Acknowledgments}

This work has been supported by NSF grants ACI-1663696, AST-1716436
and PHY-1748958, and NASA grant 80NSSC20K0515. This research was also supported by STFC through grant ST/T00049X/1. The authors thank the referee for comments that have improved this paper.

\facilities{We have made extensive use of NASA's
  Astrophysics Data System Bibliographic Services.}

\software{Astropy \citep{astropy:2013,astropy:2018,astropy:2022},
  GYRE\ \citep{Townsend:2013,Townsend:2018,Goldstein:2020},
  Matplotlib \citep{Hunter:2007}, MESA\ 
  \citep{Paxton:2011,Paxton:2013,Paxton:2015,Paxton:2018,Paxton:2019,Jermyn:2022}}



\section{Introduction} \label{s:intro}

The GYRE\ stellar oscillation code
\citep{Townsend:2013,Townsend:2018,Goldstein:2020} is a open-source
software instrument that solves the linear, non-radial, non-adiabatic
oscillation equations for an input stellar model. Released in 2013, it
has been productively used to study of heat-driven oscillations in
$\gamma$ Doradus and $\delta$ Scuti pulsators
\citep[e.g.,][]{Van-Reeth:2022,Murphy:2022}, slowly-pulsating B stars
\citep[e.g.,][]{Michielsen:2021}, variable sub-dwarf B stars
\citep[e.g.][]{Silvotti:2022}, pulsating pre-main sequence stars
\citep[e.g.,][]{Steindl:2021}, DBV white dwarfs
\citep[e.g.,][]{Chidester:2021}, and hypothetical `dark' stars
\citep{Rindler-Daller:2021}; to explore stochastically excited
oscillations in solar-like, subgiant and red-giant stars
\citep[e.g.,][]{Nsamba:2021,LiTanda:2020,LiTanda:2022}; and to model
oscillations of uncertain origin in red supergiant stars
\citep{Goldberg:2020}, post-outburst recurrent novae \citep{Wolf:2018}
and even gas-giant planets \citep{Mankovich:2019}.

This paper describes new functionality in GYRE\ for modeling static
and dynamic tides in binary\footnote{While our narrative focuses on
binary \emph{star} systems, it remains equally applicable to
star-planet systems.} systems. The equations governing small tidal
perturbations to one component of a binary are the linear
oscillation equations with extra terms representing the gravitational
forcing by the companion. Release 7.0 of GYRE\ implements these terms
and the supporting infrastructure necessary to solve the tidal
equations.

The view of astrophysical tides through the lens of forced
oscillations was pioneered in a pair of seminal papers by
\citet{Zahn:1970,Zahn:1975}. These papers also introduce complementary
approaches to solving the tidal equations, that we dub `modal
decomposition' (MD) and `direct solution' (DS). In MD the tidal
perturbations are decomposed as a superposition of the star's
free-oscillation modes, with weights determined from overlap integrals
between the mode eigenfunctions and the tidal force field. In DS the
two-point boundary value problem (BVP) posed by the tidal equations is
solved directly using a standard approach such as shooting or
relaxation. Examples of MD are given by \cite{Kumar:1995},
\citet{Lai:1997}, \citet{Schenk:2001}, \citet{Arras:2003} and \citet{Fuller:2012}; and of DS by
\citet{Savonije:1983,Savonije:1984}, \citet{Pfahl:2008}, and
\citet{Valsecchi:2013}. The study by \citet[][hereafter
  B12]{Burkart:2012} is noteworthy in that it adopts \emph{both}
approaches, although no direct comparison is made between them (a
lacuna that appears to extend into the wider literature).

The following section lays out the theoretical foundations for our
treatment of tides. Section~\ref{s:implement} describes the
modifications to GYRE\ to implement this formalism via a DS
methodology, and Section~\ref{s:calcs} presents illustrative calculations
focused on selected problems in the published
literature. Section~\ref{s:discuss} summarizes the paper, discusses
potential applications of the new GYRE\ functionality, and outlines
future improvements.


\section{Theoretical Formalism} \label{s:formalism}

Rather an exhaustive derivation of all equations, we opt to focus on
the key expressions that support our narrative and define the choices
(e.g., normalizations, sign conventions) dictated by the existing
numerical framework of GYRE. For more-detailed exposition,
we refer the reader to the papers by \citet{Polfliet:1990},
\citet{Smeyers:1991}, \citet{Smeyers:1998}, and
\citet{Willems:2003,Willems:2010}.

\subsection{Binary Configuration} \label{s:formalism-config}

Consider a binary system comprising a primary star of mass $M$
and photospheric radius $R$, together with a secondary star of
mass $q M$. To model the tides raised on the primary by the
secondary, we adopt a non-rotating reference frame with the primary's
center-of-mass fixed at the origin, the orbit of the secondary lying
in the Cartesian $xy$-plane, and the line of apsides coinciding with
the $x$-axis. The position vector of the secondary at time $t$ is then
given by
\begin{equation} \label{e:vsec}
  \mathbf{r}_{2}(t) = \frac{a(1 - e^{2})}{1 + e \cos\upsilon} \left( \mathbf{e}_{x} \cos
  \upsilon + \mathbf{e}_{y} \sin \upsilon \right),
\end{equation}
where $\mathbf{e}_{x}$ and $\mathbf{e}_{y}$ are the unit vectors along the $x$ and $y$
axes, respectively, $a$ is the orbital semi-major axis, and $e$ the
eccentricity. The true anomaly $\upsilon$ is linked to $t$ via Kepler's
equation
\begin{equation} \label{e:kepler}
  \mathcal{E} - e \sin\mathcal{E} = \mathcal{M}
\end{equation}
and the auxiliary relations
\begin{gather}
  \mathcal{M} = \Omega_{\rm orb} (t - t_{0}), \\
  \tan \frac{\upsilon}{2} = \sqrt{\frac{1-e}{1+e}} \tan \frac{\mathcal{E}}{2},
\end{gather}
where $t_{0}$ is a time of periastron passage, and $\mathcal{M}$ and
$\mathcal{E}$ are the mean and eccentric anomalies, respectively.  The
orbital angular frequency $\Omega_{\rm orb}$ is given by Kepler's third law,
\begin{equation}
  G M (1 + q) = a^{3} \Omega_{\rm orb}^{2},
\end{equation}
with $G$ the gravitational constant.

\subsection{Tidal Potential} \label{s:formalism-pot}

Tides are raised on the primary star by the forces arising from the
gravitational potential $\Phi_{2}$ of the secondary, which at position
vector $\mathbf{r}$ and time $t$ is
\begin{equation}
  \Phi_{2}(\mathbf{r};t) = - \frac{q G M}{|\mathbf{r}_{2} - \mathbf{r}|}.
\end{equation}
Using a multipolar expansion in space and a Fourier-series expansion
in time, this expression can be recast as
\begin{multline} \label{e:pot-expand}
  \Phi_{2}(\mathbf{r};t) =
  - \frac{q G M}{r_{2}}
  - \frac{q G M}{r_{2}^{2}} r \sin\vartheta \cos(\varphi - \upsilon) + \mbox{} \\
  \pottide(\mathbf{r};t)
  \end{multline}
where $(r,\vartheta,\varphi)$ are the spherical-polar radius,
colatitude and azimuth coordinates corresponding to $\mathbf{r}$. The first
term on the right-hand side is constant and therefore does not
generate a force. The second term produces a spatially uniform force
directed from the primary star toward the secondary, and precisely
cancels the inertial force arising from the orbital motion of the
primary center-of-mass about the system center-of-mass. The third term
represents the tidal part of the secondary potential, and is expressed
as a superposition
\begin{equation} \label{e:pottide}
  \pottide(\mathbf{r};t) =
  \sum_{\ell=2}^{\infty} \sum_{m=-\ell}^{\ell} \sum_{k=-\infty}^{\infty}
  \pottide[\ell,m,k] (\mathbf{r};t)
\end{equation}
of partial tidal potentials defined by
\begin{multline} \label{e:pottide-part}
  \pottide[\ell,m,k](\mathbf{r};t) =
  - \varepsilon_{\rm T} \,
  \frac{G M}{R} \,
  \cbar[\ell,m,k]
  \left( \frac{r}{R} \right)^{\ell} \times \mbox{} \\
  Y^{m}_{\ell}(\vartheta,\varphi) \,
  \ee{- \mathrm{i} k \mathcal{M}}.
\end{multline}
Here,
\begin{equation}
  \varepsilon_{\rm T} \equiv
  \left( \frac{R}{a} \right)^{3} \, q =
  \frac{\Omega_{\rm orb}^{2} R^{3}}{G M} \, \frac{q}{1 + q}
\end{equation}
is a dimensionless parameter that quantifies the overall strength of
the tidal forcing, $\cbar[\ell,m,k]$ is an tidal expansion coefficient
(Appendix~\ref{a:exp-coeff}), and $Y^{m}_{\ell}$ is a spherical
harmonic (Appendix~\ref{a:sph-harm}).

\subsection{Tidal Response} \label{s:formalism-resp}

Appendix~\ref{a:lin-eqns} summarizes the set of linearized equations
governing the response of the primary star to the tidal potential
$\pottide$. Based on the form~(\ref{e:pottide-part}) of the partial
potentials, solutions to these equations take the form
\begin{multline} \label{e:sol-xi}
  \boldsymbol{\xi}(\mathbf{r};t) =
  \sum_{\ell,m,k}
  \left[
    \txir[\ell,m,k](r) \, \mathbf{e}_{r} + \phantom{\left( \mathbf{e}_{\vartheta} \, \pderiv{}{\vartheta} \right)}
  \right. \\
  \left.
    \txih[\ell,m,k](r) \left( \mathbf{e}_{\vartheta} \, \pderiv{}{\vartheta} + \frac{\mathbf{e}_{\varphi}}{\sin\vartheta} \, \pderiv{}{\varphi}\right)
  \right] \,
  \angtime[\ell,m,k](\vartheta,\varphi;t)
\end{multline}
for the displacement perturbation vector $\boldsymbol{\xi}$, and
\begin{equation} \label{e:sol-f}
  f'(\mathbf{r};t) =
  \sum_{\ell,m,k}
  \tf[\ell,m,k]'(r) \, \angtime(\vartheta,\varphi;t)
\end{equation}
for the Eulerian ($f'$) perturbation to a scalar variable $f$ (the
corresponding Lagrangian perturbation $\delta f$ follows from
equation~\ref{e:lag-eul}). In these expressions, the notation
$\sum_{\ell,m,k}$ abbreviates the triple sum appearing in
equation~(\ref{e:pottide}), while $\mathbf{e}_{r}$, $\mathbf{e}_{\vartheta}$ and $\mathbf{e}_{\varphi}$ are the
unit basis vectors in the radial, polar, and azimuthal directions,
respectively. The functions
\begin{equation} \label{e:angtime}
  \angtime[\ell,m,k](\vartheta,\varphi;t) \equiv Y^{m}_{\ell}(\vartheta,\varphi) \ee{-\mathrm{i} k \mathcal{M}}
\end{equation}
describe the angular and time dependence of the response, while the
functions with tilde accents ($\txir,\txih,\tf'$) encapsulate the
radial dependence. The latter are found as solutions to a system of
tidal equations summarized in
Appendix~\ref{a:tidal-eqns}. Importantly, the set of radial functions
for a given combination of indices $\{\ell,m,k\}$ can be determined
independently of any other combination.

\section{Implementation in GYRE} \label{s:implement}

To implement the tidal equations~(\ref{e:tidal-mass}--\ref{e:tidal-diff}) in GYRE, which follows a DS methodology, we first
transform to a dimensionless independent variable $x \equiv r/R$
and a set of dimensionless dependent variables
\begin{equation}
  \begin{aligned}
    y_{1} &\equiv x^{2 - \ell} \frac{\txir}{r}, \\
    y_{2} &\equiv x^{2 - \ell} \frac{\tP'}{\rho g r}, \\
    y_{3} &\equiv x^{2 - \ell} \frac{\tpottot'}{gr}, \\
    y_{4} &\equiv x^{2 - \ell} \frac{1}{g} \deriv{\tpottot'}{r}, \\
    y_{5} &\equiv x^{2 - \ell} \frac{\delta \tS}{c_{P}}, \\
    y_{6} &\equiv x^{-1 - \ell} \frac{\delta \tLrad}{L}.
  \end{aligned}
\end{equation}
Here $g \equiv \sderiv{\Phi}{r}$ is the scalar gravity, and the other symbols are defined in Appendix~\ref{a:lin-eqns} (for notational simplicity we neglect the $\ell,m,k$ subscripts from perturbed quantities). With these transformations, the differential equations and boundary conditions
governing $y_{1},\ldots,y_{6}$ can be written in a form almost
identical to the linear, non-radial, non-adiabatic free oscillation equations detailed in
Appendix B2 of \citet{Townsend:2018}. The only differences are that
the interpretation of the $y_{3}$ and $y_{4}$ variables is altered;
the outer boundary condition governing the gravitational potential
acquires an inhomogeneous term on the right-hand side
\begin{equation} \label{e:inhom-bc}
 U y_{1} +  (\ell + 1) y_{3} + y_{4} = (2 \ell + 1) y_{\rm T},
\end{equation}
where $U$ is the usual homology invariant \citep[e.g.,][]{Kippenhahn:2013} and
\begin{equation} \label{e:y-T}
  y_{\rm T} \equiv x^{2 - \ell} \frac{\tpottide[\ell,m,k]}{gr};
\end{equation}
and the dimensionless oscillation frequency is defined by
\begin{equation} \label{e:omega}
  \omega \equiv \sqrt{\frac{R^{3}}{G M}} \sigma_{m,k},
\end{equation}
where $\sigma_{m,k}$ is the rotating-frame frequency defined in
equation~(\ref{e:sigmk}).

GYRE\ uses a multiple shooting algorithm to discretize BVPs on a
radial grid $x = x_{1}, \ldots, x_{N}$. For a system of of $n$
differential equations and $n$ boundary conditions (in the present
case, $n=6$) this leads to a corresponding system of linear algebraic
equations with the form
\begin{equation} \label{e:lin-sys}
  \mathsf{\mathbf{S}} \mathbf{u} = \mathbf{b},
\end{equation}
where $\mathsf{\mathbf{S}} \in \mathbb{C}^{N n \times N n}$ is a block staircase
matrix and the solution vector $\mathbf{u} \in \mathbb{C}^{N n}$ contains the
dependent variables $y_{1},\ldots,y_{n}$ evaluated at successive grid
points. For free-oscillation problems the right-hand side vector $\mathbf{b}
\in \mathbb{C}^{N n}$ is identically zero, and so the linear
system~(\ref{e:lin-sys}) is homogeneous; non-trivial solutions exist
only when the determinant of $\mathsf{\mathbf{S}}$ vanishes, a condition that
ultimately determines the dimensionless eigenfrequencies $\omega$ of
the star. However, for the tidal problem considered here, $\mathbf{b}$ is
non-zero due to the inhomogeneous term appearing on the right-hand
side of the boundary condition~(\ref{e:inhom-bc}). The linear system
can then be solved for any choice of $\omega$ --- that is, the
dimensionless frequency is an \emph{input} rather than an output.

The reordered flow of execution ($\omega$ as input rather than output)
motivates the decision to provide two separate executables in release
7.0 of GYRE: \texttt{gyre}\ for modeling free stellar oscillations (the same
as in previous releases), and \texttt{gyre\_tides}\ for modeling stellar tides
(new to this release). A full description of these programs, including
their input parameters and output data, is provided on the
GYRE\ documentation
site\footnote{\url{https://gyre.readthedocs.io/en/stable/}}


\section{Example Calculations} \label{s:calcs}

\subsection{Surface Perturbations in a KOI-54 Model} \label{s:calcs-surf}

\begin{table}
  \begin{centering}
  \begin{tabular}{cccccc}
    Model  & $M$ (\si{M_{\odot}}) & $R$ (\si{R_{\odot}}) & $T_{\rm eff}$ (\si{\kelvin}) & $Z_{\rm i}$ & $X_{\rm c}$ \\ \hline
    KOI-54 & 2.32 & 2.19 & \num{9400}  & 0.0328 & 0.487 \\
    B-star & 5.00 & 2.80 & \num{16700} & 0.0200 & 0.658 \\
  \end{tabular}
  \caption{Fundamental parameters for the two stellar models
    discussed in the text: mass $M$, radius $R$, effective
    temperature $T_{\rm eff}$, initial metal mass fraction $Z_{\rm i}$ and center hydrogen mass fraction $X_{\rm c}$.} \label{t:models}
  \end{centering}
\end{table}

\begin{figure}
  \includegraphics{fig_prop_diag_koi54}
  \caption{Propagation diagram for the KOI-54 model discussed in the
    text, plotting the Brunt-V\"ais\"al\"a\ ($N$) and Lamb ($S_{\ell}$) frequencies as a
    function of pressure coordinate for $\ell=2$ modes. The hatched
    regions indicate the acoustic wave ($\sigma > N,S_{\ell}$) and
    internal gravity wave q($\sigma < N,S_{\ell}$) propagation
    regions. The sloped dashed line shows the frequency corresponding
to the local thermal timescale $t_{\rm thm}$ (equation~\ref{e:t-thm}), while the upper (lower) horizontal dashed line corresponds to a frequency $\sigma/\Omega_{\rm orb} = 30$ (15).}
  \label{f:prop-diag-koi54}
\end{figure}

KOI-54 is a highly-eccentric ($e \approx 0.8$), near face-on binary
system comprising a pair of A-type stars. Discovered by
\citet{Welsh:2011} in Kepler observations, it serves as the archetype
of the heartbeat class of periodic variables
\citep{Thompson:2012}. B12 and \citet{Fuller:2012} each present
initial attempts to model the light curve of KOI-54 in terms of
contributions from the equilibrium tide, dynamical tides and stellar
irradiation. Here we undertake a calculation to reproduce Fig.~6 of
B12, which illustrates how a model for the KOI-54 primary responds to
forcing by a single partial tidal potential.

Guided by the parameters given in Table~1 of B12, we use release
r22.11.1 of the MESA\ software instrument
\citep{Paxton:2011,Paxton:2013,Paxton:2015,Paxton:2018,Paxton:2019,Jermyn:2022}
to evolve a \SI{2.32}{M_{\odot}} model from the zero-age main sequence
(ZAMS) until its photospheric radius has expanded to
\SI{2.19}{R_{\odot}}. Input files for this and subsequent
MESA\ calculations are available through Zenodo at
\dataset[10.5281/zenodo.7489814]{\doi{10.5281/zenodo.7489814}}. The
growth of the core is followed using the convective premixing
algorithm \citep{Paxton:2019}, but rotation is neglected. We avoid any
smoothing of the Brunt-V\"ais\"al\"a\ frequency profile, as this can introduce small
but consequential departures from mass conservation. Fundamental
parameters for the model are summarized in Table~\ref{t:models}, and
its propagation diagram is plotted in Fig.~\ref{f:prop-diag-koi54}.

\begin{figure*}
  \includegraphics{fig_scan_sigma}
  \caption{Plot of the radial displacement perturbation $\txir$ and
    the Lagrangian radial flux perturbation $\delta \tFradr$ at the
    surface of the KOI-54 model, as a function of forcing frequency
    $\sigma$. Selected peaks are labeled by the resonant mode's
    classification within the Eckart-Osacki-Scuflaire scheme.}
  \label{f:scan-sigma}
\end{figure*}

Figure 6 of B12 plots the modulus of $\txir$ and $\delta \tFradr$ as a
function of $\sigma/\Omega_{\rm orb}$, for a single partial potential with
$\ell=2$. The choices of $m$ and $k$ are left undetermined because B12
neglect rotation in their figure, and treat $\sigma$ as a free
parameter rather than being constrained by
equation~(\ref{e:sigmk}). They normalize the strength of the partial
potential in a manner equivalent to setting $y_{\rm T} = 1$ in
equation~(\ref{e:inhom-bc}). We use \texttt{gyre\_tides}\ to repeat these steps
for our KOI-54 model, plotting the results in Fig.~\ref{f:scan-sigma}.

This figure shows qualitative agreement with Fig.~6 of B12. For
$\sigma/\Omega_{\rm orb} \gtrsim 30$ the surface perturbations exhibit distinct
peaks, corresponding to resonances with the star's $\ell=2$
free-oscillation modes. Selected peaks are labeled with the resonant
mode's classification within the Eckart-Osaki-Scuflaire scheme
\citep[e.g.,][]{Unno:1989}. For $\sigma/\Omega_{\rm orb} \lesssim 30$ the peaks
merge together and dissolve, because the periods of the resonant modes
become appreciably shorter than the local thermal timescale
\begin{equation} \label{e:t-thm}
  t_{\rm thm}(r) \equiv \int_{r}^{R} \frac{4 \pi r^{2} \rho T c_{P}}{L} \diff{r}
\end{equation}
in the outer parts of the main mode-trapping cavity (see
Fig.~\ref{f:prop-diag-koi54}), resulting in significant
radiative-diffusion damping that broadens and suppresses the
resonances. Eventually, for $\sigma/\Omega_{\rm orb} \lesssim 15$ the surface
perturbations reach the limits $\txir/R \rightarrow 1$ and $\delta
\tFradr/F_{{\rm rad},r} \rightarrow 4$ corresponding to the $\ell=2$
equilibrium tide (see Section~6.2 of B12 for a discussion of these
limits).

On closer inspection, some differences between the two figures are
apparent. While the KOI-54 stellar models in B12 and the present work
are similar (in particular, sharing the same $M$ and $R$),
they are not identical; therefore, the locations and heights of the
resonance peaks are not the same in each figure. More significantly,
over the range $5 \lesssim \sigma/\Omega_{\rm orb} \lesssim 30$ the flux
perturbation behaves much more smoothly in Fig.~\ref{f:scan-sigma}
than in Fig.~6 of B12. The reason for this discrepancy is not obvious,
but we speculate that it may be linked to differences in the
near-surface convection of the models. The propagation diagram for our
KOI-54 model (Fig.~\ref{f:prop-diag-koi54}) shows a He\textsc{ii}
convection zone at $\log(P/\si{\dyne\per\cm\squared}) \approx 5.5$,
but this zone is absent in B12's model (cf. their Fig.~1).


\subsection{Direct Solution versus Mode Decomposition in a KOI-54 Model} \label{s:calcs-ds-vs-md}

\begin{figure*}
  \includegraphics{fig_scan_sigma_zoom}
  \caption{Zoom-in on the four labeled resonance peaks in
    Fig.~\ref{f:scan-sigma}, plotting the complex amplitude (upper
    panels; rescaled to have a maximum value of unity) and phase
    (lower panels) of $\txir/R$ as a function of normalized detuning
    parameter $\epsilon$ (equation~\ref{e:detune}). Separate curves are
    shown for the DS (\texttt{gyre\_tides}) and MD approaches.}
  \label{f:scan-sigma-zoom}
\end{figure*}

\begin{figure}
  \includegraphics{fig_eigfuncs}
  \caption{Real (upper panel) and imaginary (lower panel) parts of the
    wavefunctions associated with the radial displacement perturbation
    $\txir$ at the peak of the $\gmode{35}$ resonance, plotted as a
    function of pressure coordinate. The separate curves correspond to
    the DS (\texttt{gyre\_tides}) and MD approaches.}
  \label{f:eigfuncs}
\end{figure}

As a validation of the results presented in Fig.~\ref{f:scan-sigma},
we re-calculate the surface perturbations using the MD approach. We
follow the formalism laid out in Section 3.2 of B12, although adopting
non-adiabatic eigenfrequencies $\sigma_{n,\ell}$ and damping rates $\gamma_{n,\ell}$
provided by \texttt{gyre}\ to evaluate the Lorentzian factor
$\Delta_{n,\ell,m,k}$ (their equation 13; here, $n$ is the mode radial
order). To evaluate the overlap integrals $Q_{n,\ell}$ that weight the
contribution of each free-oscillation mode in the MD superposition, we
use the third expression of equation~(9) in B12; we find that the the
first expression yields unreliable values when $|n| \gtrsim 20$,
because the integrand is highly oscillatory and suffers from
significant cancellation.

Figure~\ref{f:scan-sigma-zoom} zooms in on the four labeled
resonance peaks from Fig.~\ref{f:scan-sigma}, plotting the
complex amplitude and phase of $\txir$ as a function of normalized
detuning parameter
\begin{equation} \label{e:detune}
  \bar{\epsilon} \equiv \frac{\sigma - \sigma_{n,\ell}}{|\gamma_{n,\ell}|},
\end{equation}
for the two approaches.  DS (i.e., \texttt{gyre\_tides}) and MD agree at higher
forcing frequencies (right-hand panels), but show mismatches toward
lower frequencies (left-hand panels).

To delve further into these discrepancies, Fig.~\ref{f:eigfuncs} plots
the $\txir$ wavefunction evaluated at the peak ($\bar{\detune}=0$) of the
$\gmode{35}$ resonance. The imaginary part of the wavefunction is
spatially oscillatory because it is dominated by the dynamical tide,
comprising the resonantly forced oscillation mode. The real part is
non-oscillatory and corresponds to the equilibrium tide, comprising
the superposition of the other, non-resonant $\ell=2$ modes.

The discrepancies in the wavefunction are restricted to the outer
parts of the stellar envelope. For $\operatorname{Re}(\txir)$, the DS and MD
curves begin to diverge at $\log (P/\si{\dyne\per\cm\squared}) \approx
8$ (corresponding to $r/R \gtrsim 0.96$), while for $\operatorname{Im}(\txir)$
the divergence begins further out at $\log
(P/\si{\dyne\per\cm\squared}) \approx 6$ ($r/R \approx 0.99$). We
hypothesize that these divergences arise because $\sigma \ll
2\pi/t_{\rm thm}$ in these superficial layers, leading to significant
non-adiabaticity that MD is unable to correctly reproduce (see Section
6.2 of B12; also, Section 6 of \citealp{Fuller:2017}).


\subsection{Circularization in a B-star Model} \label{s:calcs-circ}

\begin{figure}
  \includegraphics{fig_prop_diag_5msun}
  \caption{Propagation diagram for the B-star model introduced in the
    text (cf. Fig.~\ref{f:prop-diag-koi54}). The dip in the
    Brunt-V\"ais\"al\"a\ frequency around $\log (P/\si{\dyne\per\cm\squared})
    \approx 7.5$ is caused by the iron opacity bump responsible for
    the overstability of the $\gmode{9}-\gmode{17}$ modes.}
  \label{f:prop-diag-5msun}
\end{figure}

\begin{figure*}
  \includegraphics{fig_secular}
  \caption{The circularization timescale $t_{e}$ plotted as a
    function of orbital period $P_{\rm orb}$, for the binary system with the
    B-star primary. Resonances are labeled beneath with the mode
    classification, and in parentheses the harmonic index $k$ and the
    sense of propagation in the co-rotating frame (p=prograde,
    r=retrograde). If the resonance is with an overstable mode, then
    the label border is bolded.}
  \label{f:secular}
\end{figure*}

\begin{figure}
  \includegraphics{fig_eigfreqs}
  \caption{Damping rates $\gamma_{n,\ell}$ plotted against
    eigenfrequencies $\sigma_{n,\ell}$ for $\gmode{1}$--$\gmode{20}$ modes of the
    B-star model, as calculated using \texttt{gyre}\ and as tabulated
    by W03. Open symbols indicate stable ($\gamma_{n,\ell} > 0$) modes, and
    filled symbols overstable ($\gamma_{n,\ell} < 0$) modes.}
  \label{f:eigfreqs}
\end{figure}

\citet[][hereafter W03]{Willems:2003} explore the secular orbital
changes due to tides in a binary system comprising a \SI{5}{M_{\odot}}
B-star primary and a \SI{1.4}{M_{\odot}} secondary on an $e=0.5$
orbit. They adopt the MD approach, but include only a single term at a
time in the modal superposition. The calculations illustrated in
Fig.~3 of W03 are used by \citet{Valsecchi:2013} to benchmark their
CAFein\ code, motivating us to do likewise here with \texttt{gyre\_tides}.

We use MESA\ as before to evolve a \SI{5}{M_{\odot}} model from the ZAMS
until its photospheric radius has grown to match the \SI{2.80}{R_{\odot}}
of W03's model. Fundamental parameters for this model are
summarized in Table~\ref{t:models}, and its propagation diagram is
plotted in Fig.~\ref{f:prop-diag-5msun}. Then, we apply
\texttt{gyre\_tides}\ to evaluate the star's response to the dominant
contributions in the tidal potential~(\ref{e:pottide}), comprising the
$\ell=2$ terms with $|m| = 2$ and $m k > 0$. We further restrict the
summation over $k$ to terms with a magnitude at least $10^{-12}$ times
that of the largest term.  As in W03, we assume the stellar angular
rotation frequency is equal to the periastron angular velocity of the
secondary,
\begin{equation}
  \Omega_{\rm peri} = \Omega_{\rm orb} \, \sqrt{\frac{(1+e)}{(1-e)^{3}}}.
\end{equation}
Based on this configuration, Fig.~\ref{f:secular} plots the circularization
timescale
\begin{equation}
  t_{e} \equiv - \left[ \frac{1}{e} \left( \deriv{e}{t} \right) \right]^{-1}_{\rm sec}
\end{equation}
as a function of orbital period $P_{\rm orb} \equiv 2\pi/\Omega_{\rm orb}$ over a pair
of intervals\footnote{These correspond to the short- and long-period
limits in Fig.~3 of W03.}. The secular rate-of-change of eccentricity
is evaluated via
\begin{multline} \label{e:secular-edot}
  \left( \deriv{e}{t} \right)_{\rm sec} = 4 \,
  \Omega_{\rm orb} \,
  q
  \sum_{\ell,m,k \geq 0}
  \left( \frac{R}{a} \right)^{\ell+3} \,
  \left( \frac{r_{\rm s}}{R} \right)^{\ell+1} \\
  \mbox{} \times \kap[\ell,m,k] \,
  \operatorname{Im}(\Fbar[\ell,m,k]) \,
  \Gbar[\ell,m,k]{3};
\end{multline}
this comes from equation~(55) of \citet{Willems:2010}, with
$F_{\ell,m,k}$ replaced by $\Fbar[\ell,m,k] \equiv F_{\ell,m,-k}$ and
$G^{(3)}_{\ell,m,k}$ by $\Gbar[\ell,m,k]{3} \equiv G^{(3)}_{\ell,m,-k}$ to account for
the differing sign convention in the assumed time dependence of
partial tides (see equation~\ref{e:angtime}). Note that the summation
is now restricted to $k \geq 0$.

The $|t_{e}|$ data plotted in the figure vary smoothly with $P_{\rm orb}$
between a series of sharp extrema. The maxima (more correctly,
singularities) arise when $(\sderiv{e}{t})_{\rm sec}$ passes through
zero. The minima arise from resonances with the star's $\ell=2$
free-oscillation modes, similar to the peaks seen in
Fig.~\ref{f:scan-sigma}. However, a key difference here is that the
star is being forced with a superposition~(\ref{e:pottide}) of partial
tidal potentials, rather than a single one as before. The criterion
for resonance $\sigma_{m,k} \approx \sigma_{n,\ell}$ can be satisfied for many
different values of $k$, leading to multiple resonances with the same
mode. This can be seen in the figure; for instance, the upper panel
shows resonances between the $f$ mode and the partial tidal potentials with
$k=31,\ldots,35$.

In the vicinity of some of the resonances shown in the figure, $t_{e}
< 0$ (blue) indicates that the tide acts to increase the eccentricity
of the orbit, driving it further away from circular. This behavior is
an instance of the `inverse tides' phenomenon discussed by
\citet{Fuller:2021}, and occurs when the summation in
equation~(\ref{e:secular-edot}) is dominated by a single,
\emph{positive} term. There are four distinct configurations that lead
to this outcome:
\begin{enumerate}[label=\Roman*.]
\item $\Gbar[\ell,m,k]{3} > 0$ and $\operatorname{Im}(\Fbar[\ell,m,k]) > 0$, the latter because
  \begin{enumerate}
  \item the resonant mode is prograde in the co-rotating frame ($\sigma_{n,\ell}/m > 0$)
    and stable ($\gamma_{n,\ell} > 0$); or \label{l:pos-pro-s}
  \item the resonant mode is retrograde in the co-rotating frame
    ($\sigma_{n,\ell}/m < 0$) and overstable ($\gamma_{n,\ell} < 0$). \label{l:pos-ret-o}
  \end{enumerate}
\item $\Gbar[\ell,m,k]{3} < 0$ and $\operatorname{Im}(\Fbar[\ell,m,k]) < 0$, the latter because
  \begin{enumerate}
  \item the resonant mode is prograde in the co-rotating frame ($\sigma_{n,\ell}/m > 0$)
    and overstable ($\gamma_{n,\ell} < 0$); or \label{l:neg-pro-o}
  \item the resonant mode is retrograde in the co-rotating frame
    ($\sigma_{n,\ell}/m < 0$) and stable ($\gamma_{n,\ell} > 0$). \label{l:neg-ret-s}

  \end{enumerate}
\end{enumerate}
All of the $t_{e} < 0$ resonances seen in Fig.~\ref{f:secular} are
instances of cases~(\ref{l:pos-ret-o}) or~(\ref{l:neg-pro-o}), and
therefore involve overstable modes. The overstability is caused by the
iron-bump opacity mechanism responsible for the slowly pulsating B
(SPB) stars \citep[e.g.,][]{Gautschy:1993,Dziembowski:1993}; in the
B-star model considered here, which falls well inside the SPB
instability strip \citep[e.g.,][]{Pamyatnykh:1999,Paxton:2015}, this
mechanism excites the $\ell=2$ $\gmode{9}$--$\gmode{17}$ modes.

Comparing Fig.~\ref{f:secular} against Fig.~3 of W03 reveals some
important differences. The latter shows numerous gaps and
discontinuities in $t_{e}$, that appear to arise because W03 only
allow a given mode to contribute toward the MD superposition when its
detuning parameter (equation~\ref{e:detune}) satisfies $0.1\,\varepsilon_{\rm T}
\leq |\bar{\detune} \gamma_{n,\ell}/ \sigma_{n,\ell}| \leq 10\,\varepsilon_{\rm T}$. The lower bound on
$\bar{\detune}$ means that the central parts of each resonance are omitted,
and so Fig.~3 of W03 does not fully reveal how small $t_{e}$ can
become when close to a resonance.

At the short-period limit of the range shown in the figures, there is
also disagreement between the typical magnitude of $t_{e}$ between
the resonances; Fig.~\ref{f:secular} shows an inter-resonance
$|t_{e}| \approx \SI{e9}{\year}$ , whereas for W03 it is 2--4 orders
of magnitude shorter. This is likely a consequence of W03
over-estimating the damping for the $\gmode{1}$--$\gmode{4}$ modes, which
dominate the tidal response at short orbital
periods. Figure~\ref{f:eigfreqs} plots eigenfrequencies $\sigma_{n,\ell}$ and
damping rates $\gamma_{n,\ell}$ for the $\gmode{1}$--$\gmode{20}$ modes, as calculated
using \texttt{gyre}\ and as obtained from Table~1 of W03. The data for the
$\gmode{5}$--$\gmode{16}$ modes are in reasonable agreement, especially given
that the underlying stellar models do not have the exact same internal
structure. However, the W03 damping rates for the $\gmode{1}$--$\gmode{4}$ modes
are four-to-five orders of magnitude larger than the \texttt{gyre}\ ones.


\subsection{Pseudo-Synchronization in a KOI-54 Model} \label{s:calcs-sync}

\begin{figure}
  \includegraphics{fig_torque}
  \caption{Secular torque $\torque_{\rm sec}$ plotted against stellar
    angular rotation frequency $\Omega_{\rm rot}$ for the KOI-54 primary
    model. The dashed line shows the equilibrium-tide torque extracted
    from Fig. 4 of B12.}
  \label{f:torque}
\end{figure}

B12 explore how tides can modify the rotation of the primary star in
the KOI-54 system, by evaluating the secular tidal torque $\torque_{\rm sec}$ on
the primary star as a function of the star's rotation rate. Using an
MD approach (see their Appendix C), they consider contributions toward
the torque from the $\ell=2$ partial tidal potentials. Their Fig. 4
shows a smoothly varying $\torque_{\rm sec}$ punctuated by many narrow peaks due
to modal resonances. The smooth part corresponds to the torque from
the equilibrium tide, and passes through zero at the
pseudo-synchronous angular frequency
\begin{equation}
  \Omega_{\rm ps} \equiv \Omega_{\rm orb} \frac{1 + (15/2)e^{2} + (45/8)e^{4} + (5/16)e^{6}}{[1 + 3e^{2} + (3/8)e^4](1 - e)^{3/2}}
\end{equation}
first derived by \citet{Hut:1981} in his theory of tides in the weak
friction limit. For the orbital parameters of KOI-54, the
pseudo-synchronous frequency is $\Omega_{\rm ps}/(2\pi) = \SI{0.395}{\per\day}$.

To repeat this calculation, we use \texttt{gyre\_tides}\ to evaluate the
response of the KOI-54 model (Section~\ref{s:calcs-surf}) to the
$\ell=2$ terms in the tidal potential~(Equation~\ref{e:pottide}) for
\num{25000} rotation frequencies in the interval $\SI{0}{\per\day}
\leq \Omega_{\rm rot}/(2\pi) \leq \SI{1}{\per\day}$. As in the preceding section,
we restrict the summation over $k$ to terms with a magnitude at least
$10^{-12}$ times that of the largest term. Then, we evaluate the tidal
torque via
\begin{multline}
  \torque_{\rm sec} = 
  4 \Omega_{\rm orb} \sqrt{\frac{G M^{3} q^{2}}{1 + q}} \, q \, a^{1/2}
  \sum_{\ell,m,k \geq 0}
  \left( \frac{R}{a} \right)^{\ell+3} \,
  \left( \frac{r_{\rm s}}{R} \right)^{\ell+1} \\
  \times \kap[\ell,m,k] \,
  \operatorname{Im}(\Fbar[\ell,m,k]) \,
  \Gbar[\ell,m,k]{4};
\end{multline}
this comes from equation~(63) of \citet{Willems:2010}, with
$G^{(4)}_{\ell,m,k}$ replaced by $\Gbar[\ell,m,k]{4} \equiv
G^{(4)}_{\ell,m,-k}$.

Fig.~\ref{f:torque} plots $\torque_{\rm sec}$ as a function of $\Omega_{\rm rot}$, using
discrete points rather than a continuous line because we are
significantly undersampling the dense forest of resonances. Also
plotted for comparison is the smooth (equilibrium tide) part of the
torque extracted from Fig.~4 of B12.  Clearly, there are some
significant differences between the two figures. Ours shows a torque
that's generally positive torque at small rotation frequencies, and
negative at high frequencies; however, the switch-over point is not
nearly as sharply defined as in B12, and occurs at a frequency
$\Omega_{\rm rot}/(2\pi) \approx \SI{0.5}{\per\day}$ around 25\% higher than
$\Omega_{\rm ps}$. Moreover, on either side of the switch-over, the lower
envelope of our torque values is around two orders of magnitude
smaller than the B12 curve.

Exploratory calculations indicate that these differences are not a
result of inaccurate overlap integrals (as was the case in
Section~\ref{s:calcs-ds-vs-md}), but rather due to a genuine
incompatibility between the DS and MD approaches. While a detailed
investigation of this problem is beyond the scope of the present
paper, we believe the fault lies with MD's over-estimation of
radiative dissipation for the equilibrium tide. If this hypothesis is correct, an immediate corollary is that pseudo-synchronization as envisaged by \citet{Hut:1981} does not operate for stars in which radiative dissipation dominates the tidal damping.


\section{Summary and Discussion} \label{s:discuss}

To briefly summarize the preceding sections: we establish the
theoretical foundations for our tides treatment
(Section~\ref{s:formalism}), describe modifications to GYRE\ to
implement tides (Section~\ref{s:implement}), and then apply the new
\texttt{gyre\_tides}\ executable to selected problems (Section~\ref{s:calcs}).
These example calculations uncover disagreements between the DS and MS
approaches, arising for a variety of reasons --- from numerical
inaccuracies in overlap integrals (Section~\ref{s:calcs-ds-vs-md}),
through to what appears to be an incompatibility between the
approaches. We look forward to future opportunities to investigate
these disagreements.

We also plan a number of enhancements to \texttt{gyre\_tides}\ that will extend
its capabilities. Key milestones include adding the ability to model
systems with spin-orbit misalignments \citep[e.g., following the
  formalism described by][]{Fuller:2017}; the incorporation of
additional damping mechanisms beyond radiative diffusion \citep[for
  instance, turbulent viscosity within convection zones; see
][]{Willems:2010}; and treating the effects of the Coriolis force,
which was neglected in deriving the linearized equations
(Appendix~\ref{a:lin-eqns}). Partial treatment of the Coriolis force
is already included in the main \texttt{gyre}\ executable, via the
traditional approximation of rotation (TAR; see, e.g.,
\citealp{Bildsten:1996,Lee:1997,Townsend:2003}). However, in its
current form this implementation cannot be used in \texttt{gyre\_tides},
because the angular operator appearing in the linearized continuity
equation within the TAR \citep[see equation~6 of][]{Bildsten:1996}
does not commute with the angular part of the Laplacian operator
appearing in the linearized Poisson equation~(\ref{e:lin-poisson}).

GYRE\ is not the first software package that adopts the DS approach
to model stellar tides; \citet{Pfahl:2008} and \citet{Valsecchi:2013}
describe functionally similar codes. Although the former authors' code
has never been made publicly available, the latters' CAFein\ code is
accessible on
GitHub\footnote{https://github.com/FrancescaV/CAFein}. After fixing a
number of bugs in CAFein\ (for instance, relating to misinterpreting
cell-centered quantities in MESA\ models as face-centered), we have
undertaken exploratory calculations comparing it against \texttt{gyre\_tides},
and find the two codes are in general agreement. Given that
CAFein\ is unmaintained, we decided that more-detailed comparison
would not be a worthwhile exercise.

It is our hope that \texttt{gyre\_tides}\ will provide a standardized and
well-supported community tool for simulating tides of spherical stars
within the linear limit. Specific areas where we anticipate productive
applications include modeling the many heartbeat systems discovered by
Kepler \citep[e.g.,][]{Thompson:2012}, TESS
\citep[e.g.,][]{Kolaczek:2021} and OGLE \citep{Wrona:2022};
investigating why most of these systems rotate faster than the
pseudo-synchronous rate \citep{Zimmerman:2017}; and exploring orbital
and rotational evolution in more-general star-star and star-planet
systems. These latter activities will initially be restricted to cases
where radiative diffusion dominates the tidal damping \citep[e.g., the
  $\gamma$ Doradus stars considered by ][]{LiGang:2020}; but with the
planned addition of other damping mechanisms, they can be extended
more broadly.


\section*{Acknowledgments}

This work has been supported by NSF grants ACI-1663696, AST-1716436
and PHY-1748958, and NASA grant 80NSSC20K0515. This research was also supported by STFC through grant ST/T00049X/1. The authors thank the referee for comments that have improved this paper.

\facilities{We have made extensive use of NASA's
  Astrophysics Data System Bibliographic Services.}

\software{Astropy \citep{astropy:2013,astropy:2018,astropy:2022},
  GYRE\ \citep{Townsend:2013,Townsend:2018,Goldstein:2020},
  Matplotlib \citep{Hunter:2007}, MESA\ 
  \citep{Paxton:2011,Paxton:2013,Paxton:2015,Paxton:2018,Paxton:2019,Jermyn:2022}}



