Dynamical resonance locking in tidally interacting binary systems (original) (raw)

Abstract

We examine the dynamics of resonance locking in detached, tidally interacting binary systems. In a resonance lock, a given stellar or planetary mode is trapped in a highly resonant state for an extended period of time, during which the spin and orbital frequencies vary in concert to maintain the resonance. This phenomenon is qualitatively similar to resonance capture in planetary dynamics. We show that resonance locks can accelerate the course of tidal evolution in eccentric systems and also efficiently couple spin and orbital evolution in circular binaries. Previous analyses of resonance locking have not treated the mode amplitude as a fully dynamical variable, but rather assumed the adiabatic (i.e. Lorentzian) approximation valid only in the limit of relatively strong mode damping. We relax this approximation, analytically derive conditions under which the fixed point associated with resonance locking is stable, and further check these analytic results using numerical integrations of the coupled mode, spin, and orbital evolution equations. These show that resonance locking can sometimes take the form of complex limit cycles or even chaotic trajectories. We provide simple analytic formulae that define the binary and mode parameter regimes in which resonance locks of some kind occur (stable, limit cycle, or chaotic). We briefly discuss the astrophysical implications of our results for white dwarf and neutron star binaries as well as eccentric stellar binaries.

1 INTRODUCTION

The most frequent treatment of tidal effects in detached binaries relies on the weak-friction theory (Murray & Dermott 1999), which considers only the large-scale ‘equilibrium tide’ (i.e. filling of the Roche potential) with dissipation parametrized by the tidal quality factor |$\mathcal {Q}$| (Goldreich & Soter 1966). Such a treatment fails to account for the resonant excitation of internal stellar waves with intrinsic frequencies comparable to the tidal forcing – the so-called ‘dynamical tide’ (Zahn 1975).

Many studies have considered the dynamical tide in different astrophysical contexts. There are two possible regimes: when the character of the excited modes is that of radially travelling waves, or when they represent standing waves. A travelling wave occurs when reflection is prohibited by a strong dissipative process at some radial location in the star or planet in question – e.g. rapid linear dissipation near the surface or non-linear wave breaking.

The standing wave regime is the subject of this paper. This is applicable when a wave's amplitude can be built up by many reflections. Existing calculations in this regime have used at least one of the following two approximations: (1) tides do not backreact on the spin of the star or planet in question (Lai 1994; Rathore, Blandford & Broderick 2005; Fuller & Lai 2011), or (2) the mode amplitude is not treated as a dynamical variable, and instead has its amplitude set by the adiabatic approximation discussed in Section 6.1.1 (Witte & Savonije 1999; Fuller & Lai 2012a; Burkart et al. 2013).

In this paper, we are interested in understanding the phenomenon of resonance locking, in which the orbital and spin frequencies vary in concert so as to hold the Doppler-shifted tidal forcing frequency _k_Ωorb − _m_Ωspin constant (Witte & Savonije 1999). Resonance locking is analogous to the phenomenon of capture into resonance in planetary dynamics (Goldreich 1965; Goldreich & Peale 1968); we provide a comparison in Section 10. Resonance locks can accelerate the course of tidal evolution, as we will show in Section 8. Moreover, recent studies (Burkart et al. 2012; Fuller & Lai 2012a) have proposed that resonance locks may have been observed in the Kepler system KOI-54 (Welsh et al. 2011), although this has been contested by O'Leary & Burkart (2014).

Since resonance locking involves a changing spin frequency, clearly it cannot occur under approximation (1) noted above. The domain of validity of approximation (2) is given in Section 6.1.1 (see also Burkart et al. 2013). In this paper, we drop both of the above assumptions and examine resonance locks accounting for a dynamically evolving mode amplitude coupled to both the orbital and spin evolution. Our aim is to investigate the general dynamical properties of resonance locking, rather than to focus on a specific astrophysical application. Our key questions concern determining when resonance locks can occur and under what conditions they are dynamically stable.

This paper is structured as follows. We first describe the essential idea behind resonance locking in Section 2, and enumerate the approximations we make in order to limit the complexity of our analysis in Section 3. We then develop evolution equations for a single stellar or planetary eigenmode in Section 4.1, and determine the implied backreaction upon the binary orbit and stellar or planetary spin in Section 4.2. We establish the existence and assess the stability of fixed points in the evolution equations associated with resonance locking in Section 5. We describe two analytic approximations that a mode's amplitude follows in certain limits in Section 6.1, and then present example numerical integrations of resonance locks in Section 6.2. In Section 6.3, we discuss the possibility of chaos during resonance locking. We numerically and analytically determine the parameter regimes that lead to resonance locking in Section 7, and show that resonance locks can accelerate tidal evolution in Section 8. We apply our results to two example astrophysical systems – inspiralling compact object binaries and eccentric stellar binaries – in Section 9. We then conclude in Section 10.

2 BASIC IDEA

We first explain the essential mechanism behind resonance locking by considering the example situation of a circular white dwarf binary inspiralling due to the emission of gravitational waves (Burkart et al. 2013). Focusing on a particular white dwarf, and shifting to a frame of reference corotating with the white dwarf's spin, the tidal forcing frequency is σ = m(Ωorb − Ωspin), where m is the azimuthal spherical harmonic index and we temporarily assume Ωorb ≫ Ωspin. Due to the influence of gravitational waves, Ωorb gradually increases, and thus so does σ. As such, σ sweeps towards resonance with the nearest normal mode, and this mode gains energy as it becomes increasingly resonant (Rathore et al. 2005).

Along with energy, however, comes angular momentum (for m ≠ 0 modes). As the mode then damps, this angular momentum is transferred to the background rotation, increasing Ωspin and consequently decreasing σ. Thus if the mode is capable of achieving a sufficient amplitude, fixed points can exist where σ (but not Ωorb or Ωspin individually) is held constant by tidal synchronization balancing orbital decay by gravitational waves, as illustrated in Fig. 1. This is the idea behind a tidal resonance lock.

Figure 1.

Diagram illustrating the basic mechanism behind tidal resonance locking. Forcing frequency drift caused for example by gravitational waves causes the tidal forcing frequency σ to advance towards the right. The influence of tidal synchronization, on the other hand, becomes stronger as resonance approaches, and tends to push the forcing frequency to the left. This creates a fixed point on the resonant mode's Lorentzian amplitude profile. (The mirror image of this situation is also possible.)

The properties of such fixed points corresponding to resonance locking clearly depend on the rate of externally driven orbital evolution and the strength of tidal coupling to the mode in question. Furthermore, the mode's damping rate influences both its maximum achievable amplitude as well as the rate at which it dissipates angular momentum into the background rotation profile. Lastly, since resonance locking involves a balance between orbital and spin evolution, the ratio of the associated moments of inertia – roughly _MR_2/μ_a_2 where μ is the reduced mass and a is the semi-major axis – also plays a key role. We will see in subsequent sections that dimensionless parameters corresponding to these four quantities will entirely determine the dynamics of resonance locking.

The ideas we have presented here in the context of inspiralling white dwarf binaries carry over to the more general scenario where some generic physical process causing the forcing frequency σ or the eigenmode frequency ω to evolve in one direction – e.g. gravitational waves, magnetic braking, the equilibrium tide, stellar evolution, etc. – is balanced by the resonant tidal interaction with a planetary or stellar normal mode causing σ to evolve in the opposite direction.

3 ESSENTIAL ASSUMPTIONS

Throughout this work, we invoke the following principal assumptions.

• As discussed in Section 1, we assume that the dynamical tide is composed of standing waves, and thus that linear dissipation (provided e.g. by radiative diffusion) is not strong enough to prohibit wave reflection.
• As a mode's amplitude grows, non-linear processes can become important. For example, gravity waves break and catastrophically dissipate when their vertical displacement becomes comparable to their vertical wavelength (Lindzen 1966). Moreover, global parametric resonances can occur at smaller amplitudes, which transfer energy to a pair of daughter modes (e.g. Weinberg et al. 2012).
  We completely neglect non-linear hydrodynamical phenomena in this work. The applicability of this approximation depends on the particular star or planet's structure and the binary's parameters. For example, although Witte & Savonije (2002) considered resonance locking in solar-type binaries, Goodman & Dickson (1998) showed that the lack of a convective core in such stars allows geometric focusing to progressively amplify the local amplitude of an inward-propagating gravity wave to the point where wave breaking occurs over a wide range of binary parameters (see also Barker & Ogilvie 2010). This precludes the establishment of the standing waves required for a resonance lock.
• We assume that the star or planet in question rotates as a rigid body, even though tidal torques are typically applied to the background rotation profile in a thin spherical shell where mode damping is strongest (Goldreich & Nicholson 1989b). We thus rely upon the action of an efficient angular momentum transport process. Whether such a process is available strongly depends on the application in question; see e.g. Burkart et al. (2013) and accompanying references. Without solid-body rotation, critical layers can develop where angular momentum is deposited, which would provide rapid, local dissipation and hence violate assumption (i).
• We take the overall strength of a mode's tidal forcing to be constant, and account for orbital evolution only insofar as it affects the forcing frequency. We thus ignore, for example, that the magnitude of the tidal force depends on the binary's semi-major axis, which necessarily must evolve if Ωorb is to change. We determine when this assumption is valid in Section 5.1; it would never be appropriate for long-time-scale simulations that follow the evolution of many resonance locks (which are not performed in this paper).
• Lastly, to simplify our analysis, we only consider binary systems where the spin and orbital angular momenta are aligned. Generalization to spin-orbit-misaligned systems is straightforward (Lai 2012), and ultimately only introduces a Wigner |$\mathcal {D}$|-matrix element into our definition of the tidal coupling coefficient U in Section 4.1.

4 FORMALISM

4.1 Mode amplitude evolution

We work in a frame of reference corotating with the spin of the body in question, where the rotation axis lies along the |$\hat{z}$| direction. Following Schenk et al. (2002), we invoke a phase space eigenmode decomposition of the tidal response (Dyson & Schutz 1979). The standard stellar mode inner product for arbitrary vector fields |$\boldsymbol {\eta }$| and |$\boldsymbol {\eta }^{\prime }$| is

\begin{equation} \langle \boldsymbol {\eta },\boldsymbol {\eta }^{\prime }\rangle = \int \boldsymbol {\eta }^*\cdot \boldsymbol {\eta }^{\prime } \rho \,\mathrm{d}V, \end{equation}

(1)

and the antiHermitian Coriolis force operator is defined by |$B\,\boldsymbol {\eta } = 2\Omega _\mathrm{spin}\hat{z}\times \boldsymbol {\eta }$|⁠.

We consider the resonant interaction between a time-varying tidal potential evaluated on a Keplerian orbit and a single linear mode of our expansion. The mode has spherical harmonic indices l and m, and has complex amplitude q. The differential equation describing the mode's linear evolution is (Schenk et al. 2002)

\begin{equation} \dot{q} + (i\omega +\gamma ) q = \frac{i\omega }{\varepsilon }\left\langle \boldsymbol {\xi }, \boldsymbol {a}_\mathrm{tide}\right\rangle , \end{equation}

