Heartbeat Stars, Tidally Excited Oscillations, and Resonance Locking

TL;DR

This paper develops advanced methods to analyze tidally excited oscillations in heartbeat stars, enhancing understanding of tidal interactions and resonance locking in eccentric binary systems.

Contribution

It introduces comprehensive computational techniques for modeling tidally excited oscillations, including effects of rotation, non-adiabaticity, and spin-orbit misalignment, with applications to resonance locking analysis.

Findings

Methods improve modeling of tidal oscillations in eccentric binaries.

Hot stars with surface temperatures above 6500 K show more visible oscillations.

Statistical framework helps distinguish chance resonances from resonance locking.

Abstract

Heartbeat stars are eccentric binary stars in short period orbits whose light curves are shaped by tidal distortion, reflection, and Doppler beaming. Some heartbeat stars exhibit tidally excited oscillations and present new opportunities for understanding the physics of tidal dissipation within stars. We present detailed methods to compute the forced amplitudes, frequencies, and phases of tidally excited oscillations in eccentric binary systems. Our methods i) factor out the equilibrium tide for easier comparison with observations, ii) account for rotation using the traditional approximation, iii) incorporate non-adiabatic effects to reliably compute surface luminosity perturbations, iv) allow for spin-orbit misalignment, and v) correctly sum over contributions from many oscillation modes. We also discuss why tidally excited oscillations are more visible in hot stars with surface temperatures $T \! \gtrsim \! 6500 \, {\rm K}$, and we derive some basic probability theory that can be used to compare models with data in a statistical manner. Application of this theory to heartbeat systems can be used to determine whether observed tidally excited oscillations can be explained by chance resonances with stellar oscillation modes, or whether a resonance locking process is operating.

Equations

$$\frac{\Delta L_N}{L}\simeq A_N\sin (N\Omega t+{\Delta }_N)$$

$$A_{N}=\epsilon_{l}X_{Nm}V_{lm}\big{|}Q_{\alpha}L_{\alpha}\big{|}\frac{\omega_{Nm}}{\sqrt{(\omega_{\alpha}-\omega_{Nm})^{2}+\gamma_{\alpha}^{2}}}\,.$$

$$\epsilon_{l}=\frac{M^{\prime}}{M}\bigg{(}\frac{R}{a}\bigg{)}^{l+1}\,.$$

$$Q_{\alpha}=\frac{\langle\xi_{\alpha}|{\mbox{\boldmath$\nabla$}}(r^{l}Y_{lm})\rangle}{\omega_{\alpha}^{2}}.$$

$$F_{Nm}=\frac{1}{\pi}\int^{\pi}_{0}d\Psi\frac{\cos\big{[}N\big{(}\Psi-e\sin\Psi\big{)}-m\phi(t)\big{]}}{\big{(}1-e\cos\Psi\big{)}^{l}}\,,$$

$${\omega }_{Nm}=N\Omega -m{\Omega }_s,$$

$${\Delta }_N={\delta }_{\alpha }-m{\varphi }_s$$

$$\delta_{\alpha}\simeq{\rm atan2}\bigg{[}\frac{(\omega_{\alpha m}-\omega_{Nm})}{-\gamma_{\alpha m}}\bigg{]}\,.$$

$$\delta \rho \sim \frac{\rho r^2}{g{H}^2}{\omega }_{\alpha }^2{\xi }_r,$$

$$\delta \rho \sim \frac{\rho r^2}{g{H}^{2}N}{\omega }_{\alpha }^3{\xi }_{\perp }.$$

$$\delta\Phi\sim\frac{-4\pi G}{2l+1}\frac{\rho r^{l+4}}{gHNR^{l+1}}\omega_{\alpha}^{3}\xi_{\perp}\bigg{|}_{r_{g}}\,.$$

$$\frac{\Delta F}{F}$$

$$\sim \sqrt{l(l+1)}\frac{{\xi }_{\perp }}{r}.$$

$$\gamma_{\alpha}\sim\frac{\int k_{r}^{2}K|{\mbox{\boldmath$\xi$}}|^{2}dr}{\int|{\mbox{\boldmath$\xi$}}|^{2}dr}\,,$$

$$\delta {\omega }_N=|{\omega }_{Nm}-{\omega }_{\alpha }|.$$

$$\delta {\omega }_{N,\min }=0$$

$$\delta {\omega }_{N,\max }=\Delta {\omega }_{\alpha }/2$$

$$\delta {\omega }_{N,\mathrm{mean}}=\Delta {\omega }_{\alpha }/4.$$

$${\mathcal{L}}_N={ϵ}_l{X}_{Nm}{V}_{lm}|Q_{\alpha }{L}_{\alpha }|,$$

$$\bigg{|}\frac{\Delta L_{N}}{L}\bigg{|}_{\rm med}\simeq\bigg{|}4\mathcal{L}_{N}\frac{\omega_{Nm}}{\Delta\omega_{\alpha}}\bigg{|}.$$

$$\bigg{|}\frac{\Delta L_{N}}{L}\bigg{|}_{\rm max}\simeq\bigg{|}\mathcal{L}_{N}\frac{\omega_{Nm}}{\gamma_{\alpha}}\bigg{|}.$$

$$\displaystyle{\rm P}\bigg{[}\frac{\Delta L_{N}}{L}>\frac{\Delta L_{c}}{L}\bigg{]}$$

$${\rm P}\bigg{[}\frac{\Delta L_{N}}{L}>\frac{\Delta L_{c}}{L}\bigg{]}={\rm P}\bigg{[}\delta\omega_{N}<\sqrt{Z_{\alpha}}\frac{\Delta\omega_{\alpha}}{2}\bigg{]}$$

$$Z_{\alpha}=\bigg{(}\frac{\Delta L_{c}}{L}\bigg{)}^{-2}\bigg{(}2\mathcal{L}_{N}\frac{\omega_{Nm}}{\Delta\omega_{\alpha}}\bigg{)}^{2}-\bigg{(}\frac{2\gamma_{\alpha}}{\Delta\omega_{\alpha}}\bigg{)}^{2}.$$

$$\displaystyle{\rm P}\bigg{[}\frac{\Delta L_{N}}{L}>\frac{\Delta L_{c}}{L}\bigg{]}$$

$$=0\mathrm{for}{Z}_{\alpha }<0$$

$$=1\mathrm{for}{Z}_{\alpha }>1.$$

$${\rm Num}_{\rm ex}\bigg{[}\frac{\Delta L_{N}}{L}>\frac{\Delta L_{c}}{L}\bigg{]}=\sum_{N_{\rm min}}^{N_{\rm max}}{\rm P}\bigg{[}\frac{\Delta L_{N}}{L}>\frac{\Delta L_{c}}{L}\bigg{]}.$$

$$\delta\omega_{N}=\sqrt{\omega_{Nm}^{2}\big{(}\mathcal{L}_{N}/A_{N}\big{)}^{2}-\gamma_{\alpha}^{2}}\,.$$

$$d\delta\omega_{N}=\frac{\omega_{Nm}\big{(}\mathcal{L}_{N}/A_{N}\big{)}^{2}}{\sqrt{\big{(}\mathcal{L}_{N}/A_{N}\big{)}^{2}-\big{(}\gamma_{\alpha}/\omega_{Nm}\big{)}^{2}}}\frac{dA_{N}}{A_{N}}\,.$$

$$\frac{dn}{d\delta {\omega }_N}=\frac{2}{\Delta {\omega }_{\alpha }}.$$

$$\frac{dn}{dA_{N}}=\frac{2\omega_{Nm}}{\Delta\omega_{\alpha}}\frac{\big{(}\mathcal{L}_{N}/A_{N}\big{)}^{2}}{\sqrt{\big{(}\mathcal{L}_{N}/A_{N}\big{)}^{2}-\big{(}\gamma_{\alpha}/\omega_{Nm}\big{)}^{2}}}\frac{1}{A_{N}}\,.$$

Full text

1 Introduction

Heartbeat stars are a growing class of eccentric ($e\gtrsim 0.3$) binary stars in short period ($1\,{\rm day}\lesssim P\lesssim 1\,{\rm yr}$) orbits, whose light curves are shaped by tidal distortion, reflection, and Doppler beaming. These effects are most prominent near periastron, combining to generate the characteristic “heartbeat” signature (whose shape resembles an EKG diagram) that is their namesake. Many heartbeat stars oscillate throughout their orbit due to the tidal excitation of stellar oscillation modes. The smoking gun signature of tidally excited oscillations (TEOs) is that they occur at exact integer multiples of the orbital frequency. Typical photometric variations of $\Delta L/L\lesssim 10^{-3}$ and characteristic time scales of days have hindered ground-based observations, and only a few heartbeat stars had been discovered in the pre-Kepler era (see De Cat et al. 2000; Willems & Aerts 2002; Handler et al. 2002; Maceroni et al. 2009).

Kepler has revolutionized studies of heartbeat stars by providing continuous, high precision photometry throughout multiple stellar orbits. The prototypical heartbeat star, KOI-54, has been examined in a series of papers (Welsh et al. 2011; Fuller & Lai 2012a; Burkart et al. 2012; O’Leary & Burkart 2014). It consists of two eccentric A-type stars in a highly eccentric ($e=0.83$), 42 day orbit. KOI-54 exhibits dozens of TEOs, the largest of which occur at exactly 90 and 91 times the orbital frequency. More recently, Thompson et al. (2012) presented 17 additional heartbeat stars with a variety of stellar components, and Beck et al. (2014) examined 15 sub-giant/red giant heartbeat systems. Additional heartbeat systems include those analyzed in Hambleton et al. (2013), Maceroni et al. (2014), and Hambleton et al. (2016). The radial velocity curves of many heartbeat systems have been measured in Smullen & Kobulnicky (2015); Shporer et al. (2016); Dimitrov et al. (2017), confirming that heartbeat stars are generally near the upper envelope of the binary eccentricity distribution. Currently, over 150 heartbeat stars have been identified in Kepler data Kirk et al. (2016), although few have been analyzed in detail.

Despite the fact that many heartbeat stars pulsate, their asteroseismic utility is currently somewhat limited. One reason is that TEOs pulsate at integer harmonics of the orbital frequency (rather than at stellar oscillation mode frequencies), with the largest amplitude TEOs occurring via resonances between orbital harmonics and mode frequencies. Although it may be possible to use this information to perform “tidal asteroseismology” (Burkart et al., 2012), this procedure is difficult. Nonetheless, TEOs differ from self-excited and stochastically excited pulsations because their amplitudes and phases can be straightforwardly predicted from linear theory and compared with observations. Moreover, TEOs allow for the observation of oscillation modes that are normally undetectable, allowing for new opportunities to constrain pulsation physics.

Heartbeat stars also offer an unprecedented opportunity to study tidal interactions between stars, and systems displaying TEOs are especially useful. Given the properties of the host star, the frequency and amplitude of a TEO can be used to identify the stellar eigenmode responsible for the oscillation and the energy contained in the pulsation. With an estimate of the mode damping rate, one can then calculate a tidal dissipation rate and orbital circularization time scale. Many heartbeat stars contain A-F type stellar components, which contain neither a thick convective envelope conducive to damping of the equilibrium tide, nor a large convective core conducive to excitation of the dynamical tide. Hence, A-F stars are not well described by commonly used tidal theories (e.g., Zahn 1977). Instead, tidal dissipation likely occurs through excitation and damping of TEOs, and so the information contained in heartbeat stars is of great importance.