(2)

where ω is the mode's eigenfrequency, γ is its damping rate, |$\boldsymbol {\xi }$| is its Lagrangian displacement vector,

\begin{equation} \varepsilon = 2\omega ^2\langle \boldsymbol {\xi },\boldsymbol {\xi }\rangle + \omega \langle \boldsymbol {\xi },iB\boldsymbol {\xi }\rangle \end{equation}

(3)

is a normalization factor (equal to the mode energy at unit amplitude), and |$\boldsymbol {a}_\mathrm{tide}$| is the time-dependent tidal acceleration vector. We assume γ ≪ ω and, without loss of generality, we take ω > 0.

Allowing for an arbitrary eccentricity, the projection of the tidal acceleration on to the mode |$\left\langle \boldsymbol {\xi }, \boldsymbol {a}_\mathrm{tide}\right\rangle$| is proportional to

\begin{equation*} \left(\frac{a}{D}\right)^{l+1}\text{e}^{-\text{i}m\left(f-\psi _\mathrm{spin}\right)}, \end{equation*}

where D is the binary separation, a is the semi-major axis, f is the true anomaly and

\begin{equation} {\left(\begin{array}{c}\psi _{\rm orb}\\ \psi _{\rm spin}\end{array}\right)} = \int {\left(\begin{array}{c}\Omega _{\rm orb}\\ \Omega _{\rm spin}\end{array}\right)} {\rm d}t. \end{equation}

(4)

Assuming that changes in the orbital frequency and eccentricity occur on time-scales much longer than an orbital period, we can expand the dependence on D and f in a Fourier series:

\begin{equation} \left(\frac{a}{D}\right)^{l+1}\text{e}^{-\text{i}m\left(f-\psi _\mathrm{spin}\right)} \approx \sum _k X_{lm}^k \text{e}^{-\text{i}\left(k\psi _\mathrm{orb}-m\psi _\mathrm{spin}\right)}, \end{equation}

(5)

where |$X_{lm}^k(\text{e})$| is a Hansen coefficient (Appendix B), with |$X_{lm}^k(\text{e}=0)=\delta ^k_{m}$| for circular orbits. Since we are concerned with resonant mode-tide interaction, we henceforth consider only a single harmonic component of this expansion. The phase associated with this component is ψ = _k_ψorb − _m_ψspin.

Our mode amplitude equation from (2) then becomes

\begin{equation} \dot{q} + (i\omega +\gamma )q = i\omega U \text{e}^{-\text{i}\psi }. \end{equation}

(6)

The dimensionless mode-tide coupling strength associated with our harmonic component is

\begin{equation} U = \left( \frac{M^{\prime }}{M} \right)\left( \frac{R}{a} \right)^{l+1}\left( \frac{E_*}{\varepsilon } \right)\mathcal {I}_{lm}\, X_{lm}^k\, W_{lm}, \end{equation}

(7)

where _M_′ is the companion mass,

\begin{equation} \mathcal {I}_{lm} = \frac{1}{MR^l}\left\langle \boldsymbol {\xi },\boldsymbol {\nabla }\left(r^l Y_{lm}\right)\right\rangle \end{equation}

(8)

is the mode's linear tidal overlap integral (Press & Teukolsky 1977),

\begin{equation} W_{lm} = \frac{4\pi }{2l+1}Y_{lm}^*\left( \frac{\pi }{2},\,0 \right) \end{equation}

(9)

is an order-unity angular coupling coefficient, and E* = _GM_2/R is the gravitational energy scale. We also define the tidal forcing frequency to be |$\sigma = \dot{\psi } = k\Omega _\mathrm{orb}- m\Omega _\mathrm{spin}$|⁠, and take σ > 0 without loss of generality. We are free to choose the sign of σ since we are considering a complex conjugate mode pair – one member of the pair is resonant with σ while the other is resonant with −σ. Because m can possess either sign, both prograde and retrograde modes are allowed for σ > 0. Since we are considering resonant interaction, our assumption is σ ≈ ω.

4.2 Forcing frequency evolution

We will first obtain the time derivative of the stellar or planetary spin frequency implied by equation (6). The canonical angular momentum associated with the mode together with its complex conjugate is given by (Appendix C)

\begin{equation} J_\mathrm{mode}= \frac{m\varepsilon }{\omega }|q|^2. \end{equation}

(10)

Differentiating with respect to time and substituting equation (6), we find

\begin{equation} \dot{J}_\mathrm{mode}= -2\gamma J_\mathrm{mode}+ 2m\varepsilon U \mathrm{Im}\left[ qe^{i\psi } \right]. \end{equation}

(11)

The first term in equation (11) results from damping, and is imparted to the stellar or planetary spin (Goldreich & Nicholson 1989a).1 We thus set

\begin{equation} \dot{J}_\mathrm{spin}= 2\gamma J_\mathrm{mode}, \end{equation}

(12)

meaning that

\begin{equation} \dot{\Omega }_\mathrm{spin}= \frac{2\gamma J_\mathrm{mode}}{I_*} + \alpha _\mathrm{spin}, \end{equation}

(13)

where I* is the planet or star's moment of inertia, and we have incorporated an additional term αspin to account for processes that can change Ωspin other than interaction with the mode in question – e.g. the equilibrium tide, magnetic braking, etc. We take αspin to be constant.

The rate at which the orbital energy changes is given by (Weinberg et al. 2012)

\begin{equation} \dot{E}_\mathrm{orb}= -2k\Omega _\mathrm{orb}\varepsilon U \mathrm{Im}\left[ q\text{e}^{\text{i}\psi } \right]. \end{equation}

(14)

We can convert equation (14) into an expression for |$\dot{\Omega }_\mathrm{orb}$| using

\begin{equation} \frac{\dot{\Omega }_\mathrm{orb}}{\Omega _\mathrm{orb}} = -3\,\frac{\dot{E}_\mathrm{orb}}{\mu a^2 \Omega _\mathrm{orb}^2 }, \end{equation}

(15)

which yields

\begin{equation} \dot{\Omega }_\mathrm{orb}= \left(\frac{3k}{\mu a^2}\right) 2\varepsilon U\mathrm{Im}\left[ q\text{e}^{\text{i}\psi } \right] + \alpha _\mathrm{orb}, \end{equation}

(16)

where we have again included an extra term αorb to account for e.g. orbital decay by gravitational waves.

It is useful to combine equations (13) and (16) to determine the time derivative of δω = ω − σ, which can be expressed as

\begin{equation} \frac{\dot{\delta \omega }}{\omega } = -\Gamma _{\mathrm{d}r}+ \Gamma _\mathrm{tide}\left( \frac{\gamma |q|^2}{\omega U^2} -r\, \mathrm{Im}\left[\frac{q\text{e}^{\text{i}\psi }}{U}\right] \right). \end{equation}

(17)

Here, we have combined both α parameters from earlier into the ‘drift’ rate

\begin{equation} \Gamma _{\mathrm{d}r}= \frac{k\alpha _\mathrm{orb}-m(1-C)\alpha _\mathrm{spin}}{\omega } - \frac{\mathrm{\partial} \ln \omega }{\mathrm{\partial} t}, \end{equation}

(18)

where

\begin{equation} C = -\frac{1}{m}\frac{\mathrm{\partial} \omega }{\mathrm{\partial} \Omega _\mathrm{spin}} \end{equation}

(19)

allows for a rotationally-dependent corotating-frame eigenmode frequency and ∂ω/∂t accounts for changes in the eigenmode frequency due to progressive modifications of the background hydrostatic profile from e.g. stellar evolution.2 The tidal backreaction rate is

\begin{equation} \Gamma _\mathrm{tide}= \frac{2m^2U^2 (1-C)\varepsilon }{I_*\omega }, \end{equation}

(20)

and parametrizes the strength of tidal coupling to the mode in question; it is related to the rate at which the mode can synchronize the star or planet at non-resonant amplitudes (|q| ∼ |U|). Lastly, the moment of inertia ratio is

\begin{equation} r = \frac{k^2}{m^2}\frac{3I_*}{(1-C)\mu a^2}. \end{equation}

(21)

We assume that both Γtide and r are positive throughout this work, but allow Γd_r_ to possess either sign.

5 RESONANCE LOCK FIXED POINTS

5.1 Existence of fixed points

We first remove the overall oscillatory time dependence of q by changing variables to |$Q=q\text{e}^{\text{i}\psi }/U$|⁠, so that equation (6) becomes

\begin{equation} \dot{Q} + \left( i\delta \omega +\gamma +\frac{\dot{U}}{U}\right)Q=i\omega , \end{equation}

(22)

where again δω = ω − σ. For simplicity, we now assume that γ is much larger than terms contributing to |$\dot{U}/U$|⁠, such as |$\dot{a}/a$| and |$\dot{e}/e$|⁠, so that equation (22) becomes

\begin{equation} \dot{Q} + \left( i\delta \omega +\gamma \right)Q=i\omega ; \end{equation}

(23)

as we will see in Section 8, this essentially amounts to assuming γ ≫ |Γd_r_|. Equation (17) becomes

\begin{equation} \frac{\dot{\delta \omega }}{\omega } = -\Gamma _{\mathrm{d}r}+ \Gamma _\mathrm{tide}\left( \frac{\gamma }{\omega }|Q|^2 - r\,\mathrm{Im}[Q] \right). \end{equation}

(24)

Having thus eliminated direct dependence on the phase ψ, our two dynamical variables are now Q and δω. Since Q is complex, we have a third-order differential system.

Resonance locking corresponds to a fixed point in the evolution equations. We thus set time derivatives to zero in equation (23) to derive

\begin{equation} Q_{\mathrm{f}}= \frac{\omega }{\delta \omega _{\mathrm{f}}- i\gamma }, \end{equation}

(25)

and hence

\begin{equation} {\left[\begin{array}{c}\mathrm{Re}(Q_{\mathrm{f}})\\ \mathrm{Im}(Q_{\mathrm{f}})\end{array}\right]} = \frac{\omega }{\delta \omega _{\mathrm{f}}^2 + \gamma ^2}\times {\left[\begin{array}{c}\delta \omega _{\mathrm{f}}\\ \gamma \end{array}\right]}. \end{equation}

(26)

Similarly, setting |$\dot{\delta \omega }=0$| in equation (24) and using the fact that (γ/ω)|_Q_f|2 = Im(_Q_f), we have

\begin{equation} \Gamma _{\mathrm{d}r}= (1-r)\,\Gamma _\mathrm{tide}\,\mathrm{Im}(Q_{\mathrm{f}}). \end{equation}

(27)

We can use equations (26) and (27) to derive

\begin{equation} \delta\omega _{\mathrm{f}}^2 = \frac{\omega\gamma \Gamma _\mathrm{tide}(1-r)}{\Gamma _{\mathrm{d}r}}-\gamma^2. \end{equation}

(28)

So far we have not determined the sign of δωf, and indeed there is one fixed point for each sign; however, we will show in Section 5.2.2 that one is always unstable.

These fixed points exists if equations (26) and (27) can be solved for _Q_f and σf, which is possible if (assuming Γtide > 0)

\begin{equation} \frac{\gamma}{\omega \Gamma _\mathrm{tide}} < \frac{1-r}{\Gamma _{\mathrm{d}r}}; \end{equation}

(29)

this in particular requires

\begin{equation} (1-r)\Gamma _{\mathrm{d}r}>0. \end{equation}

(30)

If Γd_r_ > 0, then equation (30) reduces to

\begin{equation} |k|<|m|\left(\frac{(1-C)\mu a^2}{3I_*}\right)^{1/2}. \end{equation}

(31)

As a result, Fuller & Lai (2012a) referred to the quantity on the right-hand side of equation (31) as the ‘critical’ harmonic.

The existence of the fixed points also relies on the ‘weak-tide’ limit, which we define to be |δωf| < Δω/2, where Δω is the eigenmode frequency spacing near the mode in question. This means that the fixed point must lie within the eigenmode's domain of influence, so that the contribution of other eigenmode resonances can be legitimately neglected. If this is not the case, then resonance locking is not possible.

For example, if Γd_r_ is very small, equation (27) requires a commensurately small value of Q. However, this could be impossible to achieve in practice, since it would require a very large value of the detuning δω, allowing the possibility for a neighbouring mode to come into resonance. The actual outcome in such a situation is that the tidal interaction would dominate the dynamics, and the drift processes contributing to Γd_r_ would be irrelevant – i.e. the ‘strong-tide’ limit.

A necessary condition for the weak-tide limit, and thus for a resonance lock to be able to occur, is |δωf| ≪ ω, which evaluates to

\begin{equation} \frac{|\Gamma _{\mathrm{d}r}|}{\gamma } \gg \frac{\Gamma _\mathrm{tide}|1-r|}{\omega }. \end{equation}

(32)

We will use this as a convenient, although very liberal, proxy for the real requirement of |δωf| < Δω/2, so as to avoid handling the additional parameter Δω. Calculations of dynamical tidal evolution that use the adiabatic approximation (Section 6.1.1) with many modes – e.g. Witte & Savonije (1999); Fuller & Lai (2012a) – already naturally account for the true requirement.

5.2 Fixed point stability

5.2.1 Linearization

We will now perform a linear stability analysis about each fixed point. First, note that the presence of the non-analytic functions | · | and Im(·) in equation (24) necessitates treating the real and imaginary parts of Q separately. Thus let

\begin{equation*} \boldsymbol {\zeta } ={\left[\begin{array}{c}\mathrm{Re}(Q)\quad \mathrm{Im}(Q)\quad (\sigma -\omega )/\omega \end{array}\right]}^T, \end{equation*}

and

\begin{equation*} \frac{\text{d}\boldsymbol {\zeta }}{\text{d}t} = \boldsymbol {f}(\boldsymbol {\zeta }), \end{equation*}

where |$\boldsymbol {f}$| represents equations (23) and (24) and satisfies |$\boldsymbol {\nabla } \cdot \boldsymbol {f} = -2\gamma$|⁠. We can then derive a time-evolution equation for |$\boldsymbol {\delta \zeta } = \boldsymbol {\zeta } - \boldsymbol {\zeta }_{\mathrm{f}}$| by invoking results from Section 5.1.

Proceeding, to linear order we have

\begin{equation} \frac{\text{d}}{\text{d}t}\boldsymbol {\delta \zeta } = A \boldsymbol {\delta \zeta } + o(\boldsymbol {\delta \zeta }), \end{equation}

(33)

where

\begin{eqnarray} A=\mathcal {D}\boldsymbol {f}(\boldsymbol {\zeta }_{\mathrm{f}})={\left[\begin{array}{ccc}-\gamma &\quad \delta \omega _{\mathrm{f}}&\quad A_{1,3} \\ -\delta \omega _{\mathrm{f}}&\quad -\gamma &\quad A_{2,3} \\ A_{3,1}&\quad A_{3,2}&\quad 0 \end{array}\right]}, \end{eqnarray}

(34)

\begin{eqnarray} A_{1,3} = -\frac{\omega \beta }{\Gamma _\mathrm{tide}} \quad A_{2,3} = \frac{\omega \beta \delta \omega _{\mathrm{f}}}{\gamma \Gamma _\mathrm{tide}} \end{eqnarray}

(35)

\begin{eqnarray} A_{3,1} = -\frac{2\beta \delta \omega _{\mathrm{f}}}{\omega }\quad A_{3,2} = r\Gamma _\mathrm{tide}- \frac{2 \beta \gamma }{\omega }, \end{eqnarray}

(36)

and β = Γd_r_/(1 − r).

5.2.2 Eigenvalues

We assume in this section that equation (29) is satisfied, so that the resonance locking fixed point formally exists, and that equa-tion (32) is also satisfied, so that we are in the weak-tide limit and the fixed point physically exists. The characteristic polynomial P for the eigenvalues λ of A is

\begin{equation} P(\lambda )=\lambda ^3 + 2\gamma \lambda ^2+P_1\lambda +P_0, \end{equation}

(37)

where

\begin{equation} P_1 = \frac{\omega \gamma \Gamma _\mathrm{tide}}{\beta } - \frac{r\omega \beta \delta \omega _{\mathrm{f}}}{\gamma } \qquad P_0 = 2\omega \Gamma _{\mathrm{d}r}\delta \omega _{\mathrm{f}}. \end{equation}

(38)

A fixed point is asymptotically stable if all eigenvalues satisfy Re(λ) < 0. A standard theorem then states that, for this to occur, it is necessary (but not sufficient) that all of the coefficients of λ_i_ (_i_ ≥ 0) in equation (37) possess the same sign. We thus must have that _P_1, 0 > 0.

First, we see that

\begin{equation} +\!1=\mathrm{sign}(P_0)=\mathrm{sign}(\Gamma _{\mathrm{d}r})\,\mathrm{sign}(\delta \omega _{\mathrm{f}}). \end{equation}

(39)

Since equations (26) and (27) admit two solutions, corresponding to δωf = ±|δωf|, the criterion that _P_0 > 0 simply allows us to select the solution that could potentially be stable. Henceforth we will focus on the solution that satisfies equation (39), which we refer to as the ‘lagging’ fixed point; similarly, the ‘leading’ fixed point is the one that fails to satisfy equation (39).

Equation (39) permits the following intuitive interpretation: the forcing frequency σ is ‘pushed’ towards an eigenfrequency in the direction specified by the sign of Γd_r_, but the approaching eigenmode ‘pushes’ σ in the opposite direction, and resonance locking occurs when these ‘forces’ cancel out (Section 2). In particular, if Γd_r_ > 0, σf should be smaller than ω, meaning δωf = ω − σf > 0, consistent with equation (39).

It remains to analyse _P_1 and to ascertain when it is also positive. Moreover, since positivity of the characteristic polynomial's coefficients is only a necessary condition for stability, we must then further examine the Hurwitz matrix associated with P to establish when its leading principal minors are also positive (e.g. Gradshteyn & Ryzhik 2007). We perform this analysis in Appendix A in the limit of γ ≪ |δωf|; the result is that the lagging fixed point is stable if and only if

\begin{equation} \Gamma _{\mathrm{d}r}<0 \quad \text{or} \quad 0<\frac{\Gamma _{\mathrm{d}r}}{\gamma }<(1-r)\left(\frac{\Gamma _\mathrm{tide}}{\omega }\right)^{1/3} \end{equation}

(40)

(subject to the assumptions we have made thus far). In particular, the lagging fixed point is thus always stable for r > 1 per equation (30).

Fig. 2 shows the stability region of the lagging fixed point as a function of the damping rate γ and the frequency drift rate Γd_r_ for two example values of the backreaction rate Γtide and the moment of inertia ratio r. Stability was determined by numerically solving for the eigenvalues λ using equation (37), and is indicated by green shading, while instability is white; regions where the fixed point does not exist are shaded dark grey. Equations (29), (32) and (40), which are displayed as blue dashed, black dotted and magenta dot–dashed lines, closely correspond to the green region's boundaries. The values used in the top panel were chosen based on our white dwarf binary application in Section 9.1. In the bottom panel, we have r > 1, which by equation (30) implies Γd_r_ < 0.

Figure 2.

Resonance locking fixed point stability analysis. Regions where the fixed point is asymptotically stable are shaded green, while unstable regions are white, and regions where the fixed point does not exist are dark grey. We determined the stability regions by numerically evaluating the eigenvalues λ_i_ of the matrix A determined by linearizing the equations of motion about the fixed point (equation 34) and enforcing Re(λ_i_) < 0 together with equation (32). The analytic results in equations (29), (32) and (40) are displayed as blue dashed, black dotted and magenta dot–dashed lines. The lower-right grey triangle in both panels corresponds to where the weak-tide limit is certainly violated, and thus the fixed point does not physically exist; see Section 5.1. In the bottom panel, we have r_ > 1, which by equation (30) implies Γd_r < 0.

The instability boundary defined by equation (40) is in fact a supercritical Hopf bifurcation (Wiggins 2003), corresponding to the loss of stability of a complex conjugate eigenvalue pair. Past the bifurcation, inside the unstable region of parameter space, this unstable complex conjugate pair splits into two unstable real eigenvalues, as shown in Fig. 3. In other words, the fixed point switches from being an unstable spiral to an unstable node. Determining the parameter space manifold on which this splitting occurs will be useful in Section 7.2; we can accomplish this by setting the discriminant of equation (37) to zero (thus requiring a repeated root), yielding

\begin{equation} 0=36\gamma P_1 P_0 - 32 \gamma ^3 P_0+4\gamma ^2 P_1^2 - 4 P_1^3-27P_0^2. \end{equation}

(41)

Although this equation cannot be solved analytically, we can nonetheless determine the asymptotic dependence of Γd_r_ on γ in the limit γ → 0. Newton's polygon for equation (41) shows that this asymptotic dependence is linear: Γd_r_ = _A_γ. Substituting this into equation (41), setting γ = 0, and solving for A, we have that the lagging fixed point is a spiral if

\begin{equation} \Gamma _{\mathrm{d}r}< \gamma \left(\frac{1-r}{r^{2/3}}\right)\left( \frac{\Gamma _\mathrm{tide}}{\omega } \right)^{1/3}. \end{equation}

(42)

Figure 3.

Real parts of the lagging fixed point's three eigenvalues as functions of Γd_r_ for γ/ω = 10−8, Γtide/ω = 10−10 and r = 0.01 (see the left middle panel of Fig. 10). A complex conjugate pair exists for Γd_r_/ω ≲ 10−10 (equation 42), which loses stability at Γd_r_/ω ≈ 10−11.3 in a Hopf bifurcation (equation 40).

6 TRAJECTORIES

Here, we will show several examples of trajectories that can be produced by our dynamical equations from Section 5.1. First, however, in Section 6.1 we will discuss two different analytic approximations that the trajectories follow in certain limits.

6.1 Analytic approximations

6.1.1 Adiabatic approximation