Previous theoretical work on tidal excitation of stellar oscillations is extensive, and includes Zahn (1970, 1975, 1977); Goldreich & Nicholson (1989); Kumar et al. (1995); Lai (1997); Smeyers et al. (1998); Willems & Aerts (2002); Willems (2003); Willems et al. (2003), but very few of these works accurately calculate luminosity fluctuations produced by TEOs. One exception is Pfahl et al. (2008), which provides constructive insight on luminosity fluctuations in different types of stars. Unfortunately, the studies above could not compare observed TEOs with theoretical expectations, as TEOs were very difficult to observe in the pre-Kepler era. Nonetheless, these works contain the foundations of tidal theory upon which our work is built.

In this paper, we make detailed theoretical predictions for photometric amplitudes and phases of TEOs, and the tidal dissipation they produce. In particular, we discuss the relationship between the observed photometric amplitude of a TEO and the energy contained within the stellar pulsation, which is quite sensitive to the stellar structure. Since TEOs are typically due to resonances with gravity modes (g modes), they are usually only detectable in stars without thick convective envelopes (although see Fuller et al. 2013 for an exception). TEOs are also more likely to be observed (although not exclusively found) in high eccentricity systems, where tidal forcing occurs across a wide range of frequencies, and resonances with stellar oscillation modes are more probable. We construct some basic probability theory regarding expectations for the luminosity fluctuations produced by TEOs. Observed systems which defy theoretical expectations can then be used to improve theoretical understanding of tidal interactions in eccentric binaries. In two companion papers Hambleton et al. (2017); Fuller (2017), we compare these theories with data for KIC 8164262. We measure the parameters of the system and the frequencies and amplitudes of its TEOs, showing that most of them can be explained by chance resonances with g modes, with the exception of its highest amplitude pulsation, which can instead by explained by resonance locking.

Our paper is organized as follows. Section 2 discusses the basic theory of TEOs, and Section 3 discusses properties of g modes favorable to producing observable TEOs. Section 4 derives statistical properties of TEOs to be compared with observations. Resonance locking is discussed in Section 5. Tedious calculations required for precise predictions of TEO amplitudes and phases are provided in Sections 6, 7, 8, and 9. We conclude and provide discussion in Section 10.

2 Tidally Excited Oscillations

Heartbeat stars are distorted by the time variable tidal potential of the companion star. The response of the star can generally be decomposed into two components: the equilibrium tide and the dynamical tide. The equilibrium tidal distortion is simply the hydrostatic deformation of the star due to the companion, i.e., the steady state distortion that would be produced in the absence of orbital motion. The dynamical tide is the additional non-hydrostatic oscillation of the star that is produced due to the time variable nature of the tidal forcing. The equilibrium tide contributes to the “heartbeat” signature near periastron, while the dynamical tide is composed of TEOs that are visible at all orbital phases.

Figure 1 shows the phased lightcurve of a representative heartbeat star. We point out the periastron heartbeat signature (produced by the equilibrium tide, in addition to reflection and Doppler boosting), and a TEO produced by the dynamical tide. The Fourier transform reveals a characteristic comb of peaks at integer harmonics of the orbital frequency which are produced by the heartbeat distortion in the lightcurve, similar to the comb of peaks produced in the FFT of an eclipse light curve. Some peaks at orbital harmonics are much higher (or lower) amplitude than surrounding peaks, indicating that they are produced primarily by the dynamical tide. Finally, peaks that are not at orbital harmonics may be produced by heat-driven oscillation modes (see, e.g., Hambleton et al. 2013), non-linear tidal effects (Fuller & Lai 2012a; Burkart et al. 2012; Hambleton et al. 2013; Borkovits et al. 2014), or rotational variation due to spots.

2.1 Dynamical Tide

We calculate the dynamical tidal response of the star by decomposing it into the response of each stellar oscillation mode, indexed by $\alpha$. The total stellar response is found by summing over all modes. Each mode contributes to both the equilibrium tide and the dynamical tide, and thus has an equilibrium amplitude $a_{\rm eq}$ and a dynamical amplitude $a_{\rm dyn}$. Binary light curve modeling techniques generally solve for only the equilibrium tidal distortion, because they compute only the hydrostatic distortion produced by the companion. Therefore, after a binary lightcurve model is subtracted away, the amplitude of the remaining oscillations corresponds to dynamical mode amplitudes, $a_{\rm dyn}$. Our calculations below thus focus on the dynamical component of each mode. In Section 6, we formally define the equilibrium mode amplitude, describe how to compute it, and how to calculate the dynamical component of each oscillation mode.

In an eccentric binary, each stellar oscillation mode is forced at every integer harmonic $N$ of the orbital frequency. Consequently, the observable tidal response of a star in an eccentric binary consists of fluctuations at every harmonic of the orbital frequency, with each harmonic composed of a sum over all the star’s oscillation modes. In Section 7, we describe how to rigorously compute these sums, which must be done carefully to obtain accurate results.

2.2 Simple Cases

High amplitude TEOs are usually produced by a near-resonance between an orbital harmonic and a mode frequency. In this case, the luminosity fluctuation at an orbital harmonic $N$ can be approximated by a single term in the sum of equation 91. This approach is also described in Fuller & Lai (2012a); Burkart et al. (2012). Here, we provide general formulae including rotation and non-adiabatic effects. When a single resonant mode (indexed by $\alpha$) of azimuthal number $m$ dominates the dynamical tidal response at an orbital harmonic $N$, it produces a sinusoidal luminosity fluctuation of form

[TABLE]

with amplitude

[TABLE]

In equation 1, $\Omega$ is the angular orbital frequency and $t$ is time measured from $t=0$ at periastron. The dimensionless tidal forcing amplitude $\epsilon_{l}$ is determined by the mass and semi-major axis of the perturber (see equation 69),

[TABLE]

Here, $M^{\prime}$ is the mass of the companion, $M$ is mass of the primary, $R$ is its radius, and $a$ is the semi-major axis of the orbit. The $l$ subscript refers to multipole of the tidal potential which excites the mode, and typically the $l\!=\!2$ components are the most important for tidal excitation. The factor $V_{lm}=|Y_{lm}(i_{s},0)|$ (with $Y_{lm}$ a spherical harmonic) accounts for the visibility of the mode given an inclination $i_{s}$ between the star’s rotation axis and the line of sight.

$L_{\alpha}$ describes the observed luminosity fluctuation produced by a mode normalized via equation 52. The value of $L_{\alpha}$ is sensitive to the stellar model and oscillation mode frequency, and is plotted in Figure 2 for g modes in stars of different masses. $L_{\alpha}$ is small for g modes in stars with thick convective envelopes ($T_{\rm eff}\lesssim 6500\,{\rm K}$) because the modes are trapped below the surface by a thick convection zone. In the language of asteroseismology, these modes have very large inertias and are not easily excited. The g modes of hot stars ($T_{\rm eff}\gtrsim 6700\,{\rm K}$) propagate close to the surface and produce larger luminosity perturbations, hence TEOs are more visible in these stars. The small dips along each curve are created by modes trapped near the convective core, an effect partially washed out by adding small molecular diffusivity to a stellar model. The value of $L_{\alpha}$ for low frequency modes (${\rm f}_{\alpha}<0.8\,{\rm d}^{-1}$) in stars with $T_{\rm eff}\approx 7000\,{\rm K}$ is somewhat affected by the treatment of the convective flux perturbation (see Appendix A), and should be treated cautiously. The value of $L_{\alpha}$ is determined primarily by temperature perturbations produced by the mode near the stellar photosphere, thus, a non-adiabatic mode calculation is required for an accurate estimate of $L_{\alpha}$.

The dimensionless number $Q_{\alpha}$ describes the spatial coupling between the stellar oscillation mode $\alpha$ and the tidal potential. It is given by

[TABLE]

with all quantities in dimensionless units (i.e., mass in units of $M$, length in units of $R$, time in units of $\sqrt{R^{3}/GM}$). $Q_{\alpha}$ can also be expressed in terms of the the surface gravitational potential perturbation associated with each mode (see Section 9). Figure 2 shows values of $Q_{\alpha}$ for g modes in different stars. Typically $Q_{\alpha}$ is small, of order $10^{-6}\!-\!10^{-2}$ for gravito-inertial modes using our normalization condition (see Section 6).

The dimensionless number $X_{Nm}$ describes the strength of the tidal forcing at an orbital harmonic $N$ given an orbit of eccentricity $e$ and is defined in equation 105.​​111$X_{Nm}$ is similar to a Hansen coefficient, $X^{k}_{lm}$ (see Burkart et al. 2012), or a Fourier coefficient, $c_{l,m,k}$, Willems (2003). Circular and spin-aligned orbits have non-zero values of $X_{Nm}$ only when $N=m$. Misaligned orbits have non-zero values of $X_{Nm}$ for other values of $m$, while higher eccentricity orbits have stronger forcing at higher values of $N$. For aligned spin and orbit, $X_{Nm}=W_{lm}F_{Nm}$, where

[TABLE]

and $W_{lm}$ is a constant of order unity defined below equation 69. In equation 5, $\Psi$ is the eccentric anomaly, and $\phi(t)$ is the true anomaly. Figure 3 shows $F_{Nm}$ for different orbital eccentricities, for $m=0$ and $m=2$. $F_{Nm}$ peaks at larger orbital harmonics for higher eccentricity, and falls off exponentially at very large $N$, preventing tidal excitation of f modes and p modes. The value of $F_{Nm}$ is very small for $m=-2$, while its value is unimportant for $m=\pm 1$ because $W_{lm}=0$ for $l=2$, $m=\pm 1$.

The last term in equation 2 is the resonant detuning factor and accounts for the effects of resonant mode excitation. Here, $\omega_{\alpha}$ is the stellar oscillation mode frequency in the star’s rotating frame, $\gamma_{\alpha}$ is its growth rate, and $\omega_{Nm}$ is the rotating frame tidal forcing frequency,

[TABLE]

where $\Omega_{s}$ is the angular spin frequency. For weakly damped modes, the resonant detuning factor obtains very large values when $\omega_{\alpha}\simeq\omega_{Nm}$, i.e., near a resonance between the mode frequency and the tidal forcing frequency. We define modes to have time dependence $\propto e^{-i(\omega_{\alpha}+i\gamma_{\alpha})t}$ such that negative values of $\gamma_{\alpha}$ correspond to damped modes, and positive values of $m$ correspond to prograde modes.

Finally, for aligned spin and orbit, the phase of the luminosity fluctuation relative to periastron is

[TABLE]

where $\delta_{\alpha}$ is defined in equation 87, and $\phi_{s}$ is the azimuthal angle of the observer at periastron. In the misaligned case, the phase is given in equation 8.1. If the mode is nearly adiabatic, then the phase shift is

[TABLE]