In this paper, we define the adiabatic approximation to be the situation where the mode amplitude can instantaneously adjust to a changing forcing frequency, and thus we can set the |$\dot{Q}$| term in equation (23) to zero and assume σ is constant. Equation (23) can then be solved exactly:

\begin{equation} Q_\mathrm{ad}= \frac{\omega }{\delta \omega -i\gamma }. \end{equation}

(43)

This approximation, also referred to as the Lorentzian approximation due to the form of equation (43), is frequently employed in the literature.

The domain of validity of the adiabatic approximation can be determined by comparing the maximum possible mode growth rate to the growth rate implied by equation (43); when the latter exceeds the former, the adiabatic approximation is no longer valid. We can determine the maximum rate at which a mode amplitude Q can grow by providing a perfect resonance to equation (23), i.e. by setting σ = ω. Dropping the damping term and setting |$\dot{\sigma }=0$|⁠, the particular solution is

\begin{equation} Q(t) =i \omega t, \end{equation}

(44)

which implies that |$|\dot{Q}|_{\mathrm{m}ax}=\omega$|⁠. Next, in order to estimate the time derivative of _Q_ad, we take |$\dot{\sigma } \approx \Gamma _{\mathrm{d}r}\omega$|⁠, which gives us

\begin{equation} |\dot{Q}_\mathrm{ad}| \approx |\Gamma _{\mathrm{d}r}|\left( \frac{\omega ^2}{\delta \omega ^2+\gamma ^2}\right). \end{equation}

(45)

We then equate |$|\dot{Q}_\mathrm{ad}|$| and |$|\dot{Q}|_{\mathrm{m}ax}$| and solve for δω, finding

\begin{equation} \delta \omega _\mathrm{ad}^2 \approx |\Gamma _{\mathrm{d}r}|\omega -\gamma ^2; \end{equation}

(46)

if |Γd_r_|ω < γ2, then the adiabatic approximation is always valid.

6.1.2 No-backreaction approximation

Although directly solving equations (23) and (24) outside the adiabatic limit requires numerical integration, we can nonetheless produce an approximate analytic expression for Q(t) in the limit that backreaction of the mode on the tidal forcing frequency σ is unimportant (Reisenegger & Goldreich 1994; Rathore et al. 2005). This approximation subsumes the adiabatic approximation, but is also more complicated.

Since we are already assuming that the mode damping rate is weak enough for the resonance locking fixed point to exist (equa-tion 29), we can simply take γ → 0. Subject to this simplification, we can solve equation (6), yielding

\begin{equation} q(t) \approx i\omega U\text{e}^{-\text{i}\omega t}\int _{t_0}^t \text{e}^{\text{i}\omega t-\text{i}\psi }\text{d}t, \end{equation}

(47)

with t_0 ≪ −(ωΓd_r)−1/2. If we then approximate ψ as

\begin{equation*} \psi \approx \psi _0+\omega t + \omega \Gamma _{\mathrm{d}r}t^2/2, \end{equation*}

where we have set resonance to occur (i.e. |$\dot{\psi }=\omega$|⁠) at t = 0, then equation (47) becomes a closed-form solution to the equations of motion.

Since the integral in equation (47) approaches a constant for t ≫ (ωΓd_r_)−1/2, and since we wish to estimate the maximum value of |Q|, we can simply extend the domain of integration to (−∞, +∞), yielding (using e.g. the method of stationary phase)

\begin{equation} q(t) \approx -(1-i)U\sqrt{\frac{\pi \omega }{\Gamma _{\mathrm{d}r}}}\text{e}^{-\text{i}\omega t-\text{i}\psi _0} \end{equation}

(48)

for t ≫ (ωΓd_r_)−1/2. We thus find that

\begin{equation} |Q|_{\mathrm{m}ax}\approx \sqrt{\frac{2\pi \omega }{|\Gamma _{\mathrm{d}r}|}}. \end{equation}

(49)

We plot this maximal value of |Q| in the middle panels of Figs 4–7 as a dash–dotted black line.

Figure 4.

Numerical integration of mode and orbital evolution equations for the case of resonance locking into a stable fixed point. Parameters used were Γd_r_/ω = 10−9, γ/ω = 10−4.5, Γtide/ω = 10−10 and r = 0.5. The adiabatic approximation is shown in the top panel as a dashed black line (Section 6.1.1), while the actual system trajectory is purple. The red circle shows the lagging fixed point (Section 5). Individual time series for the mode amplitude Q and the forcing frequency σ are shown in the bottom two panels. The dash–dotted black line in the mode amplitude panel shows equation (49), which gives the maximum amplitude attainable under the no-backreaction approximation (Section 6.1.2). The lagging fixed point's eigenvalues are λ/ω ∈ {(−0.031 ± 1.2 i) × 10−3, −1.6 × 10−6}.

Figure 5.

Resonance lock limit cycle for a case in which the fixed points exist (red and black circles) but are linearly unstable. Parameters used were Γd_r_/ω = 10−7.7, γ/ω = 10−4.5, Γtide/ω = 10−10 and r = 0.5. Conventions used are the same as in Fig. 4. Colour shows time, ranging from purple at t = 0 to light blue. The red & black circles show the lagging and leading fixed points, respectively (Section 5). System sweeps through resonance without initially being captured. However, the system then oscillates through resonance several times, pumping up the mode's amplitude. Eventually, the oscillation ceases and the mode's angular momentum discharges into rotation causing the system to travel back away from resonance and start over. The lagging fixed point's eigenvalues are λ/ω ∈ {−3.9 × 10−4, (1.6 ± 0.49 i) × 10−4}.

Figure 6.

Failed resonance lock. Parameters used were Γd_r_/ω = 10−7, γ/ω = 10−4.5, Γtide/ω = 10−10 and r = 0.5. Conventions used are the same as in Fig. 4. Colour shows time, ranging from purple at t = 0 to light blue. The unstable, repulsive fixed points (red and black circles) suppress the maximum attainable mode amplitude below the dash-dotted black line in the second panel, which shows the maximum amplitude that would be achieved without backreaction (Section 6.1.2). This suppression is severe enough to prevent resonance locking from occurring. The lagging fixed point's eigenvalues are λ/ω ∈ {−6.7 × 10−4, 5.4 × 10−4, 6.8 × 10−5}.

Figure 7.

Chaotic resonance lock around an unstable fixed point (red circle). Parameters used were Γd_r_/ω = 10−7.4, γ/ω = 10−4, Γtide/ω = 10−10 and r = 0.5. Conventions used are the same as in Fig. 4. Colour shows time, ranging from purple at t = 0 to light blue. This situation is similar to that depicted in Fig. 5, but in this case backreaction is so significant that the system never fully returns to the adiabatic approximation once deviating from it near the fixed point. Chaotic orbits are instead executed around the fixed point. The lagging fixed point's eigenvalues are λ/ω ∈ {(0.98 ± 2.4 i) × 10−4, −4.0 × 10−4}.

6.2 Numerical results

Figs 4–7 show full numerical solutions to equations (23) and (24) for several different choices of our four parameters γ, Γd_r_, Γtide and r. In each case, the mode amplitude is initialized to the adiabatic approximation in the regime where it should be valid (Section 6.1.1).

First, Fig. 4 gives a simple example of resonance locking into a stable fixed point. Next, we hold Γtide, γ and r constant, but increase Γd_r_ by a factor of ∼10, thus making σ sweep towards resonance more quickly. In this case the fixed point is no longer stable, since equation (40) is no longer satisfied, and Fig. 5 shows the limit cycle resonance lock that then occurs. The system passes through resonance in between the two unstable fixed points, and then oscillates back and forth through resonance. This oscillation pumps the mode amplitude up as the system is repelled by the fixed points. Eventually, the oscillation ceases and the mode's angular momentum discharges into rotation causing the system to travel back away from resonance and decay back on to the adiabatic approximation. The cycle then begins again. This limit cycle is in fact precisely the stable periodic orbit generated by the supercritical Hopf bifurcation (Section 5.2.2).

Fig. 6 shows the resulting evolution again holding all parameters constant except for Γd_r_, which we increase by another factor of ∼10. The resonance locking fixed point is now sufficiently unstable that it suppresses the mode amplitude's growth and prevents resonance locking from occurring. Near resonance the mode amplitude still grows appreciably, but after the lock fails to hold, damping causes the mode amplitude to decay exponentially. We have found the linear character of the lagging fixed point to be the principal distinguishing factor between a limit cycle occurring and a failed resonance capture due to fixed point suppression. Specifically, we find that if the lagging fixed point is an unstable node, meaning all its eigenvalues are real (Fig. 3), then it acts to suppress the mode amplitude and leads to a failed capture. If instead the lagging fixed point is an unstable spiral, meaning its eigenvalues contain an unstable complex conjugate pair (Section 5.2.2), then it is able to ‘pump’ a trajectory to high amplitude and allow the formation of a limit cycle. This distinction, which is valid within a factor of ∼3 in Γd_r_, will be critical in Section 7.2.

6.3 Chaos

Lastly, Fig. 7 shows a chaotic trajectory. This situation is very similar to the limit cycle evolution shown in Fig. 5, in that the resonance locking fixed point is unstable but Γd_r_ is not so large that resonance locking does not occur altogether; the essential difference is that the mode amplitude profile resulting from the adiabatic approximation (dashed line) is capped by damping not far above the fixed points, unlike in Fig. 5 where it ascends much higher. This corresponds to the fact that the choice of parameters for Fig. 7 lies close (logarithmically speaking) to the bifurcation manifold in parameter space where equation (29) ceases to be satisfied and the fixed points no longer exist. This appears to be a key ingredient for chaos, as we will explain below, which is why we have changed γ to a larger value than that used in Figs 4–6 (so that equation 29 is closer to not being satisfied).

Since the fixed points are so close to the peak of the adiabatic profile, the pumping process that occurs due to repulsion from the fixed points cannot allow the mode to acquire a very large amplitude. As a result, when the mode's angular momentum eventually drains into the background spin, the system does not retreat back from resonance very far, and has little time to decay back on to the adiabatic approximation before resonance is reached again. The initial condition upon entering resonance is consequently somewhat different each cycle, leading to the potential for chaos.

We show a three-dimensional projection of the orbit from Fig. 7 in Fig. 8. The lagging fixed point is shown by a small red sphere, while its unstable plane corresponding to eigenvalues λ/ω = (0.98 ± 2.4 i) × 10−4 is also displayed.

Figure 8.

Three-dimensional projection of the integration from Fig. 7. Time is indicated by colour, ranging from dark purple at early times to light blue. The unstable plane corresponding to eigenvalues (0.98 ± 2.4 i) × 10−4 is shown, centred on the fixed point (red sphere). Each cycle, the system begins near the fixed point, but is then ejected along the unstable manifold. Non-linear terms cause the system to decay back on to the adiabatic solution; this motion comprises the spiral structure on the left. The adiabatic solution then transports the system near to the fixed point, and the cycle begins again with perturbed initial conditions.

We now present numerical evidence that the path depicted in Fig. 8 follows a strange attractor. We emphasize that our evidence is not rigorous. A strange attractor of a dynamical system, also known as an attracting chaotic invariant set, is a set that (Wiggins 2003):

• is compact,
• is invariant under the dynamical equations,
• is attracting,
• has sensitive dependence on initial conditions, and
• is topologically transitive.

Fig. 8 appears to begin to trace out a bounded, attracting, invariant set, which we denote Λ; this addresses conditions (i)–(iii). In Fig. 9, we estimate the largest Lyapunov exponent of the trajectory in Fig. 8 to be 6 × 10−5ω. Since this is positive, trajectories that begin together deviate exponentially as time passes. This then points towards condition (iv) being satisfied. Lastly, since colour indicates time in Fig. 8, the fact that dark purple (early times) is tightly and randomly interwound with light blue (late times) leads us to believe that condition (v) is likely also satisfied. Again, we have presented only suggestive evidence; further study is required to fully address the presence of a strange attractor.

Figure 9.

Estimation of largest Lyapunov exponent for the chaotic trajectory in Figs 7 and 8. We initialized two numerical integrations of our dynamical equations on the adiabatic solution (Section 6.1.1) with slightly perturbed initial detuning frequencies: δω0 = 0.95δωf and 0.950 001δωf. The blue curve shows the norm of the difference between the resulting values of the reduced mode amplitude Q as a function of time. The trajectories initially deviate exponentially, demonstrating chaos, with a rough functional form of e t/τ for τ ≈ 2 × 104/ω (dashed magenta line). We thus estimate the largest Lyapunov exponent (Wiggins 2003) for the trajectories to be ≈1/τ ≈ 6 × 10−5ω > 0, which is close to the damping rate of γ = 10−4ω.

In addition, we note that our preliminary investigations show that the chaos present results from the Hopf orbit undergoing a period-doubling cascade, and is similar in several ways to the Rössler attractor (Rössler 1976). This requires further study.

7 ACHIEVING RESONANCE LOCKS

7.1 Numerical results

In order to ascertain the conditions that lead to resonance locks, we performed sample integrations of equations (23) and (24) numerically. We initialized each integration with

\begin{equation*} \delta \omega _0 = 10\times \max \left(\delta \omega _{\mathrm{f}},\delta \omega _\mathrm{ad}, \gamma \right) \end{equation*}

and _Q_0 set by the adiabatic approximation from equation (43). We then performed each integration until t_1 = 2δω0/ωΓd_r. We determined that a resonance lock occurred if, assuming r < 1,

\begin{equation*} \frac{\min \delta \omega }{\delta \omega _0} > -0.5, \end{equation*}

where min δω is the minimum value of δω attained over the final 10 per cent of integration.3 A similar formula was used for r > 1, but accounting for the fact that resonance is then approached from the left (in terms of δω; see Section 5.1). Note that these conditions account for stable, limit cycle, and chaotic forms of resonance locking (Section 6.2).

Fig. 10 shows our results. Light blue regions indicate that a resonance lock did occur, while white indicates the reverse, and dark grey indicates that the fixed point does not exist (Section 5.1). The green lines are the analytic formula from equation (54) below for the boundary between successful and failed resonance locking. The thin purple lines are the equivalent condition for resonance locks to occur when r > 1, from equation (53). Both analytic approximations, which we will develop in the next section, show good agreement with numerical results.

Figure 10.

Resonance lock regimes based on numerical solution of mode, spin and orbit evolution equations from Section 5.1. Light blue shading indicates that resonance locking occurs (including stable fixed point locks, limit cycles and chaotic locks, as in Figs 4, 5 and 7), while white indicates the reverse (as in Fig. 6). Above each dashed blue line, the resonance locking fixed point does not exist (Section 5.1; equation 29). The green and purple lines correspond to our analytic formulae for the resonance locking regime (Section 7.2; equations 54 and 53, respectively). The dash–dotted magenta line indicates the upper boundary of the domain of stability of the resonance lock fixed point (Section 5.2; equation 40), while the dotted black line shows where the weak-tide limit is certainly violated (Section 5.1; equation 32). Limit cycles and chaotic orbits (as in Figs 5 and 7) preferentially occur in the regions between the green and magenta lines (where the fixed point is unstable but locks still occur).

7.2 Analytic approximations

We now seek to obtain an analytical understanding of our numerical results in Fig. 10. We will thus attempt to assemble a set of analytic approximations to determine what values of our four parameters γ, Γd_r_, Γtide and r (Section 4) lead to resonance locks, and what values do not. We define resonance locking in this context to be any behaviour such that σ does not increase without bound as t → ∞; this definition comprises locking into a stable fixed point (Fig. 4), limit cycles about an unstable fixed point (Fig. 5), and chaotic behaviour like in Fig. 7, but does not include the behaviour seen e.g. in Fig. 6.

We see by inspecting equations (25) and (43) that the adiabatic solution exactly passes through the resonance locking fixed point. As a result, a sufficient condition for a resonance lock to be achieved is |δωf| > |δωad|, which evaluates to

\begin{equation} \Gamma _{\mathrm{d}r}^2 \lesssim (1-r)|\Gamma _\mathrm{tide}| \gamma , \end{equation}

(50)

together with fixed point stability (Section 5.2; equation 40). However, this is a very conservative estimate of the resonance locking regime. To develop a set of necessary and sufficient criteria for resonance locks, recall that the time derivative of δω is given by (equation 24)

\begin{equation} \dot{\delta \omega }/\omega = -\Gamma _{\mathrm{d}r}+ \Gamma _\mathrm{tide}(g_1-rg_2), \end{equation}

(51)

where _g_1 and _g_2 are

\begin{equation} g_1 = \frac{\gamma }{\omega }|Q|^2 \qquad g_2=\mathrm{Im}[Q]. \end{equation}

(52)

First, assume r ≫ 1. By equation (30), we see that Γd_r_ < 0 and thus that the system approaches resonance from the left (in terms of δω; note that in Figs 4–7 the abscissa is σ − ω = −δω). Since the lagging fixed point is always stable in this situation (Section 5.2.2), we simply need the resonance passage to provide sufficient amplitude to reach the fixed point in order for a stable resonance lock to take hold. Invoking the no-backreaction approximation results from Section 6.1.2 to substitute a maximum value of Im[_Q_] into equation (51), dropping _g_1, and setting |$\dot{\delta \omega }\sim 0$| leads to the following condition for resonance locking:

\begin{equation} -\frac{\Gamma _{\mathrm{d}r}}{\Gamma _\mathrm{tide}} > 3r\sqrt{-\frac{\pi \omega }{\Gamma _{\mathrm{d}r}}},\qquad (r\gg 1). \end{equation}

(53)

Although our analysis is strictly valid only for r much larger than unity, we find it to work well even for r ≳ 1, as can be seen in the bottom row of Fig. 10 (where r = 1.5). We have inserted an additional factor of 3 on the right-hand side of equation (53) to match our numerical results.

Next, assume r ≪ 1. By equation (30) this means Γd_r_ > 0, so the system approaches the resonance locking fixed point from the right (in terms of δω). Here, however, the fixed point is not always stable, as we found in Section 5.2.2. In Section 6.2, we argued that resonance locks in the form of limit cycles or chaotic trajectories could occur when the fixed point was an unstable spiral (with a complex conjugate pair of eigenvalues), but that resonance capture failed when the fixed point was an unstable node (with all real eigenvalues); this was true within a factor of ∼3 in terms of the value of Γd_r_. Thus a necessary condition for resonance locks (to within a factor of ∼3) is equation (42), which specifies where the fixed point is a spiral. Next, we drop _g_2 and again substitute the no-backreaction approximation results from Section 6.1.2 to find that resonance passage can deliver a system to the resonance locking fixed point if

\begin{equation*} \frac{\gamma }{\Gamma _{\mathrm{d}r}} > \frac{\Gamma _{\mathrm{d}r}}{2\pi \Gamma _\mathrm{tide}}. \end{equation*}

Augmented with equation (42), this approximately becomes the following condition for resonance locking:

\begin{equation} \frac{\gamma }{\Gamma _{\mathrm{d}r}} > \frac{\Gamma _{\mathrm{d}r}}{6\pi \Gamma _\mathrm{tide}} + \frac{1}{3}\left(\frac{r^{2/3}}{1-r}\right)\left( \frac{\omega }{\Gamma _\mathrm{tide}} \right)^{1/3}\!\!\!,\qquad (r\ll 1). \end{equation}

(54)

Similar to our formula for r > 1, our analysis is valid only for r very close to zero; however, we again find it to work well even for r ≲ 1, as can be seen in the right-hand panels of the top two rows of Fig. 10 (where r = 0.5). We have inserted an additional factor of 1/3 on the right-hand side of equation (54) to match our numerical results.

8 TIDAL EVOLUTION DURING RESONANCE LOCKS

8.1 Accelerating tidal evolution

Here, we generalize the energetic arguments made in Burkart et al. (2013) to estimate the orbital and spin evolution during a resonance lock. During a lock, the reduced mode amplitude Q is roughly given by its value at the fixed point, i.e. equation (27); note, however, that this approximation is very crude for limit cycles and chaotic orbits (e.g. Figs 5 and 7) since in those cases the real and imaginary parts of Q have a more complicated dependence on time. Using equa-tions (16) and (27) together with our definitions of Γd_r_, Γtide and r from equations (18)–(21), we can derive

\begin{equation} \dot{\Omega }_\mathrm{orb}= \left( \frac{1}{1-r} \right)\left[ \alpha _\mathrm{orb}- \frac{r}{k}\left( m \alpha _\mathrm{spin}+\frac{\mathrm{\partial} \omega }{\mathrm{\partial} t}\right)\right], \end{equation}

(55)

where again αorb and αspin represent the contributions to |$\dot{\Omega }_\mathrm{orb}$| and |$\dot{\Omega }_\mathrm{spin}$| from slowly varying processes other than resonant interaction with the normal mode in question (see equations 13 and 16). This can be converted into an energy transfer rate by equation (15). Performing a similar derivation for the spin frequency, we have

\begin{equation} \dot{\Omega }_\mathrm{spin}= \left( \frac{r}{1-r} \right)\left[ - \alpha _\mathrm{spin}+ \frac{1}{rm}\left( k \alpha _\mathrm{orb}-\frac{\mathrm{\partial} \omega }{\mathrm{\partial} t}\right)\right]. \end{equation}

(56)

Examining equations (55) and (56), we see that a resonance lock acts to accelerate the orbital and spin evolution given by the non-resonant processes contributing to αorb and αspin, which are due e.g. to gravitational wave orbital decay, the equilibrium tide, etc. Moreover, since the time derivative of the eccentricity is simply a linear combination of |$\dot{\Omega }_\mathrm{orb}$| and |$\dot{\Omega }_\mathrm{spin}$| (Witte & Savonije 1999), resonance locking also accelerates circularization. The degree of acceleration depends on how close the moment of inertia ratio r is to unity; we estimate under what conditions r ∼ 1 in Section 8.2. This acceleration of tidal evolution is what led Witte & Savonije (2002) to conclude that resonance locks solve the solar-binary problem (Meibom, Mathieu & Stassun 2006), although they neglected essential non-linear effects that obviate their results (see Section 3).