If the mode is not extremely close to resonance, i.e., $(\omega_{\alpha}-\omega_{Nm})\gg\gamma_{\alpha}$, then $\delta_{\alpha}\simeq\pm\pi/2$. This is the regime examined in O’Leary & Burkart (2014), in which one may easily compare observed mode amplitudes/phases to theoretical expectations. Note that for a single mode to dominate the tidal response, we must be in the regime where the typical mode spacing, $\Delta\omega_{\alpha m}$, is larger than the typical mode damping rates $\gamma_{\alpha m}$. This criterion is violated for low frequency g-modes with large damping rates, i.e., modes in the traveling wave regime.

3 Gravity Mode Properties

Here we examine the g mode properties shown in Figure 2, and discuss the underlying physics and resulting scaling with mode frequency in different types of stars.

The scaling of $Q_{\alpha}$ with frequency can be understood as follows. The numerator of equation 4 is equal to an integral over the density perturabation throughout the star (see equation 9), which is dominated by the integral over the non-oscillatory evanescent tail of the g-mode (see also Goldreich & Wu 1999 and Goodman & Dickson 1998). For g modes in the WKB limit, the Eulerian density perturbation scales as (Luan et al., 2017)

[TABLE]

where $g$ is gravity and $H$ is the scale height. For WKB g modes, $\xi_{r}\sim(\omega_{\alpha}/N)\xi_{\perp}$, and we have

[TABLE]

The length scale of the evanescent tail is $\sim\!H$ and is independent of frequency. Hence, following Luan et al. (2017), the surface gravitational potential perturbation is

[TABLE]

This quantity should be evaluated at the radial location $r_{g}$ just inside the edge of the g mode propagation cavity, for the evanescent tail that dominates the coupling. For the mode frequencies and stars considered here, this is typically the outer edge of the g mode cavity where $\omega_{\alpha}=N$ or $\omega_{\alpha}=L_{l}$. For low frequency modes in more massive stars, the location of $r_{g}$ may be just outside the convective core, though we find this is only the case for low frequency modes in the $2.5\,M_{\odot}$ star in this study.

For our normalization condition, the magnitude of $\xi_{\perp}$ at a given radius is nearly independent of $\omega_{\alpha}$. If the mode’s turning point $r_{g}$ is independent of frequency, we find $\delta\Phi\propto\omega_{\alpha}^{3}$ and to good approximation $Q_{\alpha}\propto\omega_{\alpha}$ (a better scaling taking into account the evanescent wave function is $Q_{\alpha}\propto\omega^{17/6}$, see Weinberg et al. 2012, but note their different normalization convention), as can be seen for cool stars in Figure 2. In cool stars, the location of $r_{g}$ is the base of the convection zone, nearly independent of mode frequency. Figure 4 shows a propagation and g mode eigenfunction in a $1.4\,M_{\odot}$ where this is this case. In hot stars, the surface convection zone is much thinner, and for modes of $f_{\alpha}\gtrsim 0.6\,{\rm d}^{-1}$, the upper edge of the g mode cavity is determined by the location where $\omega_{\alpha}=L_{l}$ rather than the location of the surface convective zone. For this reason, the scaling of $Q_{\alpha}$ with frequency is steeper, because lower frequency modes have their outer turning point closer to the surface where the density is lower and the mode cannot couple effectively to the tidal potential.

The behavior of $L_{\alpha}$ in Figure 2 can also be understood at a basic level, and is discussed in detail of Pfahl et al. (2008). Following Appendix A2 of Luan et al. (2017), the Lagrangian perturbation to the energy flux for g modes in the WKB limit is of order

[TABLE]