The presence of αspin in the evolution equation for |$\dot{\Omega }_\mathrm{orb}$| (and αorb in the equation for |$\dot{\Omega }_\mathrm{spin}$|⁠) implies that a resonance lock efficiently couples orbital and spin evolution together, as well as to stellar evolution through the rate of change of the eigenfrequency ∂ω/∂t. For example, if gravitational waves in an inspiralling white dwarf binary cause orbital decay, a resonance lock will cause tidal synchronization to occur on a gravitational wave time-scale (Section 9.1). Similarly, resonance locking can cause stellar spin-down by magnetic braking or eigenfrequency evolution due to tidal heating to backreact on the binary orbit.

8.2 Conditions for rapid tidal evolution

Here, we estimate the ‘typical’ value of r (defined in equation 21) expected to occur in a given binary, so as to assess whether resonance locks significantly accelerate tidal evolution (Section 8.1). First, consider a binary in a circular, spin-aligned orbit. For the lowest order l = 2 spherical harmonic of the tidal potential, only one temporal Fourier component of the tidal forcing exists: k = m = 2 (Section 4.1). With only a single forcing component, it is thus generically unlikely to find r ∼ 1. The exception is when considering a star whose companion's mass is ≪M, in which case r may be ∼1 even for a circular orbit.

For an eccentric orbit, however, there is significant power at many harmonics k of the orbital frequency. In this case, a resonance lock persists until it is disrupted when another mode, driven by a different Fourier component, also comes into resonance and upsets the balance of spin and orbital evolution theretofore enforcing |$\dot{\delta \omega }=0$|⁠. Whether the existing resonance lock can withstand a second resonance passage depends on how ‘robust’ it is; we can estimate the degree of ‘robustness’ using the maximum amplitude achieved under the no-backreaction approximation from Section 6.1.2:4

\begin{equation} |q|_{\mathrm{max}} = |U|\sqrt{\frac{2\pi \omega }{|\Gamma _{\mathrm{d}r}|}}. \end{equation}

(57)

If we equate ω ≈ _k_Ωorb, having taken |_m_Ωspin| ≪ |_k_Ωorb| for simplicity, then |q|max depends on the harmonic index k as

\begin{equation} |q|_{\mathrm{max}} \mathop {\propto }\limits _{\sim }k^b X_{lm}^k, \end{equation}

(58)

where X is a Hansen coefficient and b > 0 is a constant that accounts for power-law dependences on the tidal overlap integral and other mode-dependent quantities entering into equation (57) (Section 4.1). Using our scaling derived in Appendix B, equation (58) becomes

\begin{equation} |q|_{\mathrm{m}ax}\mathop {\propto }\limits _{\sim }k^{b-1/2} \exp \left[ -k g(e) \right], \end{equation}

(59)

where the full form of |$g(e)\approx (1-\text{e}^2)^{3/2}/3$| is given in equa-tion (B4).

The longest lived resonance locks in an eccentric binary will be those with the largest values of |q|m_ax_. Thus, we can estimate the ‘typical’ value of r by finding the value of k that maximizes |q|m_ax_. This is

\begin{equation} \mathrm{argmax}_k |q|_{\mathrm{max}} \approx \frac{b-1/2}{g(e)}. \end{equation}

(60)

In order to produce a simple order-of-magnitude estimate, we take b − 1/2 ∼ 1 to find that the most robust resonance locks will have r ∼ 1, and thus significantly increase the rate of tidal evolution, if5

\begin{equation} 1-\text{e}^2 \sim \left(10\frac{I_*}{\mu a^2}\right)^{1/3}. \end{equation}

(61)

When this condition is not satisfied, then either the longest lived resonance locks will not provide much acceleration of tidal evolution, or there will be no long-lived locks at all.

9 ASTROPHYSICAL APPLICATIONS

In Section 9.1 and 9.2, we determine where tidal resonance locks may be able to occur by applying the criteria we have developed – equations (29) and (54) – to inspiralling compact object binaries and eccentric stellar binaries.

9.1 Inspiralling compact object binaries

We consider the case of a circular, spin-orbit-aligned binary consisting of either two white dwarfs or two neutron stars. The equivalent of the drift rate Γd_r_ in this case is

\begin{equation} \Gamma _{\mathrm{d}r}=\frac{m\Omega _\mathrm{orb}}{\omega t_{\mathrm{g}w}}, \end{equation}

(62)

where t_g_w is the gravitational wave orbital decay time given by (Peters 1964)

\begin{equation} t_{\mathrm{g}w}= \omega _*^{-1}\frac{5}{96}\frac{\left( 1+M^{\prime }/M \right)^{1/3}}{M^{\prime }/M}\beta _*^{-5}\left( \frac{\omega _*}{\Omega _\mathrm{orb}} \right)^{8/3}. \end{equation}

(63)

Here, |$\omega _*^2=GM/R^3$| is the dynamical frequency and |$\beta _*^2=GM/Rc^2$| is the ratio of the escape velocity to the speed of light. We consider only the l = |m| = 2 component of the tidal response.

For a fiducial double-white dwarf binary, we use the |$0.6\,\text{M}_{\odot }$|⁠, T_e_ff = 5500 K carbon/oxygen white dwarf model described in Burkart et al. (2013), with a radius of |$R=0.013\,\text{R}_{\odot }$|⁠, and a moment of inertia of I* = 0.18_MR_2. For the damping rate γ and tidal overlap integral |$\mathcal {I}$|⁠, we use the following scaling relations listed in table 3 of Burkart et al. (2013):

\begin{eqnarray*} \mathcal {I}&\sim & 27\times \left(\displaystyle\frac{\sigma }{\omega _*}\right)^{3.69}\\ \gamma &\sim & 2.9\times 10^{-14}\omega _*\times \left(\displaystyle\frac{\sigma }{\omega _*}\right)^{-1.88}, \end{eqnarray*}

where our mode normalization convention is E* = ε and we are neglecting rotational modifications of the stellar eigenmodes (in other words, setting the Coriolis force operator B from Section 4.1 to zero).

For a fiducial double-neutron star binary, we use the |$M=1.4\,\text{M}_{\odot }$|⁠, R = 12 km cold neutron star model employed in Weinberg, Arras & Burkart (2013), which assumed the Skyrme–Lyon equation of state (Chabanat et al. 1998; Steiner & Watts 2009). We assume that I* = 0.18_MR_2, as with our white dwarf model. Weinberg et al. (2013) give the following scaling relations for l = 2 g modes:

\begin{eqnarray*} \mathcal {I}&\sim & 0.3\times \left(\displaystyle\frac{\sigma }{\omega _*}\right)^{2}\\ \gamma &\sim & 4\times 10^{-8}\times \left(\displaystyle\frac{\sigma }{\omega _*}\right)^{-2}T_8^{-2}\ \mathrm{Hz}, \end{eqnarray*}

where _T_8 is the core temperature in units of 108 K.

Fig. 11 shows our result for the white dwarf case. The yellow region shows where the resonance locking fixed point exists (equa-tion 29) but where resonance locks are nonetheless not possible according to our numerical and analytic results in Section 7. The orange region shows the complete resonance locking regime derived from our numerical results (equation 54). The dashed black line in the top panel is a simple, schematic system trajectory for a double-white dwarf binary. The system begins in the upper right with a long orbital period and a small rotation rate. Once orbital decay causes the system to reach the orange region, a resonance lock occurs and the forcing frequency σ is held approximately constant. This was already demonstrated in Burkart et al. (2013) using the adiabatic approximation (Section 6.1.1), which is valid for determining when the resonance locking fixed points exist.

Figure 11.

Resonance locking regimes for a fiducial white dwarf binary. Yellow shading indicates that the resonance locking fixed points exists (equation 29), but that resonance locks are nonetheless not possible. Orange shading indicates the complete resonance locking regime, derived from our numerical results in Section 7.1 (equation 54). The black dashed line shows a schematic system trajectory, assuming the initial spin is negligible. The dot–dashed magenta line shows where the lagging fixed point becomes unstable (equation 40).

Eventually, however, the system exits the orange region. At this point resonance locking is no longer possible because of the short gravitational wave inspiral time. This novel prediction comes from our analysis in this work.6 Note that this prediction neglects non-linear hydrodynamical phenomena (Section 3): in reality, the wave amplitude eventually becomes large enough to cause wave breaking, as shown in Burkart et al. (2013), Fuller & Lai (2012b). This may be a more stringent constraint on the existence of resonance locks in many close white dwarf binaries.

In the neutron star case, the gravitational wave time is much shorter than for white dwarf binaries at comparable values of R/a, since neutron stars are much more relativistic objects. Resonance passage thus happens very quickly, which prevents modes from reaching amplitudes large enough to allow locking. As a result, resonance locks are never possible, even though the resonance locking fixed point exists for _P_orb ≲ 50 ms.

To get a sense of the degree to which equation (54) fails to be satisfied for neutron star binaries, we first compute the following quantities with Ωspin = 0:

\begin{eqnarray*} &&\!\!\!\!&&{\frac{\gamma }{\omega } = 10^{-7}T_8^{-2} \left( \frac{P_\mathrm{orb}}{50\ \mathrm{ms}} \right)^3 \quad \frac{\Gamma _{\mathrm{d}r}}{\omega }=6\times 10^{-5}\left( \frac{P_\mathrm{orb}}{50\ \mathrm{ms}} \right)^{-5/3}}\\ &&\!\!\!\!\frac{\Gamma _\mathrm{tide}}{\omega } = 10^{-11}\left( \frac{P_\mathrm{orb}}{50\ \mathrm{ms}} \right)^{-6}\quad r = 0.002\left( \frac{P_\mathrm{orb}}{50\ \mathrm{ms}} \right)^{-4/3}. \end{eqnarray*}

We thus set r ≈ 0. Substituting the remaining values into equa-tion (54) and simplifying, we have

\begin{equation} 10\,T_8^2 < 10^{-7}\left( \frac{P_\mathrm{orb}}{50\ \mathrm{ms}} \right)^{1/3}. \end{equation}

(64)

This shows that resonance locking fails to occur by ∼8 orders of magnitude at a wide range of orbital periods. The conclusion that resonance locks cannot occur in neutron star binaries is thus very robust.

9.2 Eccentric binaries

In this section, we estimate whether resonance locking can occur in eccentric stellar binaries (Witte & Savonije 1999). For a fiducial system, we take parameters consistent with KOI-54 (Welsh et al. 2011), where it has been recently suggested that one or more of the observed tidally excited pulsations may be the signatures of resonance locks (Burkart et al. 2012; Fuller & Lai 2012a; see however O'Leary & Burkart 2014). This system consists of two similar A stars with |$M\approx 2.3\,\text{M}_{\odot }$|⁠, |$R\approx 2.2\,\text{R}_{\odot }$| and T_e_ff ≈ 8500 K. The binary's orbital parameters are e = 0.83 and _P_orb = 43 d. Burkart et al. (2012) estimated damping rates for such stars to be

\begin{equation} \gamma \sim 0.1 \left( \frac{\omega }{\omega _*} \right)^{-4}\ \mathrm{Myr}^{-1}. \end{equation}

(65)

For overlap integrals |$\mathcal {I}$|⁠, we take the following scaling derived for A stars from Burkart et al. (2012):

\begin{equation} \mathcal {I}\sim 10^{-5}\left( \frac{\omega }{\omega _*} \right)^{11/6}. \end{equation}

(66)

The drift rate Γd_r_ in this case comes from the equilibrium tide's influence on each star's spin and on the overall orbital frequency. We account only for the equilibrium tide's effect on the orbital frequency for simplicity. Parametrizing the equilibrium tide's energy transfer rate by its quality factor |$\mathcal {Q}_{\mathrm{eq}}$|⁠, so that (Goldreich & Soter 1966)

\begin{equation} |\dot{E}_{\mathrm{eq}}|\sim \frac{E_{\mathrm{tide}}\Omega _\mathrm{orb}}{\mathcal {Q}_{\mathrm{eq}}}, \end{equation}

(67)

we have the approximate formula

\begin{equation} \Gamma _{\mathrm{d}r}\sim \frac{E_{\mathrm{tide}}}{\mathcal {Q}_{\mathrm{eq}}\mu a^2}\frac{k}{\omega }, \end{equation}

(68)

where the energy contained in the equilibrium tide is roughly _E_tide ∼ λ(M_′/M)2(R/a)6_E*, and the constant λ (related to the apsidal motion constant) is ∼3 × 10−3 for an A star.

Witte & Savonije (1999) invoked the adiabatic approximation (Section 6.1.1) to show that resonance locking could occur in various fiducial eccentric binary systems; Burkart et al. (2012) performed a similar analysis for KOI-54. However, as we have established in this work, this only means that the resonance locking fixed point exists, and not necessarily that resonance locking actually occurs. For the latter, we need to apply our criterion from equation (54).

As in Section 9.1, we proceed to compute the values of our four parameters that affect the possibility of resonance locks in the current situation. We fix r, but assume 0 ≪ r < 1. We also assume that the star is non-rotating. We then find7

\begin{eqnarray*} \frac{\gamma }{\omega } &=& 10^{-11} \left( \displaystyle\frac{P_\mathrm{orb}}{40\ \mathrm{d}} \right)^{5/3}\\ \frac{\Gamma _{\mathrm{d}r}}{\omega } &=& 3\times 10^{-21}\left( \displaystyle\frac{P_\mathrm{orb}}{40\ \mathrm{d}} \right)^{-4} \left( \displaystyle\frac{\mathcal {Q}_{\mathrm{e}q}}{10^8} \right)^{-1}\\ \frac{\Gamma _\mathrm{tide}}{\omega } &=& 6\times 10^{-18} \left(\displaystyle\frac{P_\mathrm{orb}}{40\ \mathrm{d}}\right)^{-4.6}. \end{eqnarray*}

Substituting into equation (54), noting that the first term on the right-hand side of equation (54) is much smaller than the second in this case (unlike in Section 9.1), we find that locks are present if

\begin{equation} 1 > 10^{-4}\left( \frac{1-r}{0.1} \right)^{-1}\left( \frac{P_\mathrm{orb}}{40\ \mathrm{d}} \right)^{4.1}\left( \frac{\mathcal {Q}_{\mathrm{e}q}}{10^8} \right)^{-1}. \end{equation}

(69)

It thus appears that resonance locks are indeed possible in eccentric binaries, subject to the validity of the assumptions enumerated in Section 3 (e.g. solid-body rotation).

10 CONCLUSION

We have studied tidally induced resonance locking in close (but detached) binary systems. In a resonance lock, the detuning frequency δω = ω − σ between a stellar or planetary eigenmode frequency ω and a particular Fourier harmonic of the tidal driving frequency σ = _k_Ωorb − _m_Ωspin is held constant (Section 2; Witte & Savonije 1999). This happens when a slowly varying physical process causing δω to evolve in one direction is balanced by resonant interaction with the eigenmode in question causing δω to evolve in the reverse direction. The slowly varying process could be, e.g. orbital decay due to gravitational waves causing Ωorb to increase, magnetic braking causing Ωspin to decrease, stellar evolution altering ω, or non-resonant components of the tidal response (the ‘equilibrium tide’) affecting both Ωorb and Ωspin simultaneously.

Our primary goal has been to understand the dynamical properties and stability of resonance locks without relying on simplifying approximations for the mode amplitude evolution used in previous calculations. We defer detailed implications of these results to future papers. We have derived a novel set of equations allowing for a dynamically evolving mode amplitude coupled to the evolution of both Ωspin and Ωorb (Section 4). In particular, we do not assume that the mode amplitude is given by a Lorentzian profile resulting from the adiabatic approximation (Section 6.1.1) used in previous work, but instead solve the fully time-dependent mode amplitude equation.

In Section 5, we analysed the stability of the dynamical fixed points associated with resonance locks. We analytically derived when such fixed points exist (equation 29); there are either two fixed points or none for a given eigenmode. Although one of these equilibria is always unstable, the other can be stable when certain restrictions on binary and mode parameters are met (equation 40 in Section 5.2). One of the important conclusions of this analysis is that resonance locks can exist and be stable even when the adiabatic approximation for the mode amplitude evolution is invalid (which happens, e.g. in the limit of moderately weak damping).

In Section 6.2, we analysed the properties of resonance locks using direct numerical integration of our dynamical equations. In the simplest case in which a resonance lock fixed point exists and is stable, two possibilities arise: either a resonance passage is able to pump the mode's amplitude up sufficiently high to reach the fixed point and be captured into it, creating the resonance lock (Fig. 4), or the system instead sweeps through resonance without locking.

The more interesting situation is when both fixed points are unstable. In this case, we showed that resonance locking can nonetheless occur in some cases in a time-averaged sense. In these situations, the mode amplitude and detuning frequency δω execute limit cycles or even chaotic trajectories around the fixed points (see Figs 5 and 7). We presented evidence in Section 6.3 suggesting that resonance locking may in fact correspond to a strange attractor for certain parameter values; see Fig. 8.

In order to determine when resonance locking of some kind occurs (either stable, limit cycle, or chaotic), we performed numerical integrations over wide ranges of parameter values in Section 7.1. Using analytic approximations from Section 6.1, we then developed approximate analytic formulae that explain our numerical results and define the binary and mode parameter regimes in which resonance locks of some kind occur. The key results are equations (53) and (54) and Fig. 10. Future studies of tidal evolution that do not include our full set of coupled mode-spin-orbit evolution equations can nonetheless utilize our results to assess whether resonance locks can occur.

One of the interesting consequences of resonance locks highlighted by Witte & Savonije (1999) and Witte & Savonije (2002) using numerical simulations in the adiabatic approximation is that locks can produce a significant speed up of orbital and spin evolution. We have explained this analytically in Section 8.1. In particular, we have demonstrated that resonance locks generically act to produce orbital and spin evolution on a time-scale that is somewhat shorter than the slowly varying physical process whose influence drives the system into a lock. The magnitude of this acceleration depends on the effective moment of inertia ratio r defined in equa-tion (21) and is large for r ∼ 1. The latter condition can only be satisfied in eccentric binaries or binaries with high-mass ratios. For the case of an eccentric orbit, we derived a rough condition for significant acceleration of the rate of orbital and spin evolution in Section 8.2.

To give a rough sense of the possible application of our results, we applied them to three sample astrophysical systems in Section 9: inspiralling white dwarf and neutron star binaries in Section 9.1, and eccentric binaries with early-type stars in Section 9.2. As has been argued previously using the adiabatic approximation for mode amplitudes, resonance locks are likely very common in white dwarf binaries and eccentric stellar binaries. They cannot, however, occur in neutron star binaries since orbital decay by gravitational wave emission is too rapid. A future application that may be of considerable interest is tidal circularization during high-eccentricity migration of hot Jupiters.

The theory of resonance locking that we have developed bears some similarity to resonance capture in planetary dynamics (Murray & Dermott 1999). In the case of both mean-motion resonances (Goldreich 1965) and spin-orbit resonances (Goldreich & Peale 1968), the generic equation governing the evolution of the relevant angle Ψ towards resonance is

\begin{equation} \ddot{\Psi } =- F\sin \Psi + G. \end{equation}

(70)

For mean-motion resonances, Ψ defines the angle between the mean anomalies of two orbiting bodies, while for spin-orbit resonances, Ψ is related to the difference between a body's mean and true anomalies (relevant only for eccentric orbits). In both cases, G provides a frequency drift term, resulting from orbital decay for mean-motion resonances and from a net tidal torque for spin-orbit resonances.

equation (70) can be compared to our equation describing the evolution of the detuning frequency δω = ω − σ in equation (24). Both specify the second time derivative of a resonance angle, and both contain a frequency drift term describing how resonance is approached. The essential difference is that in place of the pendulum restoring force present in equation (70), tidal resonance locking instead contains two terms providing the complicated interaction with a stellar or planetary normal mode. Thus, although there are qualitative similarities between resonance capture in planetary dynamics and resonance locking, no formal mathematical analogy exists.

It is important to reiterate that resonance locks can only occur under a specific set of conditions (Section 3). They are not relevant to all close binaries. In particular, the dynamical tide must be composed of global radial standing waves, with damping times much longer than radial wave travel times. For this reason, an efficient angular momentum transport process must maintain approximate solid body rotation; if not, critical layers may develop where mode angular momentum is deposited into the background rotation profile, which would lead to efficient local wave damping. In addition, we have restricted our analysis to linear perturbation theory. In practice, this represents a restriction on the maximum mode amplitudes that are allowed, since non-linear instabilities can act on large-amplitude waves. For example, in the case of binaries containing solar-type stars with radiative cores, wave breaking in the core likely prohibits the establishment of global standing waves (Goodman & Dickson 1998), thus also precluding resonance locks from developing.

The authors are pleased to thank Edgar Knobloch, Keaton Burns, Eugene Chiang, and Scott Tremaine for useful discussions and guidance. JB is an NSF Graduate Research Fellow. EQ was supported by a Simons Investigator award from the Simons Foundation, the David and Lucile Packard Foundation, and the Thomas Alison Schneider Chair in Physics at UC Berkeley. PA is an Alfred P. Sloan Fellow, and received support from the Fund for Excellence in Science and Technology from the University of Virginia.

1

It can be shown that the second term in equation (11) is exactly the angular momentum transfer rate from the orbit; see e.g. Weinberg et al. (2012). This justifies attributing the first term to changes in the spin.

2

Tidal heating could contribute to the |q|2 term in equation (17), due to heat deposited by the mode in question affecting the background star or planet; we neglect this for simplicity.

3

Numerous other potential resonance lock criteria exist; however, we have found this to be the most reliable.

4

Our results are insensitive to the expression by which a resonance lock's ‘robustness’ is quantified, so long as it is proportional to U.

5

Note that whether r > 1 or r < 1 is determined by the sign of Γd_r_, which is independent of our estimates here; see Section 5.1.

6

Burkart et al. (2013) did assess where the adiabatic approximation (which was inaccurately referred to as the ‘secular approximation’) became invalid; however, this is an incomplete consideration that fails to account for fixed point instability, limit cycles, etc.

7

Since we are considering an eccentric orbit, the tidal coupling coefficient U (on which Γtide depends) is additionally proportional to a Hansen coefficient, due to the Fourier expansion of the orbital motion (Section 4.1). We take this coefficient to be of order unity.

REFERENCES

,

Handbook of Mathematical Functions

,

1972

New York

Dover Press

,

MNRAS

,

2010

, vol.

404

pg.

1849

,

MNRAS

,

2012

, vol.

421

pg.

983

,

MNRAS

,

2013

, vol.

433

pg.

332

,

Nucl. Phys. A

,

1998

, vol.

635

pg.

231

,

Proc. R. Soc. A

,

1979

, vol.

368

pg.

389

,

ApJ

,

1978a

, vol.

221

pg.

937

,

ApJ

,

1978b

, vol.

222

pg.

281

,

MNRAS

,

2011

, vol.

412

pg.

1331

,

MNRAS

,

2012a

, vol.

420

pg.

3126

,

MNRAS

,

2012b

, vol.

421

pg.

426

,

MNRAS

,

1965

, vol.

130

pg.

159

,

ApJ

,

1989a

, vol.

342

pg.

1075

,

ApJ

,

1989b

, vol.

342

pg.

1079

,

ARA&A

,

1968

, vol.

6

pg.

287

,

Icarus

,

1966

, vol.

5

pg.

375

,

ApJ

,

1998

, vol.

507

pg.

938

,

Table of Integrals, Series, and Products

,

2007

2nd edn

Amsterdam

Elsevier

,

MNRAS

,

1994

, vol.

270

pg.

611

,

MNRAS

,

2012

, vol.

423

pg.

486

,

Mon. Weather Rev.

,

1966

, vol.

94

pg.

295

,

ApJ

,

2006

, vol.

653

pg.

621

,

Solar System Dynamics

,

1999

Cambridge

Cambridge Univ. Press

,

MNRAS

,

2014

, vol.

440

pg.

3036

,

Phys. Rev.

,

1964

, vol.

136

pg.

1224

,

ApJ

,

1977

, vol.

213

pg.

183

,

MNRAS

,

2005

, vol.

357

pg.

834

,

ApJ

,

1994

, vol.

426

pg.

688

,

Phys. Lett. A

,

1976

, vol.

57

pg.

397

,

Phys. Rev. D

,

2002

, vol.

65

pg.

024001

,

Phys. Rev. Lett.

,

2009

, vol.

103

pg.

181101

,

ApJ

,

2012

, vol.

751

pg.

136

,

ApJ

,

2013

, vol.

769

pg.

121

et al. ,

ApJS

,

2011

, vol.

197

pg.

4

,

Introduction to Applied Nonlinear Dynamical Systems and Chaos

,

2003

2nd edn

Berlin

Springer-Verlag

,

A&A

,

1999

, vol.

350

pg.

129

,

A&A

,

2002

, vol.

386

pg.

222

,

A&A

,

1975

, vol.

41

pg.

329

APPENDIX A: DERIVING FIXED POINT STABILITY CONDITIONS

Here, we will determine the stability region for the resonance locking fixed point described in Section 5. The characteristic polynomial P in question is given in equation (37).

First, _P_1 > 0 reduces to

\begin{equation} \delta \omega _{\mathrm{f}}< \frac{(1-r)^2\gamma ^2\Gamma _\mathrm{tide}}{r\Gamma _{\mathrm{d}r}^2}. \end{equation}

(A1)

We immediately see that this is satisfied if Γd_r_ < 0 (assuming Γtide > 0), since by equation (39) we then have δωf < 0. If instead δωf > 0, then we can take |δωf| ≫ γ to derive

\begin{equation} \Gamma _{\mathrm{d}r}<\gamma \left( \frac{1-r}{r^{2/3}} \right)\left( \frac{\Gamma _\mathrm{tide}}{\omega } \right)^{1/3}. \end{equation}

(A2)

Next, the Hurwitz matrix for P is

\begin{equation} H={\left(\begin{array}{cc}2\gamma &\quad P_0 \\ 1 &\quad P_1 \\ 2\gamma &\quad P_0 \end{array}\right)}. \end{equation}

(A3)

The remaining condition for stability of the fixed point is that the three leading principal minors of H must be positive; this formally yields three additional inequalities. However, one is γ > 0 which is always satisfied, and the other two are actually identical:

\begin{equation} 2\gamma P_1 > P_0, \end{equation}

(A4)

which expands to

\begin{equation} \gamma ^2\Gamma _\mathrm{tide}(1-r)^2 > \Gamma _{\mathrm{d}r}^2\delta \omega _{\mathrm{f}}. \end{equation}

(A5)

We now see that Γd_r_ < 0 implies asymptotic stability. If Γd_r_ > 0, then the inequality reduces to (again assuming |δωf| ≫ γ)

\begin{equation} \Gamma _{\mathrm{d}r}<\gamma \left( 1-r \right)\left( \frac{\Gamma _\mathrm{tide}}{\omega } \right)^{1/3}. \end{equation}

(A6)

Since equation (A6) is more restrictive than equation (A2) (since r < 1 for Γd_r_ > 0; see equation 30), equation (A6) is the condition for asymptotic stability when Γd_r_ > 0.

APPENDIX B: HANSEN COEFFICIENT SCALING

Here, we will determine how the Hansen coefficients |$X_{lm}^k$| scale with harmonic index k. These coefficients are defined to satisfy

\begin{equation} \left( \frac{a}{D(t)} \right)^{l+1}\text{e}^{-\text{i}mf(t)}=\sum _{k=-\infty }^\infty X_{lm}^k \text{e}^{-\text{i}k\Omega _\mathrm{orb}t}, \end{equation}

(B1)

where D is the binary separation and f is the true anomaly. For |k| ≫ |m|Ωperi/Ωorb, where Ωperi is the effective orbital frequency at periapse, we have that |$X_{lm}^k\mathop {\propto }\limits _{\sim }X_{00}^k$|⁠. We can express |$X_{lm}^k$| in general as an integral over the eccentric anomaly E (Burkart et al. 2012); for l = m = 0, this is

\begin{eqnarray} X_{00}^k &=& \frac{1}{\pi }\int _{0}^{\pi} \cos \left[ k(E-e\sin E) \right] {\rm d}E \nonumber\\ &=& J_k(ek) \end{eqnarray}

(B2)

where J is a Bessel function. We can expand J k(ek) as (Abramowitz & Stegun 1972)

\begin{equation} J_k(ek) \mathop {\propto }\limits _{\sim }\frac{1}{\sqrt{k}} \exp \left[ -k g(e)\right], \end{equation}

(B3)

where

\begin{eqnarray} g(e) &=& \displaystyle\frac{1}{2}\ln \left( \frac{1+\eta }{1-\eta } \right)-\eta \nonumber\\ &=& \displaystyle\frac{\eta ^3}{3} + \frac{\eta ^5}{5} + \frac{\eta ^7}{7} + \cdots \end{eqnarray}

(B4)

and |$\eta = \sqrt{1-\text{e}^2}$|⁠.

APPENDIX C: CANONICAL ANGULAR MOMENTUM

Here, we will derive the canonical angular momentum associated with a stellar perturbation. The Lagrangian density for a stellar perturbation is (Friedman & Schutz 1978a)

\begin{equation} \mathcal {L} = \frac{1}{2}\rho \left( |\dot{\boldsymbol {\xi }}|^2 + \dot{\boldsymbol {\xi }}\cdot B\boldsymbol {\xi } - \boldsymbol {\xi }\cdot C\boldsymbol {\xi }\right), \end{equation}

(C1)

where |$\boldsymbol {\xi }$| is the Lagrangian displacement vector, the Coriolis force operator B was defined in Section 4.1, and the Hermitian operator C is defined in e.g. Schenk et al. (2002). The z component of the canonical angular momentum is then

\begin{eqnarray} J &=& -\left\langle \mathrm{\partial} _\phi \boldsymbol {\xi },\,\displaystyle\frac{\mathrm{\partial} \mathcal {L}}{\mathrm{\partial} \dot{\boldsymbol {\xi }}}\right\rangle \nonumber\\ &=& -\left\langle \mathrm{\partial} _\phi \boldsymbol {\xi },\, \dot{\boldsymbol {\xi }}+B\boldsymbol {\xi }/2\right\rangle . \end{eqnarray}

(C2)

As in Section 4.1, we perform a phase space expansion of the Lagrangian displacement vector and its time derivative (Schenk et al. 2002), so that

\begin{equation} {\left(\begin{array}{c}\boldsymbol {\xi }\\ \dot{\boldsymbol {\xi }} \end{array}\right)} = \sum _A q_A {\left(\begin{array}{c}\boldsymbol {\xi }_A\\ -i\omega _A\boldsymbol {\xi }_A \end{array}\right)}, \end{equation}

(C3)

where A runs over all rotating-frame stellar eigenmodes and both frequency signs. The canonical angular momentum then becomes

\begin{equation} J = \frac{1}{2} \sum _{AB} q_A^* q_B^{}\, m_A^{} \left[ 2\omega _B\left\langle \boldsymbol {\xi }_A,\boldsymbol {\xi }_B\right\rangle + \left\langle \boldsymbol {\xi }_A,iB\boldsymbol {\xi }_B\right\rangle \right]. \end{equation}

(C4)

Since |$\left\langle \mathrm{\partial} _\phi \boldsymbol {\xi },\, \dot{\boldsymbol {\xi }}\right\rangle$| is real valued, we can re-express equation (C4) as

\begin{eqnarray} J = \frac{1}{2} \sum _{AB} q_A^* q_B^{} \left[ (m_A\omega _B+m_B\omega _A)\left\langle \boldsymbol {\xi }_A,\boldsymbol {\xi }_B\right\rangle + m_A\left\langle \boldsymbol {\xi }_A,iB\boldsymbol {\xi }_B\right\rangle \right]. \nonumber\\ \end{eqnarray}

(C5)

We have that |$\left\langle \boldsymbol {\xi }_A,\boldsymbol {\xi }_B\right\rangle \propto \delta _{m_A,m_B}$|⁠, so J reduces to

\begin{equation} J = \frac{1}{2} \sum _{A} \frac{\varepsilon _A m_A}{\omega _A}|q_A|^2, \end{equation}

(C6)

where we have used the orthogonality relation (Friedman & Schutz 1978b)

\begin{equation} (\omega _A+\omega _B)\left\langle \boldsymbol {\xi }_A,\boldsymbol {\xi }_B\right\rangle + \left\langle \boldsymbol {\xi }_A,iB\boldsymbol {\xi }_B\right\rangle = \delta _{AB} \frac{\varepsilon _A}{\omega _A}. \end{equation}

(C7)

© 2014 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society