with the second line following from the g mode dispersion relation. Equation 3 immediately shows $\Delta F/F\gg\xi_{r}/r$ for modes in the WKB limit, and temperature/flux perturbations are more important than surface area perturbations in creating the observed luminosity perturbation. To estimate the flux perturbation at the surface of the star, and hence the value of $L_{\alpha}$, equation 3 should be evaluated at the location where $\omega_{\alpha}t_{\rm therm}\sim 1$, where $t_{\rm therm}$ is the local thermal time, as discussed in Pfahl et al. (2008). At shallower depths, the flux perturbation is approximately constant. In cool stars, the value of $\Delta F/F$ is small where $\omega_{\alpha}t_{\rm therm}\sim 1$ because this location is in the convective region, multiple scale heights above the outer boundary of the g mode cavity, where the eigenfunction $\xi_{r}$ has become much smaller. In hot stars, for g modes with low frequency, this location can be within the g mode cavity, or just above it. Lower frequency modes have their turning point closer to the surface where the density is lower, and hence $\xi_{\perp}$ is larger. This creates larger values of $L_{\alpha}$ in hotter stars and for modes with lower frequency. For very low frequency modes ($f_{\alpha}\lesssim 0.4\,{\rm d}^{-1}$ in Figure 2, their amplitude is significantly attenuated by damping because they propagate where $\omega t_{\rm therm}\lesssim 1$, so the value of $L_{\alpha}$ decreases.

Finally the behavior of $\gamma_{\alpha}$ is relatively simple to understand. The mode damping rate is

[TABLE]

where $K$ is the thermal diffusivity and $k_{r}\approx\sqrt{l(l+1)}N/(\omega_{\alpha}r)$ is the g mode radial wavenumber. For g modes with a fixed propagation cavity, we obtain $\gamma_{\alpha}\propto\omega_{\alpha}^{-2}$ as seen for g modes in cool stars in Figure 2. In hot stars, lower frequency modes propagate closer to the surface where $K$ is larger. This creates larger mode damping rates and a steeper scaling with mode frequency. As mentioned above, very low frequency modes ($f_{\alpha}\lesssim 0.4\,{\rm d}^{-1}$) propagate into regions where they are strongly damped. In this case, they are essentially damped on a wave crossing time scale, such that $\gamma_{\alpha}\sim\int v_{g}dr\sim\Delta\omega_{\alpha}\propto\omega_{\alpha}^{2}$, and so the mode damping rate starts to decrease with decreasing frequency. Our GYRE calculations only consider radiative diffusion as a source of damping, and neglect interaction with convection, or damping due to convective viscosity. This is problematic for stars near the $\gamma$-Doradus instability strip whose instability is driven by interaction with convection (Guzik et al., 2000). Convective turbulent viscosity could potentially be important, but is likely sub-dominant to convective flux blocking effects as discussed in (Guzik et al., 2000).

We comment on the real vs. imaginary pieces of mode eigenfunctions. First, imaginary components of the mode eigenfunctions can be quite large (see Figure 4), comparable to the real components, even in somewhat cool stars where non-adiabatic effects are weaker than in hot stars. This generally causes $L_{\alpha}$ to have a significant imaginary component, and mode luminosity perturbations discussed in Section 7.1 can be significantly shifted in phase relative to the expectation for an adiabatic mode. Figures 4 and 5 show that $\xi_{\perp}^{I}=0$ at the star’s surface, but this is merely a consequence of the arbitrary choice of a phase convention $\xi_{r}^{I}=0$ at the surface. Since the outer boundary condition is $\delta P=\rho g\xi_{r}$ and $\xi_{\perp}\propto\delta P$, the boundary condition results in $\xi_{\perp}^{I}=0$ a the surface. This is not true for temperature and flux perturbations, which can have a substantial phase shift relative to the mode displacement.

4 Statistical Approach

In practice, accurately predicting amplitudes of TEOs is extremely difficult. The reason is that the mode amplitudes depend sensitively on the detuning factor, $\omega_{Nm}/\sqrt{(\omega_{\alpha}-\omega_{Nm})^{2}+\gamma_{\alpha}^{2}}\gg 1$ for nearly resonant modes. Very small changes in the stellar model that alter the mode frequencies $\omega_{\alpha}$ translate into very large changes in the detuning factor and predicted TEO amplitudes.

A useful approach is thus to examine tidal mode excitation from a statistical point of view. Given a stellar and orbital model, we may reliably calculate the values of the terms in equation 2 for modes with $\omega_{\alpha}\approx\omega_{Nm}$. We may then estimate the expected TEO amplitudes from the expected values of the mode detuning factor. Consider tidal forcing frequencies for modes of azimuthal number $m$. The forcing frequencies are uniformly distributed at intervals of $\Delta\omega_{Nm}=N\Omega$ in frequency. The stellar eigenfrequencies, in turn, are distributed with a typical spacing of $\Delta\omega_{\alpha}$ determined by the stellar structure. Typically these mode frequencies are distributed smoothly so that we may calculate a typical mode frequency separation $\Delta\omega_{\alpha}$ at frequencies where $\omega_{\alpha}\approx\omega_{Nm}$. In heartbeat star systems, we find that both regimes $\Delta\omega_{Nm}<\Delta\omega_{\alpha}$ and $\Delta\omega_{Nm}>\Delta\omega_{\alpha}$ may be realized.

For a forcing frequency $\omega_{Nm}$, we are interested only in the most resonant mode, i.e., the mode for which

[TABLE]

is minimized. The minimum possible frequency difference is

[TABLE]

while the maximum possible frequency difference between a forcing frequency and the mode closest to resonance is

[TABLE]

This is true regardless of whether $\Delta\omega_{Nm}<\Delta\omega_{\alpha}$ or whether $\Delta\omega_{Nm}>\Delta\omega_{\alpha}$. If there is no correlation between mode frequencies and forcing frequencies, the mean frequency difference is

[TABLE]

For ease in presentation, we define the non-resonant luminosity amplitude

[TABLE]

Note that some factors in this equation (e.g., $X_{Nm}$) are defined at forcing frequencies $\omega_{Nm}$, while others (e.g., $Q_{\alpha}$) are defined at mode frequencies $\omega_{\alpha}$. In practice, the values of $Q_{\alpha}$, $\Delta\omega_{\alpha}$, etc., are sufficiently smooth that one can estimate their values at a forcing frequency by interpolating between neighboring values of $\omega_{\alpha}$. This process may be more difficult in the sub-inertial regime ($\omega_{\alpha}<2\Omega_{s}$) where inertial modes make the mode spectrum very dense.

From equation 17, the most probable value of the luminosity fluctuation produced by a nearly-resonant mode is

[TABLE]

Equation 19 represents the most likely luminosity fluctuation amplitude, assuming the TEO is produced by a single near-resonant oscillation mode. The expected value of the amplitude is somewhat larger than equation 19, because of the high possible amplitudes attainable very close to resonance. The maximum possible luminosity fluctuation amplitude is

[TABLE]

We can also calculate the probability of a luminosity fluctuation $\Delta L_{N}/L$ at an orbital harmonic $N$ exceeding some critical value $\Delta L_{c}/L$, given by

[TABLE]

with $A_{N}$ given by equation 2. A little rearranging yields

[TABLE]

with

[TABLE]

Assuming no correlation between $\omega_{Nm}$ and $\omega_{\alpha}$, the probability that the frequency detuning $\delta\omega_{N}$ be less than some fraction $x$ of the maximum detuning $\Delta\omega_{\alpha}/2$ is simply equal to $x$, as long as $0<x<1$. Then we have

[TABLE]

Then, the expected number of luminosity fluctuations above some amplitude $\Delta L_{c}/L$, distributed across all orbital harmonics $N$, is

[TABLE]

In practical calculations, the sum of equation 25 must be truncated, and the value of $N_{\rm max}$ must be large enough to include all harmonics that can substantially contribute to the sum. In most cases, very low and high frequency TEOs do not substantially contribute to the sum, so the truncation is not problematic.

It is also useful to calculate the expected density of TEOs as a function of $A_{N}$. Inverting equation 1, we have

[TABLE]

Differentiation with respect to $A_{N}$ yields

[TABLE]

Now, for a given orbital harmonic N with a resonant detuning $\delta\omega_{N}$, the number of TEOs per unit detuning frequency is

[TABLE]

This can be understood from the fact that around each forcing frequency the most resonant mode always has $0<\delta\omega_{N}<\Delta\omega_{\alpha}/2$. Then we find that the probability density of modes as a function of amplitude is

[TABLE]

With a stellar/orbital model, this expression can be evaluated at each orbital harmonic N, providing a probability density of TEOs as a function of frequency and amplitude. One can then compare the theoretical TEO distribution with the observed distribution.

5 Tidal Dissipation and Resonance Locking

The angular momentum contained by a mode $\alpha$, together with its complex conjugate (c.f. Burkart et al. 2014), using our normalization condition in Section 6 is

[TABLE]

Differentiating with respect to time yields $\dot{J}_{\alpha}=2\gamma_{\alpha}J_{\alpha}$, and using $\dot{E}_{\alpha}=(N\Omega/m)\dot{J}_{\alpha}$, the orbital energy dissipation rate due to a tidally forced mode (and its complex conjugate) is

[TABLE]

Recalling that $\gamma_{\alpha}\!<\!0$ corresponds to a damped mode, $\dot{E}_{\rm orb,tide}$ may be either positive or negative, depending on the sign of $\omega_{\alpha}$, the mode frequency in the star’s rotating reference frame. Physically, this corresponds to the fact that the orbit may lose energy due to a resonance with a prograde mode (positive $\omega_{\alpha}$) but it may gain energy due to a resonance with a mode that is retrograde in the rotating frame (negative $\omega_{\alpha}$) but prograde in the inertial frame (positive $N\Omega$).

Using the relation $J_{\rm orb}=-2\sqrt{1-e^{2}}E_{\rm orb}/\Omega$, we can solve for the eccentricity evolution of the orbit. The resulting circularization time scale is

[TABLE]

where the orbital decay timescale is

[TABLE]

For a decaying orbit, $t_{\rm orb}$ is positive. In this case, the eccentricity is damped when $N\sqrt{1-e^{2}}>m$, but it is excited when $N\sqrt{1-e^{2}}<m$.

The luminosity fluctuation amplitude produced by a mode (together with its complex conjugate) is

[TABLE]

The energy dissipation rate of a mode can thus be measured from its observed amplitude,

[TABLE]

Hence, with a good stellar/orbital model such that one may calculate values of $\gamma_{\alpha}$, $L_{\alpha}$, etc., inserting observed TEO amplitudes into equation 35 allows one to calculate the tidal energy dissipation rate.

5.1 Resonance Locking

TEOs with higher amplitudes than expected from a statistical analysis (Section 4) are created by tidally forced modes unusually close to resonance. Such TEOs are good candidates for resonantly locked modes, the process of which is analyzed in detail by Witte & Savonije (1999); Fuller & Lai (2012a); Burkart et al. (2014). Here we discuss observational signatures of resonantly locked modes.

A mode is near resonance when one of its eigenfrequencies is approximately equal to a forcing frequency, $\omega_{\alpha}\simeq\omega_{Nm}$, or in the inertial frame,

[TABLE]

If the mode is locked in resonance, then this condition does not change with time, i.e.,

[TABLE]

where the dot denotes the time derivative.

Mode eigenfrequencies change with time for various reasons (e.g., stellar evolution, tidal spin-up or spin-down, magnetic braking, etc.). Similarly the orbital frequency may change due to several effects (e.g., orbital decay via gravitational waves, tidally induced orbital decay). For heartbeat stars beginning to evolve off the main sequence, we expect stellar evolution and tidal dissipation to be the dominant effects. Assuming the mode frequencies evolve only via stellar evolution and tidal spin-up, while the orbital frequency evolves due only to tidal dissipation, resonance locking requires

[TABLE]

where the first term on the left hand side accounts for the change in mode frequency due to stellar evolution and the second accounts for frequency change due to tidal torques on the star. For simplicity, we assume the star maintains rigid rotation at all times. The tidal spin term is (see Burkart et al. 2014)

[TABLE]

Here, $B_{\alpha}=(1/m)\partial\sigma_{\alpha}/\partial\Omega_{s}$, which evaluates to $B_{\alpha}\simeq 1-1/l(l+1)$ for asymptotic g-modes in the the slowly rotating limit, while $I$ is the star’s moment of inertia.

The energy dissipation rate $\dot{E}_{\rm tide}$ is the tidal energy dissipation rate in the star in the inertial frame and obeys the relation $\dot{E}_{\rm tide}=N\Omega\dot{J}_{\rm tide}$. Energy conservation requires $\dot{E}_{\rm tide}=-\dot{E}_{\rm orb,tide}$, where $E_{\rm orb}=-GMM^{\prime}/(2a)$ is the orbital energy. The orbital frequency and energy are related by

[TABLE]

so the resonance locking condition (equation 37) requires

[TABLE]

Equation 41 demonstrates a key result: if resonance locking occurs, the orbital energy loss is set by a time scale proportional to $\sigma_{\alpha}/\dot{\sigma}_{\alpha,*}$, i.e., a stellar evolution time scale. Therefore, if resonance locking occurs, the orbital decay/circularization timescales are set by the rate of stellar evolution. The term in brackets can be positive or negative, and resonantly locked modes could have either increasing or decreasing frequencies.

Inserting equation 31 into equation 41, we find the resonantly locked mode amplitude is

[TABLE]

where we have defined the mode evolution time scale due to stellar evolution,

[TABLE]

which may be either positive or negative. In equation 42,

[TABLE]

with $\kappa$ the dimensionless moment of inertia such that $I=\kappa MR^{2}$. The value of $\chi_{\alpha}$ is closely related to the critical harmonic $N_{c}$ discussed in Fuller & Lai (2012a) and the moment of inertia ratio $r$ from Burkart et al. (2014). Note that the mode amplitude is only real if the quantity in brackets in equation 42 is positive, this is equivalent to the existence of a fixed point as discussed in Burkart et al. (2014).

Equations 34 and 42 illustrate the purpose of this endeavor: if we know the properties of a system (stellar and orbital parameters, oscillation mode properties, etc.), then we can predict the luminosity fluctuation produced by a resonantly locked mode. Including the resonantly locked mode and its complex conjugate, the explected luminosity variation is

[TABLE]

where $c_{\alpha}=N\Omega/(\omega_{\alpha}\chi_{\alpha})$. Equation 42 also provides contstraints on the possible frequencies of resonantly locked modes, as only modes in limited frequency regimes will produce real values of $|a_{\alpha}|_{\rm lock}$.

The orbital evolution is easily calculated during a resonance lock. Equation 41 can be rearranged to yield

[TABLE]

When the first term of equation 44 dominates, this simplifies to

[TABLE]

which is also derived in Fuller et al. (2016). Equation 47 clearly exhibits that while resonance locking occurs, orbital decay proceeds on the timescale $t_{\alpha,{\rm ev}}$ at which the resonantly locked mode’s intrinsic frequency changes. Typically, this is comparable to the stellar evolution timescale.

Equation 42 can be easily generalized to include additional forms of stellar and orbital evolution. Incorporating these terms merely requires the replacement of $t_{\alpha,{\rm ev}}$ with a general evolution timescale $t_{\rm ev}$ given by

[TABLE]

where the second and third terms on the right hand side account for additional orbital and spin frequency changes.

We comment that resonance locks are not always stable, and Burkart et al. (2014) presents a detailed analysis of their stability. A key issue is that the above analysis assumes an “adiabatic approximation” for the mode amplitude, such that the mode amplitude is given by equation 68. This approximation neglects changes in mode energy, and results in the approximation $\dot{E}_{\rm tide}=-\dot{E}_{\rm orb}$ used above. This approximation is always valid when $\omega_{\alpha}/(\gamma_{\alpha}^{2}t_{\alpha,{\rm ev}})<1$. For typical mode frequencies $f_{\alpha}=1\,{\rm c/d}$ and evolution rates $t_{\alpha,{\rm ev}}=1\,{\rm Gyr}$, this criterion is satisfied for $\gamma_{\alpha}\!\gtrsim\!4\times 10^{-6}\,{\rm c/d}$. Comparison with Figure 2 shows that the adiabatic approximation is valid for sufficiently low frequency modes in relatively hot stars. It is less likely to be valid for higher frequency g modes and in cool stars with smaller damping rates. When the adiabatic approximation is not valid, resonance locking may still occur, but in a restricted parameter space and with more complex dynamics, as discussed in Burkart et al. (2014). In our companion paper (Fuller, 2017), we show that modes in KIC8164262 easily satisfy the “no backreaction” approximation (equation 53 of Burkart et al. 2014) even though those modes do not satisfy the strict adiabatic criterion above. The no backreaction approximation allows us to use the adiabatic approximation at the mode amplitude required to sustain a resonance lock, even though the adiabatic approximation could fail at higher amplitudes.

6 Dynamical Mode Amplitudes

The tidal response of the star can be generally separated into two components: a hydrostatic equilibrium tide and a non-hydrostatic dynamical tide. The equilibrium tide corresponds to the perturbation induced in the star by a static external gravitational potential, and causes the star to distort into an elliptical shape so that hydrostatic equilibrium is maintained. The equilibrium tide is usually responsible for tidal ellipsoidal variations in close binaries, and creates part of the near-periastron luminosity variation observed in heartbeat stars.

In realistic systems, the external gravitational potential is not static (unless the orbit is both circular and synchronized), and the instantaneous amplitude of the equilibrium tide fluctuates. However, because of the finite inertia of the fluid, the fluctuating equilibrium tide induces a dynamical tidal response which can be extremely important for dissipative tidal processes. To understand the relation between the equilibrium tidal response and the dynamical tidal response, consider the linearized momentum equation in the star’s rotating reference frame:

[TABLE]

where $U$ is the external gravitational potential, and \mathcal{C}{\mbox{\boldmath\xi}} represents internal forces from the perturbed fluid. A free oscillation mode (indexed by the label $\alpha$) of the star satisfies

[TABLE]

for a mode with time dependence $\propto e^{-i\omega_{\alpha}t}$ (see Lai & Wu 2006). This convention implies that a mode with positive frequency and $m>0$ is prograde in the rotating frame of the star, while modes with positive frequency and $m<0$ are retrograde in the rotating frame (although they may be prograde in the inertial frame).

In non-rotating stars, modes obey the orthonormality condition

[TABLE]

where the integral extends over the volume of the star. However, in rotating stars (c.f. Schenk et al. 2002; Lai & Wu 2006) modes obey the modified orthonomality condition

[TABLE]

where

[TABLE]

is a Coriolis coupling coefficient. This modified orthonomality condition is important in accurately calculating the dynamical tidal response, especially when rotation strongly modifies the modes of interest.

A general perturbation to the star must be expanded in terms of its displacement and velocity perturbations via a phase space mode expansion:

[TABLE]

where $a_{\beta}$ is the time-dependent mode amplitude. The summation runs over both positive and negative values of the mode azimuthal number $m$, and also over both positive and negative eigenfrequencies $\omega_{\beta}$. Since the mode eigenfunction contains no time-dependence, equation 54 implies that for any perturbation,

[TABLE]

although individual terms in the summations are not generally equal to one another for a forced mode. Equation 55 additionally leads to the relation

[TABLE]

To obtain the tidal response of an oscillation mode (indexed by $\alpha$), we insert the mode expansion 54 into equation 49, and take the inner product with {\mbox{\boldmath\xi}}_{\alpha}. Using the orthonormality condition 52 and the identity of equation 56, one obtains the mode amplitude equation (Schenk et al., 2002)

[TABLE]

This equation describes the evolution of the mode amplitude $a_{\alpha}$, which includes both its equilibrium tide component and its dynamical tide component. Its precise form depends on the adopted normalization condition, which determines the meaning of the mode amplitude $a_{\alpha}$.

To understand tidal dynamics, it is useful to separate the equilibrium and dynamical tidal responses. To do this, we must subtract out the equilibrium tidal response from equation 49. The equilibrium tide is defined by taking the static ($\partial/\partial t\rightarrow 0$) limit of equation 49:

[TABLE]

Equation 58 also implies that \partial{\mbox{\boldmath\xi}}_{\rm eq}/\partial t=-i\omega{\mbox{\boldmath\xi}}_{\rm eq} for a tidal potential of form $U\propto e^{-i\omega t}$.

We then separate the dynamical and equilibrium tide components, {\mbox{\boldmath\xi}}={\mbox{\boldmath\xi}}_{\rm dyn}+{\mbox{\boldmath\xi}}_{\rm eq}, which upon substitution into equation 49 yields

[TABLE]

Next, we decompose the dynamical and equilibrium responses in the same fashion as equation 54:

[TABLE]

and likewise for the equilibrium tide. Then, using the identity

[TABLE]

and taking the same steps used to derive equation 57, we find

[TABLE]

From equation 58 one can also find

[TABLE]

and so

[TABLE]

Equation 64 describes the evolution of the dynamical component of the mode amplitude $a_{\alpha,{\rm dyn}}$. It contains two forcing terms on the right-hand side, which we discuss below.

The adiabatic solution to equation 64 is obtained by assuming the magnitude of the mode amplitude changes slowly such that $\dot{a}_{\alpha,{\rm dyn}}=-i\omega a_{\alpha,{\rm dyn}}$. We allow for complex eigenfrequencies by replacing $\omega_{\alpha}$ with $\omega_{\alpha}+i\gamma_{\alpha}$, where $\gamma_{\alpha}$ is the mode growth rate which has been defined such that negative values of $\gamma_{\alpha}$ correspond to damped modes. The adiabatic mode amplitude is

[TABLE]

Here, we have assumed for simplicity that the mode is weakly damped, i.e., $\gamma_{\alpha}\ll\omega_{\alpha}$, which is always the case for detectable modes in stars.

Equation 65 describes the dynamical response of a mode $\alpha$ to tidal forcing. The first term is the same as the total mode amplitudes discussed in previous works, but multiplied by a factor $\omega/\omega_{\alpha}$. This term is largest near resonances where $\omega\simeq\omega_{\alpha}$. The multiplicative factor $\omega/\omega_{\alpha}$ does not strongly affect nearly resonant modes, but is important for cancellations between non-resonant modes. As we show below, these cancellations entail that high frequency (relative to the tidal forcing frequency) f-modes which dominate the equilibrium tide response have little contribution to the dynamical tide.

The second term in equation 65 is due to damping of the equilibrium tide. Upon summation over both a mode and its complex conjugate, it produces a term

[TABLE]

This term accounts for the lag of the equilibrium tide due to damping, generating a phase shift in the equilibrium tidal bulge. The summation of this term over all modes accounts for the classical tidal lag angle $\delta\propto 1/Q$ associated with the equilibrium tide, where $Q$ is the equilibrium tide quality factor.

Finally, recall that equation 65 is only the forced mode amplitude. The free oscillation modes of the star satisfy

[TABLE]

Thus, the stellar oscillation modes are still permitted to oscillate at their natural oscillation frequency, $\omega_{\alpha}$, at an arbitrary amplitude. Therefore, stellar oscillation modes of stars in close binaries have three simultaneous amplitudes: $a_{\rm free}$, $a_{\rm eq}$, and $a_{\rm dyn}$.

Using similar notation as Fuller & Lai (2012a), the dynamical part of the mode amplitude (computed from the first term of equation 65, in the rotating frame of the star for aligned spin and orbit) is then

[TABLE]

where

[TABLE]

describes the strength of the tidal forcing due to the component of the tidal potential with indices $l$ and $m$, $W_{22}=W_{2-2}=\sqrt{3\pi/10}$, and $W_{20}=-\sqrt{\pi/5}$. We generalize to the misaligned case in Section 8.

The formalism used above is strictly only valid for adiabatic modes, as non-adiabatic modes do not obey the same orthogonality conditions. We suspect non-adiabatic corrections to equation 64 will be of order $\gamma_{\alpha}a_{\alpha}$. For modes far from resonance, we expect these corrections to be negligible. However, for highly damped modes with $\gamma_{\alpha}\sim\Delta\omega_{\alpha}$, i.e., modes near the traveling wave regime, the corrections may be significant. There, the physical picture is a train of tidally excited gravity waves propagating into the stellar envelope, as examined for early type stars (Zahn, 1975, 1977) and white dwarfs (Fuller & Lai, 2012b; Burkart et al., 2013). In this case, there are no resonances in the tidal response, and the wave amplitude is a smoothly varying function of frequency. This situation is similar to the result of equation 65, because the amplitude of the resonant detuning term is smoothly varying and can only vary by a factor of $\sim 2$ when $\gamma_{\alpha}\sim\Delta\omega_{\alpha}$. Hence, in this circumstance, we suspect our formalism still yields approximately correct results (at the factor of 2 level), but it cannot be used for precise predictions of TEO amplitudes or phases.

7 Summation over Oscillation Modes

The forced response at a given orbital harmonic $N$ is found by summing over all modes:

[TABLE]

This sum over $\alpha$ includes all modes which significantly couple with the tidal potential. For a given value of $m$, this includes modes with both positive and negative frequencies $\omega_{\alpha}$. Note that a mode with $(\omega_{\alpha},m)$ has a physically identical counterpart with $(-\omega_{\alpha},-m)$, and both must be formally included in the summation. We will refer to the former mode as the $\omega_{\alpha,m}$ mode and the latter as the $\omega_{-\alpha,-m}$ mode. We must also include forcing at both $N$ and $-N$ since the decomposition of the tidal potential into orbital harmonics includes both positive and negative forcing frequencies. The observed response at a harmonic $|N|$ is a combination of both the positive and negative forcing frequencies. In what follows, we will consider the $m=0$ and $m=\pm 2$ components of the $l=2$ tidal potential.222As is commonly done, we ignore components of the tidal potential with $l>2$ because the strength of the potential scales as $(R/a)^{(l+1)}$. We only compute the contributions from $m=0$ and $m=\pm 2$ modes because the $m=\pm 1$ contributions vanish if the spin and orbit are aligned. However, for misaligned spins, the $m=\pm 1$ terms should be included.

Including the contributions listed above, the tidal response at a harmonic $|N|$ of the orbital frequency is, in the inertial frame, for $|m|>0$ modes

[TABLE]

Here, we have factored out the angular dependence of the eigenfunctions in the traditional approximation, {\mbox{\boldmath\xi}}_{\alpha}={\mbox{\boldmath\xi}}_{\alpha}(r)H_{\alpha m}(\theta)e^{im\phi}, where $H_{\alpha m}(\theta)$ is the Hough function associated with a mode $\alpha$ (see Section 9). We have also omitted the retrograde terms (those with dependence $\propto e^{\pm i(N\Omega t+m\phi)}$) because they are typically excited to very low amplitudes since the value of $F_{N-m}=F_{-Nm}$ is very small for any realistic value of the orbital eccentricity. Here, retrograde signifies modes that are retrograde in the inertial frame, i.e., modes whose patterns revolve in the opposite direction of the orbital motion. We may still obtain large amplitude modes that are prograde in the inertial frame, yet are retrograde in rotating frame of the star.

To simplify equation 7, we use the fact that $\omega_{Nm}=-\omega_{-N-m}$, $\omega_{\alpha m}=-\omega_{-\alpha-m}$, $F_{Nm}=F_{-N-m}$, $Q_{\alpha m}=Q_{-\alpha-m}$, $H_{\alpha m}(\theta)=H_{-\alpha-m}(\theta)$, {\mbox{\boldmath\xi}}_{\alpha m}(r)={\mbox{\boldmath\xi}}_{-\alpha-m}(r), $\gamma_{\alpha m}=\gamma_{-\alpha-m}$, and likewise for $(-\alpha,m)$ and $(\alpha,-m)$ combinations. After careful manipulation we obtain

[TABLE]

with

[TABLE]

and

[TABLE]

Here, the ${\rm atan2}$ function must be used to obtain the correct phase shift. These equations assume that the mode displacement $\xi$ is a purely real quantity, i.e., there are no non-adiabatic effects. We examine non-adiabatic effects in more detail in Section 7.1.

The first line in equation 7 represents prograde modes forced in the prograde direction, while the second term represents retrograde modes forced in the prograde direction. In some previous works, the second line has been omitted. However, it is important to retain it for two reasons. First, it can still to lead to resonant excitation because the forcing frequency $\omega_{Nm}=N\Omega-m\Omega_{s}$ is negative for $N\Omega<m\Omega_{s}$. Physically this corresponds to retrograde modes (in the rotating frame) being forced in the prograde direction (in the inertial frame). Second, terms on the first and second lines tend to cancel each other (due to their nearly opposite phase shifts) for high frequency modes with $\omega_{\alpha m},\omega_{\alpha-m}\gg\omega_{Nm}$. Consequently the stellar f modes that dominate the equilibrium tide deconstructively interfere such that their contribution to the dynamical tide is substantially reduced.

The same approach can be used for $m=0$ modes. In this case there are no terms corresponding to $(\alpha,-m)$ combinations, but the $(-N,m)$ combinations cannot be dropped because $F_{-Nm}=F_{Nm}$ for $m=0$. The result is

[TABLE]

with

[TABLE]

and

[TABLE]

Once again, the second term in equation 7 has been dropped in some previous works. In this case, it cannot lead to resonant excitation, although it must still be retained in order to obtain the deconstructive interference described above. Note that the form of equations 7-77 is the same as equations 7-74 without the retrograde terms, hence the latter can be used for $m=0$ modes.

7.1 Luminosity Fluctuations and Non-adiabatic Effects

TEOs can be observed with high-accuracy photometric data, for which the observable feature of TEOs is a disk-integrated luminosity fluctuation $\Delta L/L$. Both the amplitude and phase of these oscillations can be calculated, so a single eccentric binary can yield large numbers of observable quantities (the orbital harmonics, amplitudes, and phases of TEOs) which may be compared with tidal theories.

A complication not discussed above is that TEOs may become strongly non-adiabatic near the surface of the star where they are observed. The mode eigenfunction {\mbox{\boldmath\xi}}(r) will thus be phase-shifted relative to an adiabatic mode. This is equivalent to the modes obtaining imaginary components to their eigenfunctions. Then any mode quantity, e.g., the radial component of the displacement $\xi_{r}$, may be expressed as

[TABLE]

The mode overlap integral $Q_{\alpha}$ also obtains an imaginary component. Recomputing the sums above, we find

[TABLE]

with

[TABLE]

and

[TABLE]

Here, \big{|}Q_{\alpha m}\xi_{r,\alpha m}(r)\big{|}=\big{[}{(Q_{\alpha m}\xi_{r,\alpha m})^{\rm R}}^{2}+{(Q_{\alpha m}\xi_{r,\alpha m})^{\rm I}}^{2}\big{]}^{1/2}, and $(Q_{\alpha m}\xi_{r,\alpha m})^{\rm R}=Q_{\alpha m}^{\rm R}\xi_{r,\alpha m}^{\rm R}-Q_{\alpha m}^{\rm I}\xi_{r,\alpha m}^{\rm I}$, and $(Q_{\alpha m}\xi_{r,\alpha m})^{\rm I}=Q_{\alpha m}^{\rm R}\xi_{r,\alpha m}^{\rm I}+Q_{\alpha m}^{\rm I}\xi_{r,\alpha m}^{\rm R}$. A similar equation holds for any perturbed quantity, e.g., the perturbed flux $\Delta F/F$, but with $\xi_{r}$ replaced by $\Delta F/F$ at each point in equations 7.1-81. These equations also apply for $m=0$ modes, except for the additional retrograde terms.

7.1.1 Observed Luminosity Fluctuation

The observable luminosity fluctuations are disk integrated quantities containing three contributing effects: surface area perturbations, surface normal perturbations, and flux perturbations. These have been examined in Buta & Smith (1979); Robinson et al. (1982); Townsend (2003), but here we follow the procedure of Burkart et al. (2012). For a luminosity fluctuation due to a single mode with spherical harmonic angular dependence, the magnitude of the luminosity fluctuation is

[TABLE]

where $b_{l}$ and $c_{l}$ are limb darkening coefficients and $i_{s}$ and $\phi_{s}$ are the angular coordinates of the observer in the star’s rotating frame at $t=0$. Note that $\phi_{s}$ is related to the argument of periapsis $\omega_{\rm p}$ of the secondary via

[TABLE]

From radial velocity measurements, one typically measures the argument of periapsis of the orbit of the primary, which is shifted by $\pi$ relative to $\omega_{\rm p}$ of the secondary.

The quantity $\Delta F_{\alpha}$ is the perturbed flux. In rotating stars, modes have Hough function angular dependence and so the result will change somewhat. However, we note that any Hough function may always be decomposed in terms of spherical harmonics,

[TABLE]

where $h_{klm}$ is the angular overlap integral of equation 132. Because high $l$ spherical harmonics suffer large cancelation effects when integrated over the stellar disk, only the low $l$ components of a Hough function will contribute appreciably to the disk-averaged luminosity fluctuation of a mode. For TEOs, the $l=2$ components will dominate. Therefore, to good approximation, the luminosity fluctuation of a TEO in the traditional approximation is

[TABLE]

Using equation 85, we find that the net luminosity fluctuation at a harmonic $N$ of the orbital frequency is

[TABLE]

with

[TABLE]

and

[TABLE]

and

[TABLE]

The sum of sinusoidal oscillations in equation 7.1.1 may be written as a single sinusoidal oscillation of form

[TABLE]

where the amplitude of the oscillation is

[TABLE]

and the amplitudes $A_{i}$ are given by the factors in front of the $\sin$ terms in equation 7.1.1, i.e.,

[TABLE]

The phases $\Delta_{i}$ are given by the phases of the $\sin$ terms in equation 7.1.1,

[TABLE]

Additionally, the phase $\delta_{N}$ of the observed oscillation is

[TABLE]

The sums over the indices $i$ and $j$ in equations 91 and 94 must run over all terms that contribute to the luminosity fluctuation at orbital harmonic $N$. This includes terms in equation 7.1.1 for both $m=2$ and $m=0$ modes, which may interfere with another. In systems where both stars are pulsating, contributions from both stars must be included. In many cases, however, one star may dominate the pulsation signal, and either the $m=0$ or $|m|=2$ terms may dominate, and we may limit the sums to that value of $m$.333Of course, the $l>2$ components of the tidal potential also contribute, as do additional branches of Hough functions (corresponding to $l>2$ modes in the non-rotating limit). However, as discussed above, these contributions are small in many circumstances.

The frequency $N\Omega$, luminosity fluctuation $A_{N}$, and phase $\delta_{N}$ from equation 90 represent the observable quantities for TEOs in photometric data. Recall that our expansion of the time dependence of the tidal potential into orbital harmonics requires that $t=0$ corresponds to the time of periastron. Therefore, reliably measuring mode phases from photometric data requires a determination of both the epoch of periastron and the argument of periapsis. As discussed in Section 3, the value of $L_{\alpha}$ generally has a substantial imaginary component, whereas $\gamma_{\alpha}$ is typically much smaller than the resonant detuning, unless the mode is extremely close to resonance. Consequently, the observed phase shift will almost always be dominated by non-adiabatic effects that contribute to the imaginary part of $L_{\alpha}$.

It is apparent from the tedious summations above that the relation between observed luminosity fluctuations and stellar properties is not transparent. The values $Q_{\alpha}$, $L_{\alpha}$, $\omega_{\alpha}$, and $\gamma_{\alpha}$ must be calculated for the modes of a given stellar model, and can be somewhat sensitive to the stellar parameters. Furthermore, the values of $\epsilon_{lm}$, $F_{Nm}$, $Y_{lm}(i_{s},0)$, $\omega_{Nm}$ require good measurements of the orbital parameters (and stellar masses, radii, and spins) to be accurately calculated. Most important, the resonant detuning term $\omega_{Nm}/\sqrt{(\omega_{\alpha m}-\omega_{Nm})^{2}+\gamma_{\alpha m}^{2}}$ is extremely sensitive to the precise values of the mode frequencies $\omega_{\alpha}$. For these reasons, reliably performing tidal asteroseismology (i.e., using the mode frequencies, amplitudes, and phases to constrain stellar parameters) is difficult in practice.

8 Accounting for Spin-Orbit Misalignment

Here we expand the formalism above to allow for misalignment between the stellar spin axis and the orbital angular momentum axis. This has been discussed briefly in Lai (1997), here we revisit this calculation to predict amplitudes and phases of tidally excited modes in misaligned binaries.

We begin from the amplitude equation 64 which applies in the rotating frame of the star with $\hat{z}_{\rm spin}$ in the direction of the spin angular momentum,

[TABLE]

We have dropped the ${\rm dyn}$ subscript for simplicity. Performing the standard decomposition of the tidal potential $U$ in the rotating frame gives

[TABLE]

Here, ($D^{\prime}$,$\theta^{\prime}$,$\phi^{\prime}$) is the time-dependent location of the secondary star in the rotating frame of the primary. Using $Y_{l-m}=(-1)^{m}Y_{lm}^{*}$, each component of the tidal potential has form

[TABLE]

To express quantities in reference frames rotated relative to one another, we expand in Wigner $\mathcal{D}$ functions,

[TABLE]

Here, $\alpha$, $\beta$, and $\gamma$ are Euler angles describing the rotations necessary to transform the new orbital frame to the rotating frame, and $\theta^{\prime}_{o}$,$\phi^{\prime}_{o}$ are the location of the secondary in the new reference frame. We define the orbital frame as having $\hat{z}_{\rm orb}$ in the direction of the orbital angular momentum vector, and $\hat{x}_{\rm orb}$ in the direction of the secondary at periastron. The value of the Wigner function is

[TABLE]

where $d^{l}_{m_{o}m}(\beta)$ is an element of Wigner’s small d-matrix. We caution that the definitions of Wigner functions and Euler angles is dependent on the chosen convention for the series of rotations to transform between coordinate systems. Equation 99 uses the quantum mechanics convention used in Varshalovich et al. (1988) and adopts a z-y-z convention for the series of rotations $\alpha,\beta,\gamma$, which differ from the classical Euler angles used below that follow a z-x-z convention.

To transform the orbital frame to the rotating frame, we can rotate it about the $z_{\rm orb}$ axis by an angle $\alpha_{o}$, then around its new $x-$axis by an angle $\beta_{o}$, then about its new $z$-axis by an angle $\gamma_{o}$. We define the $\hat{x}_{\rm spin}$ of the rotating frame such that it lies in the orbital plane at time $t=0$. The angle $\alpha_{o}$ is the angle between $\hat{x}_{\rm spin}$ and $\hat{x}_{\rm orb}$ at $t=0$, and is the longitudinal spin-orbit misalignment measured in the orbital plane. The angle $\beta_{o}$ is the misalignment angle between $\hat{z}_{\rm spin}$ and $\hat{z}_{\rm orb}$, and our definition of $\hat{x}_{\rm spin}$ entails $\gamma_{o}=\Omega_{s}t$, i.e., it increases with time due to the rotation of the star. Note that this series of rotations uses a z-x-z sequence corresponding to classical Euler angles that is different from the quantum convention used in the $\mathcal{D}$ functions above. The two are equivalent if we use $\alpha=\alpha_{o}-\pi/2$, $\beta=\beta_{o}$, and $\gamma=\gamma_{o}+\pi/2$. Figure 6 shows a diagram of our reference frames.

The tidal potential in the rotating frame is thus

[TABLE]

In the orbital frame, the secondary lies in the equator at $\theta^{\prime}_{o}=\pi/2$ such that only even values of $m_{o}$ contribute to the sum. Inserting this expression into equation 95, we obtain

[TABLE]

As above, we can decompose the forcing terms into orbital harmonics to obtain the forced mode amplitude at each orbital harmonic $N$

[TABLE]

with $F_{Nm}$ defined in equation 5.

Then, as above, the forced amplitude of each mode at each orbital harmonic $N$ is

[TABLE]

For arbitrary spin-orbit misalignment, odd $m$ terms can be excited, in contrast the aligned case. The misalignment angle $\beta_{o}$ modifies the strength of the forcing for different values of $m$, while the angle $\alpha_{o}$ induces a phase shift.

To simplify equation 101, we can combine the terms in the sum into a single forcing term:

[TABLE]

with amplitude $X_{Nm}$ of

[TABLE]

and phase shift

[TABLE]

In the limit of aligned spin and orbit, $\beta_{o}\rightarrow 0$ and we can choose $\alpha_{o}=0$. Then $d^{l}_{m_{o}m}(\beta_{o})\rightarrow\delta_{m_{o}m}$, $X_{Nm}\rightarrow W_{lm}F_{Nm}$, $x_{Nm}\rightarrow m\pi/2$, and equation 103 reduces to equation 68.

8.1 Observed Pulsation Amplitude and Phase

At any instant, the observed luminosity variation due to a mode forced at orbital harmonic $N$ is

[TABLE]

To obtain the observed luminosity fluctuation we must calculate $Y_{2m}(i_{s},\phi_{s})$. To do this, we can first transform back into the orbital frame,

[TABLE]

with ($\theta_{\rm obs}$,$\phi_{\rm obs}$) the angular coordinate of the observer measured in the orbit-centered frame described in Section 8. Then we can see $\theta_{\rm obs}=i_{\rm orb}$, where $i_{\rm orb}$ is the orbital inclination to the line of sight, and $\phi_{\rm obs}=-\omega_{p}-\pi/2$. Hence,

[TABLE]

Similar to the forcing terms, we can express this visibility term as an amplitude and phase shift,

[TABLE]

with amplitude

[TABLE]

and phase shift

[TABLE]

Inserting equation 110 into equation 107 yields the observed luminosity variation due to a tidally forced mode,

[TABLE]

Note the time-dependence of the spin-term has canceled back out upon transformation to the observer’s frame, so that tidally excited modes will be observed at integer harmonics of the orbital frequency regardless of the spin-orbit misalignment.

Now, as above, we must sum over modes of different signs of $m$, $\omega_{\alpha}$, and $N$. This summation yields

[TABLE]

In the aligned case, the first line represents prograde modes being forced in the same direction as the orbital motion, while the second line represents retrograde modes being forced in the direction of orbital motion. The third line represents retrograde modes being forced opposite to the direction of orbital motion, while the fourth line represents prograde modes being forced opposite to the direction of orbital motion. In the aligned case, the value of $X_{N-m}$ is very small and the third and fourth lines are negligible. However, these terms can be important for substantial misalignment. Note that $\omega_{Nm}$ can be positive or negative, but $\omega_{N-m}>0$ and $\omega_{\alpha}>0$ for all $N$ and $m$ in these sums. The first three terms in equation 8.1 can produce resonantly excited modes, but the fourth term cannot (regardless of the spin-orbit alignment), so this term is less important. In the aligned limit, $V_{lm}\rightarrow Y_{lm}(i_{s},0)$, $v_{lm}\rightarrow m(\alpha_{o}+\omega_{p})$, $x_{Nm}-v_{lm}\rightarrow m(-\pi/2-\omega_{p})=m\phi_{s}$, and equation 8.1 reduces to equation 7.1.1.

8.2 Determining Spin-Orbit Misalignment

In practice, it can be difficult to observationally determine the spin-orbit misalignment angles $\alpha_{o}$ and $\beta_{o}$. However, they can be constrained if there is a measurement of the projected spin-orbit misalignment via the Rossiter-McLaughlin effect, or if there are measurements of rotation period, stellar radius, and $v\sin i_{s}$. We use a method similar to Fabrycky & Winn (2009) to present relations between $\alpha_{o}$, $\beta_{o}$, and observable quantities.

In addition to the spin and orbit-oriented reference frames discussed above, we now consider the observer’s reference frame with $\hat{z}_{\rm obs}$ pointing toward the observer. The $x_{\rm obs}-y_{\rm obs}$ plane is the plane of the sky, and we choose to orient this plane such that the $x_{\rm obs}-$axis lies on the intersection between the orbital plane and the plane of the sky. In this frame, the orbital axis has unit vector $\hat{n}_{\rm orb}=\sin i_{\rm orb}\hat{y}_{\rm obs}+\cos i_{\rm orb}\hat{z}_{\rm obs}$. The stellar spin axis has unit vector $\hat{n}_{s}=\sin i_{s}\sin\lambda\hat{x}_{\rm obs}+\sin i_{s}\cos\lambda\hat{y}_{\rm obs}+\cos i_{s}\hat{z}_{\rm obs}$, where the angle $\lambda$ is the projected misalignment between the spin and orbital axis.

Let us also consider a reference frame defined relative to the orbit with $z_{\rm orb^{\prime}}-$axis perpendicular to the orbital plane and $x_{\rm orb^{\prime}}-$axis in the same direction as the $x_{\rm obs}$-axis above. In this frame, $\hat{n}_{\rm orb}=\hat{z}_{\rm orb^{\prime}}$ and $\hat{n}_{s}=\sin\beta_{o}\sin(\alpha_{o}+\omega_{p})\hat{x}_{\rm orb^{\prime}}-\sin\beta_{o}\cos(\alpha_{o}+\omega_{p})\hat{y}_{\rm orb^{\prime}}+\cos\beta_{o}\hat{z}_{\rm orb^{\prime}}$. The argument of periapse $\omega_{p}$ is the angle between $\hat{x}_{\rm orb}$ and $\hat{x}_{\rm orb^{\prime}}$. The obs and orb’ frames are related by a rotation about the $\hat{x}_{\rm obs}=\hat{x}_{\rm orb^{\prime}}$ axis by an angle $i_{\rm orb}$, such that $\hat{y}_{\rm orb^{\prime}}=\cos i_{\rm orb}\hat{y}_{\rm obs}-\sin i_{\rm orb}\hat{z}_{\rm obs}$, and $\hat{z}_{\rm orb^{\prime}}=\sin i_{\rm orb}\hat{y}_{\rm obs}+\cos i_{\rm orb}\hat{z}_{\rm obs}$. Equating the unit vector of the spin axis in each frame thus requires

[TABLE]

[TABLE]

[TABLE]

It is not generally possible to determine both $\alpha_{o}$ and $\beta_{o}$ from observed quantities because the above equations are not linearly independent. However, equation 117 relates the values of $\alpha_{o}$ and $\beta_{o}$, provided that $i_{s}$, $i_{\rm orb}$, and $\omega_{p}$ can be measured,

[TABLE]

9 Modes in the Traditional Approximation

Gravity modes in rapidly rotating stars are frequently studied in the traditional approximation (see Bildsten et al. 1996; Lee & Saio 1997; Townsend 2003 and references therein), which allows one to separate the radial and angular dependence of oscillation modes. Observed mode frequencies appear to conform very well to predictions of the traditional approximation (Moravveji et al., 2016; Van Reeth et al., 2016), although it fails to capture some dynamics of sub-inertial modes in convective zones (Mathis et al., 2014). The latitudinal depedendence of modes is no longer described by associated Legendre polynomials, but rather by Hough functions. The angular dependence of a mode $\alpha$ is found by solving Laplace’s tidal equation:

[TABLE]

where $\lambda_{k}$ is an angular eigenvalue, $H_{km}$ is its associated Hough function, and the operator $\mathcal{L}$ is

[TABLE]

with $\mu=\cos(\theta)$. The parameter $q$ that determines the behavior of the Hough functions is

[TABLE]

Rotation becomes very important for $q\gtrsim 1$, and it is easily verified that the solutions to equation 119 converge to the associated Legendre polynomials as $q\rightarrow 0$. Similar to spherical harmonics, we normalize the Hough functions via

[TABLE]

with the integral taken over a spherical surface.

Because the traditional approximation only works for g-modes, which exhibit small gravitational perturbations, it is almost always combined with the Cowling approximation, in which the gravitational perturbations of the modes are ignored. However, in tidal applications, the small (but finite) gravitational potential perturbations created by a g-mode determine the value of the overlap integral $Q_{\alpha}$ (equation 4). While it is still possible to approximately calculate this overlap integral in the Cowling approximation using equation 4, this equation is susceptible to numerical inaccuracies (see discussion in Fuller & Lai 2011), and we find it unsuitable to reliably calculate overlap integrals for low frequency g-modes.

Therefore, we choose to compute mode eigenfunctions using the traditional approximation but not using the Cowling approximation. To do this, we examine the linearized Poisson equation:

[TABLE]

where $\delta\Phi(r,\theta,\phi)$ is the Eulerian gravitational potential perturbation, $\delta\rho(r,\theta,\phi)$ is the Eulerian density perturbation, and $\nabla_{\perp}^{2}$ is the angular part of the Laplacian. There arises an immediate problem with calculating gravitational potential perturbations of modes in the traditional approximation. Although $\delta\rho(r,\theta,\phi)$ has Hough function angular dependence, the quantity $\nabla_{\perp}^{2}\delta\Phi(r,\theta,\phi)$ will in general not be described by a single Hough function.444The exception to this rule occurs when $q=0$, and the Hough function dependence of $\delta\rho$ reduces to an associated Legendre polynomial. Then, the angular operator $\mathcal{L}=r^{2}\nabla_{\perp}^{2}$, and both $\nabla_{\perp}^{2}\delta\Phi(r,\theta,\phi)$ and $\delta\rho(r,\theta,\phi)$ will be proportional to associated Legendre polynomials. Consequently, the modes are not orthogonal because they will exhert gravitational torques on one another, and the radial and horizontal dependence of $\delta\Phi(r,\theta,\phi)$ are not separable.

In what follows, we ignore this issue and simply assume $\nabla_{\perp}^{2}\delta\Phi(r,\theta,\phi)$ has the same Hough function angular dependence as $\delta\rho(r,\theta,\phi)$. We believe this approximation is justified for g-modes because $\delta\Phi$ is very small and does not strongly affect the mode behavior. For higher frequency modes (f modes and p modes), the approximation is also valid because they typically have $q\ll 1$ such that their angular dependence is nearly equal to a spherical harmonic.

Letting $\delta\Phi(r,\theta,\phi)=\delta\Phi(r)H_{km}(\theta)e^{im\phi}$, equation 123 reduces to

[TABLE]

where

[TABLE]

is an effective angular wave number. For $q\rightarrow 0$, $z_{km}\rightarrow-l(l+1)$. The Hough function dependence of $\delta\Phi$ changes the boundary conditions used to calculate the mode eigenfunctions. At both $r\rightarrow 0$ and $r\rightarrow R$, equation 124 reduces to (for non-radial modes)

[TABLE]

Letting $\delta\Phi\propto r^{b}$, we find $b=-1/2\pm(1/2)\sqrt{1-4z_{km}}$. Then regularity requires that

[TABLE]

and

[TABLE]

where $b_{+}=-1/2+(1/2)\sqrt{1-4z_{km}}$ and $b_{-}=-1/2-(1/2)\sqrt{1-4z_{km}}$. As $q\rightarrow 0$, we obtain the usual dependence $b_{+}\rightarrow l$ and $b_{-}\rightarrow-(l+1)$. There are additional subtleties on the inner boundary condition for $\xi_{r}$ because the traditional approximation breaks down as $r\rightarrow 0$. In this study, we use the same approach described in Fuller & Lai (2014).

With an explicit calculation of $\delta\Phi_{\alpha}$ for each mode, we may now use equation 4 to calculate the value of $Q_{\alpha}$ for each mode. Using the continuity equation and Poisson’s equation, one can show

[TABLE]

where we have assumed $\rho(R)\simeq 0$ to obtain the second line. Then, using Green’s second identity, we have

[TABLE]

Then, since $\nabla^{2}U=0$ for a tidal potential, we have

[TABLE]

where $h_{klm}$ describes the angular overlap between a Hough function and spherical harmonic,

[TABLE]

Then we find from equation 4 that

[TABLE]

for all quantities expressed in dimensionless units ($G=M=R=1$), which reduces to $Q_{\alpha}=-(2l+1)\delta\Phi_{\alpha}(R)/(4\pi\omega_{\alpha}^{2})$ in the non-rotating limit.

10 Discussion and Conclusions

We have examined the visible luminosity oscillations produced by dynamical tides in eccentric binary systems known as heartbeat stars. The signature of dynamical tides is a stable oscillation at an exact integer multiple of the orbital frequency. Most observable tidally excited oscillations (TEOs) are produced by gravity modes (g modes) within one or both stars that are resonantly forced, i.e., a g mode frequency is nearly equal to an integer multiple of the orbital frequency. TEOs are expected to have higher amplitudes in hot stars with $T_{\rm ef}\gtrsim 6500\,{\rm K}$, because the absence of surface convection zones in hot stars allows g modes to propagate much closer to the surface and produce larger surface temperature perturbations. In principle, it is straightforward to calculate the luminosity oscillations produced by TEOs given accurate stellar and orbital parameters. We have provided precise formulae to predict TEO amplitudes, frequencies, and phases, including Coriolis forces, non-adiabatic effects, and spin-orbit misalignment.

In practice, however, the uncertainties in stellar/orbital parameters make a precise calculation difficult because mode amplitudes are extremely sensitive to the resonant detuning between stellar oscillation frequencies and tidal forcing frequencies. It is often more constructive to compare observed and expected TEOs using a statistical framework. In this approach, one can compare the number of observed TEOs exceeding a given threshold to the expected number of oscillations. One can also compare observed frequencies with the range in which TEOs are expected to be observed. This probabilistic approach assumes that tidal forcing frequencies (which occur at integer multiples of the orbital frequency) and oscillation mode frequencies are uncorrelated, i.e., the TEOs have no backreaction on the orbit. In reality, some feedback does occur because TEOs dissipate energy and cause tidal orbital evolution.

If a heartbeat star exhibits a very large amplitude TEO (e.g., KIC 8164262, see companion papers Hambleton et al. 2017; Fuller 2017) that is unlikely to stem from a chance resonance, it is a good candidate to be a resonantly locked mode. Resonantly locked modes arise from feedback between stellar evolution and TEOs, exciting a mode to large amplitude such that it increases the tidal dissipation rate, causing the orbital frequency to evolve such that the mode remains resonant. A detailed analysis of a population of heartbeat stars will determine whether resonance locking is a common phenomenon. If very few stars exhibit resonantly locked modes, this would indicate that resonance locks can be broken by some effect not accounted for here. Speculative possibilities include resonance disruption due to resonances of other modes, avoided crossings of g modes due to stellar evolution, non-linear mode coupling, or non-Keplerian orbital dynamics (e.g., precession or three-body effects).

If resonance locking does commonly occur in evolving systems, it can greatly enhance tidal dissipation rates. Moreover, it simplifies the physics of tidal dissipation, as orbital decay and spin synchronization proceed on stellar evolution time scales when resonance locking occurs. Resonance locking is not necessarily limited to eccentric binary star systems, and may operate in many astrophysical settings, including inspiraling WDs (Burkart et al., 2013), planetary moon systems (Fuller et al., 2016), and many other scenarios. We hope to investigate new possibilities in future works.

Acknowledgments

I am grateful to Susan Mullally for providing the data in Figure 1, and to the anonymous referee for a thoughtful report. I thank Rich Townsend and Zhao Guo for helpful discussions. This research was supported in part by a Lee DuBridge Fellowship at Caltech, the National Science Foundation under grants AST-1205732 and PHY-1125915, and the Gordon and Betty Moore Foundation through grant GBMF5076.

Appendix A MESA Models

Our stellar models are made using the MESA stellar evolution code (Paxton et al., 2011; Paxton et al., 2013, 2015), version 9575. Important settings include the use of convective overshoot with an exponential decline above the convective zone, and a small amount of molecular diffusivity of $D=1\,{\rm cm}^{2}{\rm s}^{-1}$. Our effective overshooting parameter of $f=0.01$ is similar to values inferred for slowly pulsating B stars (Moravveji et al., 2015, 2016) and solar-like oscillators (Deheuvels et al., 2016). Larger diffusivities smooth the compositional gradient and its contribution to the Brunt-Väisälä frequency, causing less mode trapping near the core and smoothing out the dips in Figures 2. An inlist for our models is given below.

&star_job

  pgstar_flag = .true.

/ ! end of star_job namelist

&controls

  write_pulse_data_with_profile = .true.
  pulse_data_format = ’GYRE’

  initial_mass = 1.7
  initial_z = 0.02
  use_Type2_opacities = .true.
  Zbase = 2.d-2

  cool_wind_RGB_scheme = ’Reimers’
  cool_wind_AGB_scheme = ’Blocker’
  RGB_to_AGB_wind_switch = 1d-4
  Reimers_scaling_factor = 0.5
  Blocker_scaling_factor = 0.5

  overshoot_f_above_nonburn_core = 0.01
  overshoot_f0_above_nonburn_core = 0.005
  overshoot_f_above_nonburn_shell = 0.01
  overshoot_f0_above_nonburn_shell = 0.005
  overshoot_f_below_nonburn_shell = 0.01
  overshoot_f0_below_nonburn_shell = 0.005

  overshoot_f_above_burn_h_core = 0.01
  overshoot_f0_above_burn_h_core = 0.005
  overshoot_f_above_burn_h_shell = 0.01
  overshoot_f0_above_burn_h_shell = 0.005
  overshoot_f_below_burn_h_shell = 0.01
  overshoot_f0_below_burn_h_shell = 0.005

  set_min_D_mix = .true.
  min_D_mix = 1d0

  photo_interval = 100
  profile_interval = 3
  max_num_profile_models = 3000
  history_interval = 3
  terminal_interval = 1
  write_header_frequency = 10
  max_number_backups = 500
  max_number_retries = 1000
  max_timestep = 1.15d14  ! in seconds.

  mesh_delta_coeff = 0.3
  varcontrol_target = 5.d-4

  xa_central_lower_limit_species(1) = ’h1’
  xa_central_lower_limit(1) = 0.001

/ ! end of controls namelist

After constructing stellar models, we compute non-adiabatic oscillation mode properties with the GYRE pulsation code (Townsend & Teitler, 2013), version 5.0. We employ rotation via the traditional approximation. An inlist is given below.

&model model_type = ’EVOL’ ! Obtain stellar structure from an evolutionary model file = ’profile118.data.GYRE’ ! File name of the evolutionary model file_format = ’MESA’ ! File format of the evolutionary model uniform_rot = .true. !Turn on rotation Omega_units = ’RAD_PER_SEC’ !Turn on rotation Omega_rot = 2.433e-5 /

&mode l = 2 ! Harmonic degree m = 2 /

&osc outer_bound = ’ZERO’ ! Use a zero-pressure outer mechanical boundary condition rotation_method = ’TAR’ !Use traditional approximation nonadiabatic = .TRUE. /

&num diff_scheme = ’COLLOC_GL2’ ! 2nd-order Magnus solver for initial-value integrations n_iter_max = 100 /

&scan grid_type = ’INVERSE’ ! Scan for modes using a uniform-in-period grid; best for g modes freq_min = 0.1 ! Minimum frequency to scan from freq_max = 10.0 ! Maximum frequency to scan to n_freq = 1000 ! Number of frequency points in scan grid_frame = ’COROT_O’ freq_frame = ’COROT_O’ /

&grid alpha_osc = 10 alpha_exp = 2 n_inner = 10 /

&ad_output summary_file = ’summary_ad.txt’ ! File name for summary file summary_file_format = ’TXT’ ! Format of summary file summary_item_list = ’M_star,R_star,l,n_pg,omega,E_norm,f_T,psi_T’ ! Items to appear in summary file mode_file_format = ’TXT’ ! Format of mode files mode_item_list = ’l,n_pg,omega,x,xi_r,xi_h,eul_phi,lag_L,lag_T,rho,Omega_rot’ ! Eigenfunction outputs /

&nad_output summary_file = ’summary_nad.txt’ ! File name for summary file summary_file_format = ’TXT’ ! Format of summary file summary_item_list = ’M_star,R_star,l,n_pg,omega,E_norm,f_T,psi_T’ ! Items to appear in summary file mode_template = ’mode%n.txt’ ! File-name prefix for mode files mode_file_format = ’TXT’ ! Format of mode files mode_item_list = ’l,n_pg,omega,x,xi_r,xi_h,eul_phi,lag_L,lag_T,rho,Omega_rot’ ! Eigenfunction outputs /

We implement a couple important changes when computing our stellar oscillation modes. First, we renormalize GYRE’s eigenfunctions according to equation 52. Second, we set $c_{\rm rad}^{\prime}=-3$ using the statement dc_rad = -3._WP in GYRE’s oscillation equations. The purpose of this is to eliminate terms arising from gradients in the fraction of energy carried by radiation, $c_{\rm rad}$. We find these terms become very large in the partial hydrogen convective zone near the surfaces of stars with $T_{\rm eff}\sim 7000\,{\rm K}$, due to the small scale height and sharp dependence of opacity on temperature. Consequently, the radial derivative of the luminosity perturbation eigenfunction, $\delta L/L$, becomes very large in this region. We deem this behavior unphysical, because the thermal time in these regions is much shorter than g mode pulsation frequencies. In this regime, the oscillations should be roughly isothermal (rather than adiabatic), and the luminosity perturbation should be essentially frozen in these layers, i.e., the derivative of $\delta L/L$ should be near zero.

We find that implementing the change above yields more physical eigenfunctions. The ultimate source of the problem is that GYRE does not compute changes to the convective flux perturbation, which is a well known and very difficult problem in asteroseismology (Unno et al., 1989). Our implementation effectively sets the convective luminosity perturbation equal to zero, while GYRE’s current default implementation effectively sets the convective entropy perturbation equal to zero (see equations 21.6 and 21.7 of Unno et al. 1989). Neither method is physically accurate. However, after experimenting with different approximations for convective flux/entropy perturbations, we find that many prescriptions (including the one we adopt) yield somewhat similar results, with slight differences in the surface luminosity perturbation in Figure 2. Based on these experiments, we believe our computed values of $L_{\alpha}$ are accurate at the factor of 2 level for modes with $f_{\alpha}\sim 0.5\,{\rm c/d}$ in stars with $T_{\rm eff}\sim 7000\,{\rm K}$. For higher frequency modes and stars of cooler/warmer temperatures, convective flux perturbations are less of an issue, and our values of $L_{\alpha}$ are probably significantly more reliable. GYRE’s default prescription produces very different (and in our opinion, unphysical) eigenfunctions due to the large $c_{\rm rad}^{\prime}$ term as described above. Since g modes have very little inertia in the hydrogen partial ionization region, GYRE’s oscillation mode frequencies and growth rates are likely to be insensitive to these effects. However, the surface luminosity perturbation is significantly impacted, and future work should attempt to more accurately treat convective flux perturbations in these regions.